我想使用欧拉角 x、y 和 z 旋转我的体积 CT 数据。为此，我使用 SimpleITK。我已经阅读了Jessop 博士的问题 - simpleitk-rotation-of-volumetric-data-eg-mri并且我认为我有同样的问题，即我的方向/方向不是单位矩阵。方向即为：
0.08716564279125966, 0.0, -0.9961938319005929, 0.9961938319005927, 6.633444000000004e-17, 0.08716564279125968, 0.0, -1.0, 6.12303124808918e-17
然而，Jessop 博士找到的解决方案是使用轴角方向，以便他只能围绕 z 轴旋转。我想使用欧拉角围绕所有轴旋转。我怎样才能做到这一点？
PS 我会评论 Dr.Jessops 的问题来问它，但我没有足够的声望点。
博士的代码。杰索普：
# This function is from https://github.com/rock-learning/pytransform3d/blob/7589e083a50597a75b12d745ebacaa7cc056cfbd/pytransform3d/rotations.py#L302
def matrix_from_axis_angle(a):
""" Compute rotation matrix from axis-angle.
This is called exponential map or Rodrigues' formula.
Parameters
----------
a : array-like, shape (4,)
Axis of rotation and rotation angle: (x, y, z, angle)
Returns
-------
R : array-like, shape (3, 3)
Rotation matrix
"""
ux, uy, uz, theta = a
c = np.cos(theta)
s = np.sin(theta)
ci = 1.0 - c
R = np.array([[ci * ux * ux + c,
ci * ux * uy - uz * s,
ci * ux * uz + uy * s],
[ci * uy * ux + uz * s,
ci * uy * uy + c,
ci * uy * uz - ux * s],
[ci * uz * ux - uy * s,
ci * uz * uy + ux * s,
ci * uz * uz + c],
])
# This is equivalent to
# R = (np.eye(3) * np.cos(theta) +
# (1.0 - np.cos(theta)) * a[:3, np.newaxis].dot(a[np.newaxis, :3]) +
# cross_product_matrix(a[:3]) * np.sin(theta))
return R
def resample(image, transform):
"""
This function resamples (updates) an image using a specified transform
:param image: The sitk image we are trying to transform
:param transform: An sitk transform (ex. resizing, rotation, etc.
:return: The transformed sitk image
"""
reference_image = image
interpolator = sitk.sitkLinear
default_value = 0
return sitk.Resample(image, reference_image, transform,
interpolator, default_value)
def get_center(img):
"""
This function returns the physical center point of a 3d sitk image
:param img: The sitk image we are trying to find the center of
:return: The physical center point of the image
"""
width, height, depth = img.GetSize()
return img.TransformIndexToPhysicalPoint((int(np.ceil(width/2)),
int(np.ceil(height/2)),
int(np.ceil(depth/2))))
def rotation3d(image, theta_z, show=False):
"""
This function rotates an image across each of the x, y, z axes by theta_x, theta_y, and
theta_z degrees
respectively
:param image: An sitk MRI image
:param theta_x: The amount of degrees the user wants the image rotated around the x axis
:param theta_y: The amount of degrees the user wants the image rotated around the y axis
:param theta_z: The amount of degrees the user wants the image rotated around the z axis
:param show: Boolean, whether or not the user wants to see the result of the rotation
:return: The rotated image
"""
theta_z = np.deg2rad(theta_z)
euler_transform = sitk.Euler3DTransform()
print(euler_transform.GetMatrix())
image_center = get_center(image)
euler_transform.SetCenter(image_center)
direction = image.GetDirection()
print(direction)
axis_angle = (direction[2], direction[5], direction[8], theta_z)
np_rot_mat = matrix_from_axis_angle(axis_angle)
euler_transform.SetMatrix(np_rot_mat.flatten().tolist())
resampled_image = resample(image, euler_transform)
if show:
slice_num = int(input("Enter the index of the slice you would like to see"))
plt.imshow(sitk.GetArrayFromImage(resampled_image)[slice_num])
plt.show()
return resampled_image
要从欧拉角方法中获取旋转矩阵，可以使用以下代码：
def matrix_from_euler_xyz(e):
"""Compute rotation matrix from xyz Euler angles.
Intrinsic rotations are used to create the transformation matrix
from three concatenated rotations.
The xyz convention is usually used in physics and chemistry.
Parameters
----------
e : array-like, shape (3,)
Angles for rotation around x-, y'-, and z''-axes (intrinsic rotations)
Returns
-------
R : array-like, shape (3, 3)
Rotation matrix
"""
alpha, beta, gamma = e
# We use intrinsic rotations
Qx = matrix_from_angle(0, alpha)
Qy = matrix_from_angle(1, beta)
Qz = matrix_from_angle(2, gamma)
R = Qx.dot(Qy).dot(Qz)
return R
但是，仍应包含方向。有谁知道如何做到这一点？