[Insight-users] NormalizedMutualInformationMetric
Luis Ibanez
luis.ibanez at kitware.com
Thu May 8 09:43:53 EDT 2008
Hi Giorgos,
Thanks for your detailed post.
Here are some answers to your questions:
1) No,
The NormalizedMutualInformationMetric it not restricted
to be used with the OnePlusOneEvolutionary optimizer.
It should work with any optimizer.
This metric provides the methods
GetValue()
GetDerivative()
GetValueAndDerivative()
so it can be used by optimizers that require
derivatives.
Note that derivatives in this metric are
computed by using finite differences.
2) Yes, there are other parameters:
a) The differential used for computing
derivatives by finite differences:
SetDerivativeStepLength()
if you pick it too small, then you
will get null derivatives
b) Lower and Upper bounds of the histograms.
These are *VERY* important. You may want
to study the intensity conten of your images
and select a range of intensities matching
the anatomical structures that you want to
register.
BTW: What are the modalities of the images
that you are trying to register ?
c) PaddingValue: in case you have pixel values
in your image that you want to be ignored
during the registration process.
3) In general it shouldn't be necessary to preprocess
your images,
*BUT*....
it all depends on the characteristics
of your specific images.
Can you tell us more about them ?
For example,
Could you post screenshots in a public web site ?
4) Can you describe better the "terribleness" of the results
that you are getting ?
a) do you get run time errors ?
b) do you get exceptions ?
c) do you get a transform that is too far off ?
d) could you post to the list the output of the Observer
that you connected to the optimizer ?
It will be useful to see the path that the optimizer
is taking.
Also,
the Similarity2D transform may be hard to control due to
the it capability for applying scaling to the image.
Do you actually need a Similarity transform for your work ?
You probably should start with a simpler transform while
you figure out how to fine tune all the other registration
parameters to your problem.
Please let us know,
Thanks
Luis
----------------------
Giorgos Pichis wrote:
>
>
> Hi there,
> I am trying to register two 2-D images,
> using the NormalizedMutualInformationMetric (NMI), but I have terrible
> results.
> I have modified the ImageRegistration7.cxx example, in order to use this
> metric.
> My questions are:
> 1.Does NMI Metric works only with itkOnePlusOneEvolutionaryOptimizer as
> shown in example 14?
> 2.Beside the number of histogram bins, is there any other parameter of
> the NMI,that could be modified,
> in order to get better results?
> 3.Does the images need to be preprocessed (normalizedImageFilter , and
> Smoothed (GaussianImageFilter)),
> before entering the registration process?
> Thanks in advance
> Giorgos
>
>
> /*========================================================================*/
>
> #include "itkImageRegistrationMethod.h"
> #include "itkNormalizedMutualInformationHistogramImageToImageMetric.h"
> #include "itkLinearInterpolateImageFunction.h"
> #include "itkRegularStepGradientDescentOptimizer.h"
> #include "itkImage.h"
>
> #include "itkCenteredTransformInitializer.h"
> #include "itkCenteredSimilarity2DTransform.h"
>
> #include "itkImageFileReader.h"
> #include "itkImageFileWriter.h"
>
> #include "itkResampleImageFilter.h"
> #include "itkCastImageFilter.h"
> #include "itkSubtractImageFilter.h"
> #include "itkRescaleIntensityImageFilter.h"
> #include "itkIdentityTransform.h"
>
>
> #include "itkMeanSquaresImageToImageMetric.h"
> #include "itkCorrelationCoefficientHistogramImageToImageMetric.h"
>
>
> #include "itkTranslationTransform.h"
> #include "itkNearestNeighborInterpolateImageFunction.h"
> #include "itkCheckerBoardImageFilter.h"
>
> //
> #include "itkCommand.h"
> class CommandIterationUpdate : public itk::Command
> {
> public:
> typedef CommandIterationUpdate Self;
> typedef itk::Command Superclass;
> typedef itk::SmartPointer<Self> Pointer;
> itkNewMacro( Self );
> protected:
> CommandIterationUpdate() {m_LastMetricValue = 0;};
> public:
> typedef itk::RegularStepGradientDescentOptimizer OptimizerType;
> typedef const OptimizerType * OptimizerPointer;
>
> void Execute(itk::Object *caller, const itk::EventObject & event)
> {
> Execute( (const itk::Object *)caller, event);
> }
>
> void Execute(const itk::Object * object, const itk::EventObject & event)
> {
> OptimizerPointer optimizer =
> dynamic_cast< OptimizerPointer >( object );
> if( ! itk::IterationEvent().CheckEvent( &event ) )
> {
> return;
> }
> double currentValue = optimizer->GetValue();
> if( fabs( m_LastMetricValue - currentValue ) > 1e-4 )
> {
> std::cout << optimizer->GetCurrentIteration() << " ";
> std::cout << optimizer->GetValue() << " ";
> std::cout << optimizer->GetCurrentPosition() << std::endl;
> m_LastMetricValue = currentValue;
> }
> }
> private:
> double m_LastMetricValue;
> };
>
>
> int main( int argc, char *argv[] )
> {
> if( argc < 4 )
> {
> std::cerr << "Missing Parameters " << std::endl;
> std::cerr << "Usage: " << argv[0];
> std::cerr << " fixedImageFile movingImageFile ";
> std::cerr << " outputImagefile [differenceBeforeRegistration] ";
> std::cerr << " [differenceAfterRegistration] [Checkerboard Image]";
>
> std::cerr << " [initialScaling] [initialAngle] ";
> std::cerr << " [steplength] ";
> std::cerr << std::endl;
> return EXIT_FAILURE;
> }
>
> const unsigned int Dimension = 2;
> typedef float PixelType;
> typedef float InternalPixelType;
>
>
> typedef itk::Image< PixelType, Dimension > FixedImageType;
> typedef itk::Image< PixelType, Dimension > MovingImageType;
>
> /// reading the images
>
> typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;
> typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
>
> FixedImageReaderType::Pointer fixedImageReader =
> FixedImageReaderType::New();
> MovingImageReaderType::Pointer movingImageReader =
> MovingImageReaderType::New();
>
> fixedImageReader->SetFileName( argv[1] );
> movingImageReader->SetFileName( argv[2] );
>
> typedef itk::CenteredSimilarity2DTransform< double > TransformType;
>
> typedef itk::RegularStepGradientDescentOptimizer OptimizerType;
> typedef itk::NormalizedMutualInformationHistogramImageToImageMetric<
> InternalImageType, InternalImageType>
> MetricType;
> typedef itk:: LinearInterpolateImageFunction< InternalImageType, double >
> InterpolatorType;
> typedef itk::ImageRegistrationMethod< InternalImageType,
> InternalImageType >
> RegistrationType;
>
> MetricType::Pointer metric = MetricType::New();
> OptimizerType::Pointer optimizer = OptimizerType::New();
> InterpolatorType::Pointer interpolator = InterpolatorType::New();
> RegistrationType::Pointer registration = RegistrationType::New();
>
> registration->SetMetric( metric );
> registration->SetOptimizer( optimizer );
> registration->SetInterpolator( interpolator );
>
> TransformType::Pointer transform = TransformType::New();
> registration->SetTransform( transform );
>
>
>
> /////NORMALIZED MI SETTINGS
> unsigned int numberOfHistogramBins;
> std::cout<< " Enter Number of Histogram Bins : " ;
> std::cin >> numberOfHistogramBins;
> std::cout<<std::endl;
>
> /*if( argc > 4 )
> {
> numberOfHistogramBins = atoi( argv[4] );
> std::cout << "Using " << numberOfHistogramBins << " Histogram bins"
> << std::endl;
> }*/
> MetricType::HistogramType::SizeType histogramSize;
> histogramSize[0] = numberOfHistogramBins;
> histogramSize[1] = numberOfHistogramBins;
> metric->SetHistogramSize( histogramSize );
>
>
> const unsigned int numberOfParameters =
> transform->GetNumberOfParameters();
>
> typedef MetricType::ScalesType ScalesType;
> ScalesType scales( numberOfParameters );
>
> scales.Fill( 1.0 );
>
> metric->SetDerivativeStepLengthScales(scales);
>
> //////
> ///// SETTING THE IMAGES FOR THE REGISTRATION PROCESS
>
> registration->SetFixedImage( fixedImageReader->GetOutput() );
> registration->SetMovingImage( movingImageReader->GetOutput() );
>
> fixedImageReader->Update();
> FixedImageType::RegionType fixedImageRegion =
> fixedImageReader->GetOutput()->GetBufferedRegion();
> registration->SetFixedImageRegion( fixedImageRegion );
>
>
> typedef itk::CenteredTransformInitializer<
> TransformType,
> FixedImageType,
> MovingImageType >
> TransformInitializerType;
>
> TransformInitializerType::Pointer initializer =
> TransformInitializerType::New();
>
> initializer->SetTransform( transform );
>
> initializer->SetFixedImage( fixedImageReader->GetOutput() );
> initializer->SetMovingImage( movingImageReader->GetOutput() );
>
> initializer->GeometryOn();
>
> initializer->InitializeTransform();
>
> // REMAINING TRANSFORM PARAMETERS
> double initialScale = 1.0;
>
> if( argc > 7 )
> {
> initialScale = atof( argv[7] );
> }
>
> double initialAngle = 0.0;
>
> if( argc > 8 )
> {
> initialAngle = atof( argv[8] );
> }
>
> // Software Guide : BeginCodeSnippet
> transform->SetScale( initialScale );
> transform->SetAngle( initialAngle );
> // Software Guide : EndCodeSnippet
>
> registration->SetInitialTransformParameters( transform->GetParameters() );
>
> typedef OptimizerType::ScalesType OptimizerScalesType;
> OptimizerScalesType optimizerScales( transform->GetNumberOfParameters() );
> const double translationScale = 1.0 / 100.0;
>
> optimizerScales[0] = 10.0;
> optimizerScales[1] = 1.0;
> optimizerScales[2] = translationScale;
> optimizerScales[3] = translationScale;
> optimizerScales[4] = translationScale;
> optimizerScales[5] = translationScale;
>
> optimizer->SetScales( optimizerScales );
>
>
> double steplength = 1.0;
>
> if( argc > 9 )
> {
> steplength = atof( argv[9] );
> }
>
> // Software Guide : BeginCodeSnippet
> optimizer->SetMaximumStepLength( steplength );
> optimizer->SetMinimumStepLength( 0.0001 );
> optimizer->SetNumberOfIterations( 500 );
> // Software Guide : EndCodeSnippet
>
>
> // Create the Command observer and register it with the optimizer.
> //
> CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
> optimizer->AddObserver( itk::IterationEvent(), observer );
>
>
> try
> {
> registration->StartRegistration();
> }
> catch( itk::ExceptionObject & err )
> {
> std::cerr << "ExceptionObject caught !" << std::endl;
> std::cerr << err << std::endl;
> return EXIT_FAILURE;
> }
>
> OptimizerType::ParametersType finalParameters =
> registration->GetLastTransformParameters();
>
>
> const double finalScale = finalParameters[0];
> const double finalAngle = finalParameters[1];
> const double finalRotationCenterX = finalParameters[2];
> const double finalRotationCenterY = finalParameters[3];
> const double finalTranslationX = finalParameters[4];
> const double finalTranslationY = finalParameters[5];
>
> const unsigned int numberOfIterations = optimizer->GetCurrentIteration();
>
> const double bestValue = optimizer->GetValue();
>
>
> // Print out results
> //
> const double finalAngleInDegrees = finalAngle * 45.0 / atan(1.0);
>
> std::cout << std::endl;
> std::cout << "Result = " << std::endl;
> std::cout << " Scale = " << finalScale << std::endl;
> std::cout << " Angle (radians) " << finalAngle << std::endl;
> std::cout << " Angle (degrees) " << finalAngleInDegrees << std::endl;
> std::cout << " Center X = " << finalRotationCenterX << std::endl;
> std::cout << " Center Y = " << finalRotationCenterY << std::endl;
> std::cout << " Translation X = " << finalTranslationX << std::endl;
> std::cout << " Translation Y = " << finalTranslationY << std::endl;
> std::cout << " Iterations = " << numberOfIterations << std::endl;
> std::cout << " Metric value = " << bestValue << std::endl;
>
>
> ...........................
> ............................
>
> return EXIT_SUCCESS;
> }
>
>
>
>
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users
More information about the Insight-users
mailing list