Image Registration Method1¶
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)
{
res <- function() {
msg <- paste("Optimizer iteration number ", method$GetOptimizerIteration(),
" = ", method$GetMetricValue(), " : ", method$GetOptimizerPosition(),
"\n" )
cat(msg)
}
return(res)
}
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', commandIteration(Reg))
outTx = Reg$Execute(fixed, moving)
WriteTransform(outTx, args[[3]])