SimpleITK.Euler3DTransform 和 scipy.spatial.transform.Rotation.from_euler 之间的区别？

Question

使用这两个库函数：

• SimpleITK.Euler3DTransform
• scipy.spatial.transform.Rotation.from_euler

从欧拉角创建一个简单的旋转矩阵：

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))

我得到两个不同的结果：

[[ 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]]

为什么？如何使用 scipy 计算与 SimpleITK 相同的矩阵？

Answer 1

问题是

itk.Euler3DTransform

类默认情况下按

Z @ X @ Y

顺序应用旋转矩阵乘法，按

Rotation.from_euler

顺序应用

Z @ Y @ X

函数。

请注意，这与您指定的

order

无关。您指定的顺序是指角度的顺序，而不是矩阵乘法的顺序。

如果您像示例中所示那样直接使用

itk.Euler3DTransform

，您实际上可以更改 itk 的默认行为以按

Z @ Y @ X

顺序执行矩阵乘法。

我从来没有和 sitk 一起工作过，但在理论上和根据文档，这样的事情应该有效：

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

或者，我写了一个类似于

Rotation.from_euler

但也可以选择指定旋转顺序的函数：

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]

要获得与 EulerTransform 匹配的相应旋转矩阵：

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

Answer 2

您的“订单”字符串是什么。当我用 order='xyz' 运行你的代码时，我得到了 SimpleITK 和 scipy 的旋转相同的结果。