freAccuracySOMetricThread.txx

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: freAccuracySOMetricThread.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 __freAccuracySOMetricThread_txx
#define __freAccuracySOMetricThread_txx

#include "freAccuracySOMetricThread.h"

#include "freVectorToNormImageFilter.h"
#include "freImageSampleCharacteristicsCalculator.h"
#include "freSessionProcessor.h"

#include "itkNumericTraits.h"
#include "itkSubtractImageFilter.h"
#include "itkImageAdaptor.h"
#include "itkPointSet.h"
#include "itkListSample.h"

namespace FREE
{

template <unsigned int VImageDimension>
void
AccuracySOMetricThread<VImageDimension>:: 
SetResultFieldPath(const IDPath& path)
{
  this->LockExecutionMutex();
    this->m_ResultFieldPath = path;
  this->UnlockExecutionMutex();
};

template <unsigned int VImageDimension>
const IDPath&
AccuracySOMetricThread<VImageDimension>:: 
GetResultFieldPath() const {return m_ResultFieldPath;};

template <unsigned int VImageDimension>
void
AccuracySOMetricThread<VImageDimension>:: 
SetReferenceFieldPath(const IDPath& path)
{
  this->LockExecutionMutex();
    this->m_ReferenceFieldPath = path;
  this->UnlockExecutionMutex();
};

template <unsigned int VImageDimension>
const IDPath&
AccuracySOMetricThread<VImageDimension>:: 
GetReferenceFieldPath() const {return m_ReferenceFieldPath;};

template <unsigned int VImageDimension>
void
AccuracySOMetricThread<VImageDimension>:: 
SetReferencePointsPath(const IDPath& path)
{
  this->LockExecutionMutex();
    this->m_ReferencePointsPath = path;
  this->UnlockExecutionMutex();
};

template <unsigned int VImageDimension>
const IDPath&
AccuracySOMetricThread<VImageDimension>:: 
GetReferencePointsPath() const {return m_ReferencePointsPath;};

template <unsigned int VImageDimension>
void
AccuracySOMetricThread<VImageDimension>:: 
SetMovingPointsPath(const IDPath& path)
{
  this->LockExecutionMutex();
    this->m_MovingPointsPath = path;
  this->UnlockExecutionMutex();
};

template <unsigned int VImageDimension>
const IDPath&
AccuracySOMetricThread<VImageDimension>:: 
GetMovingPointsPath() const {return m_MovingPointsPath;};

template <unsigned int VImageDimension>
void
AccuracySOMetricThread<VImageDimension>:: 
SetUseField(const bool& bUseField)
{
  this->LockExecutionMutex();
    this->m_UseField = bUseField;
  this->UnlockExecutionMutex();
};

template <unsigned int VImageDimension>
bool
AccuracySOMetricThread<VImageDimension>:: 
GetUseField() const {return m_UseField;};

template <unsigned int VImageDimension>
AccuracySOMetricThread<VImageDimension>:: 
AccuracySOMetricThread()
{
  m_UseField = true;
  m_ReferenceFieldPath.Reset();
  m_ReferencePointsPath.Reset();
  m_MovingPointsPath.Reset();
  m_ResultFieldPath.Reset();
};

template <unsigned int VImageDimension>
AccuracySOMetricThread<VImageDimension>:: 
~AccuracySOMetricThread() {};

template <unsigned int VImageDimension>
bool
AccuracySOMetricThread<VImageDimension>:: 
ProcessSetup(Setup* pAdaptationSetup) throw()
{
  bool result = false;

  ControllerCentral::AddOnProgressEvent(this->m_EvaluatedSessionProgressEvent);

  try
        {
    SessionProcessor processor;

    processor.SetSetup(pAdaptationSetup);
    processor.DefineOutput("resultField", this->m_ResultFieldPath);

    processor.InitializeSession();
        
    this->m_EvaluatedSessionProgressEvent->SetACR(processor.GetSessionInfo()->GetSessionCache());

    clock_t start = clock();

    typename TransformationFieldType::Pointer smpResult = processor.GetCastedOutput<TransformationFieldType>("resultField");
                
    clock_t stop = clock();
    double dDuration = (stop-start) / (CLOCKS_PER_SEC/10); //convert clocks to thenth second

    //evaluating registration result
                if (m_UseField)
                {
      EvaluateField(smpResult, processor.GetSessionInfo());
                }
                else
                {
                        EvaluatePoints(smpResult, processor.GetSessionInfo());
                }

                m_Results.SetDuration((unsigned long)dDuration);

    result = true;
        }
        catch( FREE::ExceptionBase & err ) 
        { 
                this->m_FailureComment = std::string(err.what());
        } 
        catch(...)
        {
                this->m_FailureComment = "unkown error";
        }

  ControllerCentral::RemoveOnProgressEvent(m_EvaluatedSessionProgressEvent);

  return result;
};

template <unsigned int VImageDimension>
void
AccuracySOMetricThread<VImageDimension>:: 
EvaluateField(TransformationFieldType* pResultField, SessionInfo* pSessionInfo)
{
  typename TransformationFieldType::Pointer smpRefField = NULL;

  if (m_ReferenceFieldPath.IsEmpty()) throwExceptionMacro("Error: Reference field path is null, but should be used to measure registration. Ensure correct initialisation of the metric.");

  try
  {
    SessionAccessor::GenericMediaPointer smpFieldReference = SessionAccessor::GetMedia(m_ReferenceFieldPath,pSessionInfo);
    smpRefField = dynamic_cast<TransformationFieldType*>(smpFieldReference.GetPointer());
  }
  catchAllNPassMacro("Error: cannot retrieve reference field. Ensure correct IDPath. Used path: "<<m_ReferenceFieldPath.ToStr());
  if (smpRefField.IsNull()) throwExceptionMacro("Error: Reference field is of wrong type. Ensure correct initialisation of the metric. Field class name: "<<smpRefField->GetNameOfClass());

  //calculating difference between reference field and final field
        typedef itk::SubtractImageFilter< TransformationFieldType, TransformationFieldType, TransformationFieldType > SubtractFilterType;
        typename SubtractFilterType::Pointer smpSubFilter = SubtractFilterType::New();
                
        //computing mean and variance of the field difference
        typedef itk::Image<double,VImageDimension> NormImageType;
        typedef VectorToNormImageFilter<TransformationFieldType, NormImageType > NormImageFilterType;
        typedef ImageSampleCharacteristicsCalculator<NormImageType> CharCalculator;
        typename NormImageFilterType::Pointer smpNormFilter = NormImageFilterType::New();

  smpSubFilter->SetInput1(pResultField);
        smpSubFilter->SetInput2(smpRefField);
        smpNormFilter->SetInput(smpSubFilter->GetOutput());
        smpNormFilter->Update();
        typename NormImageType::Pointer smpNormImage = smpNormFilter->GetOutput();

        CharCalculator calculator;
        calculator.SetImage(smpNormImage);
        calculator.Compute();

        this->m_Results.SetError(calculator.GetMean());
        this->m_Results.SetVariance(calculator.GetVariance());
        this->m_Results.SetSamplesize(pResultField->GetLargestPossibleRegion().GetNumberOfPixels());

        //checking for minimum or maximum
        typedef itk::ImageRegionConstIterator<NormImageType> NormIterator;

  double dMin = itk::NumericTraits< double >::max();
  double dMax = 0.0;

        NormIterator it(smpNormImage,pResultField->GetBufferedRegion());
        for ( it.GoToBegin(); !it.IsAtEnd(); ++it )
        {
                if (it.Get()<dMin) dMin = it.Get();
                if (it.Get()>dMax) dMax = it.Get();
        }

  this->m_Results.SetMinError(dMin);
  this->m_Results.SetMaxError(dMax);
};

template <unsigned int VImageDimension>
void
AccuracySOMetricThread<VImageDimension>:: 
EvaluatePoints(TransformationFieldType* pResultField, SessionInfo* pSessionInfo)
{
  typedef TransformationFieldType FieldType;
  typedef typename TransformationFieldType::Pointer FieldPointer;

        typedef typename ImageTypes<VImageDimension>::PointSetType PointSetType;
        typedef typename PointSetType::Pointer PointSetPointer;

  PointSetPointer smpRefPoints = 0;
  PointSetPointer smpMovingPoints = 0;

  if (m_ReferencePointsPath.IsEmpty()) throwExceptionMacro("Error: Reference points path is null, but should be used to measure registration. Ensure correct initialisation of the metric.");
  try
  {
    typename SessionAccessor::GenericMediaPointer genPoints = SessionAccessor::GetMedia(m_ReferencePointsPath,pSessionInfo);
    smpRefPoints = dynamic_cast<PointSetType*>(genPoints.GetPointer());
  }
  catchAllNPassMacro("Error: cannot retrieve reference point set Ensure correct IDPath. Used path: "<<m_ReferenceFieldPath.ToStr());

  if (m_MovingPointsPath.IsEmpty()) throwExceptionMacro("Error: Reference points path is null, but should be used to measure registration. Ensure correct initialisation of the metric.");
  try
  {
    SessionAccessor::GenericMediaPointer genPoints = SessionAccessor::GetMedia(m_MovingPointsPath,pSessionInfo);
    smpMovingPoints = dynamic_cast<PointSetType*>(genPoints.GetPointer());
  }
  catchAllNPassMacro("Error: cannot retrieve reference point set Ensure correct IDPath. Used path: "<<m_ReferenceFieldPath.ToStr());

  if (smpRefPoints->GetNumberOfPoints() != smpMovingPoints->GetNumberOfPoints()) throwExceptionMacro("Error: reference point set and moving point set must have same number of points. Reference count: " << smpRefPoints->GetNumberOfPoints() << "; moving count: " << smpMovingPoints->GetNumberOfPoints());
  
  //map each reference point into the moving domaime and calculate difference

        typename PointSetType::PointsContainer::Iterator refIter = smpRefPoints->GetPoints()->Begin();
  typename PointSetType::PointsContainer::Iterator movingIter = smpMovingPoints->GetPoints()->Begin();

  typename PointSetType::PointsContainer::Iterator endOfPoints = smpRefPoints->GetPoints()->End();

        const unsigned int iMeasurementVectorSize = 1;
  typedef itk::Vector<double,iMeasurementVectorSize> MeasurementVectorType;
  typedef itk::Statistics::ListSample< MeasurementVectorType > SampleType;

        typename SampleType::Pointer smpSample = SampleType::New();
        smpSample->SetMeasurementVectorSize(iMeasurementVectorSize);

  double dMin = itk::NumericTraits< double >::max();
  double dMax = 0.0;

  while (refIter != endOfPoints)
  {
    typename PointSetType::PointType refPoint = refIter.Value();
    typename PointSetType::PointType movingPoint = movingIter.Value();

    typename FieldType::PixelType translation;
    translation.Fill(0.0);

    typename FieldType::IndexType index;
                bool bInField =  pResultField->TransformPhysicalPointToIndex(refPoint, index);

                if (bInField)
    { //referencepoint is within the transformation field
      translation = pResultField->GetPixel(index);
            refPoint += translation;

      MeasurementVectorType mv;

                  typename PointSetType::PointType::VectorType difVector = movingPoint - refPoint;
      mv[0] = difVector.GetNorm();

                        this->m_Results.GetPointErrors().push_back(mv[0]);
      smpSample->PushBack(mv);

                  if (mv[0]<dMin) dMin = mv[0];
                  if (mv[0]>dMax) dMax = mv[0];
                }
                else
                {
                        this->m_Results.GetUnevaluatedPoints()++;
                }
                ++refIter;
                ++movingIter;
  }

        // Calculate mean and variance
        typedef itk::Statistics::CovarianceCalculator< SampleType > CovarianceAlgorithmType;

        typename CovarianceAlgorithmType::Pointer covarianceAlgorithm = CovarianceAlgorithmType::New();
        covarianceAlgorithm->SetInputSample( smpSample );
        covarianceAlgorithm->SetMean( 0 );
        covarianceAlgorithm->Update();

        double mean = covarianceAlgorithm->GetMean()->GetElement(0);
        double var = (*(covarianceAlgorithm->GetOutput()))(0,0);

        this->m_Results.SetError(mean);
        this->m_Results.SetVariance(var);
  this->m_Results.SetMinError(dMin);
  this->m_Results.SetMaxError(dMax);
};



} // end namespace FREE

#endif

---