oSensitivityGenerator.cc

Go to the documentation of this file.

/*
  Implementation of class oSensitivityGenerator

  - separators: 
  - doxygen: 
  - default initialization: 
  - CASTOR_DEBUG: 
  - CASTOR_VERBOSE: 
*/

#include "gVariables.hh"
#include "oSensitivityGenerator.hh"
#include "iDataFilePET.hh"

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

oSensitivityGenerator::oSensitivityGenerator()
{
  // Simply default all members to "empty" values
  m_pathToSensitivityImage = "";
  mp_ImageDimensionsAndQuantification = NULL;
  mp_ImageSpace = NULL;
  mp_Scanner = NULL;
  mp_ProjectorManager = NULL;
  mp_ImageConvolverManager = NULL;
  mp_DeformationManager = NULL;
  m2p_DataFile = NULL;
  mp_pathToNormalizationFileName = {};
  m3p_NormalizationDataFile = NULL;
  m_OneNormalizationFileForAllBeds = false;
  m_OneNormalizationFileForAllFrames = false;
  m_pathToAttenuationImage = "";
  m_attenuationFlag = false;
  m_forwardProjectAttenuation = false;
  m_nbAtnRespGateImages = -1;
  m_nbAtnCardGateImages = -1;
  mp_lineCounter = NULL;
  m_flagGPU = false;
  m_verbose = -1;
  m_checked = false;
  m_initialized = false;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

oSensitivityGenerator::~oSensitivityGenerator()
{
  // Delete all normalization data files if any
  if (m3p_NormalizationDataFile)
  {
    for (int bed=0; bed<mp_ImageDimensionsAndQuantification->GetNbBeds(); bed++)
    {
      if (m3p_NormalizationDataFile[bed])
      {
        for (int fr=0; fr<mp_ImageDimensionsAndQuantification->GetNbTimeFrames(); fr++)
          if (m3p_NormalizationDataFile[bed][fr]) delete m3p_NormalizationDataFile[bed][fr];
        free(m3p_NormalizationDataFile[bed]);
      }
      free(m3p_NormalizationDataFile);
    }
  }
  if (mp_lineCounter) free(mp_lineCounter);
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int oSensitivityGenerator::CheckParameters()
{
  // Check all mandatory parameters
  if (mp_ImageDimensionsAndQuantification==NULL)
  {
    Cerr("***** oSensitivityGenerator::CheckParameters() -> Image dimensions and quantification object is null !" << endl);
    return 1;
  }
  if (mp_ImageSpace==NULL)
  {
    Cerr("***** oSensitivityGenerator::CheckParameters() -> Image space object is null !" << endl);
    return 1;
  }
  if (mp_Scanner==NULL)
  {
    Cerr("***** oSensitivityGenerator::CheckParameters() -> Scanner object is null !" << endl);
    return 1;
  }
  if (mp_ProjectorManager==NULL)
  {
    Cerr("***** oSensitivityGenerator::CheckParameters() -> Projector manager object is null !" << endl);
    return 1;
  }
  if (mp_ImageConvolverManager==NULL)
  {
    Cerr("***** oSensitivityGenerator::CheckParameters() -> Convolver manager object is null !" << endl);
    return 1;
  }
  if (mp_DeformationManager==NULL)
  {
    Cerr("***** oSensitivityGenerator::CheckParameters() -> Deformation manager object is null !" << endl);
    return 1;
  }
  if (m2p_DataFile==NULL)
  {
    Cerr("***** oSensitivityGenerator::CheckParameters() -> Data file array is null !" << endl);
    return 1;
  }
  for (int b=0; b<mp_ImageDimensionsAndQuantification->GetNbBeds(); b++)
  {
    if (m2p_DataFile[b]==NULL)
    {
      Cerr("***** oSensitivityGenerator::CheckParameters() -> Data file object for bed " << b+1 << " is null !" << endl);
      return 1;
    }
  }
  if (m_verbose<0)
  {
    Cerr("***** oSensitivityGenerator::CheckParameters() -> Verbose level is negative !" << endl);
    return 1;
  }
  // In SPECT with attenuation correction
  if (m_pathToAttenuationImage!="" && m2p_DataFile[0]->GetDataType()==TYPE_SPECT)
  {
    // Check that the projection line strategy is not IMAGE_COMPUTATION_STRATEGY (this is not compatible)
    // Note that we cannot do this check in the oProjectionLine directly, because we cannot now in advance
    // that the projection method including attenuation will be used...
    if (mp_ProjectorManager->GetComputationStrategy()==IMAGE_COMPUTATION_STRATEGY)
    {
      Cerr("***** oSensitivityGenerator::CheckParameters() -> The image-computation strategy of the oProjectionLine is not compatible with SPECT attenuation correction !");
      return 1;
    }
    // Check that the projector is compatible with SPECT with attenuation correction
    if (!mp_ProjectorManager->IsBackwardOperatorCompatibleWithSPECTAttenuationCorrection())
    {
      Cerr("***** oSensitivityGenerator::CheckParameters() -> The backward projector is not compatible with SPECT attenuation correction !" << endl);
      return 1;
    }
  }
  // Now it is checked
  m_checked = true;
  // End
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int oSensitivityGenerator::Initialize()
{
  // First check that the parameters has been checked !
  if (!m_checked)
  {
    Cerr("***** oSensitivityGenerator::Initialize() -> Must call the CheckParameters() function before initialization !" << endl);
    return 1;
  }
  // Verbose
  if (m_verbose>=2) Cout("oSensitivityGenerator::Initialize() -> Start initialization" << endl);
  // For the moment, SPECT sensitivity is not implemented
  if (m2p_DataFile[0]->GetDataType()==TYPE_SPECT)
  {
    Cerr("***** oSensitivityGenerator::Initialize() -> Sensitivity for SPECT not yet implemented !" << endl);
    return 1;
  }
  // For the moment, PET sensitivity when using TOF is not possible, will soon be !
  if (m2p_DataFile[0]->GetDataType()==TYPE_PET && (dynamic_cast<iDataFilePET*>(m2p_DataFile[0]))->GetTOFInfoFlag())
  {
    Cerr("***** oSensitivityGenerator::Initialize() -> Sensitivity for PET with TOF is not yet possible, will soon be !" << endl);
    return 1;
  }
  // Initialize normalization files (must be done before attenuation files, because attenuation can be inside the normalization data file)
  if (InitializeNormalizationFiles())
  {
    Cerr("***** oSensitivityGenerator::Initialize() -> A problem occured while initializing the normalization data files !" << endl);
    return 1;
  }
  // Initialize attenuation related stuff
  if (InitializeAttenuationFiles())
  {
    Cerr("***** oSensitivityGenerator::Initialize() -> A problem occured while initializing the attenuation files !" << endl);
    return 1;
  }
  // We want to know if we will need to forward project into the attenuation image, so here the full condition
  if (
       m_attenuationFlag && // 1. we must have an attenuation image
       m2p_DataFile[0]->GetDataType()==TYPE_PET && // 2. This is only for PET
       ( m3p_NormalizationDataFile==NULL || // 3. Either there are no normalization data file (potentially including ACF)
         !(dynamic_cast<iDataFilePET*>(m3p_NormalizationDataFile[0][0]))->GetAtnCorrectionFlag() ) // 4. Or the normalization data file does not contain ACF
     )
  {
    // We use a specific boolean that says we need to forward project into the attenuation image
    m_forwardProjectAttenuation = true;
  }
  else
  {
    // Otherwise, we don't
    m_forwardProjectAttenuation = false;
  }
  // Allocate the backward image for multi-thread
  mp_ImageSpace->LMS_InstantiateBackwardImage();
  // Create sensitivity image file name
  sOutputManager* p_outputManager = sOutputManager::GetInstance();
  m_pathToSensitivityImage = p_outputManager->GetPathName() + p_outputManager->GetBaseName() + "_sensitivity";
  // Allocate line counters
  mp_lineCounter = (uint64_t*)calloc(mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection(),sizeof(uint64_t));
  // End
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int oSensitivityGenerator::InitializeNormalizationFiles()
{
  // TODO: check the consistency between all normalization files (atnflag, normfactor, etc)

  // If no normalization file name provided, then we will exit, but we do some checks before
  if (mp_pathToNormalizationFileName.size()==0)
  {
    // If the ignore-normalization-correction-flag is on, then we can exit the function
    if (mp_ImageDimensionsAndQuantification->GetIgnoreNormCorrectionFlag()) return 0;
    // In PET, throw an error if normalization correction is in the datafile and is taken into account
    if (m2p_DataFile[0]->GetDataType()==TYPE_PET)
    {
      // Cast the vDataFile into iDataFilePET
      iDataFilePET* p_pet_file = (dynamic_cast<iDataFilePET*>(m2p_DataFile[0]));
      // Check if the normalization data are in the file and we do not ignore it
      if (p_pet_file->GetNormCorrectionFlag())
      {
        Cerr("***** oSensitivityGenerator::InitializeNormalizationFiles() -> Normalization correction is included in the data file while it is not in the sensitivity computation !" << endl);
        return 1;
      }
    }
    // We can exit now
    return 0;
  }

  // Verbose
  if (m_verbose>=2) Cout("oSensitivityGenerator::InitializeNormalizationFiles() -> Initialize normalization files" << endl);

  // Check that the number of normalization files options is the same as the number of beds
  if (mp_pathToNormalizationFileName.size()!=((size_t)(mp_ImageDimensionsAndQuantification->GetNbBeds())))
  {
    // If only one is provided, then we assume that it is the same for all bed positions
    if (mp_pathToNormalizationFileName.size()==1) m_OneNormalizationFileForAllBeds = true;
    // Otherwise, throw an error
    else
    {
      Cerr("***** oSensitivityGenerator::InitializeNormalizationFiles() -> Number of normalization files options should be one or equal to the number of data files !" << endl);
      return 1;
    }
  }

  // Allocate the normalization data files array for the bed positions
  m3p_NormalizationDataFile = (vDataFile***)malloc(mp_ImageDimensionsAndQuantification->GetNbBeds()*sizeof(vDataFile**));

  // Get the scanner manager
  sScannerManager* p_ScannerManager = sScannerManager::GetInstance();

  // Loop on the number of beds
  for (int bed=0; bed<mp_ImageDimensionsAndQuantification->GetNbBeds(); bed++)
  {
    // --------------------------------------------------------------------------
    // Step 0: When only one normalization file names is provided for all beds
    // --------------------------------------------------------------------------
    if (m_OneNormalizationFileForAllBeds)
    {
      // If not the first bed
      if (bed!=0)
      {
        // We copy the pointer of the first bed into this one
        m3p_NormalizationDataFile[bed] = m3p_NormalizationDataFile[0];
        // And we continue the loop
        continue;
      }
    }
    // --------------------------------------------------------------------------
    // Step 1: Get filenames and check consistency with number of beds and frames
    // --------------------------------------------------------------------------
    // We will get the normalization file name for this bed
    vector<string> normalization_files;
    // Then we search for commas separating the file name associated to each frame
    size_t comma = mp_pathToNormalizationFileName[bed].find_first_of(",");
    // Loop to get all file names
    while (comma!=string::npos)
    {
      // Get the file name before the comma
      normalization_files.push_back(mp_pathToNormalizationFileName[bed].substr(0,comma));
      // Remove this part from the collection of file names
      mp_pathToNormalizationFileName[bed] = mp_pathToNormalizationFileName[bed].substr(comma+1);
      // Search for the next comma
      comma = mp_pathToNormalizationFileName[bed].find_first_of(",");
    }
    // Get the last file name
    normalization_files.push_back(mp_pathToNormalizationFileName[bed]);
    // Check that the number of file names found is equal to the number of frames
    if (normalization_files.size()!=((size_t)(mp_ImageDimensionsAndQuantification->GetNbTimeFrames())))
    {
      // If only one is provided, then we assume that it is the same for all frames
      if (mp_pathToNormalizationFileName.size()==1) m_OneNormalizationFileForAllFrames = true;
      // Otherwise, throw an error
      else
      {
        Cerr("***** oSensitivityGenerator::InitializeNormalizationFiles() -> Number of normalization files for bed " << bed+1 << " should be one or equal to the number of frames !" << endl);
        return 1;
      }
    }
    // If this is the case, then check for previous beds that it was so
    else if (bed!=0 && !m_OneNormalizationFileForAllFrames) 
    {
      Cerr("***** oSensitivityGenerator::InitializeNormalizationFiles() -> If only one normalization file is provided for all frames, then it must be so for all beds !" << endl);
      return 1;
    }
    // --------------------------------------------------------------------------
    // Step 2: Create the normalization data file objects
    // --------------------------------------------------------------------------
    // Allocate the normalization data file array for the frames
    m3p_NormalizationDataFile[bed] = (vDataFile**)malloc(mp_ImageDimensionsAndQuantification->GetNbTimeFrames()*sizeof(vDataFile*));
    // Loop on the number of frames
    for (int fr=0; fr<mp_ImageDimensionsAndQuantification->GetNbTimeFrames(); fr++)
    {
      // If one normalization file for all frames
      if (m_OneNormalizationFileForAllFrames)
      {
        // Then, if not the first frame
        if (fr!=0)
        {
          // We copy the pointer of the first frame into this one
          m3p_NormalizationDataFile[bed][fr] = m3p_NormalizationDataFile[bed][0];
          // And we continue the loop
          continue;
        }
      }
      // Switch on the scanner type and create the specific data files
      if (p_ScannerManager->GetScannerType() == SCANNER_PET) m3p_NormalizationDataFile[bed][fr] = new iDataFilePET();
      else
      {
        Cerr("***** oSensitivityGenerator::InitializeNormalizationFiles() -> Unknown scanner type (" << p_ScannerManager->GetScannerType() << ") or not yet implemented !" << endl);
        return 1;
      }
      m3p_NormalizationDataFile[bed][fr]->SetHeaderDataFileName(normalization_files[fr]);
      m3p_NormalizationDataFile[bed][fr]->SetBedIndex(bed);
      m3p_NormalizationDataFile[bed][fr]->SetPercentageLoad(m2p_DataFile[0]->GetPercentageLoad());
      m3p_NormalizationDataFile[bed][fr]->SetVerbose(m2p_DataFile[0]->GetVerbose());
      m3p_NormalizationDataFile[bed][fr]->SetImageDimensionsAndQuantification(mp_ImageDimensionsAndQuantification);
      if (m3p_NormalizationDataFile[bed][fr]->ReadInfoInHeader())
      {
        Cerr("***** oSensitivityGenerator::InitializeNormalizationFiles() -> A problem occurred during normalization datafile header reading !" << endl);
        return 1;
      }
      if (m3p_NormalizationDataFile[bed][fr]->CheckParameters())
      {
        Cerr("***** oSensitivityGenerator::InitializeNormalizationFiles() -> A problem occurred while checking normalization datafile parameters !" << endl);
        return 1;
      }
      if (m3p_NormalizationDataFile[bed][fr]->ComputeSizeEvent())
      {
        Cerr("***** oSensitivityGenerator::InitializeNormalizationFiles() -> A problem occurred in normalization datafile event size computation !" << endl);
        return 1;
      }
      if (m3p_NormalizationDataFile[bed][fr]->InitializeFile())
      {
        Cerr("***** oSensitivityGenerator::InitializeNormalizationFiles() -> A problem occurred in normalization datafile initialization !" << endl);
        return 1;
      }
      if (m3p_NormalizationDataFile[bed][fr]->PrepareDataFile())
      {
        Cerr("***** oSensitivityGenerator::InitializeNormalizationFiles() -> A problem occured in normalization datafile preparation !" << endl);
        return 1;
      }
    }
  }

  // End
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int oSensitivityGenerator::InitializeAttenuationFiles()
{
  // If it is asked to ignore attenuation correction, then nothing to do
  if (mp_ImageDimensionsAndQuantification->GetIgnoreAttnCorrectionFlag()) return 0;

  // If attenuation image file name is empty, then we do some checks before
  if (m_pathToAttenuationImage=="")
  {
    // Case for PET: if no normalization data file OR no ACF in the normalization data file
    if (m2p_DataFile[0]->GetDataType()==TYPE_PET && (m3p_NormalizationDataFile==NULL || !(dynamic_cast<iDataFilePET*>(m3p_NormalizationDataFile[0][0]))->GetAtnCorrectionFlag()))
    {
      // Cast the vDataFile into iDataFilePET
      iDataFilePET* p_pet_file = (dynamic_cast<iDataFilePET*>(m2p_DataFile[0]));
      // Throw an error if attenuation correction is in the datafile (at this step we know we do not ignore it)
      if (p_pet_file->GetAtnCorrectionFlag())
      {
        Cerr("***** oSensitivityGenerator::InitializeAttenuationFiles() -> Attenuation correction is included in the data file while it is not in the sensitivity computation !" << endl);
        return 1;
      }
    }
    // We can exit now
    return 0;
  }

  // In PET, if the ACF is already present in the normalization data file, then say it and ignore mumap
  if (m2p_DataFile[0]->GetDataType()==TYPE_PET && m3p_NormalizationDataFile!=NULL && (dynamic_cast<iDataFilePET*>(m3p_NormalizationDataFile[0][0]))->GetAtnCorrectionFlag())
  {
    Cout("oSensitivityGenerator::InitializeAttenuationFiles() -> Ignore provided mumap and consider ACF provided in the normalization data file" << endl);
    return 0;
  }

  // Verbose
  if (m_verbose>=2) Cout("oSensitivityGenerator::InitializeAttenuationFiles() -> Allocate and read attenuation images (assumed to be in cm-1)" << endl);

  // Allocate and read the attenuation image (this puts the image into the m4p_attenuation buffer of oImageSpace)
  if (mp_ImageSpace->InitAttenuationImage(m_pathToAttenuationImage))
  {
    Cerr("***** oSensitivityGenerator::Initialize() -> A problem occured while initializing the attenuation image into the image space !" << endl);  
    return 1;
  }

  // Allocate the foward image, where the attenuation will be copied and deformed if needed
  mp_ImageSpace->LMS_InstantiateForwardImage();

  // Copy attenuation image into forward-image buffer
  mp_ImageSpace->LMS_CopyAtnToForwardImage();

  // We do not need the m4p_attenuation buffer of oImageSpace anymore, so free it now
  mp_ImageSpace->LMS_DeallocateAttenuationImage();

  // TODO: Think about adding a field in the attenuation header to specify if the attenuation images have the same spatial resolution as the
  //       scanner in use, in order to apply a smoothing on it. For the moment, we assume it has been done beforehand.

  // TODO: for the moment, we assume that the attenuation image is in cm-1, maybe use a header field to detect that and convert it
  //       automatically.

  // Put attenuation flag on
  m_attenuationFlag = true;

  // End
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int oSensitivityGenerator::Launch()
{
  // Verbose
  if (m_verbose>=1) Cout("oSensitivityGenerator::Launch() -> Start the sensitivity computation" << endl);
  // Simply switches between the CPU and GPU versions
  #ifdef CASTOR_GPU
  if (m_flagGPU) return LaunchGPU();
  else return LaunchCPU();
  #else
  return LaunchCPU();
  #endif
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int oSensitivityGenerator::LaunchCPU()
{
  // Loop on the bed positions
  for (int bed=0; bed<mp_ImageDimensionsAndQuantification->GetNbBeds(); bed++)
  {
    // Reset line counters
    for (int th=0; th<mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection(); th++) mp_lineCounter[th] = 0;
    // If a normalization file is provided, then use them to iterate on the elements to be taken into account
    if (m3p_NormalizationDataFile!=NULL && ComputeSensitivityFromNormalizationFile(bed))
    {
      Cerr("***** oSensitivityGenerator::LaunchCPU() -> A problem occured while computing the sensitivity using a normalization file for bed " << bed+1 << " !" << endl);
      return 1;
    }
    // Otherwise, use a generic method from the scanner elements
    if (m3p_NormalizationDataFile==NULL && ComputeSensitivityFromScanner(bed))
    {
      Cerr("***** oSensitivityGenerator::LaunchCPU() -> A problem occured while computing the sensitivity from scanner elements for bed " << bed+1 << " !" << endl);
      return 1;
    }
    // Sum up the line counter into the first thread
    for (int th=1; th<mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection(); th++) mp_lineCounter[0] += mp_lineCounter[th];
    // Verbose
    if (m_verbose>=2) Cout("  --> Number of effectively projected lines: " << mp_lineCounter[0] << endl);
  }
  // Deallocate the forward image only if attenuation was used
  if (m_attenuationFlag) mp_ImageSpace->LMS_DeallocateForwardImage();
  // Reduce all results (merging parallel computation)
  for (int fr=0 ; fr<mp_ImageDimensionsAndQuantification->GetNbTimeFrames() ; fr++)
    for (int rg=0 ; rg<mp_ImageDimensionsAndQuantification->GetSensNbRespGates() ; rg++)
      for (int cg=0 ; cg<mp_ImageDimensionsAndQuantification->GetSensNbCardGates() ; cg++)
      {
        // Merge parallel results into the backward image
        mp_ImageSpace->LMS_Reduce(fr, rg, cg);
        // Perform any deformations on the backward image
        mp_DeformationManager->ApplyDeformationForSensitivityGeneration(mp_ImageSpace, BACKWARD_DEFORMATION, 0, fr, rg, cg);
        mp_DeformationManager->ApplyDeformationForSensitivityGeneration(mp_ImageSpace, BACKWARD_DEFORMATION, 1, fr, rg, cg);
        mp_DeformationManager->ApplyDeformationForSensitivityGeneration(mp_ImageSpace, BACKWARD_DEFORMATION, 2, fr, rg, cg);
      }
  // Allocate sensitivity image
  mp_ImageSpace->LMS_InstantiateSensitivityImage(); 
  // Copy backward image to sensitivity image
  mp_ImageSpace->LMS_CopyBackwardToSensitivity();
  // Deallocate backward image
  mp_ImageSpace->LMS_DeallocateBackwardImage();
  // Apply PSF if needed
  mp_ImageConvolverManager->ConvolveSensitivity(mp_ImageSpace);
  // Save sensitivity image
  mp_ImageSpace->LMS_SaveSensitivityImage(m_pathToSensitivityImage, mp_DeformationManager);
  // Deallocate sensitivity image
  mp_ImageSpace->LMS_DeallocateSensitivityImage();
  // End
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
// TODO: implement a function ComputeSensitivityFromDataFile when in histogram mode
int oSensitivityGenerator::ComputeSensitivityFromNormalizationFile(int a_bed)
{
  // Verbose
  if (m_verbose>=2)
  {
    if (mp_ImageDimensionsAndQuantification->GetNbBeds()>1) Cout("oSensitivityGenerator::ComputeSensitivityFromNormalizationFile() -> Start computation for bed " << a_bed+1 << endl);
    else Cout("oSensitivityGenerator::ComputeSensitivityFromNormalizationFile() -> Start computation" << endl);
  }

  // Initial clock
  clock_t clock_start = clock();
  time_t time_start = time(NULL);

  // Loop on frames
  for (int fr=0 ; fr<mp_ImageDimensionsAndQuantification->GetNbTimeFrames() ; fr++)
  {
    // Verbose
    if (m_verbose>=2 && mp_ImageDimensionsAndQuantification->GetNbTimeFrames()>1)
      cout << "  --> Processing frame " << fr+1 << endl;
    // Loop on respiratory gates
    for (int rg=0 ; rg<mp_ImageDimensionsAndQuantification->GetSensNbRespGates() ; rg++)
    {
      // Verbose
      if (m_verbose>=2 && mp_ImageDimensionsAndQuantification->GetSensNbRespGates()>1)
        cout << "  --> Processing respiratory gate " << rg+1 << endl;
      // Loop on cardiac gates
      for (int cg=0 ; cg<mp_ImageDimensionsAndQuantification->GetSensNbCardGates() ; cg++) // TODO: why NbCardGates and not SensNbCardGates here ?
      {
        // Verbose
        if (m_verbose>=2 && mp_ImageDimensionsAndQuantification->GetSensNbCardGates()>1)
          cout << "  --> Processing cardiac gate " << cg+1 << endl;
        // Perform any forward deformations on the forward image (only if attenuation is used)
        if (m_attenuationFlag)
        {
          mp_DeformationManager->ApplyDeformationForSensitivityGeneration(mp_ImageSpace, FORWARD_DEFORMATION, 0, fr, rg, cg);
          mp_DeformationManager->ApplyDeformationForSensitivityGeneration(mp_ImageSpace, FORWARD_DEFORMATION, 1, fr, rg, cg);
          mp_DeformationManager->ApplyDeformationForSensitivityGeneration(mp_ImageSpace, FORWARD_DEFORMATION, 2, fr, rg, cg);
        }
        // Reset the file buffer range before launching the loop over the datafile (this is done only if the percentage load is under 100)
        m3p_NormalizationDataFile[a_bed][fr]->ResetBufferRange();
        // In case a problem occurs in the parallel loop
        bool problem = false;
        // Set the number of threads for projections (right after this loop, we set back the number of threads for image computations)
        #ifdef CASTOR_OMP
        omp_set_num_threads(mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection());
        #endif
        // Compute the index limit over which we consider that threads have finished the loop
        // (to be valid, the omp loop need to use a static scheduling with a chunck size of 1)
        int64_t index_limit = m3p_NormalizationDataFile[a_bed][fr]->GetSize() - 1 - 2*((int64_t)(mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection()));
        // Launch the loop on all events
        int64_t index = 0, printing_index = 0;
        #pragma omp parallel for private(index) schedule(static, 1)
        for (index=0; index<m3p_NormalizationDataFile[a_bed][fr]->GetSize(); index++)
        {              
          // Get the thread index
          int th = 0;
          #ifdef CASTOR_OMP
          th = omp_get_thread_num();
          #endif
          // Verbose
          if (th==0 && m_verbose>=2)
          {
            if (printing_index%1000==0)
            {
              int percent = ((int)( ((FLTNB)(index)) * 100. / ((FLTNB)(m3p_NormalizationDataFile[a_bed][fr]->GetSize())) ));
              cout << "\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b      "
                   << percent << " %                    " << flush;
            }
            printing_index++;
          }
          // Get the current event for that thread index
          vEvent* event = m3p_NormalizationDataFile[a_bed][fr]->GetEventWithAscendingOrderAssumption(index, index_limit, th);
          if (event==NULL)
          {
            Cerr("***** oSensitivityGenerator::ComputeSensitivityFromNormalizationFile() -> An error occured while getting the event from index "
                 << index << " (thread " << th << ") !" << endl);
            // A problem was encountered
            problem = true;
            // We must continue here because we are inside an OpenMP loop
            continue;
          }
          // Compute the projection line
          oProjectionLine *line = mp_ProjectorManager->ComputeProjectionLine(event, th);
          if (line==NULL)
          {
            Cerr("***** oSensitivityGenerator::ComputeSensitivityFromNormalizationFile() -> A problem occured while computing the projection line !" << endl);
            // Specify that there was a problem
            problem = true;
            // We must continue here because we are inside an OpenMP loop
            continue;
          }
          // Process this line
          if (line->NotEmptyLine() && ProcessThisLine(line, event, a_bed, fr, rg, cg, th))
          {
            Cerr("***** oSensitivityGenerator::ComputeSensitivityFromNormalizationFile() -> A problem occured while processing a line !" << endl);
          }
        }
        // End of progression printing (do not log out with Cout here)
        if (m_verbose>=2) cout << "\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b"
                               << "  --> 100 %                        " << endl;
        // If a problem was encountered, then report it here
        if (problem)
        {
          Cerr("***** oSensitivityGenerator::ComputeSensitivityFromNormalizationFile() -> A problem occured inside the parallel loop over events !" << endl);
          return 1;
        }
      }
    }
  }

  // Clock total
  clock_t clock_stop = clock();
  time_t time_stop = time(NULL);
  if (mp_ImageDimensionsAndQuantification->GetNbBeds()>1)
    Cout("  --> Time spent for sensitivity generation for bed " << a_bed+1 << " | User: " << time_stop-time_start
         << " sec | CPU: " << (clock_stop-clock_start)/((FLTNB)CLOCKS_PER_SEC) << " sec" << endl);
  else
    Cout("  --> Time spent for sensitivity generation | User: " << time_stop-time_start
         << " sec | CPU: " << (clock_stop-clock_start)/((FLTNB)CLOCKS_PER_SEC) << " sec" << endl);

  // End
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int oSensitivityGenerator::ComputeSensitivityFromScanner(int a_bed)
{
  // Verbose
  if (m_verbose>=2)
  {
    if (mp_ImageDimensionsAndQuantification->GetNbBeds()>1) Cout("oSensitivityGenerator::ComputeSensitivityFromScanner() -> Start computation for bed " << a_bed+1 << endl);
    else Cout("oSensitivityGenerator::ComputeSensitivityFromScanner() -> Start computation" << endl);
  }

  // Get the scanner manager
  sScannerManager* p_scannerManager = sScannerManager::GetInstance();

  // Total number of elements
  int nb_total_elts = mp_Scanner->GetSystemNbElts();

  // Initial clock
  clock_t clock_start = clock();
  time_t time_start = time(NULL);

  // Initialize main loop start and stop values
  int64_t main_loop_start_index = 0 ;
  int64_t main_loop_stop_index = 0;

  // Check beforehand any issue with the loop start/stop values (not possible in the inner multithreaded loop)
  if (p_scannerManager->PROJ_GetModalityStopValueMainLoop()>0) main_loop_stop_index = p_scannerManager->PROJ_GetModalityStopValueMainLoop();
  else
  {
    Cerr("***** oAnalyticProjection::LaunchCPU()-> An error occured when trying to initialize main loop stop index !" << endl); 
    return 1;
  }
  if (p_scannerManager->PROJ_GetModalityStartValueInnerLoop(0)<=0)
  {
    Cerr("***** oAnalyticProjection::LaunchCPU()-> An error occured when trying to initialize inner loop start index !" << endl); 
    return 1;
  }

  // Prepare pre-computed sums of events to avoid the exponential evolution of the percentage (for PET)
  // TODO Perhaps replace this with a call to scannerManager for error management
  int64_t* progression_elts_array = new int64_t[nb_total_elts*sizeof(int64_t)];
  progression_elts_array[0] = 0;
  for (int idx_elt=1 ; idx_elt<mp_Scanner->GetSystemNbElts() ; idx_elt++) 
    progression_elts_array[idx_elt] = progression_elts_array[idx_elt-1] + (mp_Scanner->GetSystemNbElts()-idx_elt);

  // Index for progression printing
  uint64_t printing_index = 0;
  uint64_t progression_index_total = p_scannerManager->PROJ_GetProgressionFinalValue()*
                                              mp_ImageDimensionsAndQuantification->GetNbTimeFrames()*
                                              mp_ImageDimensionsAndQuantification->GetSensNbRespGates()*
                                              mp_ImageDimensionsAndQuantification->GetSensNbCardGates();
  // Loop on frames
  for (int fr=0 ; fr<mp_ImageDimensionsAndQuantification->GetNbTimeFrames() ; fr++)
  {
    // Verbose
    if (m_verbose>=2 && mp_ImageDimensionsAndQuantification->GetNbTimeFrames()>1)
      cout << "  --> Processing frame " << fr+1 << endl;
    // Loop on respiratory gates
    for (int rg=0 ; rg<mp_ImageDimensionsAndQuantification->GetSensNbRespGates() ; rg++)
    {
      // Verbose
      if (m_verbose>=2 && mp_ImageDimensionsAndQuantification->GetSensNbRespGates()>1)
        cout << "  --> Processing respiratory gate " << rg+1 << endl;
      // Loop on cardiac gates
      for (int cg=0 ; cg<mp_ImageDimensionsAndQuantification->GetSensNbCardGates() ; cg++)
      {
        // Verbose
        if (m_verbose>=2 && mp_ImageDimensionsAndQuantification->GetSensNbCardGates()>1)
          cout << "  --> Processing cardiac gate " << cg+1 << endl;
        // Perform any forward deformations on the forward image (only if attenuation is used)
        if (m_attenuationFlag)
        {
          mp_DeformationManager->ApplyDeformationForSensitivityGeneration(mp_ImageSpace, FORWARD_DEFORMATION, 0, fr, rg, cg);
          mp_DeformationManager->ApplyDeformationForSensitivityGeneration(mp_ImageSpace, FORWARD_DEFORMATION, 1, fr, rg, cg);
          mp_DeformationManager->ApplyDeformationForSensitivityGeneration(mp_ImageSpace, FORWARD_DEFORMATION, 2, fr, rg, cg);
        }
        // In case a problem occurs in the parallel loop
        bool problem = false;
        // Set the number of threads for projections (right after this loop, we set back the number of threads for image computations)
        #ifdef CASTOR_OMP
        omp_set_num_threads(mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection());
        #endif

        // Start the loop
        int64_t idx_elt1;

        #pragma omp parallel for private(idx_elt1) schedule(static, 1)
        for(idx_elt1=main_loop_start_index ; idx_elt1<main_loop_stop_index ; idx_elt1++)
        {
          // Get the thread index
          int th = 0;
          #ifdef CASTOR_OMP
          th = omp_get_thread_num();
          #endif

          // Initialize inner loop start and stop values
          int64_t inner_loop_start_index = p_scannerManager->PROJ_GetModalityStartValueInnerLoop(idx_elt1) ;
          int64_t inner_loop_stop_index = mp_Scanner->GetSystemNbElts();
          
          // Inner loop on scanner elements (idx_elt2) which are used to form LOR using the first scanner element (idx_elt1) 
          for (int64_t idx_elt2=inner_loop_start_index ; idx_elt2<inner_loop_stop_index ; idx_elt2++)
          {
            // Print progression (do not log out with Cout here)
            if (th==0)
            {
              if (m_verbose>=2 && printing_index%10000==0)
              {
               int64_t progression_index_current = p_scannerManager->PROJ_GetCurrentProgression(idx_elt1, idx_elt2, progression_elts_array, 
                                                                                                  mp_ImageDimensionsAndQuantification->GetSensNbRespGates(), 
                                                                                                  mp_ImageDimensionsAndQuantification->GetSensNbCardGates(), 
                                                                                                                                fr, rg, cg);
                int64_t percent = (int64_t) (( ((FLTNB)progression_index_current)/((FLTNB)progression_index_total) ) * ((FLTNB)100));
                cout << "\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b      "
                   << percent << " %                       ";
              }
              printing_index++;
            }

            // Allocate an event using the iDataFile
            vEvent* event = m2p_DataFile[a_bed]->PROJ_GenerateEvent(idx_elt1, idx_elt2, th);

            // Check from the scanner requirements that this LOR is allowed
            if (!mp_Scanner->IsAvailableLOR(event->GetID1(0), event->GetID2(0))) continue;

            // Generate the projection event and compute the projection line 
            oProjectionLine* line = mp_ProjectorManager->ComputeProjectionLine(event, th);
            if (line==NULL)
            {
              Cerr("***** oSensitivityGenerator::ComputeSensitivityFromScanner() -> A problem occured while computing the projection line !" << endl);
              // Specify that there was a problem
              problem = true;
              // We must continue here because we are inside an OpenMP loop
              continue;
            }

            // Process this line
            if (line->NotEmptyLine() && ProcessThisLine(line, event, a_bed, fr, rg, cg, th))
            {
              Cerr("***** oSensitivityGenerator::ComputeSensitivityFromScanner() -> A problem occured while processing a line !" << endl);
            }
          }
        }
        // End of progression printing (do not log out with Cout here)
        if (m_verbose>=2) cout << "\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b"
                               << "  --> 100 %                       " << endl;
        // If a problem was encountered, then report it here
        if (problem)
        {
          Cerr("***** oSensitivityGenerator::ComputeSensitivityFromScanner() -> A problem occured inside the parallel loop over events !" << endl);
          return 1;
        }
      }
    }
  }
  delete[] progression_elts_array;

  // Clock total
  clock_t clock_stop = clock();
  time_t time_stop = time(NULL);
  if (mp_ImageDimensionsAndQuantification->GetNbBeds()>1)
    Cout("  --> Time spent for sensitivity generation for bed " << a_bed+1 << " | User: " << time_stop-time_start
         << " sec | CPU: " << (clock_stop-clock_start)/((FLTNB)CLOCKS_PER_SEC) << " sec" << endl);
  else
    Cout("  --> Time spent for sensitivity generation | User: " << time_stop-time_start
         << " sec | CPU: " << (clock_stop-clock_start)/((FLTNB)CLOCKS_PER_SEC) << " sec" << endl);

  // End
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int oSensitivityGenerator::ProcessThisLine(oProjectionLine* ap_Line, vEvent* ap_Event, int a_bed, int a_frame, int a_respGate, int a_cardGate, int a_thread)
{
  // ----------------------------------------------------------------------------------------------
  // Part 0: get the multiplicative corrections from the vEvent and the quantification factor
  // ----------------------------------------------------------------------------------------------

  // The multiplicative corrections from the vEvent can be:
  //   1: generated from vDataFile::GenerateEvent(), this will be a default value, so this will be 1.
  //   2: taken from the normalization data file, so this will be the normalization correction factor times the attenuation correction factor
  FLTNB multiplicative_factor = ap_Event->GetMultiplicativeCorrections();
  // If null, then skip this event
  if (multiplicative_factor<=0.) return 0;
  // Get the quantification factor
  FLTNB quantification_factor = mp_ImageDimensionsAndQuantification->GetQuantificationFactor(a_bed, a_frame, a_respGate, a_cardGate);
  // If null, then skip this event
  if (quantification_factor<=0.) return 0;

  // ----------------------------------------------------------------------------------------------
  // Part 1: deal with the attenuation forward projection for PET if any
  // ----------------------------------------------------------------------------------------------

  // The ACF
  FLTNB attenuation_correction_factor = 1.;

  // Do a forward projection only if the attenuation image was provided and the ACF was not included in the normalization data file
  if ( m_forwardProjectAttenuation )
  {
    // Cumulative mu (in cm-1)
    FLTNB cumulative_mu = 0.;
    // Loop on TOF bins
    for (int b=0; b<ap_Line->GetNbTOFBins(); b++)
    {
      // Projection operation (the attenuation image in the m4p_forwardImage buffer of the oImageSpace)
      for (int vl=0 ; vl<ap_Line->GetCurrentNbVoxels(FORWARD, b) ; vl++)
        cumulative_mu += ap_Line->GetVoxelWeights(FORWARD, b, vl) * mp_ImageSpace->m4p_forwardImage[a_frame][a_respGate][a_cardGate][ap_Line->GetVoxelIndex(FORWARD, b,vl)];
    }
    // Update the ACF (converting the cumulative mu in mm-1)
    attenuation_correction_factor = max(((FLTNB)1.),exp(cumulative_mu*((FLTNB)0.1)));
  }

  // ----------------------------------------------------------------------------------------------
  // Part 2: back-projection of the LOR sensitivity
  // ----------------------------------------------------------------------------------------------

  // Compute the global sensitivity (calibration / (ACF * norm))
  FLTNB lor_sensitivity = 1. / (quantification_factor * attenuation_correction_factor * multiplicative_factor);

  // Loop on TOF bins
  for (int b=0; b<ap_Line->GetNbTOFBins(); b++)
  {
    // Backprojection operation into the backward image
    for (int vl=0 ; vl<ap_Line->GetCurrentNbVoxels(BACKWARD, b) ; vl++)
      mp_ImageSpace->m6p_backwardImage[0][a_thread][a_frame][a_respGate][a_cardGate][ap_Line->GetVoxelIndex(BACKWARD, b, vl)] += ap_Line->GetVoxelWeights(BACKWARD, b, vl)
                                                                                                                               * lor_sensitivity;
  }

  // Increment line counters
  mp_lineCounter[a_thread]++;

  // End
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

string oSensitivityGenerator::GetPathToSensitivityImage() 
{
  // Get the flag that says if we merge the dynamic files or not
  bool merge_dynamic_files = sOutputManager::GetInstance()->MergeDynImages();
  // If we merge dynamic file, or we are in a static case (no dynamic), we stay with the hdr extension
  if ( ( mp_ImageDimensionsAndQuantification->GetSensNbCardGates()==1 &&
         mp_ImageDimensionsAndQuantification->GetSensNbRespGates()==1 &&
         mp_ImageDimensionsAndQuantification->GetNbTimeFrames()==1 ) ||
       merge_dynamic_files )
    return m_pathToSensitivityImage + ".hdr";
  else
    return m_pathToSensitivityImage + ".mhd";
}