Image Registration Method 4
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', 'displaceMeth4.hdf5')
produces the text and images below.
Output Text
Text Output (click triangle to collapse)
0 = -0.33742 : (0.8883551991522953, 0.4591568796599763)
1 = -0.35088 : (1.7268176481128081, 1.004116255836898)
2 = -0.36229 : (2.4187469173877156, 1.7260815502387081)
3 = -0.37514 : (3.3203576953706744, 2.158629819250346)
4 = -0.38488 : (3.9693971620112545, 2.9193845568254146)
5 = -0.39582 : (4.580401367036123, 3.711011906935255)
6 = -0.40560 : (4.982807469043353, 4.626473172673651)
7 = -0.41584 : (5.724526806050685, 5.297183561074937)
8 = -0.42722 : (6.3011461810254055, 6.114196471868739)
9 = -0.43198 : (6.879495374940256, 6.929985788970921)
10 = -0.44823 : (7.679696390998657, 7.529717673982133)
11 = -0.47282 : (8.372140895644197, 8.251188815460566)
12 = -0.50020 : (9.027329551231393, 9.006654119488385)
13 = -0.51452 : (9.406946086345647, 9.931798046814208)
14 = -0.54356 : (9.854355564205843, 10.826127269819287)
15 = -0.57008 : (10.264665927853011, 11.738073214206194)
16 = -0.60484 : (10.776763753895718, 12.597000345319383)
17 = -0.65100 : (11.404217723671293, 13.37565400119783)
18 = -0.70046 : (11.871902871940648, 14.25954913190524)
19 = -0.77212 : (12.347454971019081, 15.139236690679432)
20 = -0.86938 : (12.764835887967642, 16.0479683178209)
21 = -1.04948 : (13.05232647713824, 17.005751780373025)
22 = -1.41999 : (12.554720245893066, 16.95688421826918)
23 = -1.23275 : (12.803933391805609, 16.97670360030731)
24 = -1.35483 : (13.053245682820764, 16.99523415747041)
25 = -1.41953 : (12.928734542524285, 17.00627842938261)
26 = -1.41252 : (12.990848387738, 16.99934155452668)
27 = -1.42587 : (13.05330474421836, 16.9970062739872)
28 = -1.41959 : (13.022106487167733, 16.998803843736308)
29 = -1.42511 : (12.99088273126701, 17.000084299691075)
30 = -1.42589 : (13.006454996306088, 16.99880165143347)
31 = -1.42614 : (12.998831820095416, 17.000511136944098)
32 = -1.42621 : (13.002653505928967, 16.999702739484835)
33 = -1.42621 : (13.00646513960582, 16.99884819997377)
34 = -1.42614 : (13.004549764853339, 16.99923034761956)
-------
itk::simple::TranslationTransform
TranslationTransform (0x57bbc6fdbe80)
RTTI typeinfo: itk::TranslationTransform<double, 2u>
Reference Count: 2
Modified Time: 13969
Debug: Off
Object Name:
Observers:
none
Offset: [13.0045, 16.9992]
Optimizer stop condition: RegularStepGradientDescentOptimizerv4: Step too small after 35 iterations. Current step (0.000976562) is less than minimum step (0.001).
Iteration: 36
Metric value: -1.426184883438515
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.Write("{0,3} = {1,10:F5} : [",
m_Method.GetOptimizerIteration(),
m_Method.GetMetricValue());
for (int i = 0; i < pos.Count; i++)
{
if (i > 0) Console.Write(", ");
Console.Write("{0:F5}", pos[i]);
}
Console.WriteLine("]");
}
}
class ImageRegistrationMethod4
{
static void Main(string[] args)
{
if (args.Length < 3)
{
Console.WriteLine("Usage: ImageRegistrationMethod4 <fixedImageFile> <movingImageFile> <outputTransformFile> [numberOfBins] [samplingPercentage]");
return;
}
Image fixedImage = SimpleITK.ReadImage(args[0], PixelIDValueEnum.sitkFloat32);
Image movingImage = SimpleITK.ReadImage(args[1], PixelIDValueEnum.sitkFloat32);
uint numberOfBins = 24;
double samplingPercentage = 0.10;
if (args.Length > 3)
{
numberOfBins = uint.Parse(args[3]);
}
if (args.Length > 4)
{
samplingPercentage = double.Parse(args[4]);
}
ImageRegistrationMethod R = new ImageRegistrationMethod();
R.SetMetricAsMattesMutualInformation(numberOfBins);
R.SetMetricSamplingPercentage(samplingPercentage);
R.SetMetricSamplingStrategy(ImageRegistrationMethod.MetricSamplingStrategyType.RANDOM);
R.SetOptimizerAsRegularStepGradientDescent(1.0, 0.001, 200);
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);
Console.WriteLine("-------");
Console.WriteLine(outTx.ToString());
Console.WriteLine("Optimizer stop condition: " + R.GetOptimizerStopConditionDescription());
Console.WriteLine(" Iteration: " + R.GetOptimizerIteration());
Console.WriteLine(" Metric value: " + R.GetMetricValue());
SimpleITK.WriteTransform(outTx, args[2]);
ResampleImageFilter resampler = new ResampleImageFilter();
resampler.SetReferenceImage(fixedImage);
resampler.SetInterpolator(InterpolatorEnum.sitkLinear);
resampler.SetDefaultPixelValue(100);
resampler.SetTransform(outTx);
Image output = resampler.Execute(movingImage);
Image simg1 = SimpleITK.Cast(SimpleITK.RescaleIntensity(fixedImage), PixelIDValueEnum.sitkUInt8);
Image simg2 = SimpleITK.Cast(SimpleITK.RescaleIntensity(output), PixelIDValueEnum.sitkUInt8);
Image cimg = SimpleITK.Compose(simg1, simg2, SimpleITK.Divide(SimpleITK.Add(simg1, simg2), 2));
// Show the composition image if SITK_NOSHOW environment variable is not set
if (Environment.GetEnvironmentVariable("SITK_NOSHOW") == null)
{
SimpleITK.Show(cimg, "ImageRegistration4 Composition");
}
}
}
}
#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::cout << "Usage: " << argv[0]
<< " <fixedImageFile> <movingImageFile> <outputTransformFile> [numberOfBins] [samplingPercentage]"
<< std::endl;
return 1;
}
sitk::Image fixed = sitk::ReadImage(argv[1], sitk::sitkFloat32);
sitk::Image moving = sitk::ReadImage(argv[2], sitk::sitkFloat32);
unsigned int numberOfBins = 24;
double samplingPercentage = 0.10;
if (argc > 4)
{
numberOfBins = atoi(argv[4]);
}
if (argc > 5)
{
samplingPercentage = atof(argv[5]);
}
sitk::ImageRegistrationMethod R;
R.SetMetricAsMattesMutualInformation(numberOfBins);
R.SetMetricSamplingPercentage(samplingPercentage);
R.SetMetricSamplingStrategy(sitk::ImageRegistrationMethod::RANDOM);
R.SetOptimizerAsRegularStepGradientDescent(1.0, 0.001, 200);
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]);
sitk::ResampleImageFilter resampler;
resampler.SetReferenceImage(fixed);
resampler.SetInterpolator(sitk::sitkLinear);
resampler.SetDefaultPixelValue(100);
resampler.SetTransform(outTx);
sitk::Image out = resampler.Execute(moving);
sitk::Image simg1 = sitk::Cast(sitk::RescaleIntensity(fixed), sitk::sitkUInt8);
sitk::Image simg2 = sitk::Cast(sitk::RescaleIntensity(out), sitk::sitkUInt8);
sitk::Image cimg = sitk::Compose(simg1, simg2, sitk::Divide(sitk::Add(simg1, simg2), 2));
// Show the composition image if SITK_NOSHOW environment variable is not set
if (getenv("SITK_NOSHOW") == nullptr)
{
sitk::Show(cimg, "ImageRegistration4 Composition");
}
return 0;
}
import org.itk.simple.*;
import java.text.DecimalFormat;
class IterationUpdate extends Command {
private ImageRegistrationMethod method;
public IterationUpdate(ImageRegistrationMethod m) {
method = m;
}
public void execute() {
VectorDouble pos = method.getOptimizerPosition();
DecimalFormat df = new DecimalFormat("0.00000");
System.out.print(String.format("%3d = %10s : [",
method.getOptimizerIteration(),
df.format(method.getMetricValue())));
for (int i = 0; i < pos.size(); i++) {
if (i > 0) System.out.print(", ");
System.out.print(df.format(pos.get(i)));
}
System.out.println("]");
}
}
public class ImageRegistrationMethod4 {
public static void main(String[] args) throws Exception {
if (args.length < 3) {
System.out.println("Usage: ImageRegistrationMethod4 <fixedImageFile> <movingImageFile> <outputTransformFile> [numberOfBins] [samplingPercentage]");
System.exit(1);
}
Image fixed = SimpleITK.readImage(args[0], PixelIDValueEnum.sitkFloat32);
Image moving = SimpleITK.readImage(args[1], PixelIDValueEnum.sitkFloat32);
int numberOfBins = 24;
double samplingPercentage = 0.10;
if (args.length > 3) {
numberOfBins = Integer.parseInt(args[3]);
}
if (args.length > 4) {
samplingPercentage = Double.parseDouble(args[4]);
}
ImageRegistrationMethod R = new ImageRegistrationMethod();
R.setMetricAsMattesMutualInformation(numberOfBins);
R.setMetricSamplingPercentage(samplingPercentage);
R.setMetricSamplingStrategy(ImageRegistrationMethod.MetricSamplingStrategyType.RANDOM);
R.setOptimizerAsRegularStepGradientDescent(1.0, 0.001, 200);
R.setInitialTransform(new TranslationTransform(fixed.getDimension()));
R.setInterpolator(InterpolatorEnum.sitkLinear);
IterationUpdate cmd = new IterationUpdate(R);
R.addCommand(EventEnum.sitkIterationEvent, cmd);
Transform outTx = R.execute(fixed, moving);
System.out.println("-------");
System.out.println(outTx.toString());
System.out.println("Optimizer stop condition: " + R.getOptimizerStopConditionDescription());
System.out.println(" Iteration: " + R.getOptimizerIteration());
System.out.println(" Metric value: " + R.getMetricValue());
SimpleITK.writeTransform(outTx, args[2]);
ResampleImageFilter resampler = new ResampleImageFilter();
resampler.setReferenceImage(fixed);
resampler.setInterpolator(InterpolatorEnum.sitkLinear);
resampler.setDefaultPixelValue(100);
resampler.setTransform(outTx);
Image output = resampler.execute(moving);
Image simg1 = SimpleITK.cast(SimpleITK.rescaleIntensity(fixed), PixelIDValueEnum.sitkUInt8);
Image simg2 = SimpleITK.cast(SimpleITK.rescaleIntensity(output), PixelIDValueEnum.sitkUInt8);
Image cimg = SimpleITK.compose(simg1, simg2, SimpleITK.divide(SimpleITK.add(simg1, simg2), 2.0));
// Show the composition image if SITK_NOSHOW environment variable is not set
if (System.getenv("SITK_NOSHOW") == null) {
SimpleITK.show(cimg, "ImageRegistration4 Composition");
}
}
}
#!/usr/bin/env python
import sys
import os
import SimpleITK as sitk
def command_iteration(method):
""" Callback invoked when the optimization has an iteration. """
print(
f"{method.GetOptimizerIteration():3} "
+ f"= {method.GetMetricValue():10.5f} "
+ f": {method.GetOptimizerPosition()}"
)
def main(args):
""" A SimpleITK example demonstrating image registration using Mattes mutual
information as the metric. """
if len(args) < 3:
print(
"Usage:",
"ImageRegistrationMethod4",
"<fixedImageFilter> <movingImageFile>",
"<outputTransformFile> <numberOfBins> <samplingPercentage>",
)
sys.exit(1)
fixed = sitk.ReadImage(args[1], sitk.sitkFloat32)
moving = sitk.ReadImage(args[2], sitk.sitkFloat32)
numberOfBins = 24
samplingPercentage = 0.10
if len(args) > 4:
numberOfBins = int(args[4])
if len(args) > 5:
samplingPercentage = float(args[5])
R = sitk.ImageRegistrationMethod()
R.SetMetricAsMattesMutualInformation(numberOfBins)
R.SetMetricSamplingPercentage(samplingPercentage, sitk.sitkWallClock)
R.SetMetricSamplingStrategy(R.RANDOM)
R.SetOptimizerAsRegularStepGradientDescent(1.0, 0.001, 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)
library(SimpleITK)
commandIteration <- function(method)
{
msg <- paste(method$GetOptimizerIteration(), "=",
method$GetMetricValue(), ":",
method$GetOptimizerPosition(), "\n" )
cat(msg)
}
args <- commandArgs( TRUE )
if (length(args) < 3) {
stop("3, 4, or 5 arguments expected - fixedImageFilter, movingImageFile, outputTransformFile [numberOfBins] [samplingPercentage]")
}
fixed <- ReadImage(args[[1]], 'sitkFloat32')
moving <- ReadImage(args[[2]], 'sitkFloat32')
numberOfBins <- 24
samplingPercentage <- 0.10
if (length(args) > 4) {
numberOfBins <- strtoi(args[[4]])
}
if (length(args) > 5) {
samplingPercentage <- as.numeric(args[[5]])
}
R <- ImageRegistrationMethod()
R$SetMetricAsMattesMutualInformation(numberOfBins)
R$SetMetricSamplingPercentage(samplingPercentage)
R$SetMetricSamplingStrategy('RANDOM')
R$SetOptimizerAsRegularStepGradientDescent(1.0,.001,200)
R$SetInitialTransform(TranslationTransform(fixed$GetDimension()))
R$SetInterpolator('sitkLinear')
R$AddCommand( 'sitkIterationEvent', function() commandIteration(R) )
outTx <- R$Execute(fixed, moving)
cat("-------\n")
outTx
cat("Optimizer stop condition:", R$GetOptimizerStopConditionDescription(), '\n')
cat("Iteration:", R$GetOptimizerIteration(), '\n')
cat("Metric value:", R$GetMetricValue(), '\n')
WriteTransform(outTx, args[[3]])