<br><div class="gmail_quote"><br>Hi there,<br>I am trying to register two 2-D images,<br>using the NormalizedMutualInformationMetric (NMI), but I have terrible results.<br>I have modified the ImageRegistration7.cxx example, in order to use this metric.<br>
My questions are:<br>1.Does NMI Metric works only with itkOnePlusOneEvolutionaryOptimizer as shown in example 14?<br>2.Beside the number of histogram bins, is there any other parameter of the NMI,that could be modified,<br>
in order to get better results?<br>3.Does the images need to be preprocessed (normalizedImageFilter , and Smoothed (GaussianImageFilter)),<br>before entering the registration process?<br>Thanks in advance<br>Giorgos<br><br>
<br>/*========================================================================*/<br><br>#include "itkImageRegistrationMethod.h"<br>#include "itkNormalizedMutualInformationHistogramImageToImageMetric.h"<br>
#include "itkLinearInterpolateImageFunction.h"<br>#include "itkRegularStepGradientDescentOptimizer.h"<br>#include "itkImage.h"<br><br>#include "itkCenteredTransformInitializer.h"<br>
#include "itkCenteredSimilarity2DTransform.h"<br><br>#include "itkImageFileReader.h"<br>#include "itkImageFileWriter.h"<br><br>#include "itkResampleImageFilter.h"<br>#include "itkCastImageFilter.h"<br>
#include "itkSubtractImageFilter.h"<br>#include "itkRescaleIntensityImageFilter.h"<br>#include "itkIdentityTransform.h"<br><br><br>#include "itkMeanSquaresImageToImageMetric.h"<br>
#include "itkCorrelationCoefficientHistogramImageToImageMetric.h"<br>
<br><br>#include "itkTranslationTransform.h"<br>#include "itkNearestNeighborInterpolateImageFunction.h"<br>#include "itkCheckerBoardImageFilter.h"<br><br>//<br>#include "itkCommand.h"<br>
class CommandIterationUpdate : public itk::Command <br>{<br>public:<br> typedef CommandIterationUpdate Self;<br> typedef itk::Command Superclass;<br> typedef itk::SmartPointer<Self> Pointer;<br>
itkNewMacro( Self );<br>
protected:<br> CommandIterationUpdate() {m_LastMetricValue = 0;};<br>public:<br> typedef itk::RegularStepGradientDescentOptimizer OptimizerType;<br> typedef const OptimizerType * OptimizerPointer;<br><br> void Execute(itk::Object *caller, const itk::EventObject & event)<br>
{<br> Execute( (const itk::Object *)caller, event);<br> }<br><br> void Execute(const itk::Object * object, const itk::EventObject & event)<br> {<br> OptimizerPointer optimizer = <br> dynamic_cast< OptimizerPointer >( object );<br>
if( ! itk::IterationEvent().CheckEvent( &event ) )<br> {<br> return;<br> }<br> double currentValue = optimizer->GetValue();<br> if( fabs( m_LastMetricValue - currentValue ) > 1e-4 )<br>
{ <br> std::cout << optimizer->GetCurrentIteration() << " ";<br> std::cout << optimizer->GetValue() << " ";<br> std::cout << optimizer->GetCurrentPosition() << std::endl;<br>
m_LastMetricValue = currentValue;<br> }<br> }<br> private:<br> double m_LastMetricValue;<br>};<br><br><br>int main( int argc, char *argv[] )<br>{<br> if( argc < 4 )<br> {<br> std::cerr << "Missing Parameters " << std::endl;<br>
std::cerr << "Usage: " << argv[0];<br> std::cerr << " fixedImageFile movingImageFile ";<br> std::cerr << " outputImagefile [differenceBeforeRegistration] ";<br>
std::cerr << " [differenceAfterRegistration] [Checkerboard Image]";<br><br> std::cerr << " [initialScaling] [initialAngle] ";<br> std::cerr << " [steplength] ";<br>
std::cerr << std::endl;<br> return EXIT_FAILURE;<br> }<br> <br> const unsigned int Dimension = 2;<br> typedef float PixelType;<br> typedef float InternalPixelType;<br> <br> <br> typedef itk::Image< PixelType, Dimension > FixedImageType;<br>
typedef itk::Image< PixelType, Dimension > MovingImageType;<br><br>/// reading the images<br><br> typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;<br> typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;<br>
<br> FixedImageReaderType::Pointer fixedImageReader = FixedImageReaderType::New();<br> MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();<br><br> fixedImageReader->SetFileName( argv[1] );<br>
movingImageReader->SetFileName( argv[2] );<br><br> typedef itk::CenteredSimilarity2DTransform< double > TransformType;<br> <br> typedef itk::RegularStepGradientDescentOptimizer OptimizerType;<br> typedef itk::NormalizedMutualInformationHistogramImageToImageMetric< InternalImageType, InternalImageType><br>
MetricType;<br> typedef itk:: LinearInterpolateImageFunction< InternalImageType, double ><br> InterpolatorType;<br> typedef itk::ImageRegistrationMethod< InternalImageType, InternalImageType ><br> RegistrationType;<br>
<br> MetricType::Pointer metric = MetricType::New();<br> OptimizerType::Pointer optimizer = OptimizerType::New();<br> InterpolatorType::Pointer interpolator = InterpolatorType::New();<br> RegistrationType::Pointer registration = RegistrationType::New();<br>
<br> registration->SetMetric( metric );<br> registration->SetOptimizer( optimizer );<br> registration->SetInterpolator( interpolator );<br><br> TransformType::Pointer transform = TransformType::New();<br>
registration->SetTransform( transform );<br><br><br><br> /////NORMALIZED MI SETTINGS<br> unsigned int numberOfHistogramBins;<br> std::cout<< " Enter Number of Histogram Bins : " ;<br> std::cin >> numberOfHistogramBins; <br>
std::cout<<std::endl;<br> <br> /*if( argc > 4 )<br> {<br> numberOfHistogramBins = atoi( argv[4] );<br> std::cout << "Using " << numberOfHistogramBins << " Histogram bins" << std::endl;<br>
}*/<br> MetricType::HistogramType::SizeType histogramSize;<br> histogramSize[0] = numberOfHistogramBins;<br> histogramSize[1] = numberOfHistogramBins;<br> metric->SetHistogramSize( histogramSize );<br> <br><br>
const unsigned int numberOfParameters = transform->GetNumberOfParameters();<br><br> typedef MetricType::ScalesType ScalesType;<br> ScalesType scales( numberOfParameters );<br><br> scales.Fill( 1.0 );<br> <br> metric->SetDerivativeStepLengthScales(scales);<br>
<br>//////<br>///// SETTING THE IMAGES FOR THE REGISTRATION PROCESS<br><br> registration->SetFixedImage( fixedImageReader->GetOutput() );<br> registration->SetMovingImage( movingImageReader->GetOutput() );<br>
<br> fixedImageReader->Update();<br> FixedImageType::RegionType fixedImageRegion =<br> fixedImageReader->GetOutput()->GetBufferedRegion();<br> registration->SetFixedImageRegion( fixedImageRegion );<br>
<br><br> typedef itk::CenteredTransformInitializer< <br> TransformType, <br> FixedImageType, <br> MovingImageType > TransformInitializerType;<br>
<br> TransformInitializerType::Pointer initializer = TransformInitializerType::New();<br><br> initializer->SetTransform( transform );<br><br> initializer->SetFixedImage( fixedImageReader->GetOutput() );<br>
initializer->SetMovingImage( movingImageReader->GetOutput() );<br><br> initializer->GeometryOn();<br><br> initializer->InitializeTransform();<br><br>// REMAINING TRANSFORM PARAMETERS<br> double initialScale = 1.0;<br>
<br> if( argc > 7 )<br> {<br> initialScale = atof( argv[7] );<br> }<br> <br> double initialAngle = 0.0;<br><br> if( argc > 8 )<br> {<br> initialAngle = atof( argv[8] );<br> }<br> <br> // Software Guide : BeginCodeSnippet<br>
transform->SetScale( initialScale );<br> transform->SetAngle( initialAngle );<br> // Software Guide : EndCodeSnippet<br><br> registration->SetInitialTransformParameters( transform->GetParameters() );<br>
<br>
typedef OptimizerType::ScalesType OptimizerScalesType;<br> OptimizerScalesType optimizerScales( transform->GetNumberOfParameters() );<br> const double translationScale = 1.0 / 100.0;<br><br> optimizerScales[0] = 10.0;<br>
optimizerScales[1] = 1.0;<br> optimizerScales[2] = translationScale;<br> optimizerScales[3] = translationScale;<br> optimizerScales[4] = translationScale;<br> optimizerScales[5] = translationScale;<br><br> optimizer->SetScales( optimizerScales );<br>
<br><br> double steplength = 1.0;<br> <br> if( argc > 9 )<br> {<br> steplength = atof( argv[9] );<br> }<br><br> // Software Guide : BeginCodeSnippet<br> optimizer->SetMaximumStepLength( steplength ); <br>
optimizer->SetMinimumStepLength( 0.0001 );<br> optimizer->SetNumberOfIterations( 500 );<br> // Software Guide : EndCodeSnippet<br><br><br> // Create the Command observer and register it with the optimizer.<br>
//<br>
CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();<br> optimizer->AddObserver( itk::IterationEvent(), observer );<br><br><br> try <br> { <br> registration->StartRegistration(); <br>
} <br>
catch( itk::ExceptionObject & err ) <br> { <br> std::cerr << "ExceptionObject caught !" << std::endl; <br> std::cerr << err << std::endl; <br> return EXIT_FAILURE;<br> } <br>
<br> OptimizerType::ParametersType finalParameters = <br> registration->GetLastTransformParameters();<br><br><br> const double finalScale = finalParameters[0];<br> const double finalAngle = finalParameters[1];<br>
const double finalRotationCenterX = finalParameters[2];<br> const double finalRotationCenterY = finalParameters[3];<br> const double finalTranslationX = finalParameters[4];<br> const double finalTranslationY = finalParameters[5];<br>
<br> const unsigned int numberOfIterations = optimizer->GetCurrentIteration();<br><br> const double bestValue = optimizer->GetValue();<br><br><br> // Print out results<br> //<br> const double finalAngleInDegrees = finalAngle * 45.0 / atan(1.0);<br>
<br> std::cout << std::endl;<br> std::cout << "Result = " << std::endl;<br> std::cout << " Scale = " << finalScale << std::endl;<br> std::cout << " Angle (radians) " << finalAngle << std::endl;<br>
std::cout << " Angle (degrees) " << finalAngleInDegrees << std::endl;<br> std::cout << " Center X = " << finalRotationCenterX << std::endl;<br> std::cout << " Center Y = " << finalRotationCenterY << std::endl;<br>
std::cout << " Translation X = " << finalTranslationX << std::endl;<br> std::cout << " Translation Y = " << finalTranslationY << std::endl;<br> std::cout << " Iterations = " << numberOfIterations << std::endl;<br>
std::cout << " Metric value = " << bestValue << std::endl;<br><br><br> ...........................<br>............................<br><br> return EXIT_SUCCESS;<br>
}<br><br><br><br>
</div><br>