A quick question concerning applying translation and rotation to volumes and then registering them to compare the results to the imposed motion:
I have a python script that uses SimpleITK to impose random motion on a series of phantom scans. The script outputs the values of translation and rotation to an excel file for later comparison. The upper limit for translation is 5mm and for rotation 0.1rad. All values are printed out and well below those maximum values.
I run a 3D Euler Rigid Registration, implemented in C++ on the reference volume and the volumes where motion was imposed. The output values for the transformation file are completely different than the imposed motion ones. Also, when I use ITK-Snap to perform a registration, the results of the registration are way over the upper limit specified in the python script.
For example:
imposed translation motion output on x-axis: 0,810653599 mm
registration result translation on x axis with ITK-Snap: 42,32 mm
Is there a resource to read on why all the values are so different? Does the confusion stem from the usage of the virtual and physical images within the registration framework and do I have to translate the values between the two domains?
Based on your description it is hard to tell what is going on. Is the ITK-SNAP registration correct? Is the SimpleITK transformation applied correctly to the images?
Let’s start with a simple setup, known, arbitrary, SimpleITK translation:
Load each image in ITK-SNAP and place the cursor at the same anatomical landmark. Open the tools->image information, info tab and see the “3D cursor coordinates, world units ITK”. Confirm that the landmark was translated as expected. Then run the ITK-SNAP registration and see if you get the expected translation in x-direction (choice of moving/fixed image will yield the expected or inverse transformation depending on the image roles).
The ITK-SNAP registration is visually absolutely correct and aligns the images as expected. This is also confirmed by running your code snippet with a volume of the datasets:
Also, checking the landmarks shows that the translation of 20 results in a shift of 10 with a voxel size of ~2 and a correct shift of 20 in “world units ITK” This sounds good to me. Would a snippet of the script which imposes motion help for further investigation?
I think we can assume that something is wrong in the python script, which imposes motion in that case, right?
Yes, based on your results I’d suspect the script used for resampling (imposing motion) isn’t doing what you want it to do. Please share any code you can so we can further help.
It’s not my script, and I am not too familiar with python and SimpleITK, but the following parts of code is where I suspect that maybe something could be wrong:
Function to impose motion:
# Function to impose predefined motion on to the 3D volumes
def impose_motion(itk_img, translation, rotation, fpath_mtn_impsd_img):
dim = itk_img.GetDimension()
translation_transform = sitk.TranslationTransform(dim, translation)
rigid_euler = sitk.Euler3DTransform()
rigid_euler.SetTranslation(translation_transform.GetOffset())
rigid_euler.SetRotation(*rotation)
mtn_impsd_img = sitk.Resample(
itk_img, itk_img, rigid_euler, sitk.sitkNearestNeighbor, 0, itk_img.GetPixelID()
)
print("Translation in function: "+
str(rigid_euler.GetTranslation()))
print("Rotation in function: "+
str(rigid_euler.GetAngleX())+
"/"+
str(rigid_euler.GetAngleY())+
"/"+
str(rigid_euler.GetAngleZ()))
sitk.WriteImage(mtn_impsd_img, fpath_mtn_impsd_img)
Code where for every file, a random value set is generated and passed to the impose motion function (nested within a for each loop):
The impose_motion is a bit strange: (1) what is the role of translation_transform, doesn’t seem to be required, use the translation parameter directly. (2) If interpolating real valued, anatomical, images, the go-to interpolation method would be sitkLinear and not sitkNearestNeighbor (use nearest neighbor if interpolating a discrete label image).
There appears to be a bug in the invocation of this function:
The itk_image is read from the file given by img_to_mtn_imps and the resampled image is written to the file by the same name, last parameter passed to impose_motion. Is overwriting the original file intentional?
Could it just be that the rotation has such an impact on the translation that the raw, generated, values are simply not valid anymore, or can’t be taken as a reference?
Because now the rotation seems to be correct? 0.3 rad should be ~17.19° around the z axis.
OK, so the problem is with the understanding of how the transforms work. By default rotation is around the world point (0,0,0), this is an arbitrary point relative to the image. It is not the origin of the image or its center. What you want is to rotate around the image center. Add the following line after setting the rotation:
transform.SetCenter(image.TransformContinuousIndexToPhysicalPoint([(sz-1)/2 for sz in image.GetSize()]))
Thank you @zivy for pointing that out. I thought SimpleITK automatically assumes the center of the rotation as the image center. I will gladly read through your provided resources!
One last question, as I just applied your suggested center-shift:
After the shift on x by 20, and rotation around z by 0.3 is applied, is there a way to confirm the registration just by the values? So for example, registering the images with ITK-Snap results in:
Having a 5.55 shift on the y-axis and 18.85 on the x-axis, instead of just 20 on x. Is it too naive or wrong to assume that the values should be the inverse ones of the imposed motion because there are more ways to transfer back the volume for a successful registration?
I am aware of other registration evaluation methods, I just ask out of curiosity and for understanding the topic a little better.
Unfortunately, there is no simple answer here. Key theoretical concepts underlying spatial transformations:
There are multiple parametric representations for the same spatial transformation and the parameterizations are not necessarily unique. Limiting ourselves to rigid 3D transformations, two ways to specify the rotation in SimpleITK are Euler angles or a unit quaternion (a.k.a. versor). Let’s assume we chose Euler angles, there are twelve possible angle combinations out of which SimpleITK supports ZXY and ZYX. Finally, using the default ZXY parameterization we have two parameterizations which represent the same rotation: T(0,0,90) = T(180,180,270) = \left[\begin{array}{ccc} 0 & -1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{array}\right]
The registration error is spatially varying, so judging whether a registration is sufficiently accurate is dependent on the spatial region of interest and the associated task.
Looking at the parameter values can serve as a sanity check but not much more. The print_transformation_differences found in the Transformations notebook may be of interest. This notebook in the SimpleITK tutorial illustrates various ways of evaluating the registration (from summary stats to histogram, to a spatial visualization).
Thank you so much for your time and explanations. It cleared up a lot of things for me - very, very helpful.
The issue with the discrepancy between imposed motion and registration result was mainly resolved by setting the center of rotation accordingly.
Also thank you for this information, I will gladly look through the notebooks. Especially print_transformation_differences sounds very helpful on a first glance.
