#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"

#include "itkResampleImageFilter.h"
#include "itkAffineTransform.h"
#include "itkNearestNeighborInterpolateImageFunction.h"

using namespace std;

int main( void)
{
  const     unsigned int   Dimension = 3;
  typedef   unsigned char  InputPixelType;
  typedef   unsigned char  OutputPixelType;
  typedef itk::Image< InputPixelType,  Dimension >   InputImageType;
  typedef itk::Image< OutputPixelType, Dimension >   OutputImageType;
  typedef itk::ImageFileReader< InputImageType  >  ReaderType;
  typedef itk::ImageFileWriter< OutputImageType >  WriterType;

  ReaderType::Pointer reader = ReaderType::New();
  WriterType::Pointer writer = WriterType::New();
 
  // a working volume
  reader->SetFileName ("brainweb1e1a10f20.mha"  );
  writer->SetFileName ("ScaleTransform_brainweb1e1a10f20.hdr");

  //// not working volume (a segmented volume)
  //reader->SetFileName("ConfiConnect_brainweb1e1a10f20_grayMatter.hdr");
  //writer->SetFileName("ScaleTransform_ConfiConnect_brainweb1e1a10f20_grayMatter.hdr");

  typedef itk::ResampleImageFilter<InputImageType,OutputImageType> FilterType;
  FilterType::Pointer filter = FilterType::New();
  
  typedef itk::AffineTransform< double, Dimension >  TransformType;  
  TransformType::Pointer transform = TransformType::New();
  filter->SetTransform( transform );
 
  typedef itk::NearestNeighborInterpolateImageFunction< InputImageType, double >  InterpolatorType; 
  InterpolatorType::Pointer interpolator = InterpolatorType::New();
  filter->SetInterpolator( interpolator );
  
  //filter->SetDefaultPixelValue( 0 );
 
  //----- set the spacing of the filter
  double spacing[ Dimension ];
  reader->Update();
  const InputImageType::SpacingType& sp = reader->GetOutput()->GetSpacing();
  spacing[0] = sp[0]; // pixel spacing in millimeters along X
  spacing[1] = sp[1]; // pixel spacing in millimeters along Y
  spacing[2] = sp[2]; // pixel spacing in millimeters along Y
  filter->SetOutputSpacing( spacing );
  cout << "Spacing = ";
  cout << sp[0] << ", " << sp[1] << ", " << sp[2] << std::endl;
  

  //----- set the origin of the filter
  double origin[ Dimension ];
  const InputImageType::PointType& orgn = reader->GetOutput()->GetOrigin();
  origin[0] = orgn[0];  // X space coordinate of origin
  origin[1] = orgn[1];  // Y space coordinate of origin
  origin[2] = orgn[2];  // Y space coordinate of origin
  filter->SetOutputOrigin( origin );
  cout << "Origin = ";
  cout << orgn[0] << ", " << orgn[1] << ", " << orgn[2] << std::endl;

  //----- set the size of the filter
  InputImageType::SizeType size, imageSize =reader->GetOutput()->GetLargestPossibleRegion().GetSize();  
  size[0] = imageSize[0];  // number of pixels along X
  size[1] = imageSize[1];  // number of pixels along Y
  size[2] = imageSize[2];  // number of pixels along Y
  filter->SetSize( size ); 
  std::cout << "ImageSize= ";
  std::cout << imageSize[0] << ", " << imageSize[1] << ", " << imageSize[2]<< std::endl;

  //// Image Translation
  //TransformType::OutputVectorType translation;
  //translation[0] = -30;  // X translation in millimeters
  //translation[1] = -50;  // Y translation in millimeters
  //translation[2] = -50;  // Y translation in millimeters
  //transform->Translate( translation );
  //filter->SetDefaultPixelValue( 100 );
  ////filter->SetDefaultPixelValue( 0 );


  //// Image Scaling
  //TransformType::OutputVectorType scale;
  //scale[0] = 1;  // X translation in millimeters
  //scale[1] = 1;  // Y translation in millimeters
  //scale[2] = 1;  // Y translation in millimeters
  //transform->Scale( scale );
  //filter->SetDefaultPixelValue( 0 );
  //filter->SetDefaultPixelValue( 100 );

 

  filter->SetInput( reader->GetOutput() );
  writer->SetInput( filter->GetOutput() );
  writer->Update();
  

  transform->Print( std::cout );

  cin.get();
  return EXIT_SUCCESS;
}