Coregistration of two rectangles using SimpleITK

Algorithms
1

Hello,
I am trying to co-register two very simple 2D png images of rectangles (white rectangle value 256 on black (0 value) background), but struggling with getting them aligned no matter what metric, optimizer or parameters I choose. Would anybody please be so kind and review the code below and give me a hint about what might be wrong with this code/parameters, etc…?
Thanks a lot!!
Images are attached.

FixedIdx
FixedIdx1890×1417 21.1 KB

MovingIdx
MovingIdx1890×1417 20.7 KB

Here is the code

import SimpleITK as sitk
import sys
import os
import matplotlib.pyplot as plt 
import matplotlib as mpl
import numpy as np

def command_iteration(method):
    print(
        f"{method.GetOptimizerIteration():3} "
        + f"= {method.GetMetricValue():10.5f} "
        + f": {method.GetOptimizerPosition()}"
    )

FixedPath =  r"D:\...path\FixedIdx.png"
MovingPath =  r"D:\...path\ImageOverlayDevel\MovingIdx.png"

fixed = sitk.ReadImage(FixedPath, sitk.sitkFloat32)
moving = sitk.ReadImage(MovingPath, sitk.sitkFloat32)

fixed = sitk.RescaleIntensity(fixed)
moving = sitk.RescaleIntensity(moving)

R = sitk.ImageRegistrationMethod()
R.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)
#R.SetOptimizerAsGradientDescent(learningRate=5000.0, numberOfIterations=100, convergenceMinimumValue=1e-7, convergenceWindowSize=10)
R.SetOptimizerAsRegularStepGradientDescent(learningRate=150.0, minStep=0.001, numberOfIterations=500, gradientMagnitudeTolerance=1e-7)

# Set initial translation offset (large enough)
initial_transform = sitk.TranslationTransform(2)
initial_transform.SetOffset((0, 0))
R.SetInitialTransform(initial_transform, inPlace=False)

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)

showimg = sitk.GetArrayFromImage(moving+fixed)

imgplot2 = plt.imshow(showimg )
2

With such a sparse image, you might need full sampling, not sub-sampling percentage.

3

Hi Dzenan,

thanks for the idea. I assume you mean I should do this?

R.SetMetricSamplingStrategy(R.NONE)
R.SetMetricSamplingPercentage(1.0)

This doesn’t change the situation, unfortunately. The final overlap still doesn’t look better…

image
image534×402 16.8 KB

1 Like
4

Your learning rate seems suspiciously large, and minimum step rather small. You should probably also increase number of iterations.

With such images, you are probably a lot better off creating a signed distance field (e.g. Maurer filter), and register those instead.

5

You may want to try using a large gaussian to increase the capture radius of the metric space. There should be a SimpleITK Notebooks which demonstrates rendering this metric space: GitHub - InsightSoftwareConsortium/SimpleITK-Notebooks: Jupyter notebooks for learning how to use SimpleITK

You could also try using an FFT based method to find the optimization translation:

2 Likes
6

Hello @Andre,

Welcome to the world of image registration! You have successfully created a corner case as your initial input. Happens to the best of us.

TL;DR

Convert your binary (i>0 or i==255 whichever way you want to binarize them) images to distance maps and register those using the mean squares similarity metric and basic registration settings, no need for any advanced settings and it will work:

import SimpleITK as sitk

FixedPath = "FixedIdx.png"
MovingPath = "MovingIdx.png"

fixed_image = sitk.ReadImage(FixedPath, sitk.sitkFloat32)
moving_image = sitk.ReadImage(MovingPath, sitk.sitkFloat32)

fixed_image_distance_map = sitk.SignedMaurerDistanceMap(
    fixed_image > 0, squaredDistance=False, useImageSpacing=True
)
moving_image_distance_map = sitk.SignedMaurerDistanceMap(
    moving_image > 0, squaredDistance=False, useImageSpacing=True
)

transform = sitk.TranslationTransform(2)
print(f"initial translation: {transform.GetOffset()}")

registration_method = sitk.ImageRegistrationMethod()

# Similarity metric settings.
registration_method.SetMetricAsMeanSquares()
registration_method.SetMetricSamplingStrategy(registration_method.RANDOM)
registration_method.SetMetricSamplingPercentage(0.01)

registration_method.SetInterpolator(sitk.sitkLinear)

# Optimizer settings.
registration_method.SetOptimizerAsGradientDescent(
    learningRate=1.0,
    numberOfIterations=100,
    convergenceMinimumValue=1e-6,
    convergenceWindowSize=10,
)

registration_method.SetInitialTransform(transform)

registration_method.Execute(fixed_image_distance_map, moving_image_distance_map)

# Always check the reason optimization terminated.
print("Final metric value: {0}".format(registration_method.GetMetricValue()))
print(
    "Optimizer's stopping condition, {0}".format(
        registration_method.GetOptimizerStopConditionDescription()
    )
)

print(f"final translation: {transform.GetOffset()}")

The details

An image is a dense regular grid sampling of a continuous signal. When the original signal occupies a small region in the image, using intensity based registration to align the two images is likely to fail.

This is particularly an issue when the intensity values in the object do not vary (two contours with single intensity value). In such cases the image-to-image registration is actually a less appropriate approach than a contour-to-contour / surface-to-surface registration. So when you have an image representation of a contour/surface you can:

  1. Extract the contour/surface and perform registration using classics such as the Iterative Closest Point algorithm (available in ITK, not in SimpleITK).
  2. Convert the sparse representation to a distance map as done above and solve it using intensity based registration.
2 Likes
7

Thanks Bradley,

Could you specify what you mean by “using a large Gaussian” ?

You mean filtering the image by large radius Gauss filter, with the aim to spread visible values of the sparse image?

Thanks heaps (sorry for a possibly dumb question)

8

Amazing!

Thanks for the code and explanation. That really helped. Works perfectly now.

Just a quick question about corner case…

I had been contemplating about simply using fiducial registration of the 4 corners, as an alternative.

Is this implemented in simpleITK ? Any example code for me to learn?

Thanks again.

9

Hello @Andre,

Yes, fiducial registration is available in SimpleITK, goes by the name LandmarkBasedTransformInitializerFilter.

As this appears to be your initial foray into registration, I highly recommend skimming the SimpleITK notebooks repository, specifically the 6* series which deal with registration and illustrate various aspects of the task and how to address them. Specifically, go over the registration initialization and the registration errors, notebooks. The concepts described there are general and will serve you well no matter what tool you use for registration.

For a concise overview of the toolkit see the tutorial.

2 Likes
10

Yes exactly. Graphic the metric value over points and the effects of performing gaussian smoothing on on or two of this images is a enlightening exercise.