#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkIdentityTransform.h"
#include "itkNearestNeighborInterpolateImageFunction.h"
#include "itkAffineTransform.h"

using namespace std;
int main( void )
{  
  typedef   unsigned char  InputPixelType;
  typedef   unsigned char  OutputPixelType;
  const     unsigned int   Dimension = 3;

  typedef itk::Image< InputPixelType,  Dimension >   InputImageType;
  typedef itk::Image< OutputPixelType, Dimension >   OutputImageType;

  typedef itk::ImageFileReader< InputImageType   >  InputReaderType;
  typedef itk::ImageFileReader< OutputImageType  >  OutputReaderType;
  typedef itk::ImageFileWriter< OutputImageType  >  WriterType;
  typedef itk::ResampleImageFilter< OutputImageType, OutputImageType >  FilterType;
  typedef itk::IdentityTransform< double, Dimension >  TransformType;
  typedef itk::NearestNeighborInterpolateImageFunction< OutputImageType, double >  InterpolatorType;
  
  InputReaderType::Pointer  InputReader  = InputReaderType::New();
  OutputReaderType::Pointer OutputReader = OutputReaderType::New();
  WriterType::Pointer writer = WriterType::New();  
  FilterType::Pointer resampler = FilterType::New();
  TransformType::Pointer transform = TransformType::New();
  InterpolatorType::Pointer interpolator = InterpolatorType::New();

  typedef itk::AffineTransform< double, Dimension >  AffineTransformType; 
  AffineTransformType::Pointer Affinetransform = AffineTransformType::New();

  //working
  InputReader->SetFileName("C:/Images/Charbag_04000200001060_080530121234.hdr");
  OutputReader->SetFileName("C:/Images/brainweb165a10f17.mha");  
  writer->SetFileName("C:/Images/Volume_IdentityTransform.hdr");

  // //not working
  //InputReader->SetFileName("C:/Images/brainweb165a10f17.mha");
  //OutputReader->SetFileName("C:/Images/ConfiConnect_brainweb165a10f17_whiteMatter.hdr");  
  //writer->SetFileName("C:/Images/Volume_IdentityTransform.hdr");

  //  get output volume params
  try 
    {
    OutputReader->Update();
    }
  catch( itk::ExceptionObject & excep )
    {
    std::cerr << "Exception catched !" << std::endl;
    std::cerr << excep << std::endl;
    }

  OutputImageType::ConstPointer OutputImage = OutputReader->GetOutput();
  OutputImageType::RegionType  region2 = OutputImage->GetLargestPossibleRegion();  
  OutputImageType::SizeType     size2  = region2.GetSize();
  OutputImageType::IndexType    start2 = region2.GetIndex();
  const OutputImageType::PointType& orgn2 = OutputImage->GetOrigin(); 
  const OutputImageType::SpacingType& sp2 = OutputImage->GetSpacing(); 

   //  get Input volume params
  try 
    {
    InputReader->Update();
    }
  catch( itk::ExceptionObject & excep )
    {
    std::cerr << "Exception catched !" << std::endl;
    std::cerr << excep << std::endl;
    }

  InputImageType::ConstPointer InputImage = InputReader->GetOutput();
  InputImageType::RegionType  region = InputImage->GetLargestPossibleRegion();  
  InputImageType::SizeType     size  = region.GetSize();
  InputImageType::IndexType    start = region.GetIndex();
  const InputImageType::PointType& orgn = InputImage->GetOrigin(); 
  const InputImageType::SpacingType& sp = InputImage->GetSpacing(); 

  AffineTransformType::OutputVectorType scale;
  scale[0] = 1;  // X scaling in millimeters
  scale[1] = 1;  // Y scaling in millimeters
  scale[2] = 1;  // Z scaling in millimeters
  //Affinetransform->Scale( scale );
  Affinetransform->SetIdentity();  


  resampler->SetDefaultPixelValue( 0 );
  resampler->SetOutputSpacing(sp);
  resampler->SetOutputOrigin(orgn);
  resampler->SetSize( size );
  resampler->SetOutputStartIndex( start );

  resampler->SetTransform( Affinetransform );
  resampler->SetInterpolator( interpolator );
  resampler->SetInput( OutputReader->GetOutput() );
  
  writer->SetInput( resampler->GetOutput() );
  try 
    {
    writer->Update();
    }
  catch( itk::ExceptionObject & excep )
    {
    std::cerr << "Exception catched !" << std::endl;
    std::cerr << excep << std::endl;
    }
 

  double Physical_ext[Dimension],Physical_ext2[Dimension];
  Physical_ext[0]=size[0]*sp[0];  Physical_ext[1]=size[1]*sp[1];  Physical_ext[2]=size[2]*sp[2];
  Physical_ext2[0]=size2[0]*sp2[0];  Physical_ext2[1]=size2[1]*sp2[1];  Physical_ext2[2]=size2[2]*sp2[2];

  cout << "Size= "; cout << size[0] << ", " <<  size[1] << ", " <<  size[2] << std::endl;
  cout << "Start= "; cout << start[0] << ", " <<  start[1] << ", " <<  start[2] << std::endl; 
  cout << "Origin= "; cout <<orgn[0] << ", " << orgn[1] << ", " << orgn[2] << std::endl;
  cout << "Spacing= "; cout <<sp[0] << ", " << sp[1] << ", " << sp[2] << std::endl;
  cout << "Physical_extent= "; cout <<Physical_ext[0]<< ", " << Physical_ext[1] << ", " << Physical_ext[2] << std::endl;
  
  cout <<std::endl; cout <<std::endl;

  cout << "Size Out= "; cout << size2[0] << ", " <<  size2[1] << ", " <<  size2[2] << std::endl;
  cout << "Start Out= "; cout << start2[0] << ", " <<  start2[1] << ", " <<  start2[2] << std::endl; 
  cout << "Origin Out= "; cout <<orgn2[0] << ", " << orgn2[1] << ", " << orgn2[2] << std::endl;
  cout << "Spacing Out= "; cout <<sp2[0] << ", " << sp2[1] << ", " << sp2[2] << std::endl;
  cout << "Physical_extent= "; cout <<Physical_ext2[0]<< ", " << Physical_ext2[1] << ", " << Physical_ext2[2] << std::endl;
  
  cout << "Code end !" << std::endl;
  cin.get();
  return 0;
}