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

**URL:** https://discourse.itk.org/t/itk-gradientdescentoptimizer-seems-not-to-return-deterministic-results-with-identical-start-conditions/1684
**Category:** Algorithms
**Created:** [March 19, 2019, 8:51am UTC](https://discourse.itk.org/t/itk-gradientdescentoptimizer-seems-not-to-return-deterministic-results-with-identical-start-conditions/1684 "2019-03-19T08:51:54Z")
**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 19, 2019, 8:51am UTC](https://discourse.itk.org/t/itk-gradientdescentoptimizer-seems-not-to-return-deterministic-results-with-identical-start-conditions/1684/1 "2019-03-19T08:51:54Z")

</div>

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???

---

<div class="post-metadata">

### Author: ![blowekamp](https://discourse.itk.org/user_avatar/discourse.itk.org/blowekamp/32/79_2.png) [@blowekamp](https://discourse.itk.org/u/blowekamp)
#### Post date: [March 19, 2019, 12:57pm UTC](https://discourse.itk.org/t/itk-gradientdescentoptimizer-seems-not-to-return-deterministic-results-with-identical-start-conditions/1684/2 "2019-03-19T12:57:02Z")

</div>

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](https://itk.org/Doxygen/html/classitk_1_1ImageRegistrationMethodv4.html#aecfd9f781b4d0d6d4198df96f0e787f9)
- 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](https://discourse.itk.org/t/reproducibility-of-image-registration/1591)