itk:GradientDescentOptimizer seems not to return deterministic results with identical start conditions

Algorithms
1

I am using the following pieceo of code which is basically copied from the ITK examples in order to perform a registration with an itk:GradientDescentOptimizer:

#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 "itkImageFileReader.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);

  typedef itk::ImageFileReader<ImageType> FileRaderTyper;
  FileRaderTyper::Pointer pFileReaderFixed = FileRaderTyper::New();
  pFileReaderFixed->SetFileName("F:/Registrierung/img/input/fixed.png");
  pFileReaderFixed->Update();
  fixedImage = pFileReaderFixed->GetOutput();

  FileRaderTyper::Pointer pFileReaderMoving = FileRaderTyper::New();
  pFileReaderMoving->SetFileName("F:/Registrierung/img/input/moving.png");
  pFileReaderMoving->Update();
  movingImage = pFileReaderMoving->GetOutput();


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

With identical start parameters I was expecting to see identical results when I start this several times, but the the results are different every time. For example the first two runs gave me

it 0, transl = [4.592627535668943e-23, 0], value=0.000764086
it 1, transl = [2.5477274229314293e-21, -1.2221097110755927e-21], value=0.0569067
Optimizer stop condition: GradientDescentOptimizer: Maximum number of iterations (2) exceeded.

Result =
Translation X = 2.54773e-21
Translation Y = -1.22211e-21
Iterations = 2
Metric value = 0.0569067
Numb. Samples = 100

and

it 0, transl = [5.543717652645973e-26, 0], value=0.00124661
it 1, transl = [1.012634919439225e-23, -3.1250277377459864e-24], value=0.0667512
Optimizer stop condition: GradientDescentOptimizer: Maximum number of iterations (2) exceeded.

Result =
Translation X = 1.01263e-23
Translation Y = -3.12503e-24
Iterations = 2
Metric value = 0.0667512
Numb. Samples = 100

What do I have to do in order to get consistent results with every run???

2

You results are very close. To put that magnitude in perspective consider the estimated size of a sub atomic quark is 1e-15 meters. For any practically driven application, this would be defined well with in acceptable tolerances to say the least.

Two things to be aware of with reproducibility of the registration:

  • Look for randomness in the components of the registration such as ImageRegistrationMethodv4:: MetricSamplingReinitializeSeed
  • Consider only using one thread to maximize reproducibility, because only within the limits floating point resolution operations are commutative.

There was also a similar recent post here: Reproducibility of image registration

2 Likes