体积数据的 SimpleITK 旋转（例如 MRI）

Question

我有一个 32x32x3（高度、宽度、深度）的图像，我试图在坐姿中围绕 z 轴旋转 45 度。然而，看起来我正在旋转的 z/深度轴是一个角度。我如何能够旋转图像，这样如果我看图像的一个切片，我会看到切片从中心旋转了 45 度？

下面是我的代码，下面是图像（第一个是原始图像，第二个是旋转尝试失败）。 此外，这些是公共图像，而不是机密数据。

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.sitkBSpline
  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_x, theta_y, 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_x = np.deg2rad(theta_x)
  theta_y = np.deg2rad(theta_y)
  theta_z = np.deg2rad(theta_z)
  euler_transform = sitk.Euler3DTransform(get_center(image), theta_x, theta_y, theta_z, (0, 0, 0))
  image_center = get_center(image)
  euler_transform.SetCenter(image_center)
  euler_transform.SetRotation(theta_x, theta_y, theta_z)
  resampled_image = resample(image, euler_transform)
  if show:
    plt.imshow(sitk.GetArrayFromImage(resampled_image)[0])
    plt.show()
  return resampled_image

if __name__ == "__main__":
  img = sitk.ReadImage("...")
  img_arr = sitk.GetArrayFromImage(img)[0] # Represents the 0th slice, since numpy swaps the first and third axes default to sitk
  plt.imshow(img_arr); plt.show()
  input("Press enter to continue...")
  rotation3d(img, 0, 0, 45, show=True)

Answer 1

根据此处提供的信息，我怀疑发生了什么。我相信您的 MRI 扫描具有非单位方向余弦矩阵。您可以通过以下方式确认：

print(img.GetDirection())

输出按行主要顺序。 当你这样做时：

img_arr = sitk.GetArrayFromImage(img)[0]

您假设方向余弦矩阵是恒等式。因此，当您抓取垂直于第三轴的切片时，它垂直于 z 轴，但并非如此（它可能很接近）。

要绕垂直于轴向图像平面的轴旋转，您需要将方向余弦矩阵的第三列作为旋转轴，并且您知道角度，它们一起定义了一个旋转矩阵 (see here for details)。

你可以这样做：

np_rot_mat = compute_rotation_matrix_from_axis_angle()
euler_transform.SetMatrix(np_rot_mat.flatten().tolist())

希望这会有所帮助。

以后的讨论，请坚持ITK discourse where you started the original discussion。

Answer 2

感谢 zivy 和 https://github.com/rock-learning/pytransform3d/blob/7589e083a50597a75b12d745ebacaa7cc056cfbd/pytransform3d/rotations.py#L302，我现在有了解决问题的方法。下面的代码可以正常工作：

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