基于病理“金标准”标记脊髓损伤慢性期MRI图像坏死区域的方法

包星星，赵璨，包书圣，饶家声，3，杨朝阳，李晓光，

1.北京航空航天大学生物与医学工程学院/生物材料与神经再生北京市重点实验室，北京100083；2.中国康复科学所康复工程研究所，北京100068；3.北京航空航天大学生物医学工程高精尖创新中心，北京100083；4.首都医科大学基础医学院/神经生物学系，北京100069

前言

脊髓损伤是一种高致残率的中枢神经性疾病，严重影响患者的正常生活，并给家庭带来沉重的经济负担［1］。在脊髓损伤急性期阶段，通过创造适宜的微环境可以实现内源性神经干细胞在脊髓损伤局部向神经元的分化，从而实现神经组织的修复［2-5］。相比于急性期的研究而言，在脊髓损伤慢性期开展修复的研究相对较少，其中最主要的原因是胶质瘢痕的存在。胶质瘢痕的形成和稳定是脊髓损伤慢性期的病理生理学特征之一，在损伤局部所形成的致密胶质瘢痕会分泌硫酸软骨素蛋白聚糖，连同其所包绕而成的囊腔共同构成物理和化学双重屏障，严重阻碍了轴突的再生与髓鞘化，从而影响神经功能的恢复［6-12］。

对脊髓损伤慢性期的相关研究而言，胶质瘢痕及其所包绕的囊腔构成的脊髓坏死区域是脊髓组织修复过程中无法绕开的部位。针对特定时间点的横断面研究可以使用免疫组化的方法，通过病理切片来评估脊髓坏死区域的位置、尺寸、边界范围等信息。然而，针对长时间的纵向连续观察研究而言，病理切片获取坏死区域信息的方法存在着明显的局限性，单个受试对象无法既提供不同时间节点的观察结果同时又给出相应的病理切片。为了解决这一矛盾，一些研究采用磁共振成像（Magnetic ResonanceImaging,MRI）的方法来进行脊髓损伤后坏死区域的检测研究，保障了纵向实验的顺利进行［13-16］。

在上述研究中无论采用何种成像模态，MRI的脊髓坏死区域大多依赖于资深放射科医生的定性观察或手动识别［17］。考虑到MRI的部分容积效应，坏死组织在MRI图像中通常具有边缘模糊的特点，这使得坏死区域的准确识别存在一定难度。此外，不同损伤位置、程度、范围、方式等导致损伤区域各异，由此形成的脊髓坏死区域也形态多样且不规则，更增添了手动识别的困难度［17］。因此，基于操作者手动勾勒得到的脊髓坏死区域往往不具有可重复性和可验证性，阻碍了对胶质瘢痕及囊腔区域的准确评估。

目前已知病理切片是识别脊髓损伤慢性期胶质瘢痕及其包绕的囊腔区域的“金标准”，因此，应用病理切片结果来指导MRI图像坏死区域的识别无疑是一个较好的思路。要实现上述工作，首先应当实现病理切片与MRI图像的对应。本研究基于病理切片与MRI图像的坏死区域空间相对位置具有一致性这一特点，通过一系列图像处理建立了两者之间的映射关系，并验证了映射后获得的MRI图像实际坏死区域与“金标准”的胶质瘢痕和囊腔范围之间的吻合度，为将来计算机自动识别与分割MRI坏死区域奠定了基础。

1 材料与方法

1.1 动物模型制备

本研究制备了9 只成年雌性Wistar 大鼠的脊髓砸伤模型，砸伤位置为脊髓的T8 节段，使用NYU 脊髓打击器（NewYork University,New York,USA）砸伤脊髓，其中撞击棒重10 g，并从50 mm 的高度自由下落完成施力。术后3个月内，实验大鼠饲养于恒温恒湿的环境中，确保充足的食物和水。术后连续5 d 注射青霉素（10 万单元/d/次），每天人工护理排尿直到实验动物恢复自主排尿功能。实验方案经过北京航空航天大学生物与医学伦理委员会批准（批准号：BM20180047）。

1.2 数据采集

本研究在大鼠脊髓损伤慢性期（术后3 个月）采集了其脊髓矢状面的T1加权（T1-weighted,T1W）和T2加权（T2-weighted,T2W）两种在体MRI结构像数据，所有数据均采自首都医科大学7.0T 小动物磁共振扫描仪（Bruker Biospec, Karlsruhe, Germany）。坏死区域在T1W 图像中主要表现为低、中等信号强度，在T2W 图像中主要表现为高信号强度。不同模态下各异的影像学表现为脊髓坏死区域的定位提供了丰富的信息。

1.3 病理切片制备

脊髓病理切片是确定脊髓坏死区域尺寸和范围的“金标准”，是建立脊髓病理切片到MRI图像映射关系的基础。本研究对实验动物脊髓矢状位的病理组织切片（10 μm）进行了胶质纤维酸性蛋白（Glial Fibrillary Acidic Protein, GFAP）染色，并通过普通光学显微镜获取了GFAP染色后的病理切片图像。图1的病理切片展示了脊髓损伤慢性期神经元坏死液化形成的空洞（即囊腔），囊腔壁由致密的胶质瘢痕构成。

1.4 方法

图1 为本研究所提出的方法示意图。为建立脊髓病理切片到MRI图像的映射关系，本研究所构建的方法包括以下几部分工作：（1）预处理：双边滤波去噪以及T1W 和T2W 图像的刚体配准；（2）脊髓分离及拉直：利用Snake模型实现MRI图像中脊髓的粗分割，剔除冗余信息，并利用双三次插值法对脊髓进行拉直；（3）图像融合：基于双树复小波变换（Dual Tree-Complex Wavelet Transform,DT-CWT）对T1W 和T2W图像进行融合；（4）脊髓MRI图像坏死区域的确定：基于病理切片与MRI图像的脊髓坏死区域空间相对位置具有一致性这一特点，建立脊髓病理切片到MRI图像的映射关系，以实现MRI图像中脊髓坏死区域的标注。

图1 脊髓坏死区域鉴别框架示意图Fig.1 Schematic diagram of spinal cord necrosis area differentiation framework

1.4.1 预处理MRI图像在获取和生成的过程中会混杂一定的噪声，这无疑会给图像的后续处理增加干扰信息。双边滤波方法在进行噪声平滑的同时能较好地保留图像的边缘信息［18］，因此本研究将该方法应用于脊髓MRI多模态图像的去噪，取得了较好的去噪效果。

T1W 和T2W 两种序列的采集存在时间差，生理因素及机器震动等均有可能导致组织发生空间位移，从而引起图像的不匹配。本文以T1W 图像为参考图像，对T2W 图像进行刚性变换以改变T2W 的空间位置，实现T1W 和T2W 图像的配准，以消除两种序列的空间位置差异。

1.4.2 脊髓分离及拉直原始脊髓MRI图像不仅存在脊髓，还包括椎体、椎间盘、肌肉以及皮肤等其他组织，为去除冗余信息，本文采用了Snake 模型从MRI图像中提取脊髓。Snake 模型在最小化能量泛函的动力下驱使曲线不断靠近目标物体边缘以获取目标物体轮廓［19］。初始轮廓曲线是通过将T1W 图像进行对数变换后手动获取的，根据曲线上点v在图像中的空间位置可确定Snake模型的能量函数：

其中，Etotal为曲线演化过程中的总能量；第一项为内部能量，保持着曲线的连续性和平滑性；第二项则为外部能量，表示灰度值或梯度等图像特征信息对曲线的能量控制，本研究取的是v点处图像的梯度信息：

由于在实际制备脊髓病理切片的过程中可能会产生任意角度的弯曲，为准确获取病理切片与MRI图像的空间对应关系，本研究采用插值效果最佳的双三次插值法对病理切片和MRI图像中的脊髓均进行了拉直处理。脊髓拉直图像中某位置的强度值，对应于距离原图像同位置像素点最近的16个像素强度值的加权平均值，各权重的计算如下［20］：

其中，x表示原图像中距离给定数据点最近的16个像素点到目标图像对应点的距离，S(x)表示16 个像素点对应权重。

1.4.3 图像融合为尽可能全面地获取脊髓坏死区域的影像信息，本研究充分考虑了脊髓T1W 和T2W 图像的信息，利用DT-CWT 方法将T1W 和T2W 图像转换到了频率域，在频域实现两种图像的融合。DTCWT 不仅保留了小波变换的时频分析特性，而且具备平移不变性、方向分析能力良好、数据冗余有限、计算效率较高以及重构效率好等特点，成为近年来应用于图像融合的热门方法之一［21］。图2 为DT-CWT 的分解示意图，通过树A 和树B 并行获取实部和虚部的小波变换系数，从而形成双树结构。每级分解通过低、高通滤波及间隔采样，可获取两个低频子图像以及±15°、±45°和±75° 这6个方向的高频子图像1, 2,…, 6)。本文首先将脊髓拉直的T1W 和T2W图像进行DT-CWT 五级分解，然后对同方向高频系数的融合取小波系数模的最大值，以使融合图像突出细节信息。低频系数的融合取其小波系数的均值，再经过逆变换可得到T1W 和T2W 融合图像（图3）。

图2 DT-CWT变换Fig.2 DT-CWT

图3 融合流程图Fig.3 Fusion flow chart

1.4.4 基于病理“金标准”标定脊髓MRI图像坏死区域的方法本研究的目的旨在构建病理切片与MRI图像的映射关系，从而依据病理切片“金标准”来客观地标记脊髓MRI图像中的坏死区域。图4 为映射关系构建的示意图，图4a 为拉直的脊髓病理图像，包含坏死区域的脊髓段已用黑色框标记并对其进行裁剪。病理图像中坏死组织的边界清晰，可手动获取坏死区域并进行二值化，将损伤区域赋值为1，其余赋值为0，得到图4b 所示二值图像。针对一张MRI图像（1 mm 层厚）对应多张病理切片（10 μm 层厚）的情况，本研究分别对多张病理切片进行二值化处理，随后分别与MRI图像建立映射关系，最后再求取总的映射关系式，以尽量减小瘢痕在MRI图像上被标记的误差。映射关系的建立是依据被拉直的病理切片和MRI图像中坏死组织边界到脊髓左侧边缘的水平距离占整个脊髓宽度的比例相等这一特点来实现的（研究中所用的距离、宽度和高度都是指相应方向的像素个数），具体公式如下：

其中，s为对应的病理切片数。对于MRI图像中脊髓的裁剪，首先根据脊髓的影像学表现，可确定图4c脊髓宽度n以及大致的损伤中心位点，再计算二值化图像的高度、宽度比例的均值可获取应裁剪的位置和高度m：

位置的标记分两种情况：若在病理图像中，瘢痕区域的数量等于1，则选取损伤区域y 轴方向上1/8、1/4、1/2、3/4 和7/8 这5 个位置，每个位置一对像素点，共10个像素点（分别对应同一y轴水平损伤边缘的两个位置）；若瘢痕的区域数量大于1，则增加3/8、5/8、1/3 和2/3 这4 个位置，共18 个像素点，最后连接这些像素点以勾勒整个瘢痕边界。

图4 病理图像与MRI图像映射关系构建示意图Fig.4 Schematic diagram of mapping relationship between pathological image and MRI image

1.4.5 评估方法传统评估方法在本研究中难以适用，这主要有两个方面的原因：（1）传统评估方法不适用于图像尺寸不一的配准情况；（2）本研究使用了二值化图像和相对位置关系两个要素，因此难以提供额外的图像纹理、灰度以及空间位置等信息供传统评估方法所使用。根据本研究的具体情况，笔者建立了一种基于形状参数的评估方法，包括圆度和偏心率两个特征，计算公式如下：

其中，PR 和PE 分别表示病理图像中胶质瘢痕及囊腔区域的圆度和偏心率，IR 和IE 则分别为MRI图像中坏死区域的圆度和偏心率。式中S为瘢痕面积，C为瘢痕周长，c为瘢痕区域焦点间的距离，a为长轴长度。

本研究通过比较病理切片区域和MRI图像中两个区域的圆度和偏心率特征来反映MRI图像坏死区域标记的准确性，计算公式如下：

其中，RR、ER 分别为病理图像中胶质瘢痕及囊腔区域与MRI图像中坏死区域的圆度、偏心率的比值。在计算时，将较小的值设置为分子，较大的值设为分母，MRI图像中所标记的实际坏死区域与病理切片所给出的“金标准”越相似，则RR 和ER 的值越接近1。当两者完全相同时，计算结果即为1。

2 实验结果与分析

2.1 脊髓分离结果

Snake 模型能量函数的外部能量大小与曲线上点的局部梯度大小呈负相关，图5a 中，脊髓正常组织与瘢痕组织的信号突变使得曲线的演化方向偏离脊髓边缘（图5b），导致最后获取的脊髓不完整，这将影响后续病理图像与MRI图像映射关系的准确建立。而图5c 中，T1W 图像通过对数变换后，降低了脊髓正常组织与瘢痕强度的对比度，最终可获取完整的脊髓（图5d）。

2.2 图像融合结果

表1 对T1W 和T2W 图像在不同DT-CWT 分解级数下融合效果进行了比较，评估的参数包括互信息（Mutual Information,MI）、熵（Entropy）以及峰值信噪比（Peak Signal to Noise Ratio,PSNR）。综合来看，在五级分解的情况下可得到最优的融合效果。由此，本研究采取DT-CWT 五级分解再融合进行后续的实验。

2.3 MRI图像中实际坏死区域标注的结果

通过建立病理切片与MRI图像的映射，本研究实现了在MRI图像上准确标记坏死区域（图6）。

另外，本研究对病理图像中的实际胶质瘢痕及其包绕的囊腔区域与MRI图像中标记的坏死区域在形状上的相似程度进行了量化分析，结果发现，这两种区域的相似度评估结果在圆度上为0.93±0.03，在偏心率上达到了0.97±0.02，表明了MRI图像中所标记的区域与病理“金标准”具有较高的吻合度。

图5 脊髓粗分割Fig.5 Coarse segmentation of spinal cord

表1 基于不同DT-CWT分解级数下融合效果统计Tab.1 Statistics offusion effect based on different DT-CWT decomposition levels

图6 MRI图像中脊髓坏死区域的标记结果Fig.6 The results of MRI in the necrotic area of spinal cord

3 总结与讨论

本研究提出一种应用病理切片来确定脊髓损伤慢性期MRI图像中实际坏死区域的方法，基于病理切片胶质瘢痕及其所包绕的囊腔区域与MRI图像的脊髓坏死区域具有空间相对位置一致性的特点，建立脊髓病理切片到MRI图像的映射关系，据此实现MRI图像中脊髓坏死区域的客观、准确标记。利用形状信息的评估方法显示MRI标记的坏死区域与病理学“金标准”的胶质瘢痕和囊腔区域之间高度相似，为今后应用图像识别技术构建计算机自动辅助分割脊髓损伤慢性期坏死区域的方法奠定基础。

后续的研究可以在两个方面继续深入开展：（1）进一步改进映射方法，采用加权求和的策略，结合主观经验，将那些与MRI图像中坏死区域的大小、形状及位置等较为相似的病理切片样本数据赋予更大的权重值，以获取带有加权系数的映射关系；（2）在已建立映射关系的基础上，评价MRI中已标记坏死区域的三维结构形态与病理学“金标准”三维结构形态之间的差异，从而更全面地评估映射关系的准确性。