ITK/Examples/Registration/ImageRegistrationMethodAffine

ImageRegistrationMethodAffine.cxx

<source lang="cpp">

1. include "itkCastImageFilter.h"
2. include "itkEllipseSpatialObject.h"
3. include "itkImage.h"
4. include "itkImageRegistrationMethod.h"
5. include "itkLinearInterpolateImageFunction.h"
6. include "itkImageFileReader.h"
7. include "itkImageFileWriter.h"
8. include "itkMeanSquaresImageToImageMetric.h"
9. include "itkRegularStepGradientDescentOptimizer.h"
10. include "itkResampleImageFilter.h"
11. include "itkRescaleIntensityImageFilter.h"
12. include "itkSpatialObjectToImageFilter.h"
13. include "itkAffineTransform.h"

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

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

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

int main(int, char *[] ) {

 //  The transform that will map the fixed image into the moving image.
 typedef itk::AffineTransform< double, Dimension > TransformType;

 //  An optimizer is required to explore the parameter space of the transform
 //  in search of optimal values of the metric.
 typedef itk::RegularStepGradientDescentOptimizer       OptimizerType;

 //  The metric will compare how well the two images match each other. Metric
 //  types are usually parameterized by the image types as it can be seen in
 //  the following type declaration.
 typedef itk::MeanSquaresImageToImageMetric<
     ImageType,
     ImageType >    MetricType;

 //  Finally, the type of the interpolator is declared. The interpolator will
 //  evaluate the intensities of the moving image at non-grid positions.
 typedef itk:: LinearInterpolateImageFunction<
     ImageType,
     double          >    InterpolatorType;

 //  The registration method type is instantiated using the types of the
 //  fixed and moving images. This class is responsible for interconnecting
 //  all the components that we have described so far.
 typedef itk::ImageRegistrationMethod<
     ImageType,
     ImageType >    RegistrationType;

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

 // Each component is now connected to the instance of the registration method.
 registration->SetMetric(        metric        );
 registration->SetOptimizer(     optimizer     );
 registration->SetTransform(     transform     );
 registration->SetInterpolator(  interpolator  );

 // Get the two images
 ImageType::Pointer  fixedImage  = ImageType::New();
 ImageType::Pointer movingImage = ImageType::New();

 CreateSphereImage(fixedImage);
 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();

 // Set the registration inputs
 registration->SetFixedImage(fixedImage);
 registration->SetMovingImage(movingImage);

 registration->SetFixedImageRegion(
   fixedImage->GetLargestPossibleRegion() );

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

 // rotation matrix
 initialParameters[0] = 1.0;  // R(0,0)
 initialParameters[1] = 0.0;  // R(0,1)
 initialParameters[2] = 0.0;  // R(1,0)
 initialParameters[3] = 1.0;  // R(1,1)

 // translation vector
 initialParameters[4] = 0.0;
 initialParameters[5] = 0.0;
 
 registration->SetInitialTransformParameters( initialParameters );

 optimizer->SetMaximumStepLength( .1 ); // If this is set too high, you will get a
 //"itk::ERROR: MeanSquaresImageToImageMetric(0xa27ce70): Too many samples map outside moving image buffer: 1818 / 10000" error

 optimizer->SetMinimumStepLength( 0.01 );

 // Set a stopping criterion
 optimizer->SetNumberOfIterations( 200 );

 // Connect an observer
 //CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
 //optimizer->AddObserver( itk::IterationEvent(), observer );

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

 //  The result of the registration process is an array of parameters that
 //  defines the spatial transformation in an unique way. This final result is
 //  obtained using the \code{GetLastTransformParameters()} method.

 ParametersType finalParameters = registration->GetLastTransformParameters();
 std::cout << "Final parameters: " << finalParameters << std::endl;
 
 //  The value of the image metric corresponding to the last set of parameters
 //  can be obtained with the \code{GetValue()} method of the optimizer.

 const double bestValue = optimizer->GetValue();

 // Print out results
 //
 std::cout << "Result = " << std::endl;
 std::cout << " Metric value  = " << bestValue          << std::endl;

 //  It is common, as the last step of a registration task, to use the
 //  resulting transform to map the moving image into the fixed image space.
 //  This is easily done with the \doxygen{ResampleImageFilter}.

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

 ResampleFilterType::Pointer resampler = ResampleFilterType::New();
 resampler->SetInput( movingImage);

 //  The Transform that is produced as output of the Registration method is
 //  also passed as input to the resampling filter. Note the use of the
 //  methods \code{GetOutput()} and \code{Get()}. This combination is needed
 //  here because the registration method acts as a filter whose output is a
 //  transform decorated in the form of a \doxygen{DataObject}. For details in
 //  this construction you may want to read the documentation of the
 //  \doxygen{DataObjectDecorator}.

 resampler->SetTransform( registration->GetOutput()->Get() );

 //  As described in Section \ref{sec:ResampleImageFilter}, the
 //  ResampleImageFilter requires additional parameters to be specified, in
 //  particular, the spacing, origin and size of the output image. The default
 //  pixel value is also set to a distinct gray level in order to highlight
 //  the regions that are mapped outside of the moving image.

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

 //  The output of the filter is passed to a writer that will store the
 //  image in a file. An \doxygen{CastImageFilter} is used to convert the
 //  pixel type of the resampled image to the final type used by the
 //  writer. The cast and writer filters are instantiated below.

 typedef unsigned char OutputPixelType;
 typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
 typedef itk::CastImageFilter<
     ImageType,
     ImageType > CastFilterType;

 WriterType::Pointer      writer =  WriterType::New();
 CastFilterType::Pointer  caster =  CastFilterType::New();
 writer->SetFileName("output.png");

 caster->SetInput( resampler->GetOutput() );
 writer->SetInput( caster->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 CreateSphereImage(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(ImageRegistrationMethodAffine)

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

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

</source>