freConstrainedMetricInterface.cxx

Go to the documentation of this file.

/*=========================================================================

  Program:   F.R.E.E. - flexible registration evaluation engine
  Version:   v.1.0.0
  Date:      $Date: 2006/09/01 12:00:00 $
  Module:    $RCSfile: freConstrainedMetricInterface.cxx,v $
  Language:  C++



  Copyright (c) 2007 Ralf o Floca (Department of Medical Informatics,
  Institute for Medical Biometry and Informatics, University of Heidelberg,
  Germany). All rights reserved.
  See FREECopyright.txt or http://www.mi.med.uni-hd.de/free/copyright.htm
  for details.

     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/

#include "freConstrainedMetricInterface.h"

#include "muParser.h"

namespace FREE
{

  
  void
        ConstrainedMetricInterface::
  ResetConstraints()
        {
                this->m_Constraints.Reset();
                this->m_BarrierSizes.clear();
        };

  
  void
        ConstrainedMetricInterface::
  AddConstraint(const SetupParameterConstraint& constraint, double dBarrierSize)
        {
                this->m_Constraints.AddElement(constraint);
                this->m_BarrierSizes.push_back(dBarrierSize);
        };


  bool
        ConstrainedMetricInterface::
	ComputeConstraintPenalties( const CnstrParametersType & parameters, CnstrDecomposedMeasureType & penalties ) const
        {
    penalties.SetSize(m_Constraints.Size());
    penalties.Fill(0.0);
    double dPenaltySum = 0.0;

                for (unsigned int iPos = 0; iPos < m_Constraints.Size(); ++iPos)
                {
      double dPenalty = ComputeConstraintPenalty(iPos,parameters);
                        penalties.SetElement(iPos, dPenalty);
      dPenaltySum += dPenalty;
                }

    if (m_FailureThreshold<=0) return true;

    return dPenaltySum<m_FailureThreshold;
        };


  double
        ConstrainedMetricInterface::
	ComputeConstraintPenalty( const unsigned int iConstraintID, const CnstrParametersType & parameters ) const
        {
                SetupParameterConstraint* pConstraint = this->m_Constraints.GetElement(iConstraintID);
                double dConstraintTermValue = 0;
                double dConstrainedValue = 0;

                if (!pConstraint) throwExceptionMacro("Error. Metric constraint specified by ID does not exist. Cannot compute penalty. ID:"<<iConstraintID);
                if ((pConstraint->GetDestinationID()<0)||(pConstraint->GetDestinationID()>=parameters.size()))
                          throwExceptionMacro("Error. Parameter specified by constraint is no part of the passed parameter values. Parameter ID:"<<pConstraint->GetDestinationID());

                dConstrainedValue = parameters[pConstraint->GetDestinationID()];

    //define variable for parser
                double* pParameterData = new double[parameters.size()];
    parameters.copy_out(pParameterData);

    mu::Parser parser;

    try
    {
                  for (unsigned int iPos = 0; iPos < parameters.size(); iPos++)
      {
        parser.DefineVar("_"+Convert::ToStr(iPos),pParameterData+iPos);
      }
      parser.SetExpr(pConstraint->GetConstraintTerm());   
    }
    catch(mu::Parser::exception_type &e)
    {
      std::string sE = e.GetMsg();
                  std::string sLocation;
                  GetClassLocationMacro( sLocation );
      FREE::LogException("muParser based exception", sE, sLocation); \
                  throw FREE::ExceptionBase( "muParser exception", "Error: "+sE, sLocation, __FILE__, __LINE__); \
    }
    catchAllNPassMacro("Unkown errer while setting parser parameter and expression.");

    try
    {
      dConstraintTermValue = parser.Eval();
    }
    catch(mu::Parser::exception_type &e)
    {
      std::string sE = e.GetMsg();
                  std::string sLocation;
                  GetClassLocationMacro( sLocation );
      FREE::LogException("muParser based exception", sE, sLocation); \
                  throw FREE::ExceptionBase( "muParser exception", "Error: "+sE, sLocation, __FILE__, __LINE__); \
    }
    catchAllNPassMacro("Error while evaluating constrained term by muParser.");

    parser.ClearVar();
                delete [] pParameterData;

    //evaluate constraint
    if (pConstraint->GetRelationType() == SetupParameterConstraint::RTEqual)
    {
                        if (dConstraintTermValue == dConstrainedValue) return 0;
                        else return this->m_MaxConstraintPenalty;
    }
    else
    {
                        //transform into an inequality against 0.
                        //basic form of inequalities are dConstrainedValue (relation) dConstraintTermValue
                        double dTransformedTermValue = dConstraintTermValue - dConstrainedValue;
                        if (pConstraint->GetRelationType() == SetupParameterConstraint::RTGreaterOrEqual)
                        { // need the form 0 <= term, but right now it would be 0 >= term, so *-1
                                dTransformedTermValue *= -1;
                        }

                        double dBarrierSize = abs((double)(this->m_BarrierSizes[iConstraintID]));

                        //Check if it would be the maximum penalty anyway.
                        if (dTransformedTermValue<=0) return this->m_MaxConstraintPenalty;

                        //Check if it is outside the barrier region
                        if (dTransformedTermValue>dBarrierSize)
                        { //the constrained value is outside the barrier region within a legal valid range.
                                return 0;
                        }

                        // the constrained value seems to violate the constraint or at least is within the
                        // barrier region, so compute the penalty

                        //calculate the barrier function value
                        double dBarrierValue = -1*log(dTransformedTermValue/dBarrierSize);

                        return std::min(dBarrierValue, this->m_MaxConstraintPenalty);
                }

                return 0;
        };

        

ConstrainedMetricInterface::
ConstrainedMetricInterface()
{
        this->m_MaxConstraintPenalty = 100000;
  this->m_FailureThreshold = 0;
        this->ResetConstraints();
};

} // end namespace FREE

---