Image Registration Method BSpline 1
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', 'BrainProtonDensitySliceBSplined10.png', 'ImageRegistrationMethodBSpline1.hdf5')
produces the text and images below.
Output Text
Text Output (click triangle to collapse)
Initial Parameters:
(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
0 = -0.83872
0 = -0.85295
0 = -0.89317
0 = -0.91369
0 = -0.91369
1 = -0.95662
1 = -0.95662
2 = -0.96795
2 = -0.96795
3 = -0.97045
3 = -0.97045
4 = -0.97769
4 = -0.97769
5 = -0.98248
5 = -0.98248
6 = -0.98273
6 = -0.98273
7 = -0.98739
7 = -0.98739
8 = -0.98956
8 = -0.98956
9 = -0.99150
9 = -0.99150
10 = -0.99169
10 = -0.99169
11 = -0.99289
11 = -0.99289
12 = -0.99322
12 = -0.99322
13 = -0.99362
13 = -0.99362
14 = -0.99409
14 = -0.99409
15 = -0.99400
15 = -0.99420
15 = -0.99420
16 = -0.99434
16 = -0.99434
17 = -0.99444
17 = -0.99444
18 = -0.99456
18 = -0.99456
19 = -0.99462
19 = -0.99462
20 = -0.99467
20 = -0.99467
21 = -0.99469
21 = -0.99469
22 = -0.99471
22 = -0.99471
23 = -0.99473
23 = -0.99473
24 = -0.99476
24 = -0.99476
25 = -0.99478
25 = -0.99478
26 = -0.99481
26 = -0.99481
27 = -0.99484
27 = -0.99484
28 = -0.99486
28 = -0.99486
29 = -0.99487
29 = -0.99487
30 = -0.99488
30 = -0.99488
-------
itk::simple::BSplineTransform
BSplineTransform (0x55db2b83e690)
RTTI typeinfo: itk::BSplineTransform<double, 2u, 3u>
Reference Count: 3
Modified Time: 2650
Debug: Off
Object Name:
Observers:
none
CoefficientImage: [ 0x55db2b7ad440, 0x55db2b7a0130 ]
TransformDomainOrigin: [-0.625, -0.625]
TransformDomainPhysicalDimensions: [221.25, 257.25]
TransformDomainDirection: 1 0
0 1
TransformDomainMeshSize: [8, 8]
GridSize: [11, 11]
GridOrigin: [-28.2812, -32.7812]
GridSpacing: [27.6562, 32.1562]
GridDirection: 1 0
0 1
Optimizer stop condition: LBFGSBOptimizerv4: Gradient tolerance reached. Gradient tolerance is 1e-05
Iteration: 31
Metric value: -0.9948792040893044
Input Images
Fixed Image |
Moving Image |
Output Image
Composition Image
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;
// use sitk's output operator for std::vector etc..
using sitk::operator<<;
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<<;
if (m_Method.GetOptimizerIteration() == 0)
{
std::cout << m_Method.ToString() << std::endl;
}
// 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(10) << m_Method.GetMetricValue() << 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 );
std::vector<unsigned int> transformDomainMeshSize(fixed.GetDimension(),8);
sitk::BSplineTransform tx = sitk::BSplineTransformInitializer(fixed, transformDomainMeshSize);
std::cout << "Initial Parameters:" << tx.GetParameters() << std::endl;
sitk::ImageRegistrationMethod R;
R.SetMetricAsCorrelation();
const double gradientConvergenceTolerance = 1e-5;
const unsigned int maximumNumberOfIterations = 100;
const unsigned int maximumNumberOfCorrections = 5;
const unsigned int maximumNumberOfFunctionEvaluations = 1000;
const double costFunctionConvergenceFactor = 1e+7;
R.SetOptimizerAsLBFGSB(gradientConvergenceTolerance,
maximumNumberOfIterations,
maximumNumberOfCorrections,
maximumNumberOfFunctionEvaluations,
costFunctionConvergenceFactor);
R.SetInitialTransform(tx, true);
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():10.5f}")
def main(args):
if len(args) < 4:
print(
"Usage:",
sys.args[0],
"<fixedImageFilter> <movingImageFile>",
"<outputTransformFile>",
)
sys.exit(1)
fixed = sitk.ReadImage(args[1], sitk.sitkFloat32)
moving = sitk.ReadImage(args[2], sitk.sitkFloat32)
transformDomainMeshSize = [8] * moving.GetDimension()
tx = sitk.BSplineTransformInitializer(fixed, transformDomainMeshSize)
print("Initial Parameters:")
print(tx.GetParameters())
R = sitk.ImageRegistrationMethod()
R.SetMetricAsCorrelation()
R.SetOptimizerAsLBFGSB(
gradientConvergenceTolerance=1e-5,
numberOfIterations=100,
maximumNumberOfCorrections=5,
maximumNumberOfFunctionEvaluations=1000,
costFunctionConvergenceFactor=1e7,
)
R.SetInitialTransform(tx, True)
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(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)
return_images = {"fixed": fixed, "moving": moving, "composition": cimg}
return return_images
if __name__ == "__main__":
return_dict = main(sys.argv)
if "SITK_NOSHOW" not in os.environ:
sitk.Show(
return_dict["composition"], "ImageRegistrationMethodBSpline1 Composition"
)
# Run with:
#
# Rscript --vanilla ImageRegistrationMethodBSpline1.R fixedImageFilter movingImageFile outputTransformFile
#
library(SimpleITK)
commandIteration <- function(method)
{
msg <- paste(method$GetOptimizerIteration(), "=",
method$GetMetricValue(), "\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')
transformDomainMeshSize <- rep(8, moving$GetDimension())
tx <- BSplineTransformInitializer(fixed, transformDomainMeshSize)
cat("Initial Parameters:\n", tx$GetParameters())
R <- ImageRegistrationMethod()
R$SetMetricAsCorrelation()
R$SetOptimizerAsLBFGSB(gradientConvergenceTolerance=1e-5,
numberOfIterations=100,
maximumNumberOfCorrections=5,
maximumNumberOfFunctionEvaluations=1000,
costFunctionConvergenceFactor=1e+7)
R$SetInitialTransform(tx, TRUE)
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]])