text: #include <iostream>
  #include <fstream>
  #include <cmath>
  #include <vector>
  #include <string>
  #include <iomanip>
  
  // VTK include
  #include <vtkPolyDataReader.h>
  #include <vtkXMLPolyDataReader.h>
  #include <vtkPolyData.h>
  #include <vtkSmartPointer.h>
  #include <vtkPointData.h>
  #include <vtkSpline.h>
  #include <vtkParametricSpline.h>
  #include <vtkCellArray.h>
  #include <vtkPoints.h>
  #include <vtkSplineRepresentation.h>
  #include <vtkMath.h>
  #include <vtkParametricFunctionSource.h>
  #include <vtkCardinalSpline.h>
  #include <vtkSCurveSpline.h>
  #include <vtkObject.h>
  #include <vtkPolyLine.h>
  
  // ITK include
  #include <itkImage.h>
  #include <itkVectorContainer.h>
  #include <itkVector.h>
  #include <itkDiffusionTensor3D.h>
  #include <itkImageAdaptor.h>
  #include <itkPoint.h>
  #include <itkImageFileReader.h>
  #include <itkImageFileWriter.h>
  #include <itkContinuousIndex.h>
  #include <itkSpatialObject.h>
  #include <itkLabelOverlapMeasuresImageFilter.h>
  #include <itkKappaStatisticImageToImageMetric.h>
  #include "itkTranslationTransform.h"
  #include "itkNearestNeighborInterpolateImageFunction.h"
  #include "itkImageRegionIterator.h"
  
  // VNL Includes
  #include <vnl/vnl_matrix.h>
  #include <vnl/vnl_vector_fixed.h>
  
  #include <CurveCompareCLP.h>
  
  using namespace std;
  
  std::vector<double> FillOverlapTable(vtkSmartPointer<vtkPolyData> Source, vtkSmartPointer<vtkPolyData> Target,string   ReferenceScalarVolume,int voxelLabel)
  {  std::cout<<"Overlap calculating..."<<std::endl;
    std::vector<double> OverlapTable;
    vtkPoints* SourcePoints = Source->GetPoints();
    vtkPoints* TargetPoints = Target->GetPoints();
    int numberPointsSource=SourcePoints->GetNumberOfPoints();
    int numberPointsTarget=TargetPoints->GetNumberOfPoints();
  
    constexpr unsigned int DIM = 3;
    typedef unsigned short ScalarPixelType;
    typedef itk::Image<ScalarPixelType, DIM> IntImageType;
    /*typedef itk::DiffusionTensor3D<double> TensorPixelType;
    typedef itk::Image<TensorPixelType, DIM> TensorImageType;
    // Setup tensor file if available
    typedef itk::ImageFileReader<TensorImageType> TensorImageReader;
    TensorImageReader::Pointer tensorreader = nullptr;
    tensorreader = TensorImageReader::New();
    tensorreader->SetFileName(tensorVolume.c_str());
    tensorreader->Update();*/
    typedef itk::ImageFileReader<IntImageType> ScalarImageReader;
    ScalarImageReader::Pointer scalarReader = nullptr;
    scalarReader = ScalarImageReader::New();
    scalarReader->SetFileName(ReferenceScalarVolume.c_str());
    scalarReader->Update();
  
    // Need to allocate two images to write into for creating
    // the fiber label maps
    IntImageType::Pointer labelimage1;//classic voxelization of Source
    IntImageType::Pointer labelimage2;//classic voxelization of Target
    IntImageType::Pointer labelimageCountF;//voxelization of Source by counting how many fibers for each voxel
    IntImageType::Pointer labelimage2CountF;//voxelization of Target by counting how many fibers for each voxel
  
    labelimage1 = IntImageType::New();
    labelimage1->SetSpacing(scalarReader->GetOutput()->GetSpacing());
    labelimage1->SetOrigin(scalarReader->GetOutput()->GetOrigin());
    labelimage1->SetDirection(scalarReader->GetOutput()->GetDirection());
    labelimage1->SetRegions(scalarReader->GetOutput()->GetLargestPossibleRegion());
    labelimage1->Allocate();
    labelimage1->FillBuffer(0);
    labelimageCountF = IntImageType::New();
    labelimageCountF->SetSpacing(scalarReader->GetOutput()->GetSpacing());
    labelimageCountF->SetOrigin(scalarReader->GetOutput()->GetOrigin());
    labelimageCountF->SetDirection(scalarReader->GetOutput()->GetDirection());
    labelimageCountF->SetRegions(scalarReader->GetOutput()->GetLargestPossibleRegion());
    labelimageCountF->Allocate();
    labelimageCountF->FillBuffer(0);
    labelimage2 = IntImageType::New();
    labelimage2->SetSpacing(scalarReader->GetOutput()->GetSpacing());
    labelimage2->SetOrigin(scalarReader->GetOutput()->GetOrigin());
    labelimage2->SetDirection(scalarReader->GetOutput()->GetDirection());
    labelimage2->SetRegions(scalarReader->GetOutput()->GetLargestPossibleRegion());
    labelimage2->Allocate();
    labelimage2->FillBuffer(0);
    labelimage2CountF = IntImageType::New();
    labelimage2CountF->SetSpacing(scalarReader->GetOutput()->GetSpacing());
    labelimage2CountF->SetOrigin(scalarReader->GetOutput()->GetOrigin());
    labelimage2CountF->SetDirection(scalarReader->GetOutput()->GetDirection());
    labelimage2CountF->SetRegions(scalarReader->GetOutput()->GetLargestPossibleRegion());
    labelimage2CountF->Allocate();
    labelimage2CountF->FillBuffer(0);
  
  
    typedef itk::Point<double, 3> PointType;
    PointType fiberpoint;
    double fiberpointtemp[3];
    for(int i=0;i<SourcePoints->GetNumberOfPoints();i++)
    {	
        SourcePoints->GetPoint( i, fiberpointtemp );
        for (int k=0;k<3;k++)
  	{
  		fiberpoint[k]=fiberpointtemp[k];
  	}
        itk::ContinuousIndex<double,3> cind;
        itk::Index<3> ind;
        labelimage1->TransformPhysicalPointToContinuousIndex(fiberpoint, cind);
        ind[0] = static_cast<long int>(itk::Math::rnd_halfinttoeven(cind[0]));
        ind[1] = static_cast<long int>(itk::Math::rnd_halfinttoeven(cind[1]));
        ind[2] = static_cast<long int>(itk::Math::rnd_halfinttoeven(cind[2]));
  	
        if(!labelimage1->GetLargestPossibleRegion().IsInside(ind))
        {
          std::cerr << "Error index: " << ind << " not in image"  << std::endl;
          std::cout << "Ignoring" << std::endl;
          //return EXIT_FAILURE;
        }
        else
        {
  	labelimageCountF->SetPixel(ind, labelimageCountF->GetPixel(ind) + 1);
  	//std::cout<<"value of pixel: "<<labelimageCountF->GetPixel(ind)<<std::endl;
          labelimage1->SetPixel(ind, voxelLabel);
        }
      }
    PointType fiberpoint2;
    double fiberpointtemp2[3];
    for(int j=0;j<TargetPoints->GetNumberOfPoints();j++)
    {	
        TargetPoints->GetPoint( j, fiberpointtemp2 );
  	 for (int k=0;k<3;k++)
  	{
  		fiberpoint2[k]=fiberpointtemp2[k];
  	}
        itk::ContinuousIndex<double,3> cind2;
        itk::Index<3> ind2;
        labelimage2->TransformPhysicalPointToContinuousIndex(fiberpoint2, cind2);
        ind2[0] = static_cast<long int>(itk::Math::rnd_halfinttoeven(cind2[0]));
        ind2[1] = static_cast<long int>(itk::Math::rnd_halfinttoeven(cind2[1]));
        ind2[2] = static_cast<long int>(itk::Math::rnd_halfinttoeven(cind2[2]));
  	
        if(!labelimage2->GetLargestPossibleRegion().IsInside(ind2))
        {
          std::cerr << "Error index: " << ind2 << " not in image"  << std::endl;
          std::cout << "Ignoring" << std::endl;
          //return EXIT_FAILURE;
        }
        else
        {
  	labelimage2CountF->SetPixel(ind2, labelimage2CountF->GetPixel(ind2) + 1);
          labelimage2->SetPixel(ind2, voxelLabel);
        }
    }
  
  ////overlap calculation////
  //With classic voxelisation
   typedef itk::LabelOverlapMeasuresImageFilter <IntImageType> LabelOverlapMeasuresImageFilterType;
   LabelOverlapMeasuresImageFilterType::Pointer  LabelOverlapMeasuresImageFilter= LabelOverlapMeasuresImageFilterType::New();
   LabelOverlapMeasuresImageFilter->SetSourceImage(labelimage1);
   LabelOverlapMeasuresImageFilter->SetTargetImage(labelimage2);
   LabelOverlapMeasuresImageFilter->Update();
  
  double total=LabelOverlapMeasuresImageFilter->GetTotalOverlap();
  double jaccard=LabelOverlapMeasuresImageFilter->GetUnionOverlap();
  double dice=LabelOverlapMeasuresImageFilter->GetMeanOverlap();
  double Volume_sim=LabelOverlapMeasuresImageFilter->GetVolumeSimilarity();
  double False_negative=LabelOverlapMeasuresImageFilter->GetFalseNegativeError();
  double False_positive=LabelOverlapMeasuresImageFilter->GetFalsePositiveError();
  OverlapTable.push_back(total);
  OverlapTable.push_back(jaccard);
  OverlapTable.push_back(dice);
  OverlapTable.push_back(Volume_sim);
  OverlapTable.push_back(False_negative);
  OverlapTable.push_back(False_positive);
  std::cout<<"total: "<<total<<" jaccard: "<<jaccard<<"dice: "<<dice<<"Volume_sim: "<<Volume_sim<<"False_negative: "<<False_negative<<"False_positive: "<<False_positive<<std::endl;
  //divide each voxel of Source voxelized by the total number of points to get the probability
  
  IntImageType::Pointer img1 = IntImageType::New();
  img1 = labelimageCountF;
  itk::ImageRegionIterator<IntImageType> img_it1 (img1, img1->GetLargestPossibleRegion());
  std::cout<<" get proba of Source "<<std::endl;
  img_it1.GoToBegin();
  //divide each voxel of Target voxelized by the total number of points to get the probability
  IntImageType::Pointer img2 = IntImageType::New();
  img2 = labelimage2CountF;
  itk::ImageRegionIterator<IntImageType> img_it2 (img2, img2->GetLargestPossibleRegion());
  std::cout<<" get proba of Target "<<std::endl;
  img_it2.GoToBegin();
   double numerator=0;
   double denominator=0;
  
  while(!img_it1.IsAtEnd() && !img_it2.IsAtEnd())
      {
      // Get the value of the current voxel
       double val1 = img_it1.Get();
       double  val2= img_it2.Get();
       //Get the proba on each voxel
       double Pa=val1/numberPointsSource;
       double Pb=val2/numberPointsTarget;
       numerator+=abs(Pa-Pb);
       denominator+=(Pa+Pb-(Pa*Pb));
       ++img_it1;
       ++img_it2;
      }
  
  std::cout<<"points in source : "<<numberPointsSource<<std::endl;
  std::cout<<"points in target : "<<numberPointsTarget<<std::endl;
  
  std::cout<<"compute POV  "<<std::endl;
  //Compute POV calculation
  double POV=1-(numerator/denominator);
  std::cout<<"POV : "<< POV<<std::endl;
  	
  OverlapTable.push_back(POV);
  return OverlapTable;
  }
  
  
  
  
  std::vector<double> FillCurvatureTable(vtkSmartPointer<vtkPolyData> SourceInterpolated, vtkSmartPointer<vtkPolyData> TargetInterpolated, double stepInterpolate)
  {	std::cout<<"Curvature metric calculating..."<<std::endl;
  	std::vector<double> CurveTable;
  	vtkPoints* SourcePointsNew = SourceInterpolated->GetPoints();
  	vtkPoints* TargetPointsNew = TargetInterpolated->GetPoints();
  	typedef itk::Vector<double, 3> VectorType;
  	std::cout<<" number of points in source: "<<SourcePointsNew->GetNumberOfPoints()<<std::endl;
  	std::cout<<" number of points in target: "<<TargetPointsNew->GetNumberOfPoints()<<std::endl;
  
  	vtkIdType nbpts=0;
          vtkIdType const *pts=0;
  	double Tbis[3],Tbis2[3];
  	double T[3],T2[3];
  	double curvature1=0;
  	double curvature2=0;
  	double dsbis=0;
  	double dsbis2=0;
  	double ds=0;
  	double ds2=0;
  	int compt=1;
  	int compt2=1;
  	vtkSmartPointer<vtkCellArray> A = vtkSmartPointer<vtkCellArray>::New();
  	A=SourceInterpolated->GetLines();
  	A->InitTraversal();
  	std::cout<<" number of fibers in source: "<<SourceInterpolated->GetNumberOfLines()<<std::endl;
  	while(SourceInterpolated->GetLines()->GetNextCell(nbpts, pts))
  	{	std::cout<<" fibers in source number: "<<compt2<<std::endl;
  		vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
   		vtkIdType ind=0;double coords[3];
          	for(int p=0; p<nbpts; p++)
          	{
              		ind=pts[p];
  			SourcePointsNew->GetPoint( ind, coords);
  			points->InsertNextPoint(coords);
  		}
  		//for each point of hte current fiber with its new points from the spline interpolation	
  		int i=0;
      		for( int k=0; k<points->GetNumberOfPoints(); k++ )
  		{	if(k==(points->GetNumberOfPoints()-1))
  			{
  				i=k-1;
  			}
  			else
  			{
  				i=k;
  			}
  			double DistanceMin=100000;
  			double SourceP[3];
  			points->GetPoint(i,SourceP);
  			double TargetPoint[3];
  			double Distance=0;
  			////for each fiber of vtk file Target////
  			vtkIdType nbpts2=0;
                          vtkIdType const *pts2=0;
  			vtkSmartPointer<vtkCellArray> P = vtkSmartPointer<vtkCellArray>::New();
  			P=TargetInterpolated->GetLines();
  			P->InitTraversal();
    			while(TargetInterpolated->GetLines()->GetNextCell(nbpts2, pts2))
  			{	
  				vtkSmartPointer<vtkPoints> points2 = vtkSmartPointer<vtkPoints>::New();
   				vtkIdType ind2=0;double coords2[3];
          			for(int p=0; p<nbpts2; p++)
          			{
              				ind2=pts2[p];
  					TargetPointsNew->GetPoint( ind2, coords2);
  					points2->InsertNextPoint(coords2);
  				}
  				//for each point of the current fiber with its new points from the spline interpolation	
  				int j=0;
        				for( int l=0; l<points2->GetNumberOfPoints(); l++ )
  				{	
  					if(l==(points2->GetNumberOfPoints()-1))
  					{
  						j=l-1;
  					}
  					else
  					{
  						j=l;
  					}
  					points2->GetPoint(j,TargetPoint);
  					
  					Distance=vtkMath::Distance2BetweenPoints(SourceP,TargetPoint );
  					if(Distance<DistanceMin )
  					{	
  						//std::cout<<"two closest points"<<SourceP[0]<<" "<<SourceP[1]<<" "<<SourceP[2]<<"and "<<TargetPoint[0]<<" "<<TargetPoint[1]<<" "<<TargetPoint[2]<<std::endl;
  	  					DistanceMin = Distance;//std::cout<<" distance : "<<DistanceMin<<std::endl;
  	  					double matchingpoint[3];matchingpoint[0]=TargetPoint[0];matchingpoint[1]=TargetPoint[1];matchingpoint[2]=TargetPoint[2];
  						double tempo2[3],tempo22[3];
  						points2->GetPoint(j+1,tempo2);
  						points2->GetPoint(j+2,tempo22);
  						ds = sqrt(vtkMath::Distance2BetweenPoints(TargetPoint,tempo2));
  						ds2= sqrt(vtkMath::Distance2BetweenPoints(tempo2,tempo22));
  						T[0] = (tempo2[0]-TargetPoint[0])/ds;
  						T[1] = (tempo2[1]-TargetPoint[1])/ds;
  						T[2] = (tempo2[2]-TargetPoint[2])/ds;
  						T2[0] = (tempo22[0]-tempo2[0])/ds2;
  						T2[1] = (tempo22[1]-tempo2[1])/ds2;
  						T2[2] = (tempo22[2]-tempo2[2])/ds2;
  						curvature2=sqrt(vtkMath::Distance2BetweenPoints(T2,T))/ds;
  
  					}
        				}
  				
      			}
  			
  			double temp2[3],temp22[3];
  			points->GetPoint(i+1,temp2);
  			points->GetPoint(i+2,temp22);
  			dsbis = sqrt(vtkMath::Distance2BetweenPoints(SourceP,temp2));
  			dsbis2= sqrt(vtkMath::Distance2BetweenPoints(temp2,temp22));
  			Tbis[0] = (temp2[0]-SourceP[0])/dsbis;
  			Tbis[1] = (temp2[1]-SourceP[1])/dsbis;
  			Tbis[2] = (temp2[2]-SourceP[2])/dsbis;
  			Tbis2[0] = (temp22[0]-temp2[0])/dsbis2;
  			Tbis2[1] = (temp22[1]-temp2[1])/dsbis2;
  			Tbis2[2] = (temp22[2]-temp2[2])/dsbis2;
  			curvature1=sqrt(vtkMath::Distance2BetweenPoints(Tbis2,Tbis))/dsbis;
  
  			double val =(curvature1-curvature2)*(curvature1-curvature2);
  			//std::cout<<"   val : "<<val<<" number : "<<compt<<std::endl;
  			compt++;
  			CurveTable.push_back(val);
  			
  		}
  	compt2++;	
  
  	}
  	std::cout<<"Done..."<<std::endl;
  	return CurveTable;
  
  	
  }
  
  
  std::vector<double> FillTangentTable(vtkSmartPointer<vtkPolyData> SourceInterpolated, vtkSmartPointer<vtkPolyData> TargetInterpolated, double stepInterpolate)
  {	std::cout<<"Tangent metric calculating..."<<std::endl;
  	std::vector<double> TangentTable;
  	std::vector<double> samples1, samples2;
  	//for each point of the current fiber with its new points from the spline interpolation	
  	vtkPoints* SourcePointsNew = SourceInterpolated->GetPoints();
  	vtkPoints* TargetPointsNew = TargetInterpolated->GetPoints();
  	
  	vtkIdType nbpts=0;
          vtkIdType const *pts=0;
  	typedef itk::Vector<double, 3> VectorType;
  	VectorType v1, v2;		
  	std::cout<<" number of points in source: "<<SourcePointsNew->GetNumberOfPoints()<<std::endl;
  	std::cout<<" number of points in target: "<<TargetPointsNew->GetNumberOfPoints()<<std::endl;
  
  	int compt=1;
  	std::cout<<" number of fibers in source: "<<SourceInterpolated->GetNumberOfLines()<<std::endl;
  	vtkSmartPointer<vtkCellArray> A = vtkSmartPointer<vtkCellArray>::New();
  	A=SourceInterpolated->GetLines();
  	A->InitTraversal();
  	while(SourceInterpolated->GetLines()->GetNextCell(nbpts, pts))
  	{	//std::cout<<" fibers in source number: "<<compt<<std::endl;
  		vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
   		vtkIdType ind=0;double coords[3];
          	for(int p=0; p<nbpts; p++)
          	{
              		ind=pts[p];
  			SourcePointsNew->GetPoint( ind, coords);
  			points->InsertNextPoint(coords);
  		}
  		
  		//std::cout<<"number of points: "<<points->GetNumberOfPoints()<<std::endl;
      		for( int i=0; i<points->GetNumberOfPoints(); i++ )
  		{
  			double DistanceMin=100000;
  			double SourceP[3];
  			points->GetPoint(i,SourceP);	
  			double TargetPoint[3];
  			double Distance=0;
  			vtkIdType nbpts2=0;
                          vtkIdType const *pts2=0;
  			vtkSmartPointer<vtkCellArray> P = vtkSmartPointer<vtkCellArray>::New();
  			P=TargetInterpolated->GetLines();
  			P->InitTraversal();
  			while(TargetInterpolated->GetLines()->GetNextCell(nbpts2, pts2))
  			{
  				
  				vtkSmartPointer<vtkPoints> points2 = vtkSmartPointer<vtkPoints>::New();
   				vtkIdType ind2=0;double coords2[3];
  				
          			for(int p=0; p<nbpts2; p++)
          			{
              				ind2=pts2[p];
  					TargetPointsNew->GetPoint( ind2, coords2);
  					points2->InsertNextPoint(coords2);
  				}
  			//for each point of the current fiber with its new points from the spline interpolation	
        				for( int j=0; j<points2->GetNumberOfPoints(); j++ )
  				{	
  					points2->GetPoint(j,TargetPoint);
  					
  					Distance=vtkMath::Distance2BetweenPoints(SourceP,TargetPoint );
  					if(Distance<DistanceMin )
  					{
  	  					DistanceMin = Distance;
  						//std::cout<<"distance min :"<<DistanceMin<<std::endl;
  	  					double matchingpoint[3];matchingpoint[0]=TargetPoint[0];matchingpoint[1]=TargetPoint[1];matchingpoint[2]=TargetPoint[2];
  						double tempo1[3],tempo2[3];
  						
  						if (j==0)////if we are at the beginning of the fiber////
  						{
  					 		points2->GetPoint(j+1,tempo2);
  							for(int k=0;k<3;k++)
  							{
  								v2[k]=tempo2[k]-TargetPoint[k];
  							}
  						}
  						if(j==points2->GetNumberOfPoints()-1)////if we are at the end of the fiber////
  						{	
  							points2->GetPoint(j-1,tempo1);
  							for(int k=0;k<3;k++)
  							{
  								v2[k]=TargetPoint[k]-tempo1[k];
  							}
  
  						}
  						if(j!=0 && j!=points2->GetNumberOfPoints()-1)
  						{	
  							points2->GetPoint(j-1,tempo1);
  							points2->GetPoint(j+1,tempo2);
  							for(int k=0;k<3;k++)
  							{
  								v2[k]=tempo2[k]-tempo1[k];
  							}
  
  						}
  
  					}
        				}
  		
  			}	
      	
  		double temp1[3],temp2[3];
  		//std::cout<<"i : "<<i<<std::endl;
  		if (i==0)////if we are at the beginning of the fiber////
  		{
  			points->GetPoint(i+1,temp2);
  			for(int k=0;k<3;k++)
  			{
  				v1[k]=temp2[k]-SourceP[k];
  			}
  		}
  		if (i==points->GetNumberOfPoints()-1)////if we are at the end of the fiber////
  		{
  			points->GetPoint(i-1,temp1);
  			for(int k=0;k<3;k++)
  			{
  				v1[k]=SourceP[k]-temp1[k];
  			}
  				
  		}
  		if(i!=0 && i!=points->GetNumberOfPoints()-1)
  		{
  			points->GetPoint(i-1,temp1);
  			points->GetPoint(i+1,temp2);
  			for(int k=0;k<3;k++)
  			{
  				v1[k]=temp2[k]-temp1[k];
  			}
  				
  		}
  
      		v1.Normalize();
      		v2.Normalize();
  		//std::cout<<v1[0]<<" "<<v1[1]<<" "<<v1[2]<<std::endl;
  		//std::cout<<v2[0]<<" "<<v2[1]<<" "<<v2[2]<<std::endl;
  
  		double val = pow ( acos ( fabs ( v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]) )*180/3.14159265 , 2.0);
  		/*if(compt==211)
  			std::cout<<"val : "<<(acos ( fabs ( v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]) )*180/3.14159265)<<std::endl;*/
  		TangentTable.push_back(val);
  			
  		}
  	
  	compt++;
  		
  	}
  	std::cout<<"done..."<<std::endl;	
  	return TangentTable;
  }
  
  
  std::vector<double> FillDistanceTable(vtkSmartPointer<vtkPolyData> Source, vtkSmartPointer<vtkPolyData> Target)
  {
  	std::cout<<"Distance (closest points) metric calculating..."<<std::endl;
  	std::vector<double> DistanceTable;
  	vtkPoints* SourcePoints = Source->GetPoints();
  	vtkPoints* TargetPoints = Target->GetPoints();
  	
  	for(int i=0; i<SourcePoints->GetNumberOfPoints(); i++)
  	{
  		double SourcePoint[3], DistanceMin = 10000;
  		SourcePoints->GetPoint(i,SourcePoint);
  		for(int j=0; j<TargetPoints->GetNumberOfPoints(); j++)
  		{
  			double TargetPoint[3], XYZ[3], Distance;
  			TargetPoints->GetPoint(j,TargetPoint);
  			for(int k=0; k<3; k++)
  				XYZ[k] = TargetPoint[k]-SourcePoint[k];
  			Distance = sqrt(XYZ[0]*XYZ[0]+XYZ[1]*XYZ[1]+XYZ[2]*XYZ[2]);
  			if(Distance<DistanceMin)
  				DistanceMin = Distance;
  		}
  		DistanceTable.push_back(DistanceMin);
  	}
  	std::cout<<"Done..."<<std::endl;
  	return DistanceTable;
  }
  
  void GetBounds(std::vector<std::vector<double> > DistanceTable, double Bounds[2])
  {
  	double Min = 100000, Max = -1;
  	for(unsigned int i=0; i<DistanceTable.size(); i++)
  	{
  		for(unsigned int j=0; j<DistanceTable[i].size(); j++)
  		{
  			if(DistanceTable[i][j]<Min)
  				Min = DistanceTable[i][j];
  			if(DistanceTable[i][j]>Max)
  				Max = DistanceTable[i][j];
  		}
  	}
  	Bounds[0] = Min;
  	Bounds[1] = Max;
  }
  
  int GetFrequency(std::vector<std::vector<double> > DistanceTable, double IntervalMin, double IntervalMax)
  {
  	int Frequency=0;
  	for(unsigned int i=0; i<DistanceTable.size(); i++)
  	{
  		for(unsigned int j=0; j<DistanceTable[i].size(); j++)
  		{
  			if(DistanceTable[i][j]<IntervalMax && DistanceTable[i][j]>=IntervalMin)
  				Frequency++;
  		}
  	}
  	return Frequency;
  }
  
  double GetMeanDistance(std::vector<double> Distance, std::vector<double> Frequency)
  {
  	double MeanDistance=0, TotalFrequency=0;
  	for(unsigned int i=0; i<Distance.size(); i++)
  	{
  		MeanDistance+=Distance[i]*Frequency[i];
  		TotalFrequency+=Frequency[i];
  	}
  	MeanDistance/=TotalFrequency;
  	
  	return MeanDistance;
  }
  
  
  std::vector<std::vector<double> > GetResultTableFrequency(std::vector<std::vector<double> >  Table,int TotalNumberOfFibers,int number_of_entries,double step)
  {
  	double Bounds[2];
  	GetBounds(Table,Bounds);
  	if(step==-1)
  		step=(Bounds[1]-Bounds[0])/number_of_entries;
  	std::cout<<Bounds[0]<<" "<<Bounds[1]<<" "<<number_of_entries<<" "<<step<<std::endl;
  	double IntervalMin = Bounds[0], IntervalMax = Bounds[0] + step;
  	std::vector<double> FrequencyVector;
  	std::vector<double> DistanceVector;
  		
  	while(IntervalMin<=Bounds[1])
  	{
  		int Frequency = GetFrequency(Table, IntervalMin, IntervalMax);
  		DistanceVector.push_back(IntervalMin);
  		FrequencyVector.push_back(Frequency);
  		IntervalMin += step;
  		IntervalMax += step;
  	}
  		
  	double Dist25=-1, Dist50=-1, Dist75=-1, Dist90=-1, Dist95=-1;
  	double CumulatedFrequency=0;
  
  	for(unsigned int i=0; i<DistanceVector.size(); i++)
  	{
  		//StatFile<<DistanceVector[i]<<","<<FrequencyVector[i]<<std::endl;
  		CumulatedFrequency+=FrequencyVector[i];
  		//std::cout<<"cumulated frequency"<<CumulatedFrequency<<"TotalNumberOfFibers"<<TotalNumberOfFibers<<std::endl;
  		if(CumulatedFrequency>=TotalNumberOfFibers*0.25 && Dist25==-1)
  		{
  			Dist25=DistanceVector[i];
  		}
  		if(CumulatedFrequency>=TotalNumberOfFibers*0.5 && Dist50==-1)
  		{
  			Dist50=DistanceVector[i];
  		}
  		if(CumulatedFrequency>=TotalNumberOfFibers*0.75 && Dist75==-1)
  		{
  			Dist75=DistanceVector[i];
  		}
  		if(CumulatedFrequency>=TotalNumberOfFibers*0.9 && Dist90==-1)
  		{
  			Dist90=DistanceVector[i];
  		}
  		if(CumulatedFrequency>=TotalNumberOfFibers*0.95 && Dist95==-1)
  		{
  			Dist95=DistanceVector[i];
  		}
  	}
  	std::vector<std::vector<double> > ResultDistanceFrequency;
  	std::vector<double> DistanceFrequency;
  	DistanceFrequency.push_back(Dist25);
  	DistanceFrequency.push_back(Dist50);
  	DistanceFrequency.push_back(Dist75);
  	DistanceFrequency.push_back(Dist90);
  	DistanceFrequency.push_back(Dist95);
  	DistanceFrequency.push_back(Bounds[0]);
  	DistanceFrequency.push_back(Bounds[1]);
  	ResultDistanceFrequency.push_back(DistanceFrequency);
  	ResultDistanceFrequency.push_back(DistanceVector);
  	ResultDistanceFrequency.push_back(FrequencyVector);
  	
  
   return ResultDistanceFrequency;
  
  }
  
  
  
  
  int main(int argc, char* argv[])
  {
  	PARSE_ARGS;
  
  	int Number_of_entries=number_of_entries;
  	double Step=step;
  	std::vector<vtkSmartPointer<vtkPolyData> > FiberTracts;
  	std::vector<std::string> Filenames;
  
  	std::cout<<"Reading VTK data..."<<std::endl;
  	
  	Filenames.push_back(vtk_input1);
  	Filenames.push_back(vtk_input2);
  	for(int i=0; i<2; i++)
  	{
  		if(Filenames[i].rfind(".vtk") != std::string::npos)
  		{
  			vtkSmartPointer<vtkPolyDataReader> reader = vtkSmartPointer<vtkPolyDataReader>::New();
  			reader->SetFileName(Filenames[i].c_str());
  			FiberTracts.push_back(reader->GetOutput());
  			reader->Update();
  		}
  		else if(Filenames[i].rfind(".vtp") != std::string::npos)
  		{
  			vtkSmartPointer<vtkXMLPolyDataReader> reader = vtkSmartPointer<vtkXMLPolyDataReader>::New();
  			reader->SetFileName(Filenames[i].c_str());
  			FiberTracts.push_back(reader->GetOutput());
  			reader->Update();
  		}
  		else
  		{
  			std::cout<<"Wrong filename : "<<Filenames[i]<<" Check file format and location."<<std::endl;
  			return 0;
  		}
  	}
  	std::cout<<Filenames[0]<<" "<<FiberTracts[0]->GetNumberOfCells()<<" fibers, "<<FiberTracts[0]->GetPoints()->GetNumberOfPoints()<<" points."<<std::endl;
  	std::cout<<Filenames[1]<<" "<<FiberTracts[1]->GetNumberOfCells()<<" fibers, "<<FiberTracts[1]->GetPoints()->GetNumberOfPoints()<<" points."<<std::endl;
  	std::cout<<"VTK Files read successfuly."<<std::endl<<std::endl;
  
  	std::cout<<"interpolation"<<std::endl<<std::endl;
  	vtkSmartPointer<vtkPolyData> Source = vtkSmartPointer<vtkPolyData>::New();
  	Source=FiberTracts[0];
  	vtkSmartPointer<vtkPolyData> Target= vtkSmartPointer<vtkPolyData>::New();
  	Target=FiberTracts[1];
  	vtkPoints* SourcePoints = Source->GetPoints();
  	vtkPoints* TargetPoints = Target->GetPoints();
  	std::cout<<"points "<<SourcePoints->GetNumberOfPoints()<<std::endl;
  	std::cout<<"points "<<TargetPoints->GetNumberOfPoints()<<std::endl;
  	//new vtk with spline interpolation
  	vtkSmartPointer<vtkPolyData> SourceInterpolated = vtkSmartPointer<vtkPolyData>::New();
  	vtkSmartPointer<vtkPolyData> TargetInterpolated = vtkSmartPointer<vtkPolyData>::New();
  	typedef itk::Vector<double, 3> VectorType;
  	
  	vtkIdType nbpts=0;
          vtkIdType const *pts=0;
  	VectorType v1, v2;
  	////for each fiber of vtk file Source////
  	vtkSmartPointer <vtkCellArray> cells = vtkSmartPointer <vtkCellArray>::New();
  	int countID=0;
  	vtkPoints* polypoints = vtkPoints::New();
    	while(Source->GetLines()->GetNextCell(nbpts, pts))
  	{	
  		//store points of the current fiber in points
  		vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
   		vtkIdType ind=0;double coords[3];
          	for(int p=0; p<nbpts; p++)
          	{
              		ind=pts[p];
  			SourcePoints->GetPoint( ind, coords);
  			points->InsertNextPoint(coords);
  		}
  		//spline interpolation
  		int numberOfInputPoints = points->GetNumberOfPoints();
     		
      		vtkCardinalSpline* aSplineX;
      		vtkCardinalSpline* aSplineY;
     		vtkCardinalSpline* aSplineZ;
     		aSplineX = vtkCardinalSpline::New();
     		aSplineY = vtkCardinalSpline::New();
      		aSplineZ = vtkCardinalSpline::New();
  		aSplineX->ClosedOff();
  		aSplineY->ClosedOff();
  		aSplineZ->ClosedOff();
  
  		/////creation of splines/////
     		 for (int i=0; i<numberOfInputPoints; i++)
     		 {
  		
      		  double x = points->GetPoint(i)[0];
      		  double y = points->GetPoint(i)[1];
       		  double z = points->GetPoint(i)[2];
            	  aSplineX->AddPoint(i, x);
            	  aSplineY->AddPoint(i, y);
            	  aSplineZ->AddPoint(i, z);
      		 }
  	
  
  		double  tRange[2];
  		aSplineX->GetParametricRange 	( tRange);
  		int compteur2=0;
      		
  		////evaluation of new points then store it in polypoints////
  		for( double t = 0 ; t <= tRange[1]; t += stepInterpolate )
      		{
        		polypoints->InsertNextPoint(aSplineX->Evaluate(t),aSplineY->Evaluate(t),aSplineZ->Evaluate(t));
  		//compteur2 gives us the number of new points of the current fiber
  		compteur2++;
      		}
  		//creation of polylines
  		vtkSmartPointer<vtkPolyLine> polyline = vtkSmartPointer<vtkPolyLine>::New();
  		polyline->GetPointIds()->SetNumberOfIds(compteur2);
  		//countID gives us the total number of points cumulated
  		countID=polypoints->GetNumberOfPoints();
  		
  		for (int i= 0; i<compteur2 ;i++)
  		{	//SETID(0->number of points in the current fiber, (total points cumulated-points in the current fiber)->total points cumulated)
  			polyline->GetPointIds()->SetId(i,i+(countID-compteur2));
  		}
  		cells->InsertNextCell(polyline);
  		aSplineX->Delete();
  		aSplineY->Delete();
  		aSplineZ->Delete();
  		
  	}
  	SourceInterpolated->SetPoints(polypoints);
  	SourceInterpolated->SetLines(cells);
  	polypoints->Delete();
  	std::cout<<"end of the first interpolation"<<std::endl<<std::endl;
  
  	vtkIdType nbpts2=0;
          vtkIdType const *pts2=0;
  	vtkSmartPointer <vtkCellArray> cells2 = vtkSmartPointer <vtkCellArray>::New();
  	int countID2=0;
  	vtkPoints* polypoints2 = vtkPoints::New();
  	while(Target->GetLines()->GetNextCell(nbpts2, pts2))
  	{	
  		//std::cout<<"step2"<<std::endl;
  		vtkSmartPointer<vtkPoints> points2 = vtkSmartPointer<vtkPoints>::New();
   		vtkIdType ind2=0;double coords2[3];
          	for(int p2=0; p2<nbpts2; p2++)
          	{
              		ind2=pts2[p2];
  			TargetPoints->GetPoint( ind2, coords2);
  			points2->InsertNextPoint(coords2);
  		}
  		//spline interpolation
  		int numberOfInputPoints2 = points2->GetNumberOfPoints();
      		vtkCardinalSpline* aSpline2X;
      		vtkCardinalSpline* aSpline2Y;
     		vtkCardinalSpline* aSpline2Z;
  
     		aSpline2X = vtkCardinalSpline::New();
     		aSpline2Y = vtkCardinalSpline::New();
      		aSpline2Z = vtkCardinalSpline::New();
  		aSpline2X->ClosedOff();
  		aSpline2Y->ClosedOff();
  		aSpline2Z->ClosedOff();
  		/////creation of splines/////
     		 for (int i=0; i<numberOfInputPoints2; i++)
     		{
  			//std::cout<<"step6"<<std::endl;
      		 	double x2 = points2->GetPoint(i)[0];
      		 	double y2 = points2->GetPoint(i)[1];
       		 	double z2 = points2->GetPoint(i)[2];
            	 	aSpline2X->AddPoint(i, x2);
            	  	aSpline2Y->AddPoint(i, y2);
            	  	aSpline2Z->AddPoint(i, z2);
      		}
  		double  tRange2[2];
  		aSpline2X->GetParametricRange 	( tRange2);
      		int compteur=0;
  		////evaluation of new points then store it in polypoints2////
  		for( double t2 = 0 ; t2 <= tRange2[1]; t2 += stepInterpolate )
      		{	
        		polypoints2->InsertNextPoint(aSpline2X->Evaluate(t2),aSpline2Y->Evaluate(t2),aSpline2Z->Evaluate(t2));
  		compteur++;
      		}
  		vtkSmartPointer<vtkPolyLine> polyline2 = vtkSmartPointer<vtkPolyLine>::New();
  		polyline2->GetPointIds()->SetNumberOfIds(compteur);
  		countID2=polypoints2->GetNumberOfPoints();
  		for (int i=0; i<compteur;i++)
  		{
  			polyline2->GetPointIds()->SetId(i,i+(countID2-compteur));
  		}
  		cells2->InsertNextCell(polyline2);
  		aSpline2X->Delete();
  		aSpline2Y->Delete();
  		aSpline2Z->Delete();
  		
  		//end
  
  	}
  	TargetInterpolated->SetPoints(polypoints2);
  	TargetInterpolated->SetLines(cells2);
  	polypoints2->Delete();
  	std::cout<<"end of the second interpolation"<<std::endl<<std::endl;
  
  	std::cout<<"Calculating..."<<std::endl;
  			
  	std::vector<std::vector<double> > DistanceTable;
  	DistanceTable.push_back(FillDistanceTable(SourceInterpolated,TargetInterpolated));
  	DistanceTable.push_back(FillDistanceTable(TargetInterpolated,SourceInterpolated));
  
  	std::vector<std::vector<double> > TangentTable;
  	TangentTable.push_back(FillTangentTable(SourceInterpolated,TargetInterpolated,stepInterpolate));
  	TangentTable.push_back(FillTangentTable(TargetInterpolated,SourceInterpolated,stepInterpolate));
  		
  	std::vector<std::vector<double> > CurveTable;
  	CurveTable.push_back(FillCurvatureTable(SourceInterpolated,TargetInterpolated,stepInterpolate));
  	CurveTable.push_back(FillCurvatureTable(TargetInterpolated,SourceInterpolated,stepInterpolate));
  	
  	int TotalNumberOfFibers=SourceInterpolated->GetPoints()->GetNumberOfPoints()+TargetInterpolated->GetPoints()->GetNumberOfPoints();
  	
  		
  	std::ofstream StatFile(output_stat_file.c_str(), std::ios::out);
  	std::cout<<"number total points"<<TotalNumberOfFibers<<std::endl;	
  	if(StatFile)
  	{
  		StatFile<<"Filename,Number of fibers,Number of points"<<std::endl;
  		StatFile<<Filenames[0]<<","<<SourceInterpolated->GetNumberOfCells()<<","<<SourceInterpolated->GetPoints()->GetNumberOfPoints()<<std::endl;
  		StatFile<<Filenames[1]<<","<<TargetInterpolated->GetNumberOfCells()<<","<<TargetInterpolated->GetPoints()->GetNumberOfPoints()<<std::endl<<std::endl;
  		std::vector<std::vector<double> > DistanceTableFrequency=GetResultTableFrequency(DistanceTable,TotalNumberOfFibers,Number_of_entries,Step);
  		std::cout<<"frequency tangent metric"<<std::endl;
  		std::vector<std::vector<double> > TangentTableFrequency=GetResultTableFrequency(TangentTable,TotalNumberOfFibers,Number_of_entries,Step);
  		std::cout<<"frequency curve metric"<<std::endl;
  		std::vector<std::vector<double> > CurveTableFrequency=GetResultTableFrequency(CurveTable,TotalNumberOfFibers,Number_of_entries,Step);
  
  		StatFile<<"Distance,Frequency closest points"<<std::endl;
  		std::vector<double> FrequencyVector1=DistanceTableFrequency[2];
  		std::vector<double> DistanceVector1=DistanceTableFrequency[1];
  		for(unsigned int i=0; i<DistanceVector1.size(); i++)
  		{
  			StatFile<<DistanceVector1[i]<<","<<FrequencyVector1[i]<<std::endl;
  		}
  		StatFile<<std::endl;
  
  		StatFile<<"Distance,Frequency of tangent metric"<<std::endl;
  		std::vector<double> FrequencyVector2=TangentTableFrequency[2];
  		std::vector<double> DistanceVector2=TangentTableFrequency[1];
  		for(unsigned int i=0; i<DistanceVector2.size(); i++)
  		{
  			StatFile<<DistanceVector2[i]<<","<<FrequencyVector2[i]<<std::endl;
  		}
  		StatFile<<std::endl;
  
  		StatFile<<"Distance,Frequency of curve metric"<<std::endl;
  		std::vector<double> FrequencyVector3=CurveTableFrequency[2];
  		std::vector<double> DistanceVector3=CurveTableFrequency[1];
  		for(unsigned int i=0; i<DistanceVector3.size(); i++)
  		{
  			StatFile<<DistanceVector3[i]<<","<<FrequencyVector3[i]<<std::endl;
  		}
  		StatFile<<std::endl;
  
  		StatFile<<"Distance Frequency Percentage,Distance"<<std::endl;
  		StatFile<<"25%,"<<DistanceTableFrequency[0][0]<<std::endl;
  		StatFile<<"50%,"<<DistanceTableFrequency[0][1]<<std::endl;
  		StatFile<<"75%,"<<DistanceTableFrequency[0][2]<<std::endl;
  		StatFile<<"90%,"<<DistanceTableFrequency[0][3]<<std::endl;
  		StatFile<<"95%,"<<DistanceTableFrequency[0][4]<<std::endl;
  		StatFile<<std::endl;
  		StatFile<<"Tangent metric Frequency Percentage,Distance"<<std::endl;
  		StatFile<<"25%,"<<TangentTableFrequency[0][0]<<std::endl;
  		StatFile<<"50%,"<<TangentTableFrequency[0][1]<<std::endl;
  		StatFile<<"75%,"<<TangentTableFrequency[0][2]<<std::endl;
  		StatFile<<"90%,"<<TangentTableFrequency[0][3]<<std::endl;
  		StatFile<<"95%,"<<TangentTableFrequency[0][4]<<std::endl;
  		StatFile<<std::endl;
  		StatFile<<"Curvature metric Frequency Percentage,Distance"<<std::endl;
  		StatFile<<"25%,"<<CurveTableFrequency[0][0]<<std::endl;
  		StatFile<<"50%,"<<CurveTableFrequency[0][1]<<std::endl;
  		StatFile<<"75%,"<<CurveTableFrequency[0][2]<<std::endl;
  		StatFile<<"90%,"<<CurveTableFrequency[0][3]<<std::endl;
  		StatFile<<"95%,"<<CurveTableFrequency[0][4]<<std::endl;
  		StatFile<<std::endl;
  		
  
  		for(unsigned int i=0; i<methods.size(); i++)
  		{	
  			if(methods[i] == "Overlap")
  			{	if( ReferenceScalarVolume == "")
      				{
        					std::cerr << "A reference volume has to be specified" << std::endl;
       					 return EXIT_FAILURE;
     				 }
  				std::vector<double> OverlapTable;
  				OverlapTable=FillOverlapTable(FiberTracts[0],FiberTracts[1],ReferenceScalarVolume,voxelLabel);
  				StatFile<<"Total volumetric overlap, "<<OverlapTable[0]<<std::endl;
  				StatFile<<"Union volumetric overlap (Jaccard coefficient), "<<OverlapTable[1]<<std::endl;
  				StatFile<<"Mean volumetric overlap (Dice coefficient), "<<OverlapTable[2]<<std::endl;
  				StatFile<<"Volume similarity, "<<OverlapTable[3]<<std::endl;
  				StatFile<<"False negative error, "<<OverlapTable[4]<<std::endl;
  				StatFile<<"False positive error, "<<OverlapTable[5]<<std::endl;
  				StatFile<<"Probabilistic overlap (POV), "<<OverlapTable[6]<<std::endl;
  			}
  
  			if(methods[i] == "Hausdorff")
  			{
  				std::cout<<"	Hausdorff..."<<std::endl;
  				
  				StatFile<<"Distance: 100%,"<<DistanceTableFrequency[0][6]<<std::endl;
  				StatFile<<"Tangent metric: 100%,"<<TangentTableFrequency[0][6]<<std::endl;
  				StatFile<<"Curve metric: 100%,"<<CurveTableFrequency[0][6]<<std::endl;
  				
  				std::cout<<"	End of Hausdorff."<<std::endl;
  			}
  			if(methods[i] == "Mean")
  			{
  				std::cout<<"	Mean..."<<std::endl;
  				
  				StatFile<<"Distance Mean,"<<GetMeanDistance(DistanceVector1,FrequencyVector1)<<std::endl;
  				StatFile<<"Tangent metric  Mean,"<<GetMeanDistance(DistanceVector2,FrequencyVector2)<<std::endl;
  				StatFile<<"Curve metric Mean,"<<GetMeanDistance(DistanceVector3,FrequencyVector3)<<std::endl;
  				
  				std::cout<<"	End of Mean."<<std::endl;
  			}
  			else if(methods[i] == "None")
  			{
  				std::cout<<"Ignoring Mean and Hausdorff methods"<<std::endl;
  			}
  			/*else
  			{
  				std::cout<<"Wrong method or syntax for argument, none option is selected : "<<methods[i]<<std::endl;
  				std::cout<<"Ignoring argument."<<std::endl;
  			}*/
  		}
  	StatFile.close();
  	}
  	else
  		std::cout<<"ERROR: Unable to save output stat file."<<std::endl;
  	
  	std::cout<<"Calculation complete."<<std::endl;
  	
  	return 0;
  }
  
  