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

我有一个 32x32x3（高度、宽度、深度）的图像，我试图在 satk 中围绕 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)

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

print(img.GetDirection())

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

img_arr = sitk.GetArrayFromImage(img)[0]

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

要绕垂直于轴向图像平面的轴旋转，您需要将方向余弦矩阵的第三列作为旋转轴，并且您知道角度，它们一起定义旋转矩阵（详情请参阅此处 https://en.wikipedia.org/wiki/Rotation_matrix#Rotation_matrix_from_axis_and_angle).

然后你可以这样做：

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

希望这可以帮助。

对于未来的讨论，请坚持ITK 讨论是您开始原始讨论的地方 https://discourse.itk.org/t/rotation-around-a-specified-axis/1875/3.