Doxygen52/html/Examples_2RegistrationITKv4_2ModelToImageRegistration1_8cxx-example.html

/*=========================================================================

 *

 *  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 : 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>;


  itkTypeMacro(IterationCallback, Superclass);

  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 ParametersType = typename Superclass::ParametersType;

  using DerivativeType = typename Superclass::DerivativeType;

  using MeasureType = typename Superclass::MeasureType;


  itkNewMacro(Self);


  itkTypeMacro(SimpleImageToSpatialObjectMetric, ImageToSpatialObjectMetric);


  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 it : m_PointList)

    {

      PointType transformedPoint = this->m_Transform->TransformPoint(it);

      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

  EllipseType::Pointer ellipse1 = EllipseType::New();

  EllipseType::Pointer ellipse2 = EllipseType::New();

  EllipseType::Pointer 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

  GroupType::Pointer 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

  SpatialObjectToImageFilterType::Pointer 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>;

  GaussianFilterType::Pointer 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>;

  RegistrationType::Pointer 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>;

  MetricType::Pointer 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>;

  InterpolatorType::Pointer 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;

  OptimizerType::Pointer 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<>;

  TransformType::Pointer 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

  itk::Statistics::NormalVariateGenerator::Pointer 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>;

  IterationCallbackType::Pointer 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

  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;

}