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

Question

我有一个 32x32x3（高、宽、深）图像，我试图在 sitk 中围绕 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 轴垂直，而事实并非如此（它可能很接近）。

要绕与轴向图像平面垂直的轴旋转，您需要将方向余弦矩阵的第三列作为您的旋转轴，并且您知道您的角度，它们一起定义了一个旋转矩阵（有关详细信息，请参见此处）。

然后你可以这样做：

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

希望这可以帮助。

对于未来的讨论，请坚持使用您开始原始讨论的 ITK 话语。