|
|
| (3 intermediate revisions by one other user not shown) |
| Line 1: |
Line 1: |
| ==MutualInformation.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 "itkImageRegistrationMethod.h"
| |
| #include "itkTranslationTransform.h"
| |
| #include "itkMutualInformationImageToImageMetric.h"
| |
| #include "itkGradientDescentOptimizer.h"
| |
| #include "itkNormalizeImageFilter.h"
| |
| #include "itkDiscreteGaussianImageFilter.h"
| |
| #include "itkResampleImageFilter.h"
| |
| #include "itkCastImageFilter.h"
| |
| #include "itkCheckerBoardImageFilter.h"
| |
| #include "itkEllipseSpatialObject.h"
| |
| #include "itkSpatialObjectToImageFilter.h"
| |
| #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>
| |
