DemonsRegistration2
Overview
This example illustrates how to use the fast symmetric forces Demons algorithm. As the name implies, unlike the classical algorithm, the forces are symmetric.
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: DemonsRegistration1.
Code
// 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::FastSymmetricForcesDemonsRegistrationFilter & 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::FastSymmetricForcesDemonsRegistrationFilter & m_Filter;
};
int
main(int argc, char * argv[])
{
if (argc < 4)
{
std::cerr << "Usage: " << argv[0]
<< " <fixedImageFilter> <movingImageFile> [initialTransformFile] <outputTransformFile>" << std::endl;
return 1;
}
sitk::Image fixed = sitk::ReadImage(argv[1]);
sitk::Image moving = sitk::ReadImage(argv[2]);
sitk::HistogramMatchingImageFilter matcher;
if (fixed.GetPixelID() == sitk::sitkUInt8 || fixed.GetPixelID() == sitk::sitkInt8)
{
matcher.SetNumberOfHistogramLevels(128);
}
else
{
matcher.SetNumberOfHistogramLevels(1024);
}
matcher.SetNumberOfMatchPoints(7);
matcher.ThresholdAtMeanIntensityOn();
moving = matcher.Execute(moving, fixed);
sitk::FastSymmetricForcesDemonsRegistrationFilter filter;
IterationUpdate cmd(filter);
filter.AddCommand(sitk::sitkIterationEvent, cmd);
filter.SetNumberOfIterations(200);
filter.SetStandardDeviations(1.0);
sitk::Image displacementField;
if (argc > 4)
{
sitk::Transform initialTransform = sitk::ReadTransform(argv[3]);
argv[3] = argv[4];
--argc;
sitk::TransformToDisplacementFieldFilter toDisplacementFilter;
toDisplacementFilter.SetReferenceImage(fixed);
displacementField = toDisplacementFilter.Execute(initialTransform);
displacementField = filter.Execute(fixed, moving, displacementField);
}
else
{
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;
}
#!/usr/bin/env python
""" A SimpleITK example demonstrating fast symmetric forces Demons image
registation. """
import sys
import os
import SimpleITK as sitk
def command_iteration(sitk_filter):
""" Callback invoked when the filter progresses. """
print(f"{sitk_filter.GetElapsedIterations():3} = {sitk_filter.GetMetric():10.5f}")
if len(sys.argv) < 4:
print(
"Usage:",
sys.argv[0],
"<fixedImageFilter> <movingImageFile>",
"[initialTransformFile] <outputTransformFile>",
)
sys.exit(1)
fixed = sitk.ReadImage(sys.argv[1])
moving = sitk.ReadImage(sys.argv[2])
matcher = sitk.HistogramMatchingImageFilter()
if fixed.GetPixelID() in (sitk.sitkUInt8, sitk.sitkInt8):
matcher.SetNumberOfHistogramLevels(128)
else:
matcher.SetNumberOfHistogramLevels(1024)
matcher.SetNumberOfMatchPoints(7)
matcher.ThresholdAtMeanIntensityOn()
moving = matcher.Execute(moving, fixed)
# The fast symmetric forces Demons Registration Filter
# Note there is a whole family of Demons Registration algorithms included in
# SimpleITK
demons = sitk.FastSymmetricForcesDemonsRegistrationFilter()
demons.SetNumberOfIterations(200)
# Standard deviation for Gaussian smoothing of displacement field
demons.SetStandardDeviations(1.0)
demons.AddCommand(sitk.sitkIterationEvent, lambda: command_iteration(demons))
if len(sys.argv) > 4:
initialTransform = sitk.ReadTransform(sys.argv[3])
sys.argv[-1] = sys.argv.pop()
toDisplacementFilter = sitk.TransformToDisplacementFieldFilter()
toDisplacementFilter.SetReferenceImage(fixed)
displacementField = toDisplacementFilter.Execute(initialTransform)
displacementField = demons.Execute(fixed, moving, displacementField)
else:
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)
cimg = sitk.Compose(simg1, simg2, simg1 // 2.0 + simg2 // 2.0)
sitk.Show(cimg, "DeformableRegistration1 Composition")
# Run with:
#
# Rscript --vanilla DemonsRegistration2.R fixedImageFilter movingImageFile [initialTransformFile] 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 (or 4) arguments expected - fixedImageFilter, movingImageFile, [initialTransformFile], outputTransformFile")
}
fixed <- ReadImage(args[[1]])
moving <- ReadImage(args[[2]])
matcher <- HistogramMatchingImageFilter()
if ( fixed$GetPixelID() %in% c( 'sitkUInt8', 'sitkInt8' ) ) {
matcher$SetNumberOfHistogramLevels(128)
} else {
matcher$SetNumberOfHistogramLevels(1024)
}
matcher$SetNumberOfMatchPoints(7)
matcher$ThresholdAtMeanIntensityOn()
moving <- matcher$Execute(moving, fixed)
# The fast symmetric forces Demons Registration Filter
# Note there is a whole family of Demons Registration algorithms included in SimpleITK
demons <- FastSymmetricForcesDemonsRegistrationFilter()
demons$SetNumberOfIterations(200)
# Standard deviation for Gaussian smoothing of displacement field
demons$SetStandardDeviations(1.0)
demons$AddCommand( 'sitkIterationEvent', function() commandIteration(demons) )
if (length(args) > 3) {
initialTransform <- ReadTransform(args[[3]])
toDisplacementFilter <- TransformToDisplacementFieldFilter()
toDisplacementFilter$SetReferenceImage(fixed)
displacementField <- toDisplacementFilter$Execute(initialTransform)
displacementField <- demons$Execute(fixed, moving, displacementField)
} else {
displacementField <- demons$Execute(fixed, moving)
}
cat("-------\n")
cat("Number Of Iterations: ",demons$GetElapsedIterations(), "\n")
cat(" RMS: ", demons$GetRMSChange(), "\n")
outTx <- DisplacementFieldTransform(displacementField)
WriteTransform(outTx, tail(args, n=1))