ITK/Examples/Registration/MutualInformation

MutualInformation.cxx

<source lang="cpp">

1. include "itkImageRegistrationMethod.h"
2. include "itkTranslationTransform.h"
3. include "itkMutualInformationImageToImageMetric.h"
4. include "itkGradientDescentOptimizer.h"
5. include "itkNormalizeImageFilter.h"
6. include "itkDiscreteGaussianImageFilter.h"
7. include "itkResampleImageFilter.h"
8. include "itkCastImageFilter.h"
9. include "itkCheckerBoardImageFilter.h"
10. include "itkEllipseSpatialObject.h"
11. include "itkSpatialObjectToImageFilter.h"
12. include "itkImageFileWriter.h"

const unsigned int Dimension = 2; typedef unsigned char PixelType;

typedef itk::Image< PixelType, Dimension > ImageType;

static void CreateEllipseImage(ImageType::Pointer image); static void CreateCircleImage(ImageType::Pointer image);

int main( int argc, char *argv[] ) {

 // Generate synthetic fixed and moving images
 ImageType::Pointer  fixedImage = ImageType::New();
 CreateCircleImage(fixedImage);
 ImageType::Pointer movingImage = ImageType::New();
 CreateEllipseImage(movingImage);

 // Write the two synthetic inputs
 typedef itk::ImageFileWriter< ImageType >  WriterType;

 WriterType::Pointer      fixedWriter =  WriterType::New();
 fixedWriter->SetFileName("fixed.png");
 fixedWriter->SetInput( fixedImage);
 fixedWriter->Update();

 WriterType::Pointer      movingWriter =  WriterType::New();
 movingWriter->SetFileName("moving.png");
 movingWriter->SetInput( movingImage);
 movingWriter->Update();

 // We use floats internally
 typedef   float                                    InternalPixelType;
 typedef itk::Image< float, 2> InternalImageType;

 // Normalize the images
 typedef itk::NormalizeImageFilter<ImageType, InternalImageType> NormalizeFilterType;

 NormalizeFilterType::Pointer fixedNormalizer = NormalizeFilterType::New();
 NormalizeFilterType::Pointer movingNormalizer = NormalizeFilterType::New();

 fixedNormalizer->SetInput(  fixedImage);
 movingNormalizer->SetInput( movingImage);

 // Smooth the normalized images
 typedef itk::DiscreteGaussianImageFilter<InternalImageType,InternalImageType> GaussianFilterType;

 GaussianFilterType::Pointer fixedSmoother  = GaussianFilterType::New();
 GaussianFilterType::Pointer movingSmoother = GaussianFilterType::New();

 fixedSmoother->SetVariance( 2.0 );
 movingSmoother->SetVariance( 2.0 );

 fixedSmoother->SetInput( fixedNormalizer->GetOutput() );
 movingSmoother->SetInput( movingNormalizer->GetOutput() );

 typedef itk::TranslationTransform< double, Dimension > TransformType;
 typedef itk::GradientDescentOptimizer                  OptimizerType;
 typedef itk::LinearInterpolateImageFunction<
                                   InternalImageType,
                                   double             > InterpolatorType;
 typedef itk::ImageRegistrationMethod<
                                   InternalImageType,
                                   InternalImageType >  RegistrationType;
 typedef itk::MutualInformationImageToImageMetric<
                                         InternalImageType,
                                         InternalImageType >    MetricType;

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

 registration->SetOptimizer(     optimizer     );
 registration->SetTransform(     transform     );
 registration->SetInterpolator(  interpolator  );

 MetricType::Pointer         metric        = MetricType::New();
 registration->SetMetric( metric  );

 //  The metric requires a number of parameters to be selected, including
 //  the standard deviation of the Gaussian kernel for the fixed image
 //  density estimate, the standard deviation of the kernel for the moving
 //  image density and the number of samples use to compute the densities
 //  and entropy values. Details on the concepts behind the computation of
 //  the metric can be found in Section
 //  \ref{sec:MutualInformationMetric}.  Experience has
 //  shown that a kernel standard deviation of $0.4$ works well for images
 //  which have been normalized to a mean of zero and unit variance.  We
 //  will follow this empirical rule in this example.

 metric->SetFixedImageStandardDeviation(  0.4 );
 metric->SetMovingImageStandardDeviation( 0.4 );

 registration->SetFixedImage(    fixedSmoother->GetOutput()    );
 registration->SetMovingImage(   movingSmoother->GetOutput()   );

 fixedNormalizer->Update();
 ImageType::RegionType fixedImageRegion =
      fixedNormalizer->GetOutput()->GetBufferedRegion();
 registration->SetFixedImageRegion( fixedImageRegion );

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

 initialParameters[0] = 0.0;  // Initial offset along X
 initialParameters[1] = 0.0;  // Initial offset along Y

 registration->SetInitialTransformParameters( initialParameters );

 //  Software Guide : BeginLatex
 //
 //  We should now define the number of spatial samples to be considered in
 //  the metric computation. Note that we were forced to postpone this setting
 //  until we had done the preprocessing of the images because the number of
 //  samples is usually defined as a fraction of the total number of pixels in
 //  the fixed image.
 //
 //  The number of spatial samples can usually be as low as $1\%$ of the total
 //  number of pixels in the fixed image. Increasing the number of samples
 //  improves the smoothness of the metric from one iteration to another and
 //  therefore helps when this metric is used in conjunction with optimizers
 //  that rely of the continuity of the metric values. The trade-off, of
 //  course, is that a larger number of samples result in longer computation
 //  times per every evaluation of the metric.
 //
 //  It has been demonstrated empirically that the number of samples is not a
 //  critical parameter for the registration process. When you start fine
 //  tuning your own registration process, you should start using high values
 //  of number of samples, for example in the range of $20\%$ to $50\%$ of the
 //  number of pixels in the fixed image. Once you have succeeded to register
 //  your images you can then reduce the number of samples progressively until
 //  you find a good compromise on the time it takes to compute one evaluation
 //  of the Metric. Note that it is not useful to have very fast evaluations
 //  of the Metric if the noise in their values results in more iterations
 //  being required by the optimizer to converge. You must then study the
 //  behavior of the metric values as the iterations progress.
 
 const unsigned int numberOfPixels = fixedImageRegion.GetNumberOfPixels();

 const unsigned int numberOfSamples =
                       static_cast< unsigned int >( numberOfPixels * 0.01 );

 metric->SetNumberOfSpatialSamples( numberOfSamples );

 //  Since larger values of mutual information indicate better matches than
 //  smaller values, we need to maximize the cost function in this example.
 //  By default the GradientDescentOptimizer class is set to minimize the
 //  value of the cost-function. It is therefore necessary to modify its
 //  default behavior by invoking the \code{MaximizeOn()} method.
 //  Additionally, we need to define the optimizer's step size using the
 //  \code{SetLearningRate()} method.

 optimizer->SetLearningRate( 15.0 );
 optimizer->SetNumberOfIterations( 200 );
 optimizer->MaximizeOn(); // We want to maximize mutual information (the default of the optimizer is to minimize)

 // Note that large values of the learning rate will make the optimizer
 // unstable. Small values, on the other hand, may result in the optimizer
 // needing too many iterations in order to walk to the extrema of the cost
 // function. The easy way of fine tuning this parameter is to start with
 // small values, probably in the range of $\{5.0,10.0\}$. Once the other
 // registration parameters have been tuned for producing convergence, you
 // may want to revisit the learning rate and start increasing its value until
 // you observe that the optimization becomes unstable.  The ideal value for
 // this parameter is the one that results in a minimum number of iterations
 // while still keeping a stable path on the parametric space of the
 // optimization. Keep in mind that this parameter is a multiplicative factor
 // applied on the gradient of the Metric. Therefore, its effect on the
 // optimizer step length is proportional to the Metric values themselves.
 // Metrics with large values will require you to use smaller values for the
 // learning rate in order to maintain a similar optimizer behavior.

 try
   {
   registration->Update();
   std::cout << "Optimizer stop condition: "
             << registration->GetOptimizer()->GetStopConditionDescription()
             << std::endl;
   }
 catch( itk::ExceptionObject & err )
   {
   std::cout << "ExceptionObject caught !" << std::endl;
   std::cout << err << std::endl;
   return EXIT_FAILURE;
   }

 ParametersType finalParameters = registration->GetLastTransformParameters();

 double TranslationAlongX = finalParameters[0];
 double TranslationAlongY = finalParameters[1];

 unsigned int numberOfIterations = optimizer->GetCurrentIteration();

 double bestValue = optimizer->GetValue();

 // Print out results
 std::cout << std::endl;
 std::cout << "Result = " << 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;
 std::cout << " Numb. Samples = " << numberOfSamples    << std::endl;

 typedef itk::ResampleImageFilter<
                           ImageType,
                           ImageType >    ResampleFilterType;

 TransformType::Pointer finalTransform = TransformType::New();

 finalTransform->SetParameters( finalParameters );
 finalTransform->SetFixedParameters( transform->GetFixedParameters() );

 ResampleFilterType::Pointer resample = ResampleFilterType::New();

 resample->SetTransform( finalTransform );
 resample->SetInput( movingImage);

 resample->SetSize(    fixedImage->GetLargestPossibleRegion().GetSize() );
 resample->SetOutputOrigin(  fixedImage->GetOrigin() );
 resample->SetOutputSpacing( fixedImage->GetSpacing() );
 resample->SetOutputDirection( fixedImage->GetDirection() );
 resample->SetDefaultPixelValue( 100 );

 WriterType::Pointer      writer =  WriterType::New();
 writer->SetFileName("output.png");
 writer->SetInput( resample->GetOutput()   );
 writer->Update();

 return EXIT_SUCCESS;

}

void CreateEllipseImage(ImageType::Pointer image) {

 typedef itk::EllipseSpatialObject< Dimension >   EllipseType;

 typedef itk::SpatialObjectToImageFilter<
   EllipseType, ImageType >   SpatialObjectToImageFilterType;

 SpatialObjectToImageFilterType::Pointer imageFilter =
   SpatialObjectToImageFilterType::New();

 ImageType::SizeType size;
 size[ 0 ] =  100;
 size[ 1 ] =  100;

 imageFilter->SetSize( size );

 ImageType::SpacingType spacing;
 spacing.Fill(1);
 imageFilter->SetSpacing(spacing);

 EllipseType::Pointer ellipse    = EllipseType::New();
 EllipseType::ArrayType radiusArray;
 radiusArray[0] = 10;
 radiusArray[1] = 20;
 ellipse->SetRadius(radiusArray);

 typedef EllipseType::TransformType                 TransformType;
 TransformType::Pointer transform = TransformType::New();
 transform->SetIdentity();

 TransformType::OutputVectorType  translation;
 TransformType::CenterType        center;

 translation[ 0 ] =  65;
 translation[ 1 ] =  45;
 transform->Translate( translation, false );

 ellipse->SetObjectToParentTransform( transform );

 imageFilter->SetInput(ellipse);

 ellipse->SetDefaultInsideValue(255);
 ellipse->SetDefaultOutsideValue(0);
 imageFilter->SetUseObjectValue( true );
 imageFilter->SetOutsideValue( 0 );

 imageFilter->Update();

 image->Graft(imageFilter->GetOutput());

}

void CreateCircleImage(ImageType::Pointer image) {

typedef itk::EllipseSpatialObject< Dimension >   EllipseType;

 typedef itk::SpatialObjectToImageFilter<
   EllipseType, ImageType >   SpatialObjectToImageFilterType;

 SpatialObjectToImageFilterType::Pointer imageFilter =
   SpatialObjectToImageFilterType::New();

 ImageType::SizeType size;
 size[ 0 ] =  100;
 size[ 1 ] =  100;

 imageFilter->SetSize( size );

 ImageType::SpacingType spacing;
 spacing.Fill(1);
 imageFilter->SetSpacing(spacing);

 EllipseType::Pointer ellipse    = EllipseType::New();
 EllipseType::ArrayType radiusArray;
 radiusArray[0] = 10;
 radiusArray[1] = 10;
 ellipse->SetRadius(radiusArray);

 typedef EllipseType::TransformType                 TransformType;
 TransformType::Pointer transform = TransformType::New();
 transform->SetIdentity();

 TransformType::OutputVectorType  translation;
 TransformType::CenterType        center;

 translation[ 0 ] =  50;
 translation[ 1 ] =  50;
 transform->Translate( translation, false );

 ellipse->SetObjectToParentTransform( transform );

 imageFilter->SetInput(ellipse);

 ellipse->SetDefaultInsideValue(255);
 ellipse->SetDefaultOutsideValue(0);
 imageFilter->SetUseObjectValue( true );
 imageFilter->SetOutsideValue( 0 );

 imageFilter->Update();

 image->Graft(imageFilter->GetOutput());

} </source>

CMakeLists.txt

<source lang="cmake"> cmake_minimum_required(VERSION 2.6)

PROJECT(MutualInformation)

FIND_PACKAGE(ITK REQUIRED) INCLUDE(${ITK_USE_FILE})

ADD_EXECUTABLE(MutualInformation MutualInformation.cxx) TARGET_LINK_LIBRARIES(MutualInformation ITKIO ITKNumerics)

</source>