/*=========================================================================
file  : ITKInterpolate.cxx
usage : ./ITKInterpolate small.vtk junk.vtk 0.98 0.98 0.98
desc. : This is an example of using ITK to perform linear interpolation of
        3D data.

        According to
        http://www.itk.org/Doxygen16/html/classitk_1_1ResampleImageFilter.html,
        ResampleImageFilter is multithreaded.

author: george j. grevera, ph.d.
date  : 7-july-2006
=========================================================================*/
#if defined(_MSC_VER)
    #pragma warning ( disable : 4786 )
#endif

#ifdef __BORLANDC__
    #define ITK_LEAN_AND_MEAN
#endif

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

int main ( int argc, char ** argv ) {
  std::cerr << "warning: "
            << "ASCII vtk input files are not read properly by itk 2.8.1."
            << std::endl;
  // verify the number of parameters in the command line
  if (argc<3) {
      std::cerr << "Usage: " << std::endl
                << argv[0] << " inputImageFile outputImageFile"
                << " [xSpace ySpace zSpace]" << std::endl
                << "    where xSpace ySpace zSpace are 3 floating point "
                << "numbers respresenting the desired output size in mm."
                << std::endl;
      return EXIT_FAILURE;
  }

  //define the dimensionality of the data as well as the data types.
  const unsigned int  Dimension = 3;

  //typedef unsigned short                            InputPixelType;
  typedef float                                     InputPixelType;
  typedef itk::Image< InputPixelType, Dimension >   InputImageType;

  typedef unsigned short      OutputPixelType;
  typedef itk::Image< OutputPixelType, Dimension >  OutputImageType;

  //create an image file reader
  typedef itk::ImageFileReader< InputImageType >  ReaderType;
  ReaderType::Pointer   reader = ReaderType::New();
  const char* const  inputFilename  = argv[1];
  reader->SetFileName( inputFilename  );
  try {
      reader->Update();
  } catch ( itk::ExceptionObject& excep ) {
      std::cerr << "Exception caught!" << std::endl;
      std::cerr << excep << std::endl;
  }
  //get the size and spacing of the input data
  const InputImageType::SpacingType&
      inputSpacing = reader->GetOutput()->GetSpacing();
  InputImageType::SizeType
      inputSize = reader->GetOutput()->GetLargestPossibleRegion().GetSize();

  //create a resample image filter to interpolate our data
  typedef itk::ResampleImageFilter< InputImageType, OutputImageType >
      ResampleFilterType;
  ResampleFilterType::Pointer  resampler = ResampleFilterType::New();
  typedef itk::IdentityTransform< double, Dimension >  TransformType;
  TransformType::Pointer  transform = TransformType::New();
  transform->SetIdentity();
  resampler->SetTransform( transform );
  //linear interpolation
  typedef itk::LinearInterpolateImageFunction< InputImageType, double >
          InterpolatorType;
  InterpolatorType::Pointer  interpolator = InterpolatorType::New();
  resampler->SetInterpolator( interpolator );
  resampler->SetDefaultPixelValue( 255 ); //for regions without source
  resampler->SetOutputOrigin( reader->GetOutput()->GetOrigin() );

  //define the spacing of the result
  OutputImageType::SpacingType  outputSpacing;
  if (argc==6) {
      char*  endptr;
      outputSpacing[0] = strtod( argv[3], &endptr );
      if (endptr==argv[3]) {
          std::cerr << "Usage: " << std::endl
                    << argv[0] << " inputImageFile outputImageFile"
                    << " [xSpace ySpace zSpace]" << std::endl;
          return EXIT_FAILURE;
      }
      outputSpacing[1] = strtod( argv[4], &endptr );
      if (endptr==argv[4]) {
          std::cerr << "Usage: " << std::endl
                    << argv[0] << " inputImageFile outputImageFile"
                    << " [xSpace ySpace zSpace]" << std::endl;
          return EXIT_FAILURE;
      }
      outputSpacing[2] = strtod( argv[5], &endptr );
      if (endptr==argv[5]) {
          std::cerr << "Usage: " << std::endl
                    << argv[0] << " inputImageFile outputImageFile"
                    << " [xSpace ySpace zSpace]" << std::endl;
          return EXIT_FAILURE;
      }
  } else {
      outputSpacing[0] = inputSpacing[0];
      outputSpacing[1] = inputSpacing[1];
      outputSpacing[2] = inputSpacing[2];
  }
  std::cout << "input spacing = (" << inputSpacing[0] << ","
            << inputSpacing[1] << "," << inputSpacing[2] << ")" << std::endl
            << "output spacing = (" << outputSpacing[0] << ","
            << outputSpacing[1] << "," << outputSpacing[2] << ")" << std::endl;
  resampler->SetOutputSpacing( outputSpacing );
  //define the size of the result
  const double  dx =  inputSize[0]      * inputSpacing[0] / outputSpacing[0];
  const double  dy =  inputSize[1]      * inputSpacing[1] / outputSpacing[1];
  const double  dz = (inputSize[2] - 1) * inputSpacing[2] / outputSpacing[2];
  OutputImageType::SizeType  outputSize;
  typedef OutputImageType::SizeType::SizeValueType  SizeValueType;
  outputSize[0] = static_cast<SizeValueType>( dx );
  outputSize[1] = static_cast<SizeValueType>( dy );
  outputSize[2] = static_cast<SizeValueType>( dz );
  resampler->SetSize( outputSize );
  resampler->SetInput( reader->GetOutput() );
  resampler->Update();

  //create the image file writer
  typedef itk::ImageFileWriter< OutputImageType >  WriterType;
  WriterType::Pointer   writer = WriterType::New();
  const char* const  outputFilename = argv[2];
  writer->SetFileName( outputFilename );
#if 1
  std::cout << "writer's input is resampler's output" << std::endl;
  writer->SetInput( resampler->GetOutput() );
#else
  std::cout << "writer's input is reader's output" << std::endl;
  writer->SetInput( reader->GetOutput() );
#endif
  //write the file
  try { 
      writer->Update(); 
  } catch( itk::ExceptionObject & err ) { 
      std::cerr << "ExceptionObject caught !" << std::endl << err << std::endl; 
      return EXIT_FAILURE;
  } 

#if 0
  InputImageType::SizeType
      inputSize = reader->GetOutput()->GetLargestPossibleRegion().GetSize();
  for (int z=0; z<inputSize[2]; z++) {
      for (int y=0; y<inputSize[1]; y++) {
          for (int x=0; x<inputSize[0]; x++) {
              InputImageType::IndexType  index;
              index[0] = x;
              index[1] = y;
              index[2] = z;
              int    p = reader->GetOutput()->GetPixel( index );
              std::cout << p << " ";
          }
      }
  }
#endif

  return EXIT_SUCCESS;
}