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

不同成像数据的直观对比

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

1
2
3
4
5
6
7
>>> 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坐标系,两者自然也无法直接对齐。而既然他们可以对应一定空间中的脑结构,那为什么会有这种差异?可以找到两个直接的原因:

  1. 它们有不同的空间分辨率。直观来看,分辨率更高的数据,对于相同体积的空间,需要更多的数据来表征,也就是有更“大”的数组。如果对数据进行降采样,数组也就会相应变小。
  2. 它们描述了大小不同的空间。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的线性变换:

nibabel中,也很容易查看不同成像数据的仿射矩阵:

1
2
3
4
5
6
7
8
9
10
>>> 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的变换:

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

通过仿射矩阵对齐IJK坐标

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

参考资料

nibabel

hcp

mdtb