Kinetic Compressive Sensing

Kinetic Compressive Sensing

Michele Scipioni,  Maria F. Santarelli, Luigi Landini, Ciprian Catana, Douglas N. Greve, Julie C. Price, and Stefano Pedemonte M. Scipioni and L. Landini are with the Department of Information Engineering, University of Pisa, Pisa, Italy (email: Santarelli and L. Landini are with the Institute of Clinical Physiology, CNR, Pisa, ItalyM. Scipioni, C. Catana, D.N. Greve, J.C. Price, and S. Pedemonte are with the Athinoula A. Martinos Center for Biomedical Imaging, Boston, MA, USA.C. Catana, D.N. Greve, J.C. Price, and S. Pedemonte are with Harvard Medical School, Boston, MA, USA.S. Pedemonte is with the Center for Clinical Data Science, MGH, Boston, MA, USA

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 post-reconstruction 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 (OS-MAP-OSL), and kinetic parameters (MAP-LM). 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 post-reconstruction 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 KCS-regularized direct method can produce spatially coherent images and parametric maps, with lower spatial noise and better tissue contrast. A GPU-based open source implementation of the algorithm is provided.

parametric images, PET, compartmental models, compressive sensing, hierarchical Bayesian model, sparsity, Markov Random Field, FDG, GPU
publicationid: pubid: 978-1-5386-2282-7/17/$31.00 © 2017 IEEE

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, ligand-protein binding interactions, response to pharmacological challenges, and gene expression.

There are two approaches for deriving kinetic parameters: region-of-interest (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 maximum-a-posteriori version of the standard Levemberg-Marquardt 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 Ordered-Subsets Maximum-A-Posteriori One-Step-Late reconstruction algorithm with Kinetics-derived Prior (OS-MAP-OSL-KP) regularization term. In this article, we describe experiments with both simulated data and real [18F]FDG PET data. We also provide a GPU-based 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 back-projection 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.

Ii-a Hierarchical Bayesian model

The overall model is composed of three parts:

  1. the model of the acquisition system;

  2. the kinetic model;

  3. 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 multi-compartmental model. The one-to-one 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 sparsity-inducing prior is introduced at this level as a Markov Random Field (MRF) with an L1-norm cost function.

Fig. 1: Hierarchical Bayesian Model. From left to right: the kinetic model treats each voxel as a realization of a hidden dynamic process; the sparsity-inducing prior acts voxel-wise on the parametric maps as an MRF with L1-norm cost function; the model of the acquisition process links image domain and measured projections.

Ii-B Image reconstruction with kinetics-derived 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 log-posterior 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 log-posterior probability:


where is the conventional Poisson likelihood, and is a positive scaling factor.

As shown by Green [OSL], this problem can be solved using the one-step-late (OSL) approach for an iterative update of the MAP estimate:


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 .


where represent the set of kinetic parameters estimated for voxel , and is the modeled TAC for voxel at time point .

To implement the aforementioned OS-MAP-OSL-KP reconstruction algorithm (2)-(3), we used the in-house developed [occiput] software, which uses GPU parallel computation for the operations of projection and back-projection to speed up the reconstruction process.

Ii-C Update of parametric maps

To update the prior in (3) after each iteration of OS-MAP-OSL-KP (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:

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

  2. 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 single-instruction-multiple-data (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 Maximum-a-Posteriori Levemberg-Marquardt (MAP-LM) algorithm for nonlinear least square optimization, based on cuBLAS batched matrix operations. An additional speed-up 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:


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 MAP-LM 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:




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 non-interacting line-site 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.

Ii-D Summary of the algorithm

Fig. 2: Algorithm flow chart: we alternate between the image reconstruction via the OS-MAP-OSL-KP algorithm, and the GPU parallel implementation of parametric mapping via MAP-LM.

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.

  1. given an updated estimate of the kinetic maps, we use it to produce a new estimate of the image time series via OS-MAP-OSL-KP algorithm;

  2. given the new estimated activity time series, we updated the parametric maps via MAP-LM optimization.

Iii Simulation

Iii-a Simulation setup

To assess the performance of the parallel GPU implementation, and the effect of the sparsity-inducing 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 2-tissues 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 non-uniform 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:

Fig. 3: Simulated geometrical phantom
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
TABLE I: Kinetic parameters used in simulation

Iii-B Simulation results

We compared the results of three different methods:

  1. indirect reconstruction (INDIRECT): OSEM and post-fitting of kinetic model;

  2. direct reconstruction (DIRECT): OS-MAP-OSL-KP an presented in (2), but without using a spatial regularization in (4);

  3. direct reconstruction with kinetic compressive sensing (KCS DIRECT): like the previous one but with the addition of the prior term in (4) defined as in (6).

Fig. 4: Simulation study: late time frame reconstruction performed using (a) OSEM reconstruction, (b) INDIRECT post-reconstruction fitting, (c) DIRECT reconstruction, and (d) KCS DIRECT reconstruction.

Fig. 4 shows the reconstruction of a late time frame. It is easy to recognize a first reduction in voxel-by-voxel 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 post-reconstruction 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.

Fig. 5: Bias/variance plot of the estimated parametric maps, for the different methods.

Iv Human Data

Iv-a Dataset

The same three methods tested by simulation (conventional indirect, direct, and KCS) were applied to a [18F]FDG brain PET scan. The data were acquired on a Siemens mMR PET-MR 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 2-tissue irreversible compartment model adopted for the simulation.

Iv-B 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 (c-f) is able to produce spatially coherent images with lower noise, above all with respect to the indirect estimate (a-d), and better tissue contrast, also compared to a direct method without the spatial KCS regularization (b-e).

It is important to take into account how in this work we decided to test the algorithm with a bi-compartmental 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 macro-parameter , derived as , also to make it easier to compare the results with other similar methods using linear models (i.e. Patlak).

Fig. 6: Human data results. Top row: the last reconstructed time frame generated by standard OSEM reconstruction (a), OS-MAP-OSL-KP direct approach, which uses the kinetic model as a guidance during the reconstruction (b), and the KCS-regularized version of direct reconstruction (c). Bottom row: the parametric map produced by each of the three tested methods.

V Conclusion

The simulation study demonstrated that the proposed method of introducing a sparsity-inducing prior in a direct reconstruction framework can help in producing high-quality images and parametric maps, that are both amenable for display and also quantitatively more accurate than what a post-reconstruction 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 OS-MAP-OSL-KP (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.


The Tesla K20 used for this research was donated by the NVIDIA Corporation.


Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description