freConsistencySOMetric.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: freConsistencySOMetric.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 __freConsistencyRegistrationMetric_txx
#define __freConsistencyRegistrationMetric_txx

#include "freConsistencySOMetric.h"
#include "freGenericSetupToImageAdaptor.h"
#include "freVectorToNormImageFilter.h"
#include "freImageSampleCharacteristicsCalculator.h"
#include "freImageAdaptationInverter.h"

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

namespace FREE
{

  template <unsigned int VImageDimension>
  typename ConsistencyRegistrationMetric<VImageDimension>::DecomposedMeasureType
        ConsistencyRegistrationMetric<VImageDimension>::
        ComputeDecomposedValue( const ParametersType & parameters ) const
        {
                if (this->m_smpAdaptations.IsNull()) throwExceptionMacro("Error. Cannot calculate metric value, adaptiation list is null.");
                if (this->m_Transform.IsNull()) throwExceptionMacro("Error. Cannot calculate metric value, no setup transform set.");

                if (this->m_Transform->GetNumberOfParameters() != parameters.size()) throwExceptionMacro("Error. Cannot calculate metric value, number of passed parameters is not equal to the parameter size of the transform. Passed parameter count: "<< parameters.size());

                m_dMinError = itk::NumericTraits< double >::max();
                m_dMaxError = 0.0;
                m_lFailedRegistrations = 0;
                m_lActAdaptation = 0;
                m_dErrorMean = 0.0;
                m_dErrorVar = 0.0;

                m_Errors.clear();
                m_Vars.clear();
                m_Samplesize.clear();

                m_Transform->SetParameters(parameters);
                m_smpSetup = m_Transform->GenerateTransformedSetup();

                for (unsigned int iAdaptID = 0; iAdaptID<m_smpAdaptations->Size(); iAdaptID++)
                {
                        if (this->fnOnNextAdaptation.IsNotNull()) this->fnOnNextAdaptation->Execute(iAdaptID,"",(void*)this);
                        this->InvokeEvent( NextAdaptationObserverEvent() );

      Adaptation::Adaptation* pAdaptation = m_smpAdaptations->GetElement(iAdaptID);
      Adaptation::Adaptation::Pointer smpReverseAdaptation = ImageAdaptationInverter::InvertAdaptation(*pAdaptation);

      Setup::Pointer smpSetup = Setup::New();
      Setup::Pointer smpReverseSetup = Setup::New();
      smpSetup->operator =(*m_smpSetup);
      smpReverseSetup->operator =(*m_smpSetup);
      
      GenericSetupToImageAdaptor::AdaptSetup(*pAdaptation, *smpSetup);
      GenericSetupToImageAdaptor::AdaptSetup(*(smpReverseAdaptation.GetPointer()), *smpReverseSetup);

                        ComputeRegistrations(smpSetup, smpReverseSetup);

      if (fnOnEvaluationDone.IsNotNull()) fnOnEvaluationDone->Execute(iAdaptID,"",(void*)this);
                  this->InvokeEvent( EvaluationDoneObserverEvent() );
    };

                unsigned long lTotalSamples = 0;
                for (unsigned int iIndex=0; iIndex<m_Errors.size(); iIndex++)
                {
                        m_dErrorMean += m_Errors[iIndex] * m_Samplesize[iIndex];
                        m_dErrorVar += m_Vars[iIndex] * (m_Samplesize[iIndex]-1);
                        lTotalSamples += m_Samplesize[iIndex];
                }

    if (lTotalSamples>0) m_dErrorMean /= lTotalSamples;
    if (lTotalSamples>1) m_dErrorVar /= lTotalSamples-1;

    DecomposedMeasureType results(5);
                results.SetElement(0,m_dErrorMean);
                results.SetElement(1,std::sqrt(m_dErrorVar));
                results.SetElement(2,m_dMinError);
                results.SetElement(3,m_dMaxError);
                results.SetElement(4,m_lFailedRegistrations);

                return results;
        };

template <unsigned int VImageDimension>
ConsistencyRegistrationMetric<VImageDimension>::
ConsistencyRegistrationMetric()
{
  //m_IterationEvent = FREE::IterationEvent<Self>::New(this,&ConsistencyRegistrationMetric<VImageDimension>::OnNextIteration);
};

template <unsigned int VImageDimension>
void
ConsistencyRegistrationMetric<VImageDimension>::
OnNextIteration(const long lIteration, Statistics* pStatistic, void* pSender)
{
        //if (fnOnEvaluatu.IsNotNull()) fnOnNextIteration->Execute(lIteration, 0, pSender);
        //this->InvokeEvent( itk::IterationEvent() );
}

template <unsigned int VImageDimension>
void
ConsistencyRegistrationMetric<VImageDimension>::
ComputeRegistrations(Setup* pSetup, Setup* pReverseSetup) const
{
 // RegistratorType registrator;
 // registrator.fnOnNextIteration= m_IterationEvent;

 // try
        //{
 //   registrator.SetSetup(pSetup);
        //  registrator.InitializeRegistration();

 //   registrator.StartRegistration(); 

 //   RegistratorType::TransformFieldPointer smpField = registrator.GetTransformationField();
 //   registrator.Reset();

 //   registrator.SetSetup(pReverseSetup);
        //  registrator.InitializeRegistration();

 //   registrator.StartRegistration(); 

 //   RegistratorType::TransformFieldPointer smpReverseField = registrator.GetTransformationField();
 //   registrator.Reset();

 //   RegistratorType::TransformFieldType::RegionType region1 = smpField->GetLargestPossibleRegion();
 //   RegistratorType::TransformFieldType::RegionType region2 = smpReverseField->GetLargestPossibleRegion();

 //   region1.Crop(smpReverseField->GetLargestPossibleRegion());
 //   region2.Crop(smpField->GetLargestPossibleRegion());

 //   smpField->SetRequestedRegion(region1);
 //   smpReverseField->SetRequestedRegion(region2);

 //   //TODO::Evaluierung so nicht richtig noch mal in Ruhe berprfen.
 //   typedef itk::AddImageFilter< RegistratorType::TransformFieldType,
 //                                RegistratorType::TransformFieldType,
 //                                RegistratorType::TransformFieldType > AddFilterType;
        //  AddFilterType::Pointer smpAddFilter = AddFilterType::New();
 //     
        //  smpAddFilter->SetInput1(smpField);
        //  smpAddFilter->SetInput2(smpReverseField);
        //  smpAddFilter->Update();
 //   smpField = smpAddFilter->GetOutput();

        //  //computing mean and variance of the field difference
        //  typedef itk::Image<double,VImageDimension> NormImageType;
        //  typedef VectorToNormImageFilter<RegistratorType::TransformFieldType,
 //                                   NormImageType > NormImageFilterType;
        //  typedef ImageSampleCharacteristicsCalculator<NormImageType> CharCalculator;

        //  NormImageFilterType::Pointer smpNormFilter = NormImageFilterType::New();
        //  smpNormFilter->SetInput(smpField);
        //  smpNormFilter->Update();
        //  NormImageType::Pointer smpNormImage = smpNormFilter->GetOutput();

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

        //  m_Errors.push_back(calculator.GetMean());
        //  m_Vars.push_back(calculator.GetVariance());
        //  m_Samplesize.push_back(smpNormImage->GetLargestPossibleRegion().GetNumberOfPixels());

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

        //  NormIterator it(smpNormImage,smpField->GetBufferedRegion());
        //  for ( it.GoToBegin(); !it.IsAtEnd(); ++it )
        //  {
        //        if (it.Get()<m_dMinError) m_dMinError = it.Get();
        //        if (it.Get()>m_dMaxError) m_dMaxError = it.Get();
        //  }
        //}
        //catch( FREE::ExceptionBase & err ) 
        //{ 
        //      if (this->fnOnEvaluationFailed.IsNotNull()) this->fnOnEvaluationFailed->Execute(m_lActAdaptation,(void*)this,std::string(err.what()));
        //      m_lFailedRegistrations++;
        //} 
        //catch(...)
        //{
        //      if (this->fnOnEvaluationFailed.IsNotNull()) this->fnOnEvaluationFailed->Execute(m_lActAdaptation,(void*)this,"unkown error");
        //      m_lFailedRegistrations++;
        //}
};

} // end namespace FREE

#endif

---