骨边界增强滤波的图割算法

石志良,范伟楠,甘梓博,袁 琼

(1.武汉理工大学机电工程学院,湖北 武汉 430070)(2.湖北第二师范学院计算机学院,湖北 武汉 430205)

医学图像具有灰度广、对比度低等特点,快速精准地分割骨组织一直是难题. 临床通常采用阈值分割和区域生长相结合的方法进行手工分割,不仅需要大量人工交互,而且需要操作人员具备深厚的专业知识,分割结果易受主观因素影响,具有不可重复性. 因此,研究一种快速、稳定的自动分割算法具有重要意义.

常见的骨自动分割算法分为机器学习和非机器学习两类. 机器学习方法,可实现令人满意的精度和自动化,但技术人员需要对特定器官进行大量标注和反复训练,过程复杂、耗时,且仅适用于单个器官. 非机器学习方法指借助图像中固有的灰度、梯度、纹理等信息,根据一定的相似性准则,将图像分割成各具特性的区域. 主要方法包括阈值分割[1]、区域生长[2-3]、水平集[4-5]、聚类[6-8]、图割[9-20]等,这类算法具有通用性的优点,适用于更广泛的分割问题,因此研究一种非机器学习的自动分割算法具有不可替代的意义.

近年来,图割作为一种非机器学习的方法,具有稳定、全局最优、可扩展性好等优点,在图像分割领域占据着重要的地位. 其中最为经典的是Boykov等[10-11]提出的基于图像强度的交互式图割算法. 然而,由于算法仅考虑图像强度,在分割骨组织时鲁棒性较差,并且采用交互式的方法增加图割约束,效率低,无法满足自动分割的需求. 为提高算法的鲁棒性和效率,许多研究人员对图割算法进行改进,先后提出了基于K聚类的图割[9]、基于双通道的图割[12]、基于形状约束的图割[13]、基于分水岭技术的交互式图割[14]、基于显著性的图割[15]、基于特征检测的图割[16]、基于图像数据映射的图割[17]、结合U-net网络的图割[18]、金字塔图割[19]、基于区域生长的图割[20]等相关算法.

虽然上述算法在一定程度上提高了图割的鲁棒性,但未能彻底解决算法的交互问题,尤其在分割对比度低的骨边界时效果较差. 因此,研究提出一个基于Hessian矩阵的骨骼边界增强滤波器,在此基础上建立新的能量函数,求解能量函数得到初步分割结果,并对初步结果进行后处理,实现骨组织的自动分割.

1 传统图割算法

图割算法的基本思想是将N维图像的分割问题转换为求解最小能量函数问题. 能量函数的解与图模型中的割线唯一对应,并将图分为目标和背景两部分. Boykov等[10]提出的能量函数为

E(A)=λ·R(A)+B(A).

(1)

式中,系数表示R(A)相对于B(A)的重要程度,A=(A1,…,Ap,…A|P|)是一组二进制向量,P是图中所有像素的集合,Ap对应P中像素p,为目标点(表示为1)或背景点(表示为0),当所有Ap已知时,向量A即为图的一种分割方式.

R(A)为区域项,定义为

(2)

式中,RP(AP)表示将像素标记为目标或背景的惩罚,分别用RP(“obj”)和RP(“bkg”)表示.

B(A)为边界项,定义为

(3)

式中,当Ap=Aq时,δ(Ap,Aq)=1,否则δ(Ap,Aq)=0;B{p,q}表示相邻像素p,q的一致性惩罚,这种一致性可以基于图像强度、梯度或其他标准.当像素p,q相似(如在图像强度上相近)时B{p,q}很大,反之接近于0.N表示由像素p及其相邻像素点构成的邻域系统.

为获得全局最优的结果,用户需要标记目标种子点和背景种子点以增加图模型的输入.通常,将区域项和边界项称之为软约束,种子点称之为硬约束.

图1为传统图割算法示意图,图1(a)为灰度图像,建立的加权无向图G=(V,E),如图1(b)所示,V为图中节点,表示为

图1 图割算法示意图

V=P∪{S,T},

(4)

式中,P为像素点集合;S,T为终端节点,分别称之为源、汇.

E为图中边,表示为

E={p,q}∪ {p,S}∪{p,T},

(5)

式中,p,q为像素点;{p,q} 表示相邻像素构成的边,称为n-links;{p,S} 表示像素点与源构成的边;{p,T} 表示像素点与汇构成的边.{p,S} ,{p,T} 统称为t-links.

边的权由软约束和硬约束共同决定,{p,q} 的权为B{p,q};{p,S} 的权分为3种情况,当p为背景种子点时其值为0,p为目标种子点时其值为1+max(B{p,q}),p为非种子像素点时其值为λ·Rp(“bkg”);类似的{p,S} 的权为,当p为背景种子点时其值为1+max(B{p,q}),p为目标种子点时其值为0,p为非种子像素点时其值为λ·Rp(“obj”).

根据上述原理,用最大流算法求解能量函数,使得t-links和n-links权值最小. 如图1(c)所示,曲线1即为所求割线,最终分割结果如图1(d)所示.

2 本文方法

2.1 骨边界增强滤波的图割模型

算法总体流程如图2所示,主要由骨边界增强滤波、能量函数建立和后处理3个部分组成.

图2 本文算法流程图

骨边界增强滤波输出结果的好坏直接决定骨边界的清晰度,进而影响图像分割质量,因此,滤波器的设计是本文第一个关键步骤. 如何将滤波结果融合到图割约束项中,建立包含图像强度与结构特征的能量函数,是本文第二个关键步骤. 此外,由于CT图像固有的部分体积效应,分割后的骨组织可能存在粘连的情况,如何设计一种后处理方法将相邻骨分离,提高分割的鲁棒性和效率,是本文第三个关键步骤. 针对这3个关键步骤对图割算法进行改进,以期提高图像的分割精度与自动化程度.

2.2 骨边界增强滤波器

骨分割难点之一是关节处骨边界对比度低,因此仅基于图像强度特征的分割方法,不能取得理想的分割效果. 为提高骨边界对比度,使用经典的边缘增强滤波器(非锐化遮蔽滤波)锐化图像,以突出CT图像中灰度突变的区域. 图3表示股骨切片图像非锐化遮蔽处理后的效果,由图可知,骨边界强度较原图得到提升,但相邻骨之间仍存在连接,对分割结果产生干扰.

图3 非锐化遮蔽效果图

为此,基于Hessian矩阵设计一种增强骨边界的滤波器,称为骨边界增强滤波.

Hessian矩阵是由多元函数的二阶偏导数构成方阵,主要用来描述函数局部曲率. 以二维图像为例,令表示图像像素关于坐标的函数,将其在泰勒展开得

(6)

则其Hessian矩阵为

(7)

在二维图像中,Hessian矩阵是二维正定矩阵,具有2个特征值和2个特征向量,其中特征值表示图像在2个特征向量所指方向上的各向异性.

利用特征向量和特征值(λ1,λ2)绘制椭圆,如图4(a)所示.当λ1=λ2时,各向同性最强,对应图4(b)中斑点结构;当λ1=1,λ2=0或λ1=0,λ2=1,各向异性最强,对应图4(c)中线性结构. Hessian矩阵特征值与图像特征对应关系见表1.

表1 二维Hessian矩阵特征值与图像特征对应关系

图4 二维图像结构特征

二维图像求二阶导的一般方法为

fxx(x,y)=f(x,y)-f(x-Δx,y)-(f(x-Δx,y)-f(x+2Δx,y)),

(8)

显然此方法只包含自身在内的3个像素信息,易受图像局部信号干扰.根据线性尺度空间理论,对一个函数求导等于该函数与高斯函数导数的卷积,其表达如下

(9)

式中,σ表示高斯函数的标准差,称为尺度,fxx(x,y;σ)表示在图像尺度σ下的二阶偏导数.高斯模板包含周围矩形范围内所有像素信息,因此能显著降低求导误差.

由于在一幅图像中,同种结构具有不同大小,如粗细不同的线,故同一尺度无法适用于整张图像.为此,采用多尺度方法,即对同一点使用不同尺度的高斯模板进行卷积,选择各向异性最强的结果作为该点的输出.

将上述理论推导至三维图像,得尺度下的三维Hessian矩阵为

(10)

式中,表示三维图像在体素v的二阶偏导数.设矩阵3个特征值为λ1,λ2,λ3(|λ3|≥|λ2|≥|λ1|).以3个特征值为半轴可得特征椭球如图5所示.

图5 特征椭球

显然,当图像各向同性强时椭球将被压缩为球体;当有两个特征值的绝对值接近1,另一个接近0时,椭球将被压缩为柱体;当有两个特征值的绝对值接近0,另一个接近1时椭球将被压缩为片状. 因此可将三维图像的特征结构分为片状、管状和斑点状3种,如图6所示,与特征值的对应关系见表2.

表2 三维Hessian矩阵特征值与图像特征对应关系

图6 理想结构

其中噪声状是斑点状的特例,其特点是特征值均为0;骨边界是由骨外层的皮质骨组成,由表可得皮质骨对应的结构为片状,因此使用Hessian矩阵提取骨边界并对其进行增强是可行的.

为量化CT图像体素v局部区域与片状结构的相似性,可定义图像在尺度下的评估函数

(11)

式中,α,β,γ为常数,R1,R2,R3表示与片状、管状、斑点状、噪声状有关的项,称之为测量因子,分别定义为

(12)

根据定义计算出不同结构对应的测量因子R1,R2,R3的值,见表3.

表3 图像结构与测量因子的对应关系

如上文所述,同一尺度无法适用于整张图像,为提高评分函数的鲁棒性,计算不同尺度下的评分最大绝对值,并将最大特征值纳入函数中,得

(13)

式中,∑为不同尺度构成的集合.

图7为股骨CT切片经骨边界增强滤波器处理效果图,由图7(b)可知,处理后相连的股骨与髋骨分离开来,骨边界清晰度度提高,证明本文提出的基于Hessian矩阵的骨边界增强滤波器是有效的.

图7 滤波处理效果图

2.3 约束项计算

2.3.1 改进区域项函数

传统的区域项函数仅基于图像强度,由于骨边界薄弱、间隙狭窄、骨小梁处强度低,不能取得很好的分割结果.为此考虑结合图像强度与评分函数,定义骨区域和背景区域.

在CT图像中图像的强度经过线性变换可以得到CT值,其单位为Hounsfield Unit(HU),CT值与图像组织的对应关系,见表4.

表4 不同组织对应HU值

由表4可得,骨组织的CT值一般在400HU以上,软组织的CT值一般不超过150HU,空气和脂肪的CT值一般小于-80HU,结合评分函数,定义估算的骨组织区域Ebone和背景区域Ebkg分别为

Ebone={x∈Ω|I(x)≥400且S(x)>0},Ebkg=lcc({x∈Ω|I(x)<-80}),

(14)

式中,表示体素的CT值,表示高强度的皮质骨,表示为排除骨间通道中的软组织体素,lcc表示输入二进制参数的最大连通分量,{x∈Ω|I(x)<-80}表示脂肪和空气. 根据式(14)对股骨切片图像进行处理,结果如图8所示,其中图8(b)为非骨部分的像素区域记为Ebkg,图8(c)为非骨部分的像素区域记为Ebone.

图8 区域估计

由区域项定义可知,区域项函数中表示将体素p标记为目标点(或背景点)的惩罚.分割骨组织时,在体素p是骨骼或在体素p是背景且时惩罚小,其他情况惩罚高.结合估算的区域,建立新的区域项函数为

(15)

新的区域项函数综合评分函数值和图像强度,相较于单一的基于图像强度的区域项函数具有更明显的特征和更高的可靠性.

2.3.2 改进边界项函数

同传统区域项函数一样,仅基于图像强度的边界项函数不能处理对比度低的骨边界,因此考虑使用评分函数设计新的边界项函数.

边界项中表示相邻体素p,q的一致性惩罚,p,q越相似(如在图像强度上相近)惩罚越大,反之则小.因此,基于评分函数提出新的边界一致性标准,当评分函数值相近时惩罚大,反之则小.定义边界项函数为

(16)

根据式(1)、(2)、(3)、(15)、(16),得能量函数E(A)为

(17)

根据上述原理对CT图建立图模型,使用最大流算法求解能量函数得出图模型中对应割线,该割线将图像分割为前景背景两部分,前景即为所求骨组织. 图9为分割股骨切片的示意图,图9(a)为原始切片;图9(b)为图割后的二值结果,图中可见,切片被分为前景、背景两部分,其中前景由两个不连通区域A和B构成,分别代表髋骨和股骨,两者边界清晰、互不交叉;图9(c)为将原图与二值图相叠加的效果图,图中可见,骨组织与软组成边界清晰,骨分割质量好,证明本文提出的基于骨边界增强滤波的图割算法是可行的.

图9 股骨切片初步分割结果

2.4 后处理

虽然上述骨边界增强滤波能够提高骨边界清晰度,但由于CT图像固有的部分体积效应,仍会导致初步分割结果中出现相邻骨组织粘连的情况,如图10所示,图10(a)是使用上文方法获得的初步二值结果,图中可见股骨区域与髋骨存在细微的连接点,导致两骨无法彻底分离.为分离相邻的骨骼,提出基于形态学和简易图割的后处理方法,将半径为R的球形元素作为卷积核对分割后的二值图像进行形态学腐蚀处理,如10(b)所示,对原图进行腐蚀操作,区域C被分割为D、E和孤岛3个区域.

图10 相邻骨分离示意图

对于骨的分离操作,目标是在边界体素最少时,将图像分割为F、G两个区域,其中,对图像做简易图割处理,步骤为

(1)添加硬约束,令腐蚀操作后的区域D和F分别为图割中的前景种子点和背景种子点;

(2)建立区域项函数,令区域项函数为

(18)

(3)为得到包含体素最少的边界,令边界项函数为;

(4)建立能量函数并使用最大流算法求解.

(5)将分离后的相邻骨赋予不同CT值作为标签.

简易图割处理结果如图10(c)所示,其中图10(b)中的孤岛区域被纳入区域D中,符合预期,证明本文提出的基于形态学与简易图割的后处理方法是有效的.

3 实验及结果

为验证本文算法,可对骨组织进行自动精确的分割,在CT影像数据集上进行实验,并与阈值分割、区域生长、传统图割进行对比.

3.1 实验数据集

实验数据集分为两组:第一组为内部数据集,来源于武汉大学中南医院,包含30例下肢CT图像,图像层内间距0.78 mm,层间间距1.25 mm,图像强度单位为HU. 以医学专家进行人工分割结果作为金标准;第二组为公开数据集,来自TCIA,名为CT Lymph Nodes,包含176例腹部CT图像.

3.2 实验环境及参数设置

实验基于Windows10系统,采用Visual Studio 2015开发平台,配置医学图像处理工具包ITK、可视化工具包VTK、QT框架和最大流库,使用C++语言开发,实现图切割的高效计算. 实验PC硬件配置为Intel(R)Core(TM)i7-8700K CPU@3.70GHz,16GB内存.

经过大量的实验得出本文算法及对比算法合适参数见表5.

表5 本文算法及对比算法主要参数

3.3 评价指标

医学图像分割领域中评价指标众多,能够从不同角度对算法的预测结果进行评估.本文采用Dice系数(Dice)、精确率(Precision)、召回率(Recall),以及F1分数(F1Score)作为评价指标.

(1)Dice系数:表示两个样本的重叠情况,设A为人工分割金标准的图像像素的集合,B分割算图像像素集合,其定义为

(19)

当A与B不相交时,Dice相似性系数为0;当A与B完全相交时,Dice相似度系数为1,即系数越大,图像分割质量越高.

(2)精确率:表示分割的骨区域中真实骨区域所占比例,其定义为

(20)

式中,TP表示骨被预测为骨,FP表示非骨被预测为骨.

(3)召回率:表示分割骨区域中的真实骨占实际骨区域的比例,其定义为

(21)

式中,FN表示骨被预测为非骨.

(4)F1分数:是精确率和召回率的调和平均数,同时考虑预测的骨区域和真实骨区域,评估结果全面,其定义为:

(22)

3.4 实验结果分析

3.4.1 骨分割结果对比

从数据集中选取5例,使用阈值分割、区域生长、传统图割、骨边界增强滤波的图割算法对股骨进行分割,分割结果分别如图11所示,其中(a)行是经裁剪后的原始CT图像,(b)行是由专业医师手工分割的实际股骨区域,作为实验的真值,(c)行是由阈值分割算法提取的股骨区域,(d)行是由区域生长算法提取的股骨区域,(e)行是由基于强度的传统图割算法提取的股骨区域,(f)行是本文算法提取的股骨区域.

图11 阈值分割、区域生长、传统图割和本文算法对股骨分割结果对比

对比可得,阈值分割算法分割的效果较差,是由于阈值分割的原理是将图像中强度在一定范围内的像素不经处理地提取出来,如四、五列所示,分割结果中不仅存在与骨组织强度相近的非骨像素,而且存在多种组织形成的不连通区域,这是导致阈值分割算法分割效果不理想的主要原因;区域生长对不连通区域的处理优于阈值分割,如第三列所示,是因为区域生长能提取图像连通区域,但其本质是基于图像强度,无法排除高强度的非骨像素,因此分割效果仍不理想;传统图割算法在骨干区域的分割效果较好,是由于骨干区域图像对比度高,但在关节区域由于第二节介绍的原因,无法准确提取股骨区域,如第四、五列所示;本文算法提取的股骨区域与实际股骨区域更为相近,经过比较可得本文算法分割结果明显优于其他方法.

从数据集中选取4例,分别对股骨、髋骨、胫骨、腓骨周围区域进行裁剪.使用4种方法分割并重建,如图12-15所示.图12为股骨重建对比图,其中图12(a)为实际股骨;图12(b)为阈值分割重建的股骨,模型存在多个不连通区域,且不能分离股骨和髋骨;图12(c)为区域生长重建的股骨,模型不连通区域减少,但仍不能分离股骨和髋骨;图12(d)为传统图割重建的股骨,模型与区域生长模型相近,不能分离股骨和髋骨;图12(e)为本文算法重建的股骨,模型与实际股骨在形状和大小十分贴近实际股骨,比较可得,本文算法对股骨的分割优于其他方法.图13为髋骨重建对比图,其中图13(a)为实际髋骨,分割难点在于髋骨区域骨组织较多,与髋骨邻近的有两侧股骨及中央骶骨,比较可得,本文算法能成功分割髋骨,分割效果优于其他方法.图14、15分别为胫骨、股骨重建对比图,比较可得本文算法对胫骨、股骨的分割效果优于其他方法.

图12 阈值分割、区域生长、传统图割和本文算法股骨分割重建对比

图13 阈值分割、区域生长、传统图割和本文算法髋骨分割重建对比

图14 阈值分割、区域生长、传统图割和本文算法胫骨分割重建对比

图15 阈值分割、区域生长、传统图割和本文算法腓骨分割重建对比

3.4.2 骨分割结果对比

为定量评价算法的可靠性,用阈值分割、区域生长、传统图割算法和本文算法分别分割骨组织,并计算4种骨用4种方法在评价指标Dice、Precision、Recall和F1 Score上的均值,并与Jcgs等[21]提出的3D U-Net比较,测试数据为TCIA中名为CT Lymph Nodes的数据集,评估结果见表6.

表6 4种方法骨分割评估结果

对比实验结果,可得阈值分割算法的精确率最低,表示实际骨区域占预测骨区域的比例较低. 区域生长算法较阈值分割算法表现较好,在精确率和F1分数上有一定的提高,但各向指标均不如本文算法. 传统图割算法的召回率高,精确率很低,这表明传统图割算法对骨区域的提取不够准确. 本文算法在保证高召回率的情况下精确率能达99.07%,对召回率和精确率的平衡较好,在综合评价指标Dice和F1 Score上有明显提升. 综合4种评估结果,本文算法的整体表现优于其他方法,且与3D U-Net算法在Dice指标上相近.

4 结论

针对骨组织的自动分割进行研究,基于Hessian矩阵提出骨边界增强滤波,增强骨边界清晰度;将滤波结果融入图割的能量函数中,提高图割算法的鲁棒性和效率;融合形态学和图割提出一种后处理方法,实现相邻骨的自动分离. 在实际应用中,该方法对骨分割处理较好,但对于其他组织的分割适应性差,因此在下一步的研究工作中,将面向多样的人体组织进行分割研究,如血管、肺、结核等,进一步提升算法的通用性.