DemonsRegistration1

Overview

This example illustrates how to use the classic Demons registration algorithm. The user supplied parameters for the algorithm are the number of iterations and the standard deviations for the Gaussian smoothing of the total displacement field. Additional methods which control regularization, total field smoothing for elastic model or update field smoothing for viscous model are available.

The underlying assumption of the demons framework is that the intensities of homologous points are equal. The example uses histogram matching to make the two images similar prior to registration. This is relevant for registration of MR images where the assumption is not valid. For other imaging modalities where the assumption is valid, such as CT, this step is not necessary. Additionally, the command design pattern is used to monitor registration progress. The resulting deformation field is written to file.

See also: DemonsRegistration2.

Code

// This example is based on ITK's DeformableRegistration2.cxx example


#include <SimpleITK.h>
#include <iostream>
#include <stdlib.h>
#include <iomanip>

namespace sitk = itk::simple;


class IterationUpdate : public sitk::Command
{
public:
  IterationUpdate(const sitk::DemonsRegistrationFilter & m)
    : m_Filter(m)
  {}

  void
  Execute() override
  {
    // use sitk's output operator for std::vector etc..
    using sitk::operator<<;

    // stash the stream state
    std::ios state(NULL);
    state.copyfmt(std::cout);
    std::cout << std::fixed << std::setfill(' ') << std::setprecision(5);
    std::cout << std::setw(3) << m_Filter.GetElapsedIterations();
    std::cout << " = " << std::setw(10) << m_Filter.GetMetric();
    std::cout << std::endl;

    std::cout.copyfmt(state);
  }

private:
  const sitk::DemonsRegistrationFilter & m_Filter;
};


int
main(int argc, char * argv[])
{

  if (argc < 4)
  {
    std::cerr << "Usage: " << argv[0] << " <fixedImageFilter> <movingImageFile> <outputTransformFile>" << std::endl;
    return 1;
  }

  sitk::Image fixed = sitk::ReadImage(argv[1], sitk::sitkFloat32);

  sitk::Image moving = sitk::ReadImage(argv[2], sitk::sitkFloat32);

  sitk::HistogramMatchingImageFilter matcher;
  matcher.SetNumberOfHistogramLevels(1024);
  matcher.SetNumberOfMatchPoints(7);
  matcher.ThresholdAtMeanIntensityOn();
  moving = matcher.Execute(moving, fixed);

  sitk::DemonsRegistrationFilter filter;

  IterationUpdate cmd(filter);
  filter.AddCommand(sitk::sitkIterationEvent, cmd);

  filter.SetNumberOfIterations(50);
  filter.SetStandardDeviations(1.0);

  sitk::Image displacementField = filter.Execute(fixed, moving);

  std::cout << "-------" << std::endl;
  std::cout << "Number Of Iterations: " << filter.GetElapsedIterations() << std::endl;
  std::cout << " RMS: " << filter.GetRMSChange() << std::endl;

  sitk::DisplacementFieldTransform outTx(displacementField);

  sitk::WriteTransform(outTx, argv[3]);

  return 0;
}