frePointSupportedImageToImageMetric.txx

Go to the documentation of this file.

/*=========================================================================
 
  Program:   F.R.E.E. - flexible registration evaluation engine
  Version:   v.0.7.2
  Date:      $Date: 2006/09/01 12:00:00 $
  Module:    $RCSfile: frePointSupportedImageToImageMetric.txx,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.
 
=========================================================================*/
#ifndef _frePointSupportedImageToImageMetric_txx
#define _frePointSupportedImageToImageMetric_txx

#include "frePointSupportedImageToImageMetric.h"
#include "itkImageRegionConstIteratorWithIndex.h"

namespace itk
{

/*
* Constructor
*/
template < class TFixedImage, class TMovingImage, class TFixedPointSet, class TMovingPointSet>
PointSupportedImageToImageMetric<TFixedImage,TMovingImage,TFixedPointSet,TMovingPointSet>
::PointSupportedImageToImageMetric()
{
    itkDebugMacro("Constructor");

    m_Alpha = 1.0;

    m_LastImageMetricValue = 0;
    m_LastPointMetricValue = 0;

    m_LastImageDerivative.Fill(0);
    m_LastPointDerivative.Fill(0);

    MovingPointSetConstPointer m_MovingPointSet;
    FixedPointSetConstPointer m_FixedPointSet;

    ImageMetricPointer m_ImageMetric;

    m_MovingPointSet    = 0; // has to be provided by the user.
    m_FixedPointSet   = 0; // has to be provided by the user.
    m_ImageMetric = 0; // has to be provided by the user.

    this->SetComputeGradient(false); // don't use the default gradient for the point calculation
}

template < class TFixedImage, class TMovingImage, class TFixedPointSet, class TMovingPointSet>
void
PointSupportedImageToImageMetric<TFixedImage,TMovingImage,TFixedPointSet,TMovingPointSet>
::Initialize(void) throw ( ExceptionObject )
{
    if( !m_ImageMetric )
    {
        itkExceptionMacro(<<"ImageMetric is not present");
    }

    m_ImageMetric->SetMovingImage(this->GetMovingImage());
    m_ImageMetric->SetFixedImage(this->GetFixedImage());
    m_ImageMetric->SetTransform(this->m_Transform);
    m_ImageMetric->SetInterpolator(this->m_Interpolator);
    m_ImageMetric->SetFixedImageRegion(this->GetFixedImageRegion());
    m_ImageMetric->SetMovingImageMask(this->m_MovingImageMask);
    m_ImageMetric->SetFixedImageMask(this->m_FixedImageMask);

    Superclass::Initialize();
    m_ImageMetric->Initialize();
}

/*
* Get the match Measure
*/
template < class TFixedImage, class TMovingImage, class TFixedPointSet, class TMovingPointSet>
typename PointSupportedImageToImageMetric<TFixedImage,TMovingImage,TFixedPointSet,TMovingPointSet>::MeasureType
PointSupportedImageToImageMetric<TFixedImage,TMovingImage,TFixedPointSet,TMovingPointSet>
::GetValue( const TransformParametersType & parameters ) const
{
    itkDebugMacro("GetValue( " << parameters << " ) ");

    if( m_ImageMetric.IsNull())
    {
        itkExceptionMacro( << "Image metric has not been assigned" );
    }

    m_LastImageMetricValue = m_ImageMetric->GetValue(parameters);
    m_LastPointMetricValue = this->GetPointValue(parameters);

    m_LastPointMetricValue *= m_Alpha;
    MeasureType measure = m_LastImageMetricValue + m_LastPointMetricValue;

    return measure;
}


/*
* Get the Derivative Measure
*/
template < class TFixedImage, class TMovingImage, class TFixedPointSet, class TMovingPointSet>
void
PointSupportedImageToImageMetric<TFixedImage,TMovingImage,TFixedPointSet,TMovingPointSet>
::GetDerivative( const TransformParametersType & parameters,
                 DerivativeType & derivative  ) const
{
    MeasureType value;
    // call the combined version
    this->GetValueAndDerivative( parameters, value, derivative );
}


/*
* Get both the match Measure and theDerivative Measure 
*/
template < class TFixedImage, class TMovingImage, class TFixedPointSet, class TMovingPointSet>
void
PointSupportedImageToImageMetric<TFixedImage,TMovingImage,TFixedPointSet,TMovingPointSet>
::GetValueAndDerivative(const TransformParametersType & parameters,
                        MeasureType & value, DerivativeType  & derivative) const
{
    itkDebugMacro("GetValueAndDerivative( " << parameters << " ) ");

    if( m_ImageMetric.IsNull())
    {
        itkExceptionMacro( << "Image metric has not been assigned" );
    }

    m_ImageMetric->GetValueAndDerivative(parameters, m_LastImageMetricValue, m_LastImageDerivative);
    this->GetPointValueAndDerivative(parameters, m_LastPointMetricValue, m_LastPointDerivative);
    
    m_LastPointMetricValue *= m_Alpha;
    value = m_LastImageMetricValue + m_LastPointMetricValue;

    derivative = m_LastImageDerivative;

    for (unsigned int i=0; i<derivative.Size(); i++)
    {
      m_LastPointDerivative[i] *= m_Alpha;
      derivative += m_LastPointDerivative[i];
    }
}

/*
* Get the match Measure
*/
template < class TFixedImage, class TMovingImage, class TFixedPointSet, class TMovingPointSet>
typename PointSupportedImageToImageMetric<TFixedImage,TMovingImage,TFixedPointSet,TMovingPointSet>::MeasureType
PointSupportedImageToImageMetric<TFixedImage,TMovingImage,TFixedPointSet,TMovingPointSet>
::GetPointValue( const TransformParametersType & parameters ) const
{
  itkDebugMacro("GetPointValue( " << parameters << " ) ");

  if( m_MovingPointSet.IsNull() || m_FixedPointSet.IsNull())
  {
    return 0;
  }

  if( m_MovingPointSet->GetNumberOfPoints()!=m_FixedPointSet->GetNumberOfPoints())
  {
    itkExceptionMacro( << "Point sets have not the same size. Moving point set: "<<m_MovingPointSet->GetNumberOfPoints()<<"; fixed point set: "<< m_FixedPointSet->GetNumberOfPoints() );
  }

  typedef typename FixedPointSetType::PointsContainer::ConstIterator  FixedPointIterator;
  typedef typename MovingPointSetType::PointsContainer::ConstIterator  MovingPointIterator;

  MovingPointIterator movingPointItr = m_MovingPointSet->GetPoints()->Begin();
  MovingPointIterator movingPointEnd = m_MovingPointSet->GetPoints()->End();

  FixedPointIterator fixedPointItr = m_FixedPointSet->GetPoints()->Begin();

  MeasureType measure = NumericTraits< MeasureType >::Zero;

  this->SetTransformParameters( parameters );

  while( movingPointItr != movingPointEnd )
  {
    InputPointType  fixedPoint;
    fixedPoint.CastFrom( fixedPointItr.Value() );
    InputPointType  movingPoint;
    movingPoint.CastFrom( movingPointItr.Value() );
   
    //Because we use the same transform as the image metric the transform is
    //a back transform, so it maps from fixed to moving
    typename Superclass::OutputPointType transformedPoint = 
    this->m_Transform->TransformPoint( fixedPoint);

    typename InputPointType::VectorType distance = movingPoint - transformedPoint;

    measure += distance.GetSquaredNorm();

    ++movingPointItr;
    ++fixedPointItr;
  }

  if (m_MovingPointSet->GetNumberOfPoints() >0)
  {
    measure /= m_MovingPointSet->GetNumberOfPoints ();
  }

  return measure;
}


/*
* Get the Derivative Measure
*/
template < class TFixedImage, class TMovingImage, class TFixedPointSet, class TMovingPointSet>
void
PointSupportedImageToImageMetric<TFixedImage,TMovingImage,TFixedPointSet,TMovingPointSet>
::GetPointDerivative( const TransformParametersType & parameters,
                 DerivativeType & derivative  ) const
{
    MeasureType value;
    // call the combined version
    this->GetPointValueAndDerivative( parameters, value, derivative );
}


/*
* Get both the match Measure and theDerivative Measure 
*/
template < class TFixedImage, class TMovingImage, class TFixedPointSet, class TMovingPointSet>
void
PointSupportedImageToImageMetric<TFixedImage,TMovingImage,TFixedPointSet,TMovingPointSet>
::GetPointValueAndDerivative(const TransformParametersType & parameters,
                        MeasureType & value, DerivativeType  & derivative) const
{
  itkDebugMacro("GetPointValueAndDerivative( " << parameters << " ) ");

  const unsigned int ParametersDimension = this->GetNumberOfParameters();
  derivative = DerivativeType( ParametersDimension );
  derivative.Fill( NumericTraits<ITK_TYPENAME DerivativeType::ValueType>::Zero );

  value = NumericTraits< MeasureType >::Zero;

  if( m_MovingPointSet.IsNotNull() && m_FixedPointSet.IsNotNull())
  {
    if( m_MovingPointSet->GetNumberOfPoints()!=m_FixedPointSet->GetNumberOfPoints())
    {
      itkExceptionMacro( << "Point sets have not the same size. Moving point set: "<<m_MovingPointSet->GetNumberOfPoints()<<"; fixed point set: "<< m_FixedPointSet->GetNumberOfPoints() );
    }

    typedef typename FixedPointSetType::PointsContainer::ConstIterator  FixedPointIterator;
    typedef typename MovingPointSetType::PointsContainer::ConstIterator  MovingPointIterator;

    MovingPointIterator movingPointItr = m_MovingPointSet->GetPoints()->Begin();
    MovingPointIterator movingPointEnd = m_MovingPointSet->GetPoints()->End();

    FixedPointIterator fixedPointItr = m_FixedPointSet->GetPoints()->Begin();

    this->SetTransformParameters( parameters );

    while( movingPointItr != movingPointEnd )
    {
      InputPointType  fixedPoint;
      fixedPoint.CastFrom( fixedPointItr.Value() );
      InputPointType  movingPoint;
      movingPoint.CastFrom( movingPointItr.Value() );
     
      //Because we use the same transform as the image metric the transform is
      //a back transform, so it maps from fixed to moving
      typename Superclass::OutputPointType transformedPoint = 
      this->m_Transform->TransformPoint( fixedPoint);

      typename InputPointType::VectorType distance = movingPoint - transformedPoint;

      const TransformJacobianType & jacobian =
                  this->m_Transform->GetJacobian( fixedPoint );

      for(unsigned int par=0; par<ParametersDimension; par++)
      {
          RealType sum = NumericTraits< RealType >::Zero;

          for(unsigned int dim=0; dim<FixedPointSetType::PointDimension; dim++)
          {
              sum += -2.0 * distance[dim] * jacobian( dim, par );
          }

          derivative[par] += sum;
      }

      value += distance.GetSquaredNorm();

      ++movingPointItr;
      ++fixedPointItr;
    }

    if (m_MovingPointSet->GetNumberOfPoints()>0)
    {
      value /= m_MovingPointSet->GetNumberOfPoints();

      for(unsigned int par=0; par<ParametersDimension; par++)
      {
        derivative[par] /= m_MovingPointSet->GetNumberOfPoints();
      }
    }
  }
}

} // end namespace itk


#endif

---