SimpleITK的差异.Euler3DTransform和scipy.spatial.transform. rotat

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的函数，但也有指定旋转顺序的选项:

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]

得到与EulerTransform相匹配的旋转矩阵:

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

Answer 2

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