/**********************************
*evaluation.cxx
*Evaluation Software
*
*This program will compare the "gold standard" deformation
*with the estimated error.  
*
*
*Input:
*      deformation norm image
*      error estimation image
*
*
*output:  % of outliers
*  An outlier is defined as a point whose value is different 
*  from "gold standard" by 2mm or more.
*         
* Authors: Eric Billet, Andriy Fedorov 
* PI: Nikos Chrisochoides
***********************************/



#include "itkImage.h"
#include "itkImageRegionIterator.h"
#include "itkImageFileReader.h"


	typedef float        InputPixelType;
	typedef itk::Image< InputPixelType, 3 >    InputImageType;
	typedef itk::ImageFileReader< InputImageType  >  ReaderType;
	typedef itk::ImageRegionConstIterator< InputImageType>        IteratorTypeC;

int main( int argc, char ** argv )
{
	
  // Verify the number of parameters in the command line
  if( argc < 4 )
    {
      std::cerr << "Usage: " << std::endl;
      std::cerr << argv[0] << "[ldmap][df_norm][edges1][edges2]" << std::endl;
      return EXIT_FAILURE;
    }

  
  ReaderType::Pointer readerDF = ReaderType::New();
  readerDF->SetFileName(argv[2]);
  readerDF->Update();
  ReaderType::Pointer reader_ld_map = ReaderType::New();
  reader_ld_map->SetFileName(argv[1]);
  reader_ld_map->Update();
  ReaderType::Pointer reader_edges2 = ReaderType::New();
  reader_edges2->SetFileName(argv[4]);
  reader_edges2->Update();
  ReaderType::Pointer reader_edges1 = ReaderType::New();
  reader_edges1->SetFileName(argv[3]);
  reader_edges1->Update();



  IteratorTypeC F_edgeIt(reader_edges1->GetOutput(), reader_edges1->GetOutput()->GetRequestedRegion());
  IteratorTypeC G_edgeIt(reader_edges2->GetOutput(), reader_edges2->GetOutput()->GetRequestedRegion()); 
  IteratorTypeC LDIt(reader_ld_map->GetOutput(), reader_ld_map->GetOutput()->GetRequestedRegion());
  IteratorTypeC DFmapIt(readerDF->GetOutput(), readerDF->GetOutput()->GetRequestedRegion());

  int outlier=0,count=0,under_outlier=0,under=0;
		
   //iterate through images
  for(F_edgeIt.GoToBegin(),G_edgeIt.GoToBegin(), LDIt.GoToBegin(), DFmapIt.GoToBegin(); !F_edgeIt.IsAtEnd(); ++F_edgeIt,++G_edgeIt, ++LDIt, ++DFmapIt){
    if((F_edgeIt.Get() || G_edgeIt.Get()) && DFmapIt.Get()){

      double dist = DFmapIt.Get();
      double est = LDIt.Get();
      if(est >=0){
            //difference is the accuracy of the local distance map	
	 double diff = fabs(est-fabs(dist));

	//compute outlier %
	if(diff>2){
	  outlier++;
          if(fabs(dist)>est)
	    under_outlier++;
	}
        if(fabs(dist)>est)
          under++;	
	count++;
      }
    }
  }

 std::cout<<"outlier underestimate underestimate_outliers"<<std::endl;
 std::cout<<double(outlier)/double(count)<<" ";
 std::cout<<double(under)/double(count)<<" ";
 std::cout<<double(under_outlier)/double(outlier)<<" ";

return 0;	
}