Examples/RegistrationITKv4/ModelToImageRegistration1.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
 *
 *         https://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 : BeginLatex
//
//  This example illustrates the use of the \doxygen{SpatialObject} as a
//  component of the registration framework in order to perform model based
//  registration. The current example creates a geometrical model composed of
//  several ellipses. Then, it uses the model to produce a synthetic binary
//  image of the ellipses. Next, it introduces perturbations on the position
//  and shape of the model, and finally it uses the perturbed version as the
//  input to a registration problem. A metric is defined to evaluate the
//  fitness between the geometric model and the image.
//
// Software Guide : EndLatex
 
 
//  Software Guide : BeginLatex
//
//  Let's look first at the classes required to support
//  SpatialObject. In this example we use the
//  \doxygen{EllipseSpatialObject} as the basic shape components and we use
//  the \doxygen{GroupSpatialObject} to group them together as a
//  representation of a more complex shape. Their respective headers are
//  included below.
//
//  \index{itk::EllipseSpatialObject!header}
//  \index{itk::GroupSpatialObject!header}
//
//  Software Guide : EndLatex
 
//  Software Guide : BeginCodeSnippet
#include "itkEllipseSpatialObject.h"
#include "itkGroupSpatialObject.h"
//  Software Guide : EndCodeSnippet
 
 
//  Software Guide : BeginLatex
//
//  In order to generate the initial synthetic image of the ellipses, we use
//  the \doxygen{SpatialObjectToImageFilter} that tests---for every pixel in
//  the image---whether the pixel (and hence the spatial object) is
//  \emph{inside} or \emph{outside} the geometric model.
//
//  \index{itk::Spatial\-Object\-To\-Image\-Filter!header}
//
//  Software Guide : EndLatex
 
//  Software Guide : BeginCodeSnippet
#include "itkSpatialObjectToImageFilter.h"
//  Software Guide : EndCodeSnippet
 
 
#include "itkImageToSpatialObjectRegistrationMethod.h"
 
 
//  Software Guide : BeginLatex
//
//  A metric is defined to evaluate the fitness between the
//  SpatialObject and the Image. The base class for this
//  type of metric is the \doxygen{ImageToSpatialObjectMetric}, whose header
//  is included below.
//
//  \index{itk::Image\-To\-Spatial\-Object\-Metric!header}
//
//  Software Guide : EndLatex
 
//  Software Guide : BeginCodeSnippet
#include "itkImageToSpatialObjectMetric.h"
//  Software Guide : EndCodeSnippet
 
 
//  Software Guide : BeginLatex
//
//  As in previous registration problems, we have to evaluate the image
//  intensity in non-grid positions. The
//  \doxygen{LinearInterpolateImageFunction} is used here for this purpose.
//
//  \index{itk::Linear\-Interpolate\-Image\-Function!header}
//
//  Software Guide : EndLatex
 
//  Software Guide : BeginCodeSnippet
#include "itkLinearInterpolateImageFunction.h"
//  Software Guide : EndCodeSnippet
 
 
//  Software Guide : BeginLatex
//
//  The SpatialObject is mapped from its own space into the image
//  space by using a \doxygen{Transform}. In this
//  example, we use the \doxygen{Euler2DTransform}.
//
//  Software Guide : EndLatex
 
//  Software Guide : BeginCodeSnippet
#include "itkEuler2DTransform.h"
//  Software Guide : EndCodeSnippet
 
 
//  Software Guide : BeginLatex
//
//  Registration is fundamentally an optimization problem. Here we include
//  the optimizer used to search the parameter space and identify the best
//  transformation that will map the shape model on top of the image. The
//  optimizer used in this example is the
//  \doxygen{OnePlusOneEvolutionaryOptimizer} that implements an
//  \href{http://www.aic.nrl.navy.mil/galist/}{evolutionary algorithm}.
//
//  Software Guide : EndLatex
 
//  Software Guide : BeginCodeSnippet
#include "itkOnePlusOneEvolutionaryOptimizer.h"
//  Software Guide : EndCodeSnippet
 
 
#include "itkDiscreteGaussianImageFilter.h"
#include "itkNormalVariateGenerator.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkCastImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
 
 
//  Software Guide : BeginLatex
//
//  As in previous registration examples, it is important to
//  track the evolution of the optimizer as it progresses through the
//  parameter space.  This is done by using the Command/Observer paradigm. The
//  following lines of code implement the \doxygen{Command} observer that
//  monitors the progress of the registration. The code is quite
//  similar to what we have used in previous registration examples.
//
//  \index{Model to Image Registration!Observer}
//
//  Software Guide : EndLatex
 
//  Software Guide : BeginCodeSnippet
#include "itkCommand.h"
template <class TOptimizer>
class IterationCallback : public itk::Command
{
public:
  using Self = IterationCallback;
  using Superclass = itk::Command;
  using Pointer = itk::SmartPointer<Self>;
  using ConstPointer = itk::SmartPointer<const Self>;
 
  itkOverrideGetNameOfClassMacro(IterationCallback);
  itkNewMacro(Self);
 
  using OptimizerType = TOptimizer;
 
  void
  SetOptimizer(OptimizerType * optimizer)
  {
    m_Optimizer = optimizer;
    m_Optimizer->AddObserver(itk::IterationEvent(), this);
  }
  void
  Execute(itk::Object * caller, const itk::EventObject & event) override
  {
    Execute((const itk::Object *)caller, event);
  }
 
  void
  Execute(const itk::Object *, const itk::EventObject & event) override
  {
    if (typeid(event) == typeid(itk::StartEvent))
    {
      std::cout << std::endl << "Position              Value";
      std::cout << std::endl << std::endl;
    }
    else if (typeid(event) == typeid(itk::IterationEvent))
    {
      std::cout << m_Optimizer->GetCurrentIteration() << "   ";
      std::cout << m_Optimizer->GetValue() << "   ";
      std::cout << m_Optimizer->GetCurrentPosition() << std::endl;
    }
    else if (typeid(event) == typeid(itk::EndEvent))
    {
      std::cout << std::endl << std::endl;
      std::cout << "After " << m_Optimizer->GetCurrentIteration();
      std::cout << "  iterations " << std::endl;
      std::cout << "Solution is    = " << m_Optimizer->GetCurrentPosition();
      std::cout << std::endl;
    }
  }
  //  Software Guide : EndCodeSnippet
 
protected:
  IterationCallback() = default;
  itk::WeakPointer<OptimizerType> m_Optimizer;
};
 
//  Software Guide : BeginLatex
//
//  This command will be invoked at every iteration of the optimizer and will
//  print out the current combination of transform parameters.
//
//  Software Guide : EndLatex
 
 
//  Software Guide : BeginLatex
//
//  Consider now the most critical component of this new registration
//  approach: the metric.  This component evaluates the match between the
//  SpatialObject and the Image. The
//  smoothness and regularity of the metric determine the difficulty of the
//  task assigned to the optimizer. In this case, we use a very robust
//  optimizer that should be able to find its way even in the most
//  discontinuous cost functions. The metric to be implemented should derive
//  from the ImageToSpatialObjectMetric class.
//
//  The following code implements a simple metric that computes the sum of
//  the pixels that are inside the spatial object. In fact, the metric
//  maximum is obtained when the model and the image are aligned. The metric
//  is templated over the type of the SpatialObject and the type of
//  the Image.
//
//  Software Guide : EndLatex
 
//  Software Guide : BeginCodeSnippet
template <typename TFixedImage, typename TMovingSpatialObject>
class SimpleImageToSpatialObjectMetric
  : public itk::ImageToSpatialObjectMetric<TFixedImage, TMovingSpatialObject>
{
  //  Software Guide : EndCodeSnippet
 
public:
  using Self = SimpleImageToSpatialObjectMetric;
  using Superclass =
    itk::ImageToSpatialObjectMetric<TFixedImage, TMovingSpatialObject>;
  using Pointer = itk::SmartPointer<Self>;
  using ConstPointer = itk::SmartPointer<const Self>;
 
  using PointType = itk::Point<double, 2>;
  using PointListType = std::list<PointType>;
  using MovingSpatialObjectType = TMovingSpatialObject;
  using typename Superclass::ParametersType;
  using typename Superclass::DerivativeType;
  using typename Superclass::MeasureType;
 
  itkNewMacro(Self);
 
  itkOverrideGetNameOfClassMacro(SimpleImageToSpatialObjectMetric);
 
  static constexpr unsigned int ParametricSpaceDimension = 3;
 
  void
  SetMovingSpatialObject(const MovingSpatialObjectType * object) override
  {
    if (!this->m_FixedImage)
    {
      std::cout << "Please set the image before the moving spatial object"
                << std::endl;
      return;
    }
    this->m_MovingSpatialObject = object;
    m_PointList.clear();
    using myIteratorType =
      itk::ImageRegionConstIteratorWithIndex<TFixedImage>;
    myIteratorType it(this->m_FixedImage,
                      this->m_FixedImage->GetBufferedRegion());
 
    itk::Point<double, 2> point;
 
    while (!it.IsAtEnd())
    {
      this->m_FixedImage->TransformIndexToPhysicalPoint(it.GetIndex(), point);
 
      if (this->m_MovingSpatialObject->IsInsideInWorldSpace(point, 99999))
      {
        m_PointList.push_back(point);
      }
      ++it;
    }
 
    std::cout << "Number of points in the metric = "
              << static_cast<unsigned long>(m_PointList.size()) << std::endl;
  }
 
  void
  GetDerivative(const ParametersType &, DerivativeType &) const override
  {
    return;
  }
 
  //  Software Guide : BeginLatex
  //
  //  The fundamental operation of the metric is its \code{GetValue()} method.
  //  It is in this method that the fitness value is computed. In our current
  //  example, the fitness is computed over the points of the
  //  SpatialObject. For each point, its coordinates are mapped
  //  through the transform into image space. The resulting point is used
  //  to evaluate the image and the resulting value is accumulated in a sum.
  //  Since we are not allowing scale changes, the optimal value of the sum
  //  will result when all the SpatialObject points are mapped on
  //  the white regions of the image. Note that the argument for the
  //  \code{GetValue()} method is the array of parameters of the transform.
  //
  //  \index{Image\-To\-Spatial\-Object\-Metric!GetValue()}
  //
  //  Software Guide : EndLatex
 
  //  Software Guide : BeginCodeSnippet
  MeasureType
  GetValue(const ParametersType & parameters) const override
  {
    double value;
    this->m_Transform->SetParameters(parameters);
    value = 0;
    for (auto point : m_PointList)
    {
      const PointType transformedPoint =
        this->m_Transform->TransformPoint(point);
      if (this->m_Interpolator->IsInsideBuffer(transformedPoint))
      {
        value += this->m_Interpolator->Evaluate(transformedPoint);
      }
    }
    return value;
  }
  //  Software Guide : EndCodeSnippet
 
  void
  GetValueAndDerivative(const ParametersType & parameters,
                        MeasureType &          Value,
                        DerivativeType &       Derivative) const override
  {
    Value = this->GetValue(parameters);
    this->GetDerivative(parameters, Derivative);
  }
private:
  PointListType m_PointList;
};
 
 
//  Software Guide : BeginLatex
//
//  Having defined all the registration components we are ready to put the
//  pieces together and implement the registration process.
//
//  Software Guide : EndLatex
 
 
int
main(int argc, char * argv[])
{
  if (argc > 1)
  {
    std::cerr << "Too many parameters " << std::endl;
    std::cerr << "Usage: " << argv[0] << std::endl;
  }
 
  //  Software Guide : BeginLatex
  //
  //  First we instantiate the GroupSpatialObject and
  //  EllipseSpatialObject. These two objects are parameterized by
  //  the dimension of the space. In our current example a $2D$ instantiation
  //  is created.
  //
  //  \index{Group\-Spatial\-Object!Instantiation}
  //  \index{Ellipse\-Spatial\-Object!Instantiation}
  //
  //  Software Guide : EndLatex
 
  //  Software Guide : BeginCodeSnippet
  using GroupType = itk::GroupSpatialObject<2>;
  using EllipseType = itk::EllipseSpatialObject<2>;
  //  Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  The image is instantiated in the following lines using the pixel
  //  type and the space dimension. This image uses a \code{float} pixel
  //  type since we plan to blur it in order to increase the capture radius of
  //  the optimizer. Images of real pixel type behave better under blurring
  //  than those of integer pixel type.
  //
  //  Software Guide : EndLatex
 
  //  Software Guide : BeginCodeSnippet
  using ImageType = itk::Image<float, 2>;
  //  Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  Here is where the fun begins! In the following lines we create the
  //  EllipseSpatialObjects using their \code{New()} methods, and
  //  assigning the results to SmartPointers. These lines will create
  //  three ellipses.
  //
  //  Software Guide : EndLatex
 
  //  Software Guide : BeginCodeSnippet
  auto ellipse1 = EllipseType::New();
  auto ellipse2 = EllipseType::New();
  auto ellipse3 = EllipseType::New();
  //  Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  Every class deriving from SpatialObject has particular
  //  parameters enabling the user to tailor its shape. In the case of the
  //  EllipseSpatialObject, \code{SetRadius()} is used to
  //  define the ellipse size. An additional \code{SetRadius(Array)} method
  //  allows the user to define the ellipse axes independently.
  //
  //  \index{itk::EllipseSpatialObject!SetRadius()}
  //
  //  Software Guide : EndLatex
 
  //  Software Guide : BeginCodeSnippet
  ellipse1->SetRadiusInObjectSpace(10.0);
  ellipse2->SetRadiusInObjectSpace(10.0);
  ellipse3->SetRadiusInObjectSpace(10.0);
  //  Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  The ellipses are created centered in space by default. We use the
  //  following lines of code to arrange the ellipses in a triangle.
  //  The spatial transform intrinsically associated with the object is
  //  accessed by the \code{GetTransform()} method. This transform can define
  //  a translation in space with the \code{SetOffset()} method.  We take
  //  advantage of this feature to place the ellipses at particular
  //  points in space.
  //
  //  Software Guide : EndLatex
 
  // Place each ellipse at the right position to form a triangle
 
  //  Software Guide : BeginCodeSnippet
  EllipseType::TransformType::OffsetType offset;
  offset[0] = 100.0;
  offset[1] = 40.0;
 
  ellipse1->GetModifiableObjectToParentTransform()->SetOffset(offset);
  ellipse1->Update();
 
  offset[0] = 40.0;
  offset[1] = 150.0;
  ellipse2->GetModifiableObjectToParentTransform()->SetOffset(offset);
  ellipse2->Update();
 
  offset[0] = 150.0;
  offset[1] = 150.0;
  ellipse3->GetModifiableObjectToParentTransform()->SetOffset(offset);
  ellipse3->Update();
  //  Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  Note that after a change has been made in the transform, the
  //  SpatialObject invokes the method
  //  \code{ComputeGlobalTransform()} in order to update its global
  //  transform.  The reason for doing this is that SpatialObjects
  //  can be arranged in hierarchies. It is then possible to change the
  //  position of a set of spatial objects by moving the parent of the group.
  //
  //  Software Guide : EndLatex
 
  //  Software Guide : BeginLatex
  //
  //  Now we add the three EllipseSpatialObjects to a
  //  GroupSpatialObject that will be subsequently passed on to the
  //  registration method. The GroupSpatialObject facilitates the
  //  management of the three ellipses as a higher level structure
  //  representing a complex shape. Groups can be nested any number of levels
  //  in order to represent shapes with higher detail.
  //
  //  \index{itk::GroupSpatialObject!New()}
  //  \index{itk::GroupSpatialObject!Pointer}
  //
  //  Software Guide : EndLatex
 
  //  Software Guide : BeginCodeSnippet
  auto group = GroupType::New();
  group->AddChild(ellipse1);
  group->AddChild(ellipse2);
  group->AddChild(ellipse3);
  //  Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  Having the geometric model ready, we proceed to generate the binary
  //  image representing the imprint of the space occupied by the ellipses.
  //  The SpatialObjectToImageFilter is used to that end. Note that
  //  this filter is instantiated over the spatial object used and the image
  //  type to be generated.
  //
  //
  //  \index{itk::Spatial\-Object\-To\-Image\-Filter!Instantiation}
  //
  //  Software Guide : EndLatex
 
  //  Software Guide : BeginCodeSnippet
  using SpatialObjectToImageFilterType =
    itk::SpatialObjectToImageFilter<GroupType, ImageType>;
  //  Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  With the defined type, we construct a filter using the \code{New()}
  //  method. The newly created filter is assigned to a SmartPointer.
  //
  //  \index{itk::SpatialObjectToImageFilter!New()}
  //  \index{itk::SpatialObjectToImageFilter!Pointer}
  //
  //  Software Guide : EndLatex
 
  //  Software Guide : BeginCodeSnippet
  auto imageFilter = SpatialObjectToImageFilterType::New();
  //  Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  The GroupSpatialObject is passed as input to the filter.
  //
  //  \index{itk::SpatialObjectToImageFilter!SetInput()}
  //
  //  Software Guide : EndLatex
 
  //  Software Guide : BeginCodeSnippet
  imageFilter->SetInput(group);
  //  Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  The \doxygen{SpatialObjectToImageFilter} acts as a resampling filter.
  //  Therefore it requires the user to define the size of the desired output
  //  image. This is specified with the \code{SetSize()} method.
  //
  //  \index{itk::SpatialObjectToImageFilter!SetSize()}
  //
  //  Software Guide : EndLatex
 
  //  Software Guide : BeginCodeSnippet
  ImageType::SizeType size;
  size[0] = 200;
  size[1] = 200;
  imageFilter->SetSize(size);
  //  Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  Finally we trigger the execution of the filter by calling the
  //  \code{Update()} method.
  //
  //  \index{itk::SpatialObjectToImageFilter!Update()}
  //
  //  Software Guide : EndLatex
 
  //  Software Guide : BeginCodeSnippet
  imageFilter->Update();
  //  Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  In order to obtain a smoother metric, we blur the image using a
  //  \doxygen{DiscreteGaussianImageFilter}. This extends the capture radius
  //  of the metric and produce a more continuous cost function to
  //  optimize. The following lines instantiate the Gaussian filter and
  //  create one object of this type using the \code{New()} method.
  //
  //  Software Guide : EndLatex
 
  //  Software Guide : BeginCodeSnippet
  using GaussianFilterType =
    itk::DiscreteGaussianImageFilter<ImageType, ImageType>;
  auto gaussianFilter = GaussianFilterType::New();
  //  Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  The output of the SpatialObjectToImageFilter is connected as
  //  input to the DiscreteGaussianImageFilter.
  //
  //  Software Guide : EndLatex
 
  //  Software Guide : BeginCodeSnippet
  gaussianFilter->SetInput(imageFilter->GetOutput());
  //  Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  The variance of the filter is defined as a large value in order to
  //  increase the capture radius. Finally the execution of the filter is
  //  triggered using the \code{Update()} method.
  //
  //  Software Guide : EndLatex
 
  //  Software Guide : BeginCodeSnippet
  constexpr double variance = 20;
  gaussianFilter->SetVariance(variance);
  gaussianFilter->Update();
  //  Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  Below we instantiate the type of the
  //  \doxygen{ImageToSpatialObjectRegistrationMethod} method and instantiate
  //  a registration object with the \code{New()} method. Note that the
  //  registration type is templated over the Image and the
  //  SpatialObject types. The spatial object in this case is the
  //  group of spatial objects.
  //
  //  \index{itk::Image\-To\-Spatial\-Object\-Registration\-Method!Instantiation}
  //  \index{itk::Image\-To\-Spatial\-Object\-Registration\-Method!New()}
  //  \index{itk::Image\-To\-Spatial\-Object\-Registration\-Method!Pointer}
  //
  //  Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
  using RegistrationType =
    itk::ImageToSpatialObjectRegistrationMethod<ImageType, GroupType>;
  auto registration = RegistrationType::New();
  // Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  Now we instantiate the metric that is templated over
  //  the image type and the spatial object type. As usual, the \code{New()}
  //  method is used to create an object.
  //
  //  \index{itk::Image\-To\-Spatial\-Object\-Metric!Instantiation}
  //
  //  Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
  using MetricType = SimpleImageToSpatialObjectMetric<ImageType, GroupType>;
  auto metric = MetricType::New();
  // Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  An interpolator will be needed to evaluate the image at non-grid
  //  positions. Here we instantiate a linear interpolator type.
  //
  //  Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
  using InterpolatorType =
    itk::LinearInterpolateImageFunction<ImageType, double>;
  auto interpolator = InterpolatorType::New();
  // Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  The following lines instantiate the evolutionary optimizer.
  //
  //  \index{itk::One\-Plus\-One\-Evolutionary\-Optimizer!Instantiation}
  //
  //  Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
  using OptimizerType = itk::OnePlusOneEvolutionaryOptimizer;
  auto optimizer = OptimizerType::New();
  // Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  Next, we instantiate the transform class. In this case we use the
  //  Euler2DTransform that implements a rigid transform in $2D$
  //  space.
  //
  //  Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
  using TransformType = itk::Euler2DTransform<>;
  auto transform = TransformType::New();
  // Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  Evolutionary algorithms are based on testing random variations
  //  of parameters. In order to support the computation of random values,
  //  ITK provides a family of random number generators. In this example, we
  //  use the \doxygen{NormalVariateGenerator} which generates values with a
  //  normal distribution.
  //
  //  \index{itk::NormalVariateGenerator!New()}
  //  \index{itk::NormalVariateGenerator!Pointer}
  //
  //  Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
  auto generator = itk::Statistics::NormalVariateGenerator::New();
  // Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  The random number generator must be initialized with a seed.
  //
  //  \index{itk::NormalVariateGenerator!Initialize()}
  //
  //  Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
  generator->Initialize(12345);
  // Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  The OnePlusOneEvolutionaryOptimizer is initialized by
  //  specifying the random number generator, the number of samples for the
  //  initial population and the maximum number of iterations.
  //
  //  Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
  optimizer->SetNormalVariateGenerator(generator);
  optimizer->Initialize(10);
  optimizer->SetMaximumIteration(400);
  // Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  As in previous registration examples, we take care to normalize the
  //  dynamic range of the different transform parameters. In particular, the
  //  we must compensate for the ranges of the angle and translations of the
  //  Euler2DTransform. In order to achieve this goal, we provide an array of
  //  scales to the optimizer.
  //
  //  Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
  TransformType::ParametersType parametersScale;
  parametersScale.set_size(3);
  parametersScale[0] = 1000; // angle scale
 
  for (unsigned int i = 1; i < 3; ++i)
  {
    parametersScale[i] = 2; // offset scale
  }
  optimizer->SetScales(parametersScale);
  // Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  Here we instantiate the Command object that will act as an
  //  observer of the registration method and print out parameters at each
  //  iteration. Earlier, we defined this command as a class templated over
  //  the optimizer type. Once it is created with the \code{New()} method, we
  //  connect the optimizer to the command.
  //
  //  Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
  using IterationCallbackType = IterationCallback<OptimizerType>;
  auto callback = IterationCallbackType::New();
  callback->SetOptimizer(optimizer);
  // Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  All the components are plugged into the
  //  ImageToSpatialObjectRegistrationMethod object. The typical
  //  \code{Set()} methods are used here. Note the use of the
  //  \code{SetMovingSpatialObject()} method for connecting the spatial
  //  object. We provide the blurred version of the original synthetic binary
  //  image as the input image.
  //
  //  \index{itk::Image\-To\-Spatial\-Object\-Registration\-Method!SetFixedImage()}
  //  \index{itk::Image\-To\-Spatial\-Object\-Registration\-Method!SetMovingSpatialObject()}
  //  \index{itk::Image\-To\-Spatial\-Object\-Registration\-Method!SetTransform()}
  //  \index{itk::Image\-To\-Spatial\-Object\-Registration\-Method!SetInterpolator()}
  //  \index{itk::Image\-To\-Spatial\-Object\-Registration\-Method!SetOptimizer()}
  //  \index{itk::Image\-To\-Spatial\-Object\-Registration\-Method!SetMetric()}
  //
  //  Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
  registration->SetFixedImage(gaussianFilter->GetOutput());
  registration->SetMovingSpatialObject(group);
  registration->SetTransform(transform);
  registration->SetInterpolator(interpolator);
  registration->SetOptimizer(optimizer);
  registration->SetMetric(metric);
  // Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  The initial set of transform parameters is passed to the registration
  //  method using the \code{SetInitialTransformParameters()} method. Note
  //  that since our original model is already registered with the synthetic
  //  image, we introduce an artificial mis-registration in order to
  //  initialize the optimization at some point away from the optimal value.
  //
  //  Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
  TransformType::ParametersType initialParameters(
    transform->GetNumberOfParameters());
 
  initialParameters[0] = 0.2; // Angle
  initialParameters[1] = 7.0; // Offset X
  initialParameters[2] = 6.0; // Offset Y
  registration->SetInitialTransformParameters(initialParameters);
  // Software Guide : EndCodeSnippet
 
  std::cout << "Initial Parameters  : " << initialParameters << std::endl;
 
  //  Software Guide : BeginLatex
  //
  //  Due to the character of the metric used to evaluate the fitness
  //  between the spatial object and the image, we must tell the optimizer
  //  that we are interested in finding the maximum value of the metric. Some
  //  metrics associate low numeric values with good matching, while others
  //  associate high numeric values with good matching. The
  //  \code{MaximizeOn()} and \code{MaximizeOff()} methods allow the user to
  //  deal with both types of metrics.
  //
  //  \index{itk::Optimizer!MaximizeOn()}
  //  \index{itk::Optimizer!MaximizeOff()}
  //
  //  Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
  optimizer->MaximizeOn();
  // Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  Finally, we trigger the execution of the registration process with the
  //  \code{Update()} method. We place this call in a
  //  \code{try/catch} block in case any exception is thrown during the
  //  process.
  //
  //  \index{itk::Image\-To\-Spatial\-Object\-Registration\-Method!Update()}
  //
  //  Software Guide : EndLatex
 
 
  // Software Guide : BeginCodeSnippet
  try
  {
    registration->Update();
    std::cout << "Optimizer stop condition: "
              << registration->GetOptimizer()->GetStopConditionDescription()
              << std::endl;
  }
  catch (const itk::ExceptionObject & exp)
  {
    std::cerr << "Exception caught ! " << std::endl;
    std::cerr << exp << std::endl;
  }
  // Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  The set of transform parameters resulting from the registration can be
  //  recovered with the \code{GetLastTransformParameters()} method. This
  //  method returns the array of transform parameters that should be
  //  interpreted according to the implementation of each transform. In our
  //  current example, the Euler2DTransform has three parameters:
  //  the rotation angle, the translation in $x$ and the translation in $y$.
  //
  //  \index{itk::Image\-To\-Spatial\-Object\-Registration\-Method!Update()}
  //
  //  Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
  const RegistrationType::ParametersType finalParameters =
    registration->GetLastTransformParameters();
 
  std::cout << "Final Solution is : " << finalParameters << std::endl;
  // Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  // \begin{figure}
  // \center
  // \includegraphics[height=0.44\textwidth]{ModelToImageRegistrationTraceAngle}
  // \includegraphics[height=0.44\textwidth]{ModelToImageRegistrationTraceTranslations}
  // \itkcaption[SpatialObject to Image Registration results]{Plots of the
  // angle and translation parameters for a registration process between an
  // spatial object and an image.}
  // \label{fig:ModelToImageRegistrationPlots}
  // \end{figure}
  //
  //  The results are presented in
  //  Figure~\ref{fig:ModelToImageRegistrationPlots}. The left side shows the
  //  evolution of the angle parameter as a function of iteration
  //  numbers, while the right side shows the $(x,y)$ translation.
  //
  //  Software Guide : EndLatex
 
  return EXIT_SUCCESS;
}