Direct estimation of kinetic parametric images for dynamic PET.

Dynamic positron emission tomography (PET) can monitor spatiotemporal distribution of radiotracer in vivo. The spatiotemporal information can be used to estimate parametric images of radiotracer kinetics that are of physiological and biochemical interests. Direct estimation of parametric images from raw projection data allows accurate noise modeling and has been shown to offer better image quality than conventional indirect methods, which reconstruct a sequence of PET images first and then perform tracer kinetic modeling pixel-by-pixel. Direct reconstruction of parametric images has gained increasing interests with the advances in computing hardware. Many direct reconstruction algorithms have been developed for different kinetic models. In this paper we review the recent progress in the development of direct reconstruction algorithms for parametric image estimation. Algorithms for linear and nonlinear kinetic models are described and their properties are discussed.

1 Introduction

1.1 Parametric Imaging with Dynamic PET

Dynamic positron emission tomography (PET) provides the additional temporal information of tracer kinetics compared to static PET. It has been shown that tracer kinetics can be useful in tumor diagnosis such as differentiation between malignant and benign lesions 1, 2, and in therapy monitoring where the kinetic values from dynamic PET can be more predictive for assessment of response than the standard uptake value at a single time point 3, 4.

There are basically two approaches for deriving tracer kinetics from dynamic PET data: region-of-interest (ROI) kinetic modeling and parametric imaging 5-7. The ROI based approach fits a kinetic model to the average time activity curve (TAC) of a selected ROI, so it is easy to implement and has low computation cost. In contrast, parametric imaging estimates kinetic parameters for every pixel and provides the spatial distribution of kinetic parameters. It is thus more suitable to study heterogeneous tracer uptake in tissue. Parametric images have been found useful in many biological research and clinical diagnosis 8-10. However, parametric imaging is more demanding computationally and more sensitive to noise in dynamic PET data than the ROI-based kinetic modeling. Obtaining parametric images of high quality therefore raises new challenges that were not faced by conventional ROI-based methods 11.

1.2 Direct Reconstruction of Parametric Images

The typical procedure of parametric imaging is to reconstruct a sequence of emission images from the measured PET projection data first, and then to fit the TAC at each pixel to a linear or nonlinear kinetic model. To obtain an efficient estimate, the noise distribution of the reconstructed activity images should be modeled in the kinetic analysis. However, exact modeling of the noise distribution in PET images can be difficult because noise is spatially variant and object-dependent. Usually the spatial variation and correlations between pixels are simply ignored in the kinetic fitting step, which can lead to very noisy results.

Direct reconstruction has been developed to reduce noise amplification in parametric imaging. Direct reconstruction methods combine tracer kinetic modeling and emission image reconstruction into a single formula to estimate parametric images directly from the raw projection data 12, 13. Because PET projection data are well-modeled as independent Poisson random variables, direct reconstruction allows accurate compensation of noise propagation from sinogram measurements to the kinetic fitting process. It has been shown that images reconstructed by direct reconstruction methods have better bias-variance characteristics than those obtained by indirect methods for both linear models 33, 34 and nonlinear kinetic models 25, 29.

One drawback of direct reconstruction of parametric images is that the optimization algorithms are more complex than indirect methods 12, 13, due to the intertwinement of the temporal and spatial correlations. In addition, when compartment models are used, the relationship between the kinetic parameters and PET data also becomes nonlinear. Therefore, the development of efficient optimization algorithms for direct reconstruction is of high significance.

1.3 A Brief History of Direct Reconstruction

The theory of direct reconstruction of parametric images was developed right after the publication of the maximum likelihood (ML) expectation-maximization (EM) algorithm for static PET reconstruction 62. In 1984, Snyder 14 presented an ML EM algorithm to estimate compartmental parameters from time-of-flight list-mode data. Carson and Lange 15 in 1985 presented an EM algorithm for maximum likelihood reconstruction of kinetic parameters from PET projection data. However, neither of these algorithms was validated using a realistic simulation or real data at the time, possibly due to the limitation of computing hardware.

In 1995, Limber et al 16 combined a least squares reconstruction with the Levenberg-Marquardt (LM) algorithm 17 to estimate the parameters of a single exponential decay model from SPECT projection data. Separable nonlinear least squares were used by Huesman et al 18 in 1998 to estimate kinetic parameters of a one-tissue compartment model from 3D SPECT projection data. To reduce computation cost, several researchers also developed methods for direct estimation of kinetic parameters for ROIs instead of reconstruction of the whole parametric images, e.g. 18-22. These algorithms are only efficient for a small-scale problem with a limited number of ROIs. When the number of unknown parameters goes beyond a few hundreds, both the LM algorithm used by Limber et al 16 and the separable nonlinear least squares used by Huesman et al 18 become inefficient.

The development of algorithms that are suitable for large-scale direct estimation began to attract more and more attention with the recent advances in computing technology. In 2005, a parametric iterative coordinate descent (PICD) algorithm was proposed by Kamasak et al 25 for penalized likelihood estimation of parametric images of a two-tissue compartment model. To our knowledge, this was the first demonstration of direct parametric reconstruction on a dense set of voxels. Yan et al 30 proposed a new EM algorithm for direct ML reconstruction of kinetic parameters of a one-tissue compartment model. Both Kamasak's PICD and Yan's EM algorithms are specific to their respective compartment models. To simplify practical implementation, Wang and Qi in 2008 proposed generalized optimization transfer algorithms for direct penalized likelihood reconstruction of parametric images that are applicable to a wide variety of kinetic models 28, 29. The algorithms resemble the empirical iterative implementation of the indirect method that alternates between an image reconstruction update and kinetic fitting step (e.g., 26, 27), but have the advantage of guaranteed monotonic convergence to the direct estimate. Later Wang and Qi 31, 32 also proposed an EM-based optimization transfer algorithm for direct penalized likelihood reconstruction that has faster convergence rate, especially at situations with low background events (randoms and scatters).

Direct reconstruction of parametric images for linear kinetic models has also been developed due to its computational efficiency. In 1997, Matthews et al 23 used the ML EM algorithm to estimate the parametric images of a set of linear temporal basis functions. Meikle et al 24 in 1998 presented a direct reconstruction for the spectral analysis model using the non-negative least squares method. Mathematically, direct reconstruction of linear parametric images is essentially the same as dynamic image reconstruction with overlapping temporal basis functions, such as B-splines 50. To avoid pre-defined basis functions, Reader et al also proposed a method to estimate the linear coefficients and temporal basis functions simultaneously 53. To obtain physiologically relevant kinetic parameters directly, Patlak model has been incorporated into direct reconstruction by defining two temporal basis functions as the blood input function and its integral 23. In 2007, Wang et al 33 presented a maximum a posteriori (MAP) reconstruction of the Patlak parameters using a preconditioned conjugate gradient (PCG) algorithm. Tsoumpas et al 34 presented an ML EM algorithm for direct Patlak reconstruction. Li and Leahy 35, and Zhu et al 36 developed direct reconstruction of the Patlak parameters from list-mode data. Tang et al 37 used anatomical prior information to improve the Patlak reconstruction. Rahmim et al used the AB-EM algorithm 70 to allow negative values in the direct estimation for a Patlak-like graphical analysis model 71. It has been observed that the strong correlation between the two temporal basis functions in the Patlak model slows down the convergence speed of the direct reconstruction. A nested EM algorithm was developed by Wang and Qi 38 to improve the convergence rate of direct linear parametric reconstruction.

1.4 Aim of This Paper

In this paper we will provide an overview of the recent progress in the development of direct reconstruction algorithms. We will describe technical details of different algorithms and analyze their properties. Methods for direct reconstruction of linear and nonlinear parametric images as well as for joint estimation of parametric images and input function will be included. We expect the information in this review will provide useful guidance to readers for choosing a proper direct reconstruction method.

2 General Formulation of Direct Reconstruction

2.1 PET Data Model

A dynamic PET scan is often divided into multiple consecutive time frames, with each frame containing coincidence events recorded from the start of the frame till the end of the frame. The image intensity at pixel in time frame m, , is then given by

where and denote the start and end times of frame , respectively, and is the decay constant of the radiotracer. is the tracer concentration in pixel at time and is determined by a linear or nonlinear kinetic model with the parameter vector . is the total number of kinetic parameters for each pixel.

Dynamic PET measurements can be well modeled as a collection of independent Poisson random variables 39-42,

where and are the indices of detector pairs and time frames, respectively, and and are the total numbers of the detector pairs and of the time frames, respectively. The expected projection is related to the dynamic image )}, through an affine transform,

where , the th element of the system matrix , is the probability of detecting an event originated in pixel by detector pair , and is the expectation of scattered and random events at detector pair in the th frame. is the total number of image pixels.

Let and be the projection vectors and the image vector in frame . The forward model for time frame can be rewritten in a matrix-vector product form:

Let us define a matrix of kinetic parameters, , in which each column represents an image of a kinetic parameter. We use to denote its ordered vector. Similarly we also define the following matrices:

and their ordered vectors are denoted by and respectively. Then the expectation of the whole dynamic PET data can be expressed either by

or

where is a identity matrix and denotes the Kronecker product. Equation (6) is the form for practical implementation and equation (7) is more suitable for algorithm analysis.

2.2 Maximum Likelihood Reconstruction

The goal of direct reconstruction is to estimate from dynamic PET measurements. Let denote a set of dynamic PET measurements. and are the matrix and the vector that are formed by in the same fashion as and , respectively. The log Poisson likelihood function of the dynamic PET data is

where a constant term is neglected and is defined by

Maximum likelihood (ML) reconstruction finds the solution by maximizing the log likelihood function,

where denotes the feasible set of the kinetic parameters (e.g. satisfying nonnegativity or box constraints). The resulting images from ML reconstruction at convergence are often very noisy because the tomography problem is ill-conditioned. In practice, the true ML solution is seldom being sought. Rather images are regularized either by early termination before convergence or using a penalty function or image prior to encourage spatial smoothness.

2.3 Penalized Likelihood Reconstruction

Penalized likelihood (PL) reconstruction (or equivalently maximum a posteriori, MAP) regularizes the solution by incorporating a roughness penalty in the objective function to encourage spatially smooth images. PL reconstruction finds the parametric image that maximizes the penalized likelihood function as

where is a smoothness penalty and is the regularization parameter that controls the tradeoff between the resolution and noise. If is too small, the reconstructed image approaches the ML estimate and becomes very noisy; if is too large, the reconstructed image becomes very smooth and useful information can be lost.

The smoothness penalty can be applied either on the kinetic parameters or on the dynamic image , depending on the application. A penalty applied on the dynamic images can be expressed by

where is a smoothness penalty given by

and is the potential function. denotes the neighborhood of pixel ; is the weighting factor equal to the inverse distance between pixels and . A typical neighborhood includes the eight nearest pixels in 2D and 26 nearest voxels in 3D.

The basic requirement of is that it is even and non-decreasing of . A common choice in PET image reconstruction is the quadratic function 42. The disadvantage of the quadratic regularization is that it can over-smooth edges and small objects when a large is used. To preserve edges, non-quadratic penalty functions can be used. Examples includes the absolute value function , the Lange function 47, and other variations 45, 46. Nonconvex penalty functions can also be used to even enhance edges, but they are much less popular because the resulting objective function may have multiple local optima.

3 Reconstruction of Parametric Images for Linear Kinetic Models

3.1 Linear Kinetic Models

One way to represent time activity curves is to use a set of linear temporal basis functions. The tracer concentration at time t can be described by,

where is the total number of basis functions, is the th temporal basis function and is the coefficient to be estimated.

The temporal basis functions can be divided into two categories. The first category primarily focuses on efficient representation of time activity curves. Examples include B-splines 48-50 and wavelets 51, as well as those obtained from principal component analysis 52 or adaptively estimated from PET data 53. One advantage of these basis functions is that they can represent a wide variety of time activity curves. A disadvantage is that the associated coefficients are not directly related to the kinetic parameters of physiological interest. Kinetic modeling is often required to estimate kinetic parameters from the TACs after reconstruction. In comparison, the second category to be described below provides linear coefficients that are directly related to the kinetic parameters of interest.

3.1.1 Spectral Analysis

In spectral analysis the basis functions are a set of pre-determined exponential functions convolved with the blood input function 54, 55,

where denotes the rate constant of the th spectrum and '' represents the convolution operator. The exponential spectral bases are consistent with the compartmental models and the distribution volume (DV) of a reversible tracer can be directly computed from the spectral coefficients by

If the blood input function in (16) is replaced by a reference region TAC , the calculated quantity from the spectral coefficients then becomes the distribution volume ratio (DVR), which is related to binding potential (BP), a major parameter of interest in neuroreceptor studies, by

3.1.2 Patlak Model

The Patlak graphical method 56 is a linear technique which has been widely used in dynamic PET data analysis. For a tracer with an irreversible compartment, the time activity curve satisfies the following linear relationship at the steady state:

where is the time for the tracer to reach steady state. The Patlak slope represents the overall influx rate of the tracer into the irreversible compartment and has found applications in many disease studies. For example, is proportional to the glucose metabolic rate in FDG scans. In some applications, the plasma input function can also be replaced by a reference region input function .

Equation (19) can be rewritten into the following linear model:

where and the two basis functions are

3.1.3 Logan Plot and Relative Equilibrium Plot

For reversible tracers, the Logan plot can be used to estimate distribution volume or binding potential 72. The standard Logan plot equation is

However, the tissue activity curve is involved nonlinearly in the Logan plot (and also in the multilinear model by Ichise et al 73), which makes direct reconstruction more difficult. Considering the fact that the ratio between the tissue time activity and plasma concentration remains constant at relative equilibrium, Zhou et al 74 proposed replacing in the denominators by . The model is referred to as the relative equilibrium (RE) model. After a slight rearrangement, we can write the RE model in the following form:

where the cumulative time activity is expressed as a linear combination of the same two basis functions as those in the Patlak plot.

To use the RE model, all the time frames should start from so that the reconstructed dynamic images represent the cumulative time activity . One consequence is that the projection data in different frames are no longer independent, so the likelihood function in (9) needs some modification to obtain a true ML estimate. In addition, the intercept in Eq. (24) is usually negative, which makes the classic EM algorithm not applicable to the resulting direct reconstruction.

Taking the derivative on both sides of (24), Zhou et al 75 presented an alternative form of the RE model,

where is the first-order derivative of . One advantage of the modified RE model in (25) is that it does not involve the cumulative time activity. In addition, since both and are negative, we can rewrite (25) in positive terms

with and the two basis functions

The model equation (26) is now suitable for direct reconstruction using the classic EM algorithm because the coefficient and the basis function are both nonnegative when . Note that can also be replaced by a reference region TAC to calculate DVR.

3.2 Forward Projection

With a linear kinetic model in (15), the dynamic PET images can be expressed in a matrix-vector product

or equivalently,

where the th element of the temporal basis matrix is given by

Substituting the above equations into Eq. (6), we get the following forward model for the linear kinetic model:

or equivalently,

3.3 Direct Reconstruction Using a Single System Matrix

By treating as a single system matrix , any existing algorithms for static PET reconstruction can be used for the direct estimation of nonnegative linear parametric images. For example, applying the ML EM algorithm 61-63, we get the following update equation for direct linear parametric image reconstruction 23:

where the superscript denotes the th iteration and

Similar to static PET reconstruction, the convergence of the parametric EM algorithm can be very slow. Preconditioned conjugate gradient (PCG) algorithm and other accelerated algorithms 57, 58 can be used to achieve a faster convergence. One example of direct application of PCG to Patlak reconstruction can be found in 33. When negative values are present in the parametric images, such as those in the original relative equilibrium model (24), the AB-EM algorithm can be used 71.

One disadvantage of direct reconstruction using a single system matrix is that the convergence can be extremely slow when the temporal basis functions are highly correlated as in the spectral analysis and Patlak model 12. This is because the combination of the temporal correlation and high spatial dimension results in an ill-posed problem that is much worse than the static PET reconstruction. One way to solve this problem is to decouple the spatial image update and temporal parameter estimation at each iteration using the nested EM algorithm described below.

3.4 Nested EM Algorithm

The nested EM algorithm 38 at iteration first calculates an intermediate dynamic image by

and then updates the kinetic parameter estimate by another EM-like equation

where , , and is the sub-iteration number in iteration . Each full iteration of the nested EM algorithm consists of one iteration of EM-like emission image update (37) and multiple iterations of kinetic parameter estimation (38).

The nested EM algorithm was derived by using the optimization transfer principle 60. A surrogate function is constructed at each iteration,

where is given by (37). This surrogate function satisfies

Then the maximization of the original likelihood function is transferred into the maximization of the surrogate function that can be solved pixel-by-pixel,

Applying the EM algorithm to (42) results in the update equation in (38). The property of the surrogate function guarantees

and converges to the global solution.

It is easy to verify that the traditional EM algorithm in (34) is a special case of the nested EM algorithm with 38. By running multiple iterations of (38) with , the nested EM algorithm can substantially accelerate the convergence rate of the direct reconstruction of linear parametric images without affecting the overall computational time as the size of matrix is much smaller than that of the system matrix . This is clearly demonstrated by a toy example shown in Fig. 1 38. Starting from the same initial image, the nested EM takes 6 iterations to converge to the true solution, while the traditional EM requires more than 60 iterations. The nested EM algorithm can be further accelerated by considering the nested EM algorithm as an implicit preconditioner and using conjugate directions 38.

Figure 1

Isocontours of the likelihood function of a toy problem and the trajectories of the iterates of the traditional EM and the nested EM . The nested EM takes 6 iterations to converge to the final solution, while the traditional EM requires more than 60 iterations. Reprinted from 38 with permission.

4 Reconstruction of Parametric Images for Compartment Models

4.1 Compartment Modeling

While linear models are advantageous for computational efficiency, nonlinear kinetic models based on compartment modeling are more related to the well-developed biochemical kinetics 5. Under a compartment model, the total tracer concentration in tissue is

where is the fractional volume of blood, is the whole blood concentration, and represents the concentration of the th tissue compartment. The compartment concentrations are related with each other by a set of ordinary differential equations 5, 7,

where , and are the kinetic parameter matrices that are formed by the rate constants in , and denotes the system input.

Taking the Laplace transform of (45), we get

where and denote the Laplace transforms of and , respectively. Then solution of the differential equation (45) can be obtained by

where denotes the inverse Laplace transform.

For the commonly used two-tissue compartment model, we have

and

with , where are the tracer rate constants, and are the concentrations in the free and bound compartments, and is the tracer concentration in plasma. The analytical solution of is given by

where with , and “*” denotes the convolution operator.

4.2 Model-dependent Algorithms

4.2.1 PICD Algorithm

The parametric iterative coordinate descent (PICD) algorithm 25 was the first direct reconstruction algorithm implemented for a large-scale parametric reconstruction with a compartment model. It first transforms the rate constants in the compartment model into a set of auxiliary parameters and then estimates the auxiliary parameters directly from the sinogram data. For example, the TAC in a two-tissue-compartment model is rewritten as

where the auxiliary parameters are . Once the auxiliary parameters are estimated, the kinetic rate constants can be calculated by

To estimate the auxiliary parameters, a coordinate descent (CD) algorithm is used. The original log-likelihood function is approximated by its second-order Taylor expansion. At each iteration, the auxiliary parameters are updated by

where , and denote the gradient and Hessian of the negative log-likelihood function with respect to at and are given by

The first-order derivative and second-order derivative are respectively given by

The regularization term equals to with fixed at for all pixels other than pixel .

The PICD algorithm uses an iterative gradient descent algorithm to update the linear parameters in and an iterative golden section search to estimate the nonlinear components. To accelerate the convergence rate and reduce computational cost, the PICD algorithm updates the linear parameters more frequently than the nonlinear parameters.

The PICD algorithm is well suited for one- and two-tissue-compartment models. For kinetic models with more than two tissue compartments, however, the parameter transformation becomes too complicated.

4.2.2 PMOLAR-1T Algorithm

The PMOLAR-1T algorithm 30 is an EM algorithm for direct reconstruction of parametric images for the one-tissue compartment model. It is derived by introducing a new set of complete data that simultaneously decouples the pixel correlation as well as the temporal convolution. Here we introduce it using the optimization transfer framework.

First we note that the surrogate function in (40) is applicable to any kinetic models. For the one-tissue-compartment model, the tissue time activity curve can be described by

and we have

where . Here we approximate the convolution integral using the summation over uniformly sampled time points with being the time interval.

Applying the EM surrogate function to (40) one more time to decouple the temporal correlation, we can find the maximum of (40) iteratively by

where and is given by

For the one-tissue compartment model with , the optimization in (63) has the following analytical solution:

where is the inverse function of the function 30

and can be calculated by a look-up table.

The original PMOLAR-1T algorithm 30 only uses one subiteration with . Based on our experience with the nest EM algorithm for linear models, we expect that the convergence rate of the PMOLAR-1T algorithm can be accelerated by using more than one subiterations .

4.3 Model-independent Algorithms

4.3.1 Iterative NLS Algorithms

Reader et al 26 proposed an iterative nonlinear least squares (NLS) algorithm for direct parametric image reconstruction. It iteratively applies the NLS kinetic fitting following an EM image update:

where is the dynamic image given in (37). The weighting factor has to be chosen empirically. In 26 a uniform weight was used. Based on the fact that (40) resembles a Poisson log-likelihood function, Matthews et al 27 proposed the following nonuniform weighting factor

which substantially improved performance over the uniform weight.

However, neither the uniform nor the nonuniform weights provides any theoretical guarantee of convergence because the quadratic function in (68) is not a proper surrogate function. It has been observed that the uniform weight can result in non-monotonic changes in the likelihood function, although the results with the nonuniform weighting are reasonably good.

4.3.2 OT-SP Algorithm

To obtain a guaranteed convergence while maintaining the simplicity of the iterative NLS algorithm, Wang and Qi 29 proposed a quadratic surrogate function for the penalized likelihood function using optimization transfer (OT). The surrogate function of the log-likelihood is derived based on the fact that the second-order derivative of is a non-decreasing function of and is bounded when . The potential function in (14) also has a bounded, non-increasing second-order derivative of . Omitting the constants that are independent of , the surrogate function of the penalized likelihood is given by

where the optimum curvature is set to 43, 44

with and is an intermediate dynamic sinogram at iteration ,

The second term in (70) is the surrogate function for the penalty function and is given by

where is the optimum curvature of the penalty function

For example for the quadratic penalty and for the nonquadratic Lange penalty.

With , the surrogate satisfies the following condition:

Hence it minorizes the original objective function .

The optimization of the quadratic surrogate function can be solved pixel-by-pixel in a coordinate descent (CD) fashion29. However, the OT-CD algorithm is difficult to parallelize because of the sequential update.

For simultaneous update, Wang and Qi further developed an OT-SP (optimization transfer using separable paraboloids) algorithm 29 that uses separable paraboloids to decouple the correlations between pixels in at each iteration. The overall surrogate function that minorizes the original objective function at iteration n is

where the weighting factor in the separable least squares is

with . The intermediate dynamic image is given by

where is the gradient of with respect to ,

evaluated at .

Because is separable for pixels, maximization of is reduced to a pixel-wise least squares fitting optimization:

which can be solved by any existing nonlinear least squares methods (e.g. the Levenberg-Marquardt algorithm).

In summary, each iteration of the OT-SP algorithm consists of two steps: an image update in (78) and a NLS fitting in (80). It has the simplicity of the iterative NLS algorithm, but guarantees a monotonic convergence with . Compared to the PICD and OT-CD algorithms, the OT-SP algorithm can be easily parallelized. A disadvantage of the OT-SP is that it requires the background to be positive for the quadratic surrogate function to work, so the convergence rate can be slow when the level of background events is low.

4.3.3 OT-EM Algorithm

The OT-EM (optimization transfer using the EM surrogate) was derived for direct reconstruction under low background levels 31, 32. It combines the EM surrogate function in (40) for the log-likelihood function and a separable surrogate function for the penalty function. The overall surrogate function is

The optimization of the penalized likelihood function with respect to is transferred into the maximization of the surrogate function , which can be solved by the following small-scale pixel-wise optimization:

A modified Levenberg-Marquardt algorithm was developed to solve this pixel-wise penalized likelihood fitting in 32.

The OT-EM guarantees a monotonic increase in the penalized likelihood following the optimization transfer properties of . It can be much faster than the OT-SP algorithm when the level of background events (scatters and randoms) is low.

5 Joint Reconstruction of Parametric Images and Input Function

In parametric image reconstruction, the blood input function is required to be known a priori. One standard method of measuring blood input is arterial blood sampling, which is invasive and technically challenging. An alternative to arterial sampling is to derive the input function from a blood region or a reference region in a reconstructed image 65. When the region used for extraction of the input function is of small size, the image-derived input function (IDIF) can be less accurate due to partial volume and other effects. Joint estimation of the input function and parametric images can potentially improve the accuracy of the IDIF as well as the parametric images by fitting the input function and parametric image to the data simultaneously 64, 65.

5.1 General Formulation

Given a set of kinetic parameters, the time activity curve is a linear function of the input function u(t), which can be either the blood input function or a reference tissue function . Let denote the vector of parameters that defines the input function (e.g. the activity values at a set of time points). Then the joint estimation can be formulated as

where denotes the feasible set of and the penalized likelihood function defined in (12) is rewritten here as an explicit function of .

The maximization in (83) can be solved by an alternating update algorithm (e.g. in 66, 69)

5.2 Joint Estimation without an Input Region

One method for joint estimation without an input region was presented by Reader et al 67. The algorithm was developed for the linear spectral analysis model. It uses the EM algorithm to find the ML estimation of the input function and the linear coefficients from dynamic projection data. Here we introduce it using the optimization transfer principle and describe how to extend it to nonlinear kinetic models.

We use the EM surrogate function in (40) for the log-likelihood function . Here we rewrite it as an explicit function of the input parameters ,

where is given by (37) and is the same as but as an explicit function of . Since is separable for pixels when v is given, the optimization (84) is transferred to

For a linear kinetic model, both above steps can be solved by the EM algorithm with multiple sub-iterations. When only one sub-iteration is used, it is equivalent to the algorithm in 67.

For compartment models, the optimization in (87) can be solved by the modified Levenberg-Marquardt algorithm given in 32. The update of in (88) can still be solved by the EM algorithm.

5.3 Joint Estimation with an Input Region

One issue with the joint estimation without an input region is that there is an unknown scaling factor on and the linear coefficients in that cannot be determined. The problem can be solved if there is a region in the image that contains purely the input function. Let denote such a region. Then the image intensity at pixel in time frame is

The joint estimation then estimates the kinetic parameters at pixels outside and the time activity curve inside simultaneously.

Using the same paraboloidal surrogate functions for the OT-SP algorithm, the following algorithm can be obtained 68:

where and are updated in an alternating fashion. for is updated by the Levenberg-Marquardt algorithm. The estimation of in ( (90) is a linear least squares problem that can be solved easily (e.g. by QR decomposition).

6 Examples

Here we show some simulation results to demonstrate the advantage of direct reconstruction over indirect reconstruction. We simulated a 60-minute F-FDG brain scan using a brain phantom that consisted of gray matter, white matter and a small tumor inside the white matter. The time activity curve of each region was generated using a two-tissue compartment model and an analytical blood input function. Figure 2 shows the images reconstructed by an indirect method and the OT-EM algorithm and Figure 3 shows the corresponding images. Clearly the direct reconstruction results have much less noise than the indirect reconstruction results while maintaining the same spatial resolution. The improvement is similar to those observed in other studies for both linear and compartment models 25, 29, 33.

Figure 2

The true and reconstructed images by the indirect and direct algorithms .

Figure 3

True and reconstructed images by the indirect and direct algorithms .

7 Conclusion

We have provided an overview of direct estimation algorithms for parametric image reconstruction from dynamic PET data. Because of the combination of tomographic reconstruction and kinetic modeling, direct reconstruction algorithms are more complicated than static image reconstruction. Different algorithms have been developed to balance the tradeoff between simplicity in implementation and fast convergence. We hope that the information in this paper can provide guidance to readers for choosing the proper direct reconstruction algorithm in real applications. In some situations, the choice is clear. For example, the nested EM algorithm should be preferred over the traditional EM algorithm for linear parametric image reconstruction because the former converges much faster while preserving the same simplicity in implementation. In other situations, the choice may be more dependent on the user. We note that as long as the algorithm guarantees global convergence, a simple algorithm can still reach the final solution, albeit with a large number of iterations. We expect that direct parametric image reconstruction will find more and more applications in dynamic PET studies and more novel algorithms will be developed in the future.