freExhaustiveSOOptimizer.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: freExhaustiveSOOptimizer.h,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 "freExhaustiveSOOptimizer.h"
#include "itkCommand.h"
#include "itkEventObject.h"

namespace FREE
{

ExhaustiveSOOptimizer
::ExhaustiveSOOptimizer()
{

  itkDebugMacro("Constructor");
      
  m_StepLength = 1.0;
  m_CurrentIteration   =   0;
  m_CurrentValue = 0;
  m_CurrentDecomposedValue = 0;
  m_CurrentIndex.Fill(0);
  m_Stop = false;
  m_NumberOfSteps.Fill(0);
  m_ParameterAxis.Fill(-1); //Axis index equals parameter index.
        m_AxisType.Fill(0);
}

unsigned int 
ExhaustiveSOOptimizer
::GetAxisCount() const
{
        return m_NumberOfSteps.GetSize();
};

void ExhaustiveSOOptimizer::StartOptimization( void )
{
        if (m_CostFunction.IsNull()) throwExceptionMacro("Error: cannot optimize. No cost function is defined.");
  this->StartWalking();
}


void
ExhaustiveSOOptimizer
::StartWalking( void )
{

  itkDebugMacro("StartWalking");
        this->InvokeEvent( itk::StartEvent() );

  m_MinimumMetricValuePosition = this->GetInitialPosition();
  m_MaximumMetricValuePosition = this->GetInitialPosition();

        //Default values to ensure that the first step will have better values
        //and so the values will be valid after the first step.
        m_MaximumDecomposedMetricValue.Fill(itk::NumericTraits< DecomposedMeasureType::ValueType >::NonpositiveMin());
        m_MinimumDecomposedMetricValue.Fill(itk::NumericTraits< DecomposedMeasureType::ValueType >::max());
  m_MaximumMetricValue = itk::NumericTraits< MeasureType >::NonpositiveMin();
  m_MinimumMetricValue = itk::NumericTraits< MeasureType >::max();

  m_CurrentIteration          = 0;
  m_MaximumNumberOfIterations = 1;
  
  const unsigned int spaceDimension = this->GetAxisCount();

        if (spaceDimension != m_AxisType.Size()) throwExceptionMacro("Error: not all axis have a specified type. Please ensure the correct initialisation of parameter axis and axis type.");

  for (unsigned int i=0; i< spaceDimension; i++)
    {
    m_MaximumNumberOfIterations *= (2 * m_NumberOfSteps[i] + 1);
    }
    
  m_CurrentIndex.SetSize(spaceDimension);
  m_CurrentIndex.Fill(0);
  
  this->SetCurrentPosition( this->GetInitialPosition() );
  
  itkDebugMacro("Calling ResumeWalking");
  
  this->ResumeWalking();
}


void
ExhaustiveSOOptimizer
::ResumeWalking( void )
{
  itkDebugMacro("ResumeWalk");
  m_Stop = false;
 
  while( !m_Stop ) 
    {
    ParametersType currentPosition = this->GetCurrentPosition();
    
    if( m_Stop )
      {
      StopWalking();
      break;
      }

    m_CurrentValue = this->GetValue( currentPosition );
                m_CurrentDecomposedValue = this->GetCostFunction()->GetCurrentDecomposedValue(); 
    
    if (m_CurrentValue > m_MaximumMetricValue) 
      {
      m_MaximumMetricValue = m_CurrentValue;
      m_MaximumMetricValuePosition = currentPosition;
                        m_MaximumDecomposedMetricValue = m_CurrentDecomposedValue;
      }
    if (m_CurrentValue < m_MinimumMetricValue) 
      {
      m_MinimumMetricValue = m_CurrentValue;
      m_MinimumMetricValuePosition = currentPosition;
                        m_MinimumDecomposedMetricValue = m_CurrentDecomposedValue;
      }
     
    if( m_Stop )
      {
      this->StopWalking();
      break;
      }

    this->AdvanceOneStep();

    m_CurrentIteration++;

    }

        if (m_MinimumIsBest)
        {
                m_BestPosition = m_MinimumMetricValuePosition;
                m_BestValue     = m_MinimumMetricValue;
        }
        else
        {
                m_BestPosition = m_MaximumMetricValuePosition;
                m_BestValue     = m_MaximumMetricValue;
        }
}


void
ExhaustiveSOOptimizer
::StopWalking( void )
{

  itkDebugMacro("StopWalking");
  
  m_Stop = true;
  this->InvokeEvent( itk::EndEvent() );
}


bool
ExhaustiveSOOptimizer
::IsStoppable() const
{
  return true;
};

 
void
ExhaustiveSOOptimizer
::StopOptimization()
{
  this->StopWalking();
};

bool
ExhaustiveSOOptimizer
::IsResumeable() const
{
  return true;
};

void
ExhaustiveSOOptimizer
::ResumeOptimization()
{
  this->ResumeWalking();
};


void
ExhaustiveSOOptimizer
::AdvanceOneStep( void )
{ 

  itkDebugMacro("AdvanceOneStep");

  const unsigned int  spaceDimension =
    m_CostFunction->GetNumberOfParameters();

  ScalesType     scales = this->GetScales();

  // Make sure the scales have been set properly
  if (scales.size() != spaceDimension) throwExceptionMacro("The size of Scales is "<< scales.size()<< ", but the NumberOfParameters is "<< spaceDimension<< ".");

  ParametersType newPosition( spaceDimension );
  ParametersType currentPosition = this->GetCurrentPosition();

  IncrementIndex( newPosition );
  
  itkDebugMacro(<<"new position = " << newPosition );
  
  this->InvokeEvent( itk::IterationEvent() );
  this->SetCurrentPosition( newPosition );
}

void
ExhaustiveSOOptimizer
::IncrementIndex( ParametersType &newPosition ) 
{

  unsigned int idx = 0;

  const unsigned int  spaceDimension =
    m_CostFunction->GetNumberOfParameters();

  const unsigned int  axisDimension = this->GetAxisCount();

        while( idx < axisDimension )
    {
    m_CurrentIndex[idx]++;

    if( m_CurrentIndex[idx] > (2*m_NumberOfSteps[idx]))
      {
      m_CurrentIndex[idx]=0;
      idx++;
      }
    else
      {
      break;
      }
    }
      
  if( idx==axisDimension )
    {
    m_Stop = true;
    }

  for(unsigned int i=0; i<spaceDimension; i++)
    { //calculating new position
                        //Get axis index to retrieve current index and number of steps
                        long axisIndex = m_ParameterAxis[i];
                        if (axisIndex==-1) axisIndex = i;

                        //calculate position
                        if (m_AxisType[axisIndex]==0)
                        { //axis is linear
                                newPosition[i] = (m_CurrentIndex[axisIndex]-m_NumberOfSteps[axisIndex]) *
                                                                                                m_StepLength * this->GetScales()[i] + 
                                                                                                this->GetInitialPosition()[i];
                        }
                        else
                        { //axis is logarithmic
                                double dB = m_StepLength * this->GetScales()[i];
                                double dP = m_CurrentIndex[axisIndex]-m_NumberOfSteps[axisIndex];

                                newPosition[i] = pow(dB, dP) * this->GetInitialPosition()[i];
                        }
    }
}



void
ExhaustiveSOOptimizer
::PrintSelf(std::ostream& os, itk::Indent indent) const
{
  Superclass::PrintSelf(os,indent);

  os << indent << "NumberOfSteps = " << m_NumberOfSteps << std::endl;
  os << indent << "ParameterAxis = " << m_ParameterAxis << std::endl;
  os << indent << "CurrentIteration = " << m_CurrentIteration << std::endl;
  os << indent << "Stop = " << m_Stop << std::endl;
  os << indent << "CurrentParameter = " << m_CurrentParameter << std::endl;
  os << indent << "StepLength = " << m_StepLength << std::endl; 
  os << indent << "CurrentIndex = " << m_CurrentIndex << std::endl;
  os << indent << "MaximumNumberOfIterations = " << m_MaximumNumberOfIterations << std::endl;
  os << indent << "MaximumMetricValue = " << m_MaximumMetricValue << std::endl;
  os << indent << "MinimumMetricValue = " << m_MinimumMetricValue << std::endl;
  os << indent << "MinimumMetricValuePosition = " << m_MinimumMetricValuePosition << std::endl;
  os << indent << "MaximumMetricValuePosition = " << m_MaximumMetricValuePosition << std::endl;
}


} // end namespace FREE

---