CODEWITHSUNDEEP
pythonscipyeuler-anglesrotational-matricessimpleitk
Peter
Peter

Reputation: 3669

Difference between SimpleITK.Euler3DTransform and scipy.spatial.transform.Rotation.from_euler?

Using these two library functions:

to create a simple rotation matrix from Euler Angles:

import numpy as np
import SimpleITK as sitk
from scipy.spatial.transform import Rotation
from math import pi

euler_angles = [pi / 10, pi / 18, pi / 36]

sitk_matrix = sitk.Euler3DTransform((0, 0, 0), *euler_angles).GetMatrix()
sitk_matrix = np.array(sitk_matrix).reshape((3,3))
print(np.array_str(sitk_matrix, precision=3, suppress_small=True))

order = 'XYZ' # Different results for any order in ['XYZ','XZY','YZX','YXZ','ZXY','ZYX','xyz','xzy','yzx','yxz','zxy','zyx']
scipy_matrix = Rotation.from_euler(order, euler_angles).as_matrix()
print(np.array_str(scipy_matrix, precision=3, suppress_small=True))

I get two different results:

[[ 0.976 -0.083  0.2  ]
 [ 0.139  0.947 -0.288]
 [-0.165  0.309  0.937]]
[[ 0.981 -0.086  0.174]
 [ 0.136  0.943 -0.304]
 [-0.138  0.322  0.937]]

Why? How can I compute the same matrix as SimpleITK using scipy?

Upvotes: 1

Views: 894

Answers (2)

Spenhouet
Spenhouet

Reputation: 7159

The issue is that the itk.Euler3DTransform class by default applies the rotation matrix multiplications in Z @ X @ Y order and the Rotation.from_euler function in Z @ Y @ X order.

Note that this is independent of the order you specified. The order you specify refers to the order of the angles, not the order of the matrix multiplications.

If you are using the itk.Euler3DTransform directly as you showed in your example, you can actually change the default behavior for itk to perform the matrix multiplication in Z @ Y @ X order.

I never worked with sitk but in theory and based on the documentation, something like this should work:

euler_transform = sitk.Euler3DTransform((0, 0, 0), *euler_angles)
euler_transform.SetComputeZYX(True)
sitk_matrix = euler_transform.GetMatrix()

Alternatively, I wrote a function which is similar to Rotation.from_euler but has the option to specify the rotation order as well:

from typing import List
from numpy.typing import NDArray

def build_rotation_3d(radians: NDArray,
                      radians_oder: str = 'XYZ',
                      rotation_order: str = 'ZYX',
                      dims: List[str] = ['X', 'Y', 'Z']) -> NDArray:
    x_rad, y_rad, z_rad = radians[(np.searchsorted(dims, list(radians_oder)))]
    x_cos, y_cos, z_cos = np.cos([x_rad, y_rad, z_rad], dtype=np.float64)
    x_sin, y_sin, z_sin = np.sin([x_rad, y_rad, z_rad], dtype=np.float64)
    x_rot = np.asarray([
        [1.0, 0.0, 0.0, 0.0],
        [0.0, x_cos, -x_sin, 0.0],
        [0.0, x_sin, x_cos, 0.0],
        [0.0, 0.0, 0.0, 1.0],
    ])
    y_rot = np.asarray([
        [y_cos, 0.0, y_sin, 0.0],
        [0.0, 1.0, 0.0, 0.0],
        [-y_sin, 0.0, y_cos, 0.0],
        [0.0, 0.0, 0.0, 1.0],
    ])
    z_rot = np.asarray([
        [z_cos, -z_sin, 0.0, 0.0],
        [z_sin, z_cos, 0.0, 0.0],
        [0.0, 0.0, 1.0, 0.0],
        [0.0, 0.0, 0.0, 1.0],
    ])

    rotations = np.asarray([x_rot, y_rot, z_rot])[(np.searchsorted(dims, list(rotation_order)))]

    return rotations[0] @ rotations[1] @ rotations[2]

To get the respective rotation matrix which matches the EulerTransform:

rotation_matrix = build_rotation_3d(numpy.array(euler_angles),
                                    radians_oder = 'XYZ',
                                    rotation_order = 'ZXY')

Upvotes: 1

Dave Chen
Dave Chen

Reputation: 2085

What is your 'order' string. When I ran your code with order='xyz', I get the same results for SimpleITK and scipy's Rotation.

Upvotes: 0

Related Questions