/*
* 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: CTLevelSetMotionRegistration.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.5 $
* Full Id    : $Id: CTLevelSetMotionRegistration.txx,v 1.5 2006-10-31 09:00:24 urschler Exp $
*
*/


#ifndef _CT_LEVEL_SET_MOTION_REGISTRATION_TXX_
#define _CT_LEVEL_SET_MOTION_REGISTRATION_TXX_


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

#include "itkLevelSetMotionRegistrationFilter.h"


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

template< typename TPixelDataType >
void 
CTLevelSetMotionRegistration<TPixelDataType>::
LevelSetMotionCommandIterationUpdate::
Execute(const itk::Object * object, const itk::EventObject & event)
{
   typedef itk::LevelSetMotionRegistrationFilter< 
         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 
CTLevelSetMotionRegistration<TPixelDataType>::
matching( const unsigned int maxNrShrinkLevels, 
   const unsigned int iterationsParam1, 
   const unsigned int iterationsParam2 )
{  
   typedef itk::LevelSetMotionRegistrationFilter< 
      InputImageType, 
      InputImageType,
      DeformationFieldType > DeformableRegistrationFilterType;
      
   // here the multi resolution level set motion 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;
      
      // calculate intensity range of images
      const float minFixed = SmallUtilityMethods<InputImageType>::
         getMinimumVoxelValue(fixedPyramid[shrinkLevel]);
      const float maxFixed = SmallUtilityMethods<InputImageType>::
         getMaximumVoxelValue(fixedPyramid[shrinkLevel]);
      const float minMoving = SmallUtilityMethods<InputImageType>::
         getMinimumVoxelValue(movingPyramid[shrinkLevel]);
      const float maxMoving = SmallUtilityMethods<InputImageType>::
         getMaximumVoxelValue(movingPyramid[shrinkLevel]);
      
      const typename InputImageType::SpacingType movingSpacing = 
         movingPyramid[shrinkLevel]->GetSpacing();
      
      const float movingImageGradientSmoothingSigma = 2.0f *  
         vnl_math_min( vnl_math_min(movingSpacing[0], movingSpacing[1]), 
                        movingSpacing[2] );
      
      std::cout << "movingImageGradientSmoothingSigma: " << 
         movingImageGradientSmoothingSigma << std::endl;
      
      const float gradientMagnitudeEpsilon = 0.04f * 0.01f *
         ( vnl_math_max( maxFixed, maxMoving ) - 
           vnl_math_min( minFixed, minMoving ) ); 
      
      std::cout << "gradientMagnitudeEpsilon: " << gradientMagnitudeEpsilon << 
         std::endl;
      
      // setup the demons algorithm
      typename DeformableRegistrationFilterType::Pointer deformableFilter = 
         DeformableRegistrationFilterType::New();
         
      deformableFilter->SetFixedImage(  fixedPyramid[shrinkLevel] );
      deformableFilter->SetMovingImage( movingPyramid[shrinkLevel] );
      deformableFilter->SetNumberOfIterations( 
         shrinkLevel * iterationsParam1 + iterationsParam2 );
      deformableFilter->SetInPlace( true );
      deformableFilter->SetAlpha( gradientMagnitudeEpsilon );
      deformableFilter->SetIntensityDifferenceThreshold( 0.1 );
      deformableFilter->SetGradientMagnitudeThreshold( 1e-9 );
      deformableFilter->SetGradientSmoothingStandardDeviations(
         movingImageGradientSmoothingSigma );
      
      //deformableFilter->Print(std::cout);
      
      // Create the Command observer and register it with the deformable 
      // registration filter.
      typename LevelSetMotionCommandIterationUpdate::Pointer observer = 
         LevelSetMotionCommandIterationUpdate::New();
      deformableFilter->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);
               
         deformableFilter->SetInitialDeformationField( tmpDefField );
      }
      
      EXCEPTION_CHECKED_ITK_FILTER_UPDATE("perform level set motion update",
         deformableFilter->Update() );

      // get output
      deformationField = deformableFilter->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_LEVEL_SET_MOTION_REGISTRATION_TXX_