Image Registration Method 1

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

Overview

Code

using System;
using itk.simple;


namespace itk.simple.examples {


  class IterationUpdate : Command {

    private ImageRegistrationMethod m_Method;

    public IterationUpdate(ImageRegistrationMethod m) {
      m_Method=m;
    }

    public override void Execute() {
      VectorDouble pos = m_Method.GetOptimizerPosition();
      Console.WriteLine("{0:3} = {1:10.5} : [{2}, {3}]",
                        m_Method.GetOptimizerIteration(),
                        m_Method.GetMetricValue(),
                        pos[0], pos[1]);
    }
  }

  class ImageRegistrationMethod1 {

    static void Main(string[] args) {

      if ( args.Length < 3 )
        {
          Console.WriteLine("Usage: %s <fixedImageFilter> <movingImageFile> <outputTransformFile>\n", "ImageRegistrationMethod1");
         return;
        }

      ImageFileReader reader = new ImageFileReader();
      reader.SetOutputPixelType( PixelIDValueEnum.sitkFloat32 );

      reader.SetFileName(args[0]);
      Image fixedImage = reader.Execute();

      reader.SetFileName(args[1]);
      Image movingImage = reader.Execute();

      ImageRegistrationMethod R = new  ImageRegistrationMethod();
      R.SetMetricAsMeanSquares();
      double maxStep = 4.0;
      double minStep = 0.01;
      uint numberOfIterations = 200;
      double relaxationFactor = 0.5;
      R.SetOptimizerAsRegularStepGradientDescent( maxStep,
                                                  minStep,
                                                  numberOfIterations,
                                                  relaxationFactor );
      R.SetInitialTransform( new TranslationTransform( fixedImage.GetDimension() ) );
      R.SetInterpolator( InterpolatorEnum.sitkLinear );

      IterationUpdate cmd = new IterationUpdate(R);
      R.AddCommand(EventEnum.sitkIterationEvent, cmd);

      Transform outTx = R.Execute( fixedImage, movingImage );

      // System.out.println("-------");
      // System.out.println(outTx.toString());
      // System.out.format("Optimizer stop condition: %s\n", R.getOptimizerStopConditionDescription());
      // System.out.format(" Iteration: %d\n", R.getOptimizerIteration());
      // System.out.format(" Metric value: %f\n", R.getMetricValue());

      outTx.WriteTransform(args[2]);

    }

  }

}
// 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;



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

  virtual void Execute( )
    {
      // 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_Method.GetOptimizerIteration();
      std::cout << " = " << std::setw(10) << m_Method.GetMetricValue();
      std::cout << " : " << m_Method.GetOptimizerPosition() << 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 );

  sitk::ImageRegistrationMethod R;
  R.SetMetricAsMeanSquares();
  const double maxStep = 4.0;
  const double minStep = 0.01;
  const unsigned int numberOfIterations = 200;
  const double relaxationFactor = 0.5;
  R.SetOptimizerAsRegularStepGradientDescent( maxStep,
                                              minStep,
                                              numberOfIterations,
                                              relaxationFactor );
  R.SetInitialTransform( sitk::TranslationTransform( fixed.GetDimension() ) );
  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;
}
import org.itk.simple.*;


class IterationUpdate  extends Command {

  private ImageRegistrationMethod m_Method;

  public IterationUpdate(ImageRegistrationMethod m) {
    super();
    m_Method=m;
  }

  public void execute() {
    org.itk.simple.VectorDouble pos = m_Method.getOptimizerPosition();
    System.out.format("%3d = %10.5f : [%f, %f]\n",
                      m_Method.getOptimizerIteration(),
                      m_Method.getMetricValue(),
                      pos.get(0), pos.get(1));
  }
}

class ImageRegistrationMethod1 {

  public static void main(String argv[]) {

  if ( argv.length < 3 )
    {
      System.out.format( "Usage: %s <fixedImageFilter> <movingImageFile> <outputTransformFile>\n", "ImageRegistrationMethod1");
      System.exit(-1);
    }

  org.itk.simple.ImageFileReader reader = new org.itk.simple.ImageFileReader();
  reader.setOutputPixelType( PixelIDValueEnum.sitkFloat32 );

  reader.setFileName(argv[0]);
  Image fixed = reader.execute();

  reader.setFileName(argv[1]);
  Image moving = reader.execute();

  org.itk.simple.ImageRegistrationMethod R = new  org.itk.simple.ImageRegistrationMethod();
  R.setMetricAsMeanSquares();
  double maxStep = 4.0;
  double minStep = 0.01;
  int numberOfIterations = 200;
  double relaxationFactor = 0.5;
  R.setOptimizerAsRegularStepGradientDescent( maxStep,
                                              minStep,
                                              numberOfIterations,
                                              relaxationFactor );
  R.setInitialTransform( new org.itk.simple.TranslationTransform( fixed.getDimension() ) );
  R.setInterpolator( InterpolatorEnum.sitkLinear );

  IterationUpdate cmd = new IterationUpdate(R);
  R.addCommand( EventEnum.sitkIterationEvent, cmd);

  org.itk.simple.Transform outTx = R.execute( fixed, moving );

  System.out.println("-------");
  System.out.println(outTx.toString());
  System.out.format("Optimizer stop condition: %s\n", R.getOptimizerStopConditionDescription());
  System.out.format(" Iteration: %d\n", R.getOptimizerIteration());
  System.out.format(" Metric value: %f\n", R.getMetricValue());

  outTx.writeTransform(argv[2]);

  }

}
require "SimpleITK"

function command_iteration(method)
    print(method)
    pos = method:GetOptimizedPosition()
    print( string.format("%3d = %f : (%f, %f)", method:GetOptimizerIteration(),
        method:GetMetricValue(), pos[1], pos[2] ) )
end

if #arg < 3 then
    print( string.format("Usage: %s <fixedImageFilter> <movingImageFile> <outputTransformFile>", arg[0]) )
    os.exit ( 1 )
end

fixed = SimpleITK.ReadImage( arg[1], SimpleITK.sitkFloat32 )

moving = SimpleITK.ReadImage( arg[2], SimpleITK.sitkFloat32 )


R = SimpleITK.ImageRegistrationMethod()
R:SetMetricAsMeanSquares()
R:SetOptimizerAsRegularStepGradientDescent( 4.0, .01, 200 )
R:SetInitialTransform( SimpleITK.Transform( fixed:GetDimension(), SimpleITK.sitkTranslation ) )
R:SetInterpolator( SimpleITK.sitkLinear )

-- callback for progress reporting doesn't work yet in Lua
-- R:AddCommand( SimpleITK.sitkIterationEvent, command_iteration(R) )

outTx = R:Execute(fixed, moving)

print("-------")
print(outTx)
print(string.format("Optimizer stop condition: %s", R:GetOptimizerStopConditionDescription()) )
print(" Iteration: ", R:GetOptimizerIteration() )
print(" Metric value: ", R:GetMetricValue() )

SimpleITK.WriteTransform(outTx,  arg[3])
#!/usr/bin/env python

from __future__ import print_function

import SimpleITK as sitk
import sys
import os


def command_iteration(method) :
    print("{0:3} = {1:10.5f} : {2}".format(method.GetOptimizerIteration(),
                                   method.GetMetricValue(),
                                   method.GetOptimizerPosition()))

if len ( sys.argv ) < 4:
    print( "Usage: {0} <fixedImageFilter> <movingImageFile> <outputTransformFile>".format(sys.argv[0]))
    sys.exit ( 1 )


fixed = sitk.ReadImage(sys.argv[1], sitk.sitkFloat32)

moving = sitk.ReadImage(sys.argv[2], sitk.sitkFloat32)

R = sitk.ImageRegistrationMethod()
R.SetMetricAsMeanSquares()
R.SetOptimizerAsRegularStepGradientDescent(4.0, .01, 200 )
R.SetInitialTransform(sitk.TranslationTransform(fixed.GetDimension()))
R.SetInterpolator(sitk.sitkLinear)

R.AddCommand( sitk.sitkIterationEvent, lambda: command_iteration(R) )

outTx = R.Execute(fixed, moving)

print("-------")
print(outTx)
print("Optimizer stop condition: {0}".format(R.GetOptimizerStopConditionDescription()))
print(" Iteration: {0}".format(R.GetOptimizerIteration()))
print(" Metric value: {0}".format(R.GetMetricValue()))

sitk.WriteTransform(outTx,  sys.argv[3])

if ( not "SITK_NOSHOW" in os.environ ):

    resampler = sitk.ResampleImageFilter()
    resampler.SetReferenceImage(fixed);
    resampler.SetInterpolator(sitk.sitkLinear)
    resampler.SetDefaultPixelValue(100)
    resampler.SetTransform(outTx)

    out = resampler.Execute(moving)
    simg1 = sitk.Cast(sitk.RescaleIntensity(fixed), sitk.sitkUInt8)
    simg2 = sitk.Cast(sitk.RescaleIntensity(out), sitk.sitkUInt8)
    cimg = sitk.Compose(simg1, simg2, simg1//2.+simg2//2.)
    sitk.Show( cimg, "ImageRegistration1 Composition" )
# Run with:
#
# Rscript --vanilla ImageRegistrationMethod1.R fixedImage movingImage outputTransform
#


library(SimpleITK)


commandIteration <- function(method)
{
    msg <- paste("Optimizer iteration number ", method$GetOptimizerIteration(),
                 " = ", method$GetMetricValue(), " : ", method$GetOptimizerPosition(),
                 "\n" )
    cat(msg)
}

args <- commandArgs( TRUE )

if (length(args) != 3) {
    stop("3 arguments expected - FixedImage, MovingImage, TransformFilename")
}

fixed <- ReadImage(args[[1]], 'sitkFloat32')

moving <- ReadImage(args[[2]], 'sitkFloat32')


Reg = ImageRegistrationMethod()
Reg$SetMetricAsMeanSquares()
Reg$SetOptimizerAsRegularStepGradientDescent(4.0, .01, 200 )
Reg$SetInitialTransform(TranslationTransform(fixed$GetDimension()))
Reg$SetInterpolator('sitkLinear')

Reg$AddCommand('sitkIterationEvent', function() commandIteration(Reg))

outTx = Reg$Execute(fixed, moving)

WriteTransform(outTx, args[[3]])