/*
* 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: CTStandardDemonsRegistration.txx,v $
* Language    : C++
* Description : 
*
* Author     : Martin Urschler
* EMail      : urschler@icg.tu-graz.ac.at
* Date       : $Date: 2007-07-04 14:18:03 $
* Version    : $Revision: 1.13 $
* Full Id    : $Id: CTStandardDemonsRegistration.txx,v 1.13 2007-07-04 14:18:03 urschler Exp $
*
*/


#ifndef _CT_STANDARD_DEMONS_REGISTRATION_TXX_
#define _CT_STANDARD_DEMONS_REGISTRATION_TXX_


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

#include "itkDemonsRegistrationFilter.h"

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

template< typename TPixelDataType >
void 
CTStandardDemonsRegistration<TPixelDataType>::
DemonsCommandIterationUpdate::
Execute(const itk::Object * object, const itk::EventObject & event)
{
   typedef itk::DemonsRegistrationFilter< 
         InputImageType, 
         InputImageType,
         DeformationFieldType > 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 
CTStandardDemonsRegistration<TPixelDataType>::
matching( const unsigned int maxNrShrinkLevels, 
   const unsigned int iterationsParam1, 
   const unsigned int iterationsParam2 )
{  
   typedef itk::DemonsRegistrationFilter< InputImageType, 
      InputImageType, DeformationFieldType > 
         DemonsRegistrationFilterType;
      
   // 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 DemonsRegistrationFilterType::Pointer demonsDeformableFilter = 
         DemonsRegistrationFilterType::New();
         
      demonsDeformableFilter->SetFixedImage(  fixedPyramid[shrinkLevel] );
      demonsDeformableFilter->SetMovingImage( movingPyramid[shrinkLevel] );
      demonsDeformableFilter->SetNumberOfIterations( 
         shrinkLevel * iterationsParam1 + iterationsParam2 );
      demonsDeformableFilter->SetStandardDeviations( 1.0 );
      demonsDeformableFilter->SetInPlace( true );
      demonsDeformableFilter->SetUseMovingImageGradient(false);
      demonsDeformableFilter->SetIntensityDifferenceThreshold( 0.1 );
      //demonsDeformableFilter->Print(std::cout);
      
      // Create the Command observer and register it with the deformable 
      // registration filter.
      typename DemonsCommandIterationUpdate::Pointer observer = 
         DemonsCommandIterationUpdate::New();
      demonsDeformableFilter->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();
         
         tmpDefField = DownAndUpsampleUtilities::
            doubleVectorFieldByLinearInterpolation<DeformationFieldType>(
               deformationField,
               this->m_movingDownsampleModeSequence[shrinkLevel]);
         
         // 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);
      
         demonsDeformableFilter->SetInitialDeformationField( tmpDefField );
      }
      
      EXCEPTION_CHECKED_ITK_FILTER_UPDATE("perform demons update",
         demonsDeformableFilter->Update() );

      // get output
      deformationField = demonsDeformableFilter->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_STANDARD_DEMONS_REGISTRATION_TXX_