Exhaustive Optimizer

Synopsis

An optimizer that fully samples a grid on the parametric space.

Results

Note

Help Wanted Implementation of Results for sphinx examples containing this message. Reconfiguration of CMakeList.txt may be necessary. Write An Example <https://itk.org/ITKExamples/Documentation/Contribute/WriteANewExample.html>

Code

C++

#include "itkImageFileReader.h"
#include "itkEuler2DTransform.h"
#include "itkExhaustiveOptimizerv4.h"
#include "itkMeanSquaresImageToImageMetricv4.h"
#include "itkCenteredTransformInitializer.h"
#include "itkImageRegistrationMethodv4.h"
#include "itkImage.h"

//  Command observer to monitor the evolution of the registration process.
//
#include "itkCommand.h"
class CommandIterationUpdate : public itk::Command
{
public:
  using Self = CommandIterationUpdate;
  using Superclass = itk::Command;
  using Pointer = itk::SmartPointer<Self>;
  itkNewMacro(Self);

protected:
  CommandIterationUpdate() = default;

public:
  using OptimizerType = itk::ExhaustiveOptimizerv4<double>;
  using OptimizerPointer = const OptimizerType *;

  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
  {
    auto optimizer = static_cast<OptimizerPointer>(object);
    if (!(itk::IterationEvent().CheckEvent(&event)))
    {
      return;
    }
    std::cout << "Iteration: ";
    std::cout << optimizer->GetCurrentIteration() << "   ";
    std::cout << optimizer->GetCurrentIndex() << "   ";
    std::cout << optimizer->GetCurrentValue() << "   ";
    std::cout << optimizer->GetCurrentPosition() << "   ";
    std::cout << std::endl;
  }
};

int
main(int argc, char * argv[])
{
  if (argc < 3)
  {
    std::cout << "Usage: " << argv[0] << " fixedImage movingImage" << std::endl;
    return EXIT_FAILURE;
  }
  using FixedImageType = itk::Image<double, 2>;
  using MovingImageType = itk::Image<double, 2>;
  using FixedImageReaderType = itk::ImageFileReader<FixedImageType>;
  using MovingImageReaderType = itk::ImageFileReader<MovingImageType>;
  using TransformType = itk::Euler2DTransform<double>;
  using OptimizerType = itk::ExhaustiveOptimizerv4<double>;
  using MetricType = itk::MeanSquaresImageToImageMetricv4<FixedImageType, MovingImageType>;
  using TransformInitializerType = itk::CenteredTransformInitializer<TransformType, FixedImageType, MovingImageType>;
  using RegistrationType = itk::ImageRegistrationMethodv4<FixedImageType, MovingImageType, TransformType>;

  FixedImageReaderType::Pointer     fixedImageReader = FixedImageReaderType::New();
  MovingImageReaderType::Pointer    movingImageReader = MovingImageReaderType::New();
  FixedImageType::Pointer           fixedImage = FixedImageType::New();
  MovingImageType::Pointer          movingImage = MovingImageType::New();
  TransformType::Pointer            transform = TransformType::New();
  MetricType::Pointer               metric = MetricType::New();
  OptimizerType::Pointer            optimizer = OptimizerType::New();
  RegistrationType::Pointer         registration = RegistrationType::New();
  TransformInitializerType::Pointer initializer = TransformInitializerType::New();

  fixedImageReader->SetFileName(argv[1]);
  fixedImageReader->Update();
  fixedImage = fixedImageReader->GetOutput();

  movingImageReader->SetFileName(argv[2]);
  movingImageReader->Update();
  movingImage = movingImageReader->GetOutput();

  // Create the Command observer and register it with the optimizer.
  //
  CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
  optimizer->AddObserver(itk::IterationEvent(), observer);

  unsigned int             angles = 12;
  OptimizerType::StepsType steps(transform->GetNumberOfParameters());
  steps[0] = int(angles / 2);
  steps[1] = 0;
  steps[2] = 0;
  optimizer->SetNumberOfSteps(steps);

  OptimizerType::ScalesType scales(transform->GetNumberOfParameters());
  scales[0] = 2.0 * itk::Math::pi / angles;
  scales[1] = 1.0;
  scales[2] = 1.0;

  optimizer->SetScales(scales);

  initializer->SetTransform(transform);
  initializer->SetFixedImage(fixedImage);
  initializer->SetMovingImage(movingImage);
  initializer->InitializeTransform();

  // Initialize registration
  registration->SetMetric(metric);
  registration->SetOptimizer(optimizer);
  registration->SetFixedImage(fixedImage);
  registration->SetMovingImage(movingImage);
  registration->SetInitialTransform(transform);
  registration->SetNumberOfLevels(1);
  try
  {
    registration->Update();
    std::cout << "  MinimumMetricValue: " << optimizer->GetMinimumMetricValue() << std::endl;
    std::cout << "  MaximumMetricValue: " << optimizer->GetMaximumMetricValue() << std::endl;
    std::cout << "  MinimumMetricValuePosition: " << optimizer->GetMinimumMetricValuePosition() << std::endl;
    std::cout << "  MaximumMetricValuePosition: " << optimizer->GetMaximumMetricValuePosition() << std::endl;
    std::cout << "  StopConditionDescription: " << optimizer->GetStopConditionDescription() << std::endl;
  }
  catch (itk::ExceptionObject & err)
  {
    std::cerr << "ExceptionObject caught !" << std::endl;
    std::cerr << err << std::endl;
    return EXIT_FAILURE;
  }
  return EXIT_SUCCESS;
}

Classes demonstrated

class ExhaustiveOptimizer : public itk::SingleValuedNonLinearOptimizer

Optimizer that fully samples a grid on the parametric space.

This optimizer is equivalent to an exhaustive search in a discrete grid defined over the parametric space. The grid is centered on the initial position. The subdivisions of the grid along each one of the dimensions of the parametric space is defined by an array of number of steps.

A typical use is to plot the metric space to get an idea of how noisy it is. An example is given below, where it is desired to plot the metric space with respect to translations along x, y and z in a 3D registration application: Here it is assumed that the transform is Euler3DTransform.

OptimizerType::StepsType steps( m_Transform->GetNumberOfParameters() );
steps[0] = 10;
steps[1] = 10;
steps[2] = 10;
m_Optimizer->SetNumberOfSteps( steps );
m_Optimizer->SetStepLength( 2 );

The optimizer throws IterationEvents after every iteration. We use this to plot the metric space in an image as follows:

if( itk::IterationEvent().CheckEvent(& event ) )
{
  IndexType index;
  index[0] = m_Optimizer->GetCurrentIndex()[0];
  index[1] = m_Optimizer->GetCurrentIndex()[1];
  index[2] = m_Optimizer->GetCurrentIndex()[2];
  image->SetPixel( index, m_Optimizer->GetCurrentValue() );
}

The image size is expected to be 11 x 11 x 11.

If you wish to use different step lengths along each parametric axis, you can use the SetScales() method. This accepts an array, each element represents the number of subdivisions per step length. For instance scales of [0.5 1 4] along with a step length of 2 will cause the optimizer to search the metric space on a grid with x,y,z spacing of [1 2 8].

Physical dimensions of the grid are influenced by both the scales and the number of steps along each dimension, a side of the region is stepLength*(2*numberOfSteps[d]+1)*scaling[d].

ITK Sphinx Examples:

See itk::ExhaustiveOptimizer for additional documentation.