/*
* 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: CTCurvatureRegistration.txx,v $
* Language    : C++
* Description : 
*
* Author     : Martin Urschler
* EMail      : urschler@icg.tu-graz.ac.at
* Date       : $Date: 2006-10-31 09:00:24 $
* Version    : $Revision: 1.4 $
* Full Id    : $Id: CTCurvatureRegistration.txx,v 1.4 2006-10-31 09:00:24 urschler Exp $
*
*/


#ifndef _CT_CURVATURE_REGISTRATION_TXX_
#define _CT_CURVATURE_REGISTRATION_TXX_

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

#include "itkCurvatureRegistrationFilter.h"
#include "itkSymmetricForcesDemonsRegistrationFunction.h"

template< typename TPixelDataType >
void 
CTCurvatureRegistration<TPixelDataType>::
CurvatureCommandIterationUpdate::
Execute(itk::Object *caller, const itk::EventObject & event)
{
   Execute( (const itk::Object *)caller, event);
}

template< typename TPixelDataType >
void 
CTCurvatureRegistration<TPixelDataType>::
CurvatureCommandIterationUpdate::
Execute(const itk::Object * object, const itk::EventObject & event)
{
   typedef itk::SymmetricForcesDemonsRegistrationFunction<
      InputImageType,
      InputImageType,
      DeformationFieldType > RegistrationFunctionType;
      
   typedef itk::CurvatureRegistrationFilter< 
         InputImageType, 
         InputImageType,
         DeformationFieldType,
         RegistrationFunctionType > DeformableRegistrationFilterType;
         
   const DeformableRegistrationFilterType * filter = 
      dynamic_cast< const DeformableRegistrationFilterType * >( object );
   if( !(itk::IterationEvent().CheckEvent( &event )) )
   {
      return;
   }
         
   std::cout << filter->GetMetric() << std::endl;
}


template< typename TPixelDataType >
CTNonlinearRegistration::ReturnType 
CTCurvatureRegistration<TPixelDataType>::
matching( const unsigned int maxNrShrinkLevels, 
   const unsigned int iterationsParam1, 
   const unsigned int iterationsParam2 )
{  
   typedef itk::SymmetricForcesDemonsRegistrationFunction<
      InputImageType,
      InputImageType,
      DeformationFieldType > RegistrationFunctionType;
       
   typedef itk::CurvatureRegistrationFilter<
      InputImageType,
      InputImageType,
      DeformationFieldType,
      RegistrationFunctionType > CurvatureRegistrationFilterType;
   
   // here the multi resolution demons registration algorithm starts   
   
   // calculate the pyramid
   PyramidType fixedPyramid(maxNrShrinkLevels+1);
   PyramidType movingPyramid(maxNrShrinkLevels+1);
   CTNonlinearRegistration::ReturnType pyramidReturnValue =
      calculateFixedAndMovingPyramid( fixedPyramid, movingPyramid );
   if (pyramidReturnValue != CTNonlinearRegistration::OK)
   {
      return pyramidReturnValue;
   }
   
   DeformationFieldType::Pointer deformationField = 0;
   
   // main loop over shrink levels
   for (int shrinkLevel=maxNrShrinkLevels; shrinkLevel>=0; --shrinkLevel)
   {
      std::cout << "shrink level: " << shrinkLevel << std::endl;
            
      // setup the demons algorithm
      typename CurvatureRegistrationFilterType::Pointer 
         curvatureDeformableFilter = CurvatureRegistrationFilterType::New();
         
      curvatureDeformableFilter->SetFixedImage(  fixedPyramid[shrinkLevel] );
      curvatureDeformableFilter->SetMovingImage( movingPyramid[shrinkLevel] );
      curvatureDeformableFilter->SetNumberOfIterations( 
         shrinkLevel * iterationsParam1 + iterationsParam2 );
      curvatureDeformableFilter->SetTimeStep( 1 );
      curvatureDeformableFilter->SetConstraintWeight( 1 );
      curvatureDeformableFilter->SetInPlace( true );
      //curvatureDeformableFilter->Print(std::cout);
      
      // Create the Command observer and register it with the deformable 
      // registration filter.
      typename CurvatureCommandIterationUpdate::Pointer observer = 
         CurvatureCommandIterationUpdate::New();
      curvatureDeformableFilter->AddObserver( itk::IterationEvent(), observer );
      
      // prepare the initial deformation field by taking the result of a 
      // lower resolution as input
      // else case: do nothing, default displacement field is chosen at the 
      // lowest resolution
      if (shrinkLevel < static_cast<int>(maxNrShrinkLevels))
      {
         DeformationFieldType::Pointer tmpDefField = 
            DeformationFieldType::New();
         tmpDefField->SetRegions( 
            fixedPyramid[shrinkLevel]->GetLargestPossibleRegion() );
         tmpDefField->SetSpacing( 
            fixedPyramid[shrinkLevel]->GetSpacing() );
         tmpDefField->SetOrigin( 
            fixedPyramid[shrinkLevel]->GetOrigin() );
         tmpDefField->Allocate();
         
         // take the last deformation field as initial solution of this level
         // needs upsampling
         NonlinearRegistrationUtilities::upsampleDeformationField<
            DeformationFieldType, 
            DeformationFieldType>( tmpDefField, deformationField );
            
         // now we smooth the upsampled deformation field with sigma 1
         NonlinearRegistrationUtilities::
            efficientInplaceDiscreteGaussSmoothDeformationField<
               DeformationFieldType>(tmpDefField, 1, 1, 1);
               
         curvatureDeformableFilter->SetInitialDeformationField( tmpDefField );
      }
      
      EXCEPTION_CHECKED_ITK_FILTER_UPDATE(
         "perform curvature registration update",
         curvatureDeformableFilter->Update() );

      // get output
      deformationField = curvatureDeformableFilter->GetOutput();
   }
   
   // if we are after the last resolution level, warp the moving image 
   // according to the displacement field
   this->warpMovingImage( deformationField );
   
   return CTNonlinearRegistration::OK;
}

#endif // _CT_CURVATURE_REGISTRATION_TXX_