如何从ITK注册中获得变换仿射？

进行 3D MRI 扫描A, B, and C我想执行仿射（联合）配准B onto A，取配准的变换仿射矩阵并将其应用于C.

我的问题是配准变换的仿射矩阵的符号错误。也许是因为方向错误？

The 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 将图像保存在LPS方向和尼巴贝尔RAS。 所以我尝试对方向差异进行更改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.        ]])

您可以看到这些值非常接近，但只有符号错误。 事实上，如果我手动设置与“真实”矩阵中相同的符号，则变换矩阵是好的。

在 MONAI 的 ITK 加载器中，我发现建议执行以下操作以将 ITK 仿射转换为 nibabel 仿射的代码：

np.diag([-1, -1, 1, 1]) @ registration_affine

如果我使用 nibabelsornt_transform从中获取 ornt 变换的方法LPS to RAS，这返回[-1, -1, 1]并与 MONAI 的 ITK 加载器中所做的匹配。

但是将其应用于上面的仿射实际上并不会产生正确的符号（仅在平移位中）：

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')

输入数据：

• 固定_2mm.nii.gz https://github.com/InsightSoftwareConsortium/ITKElastix/files/8226904/fixed_2mm.nii.gz
• moving_2mm.nii.gz https://github.com/InsightSoftwareConsortium/ITKElastix/files/8226905/moving_2mm.nii.gz

输出数据：

• reg_itk.nii.gz https://github.com/InsightSoftwareConsortium/ITKElastix/files/8226906/reg_itk.nii.gz
• reg_monai.nii.gz https://github.com/InsightSoftwareConsortium/ITKElastix/files/8226907/reg_monai.nii.gz

封装版本：

itk-elastix==0.12.0
monai==0.8.0
nibabel==3.1.1
numpy==1.19.2

我之前确实在 GitHub 上的 ITKElastix 项目上问过这个问题#145 https://github.com/InsightSoftwareConsortium/ITKElastix/issues/145但无法解决我的问题。感谢 dzenanz 和 mstaring 试图在那里提供帮助。

经过与我的团队进行大量尝试和讨论后，我们意识到正在发生的事情。

我们已经确定了如何阅读 ITKTransformParameters前 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 方向的图像和 LPS 方向的 ITK 图像。

我们也已经知道，如果我们想要改变图像的方向，我们需要翻转相应的轴。 LPS 到 RAS 意味着翻转 L->R 和 P->A。 所以翻转前两个轴。 Flip由此表示为-1并且没有翻页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]])

请注意，这种翻转是双向的。 LPS -> RAS 和 RAS -> LPS。

如果您有 3D 图像和相应的仿射矩阵，则可以通过应用“flip_LPS_RAS”来翻转该图像中的轴。 如果你想计算该图像的新仿射，你可以这样做：

flipped = flip_LPS_RAS @ image_affine

既然我们已经奠定了基础，那么现在让我们看看我们未能弄清楚的内容。

我们知道配准变换的仿射矩阵基于 LPS 方向的图像，并且我们知道 nibabel 图像位于 RAS。 思路是，我们需要将仿射变换从 LPS 方向转换为 RAS 方向，类似于上面提到的图像重新定向。 所以我们应用了flip_LPS_RAS仿射于registration仿射。 我们的错误在于，这并没有使仿射变换成为面向 RAS 的变换。

问题是registrationaffine 期望应用于 LPS 方向的图像并输出 LPS 方向的图像。 让我们回顾一下，图像仿射具有 RAS 方向，并且registration仿射期望应用于 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 - 为我们提供一个仿射变换，在面向 RAS 上执行 ITK 配准图像。

从上面的最小代码示例来看，以下内容现在应该可以工作：

affine_transform = Affine(affine=registration_ras, image_only=True)
out_img = affine_transform(moving_image_np[np.newaxis, ...])