Как получить аффинное преобразование от регистрации ITK?
Имея 3D МРТ-сканы , , и я хочу выполнить аффинную (со)регистрацию
Bна
A, возьмем аффинную матрицу преобразования регистрации и применим ее к
C.
Моя проблема в том, что аффинная матрица регистрационного преобразования имеет неправильные знаки. Может из-за неправильной ориентации?
В
TransformParametersсодержат 12 значений, из которых первые 9 представляют собой матрицу вращения в порядке возрастания строк, а последние 3 являются значениями перемещения.
TransformParameters = [R1, R2, R3, R4, R5, R6, R7, R8, R9, Tx, Ty, Tz]
registration_affine = [[R1, R2, R3, Tx],
[R4, R5, R6, Ty],
[R7, R8, R9, Tz],
[0, 0, 0, 1 ]]
Я знаю, что ITK хранит изображения в ориентациях и nibabel в . Поэтому я попытался применить изменение относительно разницы ориентации к
transform_affineно это не сработало.
Я не могу получить тот же результат регистрации, что и ITK, ниже я покажу несколько примеров чисел и мой минимальный пример кода.
Чтобы проверить это, я применил аффинное преобразование к существующему изображению. Обратная этой матрице преобразования является истинным аффинным, которое может найти регистрация.
array([[ 1.02800583, 0.11462834, -0.11426342, -0.43383606],
[ 0.11462834, 1.02800583, -0.11426342, 0.47954143],
[-0.11426342, -0.11426342, 1.02285268, -0.20457054],
[ 0. , 0. , 0. , 1. ]])
Но аффинное построение, как объяснялось выше, дает:
array([[ 1.02757335, 0.11459412, 0.11448339, 0.23000557],
[ 0.11410441, 1.02746452, 0.11413955, -0.20848751],
[ 0.11398788, 0.11411115, 1.02255042, -0.04884404],
[ 0. , 0. , 0. , 1. ]])
Вы можете видеть, что значения довольно близки, но неверны только знаки. На самом деле, если я вручную устанавливаю те же знаки, что и в «истинной» матрице, матрица преобразования хороша.
В загрузчике ITK MONAI я нашел код, предлагающий сделать следующее для преобразования аффинного ITK в аффинный nibabel:
np.diag([-1, -1, 1, 1]) @ registration_affine
Если я использую нибабелс
ornt_transformметоды для получения преобразования ornt из
LPSк
RAS, это возвращает
[-1, -1, 1]и соответствует тому, что делается в загрузчике ITK MONAI.
Но применение этого к аффинному сверху на самом деле не дает правильных знаков (только в бите перевода):
array([[-1.02757335, -0.11459412, -0.11448339, -0.23000557],
[-0.11410441, -1.02746452, -0.11413955, 0.20848751],
[ 0.11398788, 0.11411115, 1.02255042, -0.04884404],
[ 0. , 0. , 0. , 1. ]])
Так что я немного застрял здесь.
Вот полный минимальный пример кода для запуска того, что я делаю/пытаюсь сделать. См. также примеры данных и версии пакетов ниже.
import nibabel
import numpy as np
from monai.transforms import Affine
from nibabel import Nifti1Image
import itk
# Import Images
moving_image = itk.imread('moving_2mm.nii.gz', itk.F)
fixed_image = itk.imread('fixed_2mm.nii.gz', itk.F)
# Import Default Parameter Map
parameter_object = itk.ParameterObject.New()
affine_parameter_map = parameter_object.GetDefaultParameterMap('affine', 4)
affine_parameter_map['FinalBSplineInterpolationOrder'] = ['1']
parameter_object.AddParameterMap(affine_parameter_map)
# Call registration function
result_image, result_transform_parameters = itk.elastix_registration_method(
fixed_image, moving_image, parameter_object=parameter_object)
parameter_map = result_transform_parameters.GetParameterMap(0)
transform_parameters = np.array(parameter_map['TransformParameters'], dtype=float)
itk.imwrite(result_image, 'reg_itk.nii.gz', compression=True)
# Convert ITK params to affine matrix
rotation = transform_parameters[:9].reshape(3, 3)
translation = transform_parameters[-3:][..., np.newaxis]
reg_affine: np.ndarray = np.append(rotation, translation, axis=1) # type: ignore
reg_affine = np.append(reg_affine, [[0, 0, 0, 1]], axis=0) # type: ignore
# Apply affine transform matrix via MONAI
moving_image_ni: Nifti1Image = nibabel.load('moving_2mm.nii.gz')
fixed_image_ni: Nifti1Image = nibabel.load('fixed_2mm.nii.gz')
moving_image_np: np.ndarray = moving_image_ni.get_fdata() # type: ignore
LPS = nibabel.orientations.axcodes2ornt(('L', 'P', 'S'))
RAS = nibabel.orientations.axcodes2ornt(('R', 'A', 'S'))
ornt_transform = nibabel.orientations.ornt_transform(LPS, RAS)[:, -1] # type: ignore
affine_transform = Affine(affine=np.diag([*ornt_transform, 1]) @ reg_affine, image_only=False)
out_img, out_affine = affine_transform(moving_image_np[np.newaxis, ...])
reg_monai = np.squeeze(out_img)
out = Nifti1Image(reg_monai, fixed_image_ni.affine, header=fixed_image_ni.header)
nibabel.save(out, 'reg_monai.nii.gz')
Входные данные:
Выходные данные:
Версии пакета:
itk-elastix==0.12.0
monai==0.8.0
nibabel==3.1.1
numpy==1.19.2
Я задавал этот вопрос раньше в проекте ITKElastix на GitHub #145 , но не смог решить свои проблемы. Спасибо dzenanz и mstaring, которые пытались помочь там.
2 ответа
После долгих попыток и обсуждений с моей командой мы пришли к пониманию того, что происходит.
Мы уже установили, как читать ИТК
TransformParametersс первыми 9 числами, являющимися частью матрицы поворота, прочитанной в главном порядке строки и последних трех частей матрицы перевода.
rot00, rot01, rot02, rot10, rot11, rot12, rot20, rot21, rot22, tx, ty, tz = parameter_map['TransformParameters']
affine = np.array([
[rot00, rot01, rot02, tx],
[rot10, rot11, rot12, ty],
[rot20, rot21, rot22, tz],
[ 0, 0, 0, 1],
], dtype=np.float32) # yapf: disable
Мы уже знали, что у nibabel есть изображения в ориентации RAS и ITK в ориентации LPS.
Мы также уже знали, что если мы хотим изменить ориентацию изображения, нам нужно перевернуть соответствующую ось. LPS в RAS означает переключение L->R и P->A. Итак, переворот первых двух осей. Флип настоящим представлен через
-1и никакого перелистывания. Таким образом, переворот первых двух осей можно описать как
[-1, -1, 1]. Мы можем построить матрицу аффинного преобразования для этого флипа с
np.diag([-1, -1, 1, 1])(последний
1только для расчетов). Таким образом, матрица аффинного преобразования для переключения между LPS и RAS:
flip_LPS_RAS = np.array([[-1, 0, 0, 0],
[ 0, -1, 0, 0],
[ 0, 0, 1, 0],
[ 0, 0, 0, 1]])
Обратите внимание, что этот флип работает в обе стороны. ЛПС -> РАН и РАН -> ЛПС.
Если у вас есть трехмерное изображение и соответствующая аффинная матрица, вы можете перевернуть ось в этом изображении, применив «flip_LPS_RAS». Если вы хотите рассчитать новый аффин этого изображения, вы можете сделать:
flipped = flip_LPS_RAS @ image_affine
Поскольку у нас есть основа, давайте теперь посмотрим, что нам не удалось выяснить.
Мы знали, что аффинная матрица регистрационного преобразования основана на изображении в ориентации LPS, и мы знали, что изображение Нибабеля находится в RAS. Мысль заключалась в том, что нам нужно преобразовать аффинное преобразование из ориентации LPS в ориентацию RAS, аналогичное переориентации изображения, упомянутой выше. Поэтому мы применили
flip_LPS_RASаффинное к аффинному. Мы ошиблись в том, что это не делает аффинное преобразование преобразованием, ориентированным на RAS.
Дело в том, что affine ожидает применения к изображению в LPS-ориентации и выводит изображение в LPS-ориентации. Напомним, что аффинное изображение имеет RAS-ориентацию и что
registrationaffine ожидает применения к изображению в LPS-ориентации. Теперь становится легче увидеть, что для применения преобразования регистрации к изображению в ориентации RAS нам сначала нужно изменить ориентацию изображения на LPS, а после регистрации обратно на RAS.
image -> flip_LPS_RAS -> registration_lps -> flip_LPS_RAS
Нас интересует только аффинная матрица для регистрационного преобразования, поэтому давайте проигнорируем изображение в приведенной выше цепочке преобразований. Написание этой цепочки аффинных преобразований в коде:
registration_ras = flip_LPS_RAS @ registration_lps @ flip_LPS_RAS
Это даст аффинную матрицу, которая принимает изображение, ориентированное на RAS, изменяет его на ориентацию LPS, выполняет регистрацию в ориентации LPS, а затем изменяет ориентацию обратно на RAS, что дает нам одно аффинное преобразование, которое выполняет регистрацию ITK на ориентированном на RAS. изображение.
Из приведенного выше минимального примера кода теперь должно работать следующее:
affine_transform = Affine(affine=registration_ras, image_only=True)
out_img = affine_transform(moving_image_np[np.newaxis, ...])
Взглянув на этот diff , вы можете быть более заинтересованы в старом способе сделать это. Он напрямую строит преобразование ITK из матрицы 4x4.
Но будьте осторожны, я думаю, что где-то в этом коде есть ошибка. Я добавил это недавно, и это снизило точность теста, что заставляет меня поверить, что где-то здесь есть ошибка.