Kinetic Compressive Sensing
Abstract
Parametric images provide insight into the spatial distribution of physiological parameters, but they are often extremely noisy, due to low SNR of tomographic data. Direct estimation from projections allows accurate noise modeling, improving the results of postreconstruction fitting. We propose a method, which we name kinetic compressive sensing (KCS), based on a hierarchical Bayesian model and on a novel reconstruction algorithm, that encodes sparsity of kinetic parameters. Parametric maps are reconstructed by maximizing the joint probability, with an Iterated Conditional Modes (ICM) approach, alternating the optimization of activity time series (OSMAPOSL), and kinetic parameters (MAPLM). We evaluated the proposed algorithm on a simulated dynamic phantom: a bias/variance study confirmed how direct estimates can improve the quality of parametric maps over a postreconstruction fitting, and showed how the novel sparsity prior can further reduce their variance, without affecting bias. Real FDG PET human brain data (Siemens mMR, 40min) images were also processed. Results enforced how the proposed KCSregularized direct method can produce spatially coherent images and parametric maps, with lower spatial noise and better tissue contrast. A GPUbased open source implementation of the algorithm is provided.
I Introduction
Kinetic analysis of medical imaging data is a standard approach to characterize the body’s physiological processes, allowing in vivo studies of metabolic patterns, perfusion, blood flow, ligandprotein binding interactions, response to pharmacological challenges, and gene expression.
There are two approaches for deriving kinetic parameters: regionofinterest (ROI) modeling and parametric imaging. The ROI based approach fits a kinetic model to the average time activity curve (TAC) of a selected ROI: it is easy to implement and has low computational cost. In contrast, parametric imaging estimates kinetic parameters for every voxel, providing richer kinetic data and enabling the visualization of the spatial distribution of physiological parameters. However, parametric images are often extremely noisy, as a reflection of the noisy measurement of tomographic data.
In conventional indirect kinetic analysis, the parametric maps are estimated by reconstructing a time series of tomographic images, then fitting a compartmental kinetic model. Direct estimation of parametric images [reader2014] from raw projection data, which has become of interest with the availability of more powerful computers and of GPUs, allows accurate noise modeling and has been shown to offer better image quality than conventional indirect methods. Despite the improvement introduced by direct kinetic modeling, parametric images remain still extremely noisy.
In this work, we use a novel approach for integrating information derived from kinetic modeling into the dynamic reconstruction in the form of a regularization prior, in order to reduce the effect of fitting errors on image reconstruction and convergence. Moreover, in order to extract from the noisy PET measurements spatially coherent parametric images, we propose a maximumaposteriori version of the standard LevembergMarquardt nonlinear least square optimization, where we use an Huber prior to enforce an hypothesis of sparsity of the spatial derivatives of the kinetic maps.
The method, which we name kinetic compressive sensing (KCS), is based on a hierarchical Bayesian model, and it can be described as an OrderedSubsets MaximumAPosteriori OneStepLate reconstruction algorithm with Kineticsderived Prior (OSMAPOSLKP) regularization term. In this article, we describe experiments with both simulated data and real [18F]FDG PET data. We also provide a GPUbased open source implementation of the reconstruction algorithm [gpuKMfit].
Ii Methods
The method is based on a hierarchical Bayesian model, from which we derived an optimization algorithm for the direct estimation of the kinetic parameters, which includes sparsity constraints. We provide an efficient GPU implementation of the projection and backprojection operations, as well as of the parametric maps estimation, based on the analytic evaluation of the gradient of the log joint posterior of the Bayesian model.
Iia Hierarchical Bayesian model
The overall model is composed of three parts:

the model of the acquisition system;

the kinetic model;

a model to promote sparsity of the kinetic maps’ derivatives.
Reading the graph in fig.1 from right to left, we first find the level of projection domain. This is linked to the activity image domain by the model of the acquisition system: it consists of the ordinary Poisson model, incorporating all effects of attenuation, scatter and randoms. Since we are working with dynamic 4D datasets, each voxel has a relevant time activity curve (TAC): the kinetic model encodes the assumption that the voxel intensities are noisy realizations of a hidden dynamic process, modeled in this case using a multicompartmental model. The onetoone connection between elements of the activity domain and elements of the kinetic parameters domain is the result of model fitting: each TAC is linked to a set of parameters, and so each layer of the kinetic parameters domain can be seen as a parametric map. The sparsityinducing prior is introduced at this level as a Markov Random Field (MRF) with an L1norm cost function.
IiB Image reconstruction with kineticsderived prior
In this work we chose to incorporate time information derived from kinetic modeling as a regularization term to help reducing fluctuations and noise in successive frames of our time series, which are common if we reconstruct each one of them independently.
In the general MAP framework, we want to maximize the logposterior distribution of our measured data , with respect to the image estimate , and subject to a regularization term . This usually means looking for that maximizes the logposterior probability:
(1) 
where is the conventional Poisson likelihood, and is a positive scaling factor.
As shown by Green [OSL], this problem can be solved using the onesteplate (OSL) approach for an iterative update of the MAP estimate:
(2) 
where is the voxel index, is the time frame number, a single bin of the projection data measured at time is represented as , and represents an element of the system matrix, modeling the contribution of voxel to projection bin . The OSL strategy involves calculation of the derivative of , in order to update the image vector from the old estimate.
Regarding the choice of the potential function , our proposal is specifically designed for 4D dynamic acquisitions because the information we want to introduce as a prior is aimed to suppress the temporal fluctuations (i.e. noise) in each voxel’s TAC due to the independent reconstruction of each time frame. This issue has been addressed by other direct reconstruction methods by forcing the new image update to start from the output of the kinetic modeling step computed on the previous 4D estimate. On the contrary, we thought of introducing a time constraint through a Gibbs prior as a gentler way to constrain the reconstruction, while still being as close as possible to the measurement.
To do this, we chose our potential function so to enforce similarity between correspondent time points between the reconstructed time series and the model simulation computed at the previous step. As for standard OSL approach, in order to do this we have to approximate as .
(3) 
where represent the set of kinetic parameters estimated for voxel , and is the modeled TAC for voxel at time point .
IiC Update of parametric maps
To update the prior in (3) after each iteration of OSMAPOSLKP (2), we need to voxelwise fit the kinetic model to the intermediate image update . This step is usual for most common nested direct reconstruction algorithms. In this work we focused on tackling the two most known issues of the standard approach:

the execution time needed to loop through each voxel’s TAC to fit the model;

the low SNR of the kinetic maps due to different convergence points for each voxel, when we treat each one of them as an independent optimization problem.
Fast parallel GPU implementation
To address the issue of long execution time, we chose to move the problem of fitting a non linear compartment model to every voxel in our image from CPU to GPU. GPUs are singleinstructionmultipledata (SIMD) devices. This means they are very efficient at performing the same algorithmic steps, in parallel, on large sets of data. In this work, we exploited this property of GPU computing to develop a set of CUDA kernels to evaluate the time integrals of the kinetic model and its derivatives with respect to the unknown parameters.
We also developed a CUDA implementation of a MaximumaPosteriori LevembergMarquardt (MAPLM) algorithm for nonlinear least square optimization, based on cuBLAS batched matrix operations. An additional speedup was obtained by avoiding numeric integration, deriving the analytic expression of the model equations and their derivatives, following an idea proposed in [miccai2015].
In short, the update of parametric maps follows the following equation:
(4) 
where, for the current voxel assigned to a specific GPU thread, is the estimated parameter set, is Jacobian matrix, and is the reconstructed TAC. is a potential function used as spatial regularization term for MAP fitting of the model to the current voxel, whose derivation is described in the next paragraph.
Derivation of the regularization term for MAPLM optimization
The modeling step in (4) is an iterative process of its own. Thanks to the parallel GPU implementation, we are now able to synchronize the update of the fitting for all the voxel so that, after each iteration of (4) we have an intermediate estimate of the parametric maps we can use to derive the prior term .
In an attempt to reduce parametric maps’ spatial noise, and to guide the convergence of voxels with similar kinetic to the same minima (i.e. parameter sets), we decided to impose an a priori assumption to the LM optimization based on a roughness penalty. We chose a Huber function:
(5) 
(6) 
where
(7) 
The idea is to rely on information derived from differences in voxel values computed on a local neighborhood, with a local penalty that penalizes large differences (i.e. borders) less than a quadratic penalty.
The actual filtering kernel used in (4) is derived from the Markov Random Field (MRF) prior conditioned on a noninteracting linesite model, presented in fig.1, which uses small cliques of neighboring voxels. The convolution between Huber kernel and current estimate of parametric maps, , is computed in a highly parallel fashion thanks to the development of a dedicated CUDA kernel running on GPU.
IiD Summary of the algorithm
The parametric maps are reconstructed by maximizing the joint probability of all unknown parameters in the graph in Fig.1 by means of the iterated conditional modes (ICM) algorithm [ICM]. As shown in fig.2, the algorithm consists of alternating the optimization of the activity time series and of the kinetic parameters.

given an updated estimate of the kinetic maps, we use it to produce a new estimate of the image time series via OSMAPOSLKP algorithm;

given the new estimated activity time series, we updated the parametric maps via MAPLM optimization.
Iii Simulation
Iiia Simulation setup
To assess the performance of the parallel GPU implementation, and the effect of the sparsityinducing regularization of the kinetic map estimation, we realized a Monte Carlo (MC) simulation with 100 noisy realizations using the phantom in Fig.3.
The kinetic behavior of the three main regions has been simulated using a 2tissues irreversible compartment model, using the parameters shown in table I. The square region in the center is chosen as to emulate a blood input region.
The nonuniform time scheme used for the simulation covers 40 minutes of acquisition with 24 time frames, and it matches the protocol used for the following human data study:
(8) 
Region  (%)  ()  ()  () 

Blood  1.0  0.0  0.0  0.0 
Inner  0.13  0.75  0.35  0.031 
Middle  0.116  0.62  0.30  0.026 
Outer  0.0985  0.51  0.27  0.018 
IiiB Simulation results
We compared the results of three different methods:

indirect reconstruction (INDIRECT): OSEM and postfitting of kinetic model;
Fig. 4 shows the reconstruction of a late time frame. It is easy to recognize a first reduction in voxelbyvoxel variance when the kinetic model is used to regularize the reconstruction (c), which is further reduced when the sparsity assumption of the spatial derivatives of the parametric maps is enforced (d).
The bias/variance plot in fig. 5 shows how a direct approach improves the quality of the estimate of parametric maps, with respect to the results provided by a standard indirect postreconstruction fitting, but also how the novel sparsity constraint is able to further reduce the variance of the produced parametric maps, without affecting, if not decreasing, the bias, also for a high number of reconstruction iterations.
Iv Human Data
Iva Dataset
The same three methods tested by simulation (conventional indirect, direct, and KCS) were applied to a [^{18}F]FDG brain PET scan. The data were acquired on a Siemens mMR PETMR scanner, recording a 40 minute 3D listmode dataset.
The raw data were then rebinned into sinograms, using the same time scheme described in (8) for the simulation.
Given the kinetic behavior of the radiotracer under study in brain, during the reconstruction and modeling phase we used the same 2tissue irreversible compartment model adopted for the simulation.
IvB Results
Fig.6 shows a comparison of the results obtained using the three different methods, in terms of images reconstruction (top row) and (net uptake rate) map estimate (bottom row). Both perspectives convey the same message: the proposed method (cf) is able to produce spatially coherent images with lower noise, above all with respect to the indirect estimate (ad), and better tissue contrast, also compared to a direct method without the spatial KCS regularization (be).
It is important to take into account how in this work we decided to test the algorithm with a bicompartmental model to describe tracer uptake in tissue. This means that the final output, in terms of parametric maps, is a set of maps, one for each micro parameter of the model ().
We chose here to show the map of the macroparameter , derived as , also to make it easier to compare the results with other similar methods using linear models (i.e. Patlak).
V Conclusion
The simulation study demonstrated that the proposed method of introducing a sparsityinducing prior in a direct reconstruction framework can help in producing highquality images and parametric maps, that are both amenable for display and also quantitatively more accurate than what a postreconstruction fitting and an unconstrained direct reconstruction can achieve, from a bias/variance point of view.
This method appears to be promising as a feasible approach to applying kinetic modeling to very large 4D clinical datasets with a reduced computational cost, thanks to the parallel GPU implementation based on the analytic expression of the kinetic model and its derivatives.
Future studies will focus on evaluating the performance of the proposed OSMAPOSLKP (2) reconstruction algorithm with respect to other known approaches for direct reconstruction [reader2014], with and without the spatial KCS prior presented in this work.
Moreover, the proposed KCS algorithm could be adapted to other applications, such as kinetic modeling of dynamic contrast enhanced (DCE) MRI and CT.
Acknowledgment
The Tesla K20 used for this research was donated by the NVIDIA Corporation.