基于稀疏表示和随机游走的磨玻璃型肺结节分割

李祥霞 李彬 田联房 张莉 朱文博

肺癌已经成为最高致死率的恶性肿瘤疾病之一.临床数据表明,早发现、早治疗能够有效地提高生存率和降低患者死亡率[1].肺结节是肺癌的早期影像学表现形式.目前,电子计算机断层扫描(Computed tomography,CT)是最广泛应用于肺癌的影像学诊断方法[2].传统的分割方法是医生获取病人的某一断层数据,进行手动地分割肺结节.然而,这不仅费时费力,而且往往容易受到主观因素的影响.因此,肺结节的自动分割已经成为研究的热点.对于实性肺结节的分割,国内外均有大量的文献报道.但是,对于磨玻璃型(Ground glass opacity,GGO)结节分割文献仍报道甚少.不幸地,与实性结节相比,GGO结节更具有发展为恶性肿瘤的潜力[3].因此,GGO结节的分割对于早期肺癌的筛选和诊断更具有临床应用价值.然而,GGO肺结节的分割一直是一个挑战性的任务,这是因为GGO结节具有形状不规则、边界模糊和结节内部灰度分布不均匀等特性.

随机游走(Random walker,RW)是将图像分割看成网格图上标签安排问题.首先需要用户在图像上标记合适数量的目标和背景种子点,然后计算每个未标记顶点首次到达给定类种子点的概率值.随后,将未标记顶点归类到最大概率值所对应的类别中,从而实现最终的分割.RW只需求解一个稀疏线性方程组,因此其无需迭代,操作简单快速.RW在处理弱边缘时能够取得较好的结果.然而,大部分RW算法采用了固定的邻域结构来构建网格图.然而,由于图像噪声、伪影等诸多因素的影响,这种邻域结构可能捕捉了不恰当的邻域关系,从而降低了分割的准确性.此外,RW仅仅采用了灰度信息构建加权矩阵,这对于对比度低、灰度不均匀的GGO结节较难直接分割.

为了解决这些问题,本文提出了基于稀疏表示和随机游走(SRRW)的算法对GGO肺结节进行自动分割.首先,根据图像的信息自动地选取种子点,避免人为的干预和减少操作的时间;其次,本文构造了一个新的图来替代传统的无向加权图,利用8-邻域和稀疏表示的K-最近邻(K-nearest neighbor,K-NN)算法解决了加权图对噪声的敏感性问题.最后,引入标签限制项到RW 的能量函数中,惩罚了预测标签与给定标签的不一致性,从而提高了分割的准确性.本文采用了美国肺部影像数据库LIDC来验证提出算法的准确性和有效性.实验结果表明,提出的算法极大地提高了GGO肺结节分割的有效性和准确性.

1 相关工作

在过去几十年里,众多学者已经提出了大量的肺结节分割方法[4],主要可分为两类:基于灰度的分割方法和基于形状的分割方法.基于灰度的分割方法是利用图像的灰度信息实现对肺结节的分割,例如阈值法[5],区域增长[6−7],聚类[8].这类分割算法对灰度均匀的实性结节的分割效果较好,但是对于胸膜粘连型、血管粘连型和GGO结节较难分割.基于形状的分割方法通过整合了结节的几何和形态信息来提高结节分割的准确性,例如形态学[9−10]、局部形状分析[11−12]、活动轮廓[13−15].目前,已经有一些学者对GGO结节分割进行了研究,并取得了一定的成果.例如,Zhu等[16]提出了基于马尔科夫随机场的算法进行GGO结节的分割.Li等[17]提出了改进的模糊聚类算法进行GGO结节的分割.Li等[18−19]提出了基于自适应局部区域能量函数的活动轮廓进行GGO结节的分割.近年来,一些学者已经将深度学习用于肺结节的分割.例如,Cheng等[20]提出了基于深度学习的肺结节计算机辅助诊断系统.Wang等[21]提出了中心聚焦卷积神经网络对肺结节进行分割.

基于图论的随机游走算法是近年来近年来逐渐兴起的交互式图像分割算法,它首先被Grady等[22]引入图像分割的领域.Eslami等[23]提出了导向随机游走的方法进行左心室的分割.Onoma等[24]利用相邻像素点之间的欧氏距离来决定参数的取值,并将标签的概率密度整合到RW能量函数中,从而提高了病变组织分割的准确性.Wang等[25]利用拓扑细化方法自动地选择种子点.Mi等[26]将肿瘤的增长预测信息整合到RW能量函数进行肿瘤的自动分割,从而提高了分割的性能.

稀疏表示(Sparse representation,SR)已经被引入到图构建中.例如,Fan等[27]提出了稀疏正则化半监督学习算法.Liu等[28]提出了低秩表示(Low-rank representation,LRR)进行子空间的聚类.Li等[29]提出了一种稀疏平衡图构建方法.Fu等[30]提出了张量低秩表示和基于稀疏编码的子空间聚类算法.Zhuang等[31]提出了一种新的半监督图学习方法进行数据的半监督分类.Liao等[32]提出了图正则化自动编码器进行数据的聚类.这些方法都是基于稀疏表示来构建稀疏图,并将构造的图应用到谱聚类、子空间表示、子空间聚类等.

2 基于稀疏表示和随机游走的GGO肺结节自动分割

本文提出了基于稀疏表示和随机游走(SRRW)的GGO肺结节分割的算法.算法的流程图如图1所示.

图1 SRRW算法的流程图Fig.1 The flowchart of the SRRW algorithm

3 基于稀疏表示的加权矩阵的构建

在随机游走算法中,图的构建是一个至关重要的环节.大部分RW算法采用了4-邻域或者8-邻域搜索策略来构建图.然而,以这种方式构建的图对噪声比较敏感.为了解决这个问题,本文将8-邻域和基于稀疏表示的K-NN方法结合起来,构建了一个可靠的新图.新的图能够更好地描述顶点的邻域关系,提高了对噪声的鲁棒性.它也能够更好地传递种子点的标签信息到未标记顶点,提高了分割的准确性.

3.1 基于8-邻域的加权矩阵构建

给定一个包含n=W×H个像素点的图像,其中W和H分别是图像的宽度和高度.然后,将输入图像表示为一个无向加权图G=(V,E,W),其中V是图的顶点集,E是图的边集.每个图像像素点对应着图的一个顶点i∈V,每条边eij∈E⊆V×V连接一对相邻顶点i和j.每一条边eij∈E被赋予一个非负权值wij∈W,衡量了一对相邻顶点之间的特征差异性.本文利用了Gabor函数将图像变成频域上的幅度谱,进行纹理特征的提取.图2是一个GGO结节CT图像中感兴趣区域(Region of interest,ROI)纹理的分析.图2(a)和2(b)分别是原始的CT图像和纹理图像.图2(d)和2(e)是在结节和背景的ROIs(如图2(c)所示)中像素的纹理值.正如图2(d)和2(e)所示,结节和背景像素的纹理值存在着明显的差异,因此可以将纹理作为特征来区分结节和背景.

图2 一个GGO结节CT图像纹理的分析Fig.2 Texture analysis in a CT image with GGO nodule

纹理特征提取后,对于每一个图像像素点,本文整合了灰度、纹理和空间位置坐标构造了一个特征向量,如式(1)所示.

其中,Imin和Imax分别是像素灰度值的最大值和最小值.Tmin和Tmax分别是像素纹理值的最大值和最小值.归一化的特征向量的表达式如式(4)所示.

本文采用了高斯核函数来计算一对8-邻域顶点i和j之间的加权值,表达式如式(5)所示.

3.2 基于稀疏表示的K-NN的加权矩阵构建

RW算法的分割性能在一定程度上依赖于构造的图.传统的邻域顶点选择方法主要包括两种:基于4-或者8-邻域,基于K-最近邻.以这些方式选择的邻域顶点依赖于邻域大小参数的选择,然而,最优参数的选择并不是一件容易的事.此外,这些方法是通过计算了两个顶点之间的欧氏距离来决定邻域顶点,因此对噪声比较敏感.为了解决这些问题,本文将稀疏表示(Sparse representation,SR)引入到图构建中.与传统方法构建的图相比,基于SR方法构建的图能够自适应地捕捉图的邻域结构,因此对噪声具有很强的鲁棒性.不同于Yuan等[33]采用K-NN来选择与目前顶点欧氏距离最小的K个顶点,本文采用了SR方法来选择K个邻域顶点.首先,利用稀疏表示的能量函数计算了每一个顶点的最佳稀疏重构系数,然后,利用K-NN算法将数据点的搜索限制在空间的局部邻域内,极大地减少了计算的复杂性.最后,用得到的稀疏矩阵定义了一个新的加权矩阵.

式(7)对应的增广拉格朗日方程表示如式(8)所示.

其中,Y1和Y2是拉格朗日乘子,u＞0是惩罚参数.第三项是ADMM算法加入的Frobenius范数惩罚项.更新迭代变量C,Y1,Y2,u,直到满足收敛条件,变量C的迭代公式如式(9)所示.

最终得到最优的稀疏矩阵C∗.然后,利用KNN将当前顶点的邻域搜索限制在局部邻域内,即选择K个顶点来表示当前顶点,从而得到最终的稀疏矩阵C∗.因为稀疏矩阵C∗测量了成对顶点的相似性,因此用其构造了一个新加权矩阵.加权值计算公式如式(10)所示.

其中,ESR表示当前顶点与稀疏表示的K-NN选取的邻域顶点彼此相连接组成的边集合.图3是构建的新的加权图每一个实心圆圈表示一个顶点,各实心圆圈分别表示当前的顶点,当前顶点的8-邻域顶点和基于稀疏的K-NN选择的邻域顶点.

图3 新加权图Fig.3 A new weighted graph

4 种子点的自动获取和SRRW能量函数的定义

众所周知,RW对种子点的位置和数量非常敏感.传统RW需要用户在感兴趣的区域手动地选取种子点.然而,人工选取种子点是费时和冗余的.本文提出了自动获取种子点的策略,包括两个主要的步骤:初始的肺结节分割,基于测地距离和局部搜索策略的种子点的选取.首先,自适应阈值和形态学开运算等方法进行粗略地分割肺结节,分割结果被表示为R.然后,利用测地距离(Geodesic distance,GD)[10]获取了R的一个小区域R0.本文通过设定一个合适的阈值T,计算从R的边界到结节中心的最大测地距离D.最后,将区间[D−T,D]范围内的像素点作为肺结节的种子点,如图4(b)所示.由于GGO结节存在形状不规则和边缘模糊等特性,基于测地距离产生的结节种子点易陷入均匀的区域.为了解决这个问题,一个局部搜索策略被提出,来选取更多的结节种子点.本文利用了灰度和Gabor纹理特征来定义相似性函数,指数函数被引入来强调灰度值的重要性,如式(11)所示.

其中,i是区域R的边缘上的像素点,Ni是像素点i的8-邻域像素集.F0表示初始的肺结节种子点集.如果相似性函数满足下面准则式(12),那么F0被迭代更新,即像素j被纳入到肺结节种子集F中.

其中,κ是预定义的阈值,凭经验设定为10.当全部的边界像素点被历遍,搜索停止.图4(c)展示了全部肺结节种子点F.本文利用提出的局部搜索策略来获取背景种子点.B0表示初始的背景种子点集.如果相似性函数满足下面准则式(13).

其中Rb表示从肺结节的中心到区域R的预先定义距离Tb(Tb＞M)的一个区域.当没有新的像素点被纳入到背景种子点集B时,迭代停止.图4(c)显示了全部的种子点S=F∪B.结节内部和外部标记点分别表示肺结节种子点和背景种子点.考虑到在连续断层CT图像上肺结节位置特征以及当前切片与相邻的上下切片对应的特征的相似性和相关性,本文从第k张切片分割的肺结节中找到肺结节的中心位置,然后将对应的肺结节和背景的种子点反馈到第k+1张CT切片的分割工作中,进行下一张连续CT切片的分割过程中.本文提出的算法在分割多帧连续的断层图像时无需反复选取种子点,大大减少了计算复杂性.

图4 种子点的自动获取Fig.4 Automatic acquisition of seeds

在种子点自动获取后,将获取的种子点注入到提出的SRRW算法中.本文引入标签限制项到原始RW能量函数中来提高分割的准确性.标签限制项迫使预测的种子点概率与给定的种子点概率的一致性.改进的能量函数如式(14)所示.

其中,fi和fj代表的是顶点i和对应的邻域顶点j的概率.λ是一个参数.|S|表示种子点集S中元素的个数.bi表示顶点i的指示函数,如式(15)所示.

式(14)的第一项和第二项是圆滑限制项,分别鼓励了在E8和ESR中邻域顶点之间概率的一致性.第三项是标签限制项,惩罚了在最优能量框架下预测的种子点概率与给定的种子点概率不一致性.式(14)可以表示为式(16).

将式(16)转化为矩阵形式,如式(17)所示.

5 实验结果和分析

本文采用数据库LIDC[34]验证提出算法的准确性和有效性.该数据库共收录了1018个病例.对于每一个病例的XML文件,记录了4位经验丰富的胸部放射科医生对直径大于3毫米的肺结节轮廓点的坐标.本文将4位医生勾画结果经协商达成一致所对应的结节轮廓作为金标准.本文从LIDC中随机地选取100个带有不同形状、大小和纹理的GGO结节图像.全部图像像素点的灰度值范围被归一化在0～255之间.实验硬件环境为CPU Intel(R)Xeon(R)E3-1225 v5,主频3.30GHz;内存4.00GB;显卡Intel(R)Ethernet Connection(2)1219-LM.实验软件环境为MATLAB R2016b开发平台.

5.1 SRRW算法的参数设置

为了设置合理的参数值,本文通过对参数设置不同参数值进行了对比实验.在种子点自动获取过程中,主要有四个参数:阈值T,κ,γ和Tb.参数T决定了在小区域R0像素的个数.参数κ决定了在区域R的边缘上像素点的个数.图5显示了在参数T和κ选取不同值时,结节种子点的选取情况.图5(a)和图5(b)分别是在T=1,κ=5和κ=10时结节种子点个数.图5(c)和图5(d)分别是在T=4,κ=5和κ=10时结节种子点情况.图5(e)～(h)是SRRW算法所对应的分割结果.从图5(e)中可以看出,当结节种子点远远小于背景种子点时,SRRW 算法会出现严重的过分割现象.但是,当增加κ值时,过分割现象会得到一定程度上缓解.因此,参数κ设置为10时,能够产生好的分割结果.在本文中,参数T设置范围为3毫米到35毫米,覆盖所有大小结节.在背景种子点选取过程中,参数Tb设定范围很广泛,在本文所有实验中设定为100.类似于此对比实验方法,参数γ经多次试验设定为5.

图5 T和κ取不同值时结节种子点选取和对应的分割结果Fig.5 Seed selection and the corresponding segmentation results when using differentTandκ

在能量函数定义中,式(14)中参数λ控制标签限制项的作用.当λ设定为0时,SRRW算法的标签限制项不起作用.图6显示了不同λ值时SRRW算法获取的分割结果.图6(b)～(d)分别是λ=10,λ=100和λ=1000时SRRW 算法分割结果.从图6(b)结果可以看出,结节分割可能会出现过分割现象,这是因为标签限制项起到很小的作用.随着λ的增大,标签限制项的作用会越来越大.当λ=100时可以得到更准确的分割结果,如图6(c)所示.但是,当λ取值过大时,结节分割会出现严重的分割错误,如图6(d)所示.因此,经多次实验,本文设定λ=100，执行了所有实验.

图6 不同λ时本文提出的SRRW算法分割的结果Fig.6 Segmentation results of the proposed SRRW algorithm with differentλ

5.2 定性SRRW算法分割结果的分析与对比

为了验证本文提出的局部搜索策略获取的结节种子点的有效性,本文在选取不同的结节种子点时,利用提出的SRRW算法对GGO结节进行分割.结节种子点获取方法分别是测地距离方法和测地距离、局部搜索策略联合.图7是一个典型GGO肺结节分割的过程.图7(a)是原始CT图像,图7(b)和7(c)分别显示了局部搜索策略获取的背景种子点和测地距离获取的结节种子点,图7(d)显示了测地距离、局部搜索策略联合所获取的全部结节种子点.图7(e)和7(f)分别是利用图7(c)和7(d)结节种子点时SRRW算法的分割结果.从图7(c)可以观察到,测地距离获取的结节种子点集中于GGO肺结节灰度均匀区域,导致了分割结果不理想(图7(e)所示).相反,利用图7(d)的结节种子点,SRRW算法能够获取更准确的分割结果,如图7(f)所示.实验结果表明,本文提出的局部搜索策略获取的结节种子点能够提高SRRW算法的准确性.

图7 不同的结节种子点时SRRW算法的GGO结节分割Fig.7 GGO nodule segmentation of SRRW algorithm using different nodule seeds

为了验证SRRW 算法分割的准确性,本文采用了传统 RW[22]、RWR[35]、SubRW 算法[36]与SRRW 算法进行了一个对比实验.为了公平起见,所有算法均采用了本文获取的种子点.图8(a)～(d)分别是RW、RWR、SubRW 和SRRW算法所分割结果的放大图.从图8(a)中可以观察到,RW对GGO结节进行分割时产生边界泄露问题.这是因为RW 仅采用了4-或8-邻域的搜索策略进行图构建,且仅利用灰度信息描述相邻顶点之间的相似性,这对结节内部灰度分布不均匀且变化较大的GGO结节分割来说是无效的.RWR算法获得了更好的分割结果,如图8(b)所示.这是因为RWR从一个顶点出发,以一定的概率转移到其邻域点,并且反复迭代,这使得算法能够捕捉图的整体结构信息,提高了分割的准确性.由于SubRW引入了标签先验知识,其获得了更好的分割结果,如图8(c)所示.然而,SubRW仍然未解决噪声敏感性问题.本文提出的SRRW算法的分割结果更接近于结节边界,如图8(d)所示.这是因为提出的算法联合了8-邻域和稀疏表示的K-最近邻算法,获得了更可靠的顶点邻域关系,捕捉了更多的结构信息,从而提高了肺结节分割的准确性.

图8 RW、RWR、SubRW和SRRW对GGO结节分割结果Fig.8 Segmentation results of GGO nodule by using RW、RWR、SubRW and SRRW

为了进一步验证SRRW算法的准确性,本文随机地选择4个具有不同大小和形状的GGO结节图像.图9显示了肺结节的分割结果.图9(a)是原始CT图像,图9(b)是SRRW 分割结果的放大图和金标准.从图9(b)中可以观察到,SRRW算法的分割结果视觉上几乎接近于金标准.实验结果表明了本文提出的SRRW算法能够准确分割各式各样的GGO肺结节.

图9 SRRW算法分割结果Fig.9 Segmentation results of the SRRW algorithm

为了验证SRRW算法对一个CT图像序列的分割性能,本文利用了提出的SRRW 算法对含有GGO肺结节的所有CT切片进行分割.首先,利用SRRW算法对单张肺结节切片进行分割.然后,考虑到在连续断层CT图像上肺结节位置特征以及当前切片与相邻的上下切片对应的特征的相似性,本文从单片分割的肺结节寻找肺结节的中心位置,然后选取相应的结节和背景种子点进行连续切片的分割.图10显示了提出的SRRW算法对一个CT序列图像的GGO结节的分割结果.图10(a)和10(b)分别是带有GGO肺结节的CT序列图像和对应放大的结节区域.图10(c)是提出的SRRW算法的分割结果和金标准.从图10(c)我们可以观察到,提出的算法能够准确地分割一个GGO肺结节序列图像.

图10 本文提出SRRW算法对GGO结节序列图像分割Fig.10 Segmentation of the proposed SRRW algorithm in a CT image sequence with a GGO nodule

5.3 定量SRRW算法分割结果的分析与对比

为了定量地验证本文提出SRRW算法的准确性和有效性,本文采用了Overlap作为分割的评价指标.它测量了算法的分割结果和金标准之间的重叠面积比,计算公式如式(20)所示.

其中,SA和SGT分别是算法分割结果和金标准.|SA∩SGT|是既包含SA又包含SGT的像素个数.|SA∪SGT|表示包含SA或SGT或包含两者的像素个数.在理想情况下,一个分割结果越准确,Overlap值越接近1.

首先,本文从LIDC数据库中随机地选取100个带有GGO肺结节的CT图像进行一个对比实验.为了验证提出算法能够分割不同样式的GGO肺结节,选取的100个图像尽可能包含大小各异、形状各异、纹理不同及灰度分布不均的GGO结节.本文在100个数据样本结果中随机地选取10个数据的Overlap值列举在表1.从表1中可以看出,与传统RW 算法[22]、SubRW 算法[36],本文提出的SRRW算法获得了较高Overlap平均值和较低方差值.实验结果表明,提出的SRRW算法具有较高的分割准确性和较强的鲁棒性.

为了统计提出的算法在种子点选取过程中的执行时间,本文从100个样本分割结果中随机选取10个样本的种子点选取的执行时间列举在表2.从列表2中可以看出,提出的种子点的选取在10个数据上执行时间平均为0.2136秒.实验结果表明了提出的自动种子点选取策略的有效性.

表1 SRRW、RW和SubRW算法的Overlap值对比Table 1 Comparison results of Overlap values among SRRW、RW and SubRW algorithms

表2 本文的种子点选择策略的执行时间Table 2 Execution times of seed selection strategy

5.4 定量SRRW算法与其他分割算法的对比

为了进一步定量地验证SRRW算法的有效性,本文采用了Overlap作为分割精度的评价指标,用本文提出的算法与其他分割算法进行了一个对比.表3列出了不同分割算法的Overlap平均值.Okada K等[37]获得了相对较低的Overlap平均值(0.45),这是因为该算法将肺结节模拟为椭圆形,不能很好地处理形状不规则的GGO肺结节.Kostis等[38]也获取了较低的Overlap平均值(0.57).Kuhnigk等[9]采用了形态学开运算方法分割形状规则的小肺结节,也不能很好地处理GGO肺结节.Kubota等[7]提出了扩散竞争和区域增长混合算法,获得的Overlap平均值为0.66,相比于Kuhnigk等提出算法的分割结果具有较大提高.最近,Messay等[39]提出了回归神经网络算法,获取了较高的分割性能.该算法获得的Overlap平均值为0.77.然而,算法需要人工提供控制点,达到提高分割精度的目的.本文提出的SRRW算法获取Overlap值为0.91,得到了较为满意的Overlap值.实验结果表明了SRRW算法对GGO肺结节分割的有效性.

表3 不同肺结节分割算法的Overlap值的对比Table 3 Comparisons of Overlap values among different pulmonary nodule segmentation algorithms

6 结论

本文提出了基于稀疏表示和随机游走(SRRW)算法对磨玻璃型肺结节进行分割.本文提出了自动选取结节和背景种子点方法,减少了人工干预.采用基于8-邻域和稀疏表示的K-NN方法构建了一个新加权图,避免了CT图像噪声和伪影干扰.通过整合灰度、纹理和空间距离和稀疏系数构建了一个新加权矩阵,更有效地衡量了相邻顶点之间的相似性关系.最后,引入标签限制项到RW能量函数中,提高了对GGO肺结节分割的准确性.通过定量和定性对比实验可以看出,本文提出的算法对GGO肺结节分割的准确性和有效性.下一步研究主要集中于进一步优化提出的模型,并将该算法扩展到对任意类型肺结节的分割.