oProjectorManager.cc

Go to the documentation of this file.

/*
  Implementation of class oProjectorManager

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

#include "oProjectorManager.hh"
#include "sOutputManager.hh"
#include "sAddonManager.hh"
#include "iDataFilePET.hh"

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

oProjectorManager::oProjectorManager()
{
  // Scanner and image dimensions
  mp_Scanner = NULL;
  mp_ImageDimensionsAndQuantification = NULL;
  // Data file (used to get some info about TOF bins and POI)
  mp_DataFile = NULL;
  // TOF and POI options
  m_applyTOF = false;
  m_applyPOI = false;
  m_nbTOFBins = -1;
  // Computation strategy for projection lines
  m_computationStrategy = -1;
  // Forward and backward options for projectors
  m_optionsForward = "";
  m_optionsBackward = "";
  // Common options
  m_optionsCommon = "";
  // Forward and backward projectors
  mp_SystemMatrixForward = NULL;
  mp_SystemMatrixBackward = NULL;
  mp_ProjectorForward = NULL;
  mp_ProjectorBackward = NULL;
  m_UseSystemMatrixForward = false;
  m_UseSystemMatrixBackward = false;
  m_UseProjectorForward = false;
  m_UseProjectorBackward = false;
  m_UseMatchedProjectors = false;
  // Forward and backward projection lines (as many as threads)
  m2p_ProjectionLines = NULL;
  // Verbosity
  m_verbose = 0;
  // Not checked yet
  m_checked = false;
  // Not initialized yet
  m_initialized = false;
}

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


oProjectorManager::~oProjectorManager() 
{
  // Go through the destructor only if the object was initialized
  if (m_initialized)
  {
    // Delete projection lines
    for (int th=0; th<mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection(); th++)
      if (m2p_ProjectionLines[th]) delete m2p_ProjectionLines[th];
    if (m2p_ProjectionLines) delete[] m2p_ProjectionLines;
    // Delete projectors
    if (m_UseProjectorForward) delete mp_ProjectorForward;
    if (m_UseSystemMatrixForward) delete mp_SystemMatrixForward;
    if (!m_UseMatchedProjectors)
    {
      if (m_UseProjectorBackward) delete mp_ProjectorBackward;
      if (m_UseSystemMatrixBackward) delete mp_SystemMatrixBackward;
    }
  }
}

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

bool oProjectorManager::IsForwardOperatorCompatibleWithSPECTAttenuationCorrection()
{
  // Case for a vProjector
  if (m_UseProjectorForward) return mp_ProjectorForward->GetCompatibilityWithSPECTAttenuationCorrection();
  // Case for a system matrix
  else return mp_SystemMatrixForward->GetCompatibilityWithSPECTAttenuationCorrection();
}

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

bool oProjectorManager::IsBackwardOperatorCompatibleWithSPECTAttenuationCorrection()
{
  // Case for a vProjector
  if (m_UseProjectorBackward) return mp_ProjectorBackward->GetCompatibilityWithSPECTAttenuationCorrection();
  // Case for a system matrix
  else return mp_SystemMatrixBackward->GetCompatibilityWithSPECTAttenuationCorrection();
}

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

int oProjectorManager::CheckParameters()
{
  // Check scanner
  if (mp_Scanner==NULL)
  {
    Cerr("***** oProjectorManager::CheckParameters() -> No scanner provided !" << endl);
    return 1;
  }
  // Check image dimensions
  if (mp_ImageDimensionsAndQuantification==NULL)
  {
    Cerr("***** oProjectorManager::CheckParameters() -> No image dimensions provided !" << endl);
    return 1;
  }
  // Check data file
  if (mp_DataFile==NULL)
  {
    Cerr("***** oProjectorManager::CheckParameters() -> No data file provided !" << endl);
    return 1;
  }
  // Check computation strategy
  if ( m_computationStrategy != IMAGE_COMPUTATION_STRATEGY &&
       m_computationStrategy != FIXED_LIST_COMPUTATION_STRATEGY &&
       m_computationStrategy != ADAPTATIVE_LIST_COMPUTATION_STRATEGY )
  {
    Cerr("***** oProjectorManager::CheckParameters() -> Unknown computation strategy provided !" << endl);
    return 1;
  }
  // Check forward projector options
  if (m_optionsForward=="")
  {
    Cerr("***** oProjectorManager::CheckParameters() -> No forward projector options provided !" << endl);
    return 1;
  }
  // Check backward projector options
  if (m_optionsBackward=="")
  {
    Cerr("***** oProjectorManager::CheckParameters() -> No backward projector options provided !" << endl);
    return 1;
  }
  // Check verbosity
  if (m_verbose<0)
  {
    Cerr("***** oProjectorManager::CheckParameters() -> Wrong verbosity level provided !" << endl);
    return 1;
  }
  // Normal end
  m_checked = true;
  return 0;
}

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

int oProjectorManager::Initialize()
{
  // Forbid initialization without check
  if (!m_checked)
  {
    Cerr("***** oProjectorManager::Initialize() -> Must call CheckParameters() before Initialize() !" << endl);
    return 1;
  }

  // Verbose
  if (m_verbose>=1) Cout("oProjectorManager::Initialize() -> Initialize projectors and projection lines" << endl);

  // TOF-PET specific stuff
  FLTNB tof_resolution = -1.;
  if (mp_DataFile->GetDataType()==TYPE_PET)
  {
    // Cast the datafile pointer
    iDataFilePET* p_pet_datafile = (dynamic_cast<iDataFilePET*>(mp_DataFile));
    // Case of TOF information in the datafile
    if (p_pet_datafile->GetTOFInfoFlag())
    {
      // Case where it is asked to ignore TOF info
      if (p_pet_datafile->GetIgnoreTOFFlag())
      {
        m_applyTOF = false;
        m_nbTOFBins = 1;
      }
      // Otherwise, we use TOF
      else
      {
        // TODO: set tof bin size
        // Get the number of TOF bins from the data file
        m_nbTOFBins = p_pet_datafile->GetNbBinsTOF();
        // Get the TOF resolution from the data file, if we apply TOF
        tof_resolution = p_pet_datafile->GetResolutionTOF();
      }
    }
    // Case no TOF
    else
    {
      m_applyTOF = false;
      m_nbTOFBins = 1;
    }
  }
  else m_nbTOFBins = 1;

  // Look if we must apply POI or not
  if (mp_DataFile->GetPOIInfoFlag() && !mp_DataFile->GetIgnorePOIFlag()) m_applyPOI = true;
  else m_applyPOI = false;
  // Get POI resolution from the data file
  FLTNB* poi_resolution = mp_DataFile->GetResolutionPOI();

  // Compare projector options to know if we use matched ones for forward and backward operations
  if (m_optionsForward==m_optionsBackward) m_UseMatchedProjectors = true;
  else m_UseMatchedProjectors = false;

  // Parse projector options and initialize them
  if (ParseOptionsAndInitializeProjectors())
  {
    Cerr("***** oProjectorManager::Initialize() -> A problem occured while parsing projector options and initializing it !" << endl);
    return 1;
  }

  // Initialize as many projection lines as threads
  m2p_ProjectionLines = new oProjectionLine*[mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection()];
  for (int th=0; th<mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection(); th++)
  {
    m2p_ProjectionLines[th] = new oProjectionLine();
    m2p_ProjectionLines[th]->SetThreadNumber(th);
    m2p_ProjectionLines[th]->SetMatchedProjectors(m_UseMatchedProjectors);
    m2p_ProjectionLines[th]->SetNbTOFBins(m_nbTOFBins);
    m2p_ProjectionLines[th]->SetTOFResolution(tof_resolution);
    m2p_ProjectionLines[th]->SetPOIResolution(poi_resolution);
    m2p_ProjectionLines[th]->SetImageDimensionsAndQuantification(mp_ImageDimensionsAndQuantification);
    m2p_ProjectionLines[th]->SetComputationStrategy(m_computationStrategy);
    m2p_ProjectionLines[th]->SetForwardProjector(mp_ProjectorForward);
    m2p_ProjectionLines[th]->SetBackwardProjector(mp_ProjectorBackward);
    m2p_ProjectionLines[th]->SetVerbose(m_verbose);
    if (m2p_ProjectionLines[th]->CheckParameters())
    {
      Cerr("***** oProjectorManager::Initialize() -> An error occured while checking parameters of an oProjectionLine !" << endl);
      return 1;
    }
    if (m2p_ProjectionLines[th]->Initialize())
    {
      Cerr("***** oProjectorManager::Initialize() -> An error occured while initializing an oProjectionLine !" << endl);
      return 1;
    }
  }

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

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

int oProjectorManager::ParseOptionsAndInitializeProjectors()
{
  string projector = "";
  string list_options = "";
  string file_options = "";

  // This is for the automatic initialization of the projectors
  typedef vProjector *(*maker_projector) ();

  // ---------------------------------------------------------------------------------------------------
  // Manage forward projector
  // ---------------------------------------------------------------------------------------------------

  // ______________________________________________________________________________
  // Get the projector name in the options and isolate the real projector's options

  // Search for a colon ":", this indicates that a configuration file is provided after the projector name
  size_t colon = m_optionsForward.find_first_of(":");
  size_t comma = m_optionsForward.find_first_of(",");

  // Case 1: we have a colon
  if (colon!=string::npos)
  {
    // Get the projector name before the colon
    projector = m_optionsForward.substr(0,colon);
    // Get the configuration file after the colon
    file_options = m_optionsForward.substr(colon+1);
    // List of options is empty
    list_options = "";
  }
  // Case 2: we have a comma
  else if (comma!=string::npos)
  {
    // Get the projector name before the first comma
    projector = m_optionsForward.substr(0,comma);
    // Get the list of options after the first comma
    list_options = m_optionsForward.substr(comma+1);
    // Configuration file is empty
    file_options = "";
  }
  // Case 3: no colon and no comma (a single projector name)
  else
  {
    // Get the projector name
    projector = m_optionsForward;
    // List of options is empty
    list_options = "";
    // Build the default configuration file
    file_options = sOutputManager::GetInstance()->GetPathToConfigDir() + "/projector/" + projector + ".conf";
  }

  // ______________________________________________________________________________
  // Case 1: projector is equal to keyword 'matrix', then use a system matrix
  if (projector==SYSTEM_MATRIX_KEYWORD)
  {
    mp_ProjectorForward = NULL;
    m_UseProjectorForward = false;
    m_UseSystemMatrixForward = true;
    // TODO: forbid TOF in PET with system matrix
    Cerr("***** oProjectorManager::ParseOptionsAndInitializeProjectors() -> Loading of custom system matrices is not yet implemented !" << endl);
    return 1;
  }
  // ______________________________________________________________________________
  // Case 2: on-the-fly projector
  else
  {
    // Unset system matrix
    mp_SystemMatrixForward = NULL;
    m_UseSystemMatrixForward = false;
    // Set projector on
    m_UseProjectorForward = true;
    // Get projector's listfrom addon manager
    std::map <string,maker_projector> list = sAddonManager::GetInstance()->mp_listOfProjectors;
    // Create the projector
    if (list[projector]) mp_ProjectorForward = list[projector]();
    else
    {
      Cerr("***** oProjectorManager::ParseOptionsAndInitializeProjectors() -> Projector '" << projector << "' does not exist !" << endl);
      return 1;
    }
    // Set parameters
    mp_ProjectorForward->SetScanner(mp_Scanner);
    if (mp_ProjectorForward->SetImageDimensionsAndQuantification(mp_ImageDimensionsAndQuantification))
    {
      Cerr("***** oProjectorManager::ParseOptionsAndInitializeProjectors() -> A problem occured while setting the image dimensions of the forward projector !" << endl);
      return 1;
    }
    if (mp_ProjectorForward->SetTOFAndPOIOptions(mp_DataFile->GetDataMode(),m_applyTOF,m_applyPOI))
    {
      Cerr("***** oProjectorManager::ParseOptionsAndInitializeProjectors() -> A problem occured while setting the TOF and POI options of the forward projector !" << endl);
      return 1;
    }
    mp_ProjectorForward->SetVerbose(m_verbose);
    // Provide common options list
    if (mp_ProjectorForward->ReadCommonOptionsList(m_optionsCommon))
    {
      Cerr("***** oProjectorManager::ParseOptionsAndInitializeProjectors() -> A problem occured while parsing and reading forward projector's common options !" << endl);
      return 1;
    }
    // Provide configuration file if any
    if (file_options!="" && mp_ProjectorForward->ReadConfigurationFile(file_options))
    {
      Cerr("***** oProjectorManager::ParseOptionsAndInitializeProjectors() -> A problem occured while parsing and reading forward projector's configuration file !" << endl);
      return 1;
    }
    // Provide options if any
    if (list_options!="" && mp_ProjectorForward->ReadOptionsList(list_options))
    {
      Cerr("***** oProjectorManager::ParseOptionsAndInitializeProjectors() -> A problem occured while parsing and reading forward projector's options !" << endl);
      return 1;
    }
    // Check parameters
    if (mp_ProjectorForward->CheckParameters())
    {
      Cerr("***** oProjectorManager::ParseOptionsAndInitializeProjectors() -> A problem occured while checking forward projector parameters !" << endl);
      return 1;
    }
    // Initialize the projector
    if (mp_ProjectorForward->Initialize())
    {
      Cerr("***** oProjectorManager::ParseOptionsAndInitializeProjectors() -> A problem occured while initializing the forward projector !" << endl);
      return 1;
    }
  }

  // ---------------------------------------------------------------------------------------------------
  // Manage backward projector
  // ---------------------------------------------------------------------------------------------------

  // If options are the same, then forward and backward are the same
  if (m_UseMatchedProjectors)
  {
    // In this case, matched projectors
    m_UseSystemMatrixBackward = m_UseSystemMatrixForward;
    m_UseProjectorBackward = m_UseProjectorForward;
    mp_SystemMatrixBackward = mp_SystemMatrixForward;
    mp_ProjectorBackward = mp_ProjectorForward;
  }
  // Else, unmatched projectors
  else
  {
    // ______________________________________________________________________________
    // Get the projector name in the options and isolate the real projector's options

    // Search for a colon ":", this indicates that a configuration file is provided after the projector name
    colon = m_optionsBackward.find_first_of(":");
    comma = m_optionsBackward.find_first_of(",");

    // Case 1: we have a colon
    if (colon!=string::npos)
    {
      // Get the projector name before the colon
      projector = m_optionsBackward.substr(0,colon);
      // Get the configuration file after the colon
      file_options = m_optionsBackward.substr(colon+1);
      // List of options is empty
      list_options = "";
    }
    // Case 2: we have a comma
    else if (comma!=string::npos)
    {
      // Get the projector name before the first comma
      projector = m_optionsBackward.substr(0,comma);
      // Get the list of options after the first comma
      list_options = m_optionsBackward.substr(comma+1);
      // Configuration file is empty
      file_options = "";
    }
    // Case 3: no colon and no comma (a single projector name)
    else
    {
      // Get the projector name
      projector = m_optionsBackward;
      // List of options is empty
      list_options = "";
      // Build the default configuration file
      file_options = sOutputManager::GetInstance()->GetPathToConfigDir() + "/projector/" + projector + ".conf";
    }

    // ______________________________________________________________________________
    // Case 1: projector is equal to keyword 'matrix', then use a system matrix
    if (projector==SYSTEM_MATRIX_KEYWORD)
    {
      mp_ProjectorBackward = NULL;
      m_UseProjectorBackward = false;
      m_UseSystemMatrixBackward = true;
      // TODO: forbid TOF in PET with system matrix
      Cerr("***** oProjectorManager::ParseOptionsAndInitializeProjectors() -> Loading of custom system matrices is not yet implemented !" << endl);
      return 1;
    }
    // ______________________________________________________________________________
    // Case 2: on-the-fly projector
    else
    {
      // Unset system matrix
      mp_SystemMatrixBackward = NULL;
      m_UseSystemMatrixBackward = false;
      // Set projector on
      m_UseProjectorBackward = true;
      // Get projector's listfrom addon manager
      std::map <string,maker_projector> list = sAddonManager::GetInstance()->mp_listOfProjectors;
      // Create the projector
      if (list[projector]) mp_ProjectorBackward = list[projector]();
      else
      {
        Cerr("***** oProjectorManager::ParseOptionsAndInitializeProjectors() -> Projector '" << projector << "' does not exist !" << endl);
        return 1;
      }
      // Set parameters
      mp_ProjectorBackward->SetScanner(mp_Scanner);
      if (mp_ProjectorBackward->SetImageDimensionsAndQuantification(mp_ImageDimensionsAndQuantification))
      {
        Cerr("***** oProjectorManager::ParseOptionsAndInitializeProjectors() -> A problem occured while setting the image dimensions of the backward projector !" << endl);
        return 1;
      }
      if (mp_ProjectorBackward->SetTOFAndPOIOptions(mp_DataFile->GetDataMode(),m_applyTOF,m_applyPOI))
      {
        Cerr("***** oProjectorManager::ParseOptionsAndInitializeProjectors() -> A problem occured while setting the TOF and POI of the backward projector !" << endl);
        return 1;
      }
      mp_ProjectorBackward->SetVerbose(m_verbose);
      // Provide common options list
      if (mp_ProjectorBackward->ReadCommonOptionsList(m_optionsCommon))
      {
        Cerr("***** oProjectorManager::ParseOptionsAndInitializeProjectors() -> A problem occured while parsing and reading backward projector's common options !" << endl);
        return 1;
      }
      // Provide configuration file if any
      if (file_options!="" && mp_ProjectorBackward->ReadConfigurationFile(file_options))
      {
        Cerr("***** oProjectorManager::ParseOptionsAndInitializeProjectors() -> A problem occured while parsing and reading backward projector's configuration file !" << endl);
        return 1;
      }
      // Provide options if any
      if (list_options!="" && mp_ProjectorBackward->ReadOptionsList(list_options))
      {
        Cerr("***** oProjectorManager::ParseOptionsAndInitializeProjectors() -> A problem occured while parsing and reading backward projector's options !" << endl);
        return 1;
      }
      // Check parameters
      if (mp_ProjectorBackward->CheckParameters())
      {
        Cerr("***** oProjectorManager::ParseOptionsAndInitializeProjectors() -> A problem occured while checking backward projector parameters !" << endl);
        return 1;
      }
      // Initialize the projector
      if (mp_ProjectorBackward->Initialize())
      {
        Cerr("***** oProjectorManager::ParseOptionsAndInitializeProjectors() -> A problem occured while initializing backward projector !" << endl);
        return 1;
      }
    }
  }

  // Normal end
  return 0;
}

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

oProjectionLine* oProjectorManager::ComputeProjectionLine(vEvent* ap_Event, int a_th) 
{
  #ifdef CASTOR_DEBUG
  if (!m_initialized)
  {
    Cerr("***** oProjectorManager::ComputeProjectionLine() -> Called while not initialized !" << endl);
    Exit(EXIT_DEBUG);
  }
  #endif

  #ifdef CASTOR_VERBOSE
  if (m_verbose>=10) Cout("oProjectorManager::ComputeProjectionLine() -> Compute a projection line" << endl);
  #endif

  // Get the list of indices
  uint32_t *index1 = ap_Event->GetEventID1();
  uint32_t *index2 = ap_Event->GetEventID2();
  int nb_indices = ap_Event->GetNbLines();

  // Clean the projection line
  m2p_ProjectionLines[a_th]->Reset();

  // With list-mode data, we may need POI and/or TOF measurements
  if (ap_Event->GetDataMode()==MODE_LIST)
  {
    // Set POI measurement
    if (m_applyPOI)
    {
      // For PET
      if (ap_Event->GetDataType()==TYPE_PET)
      {
        // Have to dynamic_cast the event into a iEventPETList to access the GetPOI functions
        m2p_ProjectionLines[a_th]->SetPOI1((dynamic_cast<iEventListPET*>(ap_Event))->GetPOI1());
        m2p_ProjectionLines[a_th]->SetPOI2((dynamic_cast<iEventListPET*>(ap_Event))->GetPOI2());
      }
      // For SPECT
      if (ap_Event->GetDataType()==TYPE_SPECT)
      {
        // By convention in SPECT, the second end point is in the camera (the first one being outside)
        //m2p_ProjectionLines[a_th]->SetPOI2((dynamic_cast<iEventListModeSPECT*>(ap_Event))->GetPOI());
      }
      // For transmission
      if (ap_Event->GetDataType()==TYPE_TRANSMISSION)
      {
        ;
      }
    }
    // Set TOF measurement (only for PET obviously)
    if (m_applyTOF && ap_Event->GetDataType()==TYPE_PET)
    {
      // Have to dynamic_cast the event into a iEventPETList to access the GetTOF function
      m2p_ProjectionLines[a_th]->SetTOFMeasurement((dynamic_cast<iEventListPET*>(ap_Event))->GetTOFMeasurement());
    }
  }

  // Project forward (and also compute line length)
  int return_value = 0;
  if (m_UseProjectorForward)
    return_value = mp_ProjectorForward->Project( FORWARD, m2p_ProjectionLines[a_th], index1, index2, nb_indices );
  else if (m_UseSystemMatrixForward)
    return_value = mp_SystemMatrixForward->Project( FORWARD, m2p_ProjectionLines[a_th], index1, index2, nb_indices );
  if (return_value)
  {
    Cerr("***** oProjectorManager::ComputeProjectionLine() -> A problem occured while forward projecting an event !" << endl);
    return NULL;
  }

  // Project backward
  if (!m_UseMatchedProjectors)
  {
    // Then project
    if (m_UseProjectorBackward)
      return_value = mp_ProjectorBackward->Project( BACKWARD, m2p_ProjectionLines[a_th], index1, index2, nb_indices );
    else if (m_UseSystemMatrixBackward)
      return_value = mp_SystemMatrixBackward->Project( BACKWARD, m2p_ProjectionLines[a_th], index1, index2, nb_indices);
    if (return_value)
    {
      Cerr("***** oProjectorManager::ComputeProjectionLine() -> A problem occured while backward projecting an event !" << endl);
      return NULL;
    }
  }

  // Return the line
  return m2p_ProjectionLines[a_th];
}

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