基于水平集和最大稳定极值区域的颈椎椎体分割方法

兰添才,陈　俊,张怡晨,李翠华*

(1.厦门大学信息科学与技术学院,福建 厦门 361005;2.龙岩学院信息工程学院,福建 龙岩 364000;3.龙岩市第二医院康复科,福建 龙岩 364000)

近年来,由于长时间面对电脑或手机,颈椎病的患病率逐年增高,发病年龄有明显年轻化的趋势[1].X射线、CT及磁共振成像(magnetic resonance imaging,MRI)等影像学检查已经成为颈椎病诊断的常规辅助手段,相较于CT和MRI,X射线检查因其最简便、直观、经济且有效,目前仍然是颈椎病早期诊断的主要检查方法[2].颈椎图像的精确分割可以辅助医生判断病情，量化分析病灶区域,为颈椎病诊断提供可靠的依据,在颈椎病的临床诊断中有着非常重要的研究价值.然而,颈部是人体几何特性和运动特性最为复杂的部分之一,由于患者个体的差异性、颈椎结构本身的复杂性和椎体周围组成物质的非均匀性等,导致难以实现颈椎图像的精确分割.对颈椎的分割,已有一些学者做了相关研究[3-5].如,Chen等[3]对颈椎形态进行了定量分析,但其采用人工定点分割,依赖于操作者的主观经验,且耗时长、重复性低;蓝智聪等[4]利用抛物线分别对椎体的四条边缘进行拟合,进而提取颈椎椎体,其提取过程采用人机交互;赵晓光等[5]在手工截取颈椎区域的基础上,提出基于边缘和角点的算法提取椎体左右边缘,并分别用直线和曲线拟合出椎体的上下边缘,但其鲁棒性有待提高,仅对颈椎的总体排列和其前方的咽后壁呈“直线”型的颈椎有较好的效果.

水平集方法在处理医学图像分割问题时具有良好的性能,很多学者在这方面取得了一定的研究成果[6-10].如,梁礼明等[6]将水平集模型应用在眼底图像血管的分割中取得了良好的效果;Jeffee等[7]和Lu等[8]把水平集方法应用于超声图像的子宫颈以及宫颈细胞的分割中,成功提取宫颈以及成功分割出单个宫颈细胞;Wang等[9]提出了一种基于概率图谱与形状强度水平集相结合的分割方法并用于肝脏CT图像的自动分割,取得了类似于手工分割的肝脏边界;Cheng等[10]提出了一种基于B-Snake模型的血管CT图像自动分割方法,算法可以分割出类似于手工获得的血管边界.医学图像分割通常采用多种分割方法的融合[11].水平集方法对于对比度低或者边缘比较模糊的图像分割不准确,而颈椎椎体后缘相对于前缘,明显具有对比度低、边缘模糊的特点.考虑到后缘具有亮度高、不连续等稳定的局部特征.而对于图像的局部不变特征提取,最大稳定极值区域(maximally stable extremal regions,MSER)在大部分情况下性能最佳[12].本研究在水平集分割的基础上,提出采用改进的MSER方法进一步提取颈椎椎体的后缘,从而提取出完整的颈椎椎体.

图2　自动提取颈椎区域Fig.2 Automatic extraction of cervical regions

1　算法流程及原理

针对颈椎椎体形态规律一致,后缘存在局部稳定特征的先验,文中提出一种融合水平集和MSER的颈椎椎体分割方法.算法的分割流程如图1所示.

图1　分割流程图Fig.1 The flow chart of segmentation

1.1　颈椎区域的自动提取与图像增强

颈椎侧位X射线图像来源于拍摄的头颅定位侧位片,如图2(a)所示.原图包含整个头颅的内容,而颈椎病诊断临床中所关注的区域是七小节颈椎椎体所在的区域.颈椎区域的提取可以避免对全图进行计算,提高处理速度.如何快速而准确地将颈椎区域从原图中分离出来,是整个临床和研究分析过程的首要步骤.赵晓光等[5]采用手工提取颈椎区域费时费力.根据颈椎侧位X射线图像中椎体位置的布局,本文中提出一种基于图像密集度分布的自动提取算法,根据阈值分割结果的密集度分布来自动提取颈椎区域.具体的提取步骤如下:

1) 降维处理.将图像统一从1 792像素×2 392像素降维到600像素×800 像素,降低原图分辨率.

2) 采用最大类间方差法(OTSU)对降维后的图像进行阈值分割(粗分割),分割结果如图2(b)所示.

3) 对阈值分割后的结果分别进行水平和垂直方向的投影，投影直方图如图2(c)所示.

4) 根据水平和垂直方向的投影直方图自动获取阈值,并将图像大小归一化为192像素×512像素，颈椎区域自动提取效果如图2(d)所示.

由图2(b)可知,对阈值分割后的结果进行垂直方向投影,颈椎部分投影值最高;对分割图进行水平方向投影,颈椎部分投影具有中间平稳的特性.图2(c)进一步验证了这个特征,水平方向投影的结果得到颈椎区域的水平方向的位置,垂直方向的投影结果得到颈椎区域的垂直方向位置.从图2(d)可知,提取结果完整包括七小节颈椎椎体,且缩小了范围,排除了上方颅骨和左上方下颌骨等较多的干扰.

受到采集系统和采集条件如光照等诸多因素的影响,采集到的X射线图像的质量不一致,噪声较多,这会对后续的图像处理产生不利的影响,因此有必要对粗分割后的颈椎区域进行图像增强.文中的图像增强主要包括以下2种方法:

1) 根据文献[13],对颈椎区域图像进行空域上的同态滤波以减少颈椎X射线图像拍摄时不均匀光照产生的影响.

2) 采用中值滤波对图像进行去噪，滤波算子为3×1.由先验知识可知,X射线图像中多是椒盐噪声,而且椎体的边缘信息特别是椎体后缘的高亮信息是本文中提取颈椎椎体的重要线索.中值滤波的核心思想是图像中像素灰度值设置为该点某领域窗口内所有像素灰度值的中值,该方法对消除颈椎图像的椒盐噪声非常有效,而且可以较好地保留图像的边缘细节,不破坏颈椎椎体后缘的局部高亮等稳定特征.

1.2　基于水平集的颈椎前缘轮廓提取

Chan和Vese[14]提出的C-V模型是一种典型的基于全局区域的水平集分割模型,C-V模型把目标与背景考虑为全域二聚类问题,其缺点主要是不能成功地分割亮度不均匀的图像.Caselles等[15]提出的测地线活动轮廓(GAC)模型是基于局部边界的水平集分割模型,缺点是不能成功地分割对比度低及边缘模糊的图像,容易造成边界泄露,且其存在对凹陷处难分割的问题.文献[16]表明,若用GAC模型中边界能量项代替C-V模型能量泛函中闭合边界曲线C的全弧长项,可获得更好的边界信息.假设定义域为Ω的图像I(x,y)被闭合边界曲线C划分为内部目标区inside(C)和外部背景区outside(C).融合区域和边界特征的能量项后,新的能量泛函表示如下[16]:

(1)

其中,ci(i=1,2)表示曲线内部和外部的平均灰度值,λi≥0为能量项权重系数,系数μ>0,g为边界停止函数.

只有当轮廓线C演化到目标边界时,能量泛函才能达到最小值.优化式(1),即可得到未知数c1、c2以及最终分割轮廓线C的位置.

分析图2(d)可知,经过颈椎区域的提取,颈椎椎体形态较相似,每个椎体的前缘清晰完整,亮度分布均匀,椎体之间有间隙,呈现凹凸状.因此本文中采用C-V模型和GAC模型相结合的方法提取颈椎的前缘轮廓.

1.3　基于MSER的颈椎后缘曲线提取

颈椎椎体后缘X射线图像实际为椎体后缘与两侧横突骨皮质的重叠影像,X射线成影后该区域的亮度值较大,形成局部高亮,同时团聚成具有一定面积大小、分布不连续的长条型区域.文献[12,17]的研究均表明MSER算法在稳定性、仿射不变性和对光照的适应性等方面都优于其他的区域特征提取算法.本文中从颈椎椎体后缘的局部高亮特征出发,利用MSER算法检测出高亮区域,从而实现颈椎椎体的后缘曲线提取.

1.3.1颈椎MSER区域提取

MSER算法由Matas等[18]于2002年提出,MSER数学定义为:

定义区域D的连通子集Q:对于∀p,q∈Q,都存在路径p,a1,a2,…,an,q,使得pAa1,a1Aa2,…,anAq相邻,这里ai∈Q,i=1,2,…,n;区域Q的边界∂Q:∂Q={q∈DQ:∃p∈Q:qAp},即∂Q当中的像素不属于子集Q,但与Q内至少一个像素相邻.则对于极值区域Q⊂D,∀p∈Q,∀q∈∂Q,有:

(2)

假设Q1,…，Qi-1，Qi,…为嵌套极值区域的序列,即Qi-1⊂Qi.若稳定方程:

q(i)=|Qi+△-Qi-△|/|Qi|，

(3)

其中,Qi表示灰度阈值为i时的连通区域,△为灰度阈值的微小增量,q(i)表示区域Qi的面积变化率.若i*处存在局部最小值,则极值区域Qi*即为MSER.

MSER提取过程具体包括以下步骤:1) 对给定的图像采用桶排序(bucket sort),按照亮度大小对所有像素点进行排序得到灰度值递增的像素点序列；2) 根据式(2)检测极值区域；3) 根据式(3)提取MSER；4) 由先验可知,颈椎椎体后缘多为长条形高亮区域,所以采用矩形拟合方法把问题转化为对MSER最小外接矩形的分析,再根据外接矩形的参数判断是否是颈椎椎体后缘区域.

1.3.2噪声区域分析与剔除

由于整幅图像中颈椎椎体后缘区域的像素点的像素值有着较大相似性,呈现高亮且不连续的短边缘特征,所以可认为颈椎的每一节椎体的后缘高亮区域可能对应一个MSER,因为不连续也可能对应多个MSER.但是因为颈椎背景复杂,经过粗分割后的颈椎区域还包括上方的部分颅骨、左上方部分下颌骨以及后缘周围的棘突、横突等,这些区域也可能存在高亮部分,在对图像进行MSER提取后,除了在椎体后缘所在区域提取到MSER外,在上述这些非目标区域可能也会提取到一些MSER.为了有效剔除这些噪声区域的MSER,同时又能最大限度地保留椎体后缘上的MSER,文中提出以下策略:

策略一:剔除不符合高度位置的MSER.从椎体结构的先验可知,第1,2节椎体结构特殊,故本文中主要提取第3～7节颈椎椎体,所以特征点也主要取自第3～7节颈椎椎体,第2节椎体以上的MSER作为噪声区域剔除.虽然由于个人体质、形变和拍摄角度导致椎体位置不尽相同,但是各节椎体高度大致相同,共7节椎体,本文中取整体椎体高度的2/7作为阈值剔除第1，2节椎体中的MSER.

策略二:剔除不满足横向间距的MSER.由颈椎结构先验可知,颈椎椎体后缘的纵向偏移应在±15°范围内,所以相邻椎体之间的横向距离也应在较小距离范围内.如图3(a)所示,假设(x1,y1)、(x2,y2)为相邻椎体MSER矩形的中心点坐标,则两个中心点的横向距离|x2-x1|应在一个较小的阈值范围内.按纵轴方向,通过计算相邻MSER对应矩形的中心位置的横向偏移,可剔除不满足条件的MSER.具体操作如下:首先扫描图像中经策略一处理之后的所有经过区域拟合后的矩形区域,分别计算每个矩形的中心点位置坐标,然后按照中心点位置的纵坐标大小对矩形进行排序.假设共有k个矩形Rk,其对应的中心点位置坐标为(xi,yi),i=1,2,…,k，则相邻中心点的平均距离为

(4)

图3　噪声区域剔除策略示意图Fig.3 Schematic diagram of the noise regions excluding strategy

以d0再加上一个增量ζ0作为阈值,用以剔除不符合条件的噪声MSER.根据大量实验,本文中设定ζ0为0.2 mm.

(5)

5) 结合椎体的前后缘宽度先验知识,以d1加上一个增量ζ1作为阈值,剔除不符合条件的MSER区域.根据大量实验,本文中设定ζ1为0.7 mm.

1.3.3颈椎后缘曲线提取

通过提取椎体的MSER区域并排除噪声区域后,提取所有目标区域MSER的矩形中心点,采用最小二乘法拟合出椎体后缘曲线,但由于排除噪声区域后的MSER主要集中在颈椎第3～7节椎体,而且MSER 特征形状不规则,可能会导致后缘曲线的提取存在误差.为此,本文中采用对后缘曲线进行二次拟合的方法,首先在拟合出的后缘曲线上获取曲线的顶点位置坐标,并参照临床中常用的Boden氏测量方法[19]提取第2节颈椎椎体后上缘和第7节颈椎后下缘位置上的两点坐标,根据这两点坐标对曲线进行二次拟合,得到更准确的后缘曲线.

1.4　算法流程

本文中颈椎椎体分割算法(简称本文算法)的具体流程如下:

1) 采用基于图像密集度分布的分割算法自动提取颈椎区域;

2) 对提取的颈椎区域进行图像增强;

3) 采用C-V和GAC相结合的水平集算法(CV-GAC)提取颈椎椎体前缘轮廓;

4) 采用MSER算法提取颈椎区域的MSER;

5) 对提取的前缘轮廓进行跟踪并细化处理;

6) 结合椎体结构刚性特征,通过文中提出的3种策略剔除噪声区域的MSER;

7) 采用最小二乘法对剔除噪声区域后的所有MSER进行拟合,得到椎体后缘曲线;

8) 参照Boden氏测量方法,分别提取第2节颈椎椎体后上缘(椎体后缘中偏上)和第7节颈椎后下缘位置上的两点坐标;

9) 结合人工提取的两点坐标,对后缘曲线采用最小二乘法进行二次拟合;

10) 融合前缘轮廓和后缘曲线,提取颈椎椎体区域.

2　实验结果与分析

2.1　实验环境与实验数据

本文中的图像分割实验是在Windows系统下,利用MATLAB R2014b实现的.Lenovo Thinkpad 笔记本具体配置为 CPU Intel Core(TM) i5 2.60 GHz,4 GB RAM,并采用 SPSS 22.0 软件进行统计学分析.

实验数据来源:本文中实验的所有颈椎X射线图像均由龙岩市第二医院放射科同一组专业技师完成,摄片条件相同(机器所设置的参数),且拍摄头颅侧位图像时要求受试者颈部姿势统一.实验选取颈椎病临床诊断中常见的生理曲度正常、反弓、变直以及旋转等4种类型X射线颈椎侧位图像共计170例进行算法验证.

2.2　实验结果分析

根据本文算法流程,首先要提取椎体的前缘轮廓,图4分别给出了颈椎生理曲度正常、反弓、变直以及旋转4种常见类型的X射线原图(图4(a)),CV-GAC提取的颈椎前缘轮廓效果图(图4(b)),专家组医生手工分割的颈椎前缘效果图(图4(c)).实验中基于区域的能量函数权重取常用值,即λ1=λ2=1;在实验中发现μ的取值对椎体右侧轮廓影响较大,对左侧轮廓影响较小,因椎体右侧结构导致右侧背景较复杂.而本研究关注的是椎体的左侧轮廓,实验中μ的值取为500.

从4(b)可知,本文中提出的自动提取颈椎区域的方法简单有效,提取的颈椎区域完整;采用CV-GAC方法提取的颈椎椎体前缘轮廓曲线非常贴近颈椎真实前缘,几乎与颈椎椎体的前缘重合.与图4(c)比较可知,采用CV-GAC方法提取的颈椎椎体前缘接近专家组手工分割的结果.对于存在较严重病变,如旋转扭曲的颈椎病变图像,虽然有误检,但实验中发现这种误检多发生在旋转图像的第2节椎体,是因为旋转使舌骨在该处形成重影所致,而其他椎体由于具有相对稳定的刚性结构特征,同样能够得到贴近椎体前缘的较完整轮廓曲线,这将为后续颈椎椎体后缘曲线的准确提取起关键作用.

图4　CV-GAC方法与手工分割方法提取颈椎前缘轮廓Fig.4 Front edge of cervical contour extracted by CV-GAC and manual segmentation method

其次,要提取颈椎区域的MSER并排除噪声区域的MSER,这是后缘曲线提取的关键.图5(a)是4种类型的颈椎区域MSER提取效果图；图5(b)所示为剔除噪声区域MSER后的效果图,图中被矩形框住的区域表示选中区域,没有被矩形框住的表示被剔除的噪声区域MSER.

从图5(a)可知,各类颈椎图像中第2～7节椎体后缘目标区域的绝大部分高亮MSER均能被有效检出,尽管旋转图像较特殊,仍能有效检出.但同时有较多噪声区域的MSER被检测到,且噪声区域的MSER主要集中在图像上部颅骨等处,椎体周围检测出的噪声区域的MSER较少.对于图像上部颅骨等处的噪声区域的MSER采用策略一容易剔除;而椎体周围检测出的噪声区域的MSER虽然较少,但它对后续的椎体后缘曲线拟合影响较大,必须剔除.本文中剔除噪声区域MSER的3种策略都是基于形状结构信息而提出.实验表明,采用本文中提出的3种策略能有效剔除噪声区域的MSER,同时又能完整保留椎体后缘高亮区域的MSER.剔除效果如图5(b)所示.

图5　颈椎区域的MSER提取与噪声区域的剔除Fig.5 MSER extraction of cervical reguons and exclusion of the noise regions

图6　本文算法与手工分割结果比较Fig.6 Results comparison of proposed segmentation and the manual segmentation

图6(a)是采用最小二乘法对排除噪声后剩下的所有MSER区域进行拟合,得到后缘曲线,并与前缘轮廓进行融合后的结果.为使椎体后缘拟合得更准确,便于临床诊断数据的提取,参照临床中常用的Boden氏测量方法,分别提取第2节椎体后缘中偏上位置和第7节椎体后下缘位置上的两点坐标,并对椎体后缘曲线进行二次拟合,拟合的效果如图6(b)所示.图6(c)给出了专家组医生采用手工分割方法提取的颈椎后缘曲线效果图.根据文献[19]中的诊断标准,从图6(b)可知,正常、反弓和变直3种类型的椎体后缘曲线拟合的效果都很好,与实际的颈椎椎体的后缘都非常的贴近;虽然旋转类型图像本身特殊,但从拟合的角度看也比较贴近.与6(c)比较可知,拟合结果接近专家组医生手工分割的结果.

从表1也可知,本文算法与专家组手工分割提取的曲度间不存在显著性差异(p>0.05),统计结果表明本文算法提取的后缘曲线接近专家组医生手工分割提取的结果.

为进一步验证本文算法的有效性,根据陈银海等[20]对曲度异常的颈椎分类标准,将颈椎图像分为正常(7 mm≤曲度值≤17 mm)、反弓(曲度值≤1 mm)和变直(1 mm<曲度值<7 mm)3类.并根据之前提取的曲度值,从170例中筛选出区分度较好的正常图像55例、反弓图像17例以及变直图像28例(共100例),采用SPSS 22.0再次进行统计学分析比较,统计结果如表2所示.从表2可知,采用本文算法提取的3种类型的颈椎后缘曲线的曲度值均与手工分割结果非常的接近,这进一步表明本文算法提取的后缘曲线接近专家组医生手工提取的结果,特别是对正常、反弓及变直的颈椎图像可获得理想的效果.

表1　本文算法和手动分割测得的颈椎曲度

表2　3种类型颈椎不同方法测得的颈椎曲度

注：横向比较有显著性差异(p<0.05);纵向比较无显著性差异(p>0.05).

3　结　论

颈椎侧位X射线图像是目前颈椎诊断中最常用的影像学检查方法,颈椎椎体的分割在X射线颈椎图像数据测量中起着关键的作用.但是因为测量方法受限,临床实际操作中普遍采用人工分割为主,人工分割主观性较强、效率低,对此本文中提出了基于水平集和MSER的颈椎椎体分割方法.从实验结果看,该方法能完整提取颈椎椎体,对正常、反弓及变直的颈椎椎体特别有效,尤其是后缘曲线的提取与专家手工提取的结果非常接近,可为颈椎曲度的测量、相邻椎体屈曲度的测量以及颈椎骨龄判定中的测量等实际的测量操作提供较准确且客观的数据.因此,本文算法在颈椎的辅助诊断中具有一定的临床实用价值.

参考文献：

[1]张少群,李义凯.颈椎病研究的历史沿革[J].中国康复医学杂志,2016,31(11):1273-1276.

[2]陈超.探讨不同的影像学方法诊断颈椎病的临床价值[J].检验医学与临床,2013,10(20):2694-2695.

[3]CHEN L L,LIN J X,XU T M,et al.The longitudinal sagittal growth changes of maxilla and mandible according to quantitative cervical vertebral maturation[J].Journal of Huazhong University of Science and Technology (Medical Sciences),2009,29(2):251-256.

[4]蓝智聪,陈莉莉,许向阳,等.计算机辅助颈椎分析系统识别颈椎标志点的准确性和重复性探讨[J].临床口腔医学杂志,2010,26(9):526-530.

[5]赵晓光,林久祥,王晴竹.颈椎影像计算机自动识别骨龄诊断系统的建立[EB/OL].(2013-09-03)[2017-06-17].http:∥www.paper.edu.cn/releasepaper/content/201309-17.

[6]梁礼明,黄朝林,石霏,等.融合形状先验的水平集眼底图像血管分割[J/OL].计算机学报，2016.(2016-11-25).[2017-06-17].http:∥www.cnki.net/kcms/detail/11.1826.TP.20161125.2311.002.html.

[7]JEFFEE A I B,PAHL C,ABDULJABBAR H N,et al.Cervical segmentation in ultrasound image using level-set algorithm[C]∥Advances in Biomedicine and Health Science.Johor:[s.n.],2013:25-30.

[8]LU Z,CARNEIRO G,BRADLEY A P.An improved joint optimization of multiple level set functions for the segmentation of overlapping cervical cells[J].IEEE Transactions on Image Processing,2015,24(4):1261-1272.

[9]WANG J K,CHENG Y Z,GUO C Y,et al.Shape-intensity prior level set combining probabilistic atlas and probability map constrains for automatic liver segmentation from abdominal CT images[J].International Journal for Computer Assisted Radiology and Surgery,2016,11(5):817-826.

[10]CHENG Y Z,HU X,WANG J,et al.Accurate vessel segmentation with constrained B-Snake[J].IEEE Transactions on Image Processing,2015,24(8):2440-2455.

[11]王阳萍,杜晓刚,赵庶旭,等.医学影像图像处理[M].北京:清华大学出版社,2012:69-109.

[12]MIKOLAJCZYK K,TUYTELAARS T,SCHMID C,et al.A comparison of affine region detectors[J].International Journal of Computer Vision,2005,65(1/2):43-72.

[13]闻莎,游志胜.性能优化的同态滤波空域算法[J].计算机应用研究,2000,17(3):62-65.

[14]CHAN T F,VESE L.Active contours without edges[J].IEEE Transactions on Image Processing,2001,10(2):266-277.

[15]CASELLES V,KIMMEL R,SAPIRO G.Geodesic active contours[J].International Journal of Computer Vision,1997,22(1):61-79.

[16]TAO W B,TAI X C.Multiple piecewise constant with geodesic active contours(MPC-GAC) framework for interactive image segmentation using graph cuts optimization[J].Image and Vision Computing,2011,29(8):499-508.

[17]MIKOLAJCZYK K,TUYTELAARS T,SCHMID C,et al.A comparision of affineregion detectors[J].International Journal of Computer Vision,2005,65(1):43-72.

[18]MATAS J,CHUM O,URBAN M,et al.Robust wide baseline stereo from maximally stable extremal regions [C]∥Proceedings of the British Machine Vision Conference.Cardiff:[s.n.],2002:384-393.

[19]王涛,周理乾,孙孟锟,等.6种颈椎曲度测量方法的可信度及可重复性比较[J].中国脊柱脊髓杂志,2015,25(4):323-327.

[20]陈银海,姚红华,杨忠.颈椎曲度的X线测量在颈椎病康复评定中的应用价值[J].中国康复,2007,22(3):156-158.