基于快速最小距离树的肝脏图像非刚性配准

陆雪松，刘 坤，谢勤岚

(中南民族大学 生物医学工程学院，武汉 430074)

目前，肝脏疾病在我国呈高发态势.从医学图像中准确地分割出肝脏组织是疾病诊断、监测和计算机辅助手术的基础[1，2].通常计算机断层扫描(CT)由于高信噪比、高分辨率在临床上被采用，但人体肝脏在CT成像与相邻器官灰度值接近，特别与胃、心脏边界较模糊[3](如图1).故肝病的精确诊疗依靠耗时而重复性差的人工手动分割，灰度值的计算机自动分割难以满足临床需求，而亟需发展快速准确的肝脏图像自动提取方法.

1)肝脏; 2) 心脏; 3) 胃图1 肝脏与心脏(a)和胃(b)边界模糊示例Fig.1 Examples of the blurry boundary between liver and the heart(a) or stomach (b)

近年来，一些研究者对传统分割算法如GrabCut进行改进[4]，或者加入新的约束来提高分割精度[5]；也有一些研究者们尝试了多种交互式、全自动分割算法提取CT肝脏图像.Lu等[6]提出运用拟蒙特卡洛(Quasi-Monte Carlo)算法获取感兴趣区域的种子点和设计区域生长准则的肝脏CT图像分割方法，但只在针对低灰度等级的CT图像时有效.Li等[7]针对水平集演化时间长的问题，提出了一种改进的水平集分割方法，但需要手动初始化零水平集，效率不足以满足临床需求.宋晓等[8]提出了一种改进的Fast Marching肝脏分割方法，相对于其他水平集方法提高了分割效率和精度，但对于某些肝脏边缘误分问题还无法有效解决.Afifi等[9]提出一种利用相邻切片之间的空间相关性构建图割能量函数迭代分割整个序列的方法，但需手动分割初始切片，对于低对比度图像分割效果较差.Li等[10]提出一种基于形状约束和形变图割的肝脏CT图像分割方法，将形状约束融入到传统的窄带图割算法能量函数的区域项和边界项中可减少肝脏分割中的欠分割和过分割现象，但肝脏与血管间常出现较大表观距离.王雪虎等[11]针对传统形变模型在形变过程中易陷入局部最优的问题，提出一种利用传统形变模型三角面片作为新模型顶点来初始化新形变模型的分割方法，能较为准确的分割肝脏凹陷区域，但对于肝脏边界的肿瘤病变可能有误差.Slagmolen等[12]应用了一种非刚性atlas匹配方法，通过仿射和Bspline配准寻找atlas与目标灰度图像的形变关系实现肝脏的自动分割，但分割精度受配准质量的影响.

本文针对肝脏CT图像中的分割难题，选择基于atlas的全自动分割方法，采用最小距离树的快速算法执行非刚性配准，依靠图像的局部结构特征完成肝脏的自动分割.先分别对atlas灰度图像和待分割图像进行局部特征提取，再在全局粗配与局部精配相结合的形变框架下，采用Renyiα-entropy度量多维局部特征，实施两幅图像的快速配准；最终利用配准结果形变场对atlas的二值图像进行形变获得待分割图像.

1 基于atlas的肝脏自动分割框架

基于atlas分割方法的基本思想是将图像的分割问题转化为配准问题，选取已分割的图像作为atlas，并将其灰度图像与待分割图像进行非刚性配准，寻求atlas到待配准图像的形变场，最后将非刚性配准所得形变场应用到atlas二值图像上，得到待配准图像的分割结果，其整体分割框架如图2所示.

基于atlas的分割方法其分割精度较大程度上取决于配准的质量，针对配准，本文采用一种全局粗配与局部精配相结合的形变框架.其模型可描述如下：

f(x,y,z)=faffine(x,y,z)+fbspline(x,y,z),

(1)

其中：Taffine是仿射模型，作为整个配准过程中的全局粗配；Tbspline是基于B样条的自由形变模型FFD(Free-form deformation)[13]，用于局部的精确配准，前者的配准结果被用作于后者配准的初始参数.

图2 基于atlas的分割框架图Fig.2 The framework of segmentation based on atals

对于全局粗配的仿射模型可表达如下:

(2)

其中，12个形变参数可通过梯度下降法对传统互信息度量进行迭代寻优获取.而对局部精配的基于B样条的FFD模型可表达为:

bm(v)Bn(w)Øi+l,j+m,k+n,

(3)

其中Bl表示B样条函数的第l阶基函数.

在这种整体的配准框架下，采用两幅图像基于最小距离树的联合Renyiα-entropy[14]度量其局部多维特征，其中Renyiα-entropy可以通过K-最邻近获取密度函数来构造，但最小距离树直接计算熵值[15]法不需要估计概率密度，需要计算边的数量比前者少，故采用后者.基于最小距离树的联合Renyiα-entropy度量准则是一种通过计算浮动图像M和参考图像F重叠区域的联合最小距离树来度量图像特征的一种度量方法.将两幅图像的联合最小距离树表示为Lfm，则二者的联合Renyiα-entropy度量可表示为：

(4)

同时，通过获得Renyiα-entropy度量相对于FFD模型的梯度解析表达式并采用随机梯度下降法[16]进行寻优.

2 快速最小距离树算法流程

快速最小距离树算法FMST(fast minimum spanning tree)是由Zhong等[17]提出的一种基于K-means的快速算法，即通过K-means将原有数据集

分为K个子集，并对每个子集运用精确的最小距离树算法，最后将K个精确的最小距离树连接起来得到算法结果.其算法流程可描述如下：

Step1随机选取K个聚类中心对初始的采样点数据集X进行K-means聚类，得到K个子集，S={S1,…，Si},1≤i≤K；

Step2对K个子集中的采样点分别运用传统的最小距离树算法(Prim算法)构建最小距离树，得到K个精确的最小距离树MST(Si)，1≤i≤K；

Step3计算每个子集的质心ci，1≤i≤K，并运用Prim算法构建质心的最小距离树MSTcentre；

Step4遍历MSTcentre中的每一条边eij，其端点分别为ci和cj，分别计算Si中距离cj最近的点a与Sj中距离最近的点b，连接a和b，最终得到第一个完整的近似最小距离树MST1；

3 快速最小距离树算法时间复杂度分析

在产生第1个近似最小距离树时，在K-means的过程中可能将原本距离很近的一些点分割到2个不同子集中，这些点很靠近K-means边界，但在构造子集的最小距离树时将它们分开在2个子集中考虑.原本在精确的最小距离树中一个子集中靠近边界的2个点应分别与另外1个子集中的2个点相连，但由于聚类的原因，它们在自己的子集中被Prim算法连接在了一起，误差很大.针对这种情况，快速算法的第2个最小距离树以2个子集的质心最小距离树边的中点为中心对点重新分类，保证了原本被K-means错分的点被重新分类到1个子集，再通过Prim算法找出正确的边.因此，联合2个近似最小距离树后能极大限度地保证精确的最小距离树所选用的边都包含在这2个近似最小距离树中，最后用Kruskal算法对边进行筛选确保算法结果在较小的误差范围内.算法总的时间复杂度Ttotal可以表示为：

Ttotal=TK-means+TMST+Tcentre+Tcombine+TMST′+Tend,

(5)

其中TK-means,TMST,Tcentre,Tcombine,TMST′,Tend分别代表算法6个步骤的时间复杂度.

Step1和Step2总的时间复杂度明显取决于K的取值，令其时间复杂度总和为T′，因此需要求解最优K的取值T′.假定K-means聚类所得每个子集的点的数目大致相等，则每个子集点的数目为N/K.由于快速算法主要针对大数据应用，可以忽略度量特征维度与K-means迭代次数对算法时间的影响，因此，分割初始数据集K-means聚类时间复杂度可表示为TK-means=NK，则有：Prim算法时间复杂度为O(N2)，所有子集构造最小距离树的时间复杂度可表示为TMST=K(N/K)2.最终总时间T′=TK-means+TMST=NK+k(N/K)2.

Step4中，在最后连接时，需要对K-1对边所连接的2个子集进行遍历需找最近的邻接边，时间复杂度为O[2N×(K-1)/K]，即O(N).O(N).，故Step4最终时间复杂度可表示为O(N).

最后Step6由重复的Step3,Step4及Kruskal算法完成，Step3和Step4总的时间复杂度为O(N)，2个近似最小距离树边的总和最大为2(N-1)，故最后总时间复杂度为O(NlogN).

综上所述，算法时间复杂度可表示为：

Ttotal=TK-means+TMST+Tcentre+Tcombine+TMST′+Tend=N1.5+N+N+N1.5+NlogN=O(N1.5).

(6)

由此确定算法最终时间复杂度为O(N1.5).

4 算法结果与分析

为了验证引用方法的精度和速度，本方法实现在基于ITK(Insight toolkit)框架的开源配准软件包Elastix中，实验在临床的CT图像数据集3Dircadb中进行，该数据集包含了20组临床肝部图像，图像像素尺寸范围为0.56～0.87 mm，切片厚度范围为1～4 mm，切片层数范围为91～260，每张切片大小均为512×512像素，采用leave-one-out交叉验证法.将每组图像选作为待分割图像，剩下的19组图像作为atlas图像，共有20×19次配准.配准过程中，先进行全局Affine+MI粗配，再进行局部Bspline+MI精配，最后再进行本文方法的配准.其中Bspline+MI与本文方法，x轴和y轴平滑参数为8，4，2；z轴平滑参数为4，2，1，3分辨率下网格控制点间距均为20 mm， Renyiα-entropy度量α设置为0.99.

通过配准所获得的形变场将atlas的二值图像进行形变所得到的图像，可作为待分割图像的自动分割结果，通过计算待分割图像金标准与自动分割结果的Dice相似系数来测量二者的重叠率，并通DSC值来评估配准的精度，DSC值越高表明配准的精度越高.同时，通过求取不同数量采样点，在传统Prim算法与FMST寻优过程中平均单步运行时间比较FMST在速度上的优势.

为分析全局粗配(Affine+MI)、局部精配(Bspline+MI)和本文方法的精度，表1比较了三者19×20组配准结果DSC值中值和均值.

表1 3种方法配准结果重叠率对照表Tab.1 Overlap rate comparison of registration results from three algorithms

由表1可见：Bspline+MI的配准结果较Affine+MI在精度上有明显提升，其总的中值和均值均分别提高了0.2066和0.2319；而本法较Bspline+MI在中值和均值又分别提升了0.015和0.0092.

Bspline+MI和本文方法分割结果与金标准的对比结果见图3.由图3可见：在第1列图像中本文方法对于Bspline+MI肝脏下边缘过分割现象有一定改进，第2列中本文方法对于肝脏和心脏的模糊边界有更好的分割效果，第3列和第4列中本文方法对于Bspline+MI肝脏尖端下边缘的过分割也有一定的改善.

1)金标准; 2) Bspline+MI; 3) 本文方法图3 本法和Bspline+MI与金标准的分割效果对比图Fig.3 Comparision results of segmentation effect using our proposal and Bspline+MI compared with gold standard

表2比较了在不同采样点数目的情况下，Prim算法与FMST算法寻优过程中迭代的单步平均运行时间及相同采样点下平均时间减少率.通过计算每次寻优迭代第2～6次的平均值来获得平均时间，采样点数从5000～3000每次增加5000个采样点.由表2可见：在配准中FSMT算法较传统的Prim算法，前者速度提升十分明显，并随着采样点数目的增长，单步时间的减少率逐步提高，说明FMST算法面对大数据优势更明显.

表2 Prim和FMST算法寻优单步迭代平均时间对照

Tab.2 Average time comparison of the single step iteration by Prim and FMST algorithm

不同采样点个数/103 时间t / msMSTFMST减少率/% 5 34613114696.6910138847398597.1315310985469698.4920567250692098.7825865386871398.973012519261064199.15

基于多atlas的图像分割方法，其最终精度取决于atlas筛选、图像配准及标记融合方法的选择.本文着重改良配准算法，对最终分割结果的速度与精度有所帮助.传统的atlas选择方法主要有平均atlas、概率模型atlas及单个最优atlas.对于标记融合，传统方法主要利用一种类基于投票策略加权选择融合，后又提出许多新的算法，如STAPLE、SIMPLE等.对于每组图像的算法都会得到19组精度各异的自动分割结果，通过对这19组结果进行atlas筛选及标记图像融合[18]，即可得到待分割图像稳定精确的最终分割结果.

5 结语

本文针对肝脏CT图像的分割，采用基于配准的方法，引用非刚性配准框架将仿射模型和FFD模型相结合，并加入基于最小距离树的联合Renyiα-entropy度量其局部多维特征.同时在模型构建的过程中，引入新的FMST算法来计算其度量准则，并采用随机梯度下降法寻优，进行快速自动配准.结果表明：本文所用的模型较传统的仿射模型在精度上有较大提升， 引用FMST算法的使得配准的速度有较大提升.由于FMST算法中，K-means的初始聚类中心随机选取，当采样点较小时，寻优单步迭代的时间波动较大，但平均时间优于传统算法；若采样点较大，则不明显，故本法运用于基于无需采样的配准模型图像分割中具有较大的优势.

致谢感谢IRCAD Laparoscopic Train Center提供的数据支持