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.34810 : (0.8170900112373908, 0.5765101157274526)
  1 =   -0.35882 : (1.421157461081397, 1.3734433099407186)
  2 =   -0.37026 : (1.9085426614258152, 2.246630384103656)
  3 =   -0.38151 : (2.751632636580906, 2.7844029138901716)
  4 =   -0.39196 : (3.425476052356483, 3.5232770919583345)
  5 =   -0.40176 : (4.199705580768797, 4.156182011621979)
  6 =   -0.40866 : (4.623283812373065, 5.06204153810001)
  7 =   -0.42017 : (5.352553529490031, 5.746267874699339)
  8 =   -0.43541 : (6.206567034025554, 6.266518712344945)
  9 =   -0.45311 : (7.113473512231072, 6.687850695214587)
 10 =   -0.46605 : (7.857789435425893, 7.355678372447749)
 11 =   -0.48649 : (8.508110515967779, 8.115337830496443)
 12 =   -0.51162 : (9.273949195192236, 8.758370580365172)
 13 =   -0.52496 : (9.834007869996828, 9.586823524578737)
 14 =   -0.54537 : (10.486715043369621, 10.344433810826025)
 15 =   -0.58103 : (11.173457743087132, 11.071334398516465)
 16 =   -0.61302 : (11.567265930282696, 11.990527039749924)
 17 =   -0.65810 : (12.129882374809362, 12.81724509105168)
 18 =   -0.70569 : (12.569575578930202, 13.715393122419533)
 19 =   -0.76771 : (12.825299870965091, 14.682142876214824)
 20 =   -0.85076 : (12.941540544813122, 15.675363952180973)
 21 =   -0.99125 : (13.041467617599132, 16.670358716048813)
 22 =   -1.29344 : (12.939224728448329, 17.665118180250133)
 23 =   -1.15414 : (13.080738488662645, 17.185562286041748)
 24 =   -1.35171 : (12.827738681852484, 16.754295433902588)
 25 =   -1.29620 : (13.004704288293269, 16.930883015960458)
 26 =   -1.40686 : (12.981453483350592, 17.179799467946896)
 27 =   -1.36363 : (13.004169364134526, 17.056880841367317)
 28 =   -1.40909 : (12.995202290853854, 16.93220288984108)
 29 =   -1.40701 : (13.002360131701693, 16.99429165947293)
 30 =   -1.41638 : (12.964636192925393, 17.044122919172057)
 31 =   -1.40767 : (12.98824305551436, 17.02364662705884)
 32 =   -1.41433 : (13.00625553581766, 16.998110122730267)
 33 =   -1.41637 : (12.990630541619018, 16.9980966582441)
 34 =   -1.41596 : (12.99844285176646, 16.998042193483332)
 35 =   -1.41634 : (13.006255082343644, 16.998107075367766)
 36 =   -1.41637 : (13.00234883284669, 16.9981050929318)
 37 =   -1.41640 : (12.998442667812787, 16.998079328801463)
 38 =   -1.41634 : (13.000395773195743, 16.99807057501985)
-------
itk::simple::TranslationTransform
 TranslationTransform (0x55a1a77b1c90)
   RTTI typeinfo:   itk::TranslationTransform<double, 2u>
   Reference Count: 2
   Modified Time: 13981
   Debug: Off
   Object Name: 
   Observers: 
     none
   Offset: [13.0004, 16.9981]

Optimizer stop condition: RegularStepGradientDescentOptimizerv4: Gradient magnitude tolerance met after 39 iterations. Gradient magnitude (7.39287e-05) is less than gradient magnitude tolerance (0.0001).
 Iteration: 40
 Metric value: -1.416386053019758

Input Images

_images/ImageRegistrationMethod4_fixed.png

Fixed Image
_images/ImageRegistrationMethod4_moving.png

Moving Image

Output Image

_images/ImageRegistrationMethod4_composition.png

Composition 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__":