Image Registration Method 4
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', 'BrainProtonDensitySliceShifted13x17y.png', 'displaceMeth4.hdf5')
produces the text and images below.
Output Text
Text Output (click triangle to collapse)
0 = -0.35281 : (0.6506890423009412, 0.7593443028228264)
1 = -0.36947 : (1.0847010252556577, 1.6602514005036664)
2 = -0.37753 : (1.5413893290165088, 2.549878172362538)
3 = -0.38791 : (2.3486985622960708, 3.1400068039858184)
4 = -0.39459 : (3.1407889572218046, 3.7504106098933465)
5 = -0.40522 : (3.977521702575541, 4.298022071472447)
6 = -0.41588 : (4.927096046366175, 4.611564039340527)
7 = -0.42514 : (5.639639274425346, 5.3131922515348595)
8 = -0.44618 : (6.434197584547169, 5.92038010699341)
9 = -0.46476 : (7.225182652335237, 6.5322155608012995)
10 = -0.47386 : (7.77689813389309, 7.366247950705544)
11 = -0.49074 : (8.50534836297095, 8.051346676260053)
12 = -0.50797 : (9.397941585500613, 8.502209559068245)
13 = -0.52170 : (9.807826555920268, 9.414346780666003)
14 = -0.53753 : (10.284459154041546, 10.293449372180874)
15 = -0.56091 : (10.683612086837567, 11.210333735795139)
16 = -0.59704 : (11.016609816815837, 12.153261363892502)
17 = -0.63377 : (11.555594460854683, 12.995576954090332)
18 = -0.68419 : (12.204872581720132, 13.756128018625406)
19 = -0.75776 : (12.698078997720419, 14.626040330997263)
20 = -0.84399 : (12.978351858374943, 15.585960704559847)
21 = -0.97169 : (12.920581812898137, 16.584290620892922)
22 = -1.24824 : (13.124196730308128, 17.563341674417774)
23 = -1.18551 : (12.989055484405798, 17.08195114308592)
24 = -1.38960 : (13.117124808555902, 16.598631116383388)
25 = -1.25009 : (13.018577883743488, 16.828388605075304)
26 = -1.35771 : (12.978365259237886, 17.0751333007232)
27 = -1.39052 : (13.03244981169743, 16.962439639919893)
28 = -1.39810 : (12.986266884835361, 17.004551244810383)
29 = -1.40246 : (13.016734540629988, 16.99760254346552)
30 = -1.40270 : (13.001215466998131, 16.9994188525276)
31 = -1.40331 : (13.009027909663939, 16.99944878317189)
32 = -1.40318 : (13.005122832969363, 16.9995445175543)
33 = -1.40329 : (13.001216583566658, 16.99954667773063)
34 = -1.40331 : (13.003168233452373, 16.99947078313784)
35 = -1.40332 : (13.005112439687908, 16.999657221812296)
-------
itk::simple::TranslationTransform
TranslationTransform (0x55aaa28d2310)
RTTI typeinfo: itk::TranslationTransform<double, 2u>
Reference Count: 2
Modified Time: 13972
Debug: Off
Object Name:
Observers:
none
Offset: [13.0051, 16.9997]
Optimizer stop condition: RegularStepGradientDescentOptimizerv4: Step too small after 36 iterations. Current step (0.000976562) is less than minimum step (0.001).
Iteration: 37
Metric value: -1.4032920395781017
Input Images
Fixed Image |
Moving Image |
Output Image
Code
#!/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} "
+ f": {method.GetOptimizerPosition()}"
)
def main(args):
if len(args) < 3:
print(
"Usage:",
"ImageRegistrationMethod4",
"<fixedImageFilter> <movingImageFile>",
"<outputTransformFile> <numberOfBins> <samplingPercentage>",
)
sys.exit(1)
fixed = sitk.ReadImage(args[1], sitk.sitkFloat32)
moving = sitk.ReadImage(args[2], sitk.sitkFloat32)
numberOfBins = 24
samplingPercentage = 0.10
if len(args) > 4:
numberOfBins = int(args[4])
if len(args) > 5:
samplingPercentage = float(args[5])
R = sitk.ImageRegistrationMethod()
R.SetMetricAsMattesMutualInformation(numberOfBins)
R.SetMetricSamplingPercentage(samplingPercentage, sitk.sitkWallClock)
R.SetMetricSamplingStrategy(R.RANDOM)
R.SetOptimizerAsRegularStepGradientDescent(1.0, 0.001, 200)
R.SetInitialTransform(sitk.TranslationTransform(fixed.GetDimension()))
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 {"fixed": fixed, "moving": moving, "composition": cimg}
if __name__ == "__main__":
# Run with:
#
# Rscript --vanilla ImageRegistrationMethod4.R fixedImageFilter movingImageFile outputTransformFile numberOfBins samplingPercentage
#
library(SimpleITK)
commandIteration <- function(method)
{
msg <- paste(method$GetOptimizerIteration(), "=",
method$GetMetricValue(), ":",
method$GetOptimizerPosition(), "\n" )
cat(msg)
}
args <- commandArgs( TRUE )
if (length(args) < 3) {
stop("3, 4, or 5 arguments expected - fixedImageFilter, movingImageFile, outputTransformFile [numberOfBins] [samplingPercentage]")
}
fixed <- ReadImage(args[[1]], 'sitkFloat32')
moving <- ReadImage(args[[2]], 'sitkFloat32')
numberOfBins <- 24
samplingPercentage <- 0.10
if (length(args) > 4) {
numberOfBins <- strtoi(args[[4]])
}
if (length(args) > 5) {
samplingPercentage <- as.numeric(args[[5]])
}
R <- ImageRegistrationMethod()
R$SetMetricAsMattesMutualInformation(numberOfBins)
R$SetMetricSamplingPercentage(samplingPercentage)
R$SetMetricSamplingStrategy('RANDOM')
R$SetOptimizerAsRegularStepGradientDescent(1.0,.001,200)
R$SetInitialTransform(TranslationTransform(fixed$GetDimension()))
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]])