Image Registration Method Displacement 1
Overview
Code
// 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);
if (m_Method.GetOptimizerIteration() == 0)
{
std::cout << "\tLevel: " << std::setw(3) << m_Method.GetCurrentLevel() << std::endl;
std::cout << "\tScales: " << m_Method.GetOptimizerScales() << std::endl;
}
std::cout << '#' << m_Method.GetOptimizerIteration() << std::endl;
std::cout << "\tMetric Value: " << m_Method.GetMetricValue() << std::endl;
std::cout << "\tLearning Rate: " << m_Method.GetOptimizerLearningRate() << std::endl;
if (m_Method.GetOptimizerConvergenceValue() != std::numeric_limits<double>::max())
{
std::cout << "\tConvergence Value: " << std::scientific << m_Method.GetOptimizerConvergenceValue() << std::endl;
}
std::cout.copyfmt(state);
}
private:
const sitk::ImageRegistrationMethod & m_Method;
};
class MultiResolutionIterationUpdate : public sitk::Command
{
public:
MultiResolutionIterationUpdate(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 << "\tStop Condition: " << m_Method.GetOptimizerStopConditionDescription() << std::endl;
std::cout << "============= Resolution Change =============" << 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::Transform initialTx =
sitk::CenteredTransformInitializer(fixed, moving, sitk::AffineTransform(fixed.GetDimension()));
sitk::ImageRegistrationMethod R;
{
std::vector<unsigned int> shrinkFactors = { 3, 2, 1 };
std::vector<double> smoothingSigmas = { 2.0, 1.0, 1.0 };
R.SetShrinkFactorsPerLevel(shrinkFactors);
R.SetSmoothingSigmasPerLevel(smoothingSigmas);
}
R.SetMetricAsJointHistogramMutualInformation(20);
R.MetricUseFixedImageGradientFilterOff();
R.MetricUseFixedImageGradientFilterOff();
{
double learningRate = 1.0;
unsigned int numberOfIterations = 100;
double convergenceMinimumValue = 1e-6;
unsigned int convergenceWindowSize = 10;
sitk::ImageRegistrationMethod::EstimateLearningRateType estimateLearningRate = R.EachIteration;
R.SetOptimizerAsGradientDescent(
learningRate, numberOfIterations, convergenceMinimumValue, convergenceWindowSize, estimateLearningRate);
}
R.SetOptimizerScalesFromPhysicalShift();
R.SetInitialTransform(initialTx);
R.SetInterpolator(sitk::sitkLinear);
IterationUpdate cmd(R);
R.AddCommand(sitk::sitkIterationEvent, cmd);
MultiResolutionIterationUpdate cmd2(R);
R.AddCommand(sitk::sitkMultiResolutionIterationEvent, cmd2);
sitk::Transform outTx1 = R.Execute(fixed, moving);
std::cout << "-------" << std::endl;
std::cout << outTx1.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::Image displacementField = sitk::Image(fixed.GetSize(), sitk::sitkVectorFloat64);
displacementField.CopyInformation(fixed);
sitk::DisplacementFieldTransform displacementTx(displacementField);
const double varianceForUpdateField = 0.0;
const double varianceForTotalField = 1.5;
displacementTx.SetSmoothingGaussianOnUpdate(varianceForUpdateField, varianceForTotalField);
R.SetMovingInitialTransform(outTx1);
R.SetInitialTransform(displacementTx, true);
R.SetMetricAsANTSNeighborhoodCorrelation(4);
R.MetricUseFixedImageGradientFilterOff();
R.MetricUseFixedImageGradientFilterOff();
{
std::vector<unsigned int> shrinkFactors = { 3, 2, 1 };
std::vector<double> smoothingSigmas = { 2.0, 1.0, 1.0 };
R.SetShrinkFactorsPerLevel(shrinkFactors);
R.SetSmoothingSigmasPerLevel(smoothingSigmas);
}
R.SetOptimizerScalesFromPhysicalShift();
{
double learningRate = 1.0;
unsigned int numberOfIterations = 300;
double convergenceMinimumValue = 1e-6;
unsigned int convergenceWindowSize = 10;
sitk::ImageRegistrationMethod::EstimateLearningRateType estimateLearningRate = R.EachIteration;
R.SetOptimizerAsGradientDescent(
learningRate, numberOfIterations, convergenceMinimumValue, convergenceWindowSize, estimateLearningRate);
}
R.Execute(fixed, moving);
std::cout << "-------" << std::endl;
std::cout << displacementTx.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::CompositeTransform outTx({ outTx1, displacementTx });
sitk::WriteTransform(outTx, argv[3]);
return 0;
}
#!/usr/bin/env python
""" A SimpleITK example demonstrating image registration using the
DisplacementFieldTransform and the JointHistogramMutualInformation
metric. """
import sys
import os
import SimpleITK as sitk
def command_iteration(method):
""" Callback invoked each iteration. """
if method.GetOptimizerIteration() == 0:
print(f"\tLevel: {method.GetCurrentLevel()}")
print(f"\tScales: {method.GetOptimizerScales()}")
print(f"#{method.GetOptimizerIteration()}")
print(f"\tMetric Value: {method.GetMetricValue():10.5f}")
print(f"\tLearningRate: {method.GetOptimizerLearningRate():10.5f}")
if method.GetOptimizerConvergenceValue() != sys.float_info.max:
print("\tConvergence Value: " + f"{method.GetOptimizerConvergenceValue():.5e}")
def command_multiresolution_iteration(method):
""" Callback invoked at the end of each multi-resolution level. """
print(f"\tStop Condition: {method.GetOptimizerStopConditionDescription()}")
print("============= Resolution Change =============")
if len(sys.argv) < 4:
print(
"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)
initialTx = sitk.CenteredTransformInitializer(
fixed, moving, sitk.AffineTransform(fixed.GetDimension())
)
R = sitk.ImageRegistrationMethod()
R.SetShrinkFactorsPerLevel([3, 2, 1])
R.SetSmoothingSigmasPerLevel([2, 1, 1])
R.SetMetricAsJointHistogramMutualInformation(20)
R.MetricUseFixedImageGradientFilterOff()
R.SetOptimizerAsGradientDescent(
learningRate=1.0,
numberOfIterations=100,
estimateLearningRate=R.EachIteration,
)
R.SetOptimizerScalesFromPhysicalShift()
R.SetInitialTransform(initialTx)
R.SetInterpolator(sitk.sitkLinear)
R.AddCommand(sitk.sitkIterationEvent, lambda: command_iteration(R))
R.AddCommand(
sitk.sitkMultiResolutionIterationEvent,
lambda: command_multiresolution_iteration(R),
)
outTx1 = R.Execute(fixed, moving)
print("-------")
print(outTx1)
print(f"Optimizer stop condition: {R.GetOptimizerStopConditionDescription()}")
print(f" Iteration: {R.GetOptimizerIteration()}")
print(f" Metric value: {R.GetMetricValue()}")
displacementField = sitk.Image(fixed.GetSize(), sitk.sitkVectorFloat64)
displacementField.CopyInformation(fixed)
displacementTx = sitk.DisplacementFieldTransform(displacementField)
del displacementField
displacementTx.SetSmoothingGaussianOnUpdate(
varianceForUpdateField=0.0, varianceForTotalField=1.5
)
R.SetMovingInitialTransform(outTx1)
R.SetInitialTransform(displacementTx, inPlace=True)
R.SetMetricAsANTSNeighborhoodCorrelation(4)
R.MetricUseFixedImageGradientFilterOff()
R.SetShrinkFactorsPerLevel([3, 2, 1])
R.SetSmoothingSigmasPerLevel([2, 1, 1])
R.SetOptimizerScalesFromPhysicalShift()
R.SetOptimizerAsGradientDescent(
learningRate=1,
numberOfIterations=300,
estimateLearningRate=R.EachIteration,
)
R.Execute(fixed, moving)
print("-------")
print(displacementTx)
print(f"Optimizer stop condition: {R.GetOptimizerStopConditionDescription()}")
print(f" Iteration: {R.GetOptimizerIteration()}")
print(f" Metric value: {R.GetMetricValue()}")
compositeTx = sitk.CompositeTransform([outTx1, displacementTx])
sitk.WriteTransform(compositeTx, sys.argv[3])
if "SITK_NOSHOW" not in os.environ:
sitk.Show(displacementTx.GetDisplacementField(), "Displacement Field")
resampler = sitk.ResampleImageFilter()
resampler.SetReferenceImage(fixed)
resampler.SetInterpolator(sitk.sitkLinear)
resampler.SetDefaultPixelValue(100)
resampler.SetTransform(compositeTx)
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, "ImageRegistration1 Composition")
# Run with:
#
# Rscript --vanilla ImageRegistrationMethodDisplacement1.R fixedImageFilter movingImageFile outputTransformFile
#
library(SimpleITK)
commandIteration <- function(method)
{
if (method$GetOptimizerIteration()==0) {
cat("\tLevel:", method$GetCurrentLevel(),
"\n\tScales:", method$GetOptimizerScales(), "\n")
}
msg <- paste("#", method$GetOptimizerIteration(),
"\n\tMetric Value: ", method$GetMetricValue(),
"\n\tLearning Rate: ", method$GetOptimizerLearningRate(),
'\n', sep="")
cat(msg)
}
commandMultiIteration <- function(method)
{
msg <- paste("\tStop Condition: ", method$GetOptimizerStopConditionDescription(),
"\n============= Resolution Change =============\n", sep="")
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')
initialTx <- CenteredTransformInitializer(fixed, moving, AffineTransform(fixed$GetDimension()))
R <- ImageRegistrationMethod()
R$SetShrinkFactorsPerLevel(c(3,2,1))
R$SetSmoothingSigmasPerLevel(c(2,1,1))
R$SetMetricAsJointHistogramMutualInformation(20)
R$MetricUseFixedImageGradientFilterOff()
R$SetOptimizerAsGradientDescent(learningRate=1.0,
numberOfIterations=100,
convergenceMinimumValue=1e-6,
convergenceWindowSize=10,
estimateLearningRate = 'EachIteration')
R$SetOptimizerScalesFromPhysicalShift()
R$SetInitialTransform(initialTx)
R$SetInterpolator('sitkLinear')
R$AddCommand( 'sitkIterationEvent', function() commandIteration(R) )
R$AddCommand( 'sitkMultiResolutionIterationEvent', function() commandMultiIteration(R) )
outTx1 <- R$Execute(fixed, moving)
cat("-------\n")
outTx1
cat("Optimizer stop condition:", R$GetOptimizerStopConditionDescription(), '\n')
cat("Iteration:", R$GetOptimizerIteration(), '\n')
cat("Metric value:", R$GetMetricValue(), '\n')
displacementField <- Image(fixed$GetSize(), 'sitkVectorFloat64')
displacementField$CopyInformation(fixed)
displacementTx <- DisplacementFieldTransform(displacementField)
rm(displacementField)
displacementTx$SetSmoothingGaussianOnUpdate(varianceForUpdateField=0.0,
varianceForTotalField=1.5)
R$SetMovingInitialTransform(outTx1)
R$SetInitialTransform(displacementTx, inPlace=TRUE)
R$SetMetricAsANTSNeighborhoodCorrelation(4)
R$MetricUseFixedImageGradientFilterOff()
R$SetShrinkFactorsPerLevel(c(3,2,1))
R$SetSmoothingSigmasPerLevel(c(2,1,1))
R$SetOptimizerScalesFromPhysicalShift()
R$SetOptimizerAsGradientDescent(learningRate=1,
numberOfIterations=300,
convergenceMinimumValue=1e-6,
convergenceWindowSize=10,
estimateLearningRate = 'EachIteration')
R$Execute(fixed, moving)
cat("-------\n")
displacementTx
cat("Optimizer stop condition:", R$GetOptimizerStopConditionDescription(), '\n')
cat("Iteration:", R$GetOptimizerIteration(), '\n')
cat("Metric value:", R$GetMetricValue(), '\n')
outTx <- CompositeTransform(outTx1$GetDimension())
outTx$AddTransform(outTx1)
outTx$AddTransform(displacementTx)
WriteTransform(outTx, args[[3]])