iPenaltyQuadratic.cc

Go to the documentation of this file.

#include "iPenaltyQuadratic.hh"

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

iPenaltyQuadratic::iPenaltyQuadratic() : vPenalty()
{
  // --------------------------
  // Specific member parameters
  // --------------------------
  m_penaltyID = "Quadratic";
  m_penaltyDerivativesOrder = 1;
  m_targetImagePath = "";
  mp_targetImage = NULL;
}

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

iPenaltyQuadratic::~iPenaltyQuadratic()
{
  // Note: there is no need to deallocate the images themselves as they are allocate using the
  //       miscellaneous image function from the image space, which automatically deals with
  //       memory deallocations.
}

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

void iPenaltyQuadratic::ShowHelpSpecific()
{
  cout << "This penalty is a simple quadratic error penalty, corresponding to a L2 distance between the reconstructed image" << endl;
  cout << "and a target image." << endl;
  cout << "N.B.: The penalty strength is scaled by 0.5." << endl;
  cout << "The following options can be used:" << endl;
  cout << "  target image path: to set a target image to fit. By default, it is a uniform image with zero value." << endl;
}

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

int iPenaltyQuadratic::ReadConfigurationFile(const string& a_configurationFile)
{
  string key_word = "";
  // Read the target image path to use
  key_word = "target image path";
  if (ReadDataASCIIFile(a_configurationFile, key_word, &m_targetImagePath, 1, KEYWORD_MANDATORY))
  {
    Cerr("***** iPenaltyQuadratic::ReadConfigurationFile() -> Failed to get the '" << key_word << "' keyword !" << endl);
    return 1;
  }
  // End
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
 
int iPenaltyQuadratic::ReadOptionsList(const string& a_optionsList)
{ 
  // Too complicated to do it that way
  Cerr("***** iPenaltyQuadratic::ReadOptionsList() -> Options can be specified only using a configuration file !" << endl);
  return 1;
}

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

int iPenaltyQuadratic::CheckSpecificParameters()
{
  // Warn user no target image was set
  if (m_targetImagePath == "")
  {
    Cerr("!!!!! iPenaltyQuadratic::CheckSpecificParameters() -> No target image was provided, thus Quadratic penalty will use a full zeros target image !"<< endl);
  }
  // Normal end
  return 0;
}

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

int iPenaltyQuadratic::InitializeSpecific()
{
  // Allocate the target image
  mp_targetImage = mp_ImageSpace -> AllocateMiscellaneousImage(); // Get a pointer to a newly allocated image
  // Initialize/Read the target image
  if (m_targetImagePath == "") // whole phantom to be considered, so uniform image
  {
    for (int v=0; v<mp_ImageDimensionsAndQuantification->GetNbVoxXYZ(); v++)
    {
        mp_targetImage[v] = 0.;
    }
  }
  else
  {
    // Interfile image reading (INTF_LERP_ENABLED = interpolation allowed)
    if (IntfReadImage(m_targetImagePath, mp_targetImage, mp_ImageDimensionsAndQuantification, m_verbose, INTF_LERP_ENABLED))
    {
      Cerr("***** iPenaltyQuadratic::InitializeSpecific() -> Error reading interfile image '" << m_targetImagePath << "' !" << endl);
      return 1;
    }
  }
  // Verbose
  if (m_verbose>=VERBOSE_NORMAL)
  {
    Cout("iPenaltyQuadratic::InitializeSpecific() -> Use the Quadratic optimizer" << endl);
    if (m_verbose>=VERBOSE_DETAIL)
    {
      if (m_targetImagePath == "") Cout("  --> Uniform zero target image" << endl);
      else Cout("  --> Target image from file: " << m_targetImagePath << endl);
    }
  }
  // Normal end
  return 0;
}

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

FLTNB iPenaltyQuadratic::ComputePenaltyValue(FLTNB* ap_image, INTNB a_voxel, int a_th)
{
  // Compute error between reconstructed image and target image to compute quadratic penalty term
  FLTNB difference;
  difference = ap_image[a_voxel] - mp_targetImage[a_voxel];
  FLTNB penalty = 0.5 * m_penaltyStrength * difference * difference;
  // Check for Inf, Nan, etc
  if (fpclassify(penalty) != FP_NORMAL) penalty = 0.;
  // Return result
  return penalty;
}

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

FLTNB iPenaltyQuadratic::ComputeFirstDerivative(FLTNB* ap_image, INTNB a_voxel, int a_th)
{
  // Compute error between reconstructed image and target image to compute first derivative
  FLTNB difference;
  difference = ap_image[a_voxel] - mp_targetImage[a_voxel];
  FLTNB first_derivative = m_penaltyStrength * difference;
  // Check for Inf, Nan, etc
  if (fpclassify(first_derivative) != FP_NORMAL) first_derivative = 0.;
  // Return result
  return first_derivative;
}


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

FLTNB iPenaltyQuadratic::ComputeSecondDerivative(FLTNB* ap_image, INTNB a_voxel, int a_th)
{
  FLTNB second_derivative = m_penaltyStrength;
  return second_derivative;
}