[Insight-users] NormalizedMutualInformationMetric

[Giorgos Pichis]

Hi there,
I am trying to register two 2-D images,
using the NormalizedMutualInformationMetric (NMI), but I have terrible
results.
I have modified the ImageRegistration7.cxx example, in order to use this
metric.
My questions are:
1.Does NMI Metric works only with  itkOnePlusOneEvolutionaryOptimizer as
shown in example 14?
2.Beside the number of histogram bins, is there any other parameter of the
NMI,that could be modified,
in order to get better results?
3.Does the images need to be preprocessed (normalizedImageFilter , and
Smoothed (GaussianImageFilter)),
before entering the registration process?
Thanks in advance
Giorgos


/*========================================================================*/

#include "itkImageRegistrationMethod.h"
#include "itkNormalizedMutualInformationHistogramImageToImageMetric.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkRegularStepGradientDescentOptimizer.h"
#include "itkImage.h"

#include "itkCenteredTransformInitializer.h"
#include "itkCenteredSimilarity2DTransform.h"

#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"

#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkSubtractImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkIdentityTransform.h"


#include "itkMeanSquaresImageToImageMetric.h"
#include "itkCorrelationCoefficientHistogramImageToImageMetric.h"


#include "itkTranslationTransform.h"
#include "itkNearestNeighborInterpolateImageFunction.h"
#include "itkCheckerBoardImageFilter.h"

//
#include "itkCommand.h"
class CommandIterationUpdate : public itk::Command
{
public:
  typedef  CommandIterationUpdate   Self;
  typedef  itk::Command             Superclass;
  typedef itk::SmartPointer<Self>  Pointer;
  itkNewMacro( Self );
protected:
  CommandIterationUpdate() {m_LastMetricValue = 0;};
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;
        }
        double currentValue = optimizer->GetValue();
        if( fabs( m_LastMetricValue - currentValue ) > 1e-4 )
        {
        std::cout << optimizer->GetCurrentIteration() << "   ";
        std::cout << optimizer->GetValue() << "   ";
        std::cout << optimizer->GetCurrentPosition() << std::endl;
        m_LastMetricValue = currentValue;
        }
    }
    private:
  double m_LastMetricValue;
};


int main( int argc, char *argv[] )
{
  if( argc < 4 )
    {
    std::cerr << "Missing Parameters " << std::endl;
    std::cerr << "Usage: " << argv[0];
    std::cerr << " fixedImageFile  movingImageFile ";
    std::cerr << " outputImagefile  [differenceBeforeRegistration] ";
    std::cerr << " [differenceAfterRegistration] [Checkerboard Image]";

    std::cerr << " [initialScaling] [initialAngle] ";
    std::cerr << " [steplength] ";
    std::cerr << std::endl;
    return EXIT_FAILURE;
    }

  const    unsigned int    Dimension = 2;
  typedef  float   PixelType;
  typedef  float           InternalPixelType;


  typedef itk::Image< PixelType, Dimension >  FixedImageType;
  typedef itk::Image< PixelType, Dimension >  MovingImageType;

/// reading the images

  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::CenteredSimilarity2DTransform< double > TransformType;

  typedef itk::RegularStepGradientDescentOptimizer       OptimizerType;
  typedef itk::NormalizedMutualInformationHistogramImageToImageMetric<
InternalImageType, InternalImageType>
    MetricType;
  typedef itk:: LinearInterpolateImageFunction< InternalImageType, double >
    InterpolatorType;
  typedef itk::ImageRegistrationMethod< InternalImageType, InternalImageType
>
    RegistrationType;

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

  registration->SetMetric(        metric        );
  registration->SetOptimizer(     optimizer     );
  registration->SetInterpolator(  interpolator  );

  TransformType::Pointer  transform = TransformType::New();
  registration->SetTransform( transform );



  /////NORMALIZED MI SETTINGS
   unsigned int numberOfHistogramBins;
  std::cout<< " Enter Number of Histogram Bins : " ;
  std::cin >> numberOfHistogramBins;
  std::cout<<std::endl;

    /*if( argc > 4 )
    {
    numberOfHistogramBins = atoi( argv[4] );
    std::cout << "Using " << numberOfHistogramBins << " Histogram bins" <<
std::endl;
    }*/
    MetricType::HistogramType::SizeType histogramSize;
  histogramSize[0] = numberOfHistogramBins;
  histogramSize[1] = numberOfHistogramBins;
  metric->SetHistogramSize( histogramSize );


  const unsigned int numberOfParameters =
transform->GetNumberOfParameters();

  typedef MetricType::ScalesType ScalesType;
  ScalesType scales( numberOfParameters );

  scales.Fill( 1.0 );

  metric->SetDerivativeStepLengthScales(scales);

//////
///// SETTING THE IMAGES FOR THE REGISTRATION PROCESS

  registration->SetFixedImage(    fixedImageReader->GetOutput()    );
  registration->SetMovingImage(   movingImageReader->GetOutput()   );

  fixedImageReader->Update();
  FixedImageType::RegionType fixedImageRegion =
       fixedImageReader->GetOutput()->GetBufferedRegion();
  registration->SetFixedImageRegion( fixedImageRegion );


  typedef itk::CenteredTransformInitializer<
                                    TransformType,
                               FixedImageType,
                                    MovingImageType >
TransformInitializerType;

  TransformInitializerType::Pointer initializer =
TransformInitializerType::New();

  initializer->SetTransform(   transform );

  initializer->SetFixedImage(  fixedImageReader->GetOutput() );
  initializer->SetMovingImage( movingImageReader->GetOutput() );

  initializer->GeometryOn();

  initializer->InitializeTransform();

// REMAINING TRANSFORM PARAMETERS
  double initialScale = 1.0;

  if( argc > 7 )
    {
    initialScale =  atof( argv[7] );
    }

  double initialAngle = 0.0;

  if( argc > 8 )
    {
    initialAngle =  atof( argv[8] );
    }

  // Software Guide : BeginCodeSnippet
  transform->SetScale( initialScale );
  transform->SetAngle( initialAngle );
  // Software Guide : EndCodeSnippet

  registration->SetInitialTransformParameters( transform->GetParameters() );

  typedef OptimizerType::ScalesType       OptimizerScalesType;
  OptimizerScalesType optimizerScales( transform->GetNumberOfParameters() );
  const double translationScale = 1.0 / 100.0;

  optimizerScales[0] = 10.0;
  optimizerScales[1] =  1.0;
  optimizerScales[2] =  translationScale;
  optimizerScales[3] =  translationScale;
  optimizerScales[4] =  translationScale;
  optimizerScales[5] =  translationScale;

  optimizer->SetScales( optimizerScales );


  double steplength = 1.0;

  if( argc > 9 )
    {
    steplength = atof( argv[9] );
    }

  // Software Guide : BeginCodeSnippet
  optimizer->SetMaximumStepLength( steplength );
  optimizer->SetMinimumStepLength( 0.0001 );
  optimizer->SetNumberOfIterations( 500 );
  // Software Guide : EndCodeSnippet


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


  try
    {
    registration->StartRegistration();
    }
  catch( itk::ExceptionObject & err )
    {
    std::cerr << "ExceptionObject caught !" << std::endl;
    std::cerr << err << std::endl;
    return EXIT_FAILURE;
    }

  OptimizerType::ParametersType finalParameters =
                    registration->GetLastTransformParameters();


  const double finalScale           = finalParameters[0];
  const double finalAngle           = finalParameters[1];
  const double finalRotationCenterX = finalParameters[2];
  const double finalRotationCenterY = finalParameters[3];
  const double finalTranslationX    = finalParameters[4];
  const double finalTranslationY    = finalParameters[5];

  const unsigned int numberOfIterations = optimizer->GetCurrentIteration();

  const double bestValue = optimizer->GetValue();


  // Print out results
  //
  const double finalAngleInDegrees = finalAngle * 45.0 / atan(1.0);

  std::cout << std::endl;
  std::cout << "Result = " << std::endl;
  std::cout << " Scale         = " << finalScale  << std::endl;
  std::cout << " Angle (radians) " << finalAngle  << std::endl;
  std::cout << " Angle (degrees) " << finalAngleInDegrees  << std::endl;
  std::cout << " Center X      = " << finalRotationCenterX  << std::endl;
  std::cout << " Center Y      = " << finalRotationCenterY  << std::endl;
  std::cout << " Translation X = " << finalTranslationX  << std::endl;
  std::cout << " Translation Y = " << finalTranslationY  << std::endl;
  std::cout << " Iterations    = " << numberOfIterations << std::endl;
  std::cout << " Metric value  = " << bestValue          << std::endl;


 ...........................
............................

  return EXIT_SUCCESS;
}
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20080508/277b49c4/attachment-0001.htm>