血管粘连型肺结节图像的序列分割方法

张 伟，张小龙，赵涓涓+，强 彦，唐笑先

(1.太原理工大学 计算机科学与技术学院，山西 晋中 030600；2.宾夕法尼亚州立大学 信息科学与技术学院，宾西法尼亚州 尤尼弗西蒂帕克 16802；3.山西省人民医院 PET/CT中心，山西 太原 030024)

0 引 言

肺癌是当今发病率和死亡率最高的恶性肿瘤[1]。早期肺癌的5年生存率超过70%[2]。因此，早发现早治疗就显得尤为重要。而肺部病灶中，血管粘连型结节的恶性度非常大，其早期CT影像的病变区域却难以被准确的分割，导致后续的识别和诊断出现误差。同时，为了提高肺结节的检出率，在CT扫描中常常采用超微小间隔CT扫描，因而又出现影像数据爆炸式增长与人工诊断力量严重不足的问题。因此，当前迫切需要一种既高效又能准确分割该类型结节的方法。

1 相关工作

肺结节图像的准确分割是后续特征提取和良恶性诊断的基础。因此，近年来对肺结节图像分割算法的研究日益受到人们重视。

Sun等[3]提出一种流向特征熵和测地线距离的K均值聚类算法来分割肺结节，取得较好的结节分割结果，但是随着结节半径增加，算法执行效率就越低；Wei Y等[4]先利用多尺度Hessian矩阵对疑似结节进行增强，再通过设计mean-shift聚类的核函数，最终实现结节疑似感兴趣区域(region of interest，ROI)的分割；Sivakumar S等[5]提出一种基于无监督分割模型的possibilistic-fuzzy聚类算法，提高了肺结节分割的速度和准确率；邹瑜等[6]针对自组织特征映射(self-organizing feature map，SOM)网络随着神经元增加分割效果变差的问题，提出一种将有限脉冲响应加入SOM中的方法，并通过合并聚类实现对乳腺超声检查图像的有效分割；赵涓涓等[7]利用PET图像中SUV阈值结合分水岭分割实现了孤立性肺结节图像的分割，分割效果较好。

虽然这些算法对实性结节的分割都很有效，但是血管粘连型结节的灰度值和血管的灰度值非常接近，分割结果中不能有效将血管和结节分离开。

因此，本文针对血管粘连型结节的序列CT图像，提出了一种序列分割方法。与文献中提到的其它方法进行比较，验证了本方法的有效性和通用性。

2 CT图像预处理

2.1 序列肺实质图像分割

该部分使用此前项目组提出的序列肺实质分割方法[8]对CT序列图像进行肺实质分割。所使用的方法分为4个阶段：第一，利用肺实质图像的位置特征得到肺部CT序列ROI图像；第二，利用改进的超像素序列图像分割算法对序列ROI图像进行超像素分割；第三，采用遗传算法(genetic algorithm，GA)优化的自生成神经森林(self-generating neural forest，SGNF)算法对超像素样本进行聚类；第四，根据聚类后超像素集的灰度特征和位置特征识别肺实质。

2.2 感兴趣区域提取

该部分采用Song J等[9]提出的改进的雪橇算法实现肺结节的准确定位，并以此结节为中心，提取30×30的感兴趣区域ROI，排除较大部分背景噪声干扰的同时显著提高后续图像分割的效率。

3 改进的超像素序列分割

超像素[10]是由一系列位置相邻且颜色、亮度、纹理等特征相似的像素点组成的小区域。其把具有相似特征的像素集合起来进行操作，降低了图像处理的复杂度。

本文充分考虑CT图像中血管和结节类型的特征，提出自适应相似度系数的超像素序列分割算法(adaptive similarity coefficient sequential linear iterative clustering，ACSLIC)，算法使用Lab颜色空间中的L，a，b这3种颜色信息和像素点的x，y坐标构成5维特征向量[L，a，b，x，y]T来表示CT图像的单个像素点，并在像素间距离计算中引入自适应相似度系数，设计新的适用于肺部CT图像像素间相似性衡量的距离计算方法，如式(1)～式(3)所示

(1)

(2)

ds=w*dLab+dxy

(3)

其中，ds表示系数叠加后CT图像两个像素点间的相似度；dLab表示像素与聚类中心间的颜色特征距离；dxy表示像素与聚类中心间的空间位置特征距离；k是超像素的聚类中心，i是其搜索范围内某一像素点；w是k、i像素点间的自适应相似度系数，其由式(8)所确定。

在CT图像中，假设λ1和λ2为Hessian矩阵的两个特征值，则线性结构满足条件λ1≈0，λ2≫0，圆形结构满足λ1≈λ2≫0。并且依据肺结节的类圆形特征和尺度不一特征。本文将Hessian矩阵与高斯函数相结合，通过设计多尺度的结节相似度函数来对CT图像中的所有像素进行结节相似度值的计算，从而构造自适应的相似度系数w。其很大程度上拉大了血管和结节区域的相似度，有效避免了血管和结节被分割到同一个超像素样本中，极大提高后续肺结节分割的准确率。如式(4)～式(8)所示

(4)

(5)

(6)

(7)

(8)

式(4)为2D图像的Hessian矩阵定义；式(5)和式(6)实现圆形结构的多尺度融合，其中σmin=dmin/4，σmin=dmax/4，dmin为肺结节图像的最小直径，dmax为肺结节图像的最大直径，实验中设置dmin=3和dmax=30，N取8，并通过不断的迭代尺度因子σ，选取各个尺度上响应最大的点，此时具有最好的效果；式(7)为构造的结节相似度函数，w1,w2表示R1和R2的权重系数，实验中设置w1=0.5，w2=15；式(8)中No(λ)i表示第i个像素点的结节相似度值。

提出的ACSLIC算法的优势同样还体现在如下两个方面：

(1)采用近似符合肺结节类圆形形状特征的六边形进行聚类。在像素聚类的过程中不会以为不同方向上距离不同而导致聚类效果有所差异，从而分割得到的超像素样本可以更好保留图像的边缘信息；

(2)精确定位肺结节，直接提取ROI进行操作，加快序列图像分割速度的同时消除了较大部分背景噪声的干扰；算法具体步骤如下：

算法1：改进的超像素序列分割算法(ACSLIC)

输入：肺部ROI图像序列；

输出：超像素分割后的序列图像；

(1)初始化六边形网格中的K个聚类中心Ck=[Lk,ak,bk,xk,yk]T，六边形网格间距S；

(2)执行式(4)～式(8)，构造颜色特征距离和空间位置特征距离间的自适应相似度系数；

(3)在2*S为边长的区域中，以[xk,yk]为开始点采用新构造的相似度度量函数来搜寻与之相似的所有像素点，得到K个超像素样本；

(4)计算所有超像素样本中L、a、b、x、y的均值来更新聚类中心；

(5)重复(3)～(4)步，直到残差小于固定的阈值T;

(6)对分割结果计算所有超像素样本的邻接矩阵，用于后续特征提取；

(7)输入下一张待分割的ROI图像，重复执行(2)～(6)步；

(8)重复(1)～(7)步，直到ROI图像序列分割完毕。

4 超像素样本特征提取

4.1 提取新特征

4.1.1 对比度增强直方图特征

灰度特征是肺结节区域与其余肺部组织最明显的差异。然而肺部CT图像通常因为噪声、伪影等导致成像模糊、对比度低等问题，有必要在提取灰度特征之前对ROI图像进行增强来突出像素间的差别,但是在图像增强的过程中，会导致边缘信息有所丢失。

因此，本文提出一种线性组合canny算子边缘检测图像与直方图均衡化图像的方法，在增强ROI图像对比度的同时保留了较多的图像细节信息。具体步骤如下所示。

(1)计算ROI图像各灰度级出现的概率，计算如式(9)所示

(9)

其中，nk表示ROI图像中灰度级为gk的像素点的数目，L为最大的灰度级数。

(2)通过如式(10)的变换函数，建立gk到Gk的映射关系

(10)

(3)将原ROI图像，标记为I；变换后的灰度值Gk对应的图像，标记为J1；

(4)利用canny算子提取ROI图像的边缘检测图像，标记为J2；

(5)将图像J1和J2进行线性组合，作为最终图像，标记为J；

(6)对图像J进行超像素分割，提取每个超像素样本对应的灰度直方图特征。

本文提出的方法与传统的直方图均衡化的对比如图1所示。

图1 本文方法与传统直方图均衡化的对比

4.1.2 超像素样本邻域纹理特征

对于原始ROI图像，为了更好地辨别超像素样本所处肺部组织的具体类型，比如结节类型、血管类型及背景类型。本文从突出每个超像素样本与其邻域超像素样本的差异角度出发，设计了一组新特征来量化每个超像素样本与其邻接超像素样本的纹理差异。

但是，由于：①血管粘连型结节的灰度值和血管的灰度值非常接近，不能有效将血管和结节分离开；②血管和结节区域的特征尽管很相似，但是在结节内部往往会包含其它构成成分(包含5种情况，分别为钙化，软组织，液体，脂肪和气体)，并且结节在结构上是类圆形，血管更趋向于线性结构。因此，提取有效的纹理特征对区分血管和结节具有很大用处；③但是单个超像素的纹理特征不能有效提取血管和结节的纹理特征，这是因为血管和结节的纹理变化通常来自比超像素更大面积的区域，类似于文献[11]中阐述。

因此，本文提出一种超像素样本邻域纹理特征，具体步骤如下：

图2 高斯金字塔图像与差异图像的计算

(2)由于Level间的尺度不同，因此，先通过双线性插值调整Level2，使其与Level1尺度相同；

(3)计算Level1与Level2的差异图像，标记为S1，如图2(c)(S1)所示；

(4)类似步骤(2)和步骤(3)，计算Level1与Level3、Level2与Level3以及Level2与Level4间的差异图像，分别标记为S2，S3，S4。

(5)调整S1-S4为原ROI图像大小，并对其分别执行改进后的超像素分割算法。对每个超像素样本分别提取灰度均值与灰度方差特征。计算如式(11)和式(12)所示

(11)

(12)

其中，ui(j)为图像Sj中第i个超像素样本的灰度均值；vari(j)为图像Sj中第i个超像素样本的灰度方差；因此，实验中标记单个图像中单个超像素样本对应的特征为SFi=[ui,vari]。

(6)构造当前超像素样本及其邻近左上、左下、右下、右上4个超像素样本的特征集合，标记为ASFi=[SFi,SFi1,SFi2,SFi3,SFi4]，如图3所示。从而构成最终4×2×5=40维的特征向量。

图3 超像素样本的邻域关系

4.1.3 位置信息特征

除辛北断层外，研究区内的断层级别稍低，延伸也较短，在沙三下亚段至沙四段形成了一些地堑和地垒。除X176西的南北向断层外，这些断层大都向下贯穿至沙四下-孔店组，虽然这些断层活动强度均较弱，断距10～50m，对沙四段沉积无明显的控制作用，但由于其切穿沙四下亚段的盐膏层，在断层活动期，仍可对流体起到断裂输导作用。

本文根据提取的ROI图像的中心坐标必处于肺结节区域，从而提出一种新的超像素样本位置信息特征。表示如式(13)和图4所示

(13)

其中，(xc,yc)为ROI图像的中心坐标；(xi,yi)为第i个超像素样本的中心坐标；H和W分别为ROI图像的高度和宽度。

图4 ROI图像中心与超像素样本中心的距离关系

4.2 特征数据归一化

特征数据的归一化是后续聚类的预处理步骤。归一化的主要目的是：有效避免具有更大数值范围特征的主导,并在聚类的过程中显著降低数值计算的困难,因此，本文采用如式(14)所示的方式对特征数据进行归一化

(14)

其中，F是一个原始的特征值；F′是归一化后的特征值；从而，特征向量中的每个值的取值范围都为[0,1]。

5 距离约束稀疏子空间聚类

稀疏子空间聚类(sparse sub-space clustering，SSC)[12]是一种基于谱聚类的子空间聚类方法。给定一组数据，该聚类方法先将所有数据作为字典来线性表示每个数据，接着采用表示系数矩阵的稀疏约束使数据尽可能被同一子空间的数据线性表示，并构造关于子空间系数的相似度矩阵，最后利用谱聚类等方法得到数据聚类的结果。

在本文中，提取的ROI图像是以结节为中心的，并且肺结节图像具有类圆形状特性。基于此先验知识，提出一种基于距离约束的稀疏子空间聚类方法。距离约束的引入有利于结节数据能够更好的被同一子空间内邻近的结节特征数据线性表示，进一步有效排除血管对结节聚类的影响。采用的距离约束函数如式(15)所示

(15)

其中，Di和Dj为第i和第j个超像素样本到ROI图像中心的距离。引入距离约束的稀疏子空间聚类模型，如式(16)所示

(16)

(17)

其中，u为子空间表示系数；U=[u1,u2,…,uK]∈RK×K，为表示系数矩阵；距离越大，1/wji越小，否则1/wji越大；W为相似度矩阵；因此，距离约束有利于结节数据尽可能被同一子空间邻近的最相似的数据线性表示。具体的算法描述如下：

算法2：距离约束稀疏子空间聚类

输入：待聚类图像I，分割后的超像素样本总个数K；

输出：聚类后的肺结节图像；

(1)ACSLIC分割得到K个超像素样本；

(2)对每个超像素样本提取对比度增强直方图特征、邻域纹理特征以及位置信息特征。得到N个特征矩阵；

(3)利用式(15)计算超像素间的距离约束值；

(4)利用式(16)得到子空间表示系数矩阵U；

(5)利用式(17)构造相似度矩阵W；

(6)根据W，采用归一化割法(normalized cut，Ncut)算法得到超像素样本聚类结果。

6 实验结果与分析

本文算法的实验平台是Visual Studio 2013和MATLAB R2014b，PC处理器为Intel(R) Core(M) i7-4770，主频3.40 GHz，内存8 GB。实验数据来源于LIDC-IDRI[13]和山西省人民医院，挑选80组血管粘连型CT序列图像共1468张，单张CT图像大小为512×512。实验中我们邀请到两位影像科医师做手动肺结节分割。同时为了验证本文方法的有效性，针对所有的CT序列图像，对本文方法、脉冲耦合神经网络算法(PCNN-pulse coupled neural network,PCNN)[14]、模糊c均值聚类算法(fuzzy c-means，FCM)[15]以及最新的流向特征熵和测地线距离的粘连血管型肺结节聚类分割方法(flowing entropy and geodesic distance，FEGD)[3]进行了大量的对比实验。

6.1 定性分析

实验中CT序列图像数目较多,本文只对部分分割过程进行说明，如图5所示。

同时，本文与FEGD、PCNN以及FCM算法进行了对比实验。部分结果如图6所示。

图5 本文方法对血管粘连型结节图像序列的分割结果注:(a)原始CT序列图像;(b)分割得到的序列肺实质图像;(c)提取的ROIs图像;(d)ACSLIC对ROIs图像分割的结果;(e)距离约束稀疏子空间聚类得到的肺结节图像掩膜序列;(f)本文方法的分割结果;(g)医师手动分割的结果

图6 4种方法对血管粘连型结节图像分割结果的对比 注：(a)原始CT图像；(b)对(a)局部放大的结果；(c)～(f)分别为PCNN、FCM、FEGD以及本文方法分割得到的肺结节图像掩膜；(g)医师手动分割的结果。

从结果图可以直观观察到，PCNN和FCM对血管粘连型结节的分割结果不准确，不能有效分离血管和结节，导致分割得到的结节粘连有较多的血管组织；此外，FEGD的分割效果较好，结果优于PCNN和FCM，但是分割结果存在少量的边界泄露(如图6，行1-4，行6所示)；并且FEGD在特定的CT图像中会不可避免的引入部分血管(如图6，行5所示)，这是因为FEGD不能有效处理血管横断面为圆形的血管；然而，本文方法可以有效的分离血管和结节，并且改进的超像素分割算法保留了更多的边界信息，分割结果更平滑。相比这3种典型分割算法，本文方法都能取得最好的分割效果。并且与医师手动分割的结果基本一致。

6.2 定量比较

本文通过采用变化信息、概率边缘指数以及时间复杂度来对分割方法进行定量比较。

6.2.1 变化信息

变化信息(variation of information，VoI)是计算图像像素点从一个聚类转变为另外一个聚类丢失或获得信息量的大小[16]。假设参考分割图像与实际分割结果分别用Sa和Sb表示。则可利用参考分割图像Sa的熵[17]H(Sa)、实际分割结果Sb的熵H(Sb)以及Sa和Sb的联合熵I(Sa,Sb)，得VoI计算如式(18)所示

VoI(Sa,Sb)=H(Sa)+H(Sb)-2I(Sa,Sb)

(18)

(19)

(20)

(21)

其中，P(m)表示一个像素被分割到第m个区域中的概率。VoI越小，表示实际分割结果与参考分割图像越接近，分割效果越好。

6.2.2 概率边缘指数

概率边缘指数(probabilistic rand index，PRI) 检验实际分割结果与参考结果之间的属性共生的一致性的参量[18]。假设(xi,yi)是原图S中的一个像素对，如果在Sa中标记为(ai,aj)且相同，则在Sb中标记为(bi,bj)也应该相同。PRI的计算如式(22)所示

(22)

其中，M为原始图像S中的像素个数；I为判别函数，主要作用是判断像素对是否具有相同标记。PRI 值越大，说明实际分割结果与医师手动分割结果的属性共生一致性越好，表明分割结果越好。

此外，为了验证本文方法的通用性，实验中计算了所有序列CT数据分割结果的VoI与PRI的均值。见表1。

表1 4种算法对所有血管粘连型结节图像序列的VoI与RPI均值

实验结果表明，对血管粘连型结节的CT序列图像，PCNN和FCM两种算法的分割效果基本一致；FEGD的分割效果远优于PCNN和FCM；但是，相比PCNN、FCM和FEGD这3种分割算法，本文方法都达到最低的VoI均值和最高的PRI均值，从而，从客观的定量比较结果进一步验证了本文方法的在分割这类型结节图像上的优越性。

6.2.3 时间复杂度

通过对所有的实验数据进行分析，发现一个CT序列图像包含结节的平均数量为8张。因此，表2显示了4种算法序列分割的平均处理时间。虽然对于单个CT序列，平均处理时间比FCM算法多出0.48 s，但是本文方法的分割效果远优于FCM算法的分割结果，在医学上这种用时间来换取分割的精度是值得的。

表2 4种算法分割序列CT图像的平均处理时间

7 结束语

针对血管粘连型结节的序列CT图像，本文提出一种基于超像素和稀疏子空间聚类的肺结节图像序列分割方法，能够实现对血管粘连这种困难分割结节图像的准确分割。该方法的主要贡献是：在超像素分割算法中构造了自适应相似度系数；在超像素样本特征提取中，设计提取了对比度增强直方图特征、超像素样本邻域纹理特征以及位置信息特征等多个新特征；在稀疏子空间聚类算法中引入距离约束项；这3个方面都有效分离了血管和结节。从定性分析和定量比较结果可以看出，本文能更好地分离血管和结节，边缘轮廓平滑，具有最好的序列分割结果，且与医师手动分割结果基本一致。因此，本文的序列分割方法能更完整、准确地分割出血管粘连型结节图像序列。