#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif

#include <fstream>
#include "itkVector.h" 
#include "itkPoint.h"
#include <itkPointSet.h>
#include "itkImage.h"
#include "itkImageFileWriter.h" 
#include "itkVectorLinearInterpolateImageFunction.h"
#include <itkImageRegionIterator.h> 

int main( int argc, char * argv[] )
{
  if( argc < 2 )
    {
    std::cerr << "Missing Parameters " << std::endl;
    std::cerr << "Usage: " << argv[0];
    std::cerr << " landmarksFile";
    std::cerr << "deformationField" << std::endl;
    return EXIT_FAILURE;
    }
	
	//typedef   itk::Vector< float, ImageDimension >  FieldVectorType;
	typedef   itk::Vector< float, 2 >  PixelType;
    typedef   itk::PointSet< PixelType, 2 >  PointSetType;
	//typedef   itk::Point< double, 2 >  PointType;
	typedef   TransformType::PointType      PointType;
	typedef   PointSetType::PointIdentifier  PointIdType;
    typedef   itk::Image< PixelType,  2 >   DeformationFieldType;
	typedef   itk::VectorLinearInterpolateImageFunction< DeformationFieldType, float >     TransformType;
	typedef   itk::ImageFileWriter< DeformationFieldType >  FieldWriterType;
	typedef   TransformType::OutputType OutputType;

	DeformationFieldType::SizeType size = { { 98,126 } };
    double origin [2] = { 0.5,   0.5 };
    double spacing[2] = { 1,   1};


	PointSetType::Pointer PointSet=PointSetType::New();
	
	//PointSetType::PointType   point;
    PointSetType::PixelType displacement;
	PointType p1; 
	PointType p2;

	PointSetType::Pointer sourceLandMarks = PointSetType::New();
    PointSetType::Pointer targetLandMarks = PointSetType::New();

    PointSetType::PointsContainer::Pointer sourceLandMarkContainer = 
                                   sourceLandMarks->GetPoints();
    PointSetType::PointsContainer::Pointer targetLandMarkContainer = 
                                   targetLandMarks->GetPoints();
    // Software Guide : EndCodeSnippet

    PointIdType id = itk::NumericTraits< PointIdType >::Zero;

    // Read in the list of landmarks
    std::ifstream infile;
    infile.open( argv[1] );
    while (!infile.eof())
      {
      infile >>  p1[0] >> p1[1]  >> p2[0] >> p2[1]; 

      // Software Guide : BeginCodeSnippet
      sourceLandMarkContainer->InsertElement( id, p1 );
      targetLandMarkContainer->InsertElement( id++, p2 );
      // Software Guide : EndCodeSnippet

      //for(int i=0;i<2;i++)
	  //{
		//  displacement[i]=p2[i]-p1[i];
	  //}
      } 
    infile.close();

	OutputType  point1;
    OutputType  point2;
	
	//PixelType displacement2;
	TransformType::Pointer linear = TransformType::New();
    point1=linear->Evaluate(p1);
	point2=linear->Evaluate(p2);
    for (int i=0;i<2;i++) displacement[i] = point2[i] - point1[i];
	
	DeformationFieldType::RegionType region;
    region.SetSize( size );

    DeformationFieldType::Pointer field = DeformationFieldType::New();
    field->SetRegions( region );
    field->SetOrigin( origin );
    field->SetSpacing( spacing );
    field->Allocate();

    typedef itk::ImageRegionIterator< DeformationFieldType > FieldIterator;
    FieldIterator fi( field, region );
    fi.GoToBegin();

    DeformationFieldType::IndexType index;

    

    while( ! fi.IsAtEnd() )
      {
      index = fi.GetIndex();
      fi.Set( displacement );
      ++fi;
      }
  
	  
	FieldWriterType::Pointer fieldWriter = FieldWriterType::New();
    fieldWriter->SetFileName( argv[2] );
    fieldWriter->SetInput( field );
    try
      {
      fieldWriter->Update();
      }
    catch( itk::ExceptionObject & excp )
      {
      std::cerr << "Exception thrown " << std::endl;
      std::cerr << excp << std::endl;
      return EXIT_FAILURE;
      }
}