Different Rotation Matrices between Euler3DTransform and Scipy Rotation

Beginner Questions
,
1

Hello,
I’m trying to rotate sections of my images according to certain orientations. However, I noticed some discrepancies between the rotations of the Euler3DTransform object and scipy.spatial.transform.Rotation when given the same set of Euler Angles.

The following code reproduces problem, with the output just below. Given a set of Euler angles in XYZ, the resulting rotation matrices between what sitk and scipy differ, especially with rotations around both X and Y. I used the metric (from the cited paper) to compare the rotation matrices and print them when they differ significantly. It’s notable in the first error (rotation in only X and Y), that the matrices only differ in position. Also worth mentioning that this error is not unique to using angles of 0.5 radians, I only chose it for demonstration purposes.

import numpy as np
from scipy.spatial.transform import Rotation as R
from SimpleITK import Euler3DTransform

euler_angles = [  # euler angles 'xyz' in radians
    (0, 0, 0),
    (0.5, 0, 0),
    (0, 0.5, 0),
    (0, 0, 0.5),
    (0.5, 0.5, 0),
    (0.5, 0, 0.5),
    (0, 0.5, 0.5),
    (0.5, 0.5, 0.5)
]

t = Euler3DTransform()
for ea in euler_angles:
    t.SetRotation(*ea)

    msitk = np.array(t.GetMatrix()).reshape(3,3)
    msci = R.from_euler("xyz", ea).as_matrix()

    # Provides a metric for comparing two rotation matrices, results in values
    # in the range of[0, 2√2]
    # Huynh, D. Q. (2009). Metrics for 3D Rotations: Comparison and Analysis.
    # Journal of Mathematical Imaging and Vision, 35(2), 155–164. https://doi.org/10.1007/s10851-009-0161-2
    error = np.linalg.norm(np.identity(3) - msitk @ msci.T)

    print(f"{ea}: error: {error:.4f}")
    if error > 1e-4:
        print("sitk:")
        print(msitk)
        print()
        print("scipy:")
        print(msci)

Output:

(0, 0, 0): error: 0.0000
(0.5, 0, 0): error: 0.0000
(0, 0.5, 0): error: 0.0000
(0, 0, 0.5): error: 0.0000
(0.5, 0.5, 0): error: 0.3456
sitk:
[[ 0.87758256  0.          0.47942554]
 [ 0.22984885  0.87758256 -0.42073549]
 [-0.42073549  0.47942554  0.77015115]]

scipy:
[[ 0.87758256  0.22984885  0.42073549]
 [ 0.          0.87758256 -0.47942554]
 [-0.47942554  0.42073549  0.77015115]]
(0.5, 0, 0.5): error: 0.0000
(0, 0.5, 0.5): error: 0.0000
(0.5, 0.5, 0.5): error: 0.3456
sitk:
[[ 0.65995575 -0.42073549  0.62244683]
 [ 0.62244683  0.77015115 -0.13938128]
 [-0.42073549  0.47942554  0.77015115]]

scipy:
[[ 0.77015115 -0.21902415  0.59907898]
 [ 0.42073549  0.88034656 -0.21902415]
 [-0.47942554  0.42073549  0.77015115]]

I’d label something like this as a bug somewhere, but I’d be very thankful to anyone who has a different explanation for this phenomenon.

Edit:
The correctly resulting rotation matrices from the first four sets of angles should showcase that it shouldn’t have anything to do with an incorrect order of Euler angles.

2

The resulting matrix depends on the order of rotations around axes, see documentation of SetComputeZYX(). Namely:

This class supports two of them, ZXY and ZYX. The default is ZXY.

This is the most likely reason for the discrepancy.

3

Thank you for your reply. While I noticed this, now that I think about it again, I start to see what type of error it is. If I understand correctly, while I give the transform the rotations in X, Y and Z, internally they are computed in order ZXY (or ZYX, depending on the mode), and thus result in a different rotation matrix.

However, I just tried again by using

t.ComputeZYXOn()

after initializing the object and then creating the Rotation object such as.

msci = R.from_euler("zyx", ea[::-1]).as_matrix()

Now the order of operation should be the same, but none except the zero and single-axis match. I turned on ZYX specifically so that I can simply reverse the given angles, but even by using ZXY and adjusting the order the results aren’t any better. I can’t seem to figure out how to make the two compatible.

EDIT: I’m not sure if double-posting is allowed here, so I’ll add this as an edit. This is the crux of where my confusion comes from. Please observe the following rotation matrices when given the two different modes.

t = Euler3DTransform()
ea = (0.2, 0.3, 0)
t.ComputeZYXOff()  # Default
t.SetRotation(*ea)

msitk = np.array(t.GetMatrix()).reshape(3,3)
t.ComputeZYXOn()
msitk2 = np.array(t.GetMatrix()).reshape(3,3)
msci = R.from_euler("zyx", ea[::-1]).as_matrix()

print(msitk)
print(msitk2)
print(msci)

results in:

msitk:
[[ 0.95533649  0.          0.29552021]
 [ 0.0587108   0.98006658 -0.18979606]
 [-0.28962948  0.19866933  0.93629336]]

msitk2:
[[ 0.95533649  0.0587108   0.28962948]
 [ 0.          0.98006658 -0.19866933]
 [-0.29552021  0.18979606  0.93629336]]

msci:
[[ 0.95533649  0.          0.29552021]
 [ 0.0587108   0.98006658 -0.18979606]
 [-0.28962948  0.19866933  0.93629336]]

notice how the matrix in the ZXY mode is identical to the scipy one in ZYX mode, and not the other way around.

4

This is contrary to your conclusion from the first post :slight_smile:

Maybe add t.ComputeZYXOff() to your original script? Or some other sanity checks?

“Double posting” is allowed here, and is preferred to editing existing posts. This is for convenience of those who follow the forum via email.

5

I kept wrecking my head over this, I finally found the root of the problem: I didn’t use the scipy module correctly; I mixed up intrinsic vs extrinsic Euler rotations, so when using it with ‘ZXY’ vs ‘zxy’ I do get the expected results.
In hindsight this is was a very avoidable mistake, but I think I learned quite a bit about how either module works in the process :sweat_smile:

1 Like
6

can I see the right code of this ?