Struggling to do simple 3D rigid registration

Hello,

I’m am struggling to perform a simple 3D rigid registration between 2 images. I have started with one of the tutorials, and tried to adapt it to my situation. Particularly I want to adjust it to do a rigid registration (i.e. including rotations), but I can’t get it to work.

I have taken a PET brain image, manually transformed it by [15, 15, -15] mm (no rotations for now), and then I’m trying to register it to the original image. If I set the SetInitialTransform to be TranslationTransform then it works as expected. However if I set it to be Euler3DTransform then it doesn’t seem to work and I don’t understand the results.

Here is my code:

#!/usr/bin/env python
from __future__ import print_function
import SimpleITK as sitk
import sys
import os

def command_iteration(method):
    print("{0:3} = {1:10.5g} : {2}".format(method.GetOptimizerIteration(),
                                           method.GetMetricValue(),
                                           method.GetOptimizerPosition()))

if len(sys.argv) < 4:
    print("Usage:", sys.argv[0], "<fixedImageFilter> <movingImageFile>",
          "<outputTransformFile>")
    sys.exit(1)

fixed = sitk.ReadImage(sys.argv[1], sitk.sitkFloat32)
moving = sitk.ReadImage(sys.argv[2], sitk.sitkFloat32)

# image size:   [75, 75, 73] pixels
# pixel size:   [4, 4, 2.78] mm
# image extent: [300, 300, 202.94] mm

fixed.SetSpacing((4.0, 4.0, 2.78));     # pixel size in mm
moving.SetSpacing((4.0, 4.0, 2.78));

# fixed.SetOrigin((150, 150, 101.47));
# moving.SetOrigin((150, 150, 101.47));

R = sitk.ImageRegistrationMethod()
R.SetMetricAsMeanSquares()
R.SetOptimizerAsRegularStepGradientDescent(4.0, 0.001, 50, 0.5, 1e-7)

tx = sitk.CenteredTransformInitializer(fixed, moving, sitk.Euler3DTransform(),
    sitk.CenteredTransformInitializerFilter.GEOMETRY)
R.SetInitialTransform(tx)
# 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("Optimizer stop condition: {0}"
      .format(R.GetOptimizerStopConditionDescription()))
print(" Iteration: {0}".format(R.GetOptimizerIteration()))
print(" Metric value: {0}".format(R.GetMetricValue()))

sitk.WriteTransform(outTx, sys.argv[3])

And the output is:

itk::simple::Transform
 Euler3DTransform (0x7fc363631460)
   RTTI typeinfo:   itk::Euler3DTransform<double>
   Reference Count: 3
   Modified Time: 2401
   Debug: Off
   Object Name: 
   Observers: 
     none
   Matrix: 
     0.968088 -0.235134 -0.0867007 
     0.203394 0.939297 -0.276318 
     0.14641 0.249866 0.957148 
   Offset: [49.8245, 7.73086, -55.7376]
   Center: [148, 148, 100.08]
   Translation: [1.62477, 1.19532, -1.37756]
   Inverse: 
     0.968088 0.203394 0.14641 
     -0.235134 0.939297 0.249866 
     -0.0867007 -0.276318 0.957148 
   Singular: 0
   Euler's angles: AngleX=2.88905 AngleY=-3.29338 AngleZ=3.38688
   m_ComputeZYX = 0

Optimizer stop condition: RegularStepGradientDescentOptimizerv4: Maximum number of iterations (50) exceeded.
 Iteration: 50
 Metric value: 3.75348540451e-05

Here are my questions:

  1. Is my code correct to do a rigid registration with possible translations and/or rotations, with a given pixel size?
  2. I expected “Translation” or “Offset” to be [15, 15, -15]. I actually don’t understand the difference between these two.
  3. Should I set SetOrigin((150, 150, 101.47)) for the images before the registration? (I would like to do this because in a real case I will usually have an offset center)
  4. What is “Center” in the output? I imagine it’s the estimated image center, but if I SetOrigin as above then it just gets added to this “Center”.
  5. Are the “Euler’s angles” the rotational part of the registration parameters? If so I expect them to be near 0 or 2*pi in this case, but they are near pi.

I’m pretty sure that I am just missing something trivial but I could not figure it out after searching the examples, tutorials, documentation, and this forum.

Many thanks for your help.

Hello @msb,

Before you go any farther, please go over the following resources (possibly you missed some of the relevant information):

  1. Fundamental concepts.
  2. Registration overview.
  3. Set of registration notebooks (6* series) on the SimpleITK Jupyter notebook repository.

With respect to your registration:

  1. Why are you manually modifying the image spacing after reading? This is very unusual (this information is always obtained directly from the image when it is read).
  2. You have premature termination (stopping condition is “Maximum number of iterations (50) exceeded.”).
  3. The nomenclature of “Center”, “Offset”, “Translation” that appear in the printed transformation are described in the fundamental concepts page referenced above and are specific to ITK/SimpleITK.
2 Likes

Hello @zivy,

Thank you for your reply. I reviewed those documents and links again in more detail and they did indeed help clarify some things that I missed when I went through them initially.

To answer your first question: I am manually modifying the image parameters because those parameters are not present in the header information. I saved this images from Matlab as nifti images and creating a nifti header from scratch is not trivial in Matlab, so I just edited the relevant values in Python.

I now know how the image origin should be set. I set the image parameters like this:

# image size:   [75, 75, 73] pixels
# pixel size:   [4, 4, 2.78] mm
# image extent: [300, 300, 202.94] mm
fixed.SetSpacing((4.0, 4.0, 2.78));     # pixel size in mm
moving.SetSpacing((4.0, 4.0, 2.78));
fixed.SetOrigin((-148, -148, -100.08));
moving.SetOrigin((-148, -148, -100.08));

Setting the origin to be half a pixel in from the image extent. The image direction vectors are the identity, which should be correct.

I know that I have premature termination. Below I ran it again until convergence (at 750 iterations), but the result is still incorrect. Being quite a simple registration problem I believe that I must be doing something else wrong.

I set up my registration like this:

R = sitk.ImageRegistrationMethod()
R.SetMetricAsMeanSquares()
R.SetOptimizerAsRegularStepGradientDescent(4.0, 0.01, 2000, 0.5, 1e-6)
R.SetInitialTransform(sitk.CenteredTransformInitializer(fixed, moving,
    sitk.Euler3DTransform(), sitk.CenteredTransformInitializerFilter.GEOMETRY))
R.SetInterpolator(sitk.sitkLinear)

And the result is

itk::simple::Transform
 Euler3DTransform (0x7f838ab99f50)
   RTTI typeinfo:   itk::Euler3DTransform<double>
   Reference Count: 2
   Modified Time: 5226
   Debug: Off
   Object Name: 
   Observers: 
     none
   Matrix: 
     0.994857 -0.0991812 -0.0205532 
     0.0967002 0.990413 -0.0986452 
     0.0301399 0.0961504 0.99491 
   Offset: [9.50359, 7.69879, -9.5203]
   Center: [0, 0, 0]
   Translation: [9.50359, 7.69879, -9.5203]
   Inverse: 
     0.994857 0.0967002 0.0301399 
     -0.0991812 0.990413 0.0961504 
     -0.0205532 -0.0986452 0.99491 
   Singular: 0
   Euler's angles: AngleX=3.04529 AngleY=-3.17188 AngleZ=3.2414
   m_ComputeZYX = 0

Optimizer stop condition: RegularStepGradientDescentOptimizerv4: Step too small after 751 iterations. Current step (0.0078125) is less than minimum step (0.01).
 Iteration: 752
 Metric value: 8.91753355411e-06

I believe that using sitk.CenteredTransformInitializerFilter.GEOMETRY is the right choice since this is a PET image being registered to a slightly shifted version of itself. I expect the registration to return a transformation of [15, 15, -15] and angles of near [0, 0, 0].

Can you possibly see if I am doing something wrong? Or can you suggest a simple case that I can test with?

Many thanks!

Hello @msb,

First step in debugging a registration is to visualize the initial condition after you used the CenteredTransformInitializer.

Second look at your code, it seems you are not calling R.SetOptimizerScalesFromPhysicalShift() this is a critical step if your transformation is rigid as the units of the rotation parameters and translation parameters are not commensurate (a unit change in rotation is not equivalent to a unit change in translation). If you are only doing translation then you don’t need to call the method because all optimized parameters are the same.

In your case you were describing the problem as a translation only, but in your code you use Euler3DTransform.

2 Likes

Ha! It was the call to SetOptimizerScalesFromPhysicalShift()! That has fixed it it seems.

Thanks very much for your help @zivy.

3 Likes