基于重采样改进的多图谱分割算法

夏 瑞，马 瑜，王文娜，罗宇卓，尚梦玉

(宁夏大学 物理与电子电气工程学院，宁夏 银川 750021)

0 引 言

阿尔兹海默症[1](Alzheimer’s disease，AD)的医学诊断表现为颞中回萎缩，海马萎缩是阿尔兹海默症早期诊断的标志。研究人脑核磁共振(MR)图像中海马体的分割，有助于医生了解患者海马体病变状况，同时益于医生结合分割结果给出的定位信息精确地执行手术计划。

由于人脑结构的复杂性以及脑组织结构图像低对比度等原因，传统分割方法不易得到精确分割结果。基于多图谱配准的分割方法利用多个图谱的先验知识，降低单一图谱配准产生的不确定性以及误差，得到的分割精确度更高。三维医学图像配准分割技术，得到的分割结果能为诊断提供丰富的信息，缺点是需要较长的计算时间。因此，本文提出基于重采样改进的多图谱配准方法，两次提取感兴趣区域(region of interest，ROI)，并改进标签融合算法，有效提高配准时间效率，同时保持较高配准精度。

1 基于多图谱配准的海马体分割方法

1.1 图谱图像重采样原理

传统的多图谱配准算法大多采用“粗精”混合的配准方法，先采用全局刚性配准进行图谱与目标图像的对齐，缩小两者之间的差异，再通过精细的非刚性配准得到精确的配准结果。直接用精细配准无法使得两幅图像很好地对齐，而非刚性配准过程要求图谱图像与目标图像具有相同的空间采样点和采样点空间距离，可对齐两幅图像。

由于3D医学图像具有巨大的数据量，图像之间的配准过程消耗计算时间过长，因此“粗精”混合的配准方法不能同时达到高精度、高时间效率的要求。基于以上原因，本文提出利用重采样的方法代替“粗”配准阶段，重采样过程包括输入图像、变换和校对机，图像的空间坐标通过变换进行映射以便生成一个新的图像。

一幅图像是使用离散网格对连续场的采样，重采样就是基于灰度的图像配准的本质。重采样具有和“粗”配准一样可以达到图谱图像与目标图像对齐的目的，使图谱图像和目标图像有相同的空间采样点和采样点空间距离，同时减少了“粗”配准环节所需时间。在重采样过程中，当图像从一个空间到另一个空间的映射需要在非网格位置进行操作时，校对机用于计算非网格位置的亮度值。变换用来映射点坐标，包括旋转、缩放、平移等，本文采用的变换是仿射变换。重采样从参考图像获得输出信息，根据提供的参考图像，使用参考图像的间距、原点和方向对待采样图像进行映射。进行重采样后，调整了浮动图像和参考图像采样点不一致的问题，并且图像的尺寸大小得到矫正。

重采样操作使参考图像与浮动图像具有同样的各向同性采样率，图谱图像与目标图像大小、中心达到一致，为后续进行非刚性配准做好准备。图1是一个简单重采样的原理图，重采样涉及从原始图像中提取像素的位置、插值灰度级，并将其重定位到校正图像中的近似矩阵坐标位置。内插的方法包括最近邻，双线性插值和立方卷积。如图1所示，第一行是原图像素空间采样点、采样空间大小示意图，第二行是对图像像素进行重采样，改变原点坐标且采样空间大小为原来的两倍的图示。

图1 重采样原理

1.2 微分同胚Demons配准算法

传统Demons算法，空间变换的复合运算使用向量加法来近似，然而，图像配准中用到的大多数空间变换形成的是李群，而不是向量空间。同时，传统Demons算法不能到可逆变形场，易引起局部形变折叠。而微分同胚是可逆的光滑映射，微分同胚Demons算法[2-5]，能够保证图像的拓扑结构在配准前后保持不变，光滑且连续，并防止引入形变折叠，对大小形变都适用，在缺少可用的空间变换信息时，是很好的配准框架。

微分同胚Demons算法，利用公式c←s∘exp(u)在李群中对几何变换s进行更新。微分同胚Demons算法的目标能量函数为

(1)

其中，变形场u是一个稠密速度场

(2)

微分同胚Demons算法的实现过程：

步骤2 对u进行流体正则化，即u←Kfluid*u；

步骤3 利用牛顿方法计算李群得到的exp(u)，计算c←s∘exp(u)；

步骤4 对s进行扩散正则化，s←Kdiff*s。

1.3 多分辨率配准

多分辨率配准[6]是在不同分辨率条件下用不同的比例进行配准，并且上一级的配准结果用作下一级配准的输入，直到达到最好的比例范围。多分辨率配准能够提高配准的成功率，并且在粗糙比例时消除局部噪声。多分辨率配准框架如图2所示，参考图像和浮动图像被视为两个图像金字塔，金字塔根据用户定义的收放因数来进行平滑和二次采样图像。

图2 多分辨率配准处理概念

本文将多分辨率配准思想加入到微分同胚配准过程中，既提高了配准的精确度，又降低了配准的时间。

1.4 基于非局部图像块相似性的加权选择算法

多图谱配准得到的形变结果传递到图谱相应的标记图像，即可得到相应的分割结果，再将多个图谱的分割结果进行融合得到一个综合的分割结果。当前研究较为成熟的标签融合算法有：加权选择算法、STAPLE算法、COLLATE算法等。本文以K近邻搜索算法对加权选择算法进行改进，并在融合阶段再次提取ROI，使融合算法运行更高效。

1.4.1 加权选择标签融合算法

简单加权选择算法[7,8](majority voting)考虑大多数分割结果在该采样点是否存在，如果大多数分割结果在该点的像素值为1，那么判定该点为分割结果。该算法不需要先验知识，对每个像素进行选择，假设图谱与目标图像进行配准，并将形变向量传递给标记图像后，图谱表示为Ai,Li，对于空间中的每一个位置x，目标图像的分割结果为

(3)

(4)

其中，i为图谱的个数，Ltarget(x)表示在位置x处的目标图像的标签像素。

改进后的加权选择算法，充分考虑到传统融合算法中未用到图谱灰度图像中信息这个缺陷，得到如下的融合公式

(5)

(6)

其中，权重系数Wi是目标图像与形变后图谱灰度图像的相似性，一般为归一化互信息(NMI)、归一化互相关系数(NCC)、平均误差平方和(MSD)。融合算法中引入了图谱灰度图像的先验信息，在一定程度上能够减小误差，得到更准确的融合结果。

为了更精确地应用目标图像与图谱灰度图像的先验信息，一些学者研究了基于局部块相似性的加权融合算法，其思想是，在计算相似性时，将图像中每个体素用以该体素为中心的一个图像块来代替，图谱图像块与目标图像块之间的相似性越高，则其拥有相同标签的概率越大

(7)

(8)

υ(x)是目标图像在位置x处的标签值，W(x,xs)是以x为中心的目标图像块与第s个图谱中以x为中心的图谱图像块的相似性，Lxs是第s个图谱中在位置x处的标签值，h是局部适应参数。

考虑到配准可能产生误差，为了减弱这种误差对融合结果的影响，Pierrick Coupe等[9]研究了基于非局部块相似性的加权选择算法，克服了局部块加权的缺点，通过在对应块的邻域进行搜索，增加了图像信息，来计算相似性

(9)

(10)

V是图谱在位置x处的一个搜索邻域块，Lxs,j是第s个图谱在搜素邻域V中位置j处的标签，W(x,xs,j)代表以x中心的目标图像块与以xs,j为中心的图谱图像块的相似性。

对于目标图像的同一个体素x，对应的图谱图像增加了块的搜索范围，则搜素邻域内图像块的个数变多，将相似度较低的块的权重直接置为0，可以排除不相似块对整体权重的影响。Pierrick Coupe等[9]在计算权重时，选择如下相似性公式对图谱块进行预选

(11)

其中，μ和σ是在位置x处目标图像块的均值和方差，μs,j和σs,j是第s个图谱在x的搜素邻域中位置j处的图谱图像块的均值和方差。s值越大，两个块越相似，当s小于预先设定的阈值时，舍弃与该s对应的图谱块，不进行权重计算。

1.4.2 基于K近邻搜索的加权选择融合算法

为了提高基于非局部块加权融合算法的效率，本文进行了更深入的研究，通过引入K近邻搜索算法，快速搜索图像块的k个近邻，只用k个非局部块来进行局部加权。

K近邻搜索算法是一种有监督的学习分类算法。对于一个样本，寻找与它最相近的k个邻居样本。实现K近邻搜索算法可以用kd-tree或VP-tree。

k-d树(k-dimensional树的缩写)是一种对K维空间中的实例点进行存储以方便对齐进行快速检索的树形数据结构，k-d树可以用于多维空间关键数据的搜索，例如范围搜索和最近邻搜索。

最简单的最近邻搜索，也称为穷举搜索，依次计算样本集中每个实例点到目标点的距离，然后取最小距离的那个点。当样本集较大时，这种策略非常耗时。对于大规模的高维数据空间，kd-tree是一种有效的快速搜素方法。利用kd-tree搜索最相近的k个数据点，减少了搜索全部数据点的计算量。

因此，本文将权重公式改写如下

(12)

xs,j∈K表示xs,j在x的搜索邻域V中，并且是x的k个近邻，这样这个图像块才被进行权重计算。由于非局部块加权融合算法增加了搜索范围，导致图像块的数量增加，加入K近邻搜素算法，只使用最近的k个非局部块来进行局部加权，大大减少了计算量，提高了融合算法的效率。

下面是基于K近邻搜素的非局部块加权融合算法实现过程：

步骤1 标注所有配准形变之后的图谱标记图像的连通区域，获得包含所有标注区域的一个区域块

Vo=V1∪V2∪…Vn

步骤2 根据步骤1中计算得到的立方区域块Vo的大小裁剪目标图像T得到TROI，以及图谱图像Ai,Li其,中i=1,2,…,N为参与配准的图谱数量，得到AROI,LROI；

步骤3 K近邻搜素算法对图谱块进行预选，剔除与目标图像块相差较大的图谱块，减少更多的计算量；

步骤4 构建概率图谱，获得最终的分割结果。

1.5 分割结果的评价

图像分割结果的好坏需要评价标准来进行衡量。主观评价即观察者主观视觉观察，易受环境等因素影响和观察者的主观认知影响。常用的客观评价方法为相似性测度(Dice)，计算两个二值分割图像的重叠程度，式(13)来评价目标图像的专家手动分割结果SF和浮动图像分割结果SM重叠的部分

(13)

∩表示两个集合的交集。相似性测度因子d取值范围是0到1，0表示两分割区域完全不重叠。1表示两分割区域完全重叠，即配准分割结果与金标准一致，分割结果十分精确。

2 实验结果分析

本文采用Visual Studio 2010软件平台结合ITK开源代码库进行MR人脑图像海马体的自动分割。数据来源于伦敦帝国理工学院医学脑部研究数据库[10]，共20个T1-MR灰度图像和对应的标记图像样本，灰度图像标号为a01到a20，对应的标记图像标号为a01-seg到a20-seg。其中，标记图像中标记了人脑中67个结构。

算法流程如图3所示。

图3 本文算法流程

图谱的选取在一定程度上影响着分割结果的好坏，当与目标图像差别较大时，配准结果越差，得到分割结果越偏离金标准。P.Aljabar等[11]在文章中描述了相似性测度因子和图像的数目存在如下关系

(14)

a、b是常数，满足0≤a≤1,b>0，n是图谱的数目。图谱数目的增多反而使分割结果与金标准的相似度下降。Sabuncu等[7]的研究表明，若随机选择图谱，Majority Vo-ting的表现不是单调的，且图谱数量超过10精确度开始下降。若选择最优图谱，图谱在15个左右时Majority Voting的结果变得平稳。基于上面几条结论，综合本文的实验，同时考虑标签融合的效率问题，本文选取了与目标图像互信息最高的10个样本作为图谱图像进行实验，对分割精度的影响较小。

2.1 图像预处理

2.1.1 颅骨剔除

人脑MRI中背景和非脑组织占有较大比重，容易造成配准的误差。剔除脑壳是重要的预处理工作，本文利用Stefan Bauer等[12]提出的算法剔除脑壳，克服了传统BET算法运算时间长的缺点，并可以有效剔除影响配准精确度的脑壳与背景部分。

图4 颅骨剔除结果

2.1.2 重采样并提取待分割组织ROI

本文以灰度图a04、标记图a04-seg为参考图像，对图谱进行以参考图像为基准的重采样操作。参考图像的尺寸大小为195×198×170，体素间距是0.937×0.937×0.937，图像中心为(90.89, 92.29, 79.18)。例如，当图谱a01的尺寸大小、体素间距、图像中心分别为176×198×160、0.937×0.937×0.937、(81.99, 92.29, 74.49)，则重采样后a01具有和参考图像一样的大小、间距、图像中心。

对目标图像以及图谱图像进行初次提取感兴趣区域(ROI)。如图5所示，为灰度图谱及对应的标记图像初次提取ROI的结果。

图5 以左、右海马体为中心提取ROI

2.2 配准及标记融合结果

本文以重采样结果取代刚性配准环节，并以提取的ROI作为微分同胚Demons配准对象，大幅度缩短配准时间。微分同胚Demons配准产生的一个形变场的示意图，如图6所示。

图6 待配准图像映射到目标图像的变形域

标签融合的结果即分割结果。图7为目标图像右海马体的金标准(专家手动分割结果)与本文方法分割结果比较图：第一行至第三行分别为轴向、矢状面、冠状面海马体切片图。每行左图为目标图像金标准；右图为分割结果与金标准的覆盖比较图，其中白色区域为两者重合区域，灰色边缘为两者非重合区域。从图中可以看到本文分割方法得到的海马体覆盖了大部分金标准结构，即分割结果较为精确。

图7 海马体切片比较

如图8所示，第一行是目标图像右海马体的金标准；第二行是传统方法分割结果(左图)以及本文方法分割结果(右图)的三维显示效果。可以看到在分割细节上，本文方法更接近于专家手工分割，传统方法分割结果表面过于光滑，丢失了分割对象的部分表面纹理信息。

图8 右海马体分割结果三维显示效果

如表1所列，是以a04为目标图像的海马体分割精度比较。改进的加权选择算法比简单加权选择算法精度略有提高，而本文方法得到的左右海马体的分割精度均高于简单加权选择算法、改进的加权选择算法，Dice值有10%～20%的提高。

表1 分割结果与目标图像的相似性测度

以右海马体分割过程耗时来进行分割效率分析，表2、表3中Resamp表示重采样、ROI表示提取感兴趣区域、GReg表示刚性配准、DReg表示微分同胚Demons配准、Wrap表示形变映射、Fuse 表示标签融合、Total是总时间。如表2所列，是使用本文方法分割海马体时各操作过程所耗时间，表3所列，是使用传统“粗精”配准方法分割海马体所用时间。

表2 本文方法分割右海马体时间/s

表3 传统“粗精”配准法分割右海马体时间/s

从表1、表2可以看到，分割以a04为目标图像分割右海马体时，本文方法重采样耗时仅4 s～5 s,相比于传统“粗精”混合算法中粗配准时间387 s～781 s，效率提高了95～155倍；并且本文方法提取ROI后，微分同胚Demons配准时间为2 s～3 s，而传统方法中微分同胚Demons配准用时147 s～151 s，精配准的效率提高了49～72倍；在形变域映射过程，相较于传统方法中7 s～8 s的用时，本文方法仅为0.5 s～1 s；最后，标签融合阶段，改进之后的加权选择算法比传统方法减少了40.1 s。整个多图谱分割过程，本文方法所用时间共117.7 s，而传统方法用时为6852.3 s，大大提高了分割的效率问题，实现了在两分钟之内完成人脑MR图像中海马体的分割。

3 结束语

三维医学图像配准分割是医学图像研究的热点，大多数研究者在分割精度上进行了深入的研究，并取得了较好的结果。处理三维医学图像耗时较长，人们在专注研究分割精度时往往忽略分割时间效率，甚至以牺牲时间效率来提高分割精度，很难适应于临床应用对于快、精、准的要求。基于此，本文对基于多图谱配准的海马体分割算法进行研究，针对分割效率问题，改变传统“粗精”混合配准方法，两次提取人脑海马体ROI，在重采样基础上进行微分同胚Demons配准，并改进标签融合算法。实验结果表明，本文方法在分割精度上比传统方法提高了10%～20%左右，在分割时间效率上是传统“粗精”配准分割算法速度的58～60倍，精度与效率均有提升。后续将考虑用实际数据来对本文算法进行验证，继续完善算法，以提高精度。