#ifndef __itkPhaseCorrelationImageRegistrationMethodTest_h
#define __itkPhaseCorrelationImageRegistrationMethodTest_h


#include "itkArray.h"
#include "itkPhaseCorrelationImageRegistrationMethodImageSource.h"
#include "itkResampleImageFilter.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkPhaseCorrelationImageRegistrationMethod.h"
#include "itkMaxPhaseCorrelationOptimizer.h"
#include "itkImageFileWriter.h"
#include "itkCastImageFilter.h"

/**
 *  This program tests one instantiation of the
 *  itk::PhaseCorrelationImageRegistrationMethod class with
 *  itk::PhaseCorrelationOperator and itk::MaxPhaseCorrelationOptimizer.
 *
 * \author Jakub Bican, jakub.bican@matfyz.cz, Department of Image Processing, 
 *         Institute of Information Theory and Automation, 
 *         Academy of Sciences of the Czech Republic.
 *
 */
template < unsigned int Dimension,
           typename TFixedImagePixelType,
					 typename TMovingImagePixelType >
int itkPhaseCorrelationImageRegistrationMethodTest(int argc, char* argv[] )
{

  bool pass = true;

  // Fixed Image Type
  typedef itk::Image<TFixedImagePixelType,Dimension>   FixedImageType;

  // Moving Image Type
  typedef itk::Image<TMovingImagePixelType,Dimension>  MovingImageType;

  // Size Type
  typedef typename MovingImageType::SizeType           SizeType;

  // test image source
  typedef itk::testhelper::PhaseCorrelationImageRegistrationMethodImageSource<
      typename FixedImageType::PixelType,
      typename MovingImageType::PixelType,
      Dimension >                                      ImageSourceType;

  // Resample filter
  typedef itk::ResampleImageFilter< MovingImageType, MovingImageType >
                                                       ResamplerType;

  // Interpolator for resampler
  typedef itk::LinearInterpolateImageFunction< MovingImageType, double >
                                                       InterpolatorType;


  // Registration method
  typedef itk::PhaseCorrelationImageRegistrationMethod< FixedImageType,
                                                        MovingImageType >
                                                       PCMType;

  // Operator type
  typedef itk::PhaseCorrelationOperator< PCMType >     OperatorType;

  // Optimizer type
  typedef itk::MaxPhaseCorrelationOptimizer<PCMType>   OptimizerType;

  // Transform type
  typedef typename PCMType::TransformType              TransformType;
  typedef typename TransformType::ParametersType       ParametersType;


  typename OperatorType::Pointer       pcmOperator   = OperatorType::New();
  typename OptimizerType::Pointer      pcmOptimizer  = OptimizerType::New();
  typename PCMType::Pointer            pcm           = PCMType::New();

  typename ImageSourceType::Pointer    imageSource   = ImageSourceType::New();
  typename ResamplerType::Pointer      resampler     = ResamplerType::New();
  typename InterpolatorType::Pointer   interpolator  = InterpolatorType::New();

  SizeType size;
  double spacing[ Dimension ];
  double origin[ Dimension ];
  typename MovingImageType::SizeType   newMovingSize;
  typename MovingImageType::PointType  newMovingOrigin;
  for (unsigned int i = 0; i < Dimension; i++)
    {
    size[i] = 100;
    spacing[i] = 1.0;
    origin[i] = 0.0;
    newMovingSize[i] = 100;
    newMovingOrigin[i] = 0.0;
    }

  // increase the resolution and size and crop moving image in 1st dimension
  // this tests the ability of PCM to padd the images to the same real size
  // and to resample the images to the same pixel size and spacing
  spacing[0] = 0.8;
  newMovingSize[0] = ( 100.0 / spacing[0] ) - 10;

  imageSource->GenerateImages( size );

  resampler->SetInterpolator( interpolator );
  resampler->SetDefaultPixelValue( 0 );
  resampler->SetOutputSpacing( spacing );
  resampler->SetOutputOrigin( origin );
  resampler->SetSize( newMovingSize );

  resampler->SetInput( imageSource->GetMovingImage() );

  resampler->Update();

  typename FixedImageType::ConstPointer  fixedImage    
                                            = imageSource->GetFixedImage();
  typename MovingImageType::Pointer      movingImage   
                                            = resampler->GetOutput();

  // shift the origin of the moving image in 2nd dimension
  // this tests the ability of PCM to introduce between-image origin offset
  // into the transformation (so the final parameters can be directly used to
  // resample the two images into the same coordinate system)
  // ! supposing that the input images have all origin components == 0.0 !
  newMovingOrigin[1] = 2.0;
  movingImage->SetOrigin(newMovingOrigin);


  //
  // Connect all the components required for Registratio
  //
  pcm->SetOperator(    pcmOperator  );
  pcm->SetOptimizer(   pcmOptimizer );
  pcm->SetFixedImage(  fixedImage   );
  pcm->SetMovingImage( movingImage  );


  //
  // Execute the registration.
  // This can potentially throw an exception
  //
  try
    {
    pcm->Update();
    }
  catch( itk::ExceptionObject & e )
    {
    std::cerr << e << std::endl;
    pass = false;
    }

  //
  // Get registration result and validate it.
  //
  ParametersType actualParameters    = imageSource->GetActualParameters();
  ParametersType finalParameters     = pcm->GetTransformParameters();
  ParametersType transformParameters = pcm->GetOutput()->Get()->GetParameters();

  const unsigned int numberOfParameters = actualParameters.Size();

  const double tolerance = 1.0;  // equivalent to 1 pixel.

  // Validate first two parameters (introduced by image source)
  for(unsigned int i=0; i<numberOfParameters; i++)
    {
    // the parameters are negated in order to get the inverse transformation.
    // this only works for comparing translation parameters....
    std::cout << finalParameters[i] << " == "
              << -(actualParameters[i]+newMovingOrigin[i]) << " == "
              << transformParameters[i] << std::endl;

    if(  ( vnl_math_abs ( finalParameters[i] -
                          (-(actualParameters[i]+newMovingOrigin[i]))
                        ) > tolerance ) ||
         ( vnl_math_abs ( transformParameters[i] -
                          (-(actualParameters[i]+newMovingOrigin[i]))
                        ) > tolerance ) )
      {
      std::cout << "Tolerance exceeded at component " << i << std::endl;
      pass = false;
      }
    }

  // All other parameters must be 0
  for (unsigned int i=numberOfParameters; i<Dimension; i++)
    {
    if (  ( vnl_math_abs ( finalParameters[i] ) > tolerance )
        ||
          ( vnl_math_abs ( finalParameters[i] ) > tolerance ) )
      {
      std::cout << "Tolerance exceeded at component " << i << std::endl;
      pass = false;
      }
    }

  if( !pass )
    {
    std::cout << "Test FAILED." << std::endl;
    return EXIT_FAILURE;
    }

  std::cout << "Test PASSED." << std::endl;
  return EXIT_SUCCESS;

}

#endif