iLinearSpectralModel.cc

Go to the documentation of this file.

/*
This file is part of CASToR.

    CASToR is free software: you can redistribute it and/or modify it under the
    terms of the GNU General Public License as published by the Free Software
    Foundation, either version 3 of the License, or (at your option) any later
    version.

    CASToR is distributed in the hope that it will be useful, but WITHOUT ANY
    WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
    FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
    details.

    You should have received a copy of the GNU General Public License along with
    CASToR (in file GNU_GPL.TXT). If not, see <http://www.gnu.org/licenses/>.

Copyright 2017-2019 all CASToR contributors listed below:

    --> Didier BENOIT, Claude COMTAT, Marina FILIPOVIC, Thibaut MERLIN, Mael MILLARDET, Simon STUTE, Valentin VIELZEUF

This is CASToR version 3.0.
*/

#include "iLinearSpectralModel.hh"
#include "cmath"
#include "math.h"



// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn iDynamicModelTemplate
  \brief Constructor of iDynamicModelTemplate. Simply set all data members to default values.
*/
iLinearSpectralModel::iLinearSpectralModel() : iLinearModel()
{
  // --- Parameters inherited from vDynamicModel class --- //
  
  m_nbTimeBF = -1; // Number of basis functions in the model un-initialised

  // Negative initialization of specific model values
  m_fast_exp=-1;
  m_slow_exp=-1;
  m_constant_basis_value=-1;
  m_spectral_bank_size=-1;
  mp_spectral_bank=NULL;

  // Initialize number of additional basis functions to 0;
  m_additionalBF_size = 0  ;

}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn ~iDynamicModelTemplate
  \brief Destructor of iDynamicModelTemplate
*/
iLinearSpectralModel::~iLinearSpectralModel()
{

}


// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn ShowHelp
  \brief Print out specific help about the implementation of the model 
         and its initialization
*/
void iLinearSpectralModel::ShowHelp()
{
  cout << "-- This class implements the Spectral Model : " << endl;
  cout << "-- By Adnrew Reader et al. " << endl;
  cout << "-- It is used to model radiotracer dynamic behaviour with a set of decaying exponential functions ( exp(-β*t) )  " << endl;
  cout << "-- The decaying exponential functions are logarithmically equally spaced within the selected range of β values   " << endl;
  cout << "   All decaying exponentials are convolved with the interpolated Arterial Input Curve, before being discretised to the duration of the reconstruction frames  " << endl;
  cout << endl;
  cout << "  The model must be initialized using either an ASCII file with the following keywords and information :" << endl;
  cout << "   As this class inherits from the iLinearModel class, the following parameters must be declared inside the couple of the following specific tags: " << endl;
  cout << "   - DYNAMIC FRAMING/ENDDF " << endl;
  cout << "   - The ASCII file must contain the following keywords :" << endl;
  cout << "    'AIC_input_file:'       (mandatory) The file containing the sampled Arterial Input Function " << endl;
  cout << "                               The file must contain the following information in successive lines, separated by ',' " << endl;
  cout << "                            -> AIC_number_of_points: " << endl;
  cout << "                            -> AIC_time_points: " << endl;
  cout << "                            -> AIC_data_points: " << endl;
  cout << "                            -> AIC_units: seconds " << endl;
  cout << "   'Parametric_image_init: path ' (optional) path to an interfile image to be used as initialization for the parametric images." << endl;
  cout << "   -> Optimisation_method : x       (mandatory) optimization method available options: " << endl;
  cout << "                                    x=0: Direct ( Implementation of basis functions side by system matrix in each tomographic iterative loop " << endl;
  cout << "                                    x=1: Nested EM  " << endl;
  cout << "                                    x=2: Iterative non-negative Least-Square " << endl;
  cout << "                                         (C.L. Lawson and R.J. Hanson, Solving Least Squares Problems)" << endl;
  cout << endl;
  cout << "   Parametric images will be initialized with 0.001 and 1.0 for Patlak slope and intercept by default " << endl;
  cout << "   The parametric images estimations will be written on disk for each iteration" << endl;
  cout << "    " << endl;
  cout << "   The following keywords are common to all dynamic models :" << endl;
  cout << "   'Number of iterations before image update: x' Set a number 'x' of iteration to reach before using the model to generate the images at each frames/gates" << endl;
  cout << "   (Default x ==  0) " << endl;
  cout << "   'No image update: x'                          If set to 1, the reconstructed images for the next iteration/subset are not reestimated using the model" << endl;
  cout << "   (Default x ==  0)                              (the code just performs standard independent reconstruction of each frames/gates) " << endl;
  cout << "   'No parameters update: x'                     If set to 1, the parameters / functions of the model are not estimated with the image" << endl;
  cout << "   (Default x ==  0)                              (this could be used to test The EstimateImageWithModel() function with specific user-provided parametric images) " << endl;
  cout << "   'Save parametric images : x'                  Enable (1)/Disable(0) saving parametric images on disk" << endl;
  cout << "   (Default x == 1)  " << endl;
  cout << "   'Save blacklisted voxels images : x'          Enable (1)/Disable(0) saving blacklisted voxels images on disk" << endl;
  cout << "   (Default x == 0)  " << endl;
  cout << "   " << endl;

}


// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn ReadAndCheckConfigurationFileSpecific
  \brief This function is used to read options from a configuration file.
  \return 0 if success, other value otherwise.
*/
int iLinearSpectralModel::ReadAndCheckConfigurationFileSpecific()
{
  if(m_verbose >=3) Cout("iLinearSpectralModel::ReadAndCheckConfigurationFileSpecific ..."<< endl);
  
  // ===================================================================
  // Implement here the reading of any options specific to this dynamic model 
  // (i.e : parameters or path to deformation files), through a configuration file
  // The ReadDataASCIIFile() functions could be helpful to recover data from a file
  // (check other dynamicModels for examples)
  // ===================================================================

  // Apply the Generic linear Checks for all Linear Models
  if( ReadAndCheckConfigurationFileSpecificToAllLinearModels()==1)
  {
    Cerr("***** iLinearSpectralModel::ReadAndCheckConfigurationFileSpecific -> Error while trying to read and check generic configuration for linear models " << m_fileOptions << endl);
    return 1;
  }

  // The file will be fully processed in the Initialize() function
  ifstream in_file(m_fileOptions.c_str(), ios::in);

  if(in_file)
  {
    // Number of requested basis/spectral functions
    if( ReadDataASCIIFile(m_fileOptions, "Number_spectral_functions", &m_spectral_bank_size, 1, KEYWORD_MANDATORY) == 1)
    {
      Cerr("***** iLinearModel::ReadAndCheckConfigurationFileSpecific -> Error while trying to read 'Number_spectral_functions' flag in " << m_fileOptions << endl);
      return 1;
    }

    // A basis function is added to the total - to account for constant exponential -> Arterial Input Basis Function
    m_additionalBF_size ++;

    // Is constant value basis function requested ?
    m_constant_basis_value =0 ;
    if( ReadDataASCIIFile(m_fileOptions, "Constant_basis_function", &m_constant_basis_value, 1, KEYWORD_OPTIONAL) == 1)
    {
      Cerr("***** iLinearModel::ReadAndCheckConfigurationFileSpecific -> Error while trying to read 'Constant_basis_function' flag in " << m_fileOptions << endl);
      return 1;
    }
    // Add another basis function if constant_basis value has been provided
    if (m_constant_basis_value>0) m_additionalBF_size ++;

    // Set model parameters and number of basis functions equal to the number of basis functions (including the additional)
    m_nbModelParam = m_spectral_bank_size + m_additionalBF_size ;
    m_nbTimeBF = m_spectral_bank_size + m_additionalBF_size ;

    // Parameters for spectral Model -> Lower, Upper decay rate constants
    if( ReadDataASCIIFile(m_fileOptions, "Fastest_rate", &m_fast_exp, 1, KEYWORD_MANDATORY) == 1)
    {
      Cerr("***** iLinearSpectralModel::ReadAndCheckConfigurationFileSpecific -> Error while trying to read 'Fastest_rate' flag in " << m_fileOptions << endl);
      return 1;
    }
    // Make rate units to 1/sec
    m_fast_exp/=60. ;

    if( ReadDataASCIIFile(m_fileOptions, "Slowest_rate", &m_slow_exp, 1, KEYWORD_MANDATORY) == 1)
    {
      Cerr("***** iLinearSpectralModel::ReadAndCheckConfigurationFileSpecific -> Error while trying to read 'Slowest_rate' flag in " << m_fileOptions << endl);
      return 1;
    }
    // Make rate units to 1/sec
    m_slow_exp/=60. ;

  }
  else
  {
    Cerr("***** iLinearPatlakModel::ReadAndCheckConfigurationFileSpecific -> Error while trying to read configuration file at: " << m_fileOptions << endl);
    return 1;
  }


return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn ReadAndCheckOptionsList
  \param a_optionsList : a list of parameters separated by commas
  \brief This function is used to read parameters from a string.
  \return 0 if success, other value otherwise.
*/
int iLinearSpectralModel::ReadAndCheckOptionsList(string a_listOptions)
{
  if(m_verbose >=3) Cout("iLinearSpectralModel::ReadAndCheckConfigurationFileSpecific ..."<< endl);

  // Must be initialized with configuration file
  Cerr("***** iLinearSpectralModel::ReadAndCheckOptionsList -> This model needs to use a file for its initialization. Use help for more info !" << endl);
  return 1;
  
  // ===================================================================
  // Implement here the reading of any options specific to this deformation model,
  // through a list of options separated by commas
  // The ReadStringOption() function could be helpful to parse the list of parameters in an array
  // ===================================================================

  // Just recover the string here, it will be processed in the Initialize() function
  //m_listOptions = a_listOptions;

  // Normal end
  //return 0;
}


// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn CheckSpecificParameters
  \brief This function is used to check whether all member variables
         have been correctly initialized or not.
  \return 0 if success, positive value otherwise.
*/
int iLinearSpectralModel::CheckSpecificParameters()
{
  // ===================================================================
  // Implement here checks over parameters which should be read using either
  // ReadAndCheckConfigurationFile() and ReadAndCheckOptionsList() functions
  // ===================================================================

  // Perform generic checks that apply for the Linear Models
  if(CheckSpecificParametersForAllLinearModels())
  {
    Cerr("***** iLinearPatlakModel::CheckSpecificParameters -> A problem occurred while checking specific parameters ! " << endl);
    return 1;
  }


  // Normal end
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn InitializeSpecific
  \brief This function is used to initialize the model parametric images and basis functions
  \return 0 if success, other value otherwise.
*/
int iLinearSpectralModel::InitializeSpecific()
{
  if (!m_checked)
  {
    Cerr("***** iLinearSpectralModel::InitializeSpecific() -> Must call CheckParameters functions before Initialize() !" << endl);
    return 1;
  }


  if(m_verbose >=2) Cout("iLinearSpectralModel::InitializeSpecific ..."<< endl);
  // Run generic Initialization for all Linear Models
  if (InitializeSpecificToAllLinearModels())
  {
    Cerr("***** iLinearSpectralModel::InitializeSpecific() -> Error while performing generic initialisations for linear models !" << endl);
    return 1;
  }

  // ===================================================================
  // Implement here the allocation/initialization of whatever member
  // variables specifically used by this deformation model
  // The set 
  // ===================================================================

  // --- Calculation of spectral function decay rates --------------------------------------------------------- //
  // allocate spectral bank
  HPFLTNB* mp_spectral_bank = (HPFLTNB*)malloc(m_spectral_bank_size*sizeof(HPFLTNB));
  // Calculate the spectral coefficients to be equally spaced in a log scale
  HPFLTNB ln_slow_rate = log10(m_slow_exp);
  HPFLTNB ln_fast_rate = log10(m_fast_exp);
  HPFLTNB step = (ln_fast_rate - ln_slow_rate) / (HPFLTNB) (m_spectral_bank_size - 1);
  for (int b = 0; b < m_spectral_bank_size; b++) mp_spectral_bank[b] = ln_slow_rate + step * b;
  // Get spectral coefficients back to scale from log scale
  for (int b = 0; b < m_spectral_bank_size; b++) mp_spectral_bank[b] = pow(10, mp_spectral_bank[b]);

  // --- Sample the exponential over the requested frames duration --------------------------------------------- //
  if(m_verbose >=3) Cout( " iLinearSpectralModel::InitializeSpecific() -> Sampling exponential functions" << endl);
  // Assumption of only one bed position at the moment !!
  int bed = 0;
  // Calculated total samples for interval of 0.1 seconds
  int total_samples =
      (int) ((((float) mp_ID->GetFrameTimeStopInMs(bed, mp_ID->GetNbTimeFrames() - 1)) / 1000.) * 10);

  HPFLTNB** spectral_exp = (HPFLTNB **) malloc((m_spectral_bank_size) * sizeof(HPFLTNB*));
  for (int ex = 0; ex < m_spectral_bank_size; ex++)
  {
    spectral_exp[ex] = (HPFLTNB*) malloc(total_samples * sizeof(HPFLTNB));
  }

  // Evaluating exponents over total time. Scaling time for 0.1 second intervals
  for (int ex = 0; ex < m_spectral_bank_size; ex++)
  {
    for (int i = 0; i < total_samples; i++)
    {
      spectral_exp[ex][i] = (FLTNB) exp(-0.1 * mp_spectral_bank[ex] * i);
    }
  }


  // --- Get Interpolated AIC and convolve with every spectral function --------------------------------------------- //

  // Allocate memory for the convolved spectral functions.
  HPFLTNB** ConvSpectralFunctions = (HPFLTNB **) malloc((m_spectral_bank_size) * sizeof(HPFLTNB*));
  for (int ex = 0; ex < m_spectral_bank_size; ex++)
  {
    ConvSpectralFunctions[ex] = (HPFLTNB*) malloc(total_samples * sizeof(HPFLTNB));
  }
  // Get Downsampled AIF and get it into a_AICIntrpY
  mp_ArterialInputCurve->Downsample();
  HPFLTNB *a_AICIntrpY = mp_ArterialInputCurve->GetDownsampledAIC();
  // For every function perform the convolution ! use of multithreading for each spectral function convolution
  int ex;
  #pragma omp parallel for private(ex) schedule(static, 1)
  for (ex = 0; ex < m_spectral_bank_size; ex++)
  {
    // Double loop for convolution of AICIntrpY with a sampled exponential function
    for (int i = 0; i <= total_samples; i++)
    {
      for (int j = i; j <= total_samples; j++)
      {
        ConvSpectralFunctions[ex][j] += (a_AICIntrpY[i] * spectral_exp[ex][j - i]);
      }
      // Normalise each value with the size of each step
      ConvSpectralFunctions[ex][i] *= 0.1;
      // Check and remove negative values - might occur if multiplication results to lower/higher values than double !!
      if ( ConvSpectralFunctions[ex][i] <0 )
      {
        Cout(" Negative spectral function value detected at " << i << ", with value: "<< ConvSpectralFunctions[ex][i] <<
         "  -> Setting value to zero "<< endl);
        ConvSpectralFunctions[ex][i] = 0 ;
      }
      if (std::isnan(ConvSpectralFunctions[ex][i]))
      {
        Cout(" NaN spectral function value detected at " << i << "--> Setting value to zero "<< endl);
        ConvSpectralFunctions[ex][i] = 0 ;
      }
    }
  }

  // Average functions over frames to calculate Basis Functions
  HPFLTNB RunningSum=0 ;
  HPFLTNB halfDeltaT = (0.1/2);
  // Loop over all spectral functions
  for (int ex = 0; ex < m_spectral_bank_size; ex++)
  {
    // Loop over frames to get the sum of the values within each frame
    for (int fr = 0; fr < mp_ID->GetNbTimeFrames(); fr++)
    {
      RunningSum  = 0 ;
      // iteration over all the datapoints within the frame
      for (uint32_t i = (mp_ID->GetFrameTimeStartInMs(bed, fr)/100)+1 ;i <(mp_ID->GetFrameTimeStopInMs(bed, fr)/100);i++)
      {
        RunningSum+=ConvSpectralFunctions[ex][i];
      }
      // Calculate integral using the trapezoidal method
      m2p_nestedModelTimeBasisFunctions[ex][fr] = (FLTNB) ((ConvSpectralFunctions[ex][(mp_ID->GetFrameTimeStartInMs(bed, fr) / 100)] +
                  2 * RunningSum + ConvSpectralFunctions[ex][(mp_ID->GetFrameTimeStopInMs(bed, fr) / 100)]) * halfDeltaT);
      // Average for the frame duration
      m2p_nestedModelTimeBasisFunctions[ex][fr] /=  (FLTNB)((mp_ID->GetFrameTimeStopInMs(bed, fr)) - (mp_ID->GetFrameTimeStartInMs(bed, fr)))/1000 ;
    }
  }

  // Add basis function equal to AIF
  for (int fr = 0; fr < mp_ID->GetNbTimeFrames(); fr++)
  {
    RunningSum = 0;
    // iteration over all the datapoints within the frame
    for (uint32_t i = (mp_ID->GetFrameTimeStartInMs(bed, fr)/100)+1 ;i <(mp_ID->GetFrameTimeStopInMs(bed, fr)/100);i++)
    {
      RunningSum += a_AICIntrpY[i];
    }
    // Calculate integral using the trapezoidal method
    m2p_nestedModelTimeBasisFunctions[m_spectral_bank_size][fr] = (FLTNB) ((a_AICIntrpY[(mp_ID->GetFrameTimeStartInMs(bed, fr) / 100)] +
                       2 * RunningSum + a_AICIntrpY[(mp_ID->GetFrameTimeStopInMs(bed, fr) / 100)]) * halfDeltaT);
    // Average for the frame duration
    m2p_nestedModelTimeBasisFunctions[m_spectral_bank_size][fr] /=   (FLTNB)((mp_ID->GetFrameTimeStopInMs(bed, fr)) - (mp_ID->GetFrameTimeStartInMs(bed, fr)))/1000 ;
  }

  // Add constant basis function if requested
  if (m_constant_basis_value>0)
  {
   // Add constant value as the last basis function
    for (int fr = 0; fr < mp_ID->GetNbTimeFrames(); fr++)
    {
      m2p_nestedModelTimeBasisFunctions[m_nbTimeBF-1][fr] = (FLTNB)1. ;
    }
  }

  // Print the basis functions
  ShowBasisFunctions();

  // Normal end
  m_initialized = true;
  return 0;
}