DemonsRegistration1
Overview
This example illustrates how to use the classic Demons registration algorithm. The user supplied parameters for the algorithm are the number of iterations and the standard deviations for the Gaussian smoothing of the total displacement field. Additional methods which control regularization, total field smoothing for elastic model or update field smoothing for viscous model are available.
The underlying assumption of the demons framework is that the intensities of homologous points are equal. The example uses histogram matching to make the two images similar prior to registration. This is relevant for registration of MR images where the assumption is not valid. For other imaging modalities where the assumption is valid, such as CT, this step is not necessary. Additionally, the command design pattern is used to monitor registration progress. The resulting deformation field is written to file.
See also: DemonsRegistration2.
Code
// A SimpleITK example demonstrating the classic Demons image registration.
using System;
using itk.simple;
namespace itk.simple.examples
{
public class DemonsRegistration1
{
public class IterationUpdate : Command
{
private DemonsRegistrationFilter filter;
public IterationUpdate(DemonsRegistrationFilter filter)
{
this.filter = filter;
}
public override void Execute()
{
Console.WriteLine(string.Format("{0,3} = {1,10:F5}", filter.GetElapsedIterations(), filter.GetMetric()));
}
}
public static void Main(string[] args)
{
if (args.Length < 3)
{
Console.WriteLine("Usage: DemonsRegistration1 <fixedImageFile> <movingImageFile> <outputTransformFile>");
return;
}
var fixedImage = SimpleITK.ReadImage(args[0], PixelIDValueEnum.sitkFloat32);
var movingImage = SimpleITK.ReadImage(args[1], PixelIDValueEnum.sitkFloat32);
// Perform histogram matching of the moving image to the fixed image
var matcher = new HistogramMatchingImageFilter();
matcher.SetNumberOfHistogramLevels(1024);
matcher.SetNumberOfMatchPoints(7);
matcher.ThresholdAtMeanIntensityOn();
movingImage = matcher.Execute(movingImage, fixedImage);
// The basic Demons Registration Filter
// Note there is a whole family of Demons Registration algorithms included in SimpleITK
var demons = new DemonsRegistrationFilter();
demons.SetNumberOfIterations(50);
// Standard deviation for Gaussian smoothing of displacement field
demons.SetStandardDeviations(1.0);
// Add iteration callback
var cmd = new IterationUpdate(demons);
demons.AddCommand(EventEnum.sitkIterationEvent, cmd);
var displacementField = demons.Execute(fixedImage, movingImage);
Console.WriteLine("-------");
Console.WriteLine(string.Format("Number Of Iterations: {0}", demons.GetElapsedIterations()));
Console.WriteLine(string.Format(" RMS: {0}", demons.GetRMSChange()));
// Create and write the displacement field transform
var outTx = new DisplacementFieldTransform(displacementField);
SimpleITK.WriteTransform(outTx, args[2]);
}
}
}
// This example is based on ITK's DeformableRegistration2.cxx example
#include <SimpleITK.h>
#include <iostream>
#include <stdlib.h>
#include <iomanip>
namespace sitk = itk::simple;
class IterationUpdate : public sitk::Command
{
public:
IterationUpdate(const sitk::DemonsRegistrationFilter & m)
: m_Filter(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_Filter.GetElapsedIterations();
std::cout << " = " << std::setw(10) << m_Filter.GetMetric();
std::cout << std::endl;
std::cout.copyfmt(state);
}
private:
const sitk::DemonsRegistrationFilter & m_Filter;
};
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::HistogramMatchingImageFilter matcher;
matcher.SetNumberOfHistogramLevels(1024);
matcher.SetNumberOfMatchPoints(7);
matcher.ThresholdAtMeanIntensityOn();
moving = matcher.Execute(moving, fixed);
sitk::DemonsRegistrationFilter filter;
IterationUpdate cmd(filter);
filter.AddCommand(sitk::sitkIterationEvent, cmd);
filter.SetNumberOfIterations(50);
filter.SetStandardDeviations(1.0);
sitk::Image displacementField = filter.Execute(fixed, moving);
std::cout << "-------" << std::endl;
std::cout << "Number Of Iterations: " << filter.GetElapsedIterations() << std::endl;
std::cout << " RMS: " << filter.GetRMSChange() << std::endl;
sitk::DisplacementFieldTransform outTx(displacementField);
sitk::WriteTransform(outTx, argv[3]);
return 0;
}
// A SimpleITK example demonstrating the classic Demons image registration.
import org.itk.simple.Image;
import org.itk.simple.SimpleITK;
import org.itk.simple.HistogramMatchingImageFilter;
import org.itk.simple.DemonsRegistrationFilter;
import org.itk.simple.Transform;
import org.itk.simple.PixelIDValueEnum;
import org.itk.simple.EventEnum;
import org.itk.simple.Command;
public class DemonsRegistration1 {
public static void main(String[] args) {
if (args.length < 3) {
System.out.println("Usage: DemonsRegistration1 <fixedImageFile> <movingImageFile> <outputTransformFile>");
System.exit(1);
}
// Read the fixed and moving images as float32
Image fixed = SimpleITK.readImage(args[0], PixelIDValueEnum.sitkFloat32);
Image moving = SimpleITK.readImage(args[1], PixelIDValueEnum.sitkFloat32);
// Perform histogram matching of the moving image to the fixed image
HistogramMatchingImageFilter matcher = new HistogramMatchingImageFilter();
matcher.setNumberOfHistogramLevels(1024);
matcher.setNumberOfMatchPoints(7);
matcher.thresholdAtMeanIntensityOn();
moving = matcher.execute(moving, fixed);
// The basic Demons Registration Filter
// Note there is a whole family of Demons Registration algorithms included in SimpleITK
DemonsRegistrationFilter demons = new DemonsRegistrationFilter();
demons.setNumberOfIterations(50);
// Standard deviation for Gaussian smoothing of displacement field
demons.setStandardDeviations(1.0);
// Add iteration callback
IterationUpdate cmd = new IterationUpdate(demons);
demons.addCommand(EventEnum.sitkIterationEvent, cmd);
Image displacementField = demons.execute(fixed, moving);
System.out.println("-------");
System.out.println("Number Of Iterations: " + demons.getElapsedIterations());
System.out.println(" RMS: " + demons.getRMSChange());
// Create and write the displacement field transform
Transform outTx = new org.itk.simple.DisplacementFieldTransform(displacementField);
SimpleITK.writeTransform(outTx, args[2]);
}
// Callback class for iteration events
static class IterationUpdate extends Command {
private DemonsRegistrationFilter filter;
public IterationUpdate(DemonsRegistrationFilter filter) {
super();
this.filter = filter;
}
public void execute() {
System.out.format("%3d = %10.5f\n", filter.getElapsedIterations(), filter.getMetric());
}
}
}
#!/usr/bin/env python
""" A SimpleITK example demonstrating the classic Demons image registration. """
import sys
import os
import SimpleITK as sitk
def command_iteration(sitk_filter):
""" Callback invoked when the filter processes an iteration. """
print(f"{sitk_filter.GetElapsedIterations():3} = {sitk_filter.GetMetric():10.5f}")
if len(sys.argv) < 4:
print(
f"Usage: {sys.argv[0]}"
+ " <fixedImageFilter> <movingImageFile> <outputTransformFile>"
)
sys.exit(1)
fixed = sitk.ReadImage(sys.argv[1], sitk.sitkFloat32)
moving = sitk.ReadImage(sys.argv[2], sitk.sitkFloat32)
matcher = sitk.HistogramMatchingImageFilter()
matcher.SetNumberOfHistogramLevels(1024)
matcher.SetNumberOfMatchPoints(7)
matcher.ThresholdAtMeanIntensityOn()
moving = matcher.Execute(moving, fixed)
# The basic Demons Registration Filter
# Note there is a whole family of Demons Registration algorithms included in
# SimpleITK
demons = sitk.DemonsRegistrationFilter()
demons.SetNumberOfIterations(50)
# Standard deviation for Gaussian smoothing of displacement field
demons.SetStandardDeviations(1.0)
demons.AddCommand(sitk.sitkIterationEvent, lambda: command_iteration(demons))
displacementField = demons.Execute(fixed, moving)
print("-------")
print(f"Number Of Iterations: {demons.GetElapsedIterations()}")
print(f" RMS: {demons.GetRMSChange()}")
outTx = sitk.DisplacementFieldTransform(displacementField)
sitk.WriteTransform(outTx, sys.argv[3])
if "SITK_NOSHOW" not 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)
# Use the // floor division operator so that the pixel type is
# the same for all three images which is the expectation for
# the compose filter.
cimg = sitk.Compose(simg1, simg2, simg1 // 2.0 + simg2 // 2.0)
sitk.Show(cimg, "DeformableRegistration1 Composition")
# Run with:
#
# Rscript --vanilla DemonsRegistration1.R fixedImageFilter movingImageFile outputTransformFile
#
library(SimpleITK)
commandIteration <- function(filter)
{
msg <- paste("Iteration number ", filter$GetElapsedIterations(),
" = ", filter$GetMetric(), "\n" )
cat(msg)
}
args <- commandArgs( TRUE )
if (length(args) != 3) {
stop("3 arguments expected - fixedImageFilter, movingImageFile, outputTransformFile")
}
fixed <- ReadImage(args[[1]], 'sitkFloat32')
moving <- ReadImage(args[[2]], 'sitkFloat32')
matcher <- HistogramMatchingImageFilter()
matcher$SetNumberOfHistogramLevels(1024)
matcher$SetNumberOfMatchPoints(7)
matcher$ThresholdAtMeanIntensityOn()
moving <- matcher$Execute(moving, fixed)
# The basic Demons Registration Filter
# Note there is a whole family of Demons Registration algorithms included in SimpleITK
demons <- DemonsRegistrationFilter()
demons$SetNumberOfIterations( 50 )
# Standard deviation for Gaussian smoothing of displacement field
demons$SetStandardDeviations( 1.0 )
demons$AddCommand( 'sitkIterationEvent', function() commandIteration(demons) )
displacementField <- demons$Execute( fixed, moving )
cat("-------\n")
cat("Number Of Iterations: ", demons$GetElapsedIterations(), "\n")
cat("RMS: ", demons$GetRMSChange(), "\n")
outTx <- DisplacementFieldTransform( displacementField )
WriteTransform(outTx, args[[3]])