如何从 ITK 注册中获得转换仿射？
How to get transformation affine from ITK registration?
鉴于 3D MRI 扫描 A、B 和 C，我想对 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 在 LPS 方向保存图像，在 RAS 方向保存 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. ]])
你可以看到数值很接近，只是符号不对。
事实上，如果我手动设置与“真实”矩阵中相同的符号，转换矩阵就很好。
在 MONAI 的 ITK 加载器中，我发现代码建议执行以下操作以将 ITK 仿射转换为 nibabel 仿射：
np.diag([-1, -1, 1, 1]) @ registration_affine
如果我使用 nibabels ornt_transform 方法获得从 LPS 到 RAS 的 ornt 转换，这个 returns [-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')
输入数据：
输出数据：
包版本：
itk-elastix==0.12.0
monai==0.8.0
nibabel==3.1.1
numpy==1.19.2
我之前在 GitHub #145 的 ITKElastix 项目上问过这个问题，但无法解决我的问题。感谢 dzenanz 和 mstaring 试图在那里提供帮助。
看看 this diff，您可能对旧方法更感兴趣。它直接从 4x4 矩阵构造一个 ITK 变换。
但请注意，我认为此代码中某处存在错误。我最近添加了这个，它降低了测试准确性，这让我相信其中某处存在错误。
经过大量尝试并与我的团队讨论后，我们意识到发生了什么。
我们已经确定了如何读取 ITK TransformParameters，其中前 9 个数字是按 row-major 顺序读取的旋转矩阵的一部分，而后三部分是平移矩阵。
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。
所以前两个轴的翻转。特此翻转通过 -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 的变换。
问题是 registration 仿射希望应用于 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 - 给我们一个执行 ITK 的仿射变换在面向 RAS 的图像上注册。
根据上面的最小代码示例，下面的代码现在应该可以工作了：
affine_transform = Affine(affine=registration_ras, image_only=True)
out_img = affine_transform(moving_image_np[np.newaxis, ...])
鉴于 3D MRI 扫描 A、B 和 C，我想对 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 在 LPS 方向保存图像，在 RAS 方向保存 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. ]])
你可以看到数值很接近，只是符号不对。 事实上，如果我手动设置与“真实”矩阵中相同的符号，转换矩阵就很好。
在 MONAI 的 ITK 加载器中，我发现代码建议执行以下操作以将 ITK 仿射转换为 nibabel 仿射：
np.diag([-1, -1, 1, 1]) @ registration_affine
如果我使用 nibabels ornt_transform 方法获得从 LPS 到 RAS 的 ornt 转换，这个 returns [-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')
输入数据：
输出数据：
包版本：
itk-elastix==0.12.0
monai==0.8.0
nibabel==3.1.1
numpy==1.19.2
我之前在 GitHub #145 的 ITKElastix 项目上问过这个问题，但无法解决我的问题。感谢 dzenanz 和 mstaring 试图在那里提供帮助。
看看 this diff，您可能对旧方法更感兴趣。它直接从 4x4 矩阵构造一个 ITK 变换。
但请注意，我认为此代码中某处存在错误。我最近添加了这个，它降低了测试准确性，这让我相信其中某处存在错误。
经过大量尝试并与我的团队讨论后，我们意识到发生了什么。
我们已经确定了如何读取 ITK TransformParameters，其中前 9 个数字是按 row-major 顺序读取的旋转矩阵的一部分，而后三部分是平移矩阵。
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。
所以前两个轴的翻转。特此翻转通过 -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 的变换。
问题是 registration 仿射希望应用于 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 - 给我们一个执行 ITK 的仿射变换在面向 RAS 的图像上注册。
根据上面的最小代码示例，下面的代码现在应该可以工作了：
affine_transform = Affine(affine=registration_ras, image_only=True)
out_img = affine_transform(moving_image_np[np.newaxis, ...])