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.34026 : (0.7150826245543834, 0.6990399416774515)
  1 =   -0.35205 : (1.5577459253547459, 1.23748080331802)
  2 =   -0.36277 : (2.282236348047206, 1.9267656699928226)
  3 =   -0.37449 : (3.118733209581384, 2.4747373851784427)
  4 =   -0.38634 : (4.09270414705766, 2.701410305818701)
  5 =   -0.39799 : (4.9610729066756765, 3.197329348929034)
  6 =   -0.41151 : (5.880882620608181, 3.5896940901525602)
  7 =   -0.42005 : (6.401552950306754, 4.443451905795603)
  8 =   -0.42927 : (7.094378952423701, 5.164556704552793)
  9 =   -0.44144 : (7.833930920236335, 5.837656167414308)
 10 =   -0.45652 : (8.462925318384015, 6.6150659995484515)
 11 =   -0.47360 : (8.922381978434776, 7.503266189559127)
 12 =   -0.48661 : (9.554036544928254, 8.278516154814165)
 13 =   -0.50874 : (10.275136644966015, 8.971347047403064)
 14 =   -0.52771 : (10.709832154605628, 9.87192453651722)
 15 =   -0.55834 : (11.18417103054953, 10.752266873774877)
 16 =   -0.58897 : (11.705617339015628, 11.605550969130626)
 17 =   -0.63189 : (12.333926798492271, 12.38351447928529)
 18 =   -0.67967 : (12.558802569237915, 13.357901921595356)
 19 =   -0.73032 : (12.811442650435938, 14.325462250139399)
 20 =   -0.81470 : (13.1553746522357, 15.264456807172616)
 21 =   -0.92386 : (12.984869538324672, 16.249813598456964)
 22 =   -1.12096 : (13.054342391067136, 17.24739744091424)
 23 =   -1.34008 : (12.921151561968108, 16.765463624873167)
 24 =   -1.32956 : (13.029991886723995, 16.990527579577233)
 25 =   -1.41778 : (12.910272374682103, 17.026475297675205)
 26 =   -1.39755 : (12.971922122766422, 17.016201128066914)
 27 =   -1.41755 : (13.031774258126534, 16.998201908838432)
 28 =   -1.41768 : (13.000662635503582, 17.001139499308035)
 29 =   -1.42032 : (13.013283510833627, 17.010351129505967)
 30 =   -1.41974 : (13.006269036099196, 17.006911318399233)
 31 =   -1.42020 : (12.99988527317096, 17.002407683519927)
 32 =   -1.42032 : (13.003784345977426, 17.002644369625095)
 33 =   -1.42030 : (13.001833848997947, 17.002543083897233)
-------
itk::simple::TranslationTransform
 TranslationTransform (0x55ecc01b8c90)
   RTTI typeinfo:   itk::TranslationTransform<double, 2u>
   Reference Count: 2
   Modified Time: 13966
   Debug: Off
   Object Name: 
   Observers: 
     none
   Offset: [13.0018, 17.0025]

Optimizer stop condition: RegularStepGradientDescentOptimizerv4: Step too small after 34 iterations. Current step (0.000976562) is less than minimum step (0.001).
 Iteration: 35
 Metric value: -1.4203230355054985

Input Images

Fixed Image

Moving Image

Output Image

Composition Image

Code

PythonR

#!/usr/bin/env python

""" A SimpleITK example demonstrating image registration using Mattes mutual
    information as the metric. """

import sys
import os
import SimpleITK as sitk


def command_iteration(method):
    """ Callback invoked when the optimization has an iteration. """
    print(
        f"{method.GetOptimizerIteration():3} "
        + f"= {method.GetMetricValue():10.5f} "
        + f": {method.GetOptimizerPosition()}"
    )


def main(args):
    """ A SimpleITK example demonstrating image registration using Mattes mutual
        information as the metric. """

    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)