我们正在使用ITK的注册算法，但我们只想要仿射变换矩阵，而不是直接应用注册。 在之前的问题中，我们已经解决了关于图像/变换方向的误解：如何从ITK注册中获得变换仿射矩阵？
现在，我们遇到了一个样本，当前的解决方案无法正常工作。旋转很好，但结果略微有偏移。 ITK的图像输出完美，因此我们知道注册成功了。这就是为什么我们将问题描述缩小到具体矩阵的仿射计算。
从ITK注册中，我们获取/读取以下参数：
parameter_map = result_transform_parameters.GetParameterMap(0)
rot00, rot01, rot02, rot10, rot11, rot12, rot20, rot21, rot22 = parameter_map[
'TransformParameters'][:9]
A = np.array([
[rot00, rot01, rot02, 0],
[rot10, rot11, rot12, 0],
[rot20, rot21, rot22, 0],
[ 0, 0, 0, 1],
], dtype=float) # yapf: disable
tx, ty, tz = parameter_map['TransformParameters'][9:]
t = np.array([
[1, 0, 0, tx],
[0, 1, 0, ty],
[0, 0, 1, tz],
[0, 0, 0, 1],
], dtype=float) # yapf: disable
# In world coordinates
cx, cy, cz = parameter_map['CenterOfRotationPoint']
c = np.array([
[1, 0, 0, cx],
[0, 1, 0, cy],
[0, 0, 1, cz],
[0, 0, 0, 1],
], dtype=float) # yapf: disable
ox, oy, oz = parameter_map['Origin']
o = np.array([
[1, 0, 0, ox],
[0, 1, 0, oy],
[0, 0, 1, oz],
[0, 0, 0, 1],
], dtype=float) # yapf: disable
moving_ras = moving_image.affine
其中A是方向/旋转矩阵，t是平移矩阵，c是旋转中心，而moving_ras是RAS方向中移动图像的仿射变换。
平移和方向矩阵可以合并为一个转换矩阵：
transform = t @ A
我们不确定如何将CenterOfRotationPoint因素考虑进去。
根据this，this和that的交流问题，我认为可能需要这样做：
transform = c @ transform @ np.linalg.inv(c)
最后，我们需要在RAS和LPS之间添加方向翻转：
registration = FLIPXY_44 @ transform @ FLIPXY_44
但这并没有导致正确的仿射变换。
在 ITK 文档和 GitHub 问题中，我们得到了以下公式来将上述参数应用于点：
T(x) = A ( x - c ) + (t + c)
虽然我们不能直接使用它，因为我们不想直接转换图像，而是只想计算正确的仿射变换矩阵，但可以看出公式与我们上面已经解释的内容非常相似。
我们的知识又到了死胡同。
我们注意到可能会出现问题的事情：
- 方向
- ITK对于图像和变换使用LPS方向
- Monai/Nibabel对于图像和变换使用RAS方向
- 旋转中心
- ITK提供所使用的旋转中心
- Monai隐含地假设旋转中心是图像的中心
- 世界空间与索引空间。
- 所有来自ITK的变换和点都在世界空间中。
- Monai似乎直接在图像上操作。
- (0, 0, 0)角 - ITK和Monai似乎使用相反的角落进行坐标 - 例如，在4x4x4图像中，ITK中的位置(0, 0, 0)是Monai中的位置(3, 3, 3)。
编辑：我注意到我的当前最小代码示例并不十分全面。因此，这里是一个更新。包含的仿射矩阵来自于ITK配准。为了简洁起见，省略了ITK代码。
这里提供新的测试数据（您可以通过MRIcoGL查看这些图像）：
这里是一个最简代码示例：
from pathlib import Path
import nibabel
import numpy as np
from monai.transforms.spatial.array import Affine
from monai.utils.enums import GridSampleMode, GridSamplePadMode
from nibabel import Nifti1Image
np.set_printoptions(suppress=True) # type: ignore
folder = Path('.')
FLIPXY_44 = np.diag([-1, -1, 1, 1])
# rot00, rot01, rot02, rot10, rot11, rot12, rot20, rot21, rot22 = parameter_map['TransformParameters'][:9]
A = np.array([[ 1.02380734, -0.05137566, -0.00766465, 0. ],
[ 0.01916231, 0.93276486, -0.23453097, 0. ],
[ 0.01808809, 0.2667324 , 0.94271694, 0. ],
[ 0. , 0. , 0. , 1. ]]) # yapf: disable
# tx, ty, tz = parameter_map['TransformParameters'][9:]
t = np.array([[ 1. , 0. , 0. , 1.12915465 ],
[ 0. , 1. , 0. , 11.76880151 ],
[ 0. , 0. , 1. , 41.54685788 ],
[ 0. , 0. , 0. , 1. ]]) # yapf: disable
# cx, cy, cz = parameter_map['CenterOfRotationPoint']
c = np.array([[ 1. , 0. , 0. , -0.1015625 ],
[ 0. , 1. , 0. , -24.5521698 ],
[ 0. , 0. , 1. , 0.1015625 ],
[ 0. , 0. , 0. , 1. ]]) # yapf: disable
# Moving image affine
x = np.array([[ 2. , 0. , 0. , -125.75732422],
[ 0. , 2. , 0. , -125.23828888],
[ 0. , 0. , 2. , -99.86506653],
[ 0. , 0. , 0. , 1. ]]) # yapf: disable
o = np.array([
[1., 0., 0., 126.8984375],
[0., 1., 0., 102.4478302],
[0., 0., 1., -126.8984375],
[0., 0., 0., 1.],
])
moving_ras = x
# Combine the direction and translation
transform = t @ A
# Factor in the center of rotation
# transform = c @ transform @ np.linalg.inv(c)
# Switch from LPS to RAS orientation
registration = FLIPXY_44 @ transform @ FLIPXY_44
y = np.array([[ 2. , 0. , 0. , -126.8984375 ],
[ 0. , 2. , 0. , -102.4478302 ],
[ 0. , 0. , 2. , -126.8984375 ],
[ 0. , 0. , 0. , 1. ]]) # yapf: disable
fixed_image_affine = y
moving_image_ni: Nifti1Image = nibabel.load(folder / 'real_moving.nii.gz')
moving_image_np: np.ndarray = moving_image_ni.get_fdata() # type: ignore
affine_transform = Affine(affine=registration,
image_only=True,
mode=GridSampleMode.NEAREST,
padding_mode=GridSamplePadMode.BORDER)
reg_monai = np.squeeze(affine_transform(moving_image_np[np.newaxis, ...]))
out = Nifti1Image(reg_monai, fixed_image_affine)
nibabel.save(out, folder / 'reg_monai.nii.gz')
当您执行此代码时，生成的reg_monai.nii.gz应该与real_fixed.nii.gz匹配（在位置和轮廓上 - 而不是实际内容）。
目前的结果如下所示（通过MRIcoGL查看）：
但是结果应该看起来像这样（这是直接从ITK注册输出中获取的硬编码仿射矩阵 - 这应该证明了注册已经成功，并且参数通常应该很好）：
为了完整起见，这里也提供执行ITK注册并获取上述仿射矩阵的代码：
from pathlib import Path
import itk
import numpy as np
np.set_printoptions(suppress=True) # type: ignore
folder = Path('.')
moving_image = itk.imread(str(folder / 'real_moving.nii.gz'), itk.F)
fixed_image = itk.imread(str(folder / 'real_fixed.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']
affine_parameter_map['MaximumNumberOfIterations'] = ['512']
parameter_object.AddParameterMap(affine_parameter_map)
# Call registration function
result_image, result_transform_parameters = itk.elastix_registration_method( # type: ignore
fixed_image, moving_image, parameter_object=parameter_object)
itk.imwrite(result_image, str(folder / 'real_reg_itk.nii.gz'), compression=True)
parameter_map = result_transform_parameters.GetParameterMap(0)
rot00, rot01, rot02, rot10, rot11, rot12, rot20, rot21, rot22 = parameter_map['TransformParameters'][:9]
A = np.array([
[rot00, rot01, rot02, 0],
[rot10, rot11, rot12, 0],
[rot20, rot21, rot22, 0],
[ 0, 0, 0, 1],
], dtype=float) # yapf: disable
tx, ty, tz = parameter_map['TransformParameters'][9:]
t = np.array([
[1, 0, 0, tx],
[0, 1, 0, ty],
[0, 0, 1, tz],
[0, 0, 0, 1],
], dtype=float) # yapf: disable
# In world coordinates
cx, cy, cz = parameter_map['CenterOfRotationPoint']
c = np.array([
[1, 0, 0, cx],
[0, 1, 0, cy],
[0, 0, 1, cz],
[0, 0, 0, 1],
], dtype=float) # yapf: disable
ox, oy, oz = parameter_map['Origin']
o = np.array([
[1, 0, 0, ox],
[0, 1, 0, oy],
[0, 0, 1, oz],
[0, 0, 0, 1],
], dtype=float) # yapf: disable
软件包版本：
itk-elastix==0.12.0
monai==0.8.0
nibabel==3.1.1
numpy==1.19.2
inv(c) @ t @ A @ c而不是c @ t @ A @ inv(c)，这将得到一个方向相同但偏移了[-2.51637511, -3.2500057, 13.11302812]的图像。 - BobA、t、c）描述了一种将一个图像叠加到另一个图像上的变换（将“移动”图像叠加到“固定”图像上）。因此，“transform”应该以一种方式转换“移动”图像，使其在视觉上覆盖“固定”图像。这应该在空间范围内工作，因此在将变换应用于数据矩阵后，经过变换的“移动”图像具有与“固定”图像相同的仿射变换，以将两个图像映射到世界空间中。 - Spenhouet