Image Registration Method BSpline 1

Overview

If you are not familiar with the SimpleITK registration framework we recommend that you read the registration overview before continuing with the example.

Example Run

Running the Python code with the following inputs:

main('BrainProtonDensitySliceBorder20.png', 'BrainProtonDensitySliceBSplined10.png', 'ImageRegistrationMethodBSpline1.hdf5')

produces the text and images below.

Output Text

Text Output (click triangle to collapse)
Initial Parameters:
(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
  0 =   -0.83872
  0 =   -0.85295
  0 =   -0.89317
  0 =   -0.91369
  0 =   -0.91369
  1 =   -0.95662
  1 =   -0.95662
  2 =   -0.96795
  2 =   -0.96795
  3 =   -0.97045
  3 =   -0.97045
  4 =   -0.97769
  4 =   -0.97769
  5 =   -0.98248
  5 =   -0.98248
  6 =   -0.98273
  6 =   -0.98273
  7 =   -0.98739
  7 =   -0.98739
  8 =   -0.98956
  8 =   -0.98956
  9 =   -0.99150
  9 =   -0.99150
 10 =   -0.99169
 10 =   -0.99169
 11 =   -0.99289
 11 =   -0.99289
 12 =   -0.99322
 12 =   -0.99322
 13 =   -0.99362
 13 =   -0.99362
 14 =   -0.99409
 14 =   -0.99409
 15 =   -0.99400
 15 =   -0.99420
 15 =   -0.99420
 16 =   -0.99434
 16 =   -0.99434
 17 =   -0.99444
 17 =   -0.99444
 18 =   -0.99456
 18 =   -0.99456
 19 =   -0.99462
 19 =   -0.99462
 20 =   -0.99467
 20 =   -0.99467
 21 =   -0.99469
 21 =   -0.99469
 22 =   -0.99471
 22 =   -0.99471
 23 =   -0.99473
 23 =   -0.99473
 24 =   -0.99476
 24 =   -0.99476
 25 =   -0.99478
 25 =   -0.99478
 26 =   -0.99481
 26 =   -0.99481
 27 =   -0.99484
 27 =   -0.99484
 28 =   -0.99486
 28 =   -0.99486
 29 =   -0.99487
 29 =   -0.99487
 30 =   -0.99488
 30 =   -0.99488
-------
itk::simple::BSplineTransform
 BSplineTransform (0x560a8338d5f0)
   RTTI typeinfo:   itk::BSplineTransform<double, 2u, 3u>
   Reference Count: 3
   Modified Time: 2650
   Debug: Off
   Object Name: 
   Observers: 
     none
   CoefficientImage: [ 0x560a83312870, 0x560a8330b1e0 ]
   TransformDomainOrigin: [-0.625, -0.625]
   TransformDomainPhysicalDimensions: [221.25, 257.25]
   TransformDomainDirection: 1 0
0 1

   TransformDomainMeshSize: [8, 8]
   GridSize: [11, 11]
   GridOrigin: [-28.2812, -32.7812]
   GridSpacing: [27.6562, 32.1562]
   GridDirection: 1 0
0 1


Optimizer stop condition: LBFGSBOptimizerv4: Gradient tolerance reached. Gradient tolerance is 1e-05
 Iteration: 31
 Metric value: -0.9948792040893044

Input Images

_images/ImageRegistrationMethodBSpline1_fixed.png

Fixed Image

_images/ImageRegistrationMethodBSpline1_moving.png

Moving Image

Output Image

_images/ImageRegistrationMethodBSpline1_composition.png

Composition Image

Code

// This one header will include all SimpleITK filters and external
// objects.
#include <SimpleITK.h>
#include <iostream>
#include <stdlib.h>
#include <iomanip>

namespace sitk = itk::simple;

// use sitk's output operator for std::vector etc..
using sitk::operator<<;


class IterationUpdate : public sitk::Command
{
public:
  IterationUpdate(const sitk::ImageRegistrationMethod & m)
    : m_Method(m)
  {}

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

    if (m_Method.GetOptimizerIteration() == 0)
    {
      std::cout << m_Method.ToString() << std::endl;
    }

    // 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_Method.GetOptimizerIteration();
    std::cout << " = " << std::setw(10) << m_Method.GetMetricValue() << std::endl;
    std::cout.copyfmt(state);
  }

private:
  const sitk::ImageRegistrationMethod & m_Method;
};


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);

  std::vector<unsigned int> transformDomainMeshSize(fixed.GetDimension(), 8);

  sitk::BSplineTransform tx = sitk::BSplineTransformInitializer(fixed, transformDomainMeshSize);

  std::cout << "Initial Parameters:" << tx.GetParameters() << std::endl;


  sitk::ImageRegistrationMethod R;
  R.SetMetricAsCorrelation();

  const double       gradientConvergenceTolerance = 1e-5;
  const unsigned int maximumNumberOfIterations = 100;
  const unsigned int maximumNumberOfCorrections = 5;
  const unsigned int maximumNumberOfFunctionEvaluations = 1000;
  const double       costFunctionConvergenceFactor = 1e+7;
  R.SetOptimizerAsLBFGSB(gradientConvergenceTolerance,
                         maximumNumberOfIterations,
                         maximumNumberOfCorrections,
                         maximumNumberOfFunctionEvaluations,
                         costFunctionConvergenceFactor);
  R.SetInitialTransform(tx, true);
  R.SetInterpolator(sitk::sitkLinear);

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

  sitk::Transform outTx = R.Execute(fixed, moving);

  std::cout << "-------" << std::endl;
  std::cout << outTx.ToString() << std::endl;
  std::cout << "Optimizer stop condition: " << R.GetOptimizerStopConditionDescription() << std::endl;
  std::cout << " Iteration: " << R.GetOptimizerIteration() << std::endl;
  std::cout << " Metric value: " << R.GetMetricValue() << std::endl;

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

  return 0;
}