Add weights to Patlak linear regression

documentation/release_5.0.htm

@@ -139,6 +139,12 @@ <h3>New functionality</h3>
 <li>
     Exposed <code>compute_total_ROI_values</code> and <code>ROIValues</code> via swig. Added python example <code>ROI_demo.py</code> to demonstrate usage.
 </li>
+<li>
+ <code>PatlakPlot::apply_linear_regression</code> now has the possibility to use weighted
+ linear regression, assuming that the TAC activity values (in units of <i>counts</i>)
+ Poisson distributed which is a reasonable approximation for OSEM images (although somewhat
+ dangerous for noisy data).
+</li>
 </ul>
 
 

examples/samples/PatlakPlot.par

@@ -9,6 +9,12 @@ blood data filename :=
 Time Shift := 0
 ; following settings are appropriate for STIR reconstructed images
 In total counts := 1
+; enable this for fully calibrated images from your scanner
+; (should remain zero for STIR images currently)
 In correct scale := 0
 
+; enable this if you want to use weighted linear regression, assuming
+; that each voxel and time frame is Poisson distributed
+Poisson distributed images := 0 ; defaults to 0
+
 end Patlak Plot Parameters:=

src/include/stir/modelling/PatlakPlot.h

@@ -2,6 +2,7 @@
 //
 /*
     Copyright (C) 2006 - 2011, Hammersmith Imanet Ltd
+    Copyright (C) 2021, University College London
     This file is part of STIR.
 
     SPDX-License-Identifier: Apache-2.0
@@ -14,7 +15,7 @@
   \brief Implementation of functions of class stir::PatlakPlot
 
   \author Charalampos Tsoumpas
-
+  \author Kris Thielemans
 */
 
 
@@ -43,6 +44,13 @@ START_NAMESPACE_STIR
     <i>Experimental and Graphical evaluation of blood-to-brain transfer constant from multiple-time uptake data: Generalizations,</i>
     J Cereb Blood Flow Metab 5: p. 584-90. 
 
+  The kinetic model is
+  \f[
+     C(t) =Ki*\int{C_p(t)\,dt}+V C_p(t)
+  \f]
+  with \f$C(t)\f$ the tissue TAC and \f$C_p(t)\f$ the (plasma) input function, with \f$Ki\f$ the slope
+  and \f$V\f$ the intercept. \f$t\f$ is a time-frame index (i.e. the activity/input function is
+  integrated over the time frame).
 
   \par Example .par file
   \verbatim
@@ -54,17 +62,27 @@ START_NAMESPACE_STIR
   blood data filename :=  blood_file.txt
   ; In seconds
   Time Shift := 0
-  In total counts := 1
+  ; work-around current lack of STIR unit info
+  In total counts := 0 ; defaults to 0, set to 1 for images in activity
+  ; enable this for fully calibrated images from your scanner
+  In correct scale := 0 ; defaults to 0
+
+  ; enable this if you want to use weighted linear regression, assuming
+  ; that each voxel and time frame is Poisson distributed
+  Poisson distributed images := 0 ; defaults to 0
 
   end Patlak Plot Parameters:=
   \endverbatim
 
   \warning
   - The dynamic images will be calibrated only if the calibration factor is given. 
-  - The [if_total_cnt] is set to true the Dynamic Image will have the total number of 
-    counts while if set to false it will have the total_number_of_counts/get_duration(frame_num).
-  - The dynamic images will always be in decaying counts.
-  - The plasma data is assumed to be in decaying counts.
+  - The \c if_total_cnt is set to \c true the Dynamic Image will be assumed to have the total number of
+    counts while if set to \c false it will have the \c total_number_of_counts/get_duration(frame_num).
+  - The dynamic images will always be assumed to be without decay correction.
+  - The plasma data is assumed to be without decay correction.
+
+  This class provides functionality for iterative estimation as well as "ordinary"
+  linear regression (see \c apply_linear_regression()).
 
   \todo Should be derived from LinearModels, but when non-linear models will be introduced, as well.  
 */
@@ -122,8 +140,16 @@ class PatlakPlot : public RegisteredParsingObject<PatlakPlot, KineticModel>
       get_dynamic_image_from_parametric_image(DynamicDiscretisedDensity & dyn_image,
 					      const ParametricVoxelsOnCartesianGrid & par_image) const;
 
-    //! This is the common method used to estimate the parametric images from the dynamic images.
-    /*! \todo There is currently no check if the time frame definitions from \a dyn_image are
+    //! estimate parametric images from dynamic images
+    /*! This performs Patlak Linear regression is applied to the data by minimising
+       \f[\sum_t w_t ( C(t)/C_p(t) - (Ki*\int{C_p(t)\,dt}/C_p(t)+V))^2 \f]
+       Weights can currently be chosen as either 1, or by assuming that \f$C(t)\f$
+       is Poisson distributed (potentially after converting to "counts"),
+       which is a reasonable approximation for OSEM images (although somewhat dangerous
+       for noisy data). The latter is chosen if
+       \c _assume_poisson_distribution is \c true.
+
+       \todo There is currently no check if the time frame definitions from \a dyn_image are
       the same as the ones encoded in the model.
     */
     void 
@@ -139,6 +165,7 @@ class PatlakPlot : public RegisteredParsingObject<PatlakPlot, KineticModel>
   float _time_shift;   //!< Shifts the time to fit the timing of Plasma Data with the Projection Data.
   bool _in_correct_scale; //!< Switch to scale or not the model_matrix to the correct scale, according to the appropriate scale factor.
   bool _in_total_cnt;   //!< Switch to choose the values of the model to be in total counts or in mean counts.
+  bool _assume_poisson_distribution; //!< Assume that image is Poisson distributed, and weight linear regression accordingly (see \c apply_linear_regression())
   std::string _blood_data_filename;   //!< Name of file in which the input function is stored
   PlasmaData _plasma_frame_data;    //!< Stores the plasma data into frames for brain studies
   std::string _time_frame_definition_filename;   //!< name of file to get frame definitions

src/modelling_buildblock/PatlakPlot.cxx

@@ -2,6 +2,7 @@
 //
 /*
     Copyright (C) 2006 - 2011, Hammersmith Imanet Ltd
+    Copyright (C) 2021, University College London
     This file is part of STIR.
 
     SPDX-License-Identifier: Apache-2.0
@@ -12,6 +13,7 @@
   \ingroup modelling
   \brief Implementations of inline functions of class stir::PatlakPlot
   \author Charalampos Tsoumpas
+  \author Ander Biguri
 
   \sa PatlakPlot.h, ModelMatrix.h and KineticModel.h
 
@@ -34,6 +36,7 @@ set_defaults()
   this->_time_shift=0.;
   this->_in_correct_scale=false;
   this->_in_total_cnt=false;
+  this->_assume_poisson_distribution=false;
 }
 
 const char * const 
@@ -191,7 +194,7 @@ PatlakPlot::apply_linear_regression(ParametricVoxelsOnCartesianGrid & par_image,
       frame_num<=num_frames ; ++frame_num )
     {      
       patlak_x[frame_num-1]=patlak_model_array[1][frame_num]/patlak_model_array[2][frame_num];
-      weights[frame_num-1]=1;                    
+      weights[frame_num-1]=1;  // if (_assume_poisson_distribution) handled later in the main loop.
     }   
   {  // Do linear_regression for each voxel // for k j i 
     float slope=0.F;
@@ -219,7 +222,19 @@ PatlakPlot::apply_linear_regression(ParametricVoxelsOnCartesianGrid & par_image,
                 // (remember, these are integrals over the time frame, not single values at discrete t) 
                 for ( frame_num = starting_frame; 
                       frame_num<=num_frames ; ++frame_num )
+                {
                   patlak_y[frame_num-1]=dyn_image[frame_num][k][j][i]/patlak_model_array[2][frame_num];
+                  if(this->_assume_poisson_distribution)
+                  {
+                    weights[frame_num-1]=dyn_image[frame_num][k][j][i]/(patlak_model_array[2][frame_num]*patlak_model_array[2][frame_num]);
+                    // if is zero, leave as is, zero weight?
+                    if(weights[frame_num-1]){
+                      weights[frame_num-1]=1/weights[frame_num-1];
+                    }else{
+                       weights[frame_num-1]=0.00001;
+                    }
+                  }
+                }
                 // Apply the regression to this pixel
                 linear_regression(y_intersection, slope,
                                   chi_square,
@@ -335,6 +350,7 @@ initialise_keymap()
   this->parser.add_key("Time Shift", &this->_time_shift);
   this->parser.add_key("In total counts", &this->_in_total_cnt);
   this->parser.add_key("In correct scale", &this->_in_correct_scale);
+  this->parser.add_key("Poisson distributed images", &this->_assume_poisson_distribution);
   this->parser.add_key("Time Frame Definition Filename", &this->_time_frame_definition_filename); 
   this->parser.add_stop_key("end Patlak Plot Parameters");
 }

src/modelling_utilities/apply_patlak_to_images.cxx

@@ -20,20 +20,7 @@
   apply_patlak_to_images output_parametric_image input_dynamic_image [par_file] 
   \endcode
 
-  \par
-  - The dynamic images will be calibrated only if the calibration factor is given. 
-  - The \a if_total_cnt is set to true the Dynamic Image will have the total number of 
-    counts while if set to false it will have the \a total_number_of_counts/get_duration(frame_num).
-  - The dynamic images will always be in decaying counts.
-  - The plasma data is assumed to be in decaying counts.
-
-  \sa PatlakPlot.h for the \a par_file
-
-  \note This implementation does not use wighted least squares because for Patlak Plot only the last frames are used, which they usually have the same duration and similar number of counts.
-
-  \todo Reimplement the method for image-based input function.
-
-  \todo Add to the Doxygen documentation a reference to their paper and how exactly this utility works.
+  \sa stir::PatlakPlot for more details, including assumptions and the parameter file.
 */
 
 #include "stir/CPUTimer.h"