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

<syntaxhighlight lang="cmake"> cmake_minimum_required(VERSION 3.9.5)

project(ImageRegistrationMethodAffine)

find_package(ITK REQUIRED) include(${ITK_USE_FILE}) if (ITKVtkGlue_LOADED)

 find_package(VTK REQUIRED)
 include(${VTK_USE_FILE})

endif()

add_executable(ImageRegistrationMethodAffine MACOSX_BUNDLE ImageRegistrationMethodAffine.cxx)

if( "${ITK_VERSION_MAJOR}" LESS 4 )

 target_link_libraries(ImageRegistrationMethodAffine ITKReview ${ITK_LIBRARIES})

else( "${ITK_VERSION_MAJOR}" LESS 4 )

 target_link_libraries(ImageRegistrationMethodAffine ${ITK_LIBRARIES})

endif( "${ITK_VERSION_MAJOR}" LESS 4 )

</syntaxhighlight>

Download and Build ImageRegistrationMethodAffine

Click here to download ImageRegistrationMethodAffine and its CMakeLists.txt file. Once the tarball ImageRegistrationMethodAffine.tar has been downloaded and extracted,

cd ImageRegistrationMethodAffine/build

• If ITK is installed:

cmake ..

• If ITK is not installed but compiled on your system, you will need to specify the path to your ITK build:

cmake -DITK_DIR:PATH=/home/me/itk_build ..

Build the project:

make

and run it:

./ImageRegistrationMethodAffine

WINDOWS USERS PLEASE NOTE: Be sure to add the ITK bin directory to your path. This will resolve the ITK dll's at run time.

Building All of the Examples

Many of the examples in the ITK Wiki Examples Collection require VTK. You can build all of the the examples by following these instructions. If you are a new VTK user, you may want to try the Superbuild which will build a proper ITK and VTK.

ItkVtkGlue

ITK >= 4

For examples that use QuickView (which depends on VTK), you must have built ITK with Module_ITKVtkGlue=ON.

ITK < 4

Some of the ITK Examples require VTK to display the images. If you download the entire ITK Wiki Examples Collection, the ItkVtkGlue directory will be included and configured. If you wish to just build a few examples, then you will need to download ItkVtkGlue and build it. When you run cmake it will ask you to specify the location of the ItkVtkGlue binary directory.