Image Registration Method 2
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', 'displaceMeth2.hdf5')
produces the text and images below.
Output Text
Text Output (click triangle to collapse)
0 = -0.49019 : (0.0028307206468947536, 0.0023491904331766162)
1 = -0.49030 : (0.015968780135700314, 0.01326142597839447)
2 = -0.49068 : (0.08103357252207889, 0.06751733183069583)
3 = -0.49262 : (0.4031918173866258, 0.3356312129926553)
4 = -0.50145 : (1.9891643379821866, 1.6551894820472324)
5 = -0.53055 : (10.002028651827025, 8.700643712668963)
6 = -0.71953 : (14.301161173235505, 16.02941144181758)
7 = -1.14334 : (12.934080126605934, 16.847083218629955)
8 = -1.28582 : (13.021228858202821, 16.980930809778222)
9 = -1.29232 : (13.010721376329036, 17.001152119792206)
10 = -1.29329 : (13.011578940166988, 17.000992823315222)
11 = -1.29327 : (13.01231557664238, 17.000844562218692)
12 = -1.29324 : (13.012426459622057, 17.000820248130232)
13 = -1.29323 : (13.012465052549832, 17.000811591622085)
-------
itk::simple::TranslationTransform
TranslationTransform (0x55eaddff4340)
RTTI typeinfo: itk::TranslationTransform<double, 2u>
Reference Count: 2
Modified Time: 25514
Debug: Off
Object Name:
Observers:
none
Offset: [13.0116, 17.001]
Optimizer stop condition: GradientDescentLineSearchOptimizerv4Template: Convergence checker passed at iteration 14.
Iteration: 14
Metric value: -1.2932878917340445
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 ImageRegistrationMethod2
{
static void Main(string[] args)
{
if (args.Length < 3)
{
Console.WriteLine("Usage: %s <fixedImageFilter> <movingImageFile> <outputTransformFile>\n", "ImageRegistrationMethod2");
return;
}
ImageFileReader reader = new ImageFileReader();
reader.SetOutputPixelType(PixelIDValueEnum.sitkFloat32);
reader.SetFileName(args[0]);
Image fixedImage = reader.Execute();
fixedImage = SimpleITK.Normalize(fixedImage);
SimpleITK.DiscreteGaussian(fixedImage, 2.0);
reader.SetFileName(args[1]);
Image movingImage = reader.Execute();
movingImage=SimpleITK.Normalize(movingImage);
movingImage = SimpleITK.DiscreteGaussian(movingImage, 2.0);
ImageRegistrationMethod R = new ImageRegistrationMethod();
R.SetMetricAsJointHistogramMutualInformation();
double learningRate = 1;
uint numberOfIterations = 200;
double convergenceMinimumValue = 1e-4;
uint convergenceWindowSize = 5;
R.SetOptimizerAsGradientDescentLineSearch(learningRate,
numberOfIterations,
convergenceMinimumValue,
convergenceWindowSize);
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);
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>
#include <numeric>
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(7) << 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 );
fixed = sitk::Normalize( fixed );
fixed = sitk::DiscreteGaussian( fixed, 2.0 );
sitk::Image moving = sitk::ReadImage( argv[2], sitk::sitkFloat32 );
moving = sitk::Normalize( moving );
moving = sitk::DiscreteGaussian( moving, 2.0);
sitk::ImageRegistrationMethod R;
R.SetMetricAsJointHistogramMutualInformation( );
const double learningRate = 1;
const unsigned int numberOfIterations = 200;
const double convergenceMinimumValue = 1e-4;
const unsigned int convergenceWindowSize=5;
R.SetOptimizerAsGradientDescentLineSearch ( learningRate,
numberOfIterations,
convergenceMinimumValue,
convergenceWindowSize);
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;
}
#!/usr/bin/env python
import SimpleITK as sitk
import sys
import os
def command_iteration(method):
print(
f"{method.GetOptimizerIteration():3} "
+ f" = {method.GetMetricValue():7.5f} "
+ f" : {method.GetOptimizerPosition()}"
)
def main(args):
if len(args) < 3:
print(
"Usage:",
"ImageRegistrationMethod2",
"<fixedImageFilter> <movingImageFile> <outputTransformFile>",
)
sys.exit(1)
fixed = sitk.ReadImage(args[1], sitk.sitkFloat32)
fixed = sitk.Normalize(fixed)
fixed = sitk.DiscreteGaussian(fixed, 2.0)
moving = sitk.ReadImage(args[2], sitk.sitkFloat32)
moving = sitk.Normalize(moving)
moving = sitk.DiscreteGaussian(moving, 2.0)
R = sitk.ImageRegistrationMethod()
R.SetMetricAsJointHistogramMutualInformation()
R.SetOptimizerAsGradientDescentLineSearch(
learningRate=1.0,
numberOfIterations=200,
convergenceMinimumValue=1e-5,
convergenceWindowSize=5,
)
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(1)
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 {"fixed": fixed,
"moving": moving,
"composition": cimg}
# Run with:
#
# Rscript --vanilla ImageRegistrationMethod2.R fixedImageFilter movingImageFile outputTransformFile
#
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 arguments expected - fixedImageFilter, movingImageFile, outputTransformFile")
}
pixelType <- 'sitkFloat32'
fixed <- ReadImage(args[[1]], pixelType)
fixed <- Normalize(fixed)
fixed <- DiscreteGaussian(fixed, 2.0)
moving <- ReadImage(args[[2]], pixelType)
moving <- Normalize(moving)
moving <- DiscreteGaussian(moving, 2.0)
R <- ImageRegistrationMethod()
R$SetMetricAsJointHistogramMutualInformation()
R$SetOptimizerAsGradientDescentLineSearch(learningRate=1.0,
numberOfIterations=200,
convergenceMinimumValue=1e-5,
convergenceWindowSize=5)
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]])