freAccuracySOMetric.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: freAccuracySOMetric.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 __freAccuracySOMetric_txx
#define __freAccuracySOMetric_txx

#include "freAccuracySOMetric.h"
#include "freStatistics.h"
#include "freGenericSetupToImageAdaptor.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
    AccuracySOMetric<VImageDimension>::
    SetResultFieldPath(const IDPath& path)
  {
    itkDebugMacro("setting ResultFieldPath to " << path.ToStr() );
    this->m_ResultFieldPath = path;
    this->Modified();
  };

  template <unsigned int VImageDimension>
  void
    AccuracySOMetric<VImageDimension>::
    SetReferenceFieldPath(const IDPath& path)
  {
    itkDebugMacro("setting ReferenceFieldPath to " << path.ToStr() );
    this->m_ReferenceFieldPath = path;
    this->Modified();
  };

  template <unsigned int VImageDimension>
  void
    AccuracySOMetric<VImageDimension>::
    SetReferencePointsPath(const IDPath& path)
  {
    itkDebugMacro("setting ReferencePointsPath to " << path.ToStr() );
    this->m_ReferencePointsPath = path;
    this->Modified();
  };

  template <unsigned int VImageDimension>
  void
    AccuracySOMetric<VImageDimension>::
    SetMovingPointsPath(const IDPath& path)
  {
    itkDebugMacro("setting MovingPointsPath to " << path.ToStr() );
    this->m_MovingPointsPath = path;
    this->Modified();
  };

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

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


  template <unsigned int VImageDimension>
  void
    AccuracySOMetric<VImageDimension>::
    InitializeValueComputation() const
  {
    m_dMinError = itk::NumericTraits< double >::max(); //will be changed in ComputeFieldMeasures()/ComputePointsMeasures()
    m_dMaxError = 0.0; //will be changed in ComputeFieldMeasures()/ComputePointsMeasures()
    m_lFailedProcessings = 0;  //will be set by ComputeRegistration();
    m_lCurAdaptation = 0;
    m_lUnevaluatedPoints = 0;  //will be used be EvaluateField()/EvaluatePoints()
    m_dErrorMean = 0.0; //will be changed in ComputeFieldMeasures()/ComputePointsMeasures()
    m_dErrorDev = 0.0; //will be changed in ComputeFieldMeasures()/ComputePointsMeasures()

    m_dDurMean = 0.0; //will be changed in ComputeFieldMeasures()/ComputePointsMeasures()
    m_dDurDev = 0.0; //will be changed in ComputeFieldMeasures()/ComputePointsMeasures()
    m_dMinDur = itk::NumericTraits< double >::max(); //will be changed in ComputeFieldMeasures()/ComputePointsMeasures()
    m_dMaxDur = 0.0; //will be changed in ComputeFieldMeasures()/ComputePointsMeasures()

    m_Errors.clear(); //will be used be EvaluateField()/EvaluatePoints()
    m_Vars.clear(); //will be used be EvaluateField()/EvaluatePoints()
    m_MinErrors.clear(); //will be used be EvaluateField()/EvaluatePoints()
    m_MaxErrors.clear(); //will be used be EvaluateField()/EvaluatePoints()
    m_Samplesize.clear(); //will be used be EvaluateField()
    m_Durations.clear(); //will be set by ComputeRegistration();

    m_PointErrors.clear(); //will be used be EvaluatePoints()
  };

  template <unsigned int VImageDimension>
  void
    AccuracySOMetric<VImageDimension>::
    InitializeMonitor(MonitorType& monitor) const
  {
    //dont have to initialize anything
  };

  template <unsigned int VImageDimension>
  void
    AccuracySOMetric<VImageDimension>::
    InitializeThread(ThreadType& thread) const
  {
    thread.SetUseField(this->m_UseField);
    thread.SetResultFieldPath(this->m_ResultFieldPath);
    thread.SetReferenceFieldPath(this->m_ReferenceFieldPath);
    thread.SetReferencePointsPath(this->m_ReferencePointsPath);
    thread.SetMovingPointsPath(this->m_MovingPointsPath);
  };

  template <unsigned int VImageDimension>
  typename AccuracySOMetric<VImageDimension>::DecomposedMeasureType
    AccuracySOMetric<VImageDimension>::
    ComputeMeasure(MonitorType& monitor) const
  {
    //get all the data from the monitor

    this->m_lFailedProcessings = monitor.GetFailureCount();

    typename MonitorType::AdaptationIDListType evaluations = monitor.GetListOfEvaluatedAdaptations();

    for (typename MonitorType::AdaptationIDListType::const_iterator pos = evaluations.begin(); pos != evaluations.end(); ++pos)
    { //go through all evaluated adaptations and store results;
      EvaluationResultType results;
      monitor.GetAdaptationResults(*pos,results);

      if (this->m_dMinError > results.GetMinError()) this->m_dMinError = results.GetMinError();
      if (this->m_dMaxError < results.GetMaxError()) this->m_dMaxError = results.GetMaxError();

      this->m_lUnevaluatedPoints += results.GetUnevaluatedPoints();

      this->m_Errors.push_back(results.GetError());
      this->m_Vars.push_back(results.GetVariance());
      this->m_Samplesize.push_back(results.GetSamplesize());
      this->m_Durations.push_back(results.GetDuration());

      /* error of all reference points of the current evaluation; only used if m_bUseField = false.*/
      for (typename EvaluationResultType::PointErrorListType::const_iterator pos2 = results.GetPointErrors().begin(); pos2 != results.GetPointErrors().end(); ++pos2)
      {
        this->m_PointErrors.push_back(*pos2);
      }
    }

    if (m_UseField)
    {
      return ComputeFieldMeasures();
    }
    else
    {
      return ComputePointsMeasures();
    }
  };

  template <unsigned int VImageDimension>
  typename AccuracySOMetric<VImageDimension>::DecomposedMeasureType
    AccuracySOMetric<VImageDimension>::
    ComputeFieldMeasures() const
  {
    if (m_ComputeAdaptationMean)
    {
      unsigned long lUsedVariances = 0;
      for (unsigned int iIndex=0; iIndex<m_Errors.size(); iIndex++)
      {
        m_dErrorMean += m_Errors[iIndex];
        if (m_Vars[iIndex]>=0)
        {
          m_dErrorDev += m_Vars[iIndex];
          lUsedVariances++;
        }
      }

      if (m_Errors.size()==0) m_dMinError = 0.0; //there was no test case, at least no successfully evaluated test case
      else m_dErrorMean /= m_Errors.size();
      if (lUsedVariances>0) m_dErrorDev /= lUsedVariances;
      m_dErrorDev = sqrt(m_dErrorDev);
    }
    else
    {
      unsigned long lTotalSamples = 0;
      for (unsigned int iIndex=0; iIndex<m_Errors.size(); iIndex++)
      {
        m_dErrorMean += m_Errors[iIndex] * m_Samplesize[iIndex];
        m_dErrorDev += m_Vars[iIndex] * (m_Samplesize[iIndex]-1);
        lTotalSamples += m_Samplesize[iIndex];
      }

      if (m_Errors.size()==0) m_dMinError = 0.0; //there was no test case, at least no successfully evaluated test case

      if (lTotalSamples>0) m_dErrorMean /= lTotalSamples;
      if (lTotalSamples>1) m_dErrorDev /= lTotalSamples-1;
      m_dErrorDev = sqrt(m_dErrorDev);
    }

    for (unsigned int iIndex=0; iIndex<m_Durations.size(); iIndex++)
    {
      m_dDurMean += m_Durations[iIndex];
      if (m_Durations[iIndex]< m_dMinDur) m_dMinDur = m_Durations[iIndex];
      if (m_Durations[iIndex]> m_dMaxDur) m_dMaxDur = m_Durations[iIndex];
    }

    if (m_Durations.size()==0) m_dMinDur = 0.0; //there was no test case, at least no successfully evaluated test case
    if (m_Durations.size()>0) m_dDurMean /= m_Durations.size();

    if (m_Durations.size()>1)
    {
      for (unsigned int iIndex=0; iIndex<m_Durations.size(); iIndex++)
      {
        m_dDurDev += (m_Durations[iIndex] - m_dDurMean)*(m_Durations[iIndex] - m_dDurMean);
      }
      m_dDurDev = sqrt(m_dDurDev/(m_Durations.size()-1));
    }

    DecomposedMeasureType results(10);
    results.SetElement(0,m_dErrorMean);
    results.SetElement(1,m_dErrorDev);
    results.SetElement(2,m_dMinError);
    results.SetElement(3,m_dMaxError);
    results.SetElement(4,m_lFailedProcessings);
    results.SetElement(5,m_lUnevaluatedPoints);
    results.SetElement(6,m_dDurMean);
    results.SetElement(7,m_dDurDev);
    results.SetElement(8,m_dMinDur);
    results.SetElement(9,m_dMaxDur);

    return results;
  };

  template <unsigned int VImageDimension>
  typename AccuracySOMetric<VImageDimension>::DecomposedMeasureType
    AccuracySOMetric<VImageDimension>::
    ComputePointsMeasures() const
  {
    const unsigned int iMeasurementVectorSize = 1;
    typedef itk::Vector<double,iMeasurementVectorSize> MeasurementVectorType;
    typedef itk::Statistics::ListSample< MeasurementVectorType > SampleType;

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

    MeasurementVectorType mv;

    typedef itk::Statistics::CovarianceCalculator< SampleType > CovarianceAlgorithmType;
    typename CovarianceAlgorithmType::Pointer covarianceAlgorithm = CovarianceAlgorithmType::New();

    if (m_ComputeAdaptationMean)
    {
      unsigned long lUsedVariances = 0;
      for (unsigned int iIndex=0; iIndex<m_Errors.size(); iIndex++)
      {
        m_dErrorMean += m_Errors[iIndex];
        if (m_Vars[iIndex]>=0)
        {
          m_dErrorDev += m_Vars[iIndex];
          lUsedVariances++;
        }
      }

      if (m_Errors.size()==0) m_dMinError = 0.0; //there was no test case, at least no successfully evaluated test case
      else m_dErrorMean /= m_Errors.size();
      if (lUsedVariances>0) m_dErrorDev /= lUsedVariances;
      m_dErrorDev = sqrt(m_dErrorDev);
    }
    else
    {
      if (m_PointErrors.size()>0)
      {
        //mean error, variance, min and max error
        //generate sample
        for (unsigned int iIndex=0; iIndex<m_PointErrors.size(); iIndex++)
        {
          mv[0]=m_PointErrors[iIndex];
          smpSample->PushBack(mv);
        }

        //Calculate mean and variance
        covarianceAlgorithm->SetInputSample( smpSample );
        covarianceAlgorithm->SetMean( 0 );
        covarianceAlgorithm->Update();

        m_dErrorMean = covarianceAlgorithm->GetMean()->GetElement(0);

        if (m_PointErrors.size()>1)
        {
          const CovarianceAlgorithmType::OutputType* matrix = covarianceAlgorithm->GetOutput();
          m_dErrorDev = sqrt((*matrix)(0,0));
        }
        else 
        {
          m_dErrorDev = 0.0;
        }
      }
      else
      {
        m_dMinError = 0.0; //there was no test case, at least no successfully evaluated test case
        //so set min error to 0.0, all other values are still 0.0.
      }
    }

    if (m_Durations.size()>0)
    {
      //duration, duration variance, min and max duration
      //gernerate sample
      smpSample->Clear();
      for (unsigned int iIndex=0; iIndex<m_Durations.size(); iIndex++)
      {
        mv[0]=m_Durations[iIndex];
        smpSample->PushBack(mv);
        //checking for minimum or maximum duration
        if (mv[0]<m_dMinDur) m_dMinDur = mv[0];
        if (mv[0]>m_dMaxDur) m_dMaxDur = mv[0];
      }
      //Calculate duration;
      covarianceAlgorithm->SetInputSample( smpSample );
      covarianceAlgorithm->SetMean( 0 );
      covarianceAlgorithm->Update();

      m_dDurMean = covarianceAlgorithm->GetMean()->GetElement(0);

      if (m_Durations.size()>1)
      {
        const typename CovarianceAlgorithmType::OutputType* matrix = covarianceAlgorithm->GetOutput();
        m_dDurDev = sqrt((*matrix)(0,0));
      }
      else 
      {
        m_dDurDev = 0.0;
      }
    }
    else
    {
      m_dMinDur = 0.0; //there was no test case, at least no successfully evaluated test case
      //so set min error to 0.0, all other values are still 0.0.
    }

    DecomposedMeasureType results(10);
    results.SetElement(0,m_dErrorMean);
    results.SetElement(1,m_dErrorDev);
    results.SetElement(2,m_dMinError);
    results.SetElement(3,m_dMaxError);
    results.SetElement(4,m_lFailedProcessings);
    results.SetElement(5,m_lUnevaluatedPoints);
    results.SetElement(6,m_dDurMean);
    results.SetElement(7,m_dDurDev);
    results.SetElement(8,m_dMinDur);
    results.SetElement(9,m_dMaxDur);

    return results;
  };

} // end namespace FREE

#endif

---