SimpleITK 3Deuler角度旋转体积数据
SimpleITK 3Deuler angles rotation volumetric data
我想使用欧拉角 x、y 和 z 旋转体积 CT 数据。为此，我使用 SimpleITK。我已经阅读了 的问题，我认为我有同样的问题，我的 orientation/direction 不是单位矩阵。方向即为：
0.08716564279125966, 0.0, -0.9961938319005929, 0.9961938319005927, 6.633444000000004e-17, 0.08716564279125968, 0.0, -1.0, 6.12303124808918e-17
然而，Jessop 博士找到的解决方案是使用轴角方向，这样他就只能绕 z 轴旋转。我想使用欧拉角绕所有轴旋转。我怎样才能做到这一点？
P.S。我会评论 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 执行，使用 VersorRigid3DTransform 并将其传递到 resampleitkfilter，我相信这就是这段代码在某种程度上所做的。
我建议使用 scipy.transform.rotation 模块将 versor 传入 rigid3Dtransform，然后相应地设置从输入到输出的方向。
我 'think' 这甚至适用于你的单位外方向，但与 3D 仿射变换一样，这可能会遗漏一些更聪明的数学家可以批评的东西：
r = R.from_euler('zyx', [90, 45, 30], degrees=True)
orientation = r.as_quat()
transform = sitk.VersorRigid3DTransform()
translation = sitk.TranslationTransform(3, (x,y,z))
transform.SetTranslation(translation.GetOffset())
#get rot_center from your center of image function as tuple (x,y,z)
rotation = sitk.VersorTransform([orientation[0], orientation[1], orientation[2], orientation[3]],rot_center)
transform.SetRotation(rotation.GetVersor())
transform.SetCenter(rotation.GetCenter())
size = img.GetSize()
bounds = list()
for x in [0, size[0]]:
for y in [0, size[1]]:
for z in [0, size[2]]:
bounds.append(img.TransformIndexToPhysicalPoint((x, y, z)))
# get the physical position of the bounds given transform
trans_bounds = list()
for b in bounds:
trans_bounds.append(transform.TransformPoint((b[0], b[1], b[2])))
trans_bounds = np.array(trans_bounds)
xmin = np.min(trans_bounds.T[0])
xmax = np.max(trans_bounds.T[0])
ymin = np.min(trans_bounds.T[1])
ymax = np.max(trans_bounds.T[1])
zmin = np.min(trans_bounds.T[2])
zmax = np.max(trans_bounds.T[2])
output_origin = (xmin, ymin, zmin)
res_all = (0.5, 0.5, 0.5)
res = res_all[0]
output_size = np.array(
[int(np.round((xmax - xmin) / res)), int(np.round((ymax - ymin) / res)), int(np.round((zmax - zmin) / res))])
output_size = output_size.astype(int)
resampleFilter = sitk.ResampleImageFilter()
resampleFilter.SetTransform(transform.GetInverse())
resampleFilter.SetInterpolator(sitk.sitkLinear)
resampleFilter.SetSize(output_size.tolist())
resampleFilter.SetOutputOrigin(output_origin)
resampleFilter.SetOutputSpacing(res_all)
resampleFilter.SetOutputDirection(img.GetDirection())
resampleFilter.SetDefaultPixelValue(0.0)
我解决了（借助 itk discourse
的 zivy
我的回答是这样的：
# The function which is used to rotate (and make the 3D image isotropic) using euler angles
def rotation3d(image, theta_x, theta_y, theta_z, output_spacing = None, background_value=0.0):
"""
This function rotates an image across each of the x, y, z axes by theta_x, theta_y, and theta_z degrees
respectively (euler ZXY orientation) and resamples it to be isotropic.
:param image: An sitk 3D 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 output_spacing: Scalar denoting the isotropic output image spacing. If None, then use the smallest
spacing from original image.
:return: The rotated image
"""
euler_transform = sitk.Euler3DTransform (image.TransformContinuousIndexToPhysicalPoint([(sz-1)/2.0 for sz in image.GetSize()]),
np.deg2rad(theta_x),
np.deg2rad(theta_y),
np.deg2rad(theta_z))
# compute the resampling grid for the transformed image
max_indexes = [sz-1 for sz in image.GetSize()]
extreme_indexes = list(itertools.product(*(list(zip([0]*image.GetDimension(),max_indexes)))))
extreme_points_transformed = [euler_transform.TransformPoint(image.TransformContinuousIndexToPhysicalPoint(p)) for p in extreme_indexes]
output_min_coordinates = np.min(extreme_points_transformed, axis=0)
output_max_coordinates = np.max(extreme_points_transformed, axis=0)
# isotropic ouput spacing
if output_spacing is None:
output_spacing = min(image.GetSpacing())
output_spacing = [output_spacing]*image.GetDimension()
output_origin = output_min_coordinates
output_size = [int(((omx-omn)/ospc)+0.5) for ospc, omn, omx in zip(output_spacing, output_min_coordinates, output_max_coordinates)]
output_direction = [1,0,0,0,1,0,0,0,1]
output_pixeltype = image.GetPixelIDValue()
return sitk.Resample(image,
output_size,
euler_transform.GetInverse(),
sitk.sitkLinear,
output_origin,
output_spacing,
output_direction,
background_value,
output_pixeltype)
我想使用欧拉角 x、y 和 z 旋转体积 CT 数据。为此，我使用 SimpleITK。我已经阅读了
0.08716564279125966, 0.0, -0.9961938319005929, 0.9961938319005927, 6.633444000000004e-17, 0.08716564279125968, 0.0, -1.0, 6.12303124808918e-17
然而，Jessop 博士找到的解决方案是使用轴角方向，这样他就只能绕 z 轴旋转。我想使用欧拉角绕所有轴旋转。我怎样才能做到这一点？
P.S。我会评论 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 执行，使用 VersorRigid3DTransform 并将其传递到 resampleitkfilter，我相信这就是这段代码在某种程度上所做的。
我建议使用 scipy.transform.rotation 模块将 versor 传入 rigid3Dtransform，然后相应地设置从输入到输出的方向。
我 'think' 这甚至适用于你的单位外方向，但与 3D 仿射变换一样，这可能会遗漏一些更聪明的数学家可以批评的东西：
r = R.from_euler('zyx', [90, 45, 30], degrees=True)
orientation = r.as_quat()
transform = sitk.VersorRigid3DTransform()
translation = sitk.TranslationTransform(3, (x,y,z))
transform.SetTranslation(translation.GetOffset())
#get rot_center from your center of image function as tuple (x,y,z)
rotation = sitk.VersorTransform([orientation[0], orientation[1], orientation[2], orientation[3]],rot_center)
transform.SetRotation(rotation.GetVersor())
transform.SetCenter(rotation.GetCenter())
size = img.GetSize()
bounds = list()
for x in [0, size[0]]:
for y in [0, size[1]]:
for z in [0, size[2]]:
bounds.append(img.TransformIndexToPhysicalPoint((x, y, z)))
# get the physical position of the bounds given transform
trans_bounds = list()
for b in bounds:
trans_bounds.append(transform.TransformPoint((b[0], b[1], b[2])))
trans_bounds = np.array(trans_bounds)
xmin = np.min(trans_bounds.T[0])
xmax = np.max(trans_bounds.T[0])
ymin = np.min(trans_bounds.T[1])
ymax = np.max(trans_bounds.T[1])
zmin = np.min(trans_bounds.T[2])
zmax = np.max(trans_bounds.T[2])
output_origin = (xmin, ymin, zmin)
res_all = (0.5, 0.5, 0.5)
res = res_all[0]
output_size = np.array(
[int(np.round((xmax - xmin) / res)), int(np.round((ymax - ymin) / res)), int(np.round((zmax - zmin) / res))])
output_size = output_size.astype(int)
resampleFilter = sitk.ResampleImageFilter()
resampleFilter.SetTransform(transform.GetInverse())
resampleFilter.SetInterpolator(sitk.sitkLinear)
resampleFilter.SetSize(output_size.tolist())
resampleFilter.SetOutputOrigin(output_origin)
resampleFilter.SetOutputSpacing(res_all)
resampleFilter.SetOutputDirection(img.GetDirection())
resampleFilter.SetDefaultPixelValue(0.0)
我解决了（借助 itk discourse
的 zivy我的回答是这样的：
# The function which is used to rotate (and make the 3D image isotropic) using euler angles
def rotation3d(image, theta_x, theta_y, theta_z, output_spacing = None, background_value=0.0):
"""
This function rotates an image across each of the x, y, z axes by theta_x, theta_y, and theta_z degrees
respectively (euler ZXY orientation) and resamples it to be isotropic.
:param image: An sitk 3D 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 output_spacing: Scalar denoting the isotropic output image spacing. If None, then use the smallest
spacing from original image.
:return: The rotated image
"""
euler_transform = sitk.Euler3DTransform (image.TransformContinuousIndexToPhysicalPoint([(sz-1)/2.0 for sz in image.GetSize()]),
np.deg2rad(theta_x),
np.deg2rad(theta_y),
np.deg2rad(theta_z))
# compute the resampling grid for the transformed image
max_indexes = [sz-1 for sz in image.GetSize()]
extreme_indexes = list(itertools.product(*(list(zip([0]*image.GetDimension(),max_indexes)))))
extreme_points_transformed = [euler_transform.TransformPoint(image.TransformContinuousIndexToPhysicalPoint(p)) for p in extreme_indexes]
output_min_coordinates = np.min(extreme_points_transformed, axis=0)
output_max_coordinates = np.max(extreme_points_transformed, axis=0)
# isotropic ouput spacing
if output_spacing is None:
output_spacing = min(image.GetSpacing())
output_spacing = [output_spacing]*image.GetDimension()
output_origin = output_min_coordinates
output_size = [int(((omx-omn)/ospc)+0.5) for ospc, omn, omx in zip(output_spacing, output_min_coordinates, output_max_coordinates)]
output_direction = [1,0,0,0,1,0,0,0,1]
output_pixeltype = image.GetPixelIDValue()
return sitk.Resample(image,
output_size,
euler_transform.GetInverse(),
sitk.sitkLinear,
output_origin,
output_spacing,
output_direction,
background_value,
output_pixeltype)