ITK  5.2.0
Insight Toolkit
Examples/RegistrationITKv4/MultiResImageRegistration1.cxx
/*=========================================================================
*
* Copyright NumFOCUS
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
// Software Guide : BeginCommandLineArgs
// INPUTS: {BrainT1SliceBorder20.png}
// INPUTS: {BrainProtonDensitySliceShifted13x17y.png}
// OUTPUTS: {MultiResImageRegistration1Output.png}
// ARGUMENTS: 128
// OUTPUTS: {MultiResImageRegistration1CheckerboardBefore.png}
// OUTPUTS: {MultiResImageRegistration1CheckerboardAfter.png}
// Software Guide : EndCommandLineArgs
// Software Guide : BeginLatex
//
// \index{itk::ImageRegistrationMethodv4!Multi-Resolution}
// \index{itk::ImageRegistrationMethodv4!Multi-Modality}
//
// This example illustrates the use of the
// \doxygen{ImageRegistrationMethodv4} to solve a simple
// multi-modality registration problem by a multi-resolution approach.
// Since ITKv4 registration method is designed based on a multi-resolution
// structure, a separate set of classes are no longer required to run
// the registration process of this example.
//
// This a great advantage over the previous versions of ITK, as
// in ITKv3 we had to use a different filter
// (\doxygen{MultiResolutionImageRegistrationMethod})
// to run a multi-resolution process. Also, we had to use image pyramids
// filters
// (\doxygen{MultiResolutionPyramidImageFilter}) for creating the sequence of
// downsampled images. Hence, you can see how ITKv4 framework is
// more user-friendly in more complex situations.
//
// To begin the example, we include the headers of the registration
// components we will use.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The ImageRegistrationMethodv4 solves a registration
// problem in a coarse-to-fine manner as illustrated in Figure
// \ref{fig:MultiResRegistrationConcept}. The registration is first performed
// at the coarsest level using the images at the first level of the fixed and
// moving image pyramids. The transform parameters determined by the
// registration are then used to initialize the registration at the next finer
// level using images from the second level of the pyramids. This process is
// repeated as we work up to the finest level of image resolution.
//
// \begin{figure}
// \center
// \includegraphics[width=\textwidth]{MultiResRegistrationConcept}
// \itkcaption[Conceptual representation of Multi-Resolution
// registration]{Conceptual representation of the multi-resolution
// registration process.} \label{fig:MultiResRegistrationConcept} \end{figure}
//
// Software Guide : EndLatex
// Software Guide : BeginLatex
//
// In a typical registration scenario, a user will tweak component settings
// or even swap out components between multi-resolution levels. For example,
// when optimizing at a coarse resolution, it may be possible to take more
// aggressive step sizes and have a more relaxed convergence criterion.
//
// Tweaking the components between resolution levels can be done using ITK's
// implementation of the \emph{Command/Observer} design pattern. Before
// beginning registration at each resolution level,
// where ImageRegistrationMethodv4 invokes a
// \code{MultiResolutionIterationEvent()}. The registration components can
// be changed by implementing a \doxygen{Command} which responds to the
// event. A brief description of the interaction between events and commands
// was previously presented in Section \ref{sec:MonitoringImageRegistration}.
//
// We will illustrate this mechanism by changing the parameters of the
// optimizer between each resolution level by way of a simple interface
// command. First, we include the header file of the Command class.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
#include "itkCommand.h"
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Our new interface command class is called
// \code{RegistrationInterfaceCommand}. It derives from
// Command and is templated over the
// multi-resolution registration type.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
template <typename TRegistration>
class RegistrationInterfaceCommand : public itk::Command
{
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// We then define \code{Self}, \code{Superclass}, \code{Pointer},
// \code{New()} and a constructor in a similar fashion to the
// \code{CommandIterationUpdate} class in Section
// \ref{sec:MonitoringImageRegistration}.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
public:
using Self = RegistrationInterfaceCommand;
using Superclass = itk::Command;
using Pointer = itk::SmartPointer<Self>;
itkNewMacro(Self);
protected:
RegistrationInterfaceCommand() = default;
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// For convenience, we declare types useful for converting pointers
// in the \code{Execute()} method.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
public:
using RegistrationType = TRegistration;
using RegistrationPointer = RegistrationType *;
using OptimizerPointer = OptimizerType *;
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Two arguments are passed to the \code{Execute()} method: the first
// is the pointer to the object which invoked the event and the
// second is the event that was invoked.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
void
Execute(itk::Object * object, const itk::EventObject & event) override
{
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// First we verify that the event invoked is of the right type,
// \code{itk::MultiResolutionIterationEvent()}.
// If not, we return without any further action.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
if (!(itk::MultiResolutionIterationEvent().CheckEvent(&event)))
{
return;
}
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// We then convert the input object pointer to a RegistrationPointer.
// Note that no error checking is done here to verify the
// \code{dynamic\_cast} was successful since we know the actual object
// is a registration method. Then we ask for the optimizer object
// from the registration method.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
auto registration = static_cast<RegistrationPointer>(object);
auto optimizer =
static_cast<OptimizerPointer>(registration->GetModifiableOptimizer());
// Software Guide : EndCodeSnippet
unsigned int currentLevel = registration->GetCurrentLevel();
typename RegistrationType::ShrinkFactorsPerDimensionContainerType
shrinkFactors =
registration->GetShrinkFactorsPerDimension(currentLevel);
typename RegistrationType::SmoothingSigmasArrayType smoothingSigmas =
registration->GetSmoothingSigmasPerLevel();
std::cout << "-------------------------------------" << std::endl;
std::cout << " Current level = " << currentLevel << std::endl;
std::cout << " shrink factor = " << shrinkFactors << std::endl;
std::cout << " smoothing sigma = ";
std::cout << smoothingSigmas[currentLevel] << std::endl;
std::cout << std::endl;
// Software Guide : BeginLatex
//
// If this is the first resolution level we set the learning rate
// (representing the first step size) and the minimum step length
// (representing the convergence criterion) to large values. At each
// subsequent resolution level, we will reduce the minimum step length by
// a factor of 5 in order to allow the optimizer to focus on progressively
// smaller regions. The learning rate is set up to the current step
// length. In this way, when the optimizer is reinitialized at the
// beginning of the registration process for the next level, the step
// length will simply start with the last value used for the previous
// level. This will guarantee the continuity of the path taken by the
// optimizer through the parameter space.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
if (registration->GetCurrentLevel() == 0)
{
optimizer->SetLearningRate(16.00);
optimizer->SetMinimumStepLength(2.5);
}
else
{
optimizer->SetLearningRate(optimizer->GetCurrentStepLength());
optimizer->SetMinimumStepLength(optimizer->GetMinimumStepLength() *
0.2);
}
// Software Guide : EndCodeSnippet
}
// Software Guide : BeginLatex
//
// Another version of the \code{Execute()} method accepting a \code{const}
// input object is also required since this method is defined as pure
// virtual in the base class. This version simply returns without taking
// any action.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
void
Execute(const itk::Object *, const itk::EventObject &) override
{
return;
}
};
// Software Guide : EndCodeSnippet
// The following section of code implements an observer
// that will monitor the evolution of the registration process.
//
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 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 << optimizer->GetCurrentIteration() << " ";
std::cout << optimizer->GetValue() << " ";
std::cout << optimizer->GetCurrentPosition() << " ";
std::cout << m_CumulativeIterationIndex++ << std::endl;
}
private:
unsigned int m_CumulativeIterationIndex{ 0 };
};
int
main(int argc, const char * argv[])
{
if (argc < 4)
{
std::cerr << "Missing Parameters " << std::endl;
std::cerr << "Usage: " << argv[0];
std::cerr << " fixedImageFile movingImageFile ";
std::cerr << " outputImagefile [backgroundGrayLevel]";
std::cerr << " [checkerBoardBefore] [checkerBoardAfter]";
std::cerr << " [numberOfBins] " << std::endl;
return EXIT_FAILURE;
}
constexpr unsigned int Dimension = 2;
using PixelType = float;
const std::string fixedImageFile = argv[1];
const std::string movingImageFile = argv[2];
const std::string outImagefile = argv[3];
const PixelType backgroundGrayLevel = (argc > 4) ? std::stoi(argv[4]) : 100;
const std::string checkerBoardBefore = (argc > 5) ? argv[5] : "";
const std::string checkerBoardAfter = (argc > 6) ? argv[6] : "";
const int numberOfBins = (argc > 7) ? std::stoi(argv[7]) : 0;
using FixedImageType = itk::Image<PixelType, Dimension>;
using MovingImageType = itk::Image<PixelType, Dimension>;
// Software Guide : BeginLatex
//
// The fixed and moving image types are defined as in previous
// examples. The downsampled images for different resolution levels
// are created internally by the registration method based on the
// values provided for \emph{ShrinkFactor} and \emph{SmoothingSigma}
// vectors.
//
// The types for the registration components are then derived using
// the fixed and moving image type, as in previous examples.
//
// Software Guide : EndLatex
using MetricType =
MovingImageType>;
using RegistrationType = itk::
ImageRegistrationMethodv4<FixedImageType, MovingImageType, TransformType>;
// All the components are instantiated using their \code{New()} method
// and connected to the registration object as in previous example.
//
TransformType::Pointer transform = TransformType::New();
OptimizerType::Pointer optimizer = OptimizerType::New();
MetricType::Pointer metric = MetricType::New();
RegistrationType::Pointer registration = RegistrationType::New();
registration->SetOptimizer(optimizer);
registration->SetMetric(metric);
using FixedImageReaderType = itk::ImageFileReader<FixedImageType>;
using MovingImageReaderType = itk::ImageFileReader<MovingImageType>;
FixedImageReaderType::Pointer fixedImageReader =
FixedImageReaderType::New();
MovingImageReaderType::Pointer movingImageReader =
MovingImageReaderType::New();
fixedImageReader->SetFileName(fixedImageFile);
movingImageReader->SetFileName(movingImageFile);
registration->SetFixedImage(fixedImageReader->GetOutput());
registration->SetMovingImage(movingImageReader->GetOutput());
using ParametersType = OptimizerType::ParametersType;
ParametersType initialParameters(transform->GetNumberOfParameters());
initialParameters[0] = 0.0; // Initial offset in mm along X
initialParameters[1] = 0.0; // Initial offset in mm along Y
transform->SetParameters(initialParameters);
registration->SetInitialTransform(transform);
registration->InPlaceOn();
metric->SetNumberOfHistogramBins(24);
if (argc > 7)
{
// optionally, override the values with numbers taken from the command
// line arguments.
metric->SetNumberOfHistogramBins(numberOfBins);
}
// Software Guide : BeginLatex
//
// To set the optimizer parameters, note that \emph{LearningRate}
// and \emph{MinimumStepLength} are set in the obsever at the begining
// of each resolution level. The other optimizer parameters are set
// as follows.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
optimizer->SetNumberOfIterations(200);
optimizer->SetRelaxationFactor(0.5);
// Software Guide : EndCodeSnippet
// Create the Command observer and register it with the optimizer.
//
CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
optimizer->AddObserver(itk::IterationEvent(), observer);
// Software Guide : BeginLatex
//
// We set the number of multi-resolution levels to three and set
// the corresponding shrink factor and smoothing sigma values for each
// resolution level. Using smoothing in the subsampled images in
// low-resolution levels can avoid large fluctuations in the
// metric function, which prevents the optimizer from becoming trapped in
// local minima. In this simple example we have no smoothing, and we have
// used small shrinkings for the first two resolution levels.
//
// \index{itk::Image\-Registration\-Methodv4!SetNumberOfLevels()}
// \index{itk::Image\-Registration\-Methodv4!SetShrinkFactorsPerLevel()}
// \index{itk::Image\-Registration\-Methodv4!SetSmoothingSigmasPerLevel()}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
constexpr unsigned int numberOfLevels = 3;
RegistrationType::ShrinkFactorsArrayType shrinkFactorsPerLevel;
shrinkFactorsPerLevel.SetSize(3);
shrinkFactorsPerLevel[0] = 3;
shrinkFactorsPerLevel[1] = 2;
shrinkFactorsPerLevel[2] = 1;
RegistrationType::SmoothingSigmasArrayType smoothingSigmasPerLevel;
smoothingSigmasPerLevel.SetSize(3);
smoothingSigmasPerLevel[0] = 0;
smoothingSigmasPerLevel[1] = 0;
smoothingSigmasPerLevel[2] = 0;
registration->SetNumberOfLevels(numberOfLevels);
registration->SetShrinkFactorsPerLevel(shrinkFactorsPerLevel);
registration->SetSmoothingSigmasPerLevel(smoothingSigmasPerLevel);
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Once all the registration components are in place we can create
// an instance of our interface command and connect it to the
// registration object using the \code{AddObserver()} method.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
using CommandType = RegistrationInterfaceCommand<RegistrationType>;
CommandType::Pointer command = CommandType::New();
registration->AddObserver(itk::MultiResolutionIterationEvent(), command);
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Then we trigger the registration process by calling \code{Update()}.
//
// Software Guide : EndLatex
try
{
registration->Update();
std::cout << "Optimizer stop condition: "
<< registration->GetOptimizer()->GetStopConditionDescription()
<< std::endl;
}
catch (const itk::ExceptionObject & err)
{
std::cout << "ExceptionObject caught !" << std::endl;
std::cout << err << std::endl;
return EXIT_FAILURE;
}
ParametersType finalParameters = transform->GetParameters();
double TranslationAlongX = finalParameters[0];
double TranslationAlongY = finalParameters[1];
unsigned int numberOfIterations = optimizer->GetCurrentIteration();
double bestValue = optimizer->GetValue();
// Print out results
//
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;
// Software Guide : BeginLatex
//
// Let's execute this example using the following images
//
// \begin{itemize}
// \item BrainT1SliceBorder20.png
// \item BrainProtonDensitySliceShifted13x17y.png
// \end{itemize}
//
// The output produced by the execution of the method is
//
// \begin{verbatim}
// 0 -0.316956 [11.4200, 11.2063]
// 1 -0.562048 [18.2938, 25.6545]
// 2 -0.407696 [11.3643, 21.6569]
// 3 -0.5702 [13.7244, 18.4274]
// 4 -0.803252 [11.1634, 15.3547]
//
// 0 -0.697586 [12.8778, 16.3846]
// 1 -0.901984 [13.1794, 18.3617]
// 2 -0.827423 [13.0545, 17.3695]
// 3 -0.92754 [12.8528, 16.3901]
// 4 -0.902671 [12.9426, 16.8819]
// 5 -0.941212 [13.1402, 17.3413]
//
// 0 -0.922239 [13.0364, 17.1138]
// 1 -0.930203 [12.9463, 16.8806]
// 2 -0.930959 [13.0191, 16.9822]
//
//
// Result =
// Translation X = 13.0192
// Translation Y = 16.9823
// Iterations = 4
// Metric value = -0.929237
// \end{verbatim}
//
// These values are a close match to the true misalignment of $(13,17)$
// introduced in the moving image.
//
// Software Guide : EndLatex
using ResampleFilterType =
ResampleFilterType::Pointer resample = ResampleFilterType::New();
resample->SetTransform(transform);
resample->SetInput(movingImageReader->GetOutput());
FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
resample->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
resample->SetOutputOrigin(fixedImage->GetOrigin());
resample->SetOutputSpacing(fixedImage->GetSpacing());
resample->SetOutputDirection(fixedImage->GetDirection());
resample->SetDefaultPixelValue(backgroundGrayLevel);
using OutputPixelType = unsigned char;
using OutputImageType = itk::Image<OutputPixelType, Dimension>;
using CastFilterType =
WriterType::Pointer writer = WriterType::New();
CastFilterType::Pointer caster = CastFilterType::New();
writer->SetFileName(outImagefile);
caster->SetInput(resample->GetOutput());
writer->SetInput(caster->GetOutput());
writer->Update();
//
// Generate checkerboards before and after registration
//
using CheckerBoardFilterType = itk::CheckerBoardImageFilter<FixedImageType>;
CheckerBoardFilterType::Pointer checker = CheckerBoardFilterType::New();
checker->SetInput1(fixedImage);
checker->SetInput2(resample->GetOutput());
caster->SetInput(checker->GetOutput());
writer->SetInput(caster->GetOutput());
resample->SetDefaultPixelValue(0);
// Before registration
TransformType::Pointer identityTransform = TransformType::New();
identityTransform->SetIdentity();
resample->SetTransform(identityTransform);
for (int q = 0; q < argc; ++q)
{
std::cout << q << " " << argv[q] << std::endl;
}
if (checkerBoardBefore != std::string(""))
{
writer->SetFileName(checkerBoardBefore);
writer->Update();
}
// After registration
resample->SetTransform(transform);
if (checkerBoardAfter != std::string(""))
{
writer->SetFileName(checkerBoardAfter);
writer->Update();
}
// Software Guide : BeginLatex
//
// \begin{figure}
// \center
// \includegraphics[width=0.32\textwidth]{MultiResImageRegistration1Output}
// \includegraphics[width=0.32\textwidth]{MultiResImageRegistration1CheckerboardBefore}
// \includegraphics[width=0.32\textwidth]{MultiResImageRegistration1CheckerboardAfter}
// \itkcaption[Multi-Resolution registration input images]{Mapped moving
// image (left) and composition of fixed and moving images before (center)
// and after (right) registration.}
// \label{fig:MultiResImageRegistration1Output}
// \end{figure}
//
// The result of resampling the moving image is presented in the left image
// of Figure \ref{fig:MultiResImageRegistration1Output}. The center and
// right images of the figure depict a checkerboard composite of the fixed
// and moving images before and after registration.
//
// Software Guide : EndLatex
// Software Guide : BeginLatex
//
// \begin{figure}
// \center
// \includegraphics[height=0.44\textwidth]{MultiResImageRegistration1TraceTranslations}
// \includegraphics[height=0.44\textwidth]{MultiResImageRegistration1TraceMetric}
// \itkcaption[Multi-Resolution registration output images]{Sequence of
// translations and metric values at each iteration of the optimizer.}
// \label{fig:MultiResImageRegistration1Trace}
// \end{figure}
//
// Figure \ref{fig:MultiResImageRegistration1Trace} (left) shows
// the sequence of translations followed by the optimizer as it searched
// the parameter space. The right side of the same figure shows the
// sequence of metric values computed as the optimizer searched the
// parameter space. From the trace, we can see that with the more
// aggressive optimization parameters we get quite close to the optimal
// value within 5 iterations with the remaining iterations just doing fine
// adjustments. It is interesting to compare these results with those
// of the single resolution example in Section
// \ref{sec:MultiModalityRegistrationMattes}, where 46 iterations were
// required as more conservative optimization parameters had to be used.
//
// Software Guide : EndLatex
return EXIT_SUCCESS;
}
itk::CastImageFilter
Casts input pixels to output pixel type.
Definition: itkCastImageFilter.h:104
itkRegularStepGradientDescentOptimizerv4.h
itkImageFileReader.h
itk::CheckerBoardImageFilter
Combines two images in a checkerboard pattern.
Definition: itkCheckerBoardImageFilter.h:46
itk::SmartPointer< Self >
itkCastImageFilter.h
itkImageRegistrationMethodv4.h
itkTranslationTransform.h
itk::ImageFileReader
Data source that reads image data from a single file.
Definition: itkImageFileReader.h:75
itk::RegularStepGradientDescentOptimizerv4
Regular Step Gradient descent optimizer.
Definition: itkRegularStepGradientDescentOptimizerv4.h:47
itk::Command
Superclass for callback/observer methods.
Definition: itkCommand.h:45
itkCheckerBoardImageFilter.h
itk::ImageFileWriter
Writes image data to a single file.
Definition: itkImageFileWriter.h:87
itk::Command
class ITK_FORWARD_EXPORT Command
Definition: itkObject.h:43
itk::TranslationTransform
Translation transformation of a vector space (e.g. space coordinates)
Definition: itkTranslationTransform.h:43
itk::Command::Execute
virtual void Execute(Object *caller, const EventObject &event)=0
itkImageFileWriter.h
itk::ResampleImageFilter
Resample an image via a coordinate transform.
Definition: itkResampleImageFilter.h:90
itk::Object
Base class for most ITK classes.
Definition: itkObject.h:62
itk::Image
Templated n-dimensional image class.
Definition: itkImage.h:86
itk::EventObject
Abstraction of the Events used to communicating among filters and with GUIs.
Definition: itkEventObject.h:57
itkResampleImageFilter.h
itk::GTest::TypedefsAndConstructors::Dimension2::Dimension
constexpr unsigned int Dimension
Definition: itkGTestTypedefsAndConstructors.h:44
itkCommand.h
itkMattesMutualInformationImageToImageMetricv4.h
itk::MattesMutualInformationImageToImageMetricv4
Computes the mutual information between two images to be registered using the method of Mattes et al.
Definition: itkMattesMutualInformationImageToImageMetricv4.h:103