SimpleITK Вращение объемных данных (например, МРТ)

У меня есть изображение размером 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)

Исходное изображение нулевого среза простаты

Неправильно повернутое изображение нулевого среза простаты


python image simpleitk itk
person DrJessop    schedule 16.05.2019    source источник

Ответы (2)


arrow_upward
4
arrow_downward

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

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, с которого вы начали исходное обсуждение.

person zivy    schedule 16.05.2019

arrow_upward
3
arrow_downward

Спасибо zivy и

person DrJessop    schedule 16.05.2019
comment
Интересно, могли бы вы использовать кватернионы так же, как описано в этой теме: itk-users.7.n7.nabble.com/. - person Miguel Rueda; 14.08.2020