Kinetic Compressive Sensing
Michele Scipioni Maria F. Santarelli Luigi Landini Ciprian Catana Douglas N. Greve Julie C. Price Stefano Pedemonte



























































978-1-5386-2282-7/17/$31.00  2017 IEEE












Kinetic Compressive Sensing
    Michele Scipioni, Student Member, IEEE,
        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: michele.scipioni@ing.unipi.it)
M.F. Santarelli and L. Landini are with the Institute of Clinical Physiology, CNR, Pisa, Italy
M. 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



    Accepted 23/03/2018
==========================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================






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

















§ 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 <cit.> 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 <cit.>. 







§ 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.



 §.§ 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.<ref> 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. 






 §.§ 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 s, with respect to the image estimate x, and subject to a regularization term U(x). This usually means looking for x̂ that maximizes the log-posterior probability:


    l(s|x) = logP(s|x)-β U(x) + const


where P(s|x) is the conventional Poisson likelihood, and β is a positive scaling factor.

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


    x_vm^new =  x_vm^old/∑_d p_vd - β∂ U(x)/∂ x_vm|_x_vm=x_vm^old∑_dp_vd s_dm/∑_vp_vd x_vm^old


where v is the voxel index, m is the time frame number, a single bin d of the projection data measured at time m is represented as s_dm, and p_vd represents an element of the system matrix, modeling the contribution of voxel v to projection bin d. The OSL strategy involves calculation of the derivative of U(x), in order to update the image vector x from the old estimate.

Regarding the choice of the potential function U(x), 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 U(x) so to enforce similarity between correspondent time points between the reconstructed time series x_v:^new and the model simulation F(ξ_v^old)_: computed at the previous step. As for standard OSL approach, in order to do this we have to approximate x_v:^new as x_v:^old.


    U(x) ∝[F(ξ_v^old)_m - x_vm^old]^2/2σ^2


where ξ_v^old represent the set of kinetic parameters estimated for voxel v, and F(ξ_v^old)_m is the modeled TAC for voxel v at time point m. 

To implement the aforementioned OS-MAP-OSL-KP reconstruction algorithm (<ref>)-(<ref>), we used the in-house developed Occiput.io <cit.> software, which uses GPU parallel computation for the operations of projection and back-projection to speed up the reconstruction process. 







 §.§ Update of parametric maps

To update the prior in (<ref>) after each iteration of OS-MAP-OSL-KP (<ref>), we need to voxelwise fit the kinetic model to the intermediate image update x_vm^new. 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 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 <cit.>. 

In short, the update of parametric maps follows the following equation:

    ξ_v^new = ξ_v^old + { (J_v^T J_v + λ I)^-1   [ J_v^T ( x_v - F(ξ_𝐯^𝐨𝐥𝐝)) + 
        - β∂ H(ξ_v)/∂ξ_v ] }


where, for the current voxel v assigned to a specific GPU thread, ξ_v is the estimated parameter set, J_v is Jacobian matrix, and x_v is the reconstructed TAC. H(ξ_v) 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 (<ref>) 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 (<ref>) we have an intermediate estimate of the parametric maps we can use to derive the prior term H(ξ_v).

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: 


    H(ξ) ∝
            t^2                               |t| ≤δ
    δ |t|-δ^2/2    |t| > δ



    ∂ H(ξ)/∂ξ = 
         
            -δ    t < -δ
    
            t              -δ≤ t ≤δ
    
            -δ    t > δ


where


    t = ξ_j-ξ_k    j,k ∈Ω



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 H(ξ) used in (<ref>) is derived from the Markov Random Field (MRF) prior conditioned on a non-interacting line-site model, presented in fig.<ref>, which uses small cliques of 3 × 3 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.









 §.§ Summary of the algorithm




The parametric maps are reconstructed by maximizing the joint probability of all unknown parameters in the graph in Fig.<ref> by means of the iterated conditional modes (ICM) algorithm <cit.>. As shown in fig.<ref>, 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 OS-MAP-OSL-KP algorithm;
    
  * given the new estimated activity time series, we updated the parametric maps via MAP-LM optimization. 










§ SIMULATION




 §.§ 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.<ref>. 

The kinetic behavior of the three main regions has been simulated using a 2-tissues irreversible compartment model, using the parameters shown in table <ref>. 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: 
    [12× 10s; 2× 30s; 4× 60s; 1× 120s; 4× 300s; 1× 600s]










 §.§ Simulation results

We compared the results of three different methods:

    
  * indirect reconstruction (INDIRECT): OSEM and post-fitting of kinetic model;
    
  * direct reconstruction (DIRECT): OS-MAP-OSL-KP an presented in (<ref>), but without using a spatial regularization in (<ref>);
    
  * direct reconstruction with kinetic compressive sensing (KCS DIRECT): like the previous one but with the addition of the prior term in (<ref>) defined as in (<ref>).




Fig. <ref> 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. <ref> 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.


























§ HUMAN DATA




 §.§ 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 (<ref>) 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.  



 §.§ Results


Fig.<ref> shows a comparison of the results obtained using the three different methods, in terms of images reconstruction (top row) and K_i (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 (K_1, k_2, k_3, f_v). 

We chose here to show the map of the macro-parameter K_i, derived as K_i = K_1 k_3 / (k_2+k_3), also to make it easier to compare the results with other similar methods using linear models (i.e. Patlak).
























§ 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 (<ref>) reconstruction algorithm with respect to other known approaches for direct reconstruction <cit.>, 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.


IEEEtran




















































































