/*
 * Copyright (c) ICG. All rights reserved.
 * See copyright.txt for more information.
 *
 * Institute for Computer Graphics and Vision
 * Graz, University of Technology / Austria
 *
 *
 * 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.
 *
 *
 * Project     : MIPItkProjects
 * Module      : Evaluation
 * Class       : $RCSfile: CompareDisplacementFields.cxx,v $
 * Language    : C++
 * Description : 
 *
 * Author     : Martin Urschler
 * EMail      : urschler@icg.tu-graz.ac.at
 * Date       : $Date: 2007-03-07 10:01:10 $
 * Version    : $Revision: 1.35 $
 * Full Id    : $Id: CompareDisplacementFields.cxx,v 1.35 2007-03-07 10:01:10 urschler Exp $
 *
 */


#include "commonItkTypedefs.h"
#include "VolumeIOWrapper.h"
#include "SmallUtilityMethods.h"
#include "NonlinearRegistrationUtilities.h"

#include "itkJoinImageFilter.h"
#include "itkImageToHistogramGenerator.h"
#include "itkImageRegionConstIterator.h"
#include "itkConstNeighborhoodIterator.h"
#include "itkCannyEdgeDetectionImageFilter.h"
#include "itkRecursiveGaussianImageFilter.h"

void calculateMutualInformationStatistics(
   const Int16ImageType::Pointer image1,
   const Int16ImageType::Pointer image2,
   double& nmi )
{
   // Using the \doxygen{JoinImageFilter} we use the two input images and put 
   // them together in an image of two components.

   typedef itk::JoinImageFilter< Int16ImageType, Int16ImageType > 
      JoinFilterType;

   JoinFilterType::Pointer joinFilter = JoinFilterType::New();
  
   joinFilter->SetInput1( image1 );
   joinFilter->SetInput2( image2 );
   joinFilter->Update();
   
   // We prepare now the types to be used for the computation of the Joint
   // histogram. For this purpose, we take the type of the image resulting from
   // the JoinImageFilter and use it as template argument of the
   // \doxygen{ImageToHistogramGenerator}. We then construct one by invoking the
   // \code{New()} method. 
   
   typedef JoinFilterType::OutputImageType VectorImageType;

   typedef itk::Statistics::ImageToHistogramGenerator< 
      VectorImageType >  HistogramGeneratorType;

   HistogramGeneratorType::Pointer histogramGenerator = 
      HistogramGeneratorType::New();
   
   // We pass the multiple components image as input to the histogram generator,
   // and setup the marginal scale value that will define the precision to be 
   // used for classifying values into the histogram bins.
   histogramGenerator->SetInput(  joinFilter->GetOutput()  );
   histogramGenerator->SetMarginalScale( 10.0 );
   
   // We must now define the number of bins to use for each one of the 
   // components in the joint image. For this purpose we take the SizeType from 
   // the traits of the histogram generator type. The array of number of bins is
   // passed to the generator and we can then invoke the \code{Compute()} method
   // in order to trigger the computation of the joint histogram.
   typedef HistogramGeneratorType::SizeType   SizeType;

   SizeType size;

   size[0] = 255;  // number of bins for the first  channel
   size[1] = 255;  // number of bins for the second channel

   histogramGenerator->SetNumberOfBins( size );
   histogramGenerator->Compute();
   
   // The histogram can be recovered from the generator by creating a variable
   // with the histogram type taken from the generator traits.
   typedef HistogramGeneratorType::HistogramType  HistogramType;
   const HistogramType * histogram = histogramGenerator->GetOutput();
   
   // We now walk over all the bins of the joint histogram and compute their
   // contribution to the value of the joint Entropy. For this purpose we use
   // histogram iterators, and the Begin() and End() methods. Since
   // the values returned from the histogram are measuring frequency we must
   // convert them to an estimation of probability by dividing them over the 
   // total sum of frequencies returned by the \code{GetTotalFrequency()} method
   HistogramType::ConstIterator itr = histogram->Begin();
   HistogramType::ConstIterator end = histogram->End();
 
   const double Sum = histogram->GetTotalFrequency();
   
   // We initialize to zero the variable to use for accumulating the value of 
   // the joint entropy, and then use the iterator for visiting all the bins of 
   // the joint histogram. For every bin we compute their contribution to the 
   // reduction of uncertainty. Note that in order to avoid logarithmic 
   // operations on zero values, we skip over those bins that have less than one
   // count. The entropy contribution must be computed using logarithms in base
   // two in order to be able express entropy in \textbf{bits}.
   double JointEntropy = 0.0;

   while( itr != end )
   {
      const double count = itr.GetFrequency();
      if( count > 0.0 )
      {
         const double probability = count / Sum;
         JointEntropy += - probability * log( probability ) / log( 2.0 );
      }
      ++itr;
   }
   
   std::cout << "Joint Entropy      = " << JointEntropy << " bits " 
      << std::endl;
   
   // Now that we have the value of the joint entropy we can proceed to estimate
   // the values of the entropies for each image independently. This can be done
   // by simply changing the number of bins and then recomputing the histogram.
   
   size[0] = 255;  // number of bins for the first  channel
   size[1] =   1;  // number of bins for the second channel

   histogramGenerator->SetNumberOfBins( size );
   histogramGenerator->Compute();
   
   // We initialize to zero another variable in order to start accumulating the
   // entropy contributions from every bin.
   
   itr = histogram->Begin();
   end = histogram->End();

   double Entropy1 = 0.0;

   while( itr != end )
   {
      const double count = itr.GetFrequency();
      if( count > 0.0 )
      {
         const double probability = count / Sum;
         Entropy1 += - probability * log( probability ) / log( 2.0 );
      }
       ++itr;
   }

   std::cout << "Image1 Entropy   = " << Entropy1 << " bits " << std::endl;
   
   // The same process is used for computing the entropy of the other component.
   // Simply by swapping the number of bins in the histogram.
   
   size[0] =   1;  // number of bins for the first channel
   size[1] = 255;  // number of bins for the second channel

   histogramGenerator->SetNumberOfBins( size );
   histogramGenerator->Compute();
   
   // The entropy is computed in a similar manner, just by visiting all the bins
   // on the histogram and accumulating their entropy contributions.
   
   itr = histogram->Begin();
   end = histogram->End();

   double Entropy2 = 0.0;

   while( itr != end )
   {
      const double count = itr.GetFrequency();
      if( count > 0.0 )
      {
         const double probability = count / Sum;
         Entropy2 += - probability * log( probability ) / log( 2.0 );
      }
      ++itr;
   }

   std::cout << "Image2 Entropy   = " << Entropy2 << " bits " << std::endl;
   
   double MutualInformation = Entropy1 + Entropy2 - JointEntropy;
  
   std::cout << "Mutual Information = " << MutualInformation << " bits " << 
      std::endl;
      
   // or Normalized Mutual Information, where the value of Mutual Information 
   // gets divided by the mean entropy of the input images.
   
   double NormalizedMutualInformation1 = 
                     2.0 * MutualInformation / ( Entropy1 + Entropy2 );

   std::cout << "Normalized Mutual Information 1 = " << 
      NormalizedMutualInformation1 <<  std::endl;
      
      
   nmi = NormalizedMutualInformation1;
}

double calculateEdgeOverlapMeasure(
   const Int16ImageType::Pointer image1,
   const Int16ImageType::Pointer image2,
   const int warpingDefaultValue )
{
   double edgeOverlapMeasure = 0.0;
   
   // perform canny edge detection on both images
   
   const Int16ImageType::SpacingType spacing1 = image1->GetSpacing();
   const Int16ImageType::SpacingType spacing2 = image2->GetSpacing();
   
   typedef itk::CannyEdgeDetectionImageFilter<FloatImageType,FloatImageType>
      CannyFilterType;
   
   FloatImageType::Pointer edges1 = 0;
   {
      FloatImageType::Pointer image1_float = 
         SmallUtilityMethodsCastItkVolume<Int16ImageType,FloatImageType>(image1);
      
      CannyFilterType::Pointer canny = CannyFilterType::New();
   
      canny->SetInput( image1_float );
   
      CannyFilterType::ArrayType variance;
      variance[0] = vnl_math_sqr(3.0*spacing1[0]);
      variance[1] = vnl_math_sqr(2.0*spacing1[1]);
      variance[2] = vnl_math_sqr(2.0*spacing1[2]);
      
      canny->SetVariance(variance);
       
      ITK_EXCEPTION_CHECKED( "perform canny on image 1", canny->Update(), 0.0);
      
      edges1 = canny->GetOutput(); 
   }
   
   std::cout << 
      SmallUtilityMethods<FloatImageType>::getMinimumVoxelValue(edges1) <<
      std::endl;
   
   std::cout << 
      SmallUtilityMethods<FloatImageType>::getMaximumVoxelValue(edges1) <<
      std::endl;
      
   //VolumeIOWrapper<FloatImageType>::writeITKVolume16Bit( 
   //   edges1, "D:\\temp\\edges1.hdr", "" );
   
   typedef itk::CannyEdgeDetectionImageFilter<FloatImageType,FloatImageType>
      CannyFilterType;
   
   FloatImageType::Pointer edges2 = 0;
   {
      FloatImageType::Pointer image2_float = 
         SmallUtilityMethodsCastItkVolume<Int16ImageType,FloatImageType>(image2);
      
      CannyFilterType::Pointer canny = CannyFilterType::New();
   
      canny->SetInput( image2_float );
   
      CannyFilterType::ArrayType variance;
      variance[0] = vnl_math_sqr(3.0*spacing2[0]);
      variance[1] = vnl_math_sqr(2.0*spacing2[1]);
      variance[2] = vnl_math_sqr(2.0*spacing2[2]);
      
      canny->SetVariance(variance);
       
      ITK_EXCEPTION_CHECKED( "perform canny on image 2", canny->Update(), 0.0);
      
      edges2 = canny->GetOutput(); 
   }
   
   std::cout << 
      SmallUtilityMethods<FloatImageType>::getMinimumVoxelValue(edges2) <<
      std::endl;
   
   std::cout << 
      SmallUtilityMethods<FloatImageType>::getMaximumVoxelValue(edges2) <<
      std::endl;
      
   //VolumeIOWrapper<FloatImageType>::writeITKVolume16Bit( 
   //   edges2, "D:\\temp\\edges2.hdr", "" );
   
   const unsigned int neighborhoodRadius = 1;
   
   typedef itk::ConstNeighborhoodIterator<Int16ImageType> IterType;
   itk::ConstNeighborhoodIterator<Int16ImageType>::RadiusType radius;
   radius.Fill(neighborhoodRadius);
   
   IterType origIt1(radius, image1, image1->GetLargestPossibleRegion());
   IterType origIt2(radius, image2, image2->GetLargestPossibleRegion());
   
   itk::ImageRegionConstIterator<FloatImageType> it1(
      edges1,
      edges1->GetLargestPossibleRegion() );
   itk::ImageRegionConstIterator<FloatImageType> it2(
      edges2, 
      edges2->GetLargestPossibleRegion() );
   
   unsigned long nrVoxels = 0L;
   const unsigned int totalNbhBufferSize = origIt1.GetSize()[0] * 
      origIt1.GetSize()[1] * origIt1.GetSize()[2];
      
   while (!it1.IsAtEnd())
   {
      // look at all of the original pixel values in a 3x3x3 neighborhood,
      // if one of them is the warpingDefaultValue -> continue
      
      bool skip = false;
      
      for (unsigned int nb=0; nb<totalNbhBufferSize; ++nb)
      {
         if (origIt1.GetPixel(nb) == warpingDefaultValue ||
             origIt2.GetPixel(nb) == warpingDefaultValue)
         {
            skip = true;
         }
      }
      
      if (skip == false)
      {
         const FloatImageType::PixelType value1 = it1.Get();
         const FloatImageType::PixelType value2 = it2.Get();
         
         const double difference = 
            vcl_fabs( static_cast<double>(value1) - 
                      static_cast<double>(value2) );

         edgeOverlapMeasure += difference;
         ++nrVoxels;         
      }
      
      ++it1;
      ++it2;
      ++origIt1;
      ++origIt2;
   }
   
   edgeOverlapMeasure /= static_cast<double>(nrVoxels);
   
   std::cout << "edgeOverlapMeasure: " << edgeOverlapMeasure << std::endl;
   
   return edgeOverlapMeasure;
}

double calculateRobustMaxAbsoluteDeviationIntensity(
   const Int16ImageType::Pointer image1,
   const Int16ImageType::Pointer image2,
   const int warpingDefaultValue )
{
   itk::ImageRegionConstIterator<Int16ImageType> it1(
      image1,
      image1->GetLargestPossibleRegion() );
   itk::ImageRegionConstIterator<Int16ImageType> it2(
      image2, 
      image2->GetLargestPossibleRegion() );
   
   double robustMaxDeviationIntensity = 0.0;
   
   const Int16ImageType::SizeType imageSize = 
      image1->GetLargestPossibleRegion().GetSize();
   const unsigned long vectorSize = imageSize[0]*imageSize[1]*imageSize[2];
   
   std::cout << "imageSize: " << imageSize << std::endl; 
   std::cout << "vectorSize: " << vectorSize << std::endl;
   
   {
      std::vector<int> intensityDifferences;
      intensityDifferences.reserve(vectorSize);
      
      while (!it1.IsAtEnd())
      {
         const Int16ImageType::PixelType value1 = it1.Get();
         const Int16ImageType::PixelType value2 = it2.Get();
         {  
            // we ignore those values, where a synthetic transformation has 
            // shown the warpingDefaultValue. these values can never be estimated
            // by a registration process
            if (value1 != warpingDefaultValue && 
                value2 != warpingDefaultValue)
            {
               const int difference = static_cast<int>(
                  vcl_fabs( static_cast<double>(value1) - 
                            static_cast<double>(value2) ) );
   
               intensityDifferences.push_back(difference);
            }
         }
         
         ++it1;
         ++it2;
      }
      
      // calculate the robust max of the intensity differences
      const std::vector<int>::iterator robustMaxIterator = 
         intensityDifferences.begin() + 
         (intensityDifferences.size() / 100L) * 95L;
      std::nth_element(intensityDifferences.begin(), robustMaxIterator, 
         intensityDifferences.end());
      robustMaxDeviationIntensity = *robustMaxIterator;
   }
   
   
   std::cout << "robust max intensity difference: " << 
      robustMaxDeviationIntensity << std::endl;
   
   return robustMaxDeviationIntensity;
}

double calculateMedianAbsoluteDeviationIntensity(
   const Int16ImageType::Pointer image1,
   const Int16ImageType::Pointer image2,
   const int warpingDefaultValue )
{
   itk::ImageRegionConstIterator<Int16ImageType> it1(
      image1,
      image1->GetLargestPossibleRegion() );
   itk::ImageRegionConstIterator<Int16ImageType> it2(
      image2, 
      image2->GetLargestPossibleRegion() );
   
   double medianDeviationIntensity = 0.0;
   
   const Int16ImageType::SizeType imageSize = 
      image1->GetLargestPossibleRegion().GetSize();
   const unsigned long vectorSize = imageSize[0]*imageSize[1]*imageSize[2];
   
   std::cout << "imageSize: " << imageSize << std::endl; 
   std::cout << "vectorSize: " << vectorSize << std::endl;
   
   int medianDifference = 0;
   {
      std::vector<int> intensityDifferences;
      intensityDifferences.reserve(vectorSize);
      
      while (!it1.IsAtEnd())
      {
         const Int16ImageType::PixelType value1 = it1.Get();
         const Int16ImageType::PixelType value2 = it2.Get();
         {  
            // we ignore those values, where a synthetic transformation has 
            // shown the warpingDefaultValue. these values can never be estimated
            // by a registration process
            if (value1 != warpingDefaultValue && 
                value2 != warpingDefaultValue)
            {
               const int difference = static_cast<int>(
                  vcl_fabs( static_cast<double>(value1) - 
                            static_cast<double>(value2) ) );
   
               intensityDifferences.push_back(difference);
            }
         }
         
         ++it1;
         ++it2;
      }
      
      // calculate the median of the intensity differences
      const std::vector<int>::iterator medianIterator = 
         intensityDifferences.begin() + intensityDifferences.size() / 2L;
      std::nth_element(intensityDifferences.begin(), medianIterator, 
         intensityDifferences.end());
      medianDifference = *medianIterator;
   }
   
   
   std::cout << "median intensity difference: " << medianDifference << 
      std::endl;
   
   it1.GoToBegin();
   it2.GoToBegin();
   
   
   {
      std::vector<int> intensityDifferences;
      intensityDifferences.reserve(vectorSize);
      
      while (!it1.IsAtEnd())
      {
         const Int16ImageType::PixelType value1 = it1.Get();
         const Int16ImageType::PixelType value2 = it2.Get();
         {  
            // we ignore those values, where a synthetic transformation has 
            // shown the warpingDefaultValue. these values can never be 
            // estimated by a registration process
            if (value1 != warpingDefaultValue && 
                value2 != warpingDefaultValue)
            {
               const int difference = static_cast<int>( vcl_fabs(
                  vcl_fabs( static_cast<double>(value1) - 
                            static_cast<double>(value2) ) - 
                  static_cast<double>(medianDifference) ) );
   
               intensityDifferences.push_back(difference);
            }
         }
         
         ++it1;
         ++it2;
      }
      
      // calculate the median of the intensity differences
      const std::vector<int>::iterator medianIterator = 
         intensityDifferences.begin() + intensityDifferences.size() / 2L;
      std::nth_element(intensityDifferences.begin(), medianIterator, 
         intensityDifferences.end());
      medianDeviationIntensity = *medianIterator;
   }
   
   std::cout << "median deviation: " << medianDeviationIntensity << 
      std::endl;
   
   return medianDeviationIntensity;
}

double calculateMedianAbsoluteDeviationDF(
   const Int16ImageType::Pointer image1,
   const Int16ImageType::Pointer image2,
   const DeformationFieldType::Pointer deformationField1,
   const DeformationFieldType::Pointer deformationField2,
   const int warpingDefaultValue )
{
   itk::ImageRegionConstIterator<Int16ImageType> it1(
      image1,
      image1->GetLargestPossibleRegion() );
   itk::ImageRegionConstIterator<Int16ImageType> it2(
      image2, 
      image2->GetLargestPossibleRegion() );
   itk::ImageRegionConstIterator<DeformationFieldType> defIt1(
      deformationField1,
      deformationField1->GetLargestPossibleRegion() );
   itk::ImageRegionConstIterator<DeformationFieldType> defIt2(
      deformationField2, 
      deformationField2->GetLargestPossibleRegion() );
      
   double medianDeviationDF = 0.0;
   
   const Int16ImageType::SizeType imageSize = 
      deformationField1->GetLargestPossibleRegion().GetSize();
   const unsigned long vectorSize = 
      3*imageSize[0]*imageSize[1]*imageSize[2];
   
   std::cout << "imageSize: " << imageSize << std::endl; 
   std::cout << "vectorSize: " << vectorSize << std::endl;
   
   float medianDifference = 0.0f;
   {
      std::vector<float> dfDifferences;
      dfDifferences.reserve(vectorSize);
      
      while (!it1.IsAtEnd())
      {
         const Int16ImageType::PixelType value1 = it1.Get();
         const Int16ImageType::PixelType value2 = it2.Get();
         {  
            // we ignore those values, where a synthetic transformation has 
            // shown the warpingDefaultValue. these values can never be estimated
            // by a registration process
            if (value1 != warpingDefaultValue && 
                value2 != warpingDefaultValue)
            {
               dfDifferences.push_back( vcl_fabs(defIt1.Get()[0] - 
                                                 defIt2.Get()[0]) );
               dfDifferences.push_back( vcl_fabs(defIt1.Get()[1] - 
                                                 defIt2.Get()[1]) );
               dfDifferences.push_back( vcl_fabs(defIt1.Get()[2] - 
                                                 defIt2.Get()[2]) );
            }
         }
         
         ++it1;
         ++it2;
         ++defIt1;
         ++defIt2;
      }
      
      // calculate the median of the intensity differences
      const std::vector<float>::iterator medianIterator = 
         dfDifferences.begin() + dfDifferences.size() / 2L;
      std::nth_element(dfDifferences.begin(), medianIterator, 
         dfDifferences.end());
      medianDifference = *medianIterator;
   }
   
   
   std::cout << "median displacement field difference: " << 
      medianDifference << std::endl;
   
   it1.GoToBegin();
   it2.GoToBegin();
   defIt1.GoToBegin();
   defIt2.GoToBegin();
   
   {
      std::vector<float> dfDifferences;
      dfDifferences.reserve(vectorSize);
      
      while (!it1.IsAtEnd())
      {
         const Int16ImageType::PixelType value1 = it1.Get();
         const Int16ImageType::PixelType value2 = it2.Get();
         {  
            // we ignore those values, where a synthetic transformation has 
            // shown the warpingDefaultValue. these values can never be 
            // estimated by a registration process
            if (value1 != warpingDefaultValue && 
                value2 != warpingDefaultValue)
            {
               dfDifferences.push_back( 
                  vcl_fabs( vcl_fabs(defIt1.Get()[0] - defIt2.Get()[0]) -
                            medianDifference ) );
               dfDifferences.push_back( 
                  vcl_fabs( vcl_fabs(defIt1.Get()[1] - defIt2.Get()[1]) -
                            medianDifference ) );
               dfDifferences.push_back( 
                  vcl_fabs( vcl_fabs(defIt1.Get()[2] - defIt2.Get()[2]) -
                            medianDifference ) );
            }
         }
         
         ++it1;
         ++it2;
         ++defIt1;
         ++defIt2;
      }
      
      // calculate the median of the intensity differences
      const std::vector<float>::iterator medianIterator = 
         dfDifferences.begin() + dfDifferences.size() / 2L;
      std::nth_element(dfDifferences.begin(), medianIterator, 
         dfDifferences.end());
      medianDeviationDF = *medianIterator;
   }
   
   std::cout << "median deviation displacement field: " << 
      medianDeviationDF << std::endl;
   
   return medianDeviationDF;
}

double calculateMedianAbsoluteDeviationDF(
   const Int16ImageType::Pointer image1,
   const Int16ImageType::Pointer image2,
   const DeformationFieldType::Pointer deformationField1,
   const int warpingDefaultValue )
{
   itk::ImageRegionConstIterator<Int16ImageType> it1(
      image1,
      image1->GetLargestPossibleRegion() );
   itk::ImageRegionConstIterator<Int16ImageType> it2(
      image2, 
      image2->GetLargestPossibleRegion() );
   itk::ImageRegionConstIterator<DeformationFieldType> defIt1(
      deformationField1,
      deformationField1->GetLargestPossibleRegion() );
   
      
   double medianDeviationDF = 0.0;
   
   const Int16ImageType::SizeType imageSize = 
      deformationField1->GetLargestPossibleRegion().GetSize();
   const unsigned long vectorSize = 
      3*imageSize[0]*imageSize[1]*imageSize[2];
   
   std::cout << "imageSize: " << imageSize << std::endl; 
   std::cout << "vectorSize: " << vectorSize << std::endl;
   
   float medianDifference = 0.0f;
   {
      std::vector<float> dfDifferences;
      dfDifferences.reserve(vectorSize);
      
      while (!it1.IsAtEnd())
      {
         const Int16ImageType::PixelType value1 = it1.Get();
         const Int16ImageType::PixelType value2 = it2.Get();
         {  
            // we ignore those values, where a synthetic transformation has 
            // shown the warpingDefaultValue. these values can never be estimated
            // by a registration process
            if (value1 != warpingDefaultValue && 
                value2 != warpingDefaultValue)
            {
               dfDifferences.push_back( vcl_fabs(defIt1.Get()[0] - 
                                                 0.0) );
               dfDifferences.push_back( vcl_fabs(defIt1.Get()[1] - 
                                                 0.0) );
               dfDifferences.push_back( vcl_fabs(defIt1.Get()[2] - 
                                                 0.0) );
            }
         }
         
         ++it1;
         ++it2;
         ++defIt1;
      }
      
      // calculate the median of the intensity differences
      const std::vector<float>::iterator medianIterator = 
         dfDifferences.begin() + dfDifferences.size() / 2L;
      std::nth_element(dfDifferences.begin(), medianIterator, 
         dfDifferences.end());
      medianDifference = *medianIterator;
   }
   
   
   std::cout << "median displacement field difference: " << 
      medianDifference << std::endl;
   
   it1.GoToBegin();
   it2.GoToBegin();
   defIt1.GoToBegin();
   
   {
      std::vector<float> dfDifferences;
      dfDifferences.reserve(vectorSize);
      
      while (!it1.IsAtEnd())
      {
         const Int16ImageType::PixelType value1 = it1.Get();
         const Int16ImageType::PixelType value2 = it2.Get();
         {  
            // we ignore those values, where a synthetic transformation has 
            // shown the warpingDefaultValue. these values can never be 
            // estimated by a registration process
            if (value1 != warpingDefaultValue && 
                value2 != warpingDefaultValue)
            {
               dfDifferences.push_back( 
                  vcl_fabs( vcl_fabs(defIt1.Get()[0] - 0.0) -
                            medianDifference ) );
               dfDifferences.push_back( 
                  vcl_fabs( vcl_fabs(defIt1.Get()[1] - 0.0) -
                            medianDifference ) );
               dfDifferences.push_back( 
                  vcl_fabs( vcl_fabs(defIt1.Get()[2] - 0.0) -
                            medianDifference ) );
            }
         }
         
         ++it1;
         ++it2;
         ++defIt1;
      }
      
      // calculate the median of the intensity differences
      const std::vector<float>::iterator medianIterator = 
         dfDifferences.begin() + dfDifferences.size() / 2L;
      std::nth_element(dfDifferences.begin(), medianIterator, 
         dfDifferences.end());
      medianDeviationDF = *medianIterator;
   }
   
   std::cout << "median deviation displacement field: " << 
      medianDeviationDF << std::endl;
   
   return medianDeviationDF;
}


double calculateRobustMaximumDeviationDF(
   const Int16ImageType::Pointer image1,
   const Int16ImageType::Pointer image2,
   const DeformationFieldType::Pointer deformationField1,
   const DeformationFieldType::Pointer deformationField2,
   const int warpingDefaultValue )
{
   itk::ImageRegionConstIterator<Int16ImageType> it1(
      image1,
      image1->GetLargestPossibleRegion() );
   itk::ImageRegionConstIterator<Int16ImageType> it2(
      image2, 
      image2->GetLargestPossibleRegion() );
   itk::ImageRegionConstIterator<DeformationFieldType> defIt1(
      deformationField1,
      deformationField1->GetLargestPossibleRegion() );
   itk::ImageRegionConstIterator<DeformationFieldType> defIt2(
      deformationField2, 
      deformationField2->GetLargestPossibleRegion() );
      
   double robustMaximumDeviationDF = 0.0;
   
   const Int16ImageType::SizeType imageSize = 
      deformationField1->GetLargestPossibleRegion().GetSize();
   const unsigned long vectorSize = 
      3*imageSize[0]*imageSize[1]*imageSize[2];
   
   std::cout << "imageSize: " << imageSize << std::endl; 
   std::cout << "vectorSize: " << vectorSize << std::endl;
   
   
   {
      std::vector<float> dfDifferences;
      dfDifferences.reserve(vectorSize);
      
      while (!it1.IsAtEnd())
      {
         const Int16ImageType::PixelType value1 = it1.Get();
         const Int16ImageType::PixelType value2 = it2.Get();
         {  
            // we ignore those values, where a synthetic transformation has 
            // shown the warpingDefaultValue. these values can never be estimated
            // by a registration process
            if (value1 != warpingDefaultValue && 
                value2 != warpingDefaultValue)
            {
               dfDifferences.push_back( vcl_fabs(defIt1.Get()[0] - 
                                                 defIt2.Get()[0]) );
               dfDifferences.push_back( vcl_fabs(defIt1.Get()[1] - 
                                                 defIt2.Get()[1]) );
               dfDifferences.push_back( vcl_fabs(defIt1.Get()[2] - 
                                                 defIt2.Get()[2]) );
            }
         }
         
         ++it1;
         ++it2;
         ++defIt1;
         ++defIt2;
      }
      
      // calculate the median of the intensity differences
      const std::vector<float>::iterator medianIterator = 
         dfDifferences.begin() + (dfDifferences.size() / 100L) * 95L;
      std::nth_element(dfDifferences.begin(), medianIterator, 
         dfDifferences.end());
      robustMaximumDeviationDF = *medianIterator;
   }
   
   
   std::cout << "robust Maximum displacement field difference: " << 
      robustMaximumDeviationDF << std::endl;
   
   return robustMaximumDeviationDF;
}

double calculateRobustMaximumDeviationDF(
   const Int16ImageType::Pointer image1,
   const Int16ImageType::Pointer image2,
   const DeformationFieldType::Pointer deformationField1,
   const int warpingDefaultValue )
{
   itk::ImageRegionConstIterator<Int16ImageType> it1(
      image1,
      image1->GetLargestPossibleRegion() );
   itk::ImageRegionConstIterator<Int16ImageType> it2(
      image2, 
      image2->GetLargestPossibleRegion() );
   itk::ImageRegionConstIterator<DeformationFieldType> defIt1(
      deformationField1,
      deformationField1->GetLargestPossibleRegion() );
      
   double robustMaximumDeviationDF = 0.0;
   
   const Int16ImageType::SizeType imageSize = 
      deformationField1->GetLargestPossibleRegion().GetSize();
   const unsigned long vectorSize = 
      3*imageSize[0]*imageSize[1]*imageSize[2];
   
   std::cout << "imageSize: " << imageSize << std::endl; 
   std::cout << "vectorSize: " << vectorSize << std::endl;
   
   
   {
      std::vector<float> dfDifferences;
      dfDifferences.reserve(vectorSize);
      
      while (!it1.IsAtEnd())
      {
         const Int16ImageType::PixelType value1 = it1.Get();
         const Int16ImageType::PixelType value2 = it2.Get();
         {  
            // we ignore those values, where a synthetic transformation has 
            // shown the warpingDefaultValue. these values can never be estimated
            // by a registration process
            if (value1 != warpingDefaultValue && 
                value2 != warpingDefaultValue)
            {
               dfDifferences.push_back( vcl_fabs(defIt1.Get()[0] - 
                                                 0.0) );
               dfDifferences.push_back( vcl_fabs(defIt1.Get()[1] - 
                                                 0.0) );
               dfDifferences.push_back( vcl_fabs(defIt1.Get()[2] - 
                                                 0.0) );
            }
         }
         
         ++it1;
         ++it2;
         ++defIt1;
      }
      
      // calculate the median of the intensity differences
      const std::vector<float>::iterator medianIterator = 
         dfDifferences.begin() + (dfDifferences.size() / 100L) * 95L;
      std::nth_element(dfDifferences.begin(), medianIterator, 
         dfDifferences.end());
      robustMaximumDeviationDF = *medianIterator;
   }
   
   
   std::cout << "robust Maximum displacement field difference: " << 
      robustMaximumDeviationDF << std::endl;
   
   return robustMaximumDeviationDF;
}


double calculateNrZeroCrossingsOfJacobian(
   DeformationFieldType::Pointer field )
{
   NonlinearRegistrationUtilities::
      efficientInplaceDiscreteGaussSmoothDeformationField<DeformationFieldType>(
         field, 2.0, 2.0, 2.0 );
   
   itk::ConstNeighborhoodIterator<DeformationFieldType>::RadiusType radius;
   radius.Fill(1);
   
   itk::ConstNeighborhoodIterator<DeformationFieldType> fieldIter( radius, 
      field, field->GetLargestPossibleRegion() );
   
   
   float ux = 0.0f, uy = 0.0f, uz = 0.0f;
   float vx = 0.0f, vy = 0.0f, vz = 0.0f;
   float wx = 0.0f, wy = 0.0f, wz = 0.0f;
   float jacobian = 0.0f;
   
   unsigned long nrVoxels = 0L;
   unsigned long nrNegativeJacobians = 0L;
   
   while(!fieldIter.IsAtEnd())
   {
      int dim = 0;
      ux = 0.5f * (fieldIter.GetPixel(14)[dim] - fieldIter.GetPixel(12)[dim]);
      uy = 0.5f * (fieldIter.GetPixel(16)[dim] - fieldIter.GetPixel(10)[dim]);
      uz = 0.5f * (fieldIter.GetPixel(22)[dim] - fieldIter.GetPixel(4)[dim]);
      
      dim = 1;
      vx = 0.5f * (fieldIter.GetPixel(14)[dim] - fieldIter.GetPixel(12)[dim]);
      vy = 0.5f * (fieldIter.GetPixel(16)[dim] - fieldIter.GetPixel(10)[dim]);
      vz = 0.5f * (fieldIter.GetPixel(22)[dim] - fieldIter.GetPixel(4)[dim]);
      
      dim = 2;
      wx = 0.5f * (fieldIter.GetPixel(14)[dim] - fieldIter.GetPixel(12)[dim]);
      wy = 0.5f * (fieldIter.GetPixel(16)[dim] - fieldIter.GetPixel(10)[dim]);
      wz = 0.5f * (fieldIter.GetPixel(22)[dim] - fieldIter.GetPixel(4)[dim]);
      
      jacobian =  (1.0f+ux)*(1.0f+vy)*(1.0f+wz)+uy*vz*wx + uz*vx*wy -
                  ( (1.0f+ux)*vz*wy + (1.0f+vy)*uz*wx + (1.0f+wz)*uy*vx );
      
      if (jacobian < 0.0f)
      {
         ++nrNegativeJacobians;
      }
      
      ++nrVoxels;
      
      ++fieldIter;
   }
   
   std::cout << "negative jacobian percentage: " << 
      static_cast<double>(nrNegativeJacobians) / static_cast<double>(nrVoxels) 
      << std::endl;
   
   return static_cast<double>(nrNegativeJacobians) / 
          static_cast<double>(nrVoxels);
}

int main( int argc, char *argv[] )
{
   if (!(argc == 7 || argc == 9))
   {
      std::cerr << "Usage: " << std::endl;
      std::cerr << argv[0] <<
         "  original_fixed_image  original_moving_image" <<
         "  moving_image_warped_to_fixed" <<
         "  gold_standard_deformation_field  calculated_deformation_field" <<
         "  resample_default_value" <<
         "  [intensity_result_file.csv displacement_result_file.csv]"
         << std::endl;
      return EXIT_FAILURE;
   }
   
   bool writeCSVFile = false;
   if (argc == 9 )
      writeCSVFile = true;
   
   std::string original_fixed_image_filename( argv[1] );
   std::string original_moving_image_filename( argv[2] );
   std::string unwarped_image_filename( argv[3] );
   std::string gold_standard_deformation_field_filename( argv[4] );
   std::string calculated_deformation_field_filename( argv[5] );
   
   std::string result_intensity_statistics_filename;
   std::string result_displacement_statistics_filename;
   if (writeCSVFile)
   {
      result_intensity_statistics_filename = std::string( argv[7] );
      result_displacement_statistics_filename = std::string( argv[8] );
   }
   
   
   const Int16ImageType::PixelType warpingDefaultValue = atoi( argv[6] ); 
   std::cout << "warpingDefaultValue: " << warpingDefaultValue << std::endl;


   Int16ImageType::Pointer original_fixed_image = 
      VolumeIOWrapper<Int16ImageType>::readITKVolume( 
         original_fixed_image_filename, "");

   Int16ImageType::Pointer original_moving_image = 
      VolumeIOWrapper<Int16ImageType>::readITKVolume( 
         original_moving_image_filename, "");

   Int16ImageType::Pointer unwarped_image = 
      VolumeIOWrapper<Int16ImageType>::readITKVolume( 
         unwarped_image_filename, "");
         
   
   if (original_fixed_image->GetLargestPossibleRegion().GetSize() !=
       unwarped_image->GetLargestPossibleRegion().GetSize())
   {
      std::cerr << "images have to be of the same size!" 
         << std::endl;
      return EXIT_FAILURE;
   }
   
   // read the deformation fields
   DeformationFieldType::Pointer calculated_deformation_field = 0;
   DeformationFieldType::Pointer gold_standard_deformation_field = 0;
   
   typedef itk::ImageFileReader<DeformationFieldType> 
      DeformationFieldReaderType;
   {
      DeformationFieldReaderType::Pointer reader = 
         DeformationFieldReaderType::New();
      reader->SetFileName( calculated_deformation_field_filename.c_str() );
      reader->Update();
   
      calculated_deformation_field = reader->GetOutput();
   }
   {
      DeformationFieldReaderType::Pointer reader = 
         DeformationFieldReaderType::New();
      reader->SetFileName( gold_standard_deformation_field_filename.c_str() );
      reader->Update();
   
      gold_standard_deformation_field = reader->GetOutput();
   }
   
   // check sizes
   if (calculated_deformation_field->GetLargestPossibleRegion().GetSize() !=
       gold_standard_deformation_field->GetLargestPossibleRegion().GetSize())
   {
      std::cout << "deformation fields have to be of the same size!" 
         << std::endl;
      return EXIT_FAILURE;
   }
   
   if (original_fixed_image->GetLargestPossibleRegion().GetSize() !=
       calculated_deformation_field->GetLargestPossibleRegion().GetSize())
   {
      std::cout << "images and displacement fields have to be of the same size!" 
         << std::endl;
      return EXIT_FAILURE;
   }
   
   
   // HERE THE VALIDATION STARTS
   
   // values that lie above this threshold are clamped to this value
   // this should make the evaluation more robust, since large differences
   // are mismatches, no matter how large they are
   // in addition this makes it possible to deal with resampleDefaultValue
   // in the unwarped image! if we meet one of them, we simply assign
   // this clampThreshold to the difference
   const double intensity_difference_clamp_threshold = 100.0;
   
   itk::ImageRegionConstIterator<Int16ImageType> fixed_it(
      original_fixed_image,
      original_fixed_image->GetLargestPossibleRegion() );
   itk::ImageRegionConstIterator<Int16ImageType> moving_it(
      original_moving_image, 
      original_moving_image->GetLargestPossibleRegion() );
   itk::ImageRegionConstIterator<Int16ImageType> unwarp_it(
      unwarped_image, 
      unwarped_image->GetLargestPossibleRegion() );

   itk::ImageRegionConstIterator<DeformationFieldType> def_it( 
      calculated_deformation_field, 
      calculated_deformation_field->GetLargestPossibleRegion() );
   itk::ImageRegionConstIterator<DeformationFieldType> gs_def_it( 
      gold_standard_deformation_field, 
      gold_standard_deformation_field->GetLargestPossibleRegion() );


   // the metric values
   long numberOfMeanSquaredIntensityVoxels = 0L;
   long numberOfMeanSquaredIntensityVoxelsUnwarped = 0L;
   long numberOfDisplacementVoxels = 0L;
   
   double mean_intensity_difference = 0.0;
   double std_intensity_difference = 0.0;
   double rms_intensity_difference = 0.0;
   double min_intensity_difference = vcl_numeric_limits<double>::max();
   double max_intensity_difference = vcl_numeric_limits<double>::min();
   
   double mean_intensity_difference_unwarped = 0.0;
   double std_intensity_difference_unwarped = 0.0;
   double rms_intensity_difference_unwarped = 0.0;
   double min_intensity_difference_unwarped = vcl_numeric_limits<double>::max();
   double max_intensity_difference_unwarped = vcl_numeric_limits<double>::min();
   
   // the displacement metric values
   double mean_displacement_difference_vector[3] = {0.0, 0.0, 0.0};
   double mean_absolute_displacement_difference_vector[3] = {0.0, 0.0, 0.0};
   
   //double rms_displacement = 0.0;
   //double mean_displacement = 0.0;
   //double min_displacement = vcl_numeric_limits<double>::max();
   //double max_displacement = vcl_numeric_limits<double>::min();
   
   double rms_displacement_field = 0.0;
   double max_displacement_field = vcl_numeric_limits<double>::min();

   double rms_displacement_field_ground_truth = 0.0;
   double max_displacement_field_ground_truth = vcl_numeric_limits<double>::min();
   
   VectorPixelType undefinedDef = 
      VectorPixelConstants<VectorPixelType::ValueType, 3>::undefined();
   
   while (!fixed_it.IsAtEnd())
   {
      const Int16ImageType::PixelType fixedValue = fixed_it.Get();
      const Int16ImageType::PixelType movingValue = moving_it.Get();
      {  
         // we ignore those values, where a synthetic transformation has 
         // shown the warpingDefaultValue. these values can never be estimated
         // by a registration process
         if (fixedValue != warpingDefaultValue && 
             movingValue != warpingDefaultValue)
         {
            const double difference = 
               vcl_fabs( static_cast<double>(fixedValue) - 
                         static_cast<double>(movingValue) );

            const double clamped_difference = 
               (difference > intensity_difference_clamp_threshold) ?
                  intensity_difference_clamp_threshold : difference;

            if (difference < min_intensity_difference)
               min_intensity_difference = difference;
            if (difference > max_intensity_difference)
               max_intensity_difference = difference;
               
            mean_intensity_difference += clamped_difference;
            rms_intensity_difference += vnl_math_sqr(clamped_difference);
            
            ++numberOfMeanSquaredIntensityVoxels;
         }
      }
      
      const Int16ImageType::PixelType unwarpedValue = unwarp_it.Get();
      {
         // we ignore those values, where a synthetic transformation has 
         // shown the warpingDefaultValue. these values can never be estimated
         // by a registration process
         if (fixedValue != warpingDefaultValue)
         {
            const double difference = 
               vcl_fabs( static_cast<double>(fixedValue) - 
                         static_cast<double>(unwarpedValue) );

            const double clamped_difference = 
               (unwarpedValue == warpingDefaultValue || 
                difference > intensity_difference_clamp_threshold) ?
                  intensity_difference_clamp_threshold : difference;

            if (difference < min_intensity_difference_unwarped)   
               min_intensity_difference_unwarped = difference;
            if (difference > max_intensity_difference_unwarped)
               max_intensity_difference_unwarped = difference;
            
            mean_intensity_difference_unwarped += clamped_difference;   
            rms_intensity_difference_unwarped += vnl_math_sqr(clamped_difference);
         
            ++numberOfMeanSquaredIntensityVoxelsUnwarped;
         }
      }
      
      // evaluate the displacements
      if ( movingValue != warpingDefaultValue &&
           fixedValue  != warpingDefaultValue )
      {
         VectorPixelType gs_def = gs_def_it.Get();
         VectorPixelType def = def_it.Get();
         
         
         
         if (def == undefinedDef)
         {
            std::cout << "unexpected undefined calculated deformation!" << 
               std::endl;
            def = VectorPixelConstants<VectorPixelType::ValueType,3>::zero();
         }
         if (gs_def == undefinedDef)
         {
            std::cout << "unexpected undefined gs deformation!" << std::endl;
            gs_def = VectorPixelConstants<VectorPixelType::ValueType,3>::zero();
         }
         
         if (def[0] > 1e7 || def[1] > 1e7 || def[2] > 1e7)
         {
            std::cout << "def: " << def << std::endl;
            def = VectorPixelConstants<VectorPixelType::ValueType,3>::zero();
         }
         
         mean_displacement_difference_vector[0] += (gs_def[0] - def[0]);
         mean_displacement_difference_vector[1] += (gs_def[1] - def[1]);
         mean_displacement_difference_vector[2] += (gs_def[2] - def[2]);
         
         mean_absolute_displacement_difference_vector[0] += 
            fabs(gs_def[0] - def[0]);
         mean_absolute_displacement_difference_vector[1] += 
            fabs(gs_def[1] - def[1]);
         mean_absolute_displacement_difference_vector[2] += 
            fabs(gs_def[2] - def[2]);
         
         rms_displacement_field += vnl_math_sqr(gs_def[0] - def[0]);
         rms_displacement_field += vnl_math_sqr(gs_def[1] - def[1]);
         rms_displacement_field += vnl_math_sqr(gs_def[2] - def[2]);
         
         if (fabs(gs_def[0] - def[0]) > max_displacement_field)
            max_displacement_field = fabs(gs_def[0] - def[0]);
         if (fabs(gs_def[1] - def[1]) > max_displacement_field)
            max_displacement_field = fabs(gs_def[1] - def[1]);
         if (fabs(gs_def[2] - def[2]) > max_displacement_field)
            max_displacement_field = fabs(gs_def[2] - def[2]);

         rms_displacement_field_ground_truth += vnl_math_sqr(gs_def[0] - 0.0);
         rms_displacement_field_ground_truth += vnl_math_sqr(gs_def[1] - 0.0);
         rms_displacement_field_ground_truth += vnl_math_sqr(gs_def[2] - 0.0);
         
         if (fabs(gs_def[0] - 0.0) > max_displacement_field_ground_truth)
            max_displacement_field_ground_truth = fabs(gs_def[0] - 0.0);
         if (fabs(gs_def[1] - 0.0) > max_displacement_field_ground_truth)
            max_displacement_field_ground_truth = fabs(gs_def[1] - 0.0);
         if (fabs(gs_def[2] - 0.0) > max_displacement_field_ground_truth)
            max_displacement_field_ground_truth = fabs(gs_def[2] - 0.0);
         
//          double normOfDisplacementDifference =
//             vcl_sqrt( vnl_math_sqr(gs_def[0] - def[0]) + 
//                       vnl_math_sqr(gs_def[1] - def[1]) +
//                       vnl_math_sqr(gs_def[2] - def[2]) );
                  
//          if (normOfDisplacementDifference < min_displacement)
//             min_displacement = normOfDisplacementDifference;
//          if (normOfDisplacementDifference > max_displacement)
//             max_displacement = normOfDisplacementDifference;
            
//          mean_displacement += normOfDisplacementDifference;
//          rms_displacement += vnl_math_sqr(normOfDisplacementDifference);
         
         ++numberOfDisplacementVoxels;
      }
      
      ++fixed_it;
      ++moving_it;
      ++unwarp_it;
      ++def_it;
      ++gs_def_it;
   }

   // never change the calculations of these three numbers, the std depends
   // on this order! 
   // we calculate intensity difference standard deviation using 
   // running statistics!
   mean_intensity_difference /= (double) numberOfMeanSquaredIntensityVoxels;
   std_intensity_difference = 
      sqrt(1.0 / ((double) (numberOfMeanSquaredIntensityVoxels - 1.0)) *
      (rms_intensity_difference - numberOfMeanSquaredIntensityVoxels * 
      mean_intensity_difference * mean_intensity_difference));
   rms_intensity_difference = sqrt( rms_intensity_difference / 
      (double)numberOfMeanSquaredIntensityVoxels );

   
   // never change the calculations of these three numbers, the std depends
   // on this order!
   // we calculate intensity difference standard deviation using running 
   // statistics!
   mean_intensity_difference_unwarped /= 
      (double) numberOfMeanSquaredIntensityVoxelsUnwarped;
   std_intensity_difference_unwarped = 
      sqrt(1.0 / ((double) (numberOfMeanSquaredIntensityVoxelsUnwarped - 1.0)) *
      (rms_intensity_difference_unwarped - 
         numberOfMeanSquaredIntensityVoxelsUnwarped * 
      mean_intensity_difference_unwarped * mean_intensity_difference_unwarped));
   rms_intensity_difference_unwarped = sqrt(rms_intensity_difference_unwarped / 
      (double) numberOfMeanSquaredIntensityVoxelsUnwarped );
      
   // displacement measures
      
   //mean_displacement /= (double) numberOfDisplacementVoxels;
   //rms_displacement = sqrt( rms_displacement /
   //   (double) numberOfDisplacementVoxels );

   rms_displacement_field_ground_truth = sqrt( rms_displacement_field_ground_truth / 
      (3.0 * numberOfDisplacementVoxels) );
   
   rms_displacement_field = sqrt( rms_displacement_field / 
      (3.0 * numberOfDisplacementVoxels) );
   
   mean_displacement_difference_vector[0] /= 
      (double) numberOfDisplacementVoxels;
   mean_displacement_difference_vector[1] /= 
      (double) numberOfDisplacementVoxels;
   mean_displacement_difference_vector[2] /= 
      (double) numberOfDisplacementVoxels;
   
   mean_absolute_displacement_difference_vector[0] /= 
      (double) numberOfDisplacementVoxels;
   mean_absolute_displacement_difference_vector[1] /= 
      (double) numberOfDisplacementVoxels;
   mean_absolute_displacement_difference_vector[2] /= 
      (double) numberOfDisplacementVoxels;

   
   const double medianDeviationIntensity = 
      calculateMedianAbsoluteDeviationIntensity(original_fixed_image,
         original_moving_image, warpingDefaultValue);

   const double medianDeviationIntensityUnwarped = 
      calculateMedianAbsoluteDeviationIntensity(original_fixed_image,
         unwarped_image, warpingDefaultValue);
   
   const double robustMaxDeviationIntensity = 
      calculateRobustMaxAbsoluteDeviationIntensity(original_fixed_image,
         original_moving_image, warpingDefaultValue);

   const double robustMaxDeviationIntensityUnwarped = 
      calculateRobustMaxAbsoluteDeviationIntensity(original_fixed_image,
         unwarped_image, warpingDefaultValue);
   
   // edge overlap comparison
   const double edgeOverlap = calculateEdgeOverlapMeasure(original_fixed_image,
      original_moving_image, warpingDefaultValue);
      
   const double edgeOverlapUnwarped = calculateEdgeOverlapMeasure(
      original_fixed_image, unwarped_image, warpingDefaultValue);
   
   const double medianDeviationDisplacementField1 = 
      calculateMedianAbsoluteDeviationDF(original_fixed_image,
         original_moving_image, gold_standard_deformation_field,
         warpingDefaultValue);

   const double medianDeviationDisplacementField2 = 
      calculateMedianAbsoluteDeviationDF(original_fixed_image,
         original_moving_image, gold_standard_deformation_field,
         calculated_deformation_field, warpingDefaultValue);
   
   const double robustMaximumDeviationDisplacementField1 = 
      calculateRobustMaximumDeviationDF(original_fixed_image,
         original_moving_image, gold_standard_deformation_field,
         warpingDefaultValue);

   const double robustMaximumDeviationDisplacementField2 = 
      calculateRobustMaximumDeviationDF(original_fixed_image,
         original_moving_image, gold_standard_deformation_field,
         calculated_deformation_field, warpingDefaultValue);
   
   const double nrZeroCrossings1 = 
      calculateNrZeroCrossingsOfJacobian(gold_standard_deformation_field);
   
   const double nrZeroCrossings2 = 
      calculateNrZeroCrossingsOfJacobian(calculated_deformation_field);
      
   // print statistics:
   
   std::cout << "numberOfMeanSquaredIntensityVoxels: " << 
      numberOfMeanSquaredIntensityVoxels << std::endl;
   std::cout << "Mean intensity difference: " << 
      mean_intensity_difference << std::endl;
   std::cout << "STD intensity difference: " << 
      std_intensity_difference << std::endl;
   std::cout << "RMS intensity difference: " << 
      rms_intensity_difference << std::endl;
   //std::cout << "min intensity difference: " << min_intensity_difference << 
   //   std::endl;
   std::cout << "max intensity difference: " << max_intensity_difference << 
      std::endl << std::endl;
      
   std::cout << "numberOfMeanSquaredIntensityVoxelsUnwarped: " << 
      numberOfMeanSquaredIntensityVoxelsUnwarped << std::endl;
   std::cout << "Mean intensity difference unwarped: " << 
      mean_intensity_difference_unwarped << std::endl;
   std::cout << "STD intensity difference unwarped: " << 
      std_intensity_difference_unwarped << std::endl;
   std::cout << "RMS intensity difference unwarped: " << 
      rms_intensity_difference_unwarped << std::endl;
   //std::cout << "min intensity difference unwarped: " << 
   //   min_intensity_difference_unwarped << std::endl;
   std::cout << "max intensity difference unwarped: " << 
      max_intensity_difference_unwarped << std::endl << std::endl;
   
   std::cout << "numberOfDisplacementVoxels: " << 
      numberOfDisplacementVoxels << std::endl;
      
//    std::cout << "mean displacement difference: " << 
//       mean_displacement << std::endl;
//    std::cout << "rms displacement difference: " << 
//       rms_displacement << std::endl;
//    std::cout << "min displacement difference: " << 
//       min_displacement << std::endl;
//    std::cout << "max displacement difference: " << 
//       max_displacement << std::endl << std::endl;
   
   std::cout << "rms displacement field ground truth: " << rms_displacement_field_ground_truth << 
      std::endl;
   std::cout << "max displacement field ground truth: " << max_displacement_field_ground_truth << 
      std::endl;

   std::cout << "rms displacement field: " << rms_displacement_field << 
      std::endl;
   std::cout << "max displacement field: " << max_displacement_field << 
      std::endl;
   
   std::cout << "mean displacement difference vector: " << 
      mean_displacement_difference_vector[0] << ":" << 
      mean_displacement_difference_vector[1] << ":" << 
      mean_displacement_difference_vector[2] << ":" << std::endl;
   std::cout << " -> its length: " <<
      vcl_sqrt( vnl_math_sqr(mean_displacement_difference_vector[0]) +
                vnl_math_sqr(mean_displacement_difference_vector[1]) +
                vnl_math_sqr(mean_displacement_difference_vector[2]) ) 
                << std::endl;
   std::cout << "mean absolute displacement difference vector: " << 
      mean_absolute_displacement_difference_vector[0] << ":" << 
      mean_absolute_displacement_difference_vector[1] << ":" << 
      mean_absolute_displacement_difference_vector[2] << ":" << std::endl;
   std::cout << " -> its length: " <<
      vcl_sqrt( vnl_math_sqr(mean_absolute_displacement_difference_vector[0]) +
                vnl_math_sqr(mean_absolute_displacement_difference_vector[1]) +
                vnl_math_sqr(mean_absolute_displacement_difference_vector[2]) ) 
                << std::endl;
                
                
                
   // now for some mutual information statistics
   
   double nmi = 0.0;
   double nmi_unwarped = 0.0;
   
   calculateMutualInformationStatistics( original_fixed_image, 
      original_moving_image, nmi );
   calculateMutualInformationStatistics( original_fixed_image,
      unwarped_image, nmi_unwarped );
   
   if (writeCSVFile)
   {
      
      // open the result file
      std::cout << "trying to open file " << 
         result_intensity_statistics_filename.c_str() << std::endl;
      std::ofstream result_file;
      try
      {
         result_file.open( result_intensity_statistics_filename.c_str(), 
            std::fstream::app ); 
   
         if (result_file.is_open())
         {
            result_file << original_fixed_image_filename << ";" <<
                           original_moving_image_filename << ";" <<
                           rms_intensity_difference << ";" <<
                           std_intensity_difference << ";" <<
                           min_intensity_difference << ";" <<
                           robustMaxDeviationIntensity << ";" <<
                           nmi << ";" <<
                           medianDeviationIntensity << ";" <<
                           edgeOverlap << ";" << 
                           std::endl;
      
            result_file << original_fixed_image_filename << ";" <<
                           unwarped_image_filename << ";" <<
                           rms_intensity_difference_unwarped << ";" <<
                           std_intensity_difference_unwarped << ";" <<
                           min_intensity_difference_unwarped << ";" <<
                           robustMaxDeviationIntensityUnwarped << ";" <<
                           nmi_unwarped << ";" <<
                           medianDeviationIntensityUnwarped << ";" <<
                           edgeOverlapUnwarped << ";" << 
                           std::endl;
                     
            result_file.close();
         }
         else
         {
            std::cerr << "file couldn't be opened!" << std::endl;
         }
      }
      catch( std::exception& )
      {
         std::cerr << "file io exception!" << std::endl;
      }
      
      // open the result file
      std::cout << "trying to open file " << 
         result_displacement_statistics_filename.c_str() << std::endl;
      std::ofstream result_file2;
      try
      {
         result_file2.open( result_displacement_statistics_filename.c_str(), 
            std::fstream::app ); 
   
         if (result_file2.is_open())
         {
            result_file2 << gold_standard_deformation_field_filename << ";" <<
                            "0" << ";" <<
                            rms_displacement_field_ground_truth << ";" <<
                            robustMaximumDeviationDisplacementField1 << ";" <<
                            medianDeviationDisplacementField1 << ";" << 
                            nrZeroCrossings1 << ";" << std::endl;
      
            result_file2 << gold_standard_deformation_field_filename << ";" <<
                            calculated_deformation_field_filename << ";" <<
                            rms_displacement_field << ";" <<
                            robustMaximumDeviationDisplacementField2 << ";" <<
                            medianDeviationDisplacementField2 << ";" <<
                            nrZeroCrossings2 << ";" << std::endl;
                     
            result_file2.close();
         }
         else
         {
            std::cerr << "file couldn't be opened!" << std::endl;
         }
      }
      catch( std::exception& )
      {
         std::cerr << "file io exception!" << std::endl;
      }
   }
   
   return EXIT_SUCCESS;
}