[Insight-users] Problem of using Euler2D transform in Example MultiResImageRegistration1

[Junyi Xia]

Hello,
        I am trying to using Euler2D transform in 
MultiResImageRegistration1.cxx. I replace the TranslationTransform with 
Euler2DTransform. While I registered with example images( 
BrainProtonDensitySliceBorder20.png and 
BrainProtonDensitySliceShifted13x17y.png), the output had the following 
errors:
//////////////////////////////////////////////////
ExceptionObject caught!
itk::ExecptionObject(0xed7150)
File:/xxx/..../itkMattesMutualInformationImageToImageMetric.txx
Line: 985
Description: itk::ERROR: MattesMutualInformationImageToImageMetric: Too 
many samples map outside moving image buffer: 4/10000.
/////////////////////////////////////////////////
The following error would pop up while I registered 
BrainProtonDensitySliceBorder20.png to itself.

///////////////////////////////////////////////
ExceptionObject caught!
itk::ExecptionObject
File:/xxx/..../itkMattesMutualInformationImageToImageMetric.txx
Line: 985
Description: itk::ERROR: MattesMutualInformationImageToImageMetric: Too 
many samples map outside moving image buffer: 0/10000.
////////////////////////////////////////////////

         It seems like it is not the initializtion problem since it 
failed when registered with itself. Any suggestion? Thank.

Junyi

//// Code

#include "itkMultiResolutionImageRegistrationMethod.h"

#include "itkEuler2DTransform.h"

#include "itkMattesMutualInformationImageToImageMetric.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkRegularStepGradientDescentOptimizer.h"
#include "itkMultiResolutionPyramidImageFilter.h"
#include "itkImage.h"
//#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"

#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkCheckerBoardImageFilter.h"

#include "itkTimeProbesCollectorBase.h"

#include "itkCommand.h"

*template* <*typename* TRegistration>
*class* RegistrationInterfaceCommand : *public* itk::Command 
{

*public*:
  *typedef*  RegistrationInterfaceCommand   Self;
  *typedef*  itk::Command                   Superclass;
  *typedef*  itk::SmartPointer<Self>        Pointer;
  itkNewMacro( Self );
*protected*:
  RegistrationInterfaceCommand() {};

*public*:
  *typedef*   TRegistration                              RegistrationType;
  *typedef*   RegistrationType *                         RegistrationPointer;
  *typedef*   itk::RegularStepGradientDescentOptimizer   OptimizerType;
//
  *typedef*   OptimizerType *                            OptimizerPointer;

  void Execute(itk::Object * object, const itk::EventObject & event)
  {

    *if*( !(itk::IterationEvent().CheckEvent( &event )) )
      {
      *return*;
      }

    RegistrationPointer registration =
                            *dynamic_cast*<RegistrationPointer>( object );

    OptimizerPointer optimizer = *dynamic_cast*< OptimizerPointer >( 
                       registration->GetOptimizer() );

    *if* ( registration->GetCurrentLevel() == 0 )
      {
      optimizer->SetMaximumStepLength( 16.00 );  
      optimizer->SetMinimumStepLength( 2.5 );
//
      }
    *else*
      {
       optimizer->SetMaximumStepLength( 
         optimizer->GetCurrentStepLength() );
      optimizer->SetMinimumStepLength(
                optimizer->GetMinimumStepLength() / 10.0 );

//
      }
  }

  void Execute(const itk::Object * , const itk::EventObject & )
    { *return*; }
};

*class* CommandIterationUpdate : *public* itk::Command 
{
*public*:
  *typedef*  CommandIterationUpdate   Self;
  *typedef*  itk::Command             Superclass;
  *typedef*  itk::SmartPointer<Self>  Pointer;
  itkNewMacro( Self );
*protected*:
  CommandIterationUpdate() {};
*public*:
  *typedef*   itk::RegularStepGradientDescentOptimizer     OptimizerType;
//
  *typedef*   const OptimizerType   *           OptimizerPointer;

  void Execute(itk::Object *caller, const itk::EventObject & event)
    {
      Execute( (const itk::Object *)caller, event);
    }

  void Execute(const itk::Object * object, const itk::EventObject & event)
    {
      OptimizerPointer optimizer = 
        *dynamic_cast*< OptimizerPointer >( object );
      *if*( !(itk::IterationEvent().CheckEvent( &event )) )
        {
        *return*;
        }
      std::cout << optimizer->GetCurrentIteration() << "   ";
      std::cout << optimizer->GetValue() << "   ";
      std::cout << optimizer->GetCurrentPosition() << std::endl;
    }
};


int main( int argc, char *argv[] )
{
  *if*( argc < 3 )
    {
    std::cerr << "Missing Parameters " << std::endl;
    std::cerr << "Usage: " << argv[0];
    std::cerr << " fixedImageFile  movingImageFile ";
    std::cerr << "outputImagefile [checkerBoardBefore] [checkerBoardAfter]" 
      << std::endl;
    *return* 1;
    }
  
  const    unsigned int    Dimension = 2;
  *typedef*  unsigned short  PixelType;
  
  *typedef* itk::Image< PixelType, Dimension >  FixedImageType;
  *typedef* itk::Image< PixelType, Dimension >  MovingImageType;


  *typedef*   float     InternalPixelType;
  *typedef* itk::Image< InternalPixelType, Dimension > InternalImageType;


  ///typedef itk::TranslationTransform< double, Dimension > TransformType;/
  *typedef* itk::RegularStepGradientDescentOptimizer       OptimizerType;
  *typedef* itk::Euler2DTransform< double > TransformType;
	//
  *typedef* itk::LinearInterpolateImageFunction< 
                                    InternalImageType,
                                    double             > InterpolatorType;
  *typedef* itk::MattesMutualInformationImageToImageMetric< 
                                    InternalImageType, 
                                    InternalImageType >   MetricType;
  *typedef* itk::MultiResolutionImageRegistrationMethod< 
                                    InternalImageType, 
                                    InternalImageType >   RegistrationType;

  *typedef* itk::MultiResolutionPyramidImageFilter<
                                    InternalImageType,
                                    InternalImageType >   FixedImagePyramidType;
  *typedef* itk::MultiResolutionPyramidImageFilter<
                                    InternalImageType,
                                    InternalImageType >   MovingImagePyramidType;

  TransformType::Pointer      transform     = TransformType::New();
  OptimizerType::Pointer      optimizer     = OptimizerType::New();
  InterpolatorType::Pointer   interpolator  = InterpolatorType::New();
  RegistrationType::Pointer   registration  = RegistrationType::New();
  MetricType::Pointer         metric        = MetricType::New();

  FixedImagePyramidType::Pointer fixedImagePyramid = 
      FixedImagePyramidType::New();
  MovingImagePyramidType::Pointer movingImagePyramid =
      MovingImagePyramidType::New();

  registration->SetOptimizer(     optimizer     );
  registration->SetTransform(     transform     );
  registration->SetInterpolator(  interpolator  );
  registration->SetMetric( metric  );
  registration->SetFixedImagePyramid( fixedImagePyramid );
  registration->SetMovingImagePyramid( movingImagePyramid );
  

  *typedef* itk::ImageFileReader< FixedImageType  > FixedImageReaderType;
  *typedef* itk::ImageFileReader< MovingImageType > MovingImageReaderType;

  FixedImageReaderType::Pointer  fixedImageReader  = FixedImageReaderType::New();
  MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();

  fixedImageReader->SetFileName(  argv[1] );
  movingImageReader->SetFileName( argv[2] );


 
  *typedef* itk::CastImageFilter< 
                        FixedImageType, InternalImageType > FixedCastFilterType;
  *typedef* itk::CastImageFilter< 
                        MovingImageType, InternalImageType > MovingCastFilterType;

  FixedCastFilterType::Pointer fixedCaster   = FixedCastFilterType::New();
  MovingCastFilterType::Pointer movingCaster = MovingCastFilterType::New();
 
  fixedCaster->SetInput(  fixedImageReader->GetOutput() );
  movingCaster->SetInput( movingImageReader->GetOutput() );


  registration->SetFixedImage(    fixedCaster->GetOutput()    );
  registration->SetMovingImage(   movingCaster->GetOutput()   );
 

  fixedCaster->Update();

  registration->SetFixedImageRegion( 
       fixedCaster->GetOutput()->GetBufferedRegion() );
   

  *typedef* RegistrationType::ParametersType ParametersType;
  ParametersType initialParameters( transform->GetNumberOfParameters() );

  initialParameters[0] = 0.0;  
  initialParameters[1] = 0.0;  
  initialParameters[2] = 0.0;


   registration->SetInitialTransformParameters( initialParameters );
	*typedef* OptimizerType::ScalesType OptimizerScalesType;

	OptimizerScalesType optimizerScales( transform->GetNumberOfParameters() );
	const double translationScale = 1.0 / 1000.0;

	optimizerScales[0] = 1.0;
	optimizerScales[1] = translationScale;
	optimizerScales[2] = translationScale;

	optimizer->SetScales( optimizerScales );

  registration->SetInitialTransformParameters( initialParameters);

  metric->SetNumberOfHistogramBins( 60 );
  metric->SetNumberOfSpatialSamples( 10000 );
//

  optimizer->SetNumberOfIterations( 200 );
	


  /// Create the Command observer and register it with the optimizer./
  ////
  CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
  optimizer->AddObserver( itk::IterationEvent(), observer );



  *typedef* RegistrationInterfaceCommand<RegistrationType> CommandType;
  CommandType::Pointer command = CommandType::New();
  registration->AddObserver( itk::IterationEvent(), command );
 
  registration->SetNumberOfLevels( 3 );
	
	itk::TimeProbesCollectorBase timer; 
  *try* 
    { 
	timer.Start("Registration");
    registration->StartRegistration(); 
	timer.Stop("Registration");
    } 
  *catch*( itk::ExceptionObject & err ) 
    { 
    std::cout << "ExceptionObject caught !" << std::endl; 
    std::cout << err << std::endl; 
    *return* -1;
    } 
  timer.Report();

  ParametersType finalParameters = registration->GetLastTransformParameters();
  
  /double rotationAngle = finalParameters[0];/
  /double TranslationAlongX = finalParameters[1];//  
  double TranslationAlongY = finalParameters[2];/

  
  unsigned int numberOfIterations = optimizer->GetCurrentIteration();
  
  double bestValue = optimizer->GetValue();


  /// Print out results/
  std::cout << "Result = " << std::endl;
  /std::cout << " Rotation = " << rotationAngle  << std::endl;/
  std::cout << " Translation X = " << TranslationAlongX  << std::endl;
  std::cout << " Translation Y = " << TranslationAlongY  << std::endl;
  std::cout << " Iterations    = " << numberOfIterations << std::endl;
  std::cout << " Metric value  = " << bestValue          << std::endl;

	return 0;
}