iPatlakModel.cc

Go to the documentation of this file.

/*
  Implementation of class iPatlakModel

  - separators: X
  - doxygen: X
  - default initialization: none  (require user inputs (Patlak basis functions), no sense to provide 'standard' configuration file as for optimizer/projector)
  - CASTOR_DEBUG: X
  - CASTOR_VERBOSE: X
*/

#include "iPatlakModel.hh"


// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn iPatlakModel
  \brief Constructor of iPatlakModel. Simply set all data members to default values.
*/
iPatlakModel::iPatlakModel() : vDynamicModel() 
{
  m_nbTimeBF = 2; // Two basis functions in the Patlak model
  m2p_parametricImages = NULL;
  m2p_patlakTACs = NULL;
  
  m_fileOptions = "";
  m_listOptions = "";
    
  m_savePImgFlag = true;
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn ~iPatlakModel
  \brief Destructor of iPatlakModel
*/
iPatlakModel::~iPatlakModel() 
{
  if(m_initialized)
  {
    for(int b=0 ; b<m_nbTimeBF ; b++)
    {
      if (m2p_patlakTACs[b]) delete m2p_patlakTACs[b];
      if (m2p_parametricImages[b]) delete m2p_parametricImages[b];
    }
  
    if (m2p_patlakTACs) delete[] m2p_patlakTACs;
    if (m2p_parametricImages) delete[] m2p_parametricImages;
  }
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn ShowHelp
  \brief Print out specific help about the implementation of the Patlak
         model and its initialization
*/
void iPatlakModel::ShowHelp()
{
  cout << "-- This class implements the Patlak Reference Tissue Model : " << endl;
  cout << "-- Patlak CS, Blasberg RG: Graphical evaluation of blood-to-brain transfer constants from multiple-time uptake data" << endl;
  cout << "-- J Cereb Blood Flow Metab 1985, 5(4):5 84-590." << endl;
  cout << "-- DOI http://dx.doi.org/10.1038/jcbfm.1985.87" << endl;
  cout << "-- It is used to model radiotracers which follows as 2-tissue compartment model with irreversible trapping  " << endl;
  cout << "-- The Patlak temporal basis functions are composed of the Patlak slope (integral of the reference TAC from the injection time " << endl;
  cout << "   divided by the instantaneous reference activity), and intercept (reference tissue TAC)  " << endl;
  cout << endl;
  cout << " It can be initialized using either an ASCII file or a list of option with the following keywords and information :" << endl; 
  cout << " - The ASCII file must contain the following keywords :" << endl;
  cout << "   'Patlak_functions:'      (mandatory) Enter the coefficients of Patlak plot and intercept for each time frame (tf) ";
  cout << "                                        on two successive lines, separated by ',' :" << endl;
  cout << "                            -> Patlak_functions: " << endl;
  cout << "                            -> coeff_Pplot_tf1,coeff_Pplot_tf2,...,coeff_Pplot_tfn" << endl;
  cout << "                            -> coeff_Pintc_tf1,coeff_Pintc_tf2,...,coeff_Pintc_tfn" << endl;
  cout << "   'Parametric_image_init:' (optional) path to an interfile image to be used as initialization for the parametric images" << endl;
  cout << "   'Patlak_save_images:'    (optional) boolean indicating if the parametric images should be saved (1) or not (0)" << endl;
  cout << endl;
  cout << " - The list of options must contain the coefficients of both Patlak functions separated by commas, with the following template :" << endl;
  cout << "   coeff_Pplot_tf1,coeff_Pplot_tf2,...,coeff_Pplot_tfn,";
  cout << "   coeff_Pintc_tf1,coeff_Pintc_tf2,...,coeff_Pintc_tfn "<< endl;
  cout << "   Parametric images will be initialized with 1.0 by default " << endl;
  cout << "   The parametric images estimations will be written on disk for each iteration" << endl;
  cout << "    " << endl;
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn ReadAndCheckConfigurationFile
  \param const string& a_configurationFile : ASCII file containing informations about a dynamic model
  \brief This function is used to read options from a configuration file.
  \return 0 if success, other value otherwise.
*/
int iPatlakModel::ReadAndCheckConfigurationFile(string a_fileOptions)
{
  if(m_verbose >=2) Cout("iPatlakModel::ReadAndCheckConfigurationFile ..."<< endl); 
  
  // Recover the file path here, it will be processed in the Initialize() function
  m_fileOptions = a_fileOptions;
  
  ifstream in_file(a_fileOptions.c_str(), ios::in);
  
  if(in_file)
  {
    if( ReadDataASCIIFile(a_fileOptions, "Patlak_save_images", &m_savePImgFlag, 1, KEYWORD_OPTIONAL) == 1)
    {
      Cerr("***** iPatlakModel::ReadAndCheckConfigurationFile -> Error while trying to read 'Patlak_save_images' flag in " << a_fileOptions << endl);
      return 1;
    }
  }
  else
  {
    Cerr("***** iPatlakModel::ReadAndCheckConfigurationFile -> Error while trying to read configuration file at: " << a_fileOptions << endl);
    return 1;
  }
    
  return 0;
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn ReadAndCheckOptionsList
  \param const string& 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 iPatlakModel::ReadAndCheckOptionsList(string a_listOptions)
{
  if(m_verbose >=2) Cout("iPatlakModel::ReadAndCheckOptionsList ..."<< endl); 
  
  // 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 iPatlakModel::CheckSpecificParameters()
{
  if(m_verbose >=2) Cout("iPatlakModel::CheckSpecificParameters ..."<< endl); 
  
  // Check image dimensions
  if (mp_ID==NULL)
  {
    Cerr("***** iPatlakModel::CheckParameters() -> ImageDimensions object has not been provided !" << endl);
    return 1;
  }
  
  // Check number time basis functions
  if (m_nbTimeBF<0)
  {
    Cerr("***** iPatlakModel::CheckParameters() -> Wrong number of time frame basis functions !" << endl);
    return 1;
  }
  
  // Check if we have somehow both a file and a list of options for init...
  if(m_listOptions != "" && m_fileOptions != "")
  {
    Cerr("***** iPatlakModel::Initialize -> Either a file or a list of options have to be selected to initialize the model, but not both ! " << endl);
    return 1;
  }
  
  // Check if we have no file not list of options for some reason...
  if(m_listOptions == "" && m_fileOptions == "")
  {
    Cerr("***** iPatlakModel::Initialize -> Either a file or a list of options should have been provided at this point ! " << endl);
    return 1;
  }
  
  // Check if we reconstruct gated data. Throw warning if it is the case
  if(mp_ID->GetNbRespGates()>1 || mp_ID->GetNbCardGates()>1)
  {
    Cerr("***** iPatlakModel::Initialize -> WARNING : the implemented Patlak model should not be used with gated reconstruction (parametric images will be the same for each gate)! " << endl);
    //return 1;
  }
  
  
  // Normal end
  return 0;
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn Initialize
  \brief This function is used to initialize Patlak parametric images and basis functions
  \return 0 if success, other value otherwise.
*/
int iPatlakModel::Initialize()
{
  if(m_verbose >=2) Cout("iPatlakModel::Initialize ..."<< endl); 

  // Forbid initialization without check
  if (!m_checked)
  {
    Cerr("***** oDynamicModelManager::Initialize() -> Must call CheckParameters functions before Initialize() !" << endl);
    return 1;
  }
  
  // --- Memory Allocation --- //
  
  // Allocate memory for Parametric images and functions of Patlak model
  m2p_patlakTACs = new FLTNB*[m_nbTimeBF];
  m2p_parametricImages = new FLTNB*[m_nbTimeBF];

  for(int b=0 ; b<m_nbTimeBF ; b++)
  {
    m2p_patlakTACs[b] = new FLTNB[mp_ID->GetNbTimeFrames()];
    m2p_parametricImages[b] = new FLTNB[mp_ID->GetNbVoxXYZ()];
  }


  // --- Data Initialization with a configuration file --- //

  if(m_fileOptions != "")
  {
    ifstream in_file(m_fileOptions.c_str(), ios::in);
    
    if(in_file)
    {
      // Patlak basis functions Initialization
      if( ReadDataASCIIFile(m_fileOptions,
                            "Patlak_functions",
                            m2p_patlakTACs,
                            mp_ID->GetNbTimeFrames(),
                            m_nbTimeBF,
                            KEYWORD_MANDATORY) )
      {
        Cerr("***** iPatlakModel::Initialize -> Error while trying to read Patlak functions coefficients !" << endl);
        Cerr("                                  'Patlak_functions' keyword in " << m_fileOptions << endl);
        return 1;
      }
      
      // Patlak Parametric images initialization
      string input_image = "";
      int return_value = 0;
      
      return_value = ReadDataASCIIFile(m_fileOptions,
                                        "Parametric_images_init",
                                        &input_image,
                                        1,
                                        KEYWORD_OPTIONAL);
      
      if( return_value == 0) // Image have been provided
      {
        // Read image // INTF_LERP_DISABLED = interpolation disabled for input image reading
        if( IntfReadImgDynCoeffFile(input_image,
                                    m2p_parametricImages,
                                    mp_ID,
                                    m_nbTimeBF,
                                    m_verbose,
                                    INTF_LERP_DISABLED) ) // Image have been provided
        {
          Cerr("***** iPatlakModel::Initialize -> Error while trying to read the provided initialization parametric images : " << input_image << endl);
          return 1;
        }
        
        //normal end
        return 0;
      }
      else if( return_value == 1) // Error during reading
      {
        Cerr("***** iPatlakModel::Initialize -> Error while trying to read Patlak functions coefficients !" << endl);
        Cerr("                                  'Parametric_image_init' keyword in " << m_fileOptions << endl);
        return 1;
      }
      else //(return_value >= 1 ) // Keyword not found : no initialization provided
      {
        // Standard initialization
        for(int b=0 ; b<m_nbTimeBF ; b++)
          for(int v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
            m2p_parametricImages[b][v] = 1.;
      }
      
    }
    
    else
    {
      Cerr("***** iPatlakModel::Initialize() -> Error while trying to read configuration file at: " << m_fileOptions << endl);
      return 1;
    }
  }
  
  
  // --- Data Initialization with a list of options --- //

  if(m_listOptions != "")
  {
    // We expect here the coefficients of the Patlak functions for each time point
    // Patlak slope before Patlak intercept
  
    // Allocate memory to recover the elements in one tmp vector
    FLTNB *p_coeffs = new FLTNB[m_nbTimeBF * mp_ID->GetNbTimeFrames()];
    
    // Read them
    if (ReadStringOption(m_listOptions,
                         p_coeffs,
                         m_nbTimeBF * mp_ID->GetNbTimeFrames(),
                         ",",
                         "Patlak model configuration"))
    {
      Cerr("***** iPatlakModel::Initialize() -> Failed to correctly read the list of options !" << endl);
      return 1;
    }
    
    // Affect coeffs
    for(int c=0 ; c<m_nbTimeBF * mp_ID->GetNbTimeFrames() ; c++)
    {
      int bf = int(c/mp_ID->GetNbTimeFrames()); // Patlak basis function index
      int fr = int(c%mp_ID->GetNbTimeFrames()); // Frame index
      m2p_patlakTACs[bf][fr] = p_coeffs[c];
    }
    
    // Delete the tmp vector
    delete[] p_coeffs;
    
    
    // Standard initialization for the parametric images
    for(int b=0 ; b<m_nbTimeBF ; b++)
      for(int v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
        m2p_parametricImages[b][v] = 1.;
  }

  // Display Patlak TACs
  if(m_verbose >=2)
  {
    Cout("iPatlakModel::Initialize() -> Patlak Normalized Time TAC coefficients :" << endl);
    Cout("                              ");
    for(int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++)
      Cout(m2p_patlakTACs[0][fr] << ", ");
    Cout(endl);
    Cout("iPatlakModel::Initialize() -> Patlak Reference Time TAC coefficients :" << endl);
    Cout("                              ");
    for(int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++)
      Cout(m2p_patlakTACs[1][fr] << ", ");
    Cout(endl);
  }
  
  // Normal end
  m_initialized = true;
  return 0;
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn EstimateModelParameters
  \param ap_ImageS : pointer to the ImageSpace
  \param a_ite : index of the actual iteration (not used)
  \param a_sset : index of the actual subset (not used)
  \brief Estimate Patlak parametric images
  \return 0 if success, other value otherwise.
*/
int iPatlakModel::EstimateModelParameters(oImageSpace* ap_ImageS, int a_ite, int a_sset) 
{
  if(m_verbose >=2) Cout("iPatlakModel::EstimateModelParameters ..." <<endl);

  #ifdef CASTOR_DEBUG
  if (!m_initialized)
  {
    Cerr("***** iPatlakModel::EstimateModelParameters() -> Called while not initialized !" << endl);
    Exit(EXIT_DEBUG);
  }
  #endif
  
  // Generate estimated image from coefficients and basis functions
  // (We use the backward image which is useless at this point, as temporary image
  // to estimate the image generated using the current estimation of the Patlak parametric images)
  for (int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++)
    for (int rg=0 ; rg<mp_ID->GetNbRespGates() ; rg++) 
      for (int cg=0 ; cg<mp_ID->GetNbCardGates() ; cg++)
        for (int v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
        {
          // Reset this voxel to 0
          ap_ImageS->m6p_backwardImage[0][0][fr][rg][cg][v] = 0;

          // Retrieve current estimation of image according to the current value of Patlak parametric images and basis functions
          // C(fr,v) = Patlak Parametric(v) * normalized time(fr) +
          //           Patlak Intercept(v)  * Reference tissue TAC(fr)
          for (int b=0 ; b<m_nbTimeBF ; b++)
            ap_ImageS->m6p_backwardImage[0][0][fr][rg][cg][v] += m2p_patlakTACs[b][fr] * m2p_parametricImages[b][v];
        }

  // Get correction images using the ratio of the current estimation of the image (m4p_image)
  // and the image generated with the combination of Patlak functions and parametric images (stored in ap_Image->m6p_backwardImage)
  // (Again, use backward image as temporary image to get the result)
  for (int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++)
    for (int rg=0 ; rg<mp_ID->GetNbRespGates() ; rg++) 
      for (int cg=0 ; cg<mp_ID->GetNbCardGates() ; cg++)
        for (int v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
          ap_ImageS->m6p_backwardImage[0][0][fr][rg][cg][v] = ap_ImageS->m4p_image[fr][rg][cg][v] / ap_ImageS->m6p_backwardImage[0][0][fr][rg][cg][v];


  // Loop on the 2 Patlak parametric images
  for (int b=0 ; b<m_nbTimeBF ; b++)
  {
    double temporal_basis_functions_norm = 0;
    
    // Compute normalization related to the temporal basis functions.
    for (int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++)
      temporal_basis_functions_norm += m2p_patlakTACs[b][fr];
    
    // Compute voxelwise corrections for the parametric images
    int v;
    #pragma omp parallel for private(v) schedule(guided)
    for (v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
    {
      double corr_factor = 0.;
      
      // Note : We don't consider here parametric images specific to respiratory (rg) or cardiac gates (cg)
      //        The model is applied to every voxels in the image
      for (int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++)
        for (int rg=0 ; rg<mp_ID->GetNbRespGates() ; rg++) 
          for (int cg=0 ; cg<mp_ID->GetNbCardGates() ; cg++)
            corr_factor += ap_ImageS->m6p_backwardImage[0][0][fr][rg][cg][v] * m2p_patlakTACs[b][fr];

      // Apply corrections and normalization to the vowelwise time basis functions coefficients.
      if (corr_factor > 0.) 
        m2p_parametricImages[b][v] *= corr_factor/temporal_basis_functions_norm; 
    }
  }
  
  return 0;
}
  
  
  
// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn FitModel
  \param ap_ImageS : pointer to the ImageSpace
  \param a_ite : index of the actual iteration (not used)
  \param a_sset : index of the actual subset (not used)
  \brief Estimate image using Patlak parametric images and basis functions
  \return 0 if success, other value otherwise.
*/
int iPatlakModel::FitModel(oImageSpace* ap_ImageS, int a_ite, int a_sset) 
{
  if(m_verbose >= 2) Cout("iPatlakModel::FitModel ... " <<endl);

  #ifdef CASTOR_DEBUG
  if (!m_initialized)
  {
    Cerr("***** iPatlakModel::FitModel() -> Called while not initialized !" << endl);
    Exit(EXIT_DEBUG);
  }
  #endif
  
  for (int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++)
    for (int rg=0 ; rg<mp_ID->GetNbRespGates() ; rg++) 
      for (int cg=0 ; cg<mp_ID->GetNbCardGates() ; cg++)
        for (int v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
        {
          // Reset current estimated value
          ap_ImageS->m4p_image[fr][rg][cg][v] = 0;
          
          // C(fr,v) = Patlak Parametric(v) * normalized time(fr) +
          //           Patlak Intercept(v)  * Reference tissue TAC(fr) 
          for (int b=0 ; b<m_nbTimeBF ; b++)
            ap_ImageS->m4p_image[fr][rg][cg][v] += m2p_parametricImages[b][v] * m2p_patlakTACs[b][fr];
        }
        
  return 0;
}





// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn SaveParametricImages
  \param ap_ImageS : pointer to the ImageSpace
  \param a_ite : index of the actual iteration
  \brief Write parametric images on disk if 'm_savePImgFlag' is enabled
  \todo Interfile management
  \return 0 if success, other value otherwise.
*/
int iPatlakModel::SaveParametricImages(int a_ite)
{
  if(m_verbose >=2) Cout("iPatlakModel::SaveParametricImages ..." <<endl);
  
  if(m_savePImgFlag)
  {
    // Get the output manager
    sOutputManager* p_output_manager = sOutputManager::GetInstance();  
    /*
    for(int bimg=0 ; bimg<m_nbTimeBF ; bimg++)
    {
      string data_file = p_output_manager->GetPathName() + p_output_manager->GetBaseName();
      if (a_ite >= 0) // Add a suffix for iteration
      {
        stringstream ss;
        ss << a_ite;
        data_file.append("_ite_").append(ss.str());
      }
    
      // Add a suffix for (basis functions) coefficients images
      stringstream ss;
      ss << bimg;
      data_file.append("_fbimg_").append(ss.str()).append(".bin");
      
      ofstream out_dfile;
    
      out_dfile.open(data_file.c_str(), ios::binary | ios::out);
          
      if(!out_dfile.is_open()) 
      {
        cout<<"***** iPatlakModel::SaveCoeffImages()->Failed to create output file for the dynamic image : "<< data_file.c_str() << endl;
        return 1;
      }
            
      out_dfile.write(reinterpret_cast<char*>(m2p_parametricImages[bimg]), mp_ID->GetNbVoxXYZ()*sizeof(FLTNB));
      out_dfile.close();
    }
    */
    
    // Interfile
    string path_to_image = p_output_manager->GetPathName() + p_output_manager->GetBaseName();
    
    // Add a suffix for iteration
    if (a_ite >= 0)
    {
      stringstream ss; ss << a_ite + 1;
      path_to_image.append("patlak_it").append(ss.str());
    }
    
    // Write interfile parametric image
    if(IntfWriteImgDynCoeffFile(path_to_image, 
                                m2p_parametricImages, 
                                mp_ID, 
                                m_nbTimeBF, 
                                m_verbose) )
    {
      Cerr("***** iPatlakModel::SaveParametricImages()-> Error writing Interfile of output image !" << endl);  
      return 1;
    }
    
  }
  
  return 0;
}