Understanding 3DVersor rotation output

1

Hi!

I’m working on a project where I use a Versor3DRigidTransform for registering 3D MRI images, and I’m having some issues interpreting the output from the registration results. My end goal is to get a transformation which is equivalent to a centered rotation (as in the middle of the image) and a translation, this is because I’m applying image corrections in the frequency domain (i.e. MRI k-space data).

I’ve generated test data with a 3D Shepp Logan phantom with a 45 deg rotation between the fixed and moving image. The images are generated outside of ITK and I use itk.image_from_array to create an image object, where I specify the centre of the image as the origin.

spacing = [1,1,1]
img = itk.image_from_array(np.abs(np.ascontiguousarray(img_array)))
img.SetSpacing([float(x) for x in spacing])
img.SetOrigin([-img_array.shape[i]/2*spacing[i] for i in range(3)])

I set up an initialiser with .GeometryOn() to ensure the image origin is used for the transformation.

After running the registration I resample the moving image to the fixed space and I get a perfect overlap, but when I look at the estimated registration parameters there is a discrepancy in the rotation estimate between the versor rotation and the rotation matrix which I dont understand.

The final parameters from the registration are:

Rotation: (-0.04, 0.01, -21.93) deg
Translation: (-0.39, 0.94, 0.00) mm

But when I pull out the rotation matrix and the offset using .GetMatrix() and .GetOffset() I get

[-8.49832563e-04 -9.54532133e-04 -1.09447133e-05]

[[ 7.07122107e-01  7.07090904e-01  8.82676813e-04]
 [-7.07091446e-01  7.07121348e-01  1.04216899e-03]
 [ 1.12748595e-04 -1.36107396e-03  9.99999067e-01]]

The offset is basically 0, which I was hoping since I set the origin to be at iso-centre. But the rotation matrix equates to the following rotations (if my math is correct):

Rotation (x,y,z): -0.078, -0.00646, -45.0 degrees

And this is what I would expect.

Can someone explain why there is this difference in the outputs?

Please let me know if you need more information, I could go on and on about the details but tried to make it as brief as possible.

Thanks!

2

Probably there is an error in converting radians to degrees in this message. A division by 360 instead of 180, or by 180 instead of 90? If it is in ITK, it would be good if you made a PR with the fix.

3

Thanks! Now I clearly see that there is a factor 2 between the values, I get the same when I run it with other parameters. My first thought was also then to check for a conversion error on my side but I can’t find anything in my code. I’m afraid that I don’t really know where to start looking in the code, I haven’t navigated through the ITK source code before. Would you be able to give me some pointers?

Thanks!

4

Search "Rotation:" (22 hits in 15 files of 11043 searched). But I examined the results and the format does not match what you show. I don’t think this message comes from ITK.

5

Hi again, after some discussion with a colleague and more thorough reading of the ITK manual, we realised that the output parameters of the 3D Versor transform is a 3D versor, and those angles do not represent the rotation as Euler angles which I wrongly had interpreted them as. Since versors uses unit quarterions, the factor 1/2 difference between the versor and rotation matrix angle is expected (although I have to admit my quarterion math isn’t good enough to fully grasp that).

For those who might run into a similar confusion about the parameters, there are more people out there wondering about that half-angle in quarterions (Concise description of why rotation quaternions use half the angle - Mathematics Stack Exchange)

3 Likes
6

I just wanted to add a further comment here to anyone who might pick up this post in the future. The statement I made above about the factor of 0.5 being the expected difference is wrong.

After further reading in the ITK manual and Wikipedia I’ve solved this confusion. The Versor transform returns the vector part of a unit quarterion. The real part (sometimes first or last component) is scaled internally to get a unit quarterion. Below is an example of how (to my understanding) get the intrinsic Euler angles from the parameters of the versor transform.

With a 3D rotation matrix R of an intrinsic rotation, here set to R=R_z(\alpha=15^\circ)\cdot R_y(\alpha=10^\circ)\cdot R_x(\alpha=5^\circ)

print(R)
>>> [[ 0.98106026 -0.03941355  0.18965056]
 [ 0.08583165  0.96616727 -0.24321542]
 [-0.17364818  0.254887    0.95125124]]

# Create a versor transform and set input matrix
TransformType = itk.VersorRigid3DTransform[itk.D]
ImageType = itk.Image[itk.D, 3]
transform = TransformType.New()
transform.SetMatrix(itk.matrix_from_array(R_in))

# Get the parameters and calculate the 4th component
par_out = transform.GetParameters()
q1, q2, q3 = np.array(par_out)[0:3]
q0 = np.sqrt(1 - q1**2 - q2**2 - q3**2)

A rotation matrix can be calculated from the quarterion as

R_quart = np.array([[1-2*(q2**2 + q3**2), 2*(q1*q2 - q0*q3), 2*(q0*q2+q1*q3)],
             [2*(q1*q2+q0*q3), 1-2*(q1**2+q3**2), 2*(q2*q3-q0*q1)],
             [2*(q1*q3-q0*q2), 2*(q0*q1+q2*q3), q0**2-q1**2-q2**2+q3**2]])
print(R_quart)

>>> [[ 0.98106026 -0.03941355  0.18965056]
 [ 0.08583165  0.96616727 -0.24321542]
 [-0.17364818  0.254887    0.95125124]]

Which we can see is the same as the input matrix above.

The conversion to Euler angles is then done by:

rz = np.rad2deg(np.arctan2(2*(q0*q1+q2*q3),(1-2*(q1**2+q2**2))))
ry = np.rad2deg(np.arcsin(2*(q0*q2 - q3*q1)))
rx = np.rad2deg(np.arctan2(2*(q0*q3+q1*q2), 1-2*(q2**2+q3**2)))
print("Quarterion conversion: {}".format([rx, ry, rz]))

>>> Quarterion conversion: [15.000000000000002, 10.0, 5.0]

which is exactly what we input. To further clarify, in contradiction to my post above, this is simply not twice the value of the versor

print("Quarterion in deg: {}".format(np.rad2deg([q0,q1,q2,q3])))

>>> Quarterion in deg: [56.56401433  7.22709397  5.27119322  1.81721431]

For context, I need to get the rotation around each axis as I’m doing motion correction for a long time series of images and want to plot the rotation as a function of time, so getting this conversion is key for that. My understanding is that it is better to do the conversion from quarterion to euler angles, rather than from the rotation matrix as there can be some ambiguities in that for large rotations, but I might be wrong.

Cheers!

4 Likes