如何从ITK配准中获得变换仿射？
在三维核磁共振扫描A、B和C的情况下，我想在A上执行B的仿射(Co)配准，获取注册的仿射矩阵，并将其应用于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方法将ornt从LPS转换为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')输入数据：
产出数据：
包版本：
itk-elastix==0.12.0
monai==0.8.0
nibabel==3.1.1
numpy==1.19.2我以前在ITKElastix项目GitHub https://github.com/InsightSoftwareConsortium/ITKElastix/issues/145上问过这个问题，但没能解决我的问题。多亏了dzenanz和瞪着眼睛想帮忙的人。
回答 2
Stack Overflow用户
发布于 2022-03-15 16:45:00
在与我的团队进行了大量的尝试和讨论之后，我们终于意识到了正在发生的事情。
我们已经建立了如何读取ITK 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方向有图像，在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方向的转换仿射转换成类似于上述图像重定向的RAS取向。因此，我们将flip_LPS_RAS仿射应用于registration仿射。我们所犯的错误是，这并不使仿射变换成为面向RAS的变换。
问题是，registration仿射期望应用于一个图像的脂多糖方向，并输出一个图像的脂多糖方向。让我们回顾一下，图像仿射具有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, ...])Stack Overflow用户
发布于 2022-03-11 19:38:51
看一看这种差异，您可能更感兴趣的是旧的方法。它直接从4x4矩阵构造一个ITK变换。
但是请注意，我认为在这段代码中有一个bug。我最近增加了这个，它降低了测试的准确性，这让我相信在某个地方有一个bug。
https://stackoverflow.com/questions/71441883
复制相似问题
腾讯云开发者
