/*
 * 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: ApplyDisplacementFieldTransform.cpp,v $
 * Language    : C++
 * Description : 
 *
 * Author     : Martin Urschler
 * EMail      : urschler@icg.tu-graz.ac.at
 * Date       : $Date: 2007-07-11 16:44:29 $
 * Version    : $Revision: 1.11 $
 * Full Id    : $Id: ApplyDisplacementFieldTransform.cpp,v 1.11 2007-07-11 16:44:29 urschler Exp $
 *
 */



// include files

#include "commonItkTypedefs.h"
#include "SmallUtilityMethods.h"
#include "VolumeIOWrapper.h"
#include "itkWarpImageFilter.h"

#include "itkImageRegionIteratorWithIndex.h"
#include "itkImageFileWriter.h"



int main( int argc, char *argv[] )
{
   if (argc < 5 || argc > 6)
   {
      std::cerr << "Usage: " << std::endl;
      std::cerr << argv[0] << 
         "  original_volume  output_volume  input_deformation_field" <<
         "  resample_default_value  [percentage betw. 0..100]" << 
         std::endl;
      return EXIT_FAILURE;
   }
   
   std::string input_filename( argv[1] );
   std::string input_displacement_filename( argv[3] );
   std::string output_filename( argv[2] );
   const Int16ImageType::PixelType warping_default_value = atoi( argv[4] );
   
   float percentage = 1.0f;
   if (argc == 6)
   {
      percentage = static_cast<float>(atoi(argv[5])) / 100.0f; 
   }
   
   Int16ImageType::Pointer original_image =
      VolumeIOWrapper<Int16ImageType>::readITKVolume( input_filename, "" );
   
   typedef itk::ImageFileReader<DeformationFieldType> 
      DeformationFieldReaderType;
   DeformationFieldReaderType::Pointer reader = 
      DeformationFieldReaderType::New();
   reader->SetFileName( input_displacement_filename.c_str() );
   reader->Update();
   
   DeformationFieldType::Pointer disp_field =
      reader->GetOutput();
   
   
   // set up the physical dimensions of the translated output image
   // they are identical with the input image's physical dimensions
   const Int16ImageType::SpacingType output_spacing = 
      original_image->GetSpacing();
   const Int16ImageType::PointType output_origin = 
      original_image->GetOrigin();
   // the output size will be taken from the output deformation field!
   
   
   if (percentage == 1.0f)
   {
      typedef itk::WarpImageFilter<Int16ImageType, Int16ImageType, 
         DeformationFieldType> WarperType;
      
      WarperType::Pointer warper = WarperType::New();
      
      typedef itk::LinearInterpolateImageFunction< 
                          Int16ImageType, double >  InterpolatorType;
      InterpolatorType::Pointer interpolator = InterpolatorType::New();
   
      warper->SetInterpolator( interpolator );
      warper->SetOutputSpacing( output_spacing );
      warper->SetOutputOrigin(  output_origin );
      warper->SetDeformationField( disp_field );
      warper->SetInput( original_image );
      warper->SetEdgePaddingValue( warping_default_value );
      
      warper->Update();
      
      VolumeIOWrapper<Int16ImageType>::writeITKVolume16Bit( 
            warper->GetOutput(), output_filename, "" );
   }
   else
   {
      typedef itk::LinearInterpolateImageFunction<Int16ImageType, double>
         MovingImageInterpolatorType;
      
      const unsigned int ImageDimension =
         itk::GetImageDimension<DeformationFieldType>::ImageDimension;
      
      Int16ImageType::SizeType size =
         original_image->GetLargestPossibleRegion().GetSize();
      
      ITK_IMAGE_ALLOCATE_MACRO( Int16ImageType, warpedImage,
         Int16ImageType, original_image );
      
      // warp the volume according to def
      itk::ImageRegionIteratorWithIndex<DeformationFieldType>
         defIter( disp_field, disp_field->GetLargestPossibleRegion() );
      itk::ImageRegionIterator<Int16ImageType> warpIter(
         warpedImage, warpedImage->GetLargestPossibleRegion());
         
      defIter.GoToBegin();
      warpIter.GoToBegin();
      
      DeformationFieldType::IndexType index;
      DeformationFieldType::PointType point;
      DeformationFieldType::PixelType currentDef;
   
      MovingImageInterpolatorType::Pointer movingImageInterpolator =
         MovingImageInterpolatorType::New();
      movingImageInterpolator->SetInputImage( original_image );
      
      Int16ImageType::PointType upperBorder;
      for (unsigned int dim=0; dim<ImageDimension; ++dim)
      {
         upperBorder[dim] = (size[dim]-1) * output_spacing[dim];
      }
   
      while(!defIter.IsAtEnd())
      {
         index = defIter.GetIndex();
   
         original_image->TransformIndexToPhysicalPoint(index, point);
   
         currentDef = defIter.Get();
   
         for (unsigned int dim=0; dim<ImageDimension; ++dim)
         {
            currentDef[dim] = percentage * currentDef[dim]; 
            
            point[dim] += currentDef[dim];
            // clamp point to the image borders
            point[dim] = (point[dim] < 0.0f) ? 0.0f : point[dim];
            point[dim] = (point[dim] > upperBorder[dim]) ?
                         upperBorder[dim] : point[dim];
         }
   
         const Int16ImageType::PixelType movingValue = 
            static_cast<Int16ImageType::PixelType>(
               movingImageInterpolator->Evaluate(point) );
   
         warpIter.Set( movingValue );
         
         defIter.Set( currentDef );
   
         ++defIter;
         ++warpIter;
      }
      
      VolumeIOWrapper<Int16ImageType>::writeITKVolume16Bit( 
         warpedImage, output_filename, "" );
   }
   
   std::string displacement_field_filename = output_filename;
   displacement_field_filename.replace(
      displacement_field_filename.find(".hdr"), 4, ".mha" ); 
   typedef itk::ImageFileWriter<DeformationFieldType> DeformationFieldWriter;
   DeformationFieldWriter::Pointer writer = DeformationFieldWriter::New();
   writer->SetFileName( displacement_field_filename.c_str() );
   writer->SetInput( disp_field );
   std::cout << "Writing deformation field!" << std::endl;
   writer->Update();

      
   return EXIT_SUCCESS;
}