SimpleITK Вращение изображения МРТ

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

Основываясь на информации, предоставленной здесь, у меня есть подозрение, что происходит. Я полагаю, что ваше МРТ имеет косинусную матрицу не в одном направлении. Вы можете подтвердить это с помощью:

print(img.GetDirection())

Выходные данные расположены в главном порядке. Когда вы делаете:

img_arr = sitk.GetArrayFromImage(img)[0]

Вы предполагаете, что матрица направляющих косинусов является тождественной. Таким образом, когда вы берете срез, перпендикулярный третьей оси, он перпендикулярен оси z, которой он не является (это может быть близко).

Чтобы вращаться вокруг оси, которая перпендикулярна плоскости осевого изображения, вам нужно взять в качестве оси вращения третий столбец матрицы косинуса направления, и вы знаете свой угол, вместе они определяют матрицу вращения ( см. Здесь для деталей ).

Затем вы можете сделать:

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

Надеюсь, это поможет.

Для будущих обсуждений, пожалуйста, придерживайтесь дискурса ITK, где вы начали первоначальное обсуждение .

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