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>
Jupyter Notebook#
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 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>;
auto fixedImage = FixedImageType::New();
auto movingImage = MovingImageType::New();
auto transform = TransformType::New();
auto metric = MetricType::New();
auto optimizer = OptimizerType::New();
auto registration = RegistrationType::New();
auto initializer = TransformInitializerType::New();
fixedImage = itk::ReadImage<FixedImageType>(argv[1]);
movingImage = itk::ReadImage<MovingImageType>(argv[2]);
// Create the Command observer and register it with the optimizer.
//
auto 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 (const itk::ExceptionObject & err)
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
Python#
# Python example demonstrating itk.ExhaustiveOptimizer usage
import sys
from math import pi
import itk
EXIT_SUCCESS = 0
EXIT_FAILURE = 1
def main(argv):
if len(argv) < 3:
raise Exception(f"Usage: {argv[0]} fixed_image moving_image")
fixed_image = itk.imread(argv[1], itk.F)
moving_image = itk.imread(argv[2], itk.F)
dimension = fixed_image.GetImageDimension()
FixedImageType = type(fixed_image)
MovingImageType = type(moving_image)
if dimension == 2:
transform = itk.Euler2DTransform[itk.D].New()
else:
raise Exception(f"Unsupported dimension: {dimension}")
TransformInitializerType = itk.CenteredTransformInitializer[
itk.MatrixOffsetTransformBase[itk.D, dimension, dimension],
FixedImageType,
MovingImageType,
]
initializer = TransformInitializerType.New(
Transform=transform,
FixedImage=fixed_image,
MovingImage=moving_image,
)
initializer.InitializeTransform()
metric = itk.MeanSquaresImageToImageMetricv4[FixedImageType, MovingImageType].New()
optimizer = itk.ExhaustiveOptimizerv4[itk.D].New()
def print_iteration():
print(
f"Iteration:"
f"{optimizer.GetCurrentIteration()}\t"
f"{list(optimizer.GetCurrentIndex())} \t"
f"{optimizer.GetCurrentValue():10.4f}\t"
f"{list(optimizer.GetCurrentPosition())}\t"
)
optimizer.AddObserver(itk.IterationEvent(), print_iteration)
angles = 12
optimizer.SetNumberOfSteps([int(angles / 2), 0, 0])
# Initialize scales and set back to optimizer
scales = optimizer.GetScales()
scales.SetSize(3)
scales.SetElement(0, 2.0 * pi / angles)
scales.SetElement(1, 1.0)
scales.SetElement(2, 1.0)
optimizer.SetScales(scales)
RegistrationType = itk.ImageRegistrationMethodv4[FixedImageType, MovingImageType]
registration = RegistrationType.New(
Metric=metric,
Optimizer=optimizer,
FixedImage=fixed_image,
MovingImage=moving_image,
InitialTransform=transform,
NumberOfLevels=1,
)
try:
registration.Update()
print(
f"MinimumMetricValue: {optimizer.GetMinimumMetricValue():.4f}\t"
f"MaximumMetricValue: {optimizer.GetMaximumMetricValue():.4f}\n"
f"MinimumMetricValuePosition: {list(optimizer.GetMinimumMetricValuePosition())}\t"
f"MaximumMetricValuePosition: {list(optimizer.GetMaximumMetricValuePosition())}\n"
f"StopConditionDescription: {optimizer.GetStopConditionDescription()}\t"
)
except Exception as e:
print(f"Exception caught: {e}")
return EXIT_FAILURE
return EXIT_SUCCESS
if __name__ == "__main__":
main(sys.argv)
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:
