# itk::ImageRegistrationMethod does not return the best value from all iterations

**URL:** https://discourse.itk.org/t/itk-imageregistrationmethod-does-not-return-the-best-value-from-all-iterations/1681
**Category:** Algorithms
**Created:** [March 18, 2019, 8:16pm UTC](https://discourse.itk.org/t/itk-imageregistrationmethod-does-not-return-the-best-value-from-all-iterations/1681 "2019-03-18T20:16:09Z")
**Posts on this page:** 2
**Page:** 1

<div class="post-metadata">

### Author: ![donelron](https://discourse.itk.org/letter_avatar_proxy/v4/letter/d/8e7dd6/32.png) [@donelron](https://discourse.itk.org/u/donelron)
#### Post date: [March 18, 2019, 8:16pm UTC](https://discourse.itk.org/t/itk-imageregistrationmethod-does-not-return-the-best-value-from-all-iterations/1681/1 "2019-03-18T20:16:09Z")

</div>

I am using the following lines of code (basically copied from an the ITK example [https://itk.org/Wiki/ITK/Examples/Registration/MutualInformation](https://itk.org/Wiki/ITK/Examples/Registration/MutualInformation)) to register two synthetic images with Mutual information:

```
#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"

#include "itkRescaleIntensityImageFilter.h"
#include "itkCommandIterationUpdate.h"
#include "itkMacro.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);

typedef float InternalPixelType;
typedef itk::Image< float, 2> InternalImageType;

typedef itk::ImageRegistrationMethod<
    InternalImageType,
    InternalImageType > RegistrationType;

class CommandIterationUpdate : public itk::Command
{
public:
    using Self = CommandIterationUpdate;
    using Superclass = itk::Command;
    using Pointer = itk::SmartPointer<CommandIterationUpdate>;
    itkNewMacro(CommandIterationUpdate);
protected:
    CommandIterationUpdate() = default;

    using InternalImageType = itk::Image< float, 2 >;
    using VectorPixelType = itk::Vector< float, 2 >;
    using DisplacementFieldType = itk::Image< VectorPixelType, 2 >;

    using RegistrationFilterType = itk::ImageRegistrationMethod<
        InternalImageType,
        InternalImageType >;

public:

    void Execute(itk::Object *caller, const itk::EventObject & event) override
    {
        Execute((const itk::Object *)caller, event);
    }

    void Execute(const itk::Object * object, const itk::EventObject & event) override
    {        
        const auto * optimizerMethod = static_cast<const itk::GradientDescentOptimizer *>(object);

        if (optimizerMethod == nullptr)
        {
            return;
        }
        if (!(itk::IterationEvent().CheckEvent(&event)))
        {
            return;
        }

        std::cout << "it " << optimizerMethod->GetCurrentIteration() << ", transl = " << optimizerMethod->GetCurrentPosition() 
          << ", value=" << optimizerMethod->GetValue()
          << "\n";
    }
};

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);

  // 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);

  typedef itk::TranslationTransform< double, Dimension > TransformType;
  typedef itk::GradientDescentOptimizer OptimizerType;
  typedef itk::LinearInterpolateImageFunction< InternalImageType, double > InterpolatorType;
  
  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 );

  metric->SetFixedImageStandardDeviation( 0.4 );
  metric->SetMovingImageStandardDeviation( 0.4 );

  registration->SetFixedImage( fixedNormalizer->GetOutput() );
  registration->SetMovingImage(movingNormalizer->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 );  
  const unsigned int numberOfPixels = fixedImageRegion.GetNumberOfPixels();
  const unsigned int numberOfSamples = static_cast<unsigned int>(numberOfPixels * 0.01);
  metric->SetNumberOfSpatialSamples( numberOfSamples );

  optimizer->SetLearningRate(50.0);
  optimizer->SetNumberOfIterations(20);
  optimizer->MaximizeOn(); // We want to maximize mutual information (the default of the optimizer is to minimize)

  // Connect an observer
  CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
  optimizer->AddObserver( itk::IterationEvent(), observer );

  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;

  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());
}

```

Now when I look at the output it seems that the result of the last iteration (in this example **0.000934234** ) is always used as “best value”, even though there are better options: e.g. lines 3 and 4 (ie 0.00124208 and 0.00240695). Of course I can put up some mechanism to keep track of all the results of all iterations and then at the very end return the best one. But this is probably not how it was intended to be used?!? There is a function  
virtual const ParametersType& itk::ImageRegistrationMethod\< TFixedImage, TMovingImage \>::GetLastTransformParameters ( ) const  
. However, I would have expected a function **GetBestTransformParameters()** to actually return the translation that resulted in the max. mutual information. What am I missing here?!?

it 0, transl = [3.080745490862392e-27, -2.716527225092382e-25], value=0.000612325  
it 1, transl = [5.1655689727723054e-23, 3.0754345820462986e-23], value=-0.00153213  
it 2, transl = [4.7576327557073874e-23, 3.0754345820462986e-23], value=0.00124208  
it 3, transl = [2.2864627443559803e-23, 3.0728931531586856e-23], value=0.00240695  
it 4, transl = [-4.666378034618216e-19, 2.1827462357148532e-23], value=-0.0626463  
it 5, transl = [-4.66734633157285e-19, -6.388334110244383e-23], value=-0.000920195  
it 6, transl = [-4.654328593482208e-19, -6.587315526432373e-23], value=0.0228237  
it 7, transl = [-4.659618853553521e-19, 1.434332129935076e-21], value=-0.0218188  
it 8, transl = [-4.659649774797798e-19, 1.4587940542757805e-21], value=-0.14773  
it 9, transl = [-4.660465615278749e-19, 1.4587889196999623e-21], value=0.00145463  
it 10, transl = [-4.661200248067382e-19, 1.4587889196999623e-21], value=-0.00058087  
it 11, transl = [-4.658168507396112e-19, 1.4685938068876238e-21], value=0.00610458  
it 12, transl = [-4.658170963459659e-19, 1.4661027965906913e-21], value=6.32885e-05  
it 13, transl = [-4.657769902544275e-19, 1.4630956262500453e-21], value=0.0076189  
it 14, transl = [-4.657769902544275e-19, 1.4630956262500453e-21], value=-0.076171  
it 15, transl = [-4.657769587250847e-19, 1.4630956262500453e-21], value=0.00135417  
it 16, transl = [-4.657769587250847e-19, 1.4629990013936299e-21], value=0.000306046  
it 17, transl = [-4.658379842485591e-19, 1.4629990013936299e-21], value=-0.00570804  
it 18, transl = [-9.334279237000942e-19, 1.6765176189990196e-21], value=-0.196626  
it 19, transl = [-9.334320099934345e-19, 1.6765106878293262e-21], value= **0.000934234**  
Optimizer stop condition: GradientDescentOptimizer: Maximum number of iterations (20) exceeded.

Result =  
Translation X = -9.33432e-19  
Translation Y = 1.67651e-21  
Iterations = 20  
Metric value = **0.000934234**  
Numb. Samples = 100

---

<div class="post-metadata">

### Author: ![dzenanz](https://discourse.itk.org/user_avatar/discourse.itk.org/dzenanz/32/1093_2.png) [@dzenanz](https://discourse.itk.org/u/dzenanz)
#### Post date: [March 19, 2019, 1:28pm UTC](https://discourse.itk.org/t/itk-imageregistrationmethod-does-not-return-the-best-value-from-all-iterations/1681/2 "2019-03-19T13:28:26Z")

</div>

Try calling [optimizer-\>SetReturnBestParametersAndValue(true);](https://itk.org/Doxygen/html/classitk_1_1GradientDescentOptimizerv4Template.html#a126e05fd40f1527f48929ac89c593c89)