|
|
| (11 intermediate revisions by 2 users not shown) |
| Line 1: |
Line 1: |
| ==ImageRegistrationMethodBSpline.cxx==
| | {{warning|1=The media wiki content on this page is no longer maintained. The examples presented on the https://itk.org/Wiki/* pages likely require ITK version 4.13 or earlier releases. In many cases, the examples on this page no longer conform to the best practices for modern ITK versions.}} |
| <source lang="cpp">
| |
| #include "itkCastImageFilter.h"
| |
| #include "itkSpatialObjectToImageFilter.h"
| |
| #include "itkEllipseSpatialObject.h"
| |
| #include "itkImage.h"
| |
| #include "itkImageRegistrationMethod.h"
| |
| #include "itkLinearInterpolateImageFunction.h"
| |
| #include "itkImageFileReader.h"
| |
| #include "itkImageFileWriter.h"
| |
| #include "itkMeanSquaresImageToImageMetric.h"
| |
| #include "itkLBFGSOptimizer.h"
| |
| #include "itkResampleImageFilter.h"
| |
| #include "itkRescaleIntensityImageFilter.h"
| |
| #include "itkSpatialObjectToImageFilter.h"
| |
| #include "itkBSplineDeformableTransform.h"
| |
| | |
| const unsigned int Dimension = 2;
| |
| typedef unsigned char PixelType;
| |
| | |
| typedef itk::Image< PixelType, Dimension > ImageType;
| |
| | |
| void CreateEllipseImage(ImageType::Pointer image);
| |
| void CreateSphereImage(ImageType::Pointer image);
| |
| | |
| int main(int, char *[] )
| |
| { | |
| // The transform that will map the fixed image into the moving image.
| |
| const unsigned int SpaceDimension = Dimension;
| |
| const unsigned int SplineOrder = 3;
| |
| typedef double CoordinateRepType;
| |
| | |
| typedef itk::BSplineDeformableTransform<
| |
| CoordinateRepType,
| |
| SpaceDimension,
| |
| SplineOrder > TransformType;
| |
|
| |
| // An optimizer is required to explore the parameter space of the transform
| |
| // in search of optimal values of the metric.
| |
| typedef itk::LBFGSOptimizer 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() );
| |
| | |
| initialParameters.Fill(0.0);
| |
| | |
| registration->SetInitialTransformParameters( initialParameters );
| |
| | |
| optimizer->SetGradientConvergenceTolerance( 0.05 );
| |
| optimizer->SetLineSearchAccuracy( 0.9 );
| |
| //optimizer->SetDefaultStepLength( 1.5 );
| |
| optimizer->SetDefaultStepLength( .5 );
| |
| optimizer->TraceOn();
| |
| optimizer->SetMaximumNumberOfFunctionEvaluations( 1000 );
| |
| | |
| 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(BSpline)
| |
| | |
| FIND_PACKAGE(ITK REQUIRED)
| |
| INCLUDE(${ITK_USE_FILE})
| |
| | |
| ADD_EXECUTABLE(BSpline BSpline.cxx)
| |
| TARGET_LINK_LIBRARIES(BSpline ITKIO ITKNumerics)
| |
| | |
| | |
| </source>
| |
