不同nifti成像数据间IJK坐标的对齐

不同的神经影像学数据有不同的IJK坐标，这里以Python库nibabel为例，简单介绍了离散数据（ROI）的坐标变换的基本方法。

不同成像数据的直观对比

经过预处理的数据神经影像学数据，一般可以视为三维或者四维的数组，前三维相当于空间上的三维，而第四维可能是功能成像的时间或弥散张量成像的不同磁感应强度。当然，这里主要考虑以体素/小立方体为基本单元的nifti数据。

>>> import nibabel
>>> mdtb = nibabel.load('mdtb_mni.nii')
>>> mdtb.shape
(141, 95, 87)
>>> hcp = nibabel.load('S1200_AverageT1w_restore.nii.gz')
>>> hcp.shape
(260, 311, 260)

这里举出了两个例子，一个是基于mdtb的小脑ROI文件，一个是来自HCP的被试平均结构像问及爱你。mdtb_mni.nii主要储存了一个形状为(141, 95, 87)的数组，每个索引对应的取值是该体素的ROI编号，而S1200_AverageT1w_restore.nii.gz主要储存了一个形状为(260, 311, 260)的数组，每个索引对应的取值是该体素在t1加权成像中的数值大小。在下文中，我们用IJK坐标指代数组的“坐标”。

在实际应用中，一个现实的问题是，我们可能希望根据一个ROI文件提供的信息，在另一个实际成像数据中找到对应的体素。然而，如同上面的示例文件，他们可能具有不同的IJK坐标系，两者自然也无法直接对齐。而既然他们可以对应一定空间中的脑结构，那为什么会有这种差异？可以找到两个直接的原因：

它们有不同的空间分辨率。直观来看，分辨率更高的数据，对于相同体积的空间，需要更多的数据来表征，也就是有更“大”的数组。如果对数据进行降采样，数组也就会相应变小。
它们描述了大小不同的空间。mdtb提供的是小脑的ROI分区，所以本身并不需要存储在小脑范围以外的数据，也就是具有更小的FOV，而S1200_AverageT1w_restore.nii.gz是全脑的结构，自然描述了更大的空间范围。

桥梁：Voxel Space与仿射变换

尽管如此，我们把这两个文件放进可视化工具wb_view里的时候会发现，可视化工具能够轻松地将两个文件各自的数据在空间中对应起来，如下图所示。

（图）

随意选中一个点，我们可以从workbench提供的信息中看出端倪：两个数组的IJK坐标不同，但如果两个数据中存在空间上对应的索引，那么这两个索引会对应Voxel Space中相同的坐标。也就是说，Voxel Space是他们的共同语言。那么，IJK坐标又是怎么和Voxel Space对应起来的呢？答案就在nifti数据头文件的仿射矩阵里面，nibabel的使用教程也提供了相应的介绍。仿射矩阵描述了从IJK到Voxel Space的线性变换：

$P_{Voxel} = AP_{IJK}$

在nibabel中，也很容易查看不同成像数据的仿射矩阵：

>>> mdtb.affine
array([[   1.,    0.,    0.,  -70.],
       [   0.,    1.,    0., -100.],
       [   0.,    0.,    1.,  -75.],
       [   0.,    0.,    0.,    1.]])
>>> hcp.affine
array([[  -0.69999999,    0.        ,    0.        ,   90.        ],
       [   0.        ,    0.69999999,    0.        , -126.        ],
       [   0.        ,    0.        ,    0.69999999,  -72.        ],
       [   0.        ,    0.        ,    0.        ,    1.        ]])

其中A一般是一个四阶方阵，对角线上的前三个值表示了各个轴上的整体伸缩，第四列则用于表示坐标的平移。IJK坐标是三维的，如果要直接左乘仿射矩阵，需要补充一个第四维的1。在实际计算中，我们一般单独取仿射矩阵的前三行、前三列与IJK坐标相乘，最后再加上第四列前三行的偏移量。另外，在相同的Voxel Space里面，前三行/列的交叉项一般是0，也就是说不同数据中的IJK空间只需要通过伸缩、平移就可以互相转换。S1200_AverageT1w_restore.nii.gz的伸缩项要小一些，相当于用更多的点表示相同体积的空间，所以空间分辨率也就高一些。

对于mdtb的仿射矩阵A，可以验证上图中我们选取的点从IJK坐标到Voxel Space的变换：

$\begin{aligned} P_{Voxel} = A_{(i, j) \in \{1, 2, 3\} \times \{1, 2, 3\} }P_{IJK} + A_{Offset} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ \end{bmatrix} \times \begin{bmatrix} ? \\ ? \\ ? \\ \end{bmatrix} + \begin{bmatrix} -70 \\ -100 \\ -75 \\ \end{bmatrix} = \begin{bmatrix} ? \\ ? \\ ? \\ \end{bmatrix} \end{aligned}$

这就得到了上图中Voxel Space的坐标。

通过仿射矩阵对齐IJK坐标

明白了两个成像数据间的桥梁之后，就可以回到最初的问题了：如何使记录体素ROI编号数据的IJK坐标与结构像的IJK坐标对齐，从而确认结构像中的某一个点所对应的ROI编号？

参考资料

nibabel

hcp

mdtb