/*
 * Copyright (c) ICG. All rights reserved.
 * See copyright.txt for more information.
 *
 * Institute for Computer Graphics and Vision
 * Graz, University of Technology / Austria
 *
 *
 * This software is distributed WITHOUT ANY WARRANTY; without even
 * the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
 * PURPOSE.  See the above copyright notices for more information.
 *
 *
 * Project     : MIPItkProjects
 * Module      : Evaluation
 * Class       : $RCSfile: ApplyUniformPeriodicTransform.cxx,v $
 * Language    : C++
 * Description : 
 *
 * Author     : Martin Urschler
 * EMail      : urschler@icg.tu-graz.ac.at
 * Date       : $Date: 2007-06-26 15:35:29 $
 * Version    : $Revision: 1.2 $
 * Full Id    : $Id: ApplyUniformPeriodicTransform.cxx,v 1.2 2007-06-26 15:35:29 urschler Exp $
 *
 */

#include "commonItkTypedefs.h"
#include "SmallUtilityMethods.h"
#include "VolumeIOWrapper.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkImageFileWriter.h"




int main( int argc, char * argv[] )
{
   if( argc != 6 )
   {
       std::cerr << "Missing Parameters " << std::endl;
       std::cerr << "Usage: " << argv[0];
       std::cerr << "  fixedImage  warpedImage  amplitude  phase" <<
          "  resampleDefaultValue" << std::endl;
       return EXIT_FAILURE;
   }
    
   std::string inputFilename( argv[1] );
   std::string outputFilename( argv[2] );
   const int amplitude = atoi( argv[3] );
   const int phase = atoi(argv[4]);
   const Int16ImageType::PixelType resampleDefaultValue = atoi( argv[5] );
   
   
   
   Int16ImageType::Pointer originalImage =
      VolumeIOWrapper<Int16ImageType>::readITKVolume( inputFilename, "" );

   Int16ImageType::SizeType originalImageSize = 
      originalImage->GetLargestPossibleRegion().GetSize();

   Int16ImageType::SpacingType spacing = 
      originalImage->GetSpacing();
      
   typedef itk::LinearInterpolateImageFunction<Int16ImageType, double>
      InterpolatorType;
   InterpolatorType::Pointer interpolator = InterpolatorType::New();
   interpolator->SetInputImage(originalImage);
   
   Int16ImageType::Pointer outputImage = 
      SmallUtilityMethods<Int16ImageType>::createNewImage( originalImage );
   itk::ImageRegionIteratorWithIndex<Int16ImageType> outIter(
      outputImage, outputImage->GetLargestPossibleRegion());
   
   // create a result deformation field
   DeformationFieldType::Pointer deformationField = 
      SmallUtilityMethods<Int16ImageType>::createNewZeroedDeformationField( 
         outputImage );
   
   itk::ImageRegionIterator<DeformationFieldType> fieldIter( deformationField, 
      deformationField->GetLargestPossibleRegion() );
   
   
   Int16ImageType::IndexType outputIndex;
   Int16ImageType::PointType outputPoint;
   Int16ImageType::PointType inputResampleLocation;
   VectorPixelType deformation;
   
   while ( !outIter.IsAtEnd() )
   {
      outputIndex = outIter.GetIndex();
      outputImage->TransformIndexToPhysicalPoint( outputIndex, 
         outputPoint );
      
      // Compute corresponding displaced point 
      inputResampleLocation[0] = outputPoint[0] + 
         (amplitude*spacing[0]) * vcl_sin(outputPoint[0] * vnl_math::pi / (phase*spacing[0]));
      inputResampleLocation[1] = outputPoint[1] + 
         (amplitude*spacing[1]) * vcl_sin(outputPoint[1] * vnl_math::pi / (phase*spacing[1]));
      inputResampleLocation[2] = outputPoint[2]; 
         (amplitude*spacing[2]) * vcl_sin(outputPoint[2] * vnl_math::pi / (phase*spacing[2]));
      
      //std::cout << outputPoint << " -> " << inputResampleLocation << std::endl;
      
      if (interpolator->IsInsideBuffer(inputResampleLocation))
      {
         outIter.Set( ROUND_DBL<Int16ImageType::PixelType>(
            interpolator->Evaluate(inputResampleLocation)) );
            
         VectorPixelType deformation = inputResampleLocation - outputPoint;
         fieldIter.Set( deformation );
      }
      else
      {
         outIter.Set( resampleDefaultValue );
         fieldIter.Set( 
            VectorPixelConstants<VectorPixelType::ValueType, 3>::undefined() );
      }
      
      ++outIter;
      ++fieldIter;
   }
   
   // file output
   
   VolumeIOWrapper<Int16ImageType>::writeITKVolume16Bit( 
      outputImage, outputFilename, "" );
   
   
   std::string displacementFieldFilename = outputFilename;
   displacementFieldFilename.replace(
      displacementFieldFilename.find(".hdr"), 4, ".mha" ); 
   typedef itk::ImageFileWriter<DeformationFieldType> DeformationFieldWriter;
   DeformationFieldWriter::Pointer writer = DeformationFieldWriter::New();
   writer->SetFileName( displacementFieldFilename.c_str() );
   writer->SetInput( deformationField );
   std::cout << "Writing deformation field" << 
      displacementFieldFilename.c_str() << std::endl;
   writer->Update();

   return EXIT_SUCCESS;
}