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:
- Is my code correct to do a rigid registration with possible translations and/or rotations, with a given pixel size?
- I expected “Translation” or “Offset” to be [15, 15, -15]. I actually don’t understand the difference between these two.
- 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)
- 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”.
- 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.