iOptimizerMLMuMap.cc

Go to the documentation of this file.

#include "iOptimizerMLMuMap.hh"
#include "sOutputManager.hh"
#include "iDataFilePET.hh"

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

iOptimizerMLMuMap::iOptimizerMLMuMap() : vOptimizer()
{
  // ---------------------------
  // Mandatory member parameters
  // ---------------------------

  // Initial value at 1
  m_initialValue = 1.;
  // Only one backward image for MLMuMap
  m_nbBackwardImages = 1;
  // MLMuMap is compatible only with histogram data
  m_listmodeCompatibility = false;
  m_histogramCompatibility = true;
  // MLMuMap is compatible only with log-converted ACF from PET histogram data
  m_emissionCompatibility = true;
  m_transmissionCompatibility = false;

  // --------------------------
  // Specific member parameters
  // --------------------------

  m_dataSpaceDenominatorThreshold = -1.;
  m_minimumImageUpdateFactor = -1.;
  m_maximumImageUpdateFactor = -1.;
}

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

iOptimizerMLMuMap::~iOptimizerMLMuMap()
{
}

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

void iOptimizerMLMuMap::ShowHelpSpecific()
{
  cout << "This optimizer is a modified MLEM algorithm where the image voxel values are attenuation factors" << endl;
  cout << "and the data are taken as the logarithm of ACF values associated to PET histogram data." << endl;
  cout << "It is numerically implemented in the multiplicative form (as opposed to the gradient form)." << endl;
  cout << "It truncates ACF lower than 1 to 1 to satisfy the positivity constraint." << endl;
  cout << "Note that although physicists are used to attenuation values in 1/cm, the reconstructed image units" << endl;
  cout << "are 1/mm as the whole castor code is using mm." << endl;
  cout << "Subsets can be used." << endl;
  cout << "The following options can be used (in this particular order when provided as a list):" << endl;
  cout << "  initial image value: to set the uniform voxel value for the initial image" << endl;
  cout << "  denominator threshold: to set the threshold of the data space denominator under which the ratio is set to 1" << endl;
  cout << "  minimum image update: to set the minimum of the image update factor under which it stays constant (0 or a negative value" << endl;
  cout << "                        means no minimum thus allowing a 0 update)" << endl;
  cout << "  maximum image update: to set the maximum of the image update factor over which it stays constant (0 or a negative value means" << endl;
  cout << "                        no maximum)" << endl;
}

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

int iOptimizerMLMuMap::ReadConfigurationFile(const string& a_configurationFile)
{
  string key_word = "";
  // Read the initial image value option
  key_word = "initial image value";
  if (ReadDataASCIIFile(a_configurationFile, key_word, &m_initialValue, 1, KEYWORD_MANDATORY))
  {
    Cerr("***** iOptimizerMLMuMap::ReadConfigurationFile() -> Failed to get the '" << key_word << "' keyword !" << endl);
    return 1;
  }
  // Read the denominator threshold option
  key_word = "denominator threshold";
  if (ReadDataASCIIFile(a_configurationFile, key_word, &m_dataSpaceDenominatorThreshold, 1, KEYWORD_MANDATORY))
  {
    Cerr("***** iOptimizerMLMuMap::ReadConfigurationFile() -> Failed to get the '" << key_word << "' keyword !" << endl);
    return 1;
  }
  // Read the minimum image update option
  key_word = "minimum image update";
  if (ReadDataASCIIFile(a_configurationFile, key_word, &m_minimumImageUpdateFactor, 1, KEYWORD_MANDATORY))
  {
    Cerr("***** iOptimizerMLMuMap::ReadConfigurationFile() -> Failed to get the '" << key_word << "' keyword !" << endl);
    return 1;
  }
  // Read the maximum image update option
  key_word = "maximum image update";
  if (ReadDataASCIIFile(a_configurationFile, key_word, &m_maximumImageUpdateFactor, 1, KEYWORD_MANDATORY))
  {
    Cerr("***** iOptimizerMLMuMap::ReadConfigurationFile() -> Failed to get the '" << key_word << "' keyword !" << endl);
    return 1;
  }
  // Normal end
  return 0;
}

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

int iOptimizerMLMuMap::ReadOptionsList(const string& a_optionsList)
{
  // There are 4 floating point variables as options
  const int nb_options = 4;
  FLTNB options[nb_options];
  // Read them
  if (ReadStringOption(a_optionsList, options, nb_options, ",", "MLEM configuration"))
  {
    Cerr("***** iOptimizerMLMuMap::ReadOptionsList() -> Failed to correctly read the list of options !" << endl);
    return 1;
  }
  // Affect options
  m_initialValue = options[0];
  m_dataSpaceDenominatorThreshold = options[1];
  m_minimumImageUpdateFactor = options[2];
  m_maximumImageUpdateFactor = options[3];
  // Normal end
  return 0;
}

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

int iOptimizerMLMuMap::CheckSpecificParameters()
{
  // Check that initial image value is strictly positive
  if (m_initialValue<=0.)
  {
    Cerr("***** iOptimizerMLMuMap->CheckSpecificParameters() -> Provided initial image value (" << m_initialValue << ") must be strictly positive !" << endl);
    return 1;
  }
  // Check that denominator threshold value is strictly positive
  if (m_dataSpaceDenominatorThreshold<=0.)
  {
    Cerr("***** iOptimizerMLMuMap->CheckSpecificParameters() -> Provided data space denominator threshold (" << m_dataSpaceDenominatorThreshold << ") must be strictly positive !" << endl);
    return 1;
  }
  // Check that maximum image update factor is higher than the minimum
  if (m_minimumImageUpdateFactor>0. && m_maximumImageUpdateFactor>0. && m_maximumImageUpdateFactor<m_minimumImageUpdateFactor)
  {
    Cerr("***** iOptimizerMLMuMap->CheckSpecificParameters() -> Provided minimum/maximum (" << m_minimumImageUpdateFactor << "/" << m_maximumImageUpdateFactor << " are inconsistent !" << endl);
    return 1;
  }
  // Can only deal with PET histogram data where ACF values are provideed
  if (m_dataType!=TYPE_PET || m_dataMode!=MODE_HISTOGRAM)
  {
    Cerr("***** iOptimizerMLMuMap->CheckSpecificParameters() -> Can only reconstruct ACF from PET histogram data !" << endl);
    return 1;
  }
  // Get the datafile and check that ACF are provided
  if (!(dynamic_cast<iDataFilePET*>(mp_DataFile))->GetAtnCorrectionFlag())
  {
    Cerr("***** iOptimizerMLMuMap->CheckSpecificParameters() -> ACF must be provided into the PET histogram data !" << endl);
    return 1;
  }
  // Normal end
  return 0;
}

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

int iOptimizerMLMuMap::InitializeSpecific()
{
  // In order to get quantitative attenuation values in the image space, we have to disable all PET related corrections.
  // We cannot disable attenuation correction though. Because the effect of ignoring a correction from the datafile means
  // that when reading the value in the datafile, the value is replaced by a neutral one (1 for ACF), so this would mean
  // that we cannot have access to the ACF anymore within the optimizer. A consequence is that the attenuation will be
  // automatically integrated into the system matrix and we have to remove it manually in the computation of the sensitivity
  // at least. No need to do it in forward and backward projections as the multiplicative factor of the system matrix cancel
  // out in the computation of the ratio (also because there is no scatter or random here).
  mp_ImageDimensionsAndQuantification->SetIgnoreDecaCorrectionFlag(true);
  mp_ImageDimensionsAndQuantification->SetIgnoreBratCorrectionFlag(true);
  mp_ImageDimensionsAndQuantification->SetIgnoreFdurCorrectionFlag(true);
  mp_ImageDimensionsAndQuantification->SetIgnoreCaliCorrectionFlag(true);
  mp_ImageDimensionsAndQuantification->SetIgnoreNormCorrectionFlag(true);
  mp_ImageDimensionsAndQuantification->SetIgnoreRandCorrectionFlag(true);
  mp_ImageDimensionsAndQuantification->SetIgnoreScatCorrectionFlag(true);
  // Have to do it also inside the datafile
  (dynamic_cast<iDataFilePET*>(mp_DataFile))->SetIgnoreNormCorrectionFlag(true);
  (dynamic_cast<iDataFilePET*>(mp_DataFile))->SetIgnoreRandCorrectionFlag(true);
  (dynamic_cast<iDataFilePET*>(mp_DataFile))->SetIgnoreScatCorrectionFlag(true);
  if (mp_ImageDimensionsAndQuantification->ResetQuantificationFactors())
  {
    Cerr("***** oOptimizerMLMuMap->InitializeSpecific() -> A problem occured while reseting the quantification factors !" << endl);
    return 1;
  }
  // Verbose
  if (m_verbose>=VERBOSE_NORMAL)
  {
    Cout("iOptimizerMLMuMap::InitializeSpecific() -> Use the MLEM optimizer" << endl);
    if (m_verbose>=VERBOSE_DETAIL)
    {
      Cout("  --> Initial image value: " << m_initialValue << endl);
      Cout("  --> Data space denominator threshold: " << m_dataSpaceDenominatorThreshold << endl);
      if (m_minimumImageUpdateFactor>0.) Cout("  --> Minimum image update factor: " << m_minimumImageUpdateFactor << endl);
      else Cerr("!!!!! The minimum update value is not set, if using subsets, voxels could be trapped in 0 value causing some negative bias !" << endl);
      if (m_maximumImageUpdateFactor>0.) Cout("  --> Maximum image update factor: " << m_maximumImageUpdateFactor << endl);
    }
  }
  // Normal end
  return 0;
}

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

int iOptimizerMLMuMap::SensitivitySpecificOperations( FLTNB a_data, FLTNB a_forwardModel, FLTNB* ap_weight,
                                                    FLTNB a_multiplicativeCorrections, FLTNB a_additiveCorrections, FLTNB a_blankValue,
                                                    FLTNB a_quantificationFactor, oProjectionLine* ap_Line )
{
  // We multiply by the multiplicative_corrections as they are already in the system matrix and we want to remove them
  *ap_weight = a_multiplicativeCorrections;
  // That's all
  return 0;
}

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

int iOptimizerMLMuMap::DataSpaceSpecificOperations( FLTNB a_data, FLTNB a_forwardModel, FLTNB* ap_backwardValues,
                                                  FLTNB a_multiplicativeCorrections, FLTNB a_additiveCorrections, FLTNB a_blankValue,
                                                  FLTNB a_quantificationFactor, oProjectionLine* ap_Line )
{
  // Do nothing as it is not used
  return 0;
}

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

int iOptimizerMLMuMap::DataStep5ComputeCorrections( oProjectionLine* ap_Line, vEvent* ap_Event,
                                                  int a_bed, int a_timeFrame, int a_respGate, int a_cardGate,
                                                  int a_th )
{
  // Note: This algorithm only works for PET data and only ACF are used, so we ignore TOF and overload this function.
  //       We cast the event to PET in order to specifically get the ACF value.
  
  // Default to NO TOF
  int NOTOF = 0;
  ap_Line->SetCurrentTOFBin(NOTOF);
  // Get the ACF
  FLTNB acf = (dynamic_cast<iEventPET*>(ap_Event))->GetAtnCorrFactor();
  // Get the forward model
  FLTNB forward = m2p_forwardValues[a_th][NOTOF];
  // Truncate acf to 1 if lower
  if (acf<1.) acf = 1.;
  // Take the log
  acf = log(acf);
  // If the foward model is too close to zero, then ignore this data (set to 1 and not 0 because this line is taken into account in the sensitivity)
  if (forward>m_dataSpaceDenominatorThreshold) *m3p_backwardValues[a_th][NOTOF] = acf / forward;
  else *m3p_backwardValues[a_th][NOTOF] = 1.;
  // End
  return 0;
}

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

int iOptimizerMLMuMap::ImageSpaceSpecificOperations( FLTNB a_currentImageValue, FLTNB* ap_newImageValue,
                                                   FLTNB a_sensitivity, FLTNB* ap_correctionValues,
                                                   INTNB a_voxel, int a_tbf, int a_rbf, int a_cbf )
{
  // Compute image update factor
  FLTNB image_update_factor = *ap_correctionValues / a_sensitivity;
  // Apply minimum image update factor
  if ( m_minimumImageUpdateFactor > 0. && image_update_factor < m_minimumImageUpdateFactor ) image_update_factor = m_minimumImageUpdateFactor;
  // Apply maximum image update factor
  if ( m_maximumImageUpdateFactor > 0. && image_update_factor > m_maximumImageUpdateFactor ) image_update_factor = m_maximumImageUpdateFactor;
  // Update image
  *ap_newImageValue = a_currentImageValue * image_update_factor;
  // End
  return 0;
}

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