【发布时间】：2021-09-08 15:20:53
【问题描述】：
我想使用欧拉角 x、y 和 z 旋转我的体积 CT 数据。为此，我使用 SimpleITK。我已经阅读了来自Dr. Jessop - simpleitk-rotation-of-volumetric-data-e-g-mri 的问题，并且我认为我有同样的问题，即我的方向/方向不是单位矩阵。方向即为：
0.08716564279125966, 0.0, -0.9961938319005929, 0.9961938319005927, 6.633444000000004e-17, 0.08716564279125968, 0.0, -1.0, 6.12303124808918e-17
但是，Jessop 博士找到的解决方案是使用轴角方向，以便他只能围绕 z 轴旋转。我想使用欧拉角围绕所有轴旋转。我怎样才能做到这一点？
附：我会评论 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
但是，方向仍应包含在内。有人知道怎么做吗？
【问题讨论】：
-
这可以使用 SimpleITK 使用 Euler3DTransform 并将其传递给 resampleitkfilter 来执行，我相信这就是这段代码正在做的事情。添加方向 -> GetDirection() 将返回一个单位方向矩阵/余弦。即单位矩阵。
标签： python orientation direction euler-angles simpleitk