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)
{}
virtual void Execute( )
{
// 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)
{}
virtual void Execute( )
{
// 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;
shrinkFactors.push_back(3);
shrinkFactors.push_back(2);
shrinkFactors.push_back(1);
std::vector<double> smoothingSigmas;
smoothingSigmas.push_back(2.0);
smoothingSigmas.push_back(1.0);
smoothingSigmas.push_back(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, true);
R.SetInterpolator(sitk::sitkLinear);
IterationUpdate cmd(R);
R.AddCommand( sitk::sitkIterationEvent, cmd);
MultiResolutionIterationUpdate cmd2(R);
R.AddCommand( sitk::sitkMultiResolutionIterationEvent, cmd2);
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::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(outTx);
R.SetInitialTransform(displacementTx, true);
R.SetMetricAsANTSNeighborhoodCorrelation(4);
R.MetricUseFixedImageGradientFilterOff();
R.MetricUseFixedImageGradientFilterOff();
{
std::vector<unsigned int> shrinkFactors;
shrinkFactors.push_back(3);
shrinkFactors.push_back(2);
shrinkFactors.push_back(1);
std::vector<double> smoothingSigmas;
smoothingSigmas.push_back(2.0);
smoothingSigmas.push_back(1.0);
smoothingSigmas.push_back(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
);
}
outTx.AddTransform( 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
from __future__ import print_function
import SimpleITK as sitk
import sys
import os
def command_iteration(method) :
if (method.GetOptimizerIteration() == 0):
print("\tLevel: {0}".format(method.GetCurrentLevel()))
print("\tScales: {0}".format(method.GetOptimizerScales()))
print("#{0}".format(method.GetOptimizerIteration()))
print("\tMetric Value: {0:10.5f}".format( method.GetMetricValue()))
print("\tLearningRate: {0:10.5f}".format(method.GetOptimizerLearningRate()))
if (method.GetOptimizerConvergenceValue() != sys.float_info.max):
print("\tConvergence Value: {0:.5e}".format(method.GetOptimizerConvergenceValue()))
def command_multiresolution_iteration(method):
print("\tStop Condition: {0}".format(method.GetOptimizerStopConditionDescription()))
print("============= Resolution Change =============")
if len ( sys.argv ) < 4:
print( "Usage: {0} <fixedImageFilter> <movingImageFile> <outputTransformFile>".format(sys.argv[0]))
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,inPlace=True)
R.SetInterpolator(sitk.sitkLinear)
R.AddCommand( sitk.sitkIterationEvent, lambda: command_iteration(R) )
R.AddCommand( sitk.sitkMultiResolutionIterationEvent, lambda: command_multiresolution_iteration(R) )
outTx = R.Execute(fixed, moving)
print("-------")
print(outTx)
print("Optimizer stop condition: {0}".format(R.GetOptimizerStopConditionDescription()))
print(" Iteration: {0}".format(R.GetOptimizerIteration()))
print(" Metric value: {0}".format(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(outTx)
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)
outTx.AddTransform( R.Execute(fixed, moving) )
print("-------")
print(outTx)
print("Optimizer stop condition: {0}".format(R.GetOptimizerStopConditionDescription()))
print(" Iteration: {0}".format(R.GetOptimizerIteration()))
print(" Metric value: {0}".format(R.GetMetricValue()))
sitk.WriteTransform(outTx, sys.argv[3])
if ( not "SITK_NOSHOW" in os.environ ):
sitk.Show(displacementTx.GetDisplacementField(), "Displacement Field")
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.+simg2//2.)
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, inPlace=TRUE)
R$SetInterpolator('sitkLinear')
R$AddCommand( 'sitkIterationEvent', function() commandIteration(R) )
R$AddCommand( 'sitkMultiResolutionIterationEvent', function() commandMultiIteration(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')
displacementField <- Image(fixed$GetSize(), 'sitkVectorFloat64')
displacementField$CopyInformation(fixed)
displacementTx <- DisplacementFieldTransform(displacementField)
rm(displacementField)
displacementTx$SetSmoothingGaussianOnUpdate(varianceForUpdateField=0.0,
varianceForTotalField=1.5)
R$SetMovingInitialTransform(outTx)
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')
outTx$AddTransform( 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]])