Image Registration Method 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', 'BrainProtonDensitySliceShifted13x17y.png', 'displaceMeth1.hdf5')
produces the text and images below.
Output Text
Text Output (click triangle to collapse)
0 = 4499.45303 : (2.928695951245586, 2.7244705953923805)
1 = 3860.83626 : (6.135143776902402, 5.115849348610004)
2 = 3508.02147 : (8.822660051952475, 8.078492808653918)
3 = 3117.30771 : (10.968558473732326, 11.454158663474674)
4 = 2125.42619 : (13.105290365964755, 14.835634202454191)
5 = 911.30826 : (12.75173580401588, 18.819978461140323)
6 = 741.41676 : (13.139053510563274, 16.857840597942413)
7 = 16.89176 : (12.356787624301035, 17.480785285045815)
8 = 233.71446 : (12.79212443526829, 17.234854683011704)
9 = 39.80270 : (13.167510875734614, 16.904574468172815)
10 = 16.57314 : (12.938831371165355, 17.005597654570586)
11 = 1.68763 : (13.063495692092735, 16.996443033457986)
12 = 1.79437 : (13.001061362657559, 16.999307384689935)
13 = 0.00076 : (12.945418587211314, 17.0277701944711)
14 = 1.74802 : (12.974454390534774, 17.01621663980765)
15 = 0.43025 : (13.002439510423766, 17.002309966416835)
16 = 0.00532 : (12.989877586882951, 16.99301810428082)
-------
itk::simple::TranslationTransform
TranslationTransform (0x5639ff927c90)
RTTI typeinfo: itk::TranslationTransform<double, 2u>
Reference Count: 2
Modified Time: 2526
Debug: Off
Object Name:
Observers:
none
Offset: [12.9899, 16.993]
Optimizer stop condition: RegularStepGradientDescentOptimizerv4: Step too small after 17 iterations. Current step (0.0078125) is less than minimum step (0.01).
Iteration: 18
Metric value: 0.07213459794461137
Input Images
Fixed Image |
Moving Image |
Output Image
Composition Image
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)
{}
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_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
import sys
import os
import SimpleITK as sitk
def command_iteration(method):
""" Callback invoked when the optimization process is performing an iteration. """
print(
f"{method.GetOptimizerIteration():3} "
+ f"= {method.GetMetricValue():10.5f} "
+ f": {method.GetOptimizerPosition()}"
)
def main(args):
""" A basic SimpleITK image registration example. """
if len(args) < 3:
print(
"Usage:",
"ImageRegistrationMethod1",
"<fixedImageFilter> <movingImageFile>",
"<outputTransformFile>",
)
sys.exit(1)
fixed = sitk.ReadImage(args[1], sitk.sitkFloat32)
moving = sitk.ReadImage(args[2], sitk.sitkFloat32)
R = sitk.ImageRegistrationMethod()
R.SetMetricAsMeanSquares()
R.SetOptimizerAsRegularStepGradientDescent(4.0, 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(f"Optimizer stop condition: {R.GetOptimizerStopConditionDescription()}")
print(f" Iteration: {R.GetOptimizerIteration()}")
print(f" Metric value: {R.GetMetricValue()}")
sitk.WriteTransform(outTx, args[3])
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.0 + simg2 // 2.0)
return_images = {"fixed": fixed, "moving": moving, "composition": cimg}
# 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]])