/*
* Copyright (c) ICG. All rights reserved.
*
* 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     : projects
* Module      : MIPViewer - plugins
* Class       : $RCSfile: CTFastBlockMatchingRegistration.txx,v $
* Language    : C++
* Description : 
*
* Author     : Martin Urschler
* EMail      : urschler@icg.tu-graz.ac.at
* Date       : $Date: 2007-07-11 13:54:43 $
* Version    : $Revision: 1.6 $
* Full Id    : $Id: CTFastBlockMatchingRegistration.txx,v 1.6 2007-07-11 13:54:43 urschler Exp $
*
*/


#ifndef _CTFASTBLOCKMATCHINGREGISTRATION_TXX_
#define _CTFASTBLOCKMATCHINGREGISTRATION_TXX_

#include "SmallUtilityMethods.h"
#include "NonlinearRegistrationUtilities.h"
#include "VolumeIOWrapper.h"

#include "itkLinearInterpolateImageFunction.h"
#include "itkImageFileWriter.h"
#include "insight-journal/suarez/itkVectorDivideImageFilter.h"
#include "insight-journal/suarez/itkNanWarpImageFilter.h"
#include "insight-journal/suarez/itkSuarezBlockMatchingRegistration.h"




template< typename TPixelDataType >
CTNonlinearRegistration::ReturnType 
CTFastBlockMatchingRegistration<TPixelDataType>::
matching(
   const float regularization,
   const float similarity,
   const unsigned int levelsNotToCompute )
{
   typedef InputImageType ScalarImageType;
   typedef DeformationFieldType VectorImageType;
   
   
   typedef itk::SuarezBlockMatchingRegistration<ScalarImageType,
      VectorImageType> FilterType;
   typename FilterType::Pointer filter = FilterType::New();
   filter->SetFixedImage( this->m_fixedImage );
   filter->SetMovingImage( this->m_movingImage );
   filter->SetVarianceRegularization( regularization );
   filter->SetVarianceSimilarity( similarity );
   filter->SetNumberOfLevelsNotToCompute( levelsNotToCompute );

   typedef itk::LinearInterpolateImageFunction<ScalarImageType, double>
      LinearInterpolatorType; // TODO: can't take 'double' off with gcc
   typename LinearInterpolatorType::Pointer linearInterpolator =
      LinearInterpolatorType::New();

   typedef itk::NanWarpImageFilter
      <ScalarImageType, ScalarImageType, VectorImageType> WarpType;
   typename WarpType::Pointer warp = WarpType::New();

   EXCEPTION_CHECKED_ITK_FILTER_UPDATE( "registration",
      filter->Update() );
  
   VectorImageType::Pointer deformation = filter->GetOutput();

   warp->SetInterpolator( linearInterpolator );
   warp->SetEdgePaddingValue(
      itk::NumericTraits<typename WarpType::PixelType>::quiet_NaN() );
   warp->SetOutputSpacing( deformation->GetSpacing() );
   warp->SetOutputOrigin( deformation->GetOrigin() );
   warp->SetDeformationField( deformation );
   warp->SetInput( this->m_movingImage );
   warp->SetEdgePaddingValue( this->m_resampleDefaultValue );
   
   EXCEPTION_CHECKED_ITK_FILTER_UPDATE( "final image warping",  
      warp->Update() );
   
   typename InputImageType::Pointer warpedImage = warp->GetOutput();

   // finally copy the warped output image back into the moving image
   ITK_IMAGE_COPY_MACRO_TN( InputImageType, warpedImage, 
                            InputImageType, this->m_movingImage );
   
   // if necessary write the displacement field
   if (this->m_displacementFieldFilename != std::string(""))
   {
      typedef itk::ImageFileWriter<DeformationFieldType> DefFieldWriter;
      DefFieldWriter::Pointer def_field_writer = DefFieldWriter::New();
      def_field_writer->SetFileName(this->m_displacementFieldFilename.c_str());
      def_field_writer->SetInput( deformation );
      std::cout << "writing displacement field: " << 
         this->m_displacementFieldFilename.c_str() << std::endl;
      def_field_writer->Update();
   }
   
   return CTNonlinearRegistration::OK;
}





#endif /*CTFASTBLOCKMATCHINGREGISTRATION_TXX_*/