DICOM: 从参考坐标系到图像方位
Published:
本文介绍如何从 DICOM 格式文件获取图像扫描方位，并将位于参考坐标系（Reference coordinate system）的点转换到图像坐标系（the coordinate system of the image）。
背景介绍
最近在处理一批人体前臂的核磁共振（MRI）医学影像数据，包含两个序列：1. T1W 序列，用作解剖学参考；2. DTI 序列，用于追踪肌纤维。
我们的目的是：
- 在 T1W 的截面数据中划分肌肉轮廓，作为目标区域（Region of interest, ROI），并沿前臂长轴方向插值形成该肌肉的三维实体（类似于建模软件中的放样功能）；
- 利用 DTI 数据在 ROI 内追踪肌纤维走势；
- 找到位于实体内的肌纤维路径。
其中肌肉轮廓划分由 3D Slicer 软件实现，建模得到的肌肉实体导出为 stl 文件，导入 Matlab 做后处理；肌纤维追踪在 DSI Studio 中完成，生成的肌纤维路径保存在 .mat 文件中；最后判断肌纤维路径上的每个采样点是否位于三角形网格面片构成的实体内，找到真正属于这个肌肉的肌纤维点。
目前遇到问题是：将 3D Slicer 建模得到的实体、fiber 路径导入 MATLAB 后，发现两者之间存在夹角（图1左）。事实上，我们把肌肉实体导入 DSI Studio 中发现两者之间是没有夹角的（图1右）。这就使后处理流程的最后一步无法进行。
图1：校正前肌肉实体图示（左）与DSI Studio中实体示意图
为了展示效果，肌肉实体在图1中分别显示为蓝色（左）和黄色（右）。
思前想后，终于搞懂了！是坐标系转换的问题。具体细节记录在以下段落中。
必要的基础知识
我们首先来讲 DICOM 中定义的 “Reference Coordinates System (RCS)” 是什么。RCS 并不同于我们熟知的全局坐标系、局部坐标系等概念，RCS 是定义在受试者本身的。可以这样理解：我们想象受试者伸平双臂，从右手到左手方向便是 x 轴正方向；从前胸到后背的方向是 y 轴正方向；根据右手法则，从脚掌到头顶的方向是 z 轴正方向（图3左 Rxyz 坐标系所示）。
了解了 RCS 坐标系后，我们再来考虑图像坐标系（the coordinate system of the image）。当我们拿到一个影像序列时，通常会利用一些常用软件（比如 ImageJ）或设备自带的软件读取这些图像的头文件。其中代号为 0020 的模块（Image Plane Module）代表了图像位置信息（图2）。
图2：DICOM 图像位置模块标签
- (0020,0037) - Image Orientation (Patient) 是 1*6 的向量，前三个数字代表了图像 x 方向（row 方向）的基在 RCS 坐标系下的投影值，即为 Xx, Xy, Xz；后三个数字代表了图像 y 方向（column 方向）的基在 RCS 坐标系下的投影，即 Yx, Yy, Yz。
- (0020,0032) - Image Position (Patient) 是 1*3 的向量，代表了图像左上角第一个像素点中心在 RCS 坐标系下的坐标。
图3：两种受试者姿态示意图
看到这里是不是有些晕？别急，我们看图3。图三左侧的受试者姿态非常好，胳膊腿平行于躯干。当我们以矢状方向（sagittal direction, 图中黄色箭头）扫描时，我们得到的图像是平行于 Ry, Rz 构成的平面的，这时 (0020, 0037) 模块的值应当是：[0 1 0 0 0 -1]。但如果由于拍摄需求，或者受试者躺偏了，变成了图三右侧的姿态，我们仍然希望针对前臂拍摄矢状方向的序列，这时序列前进的方向不再是 Rx 方向了，而是 Ix 方向，得到的图像与 Ry, Rz 构成的平面形成了一定夹角。
回过头看图1。图1左侧肌肉实体的坐标显然是在 RCS(Rxyz) 坐标系，而肌纤维的坐标却是在 Ixyz 坐标系,所以存在一定夹角是必然的了。
坐标转换过程
下面我们来看如何做 Rxyz -> Ixyz 的坐标转换。 以下是本次数据集的一些信息：
IOP = [-0.0112, 0.9999, 4.623e-08;
-0.1732, -0.0021, -0.9849]';
其中 IOP 代表 Patient orientation。
通过叉积运算计算 Iz 方向：
IOP(:,3) = cross(IOP(:,1), IOP(:,2));
这样便得到了初步的旋转矩阵。
但观察图一左右，我们发现 IOP 矩阵给出的 x, y, z 并不是我们期望的 x, y, z。事实上，我们希望：Iz = -IOP(:,3); Iy = IOP(:,1); Iz = -IOP(:,2)。具体细节这里不做赘述。因此有：
TransM = [-IOP(:,3), IOP(:,1), -IOP(:,2)];
P_i = TransM\P_r;
其中 P_r 代表了图1左侧肌肉实体每个网格面片顶点的坐标（在 Rxyz 坐标系下），P_i 转换到 Ixyz 坐标系下的坐标。通过这一步骤就完成校正啦！校正后的图像如图4所示。
图4：校正后肌肉实体及肌纤维示意图
当然如果要将图片在 RCS 中的位置信息也加上，可以使用4维的旋转矩阵，即对 IOP 加一行 [0,0,0,1]，最后一列前三个数字为 (0020, 0032) 的数字即可。
希望大家能有所收获！
参考
http://dicomiseasy.blogspot.com/2013/06/getting-oriented-using-image-plane.html