基金项目:国家自然科学基金(61373077); 福建省自然科学基金(2015J01587); 福建省科技厅资助高校项目(JK2010056); 福建省教育厅项目(JB10160)
通信作者:chli@xmu.edu.cn
(1.厦门大学信息科学与技术学院,福建 厦门 361005; 2.龙岩学院信息工程学院, 福建 龙岩 364000; 3.龙岩市第二医院康复科,福建 龙岩 364000)
(1.College of Information Science and Engineering,Xiamen University,Xiamen 361005,China; 2.College of Information Engineering,Longyan University,Longyan 364000,China; 3.Department of Rehabilitation,Longyan Second Hospital,Longyan 364000,China)
cervical segmentation; level set; maximally stable extremal regions(MSER); least square method; structure features
DOI: 10.6043/j.issn.0438-0479.201706016
颈椎椎体的分割在颈椎图像处理中起着关键的作用,是颈椎病灶确定和辅助诊断的重要基础.针对颈椎椎体边缘特征复杂的特点,提出一种基于水平集和最大稳定极值区域(maximally stable extremal regions,MSER)融合的颈椎椎体分割方法.首先采用基于图像密集度分布的图像分割方法对图像进行粗分割,自动提取颈椎区域; 然后采用改进的水平集方法提取出颈椎椎体的前缘轮廓; 根据颈椎椎体后缘的局部稳定特征,采用改进的MSER方法提取出椎体的后缘高亮区域,并结合椎体结构特征,采用最小二乘法拟合出椎体的后缘曲线; 最后融合颈椎椎体前缘轮廓与后缘曲线,从而提取完整的颈椎椎体.实验结果表明,该方法能有效地分割和提取颈椎椎体,提取的后缘曲线接近专家手工提取的结果,可以为颈椎病的临床诊断提供更客观的诊断依据.
Segmentation of cervical centrum plays a key part in cervical image processing,and is an important basis for confirming cervical lesions and auxiliary diagnosis.Targeting on the complexity of cervical centrum margin,this paper presents cervical centrum segmentation based on level set and combination of maximally stable extremal regions(MSER).First,rough segmentation of image will be undergone based on image distribution density to extract cervical regions automatically.Then,the improved level set is adopted to extract lip sketch of the cervical centrum.In line with partial stable features of cervical centrum trailing edge,the improved MSER is adopted to extract highlighted areas in rip of centrum.At the same time,with combination of structure features of the centrum,least square method is taken for centrum rip curve fitting.Finally,cervical centrum lip sketch and trailing curve are combined for extraction of completed cervical centrum.Experiment results show that methods used in this paper can efficiently segment and extract cervical centrum,and the result of trailing curve is similar to that by manual segmentation,and can offer more objective diagnostic basis for clinical diagnosis of cervical spondylosis.
近年来,由于长时间面对电脑或手机,颈椎病的患病率逐年增高,发病年龄有明显年轻化的趋势[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方法进一步提取颈椎椎体的后缘,从而提取出完整的颈椎椎体.
针对颈椎椎体形态规律一致,后缘存在局部稳定特征的先验,文中提出一种融合水平集和MSER的颈椎椎体分割方法.算法的分割流程如图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射线图像中多是椒盐噪声,而且椎体的边缘信息特别是椎体后缘的高亮信息是本文中提取颈椎椎体的重要线索.中值滤波的核心思想是图像中像素灰度值设置为该点某领域窗口内所有像素灰度值的中值,该方法对消除颈椎图像的椒盐噪声非常有效,而且可以较好地保留图像的边缘细节,不破坏颈椎椎体后缘的局部高亮等稳定特征.
Chan和Vese[14]提出的C-V模型是一种典型的基于全局区域的水平集分割模型,C-V模型把目标与背景考虑为全域二聚类问题,其缺点主要是不能成功地分割亮度不均匀的图像.Caselles等[15]提出的测地线活动轮廓(GAC)模型是基于局部边界的水平集分割模型,缺点是不能成功地分割对比度低及边缘模糊的图像,容易造成边界泄露,且其存在对凹陷处难分割的问题.文献[16]表明,若用GAC模型中边界能量项代替C-V模型能量泛函中闭合边界曲线C的全弧长项,可获得更好的边界信息.假设定义域为Ω的图像I(x,y)被闭合边界曲线C划分为内部目标区inside(C)和外部背景区outside(C).融合区域和边界特征的能量项后,新的能量泛函表示如下[16]:
ECV-GAC(c1,c2,C)=λ1inside(C)(I-c1)2dxdy+
λ2outside(C)(I-c2)2dxdy+μ∮Cg(C)ds.(1)
其中,ci(i=1,2)表示曲线内部和外部的平均灰度值,λi≥0为能量项权重系数,系数μ>0,g为边界停止函数.
只有当轮廓线C演化到目标边界时,能量泛函才能达到最小值.优化式(1),即可得到未知数c1、c2以及最终分割轮廓线C的位置.
分析图2(d)可知,经过颈椎区域的提取,颈椎椎体形态较相似,每个椎体的前缘清晰完整,亮度分布均匀,椎体之间有间隙,呈现凹凸状.因此本文中采用C-V模型和GAC模型相结合的方法提取颈椎的前缘轮廓.
颈椎椎体后缘X射线图像实际为椎体后缘与两侧横突骨皮质的重叠影像,X射线成影后该区域的亮度值较大,形成局部高亮,同时团聚成具有一定面积大小、分布不连续的长条型区域.文献[12,17]的研究均表明MSER算法在稳定性、仿射不变性和对光照的适应性等方面都优于其他的区域特征提取算法.本文中从颈椎椎体后缘的局部高亮特征出发,利用MSER算法检测出高亮区域,从而实现颈椎椎体的后缘曲线提取.
MSER算法由Matas等[18]于2002年提出,MSER数学定义为:
对灰度图像I:区域DZ2→S,其中S={0,1,…,255},S是全序的,且满足非对称性、自反性和传递性,定义邻域关系:AD×D,即若∑di=1|pi-qi|≤1,则点p,q∈D相邻,并记为pAq,本文中取4邻域.
定义区域D的连通子集Q:对于p,q∈Q,都存在路径p,a1,a2,…,an,q,使得pAa1,a1Aa2,…,anAq相邻,这里ai∈Q,i=1,2,…,n; 区域Q的边界Q:Q={q∈D\Q:p∈Q:qAp},即Q当中的像素不属于子集Q,但与Q内至少一个像素相邻.则对于极值区域QD,p∈Q,q∈Q,有:
Q={MSER+(最大极值区域),I(p)>I(q),
MSER-(最小极值区域),I(p)<I(q).(2)
假设Q1,…,Qi-1,Qi,…为嵌套极值区域的序列,即Qi-1Qi.若稳定方程:
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最小外接矩形的分析,再根据外接矩形的参数判断是否是颈椎椎体后缘区域.
由于整幅图像中颈椎椎体后缘区域的像素点的像素值有着较大相似性,呈现高亮且不连续的短边缘特征,所以可认为颈椎的每一节椎体的后缘高亮区域可能对应一个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,则相邻中心点的平均距离为
d0=(∑k-1i=1|xi+1-xi|)/k.(4)
以d0再加上一个增量ζ0作为阈值,用以剔除不符合条件的噪声MSER.根据大量实验,本文中设定ζ0为0.2 mm.
策略三:剔除与椎体前缘距离不相符的MSER.经过策略一和策略二,绝大部分噪声区域的MSER被排除,但仍存在个别噪声MSER且会对后续处理产生较大干扰.进一步分析椎体结构特征,椎体的前后缘接近平行,椎体前后缘之间的距离应基本相同.如图3(b)所示,假设点(x1,y1)、(x2,y2)为后缘上的MSER对应矩形中心点坐标,根据椎体的刚性特征,则它们与对应的椎体前缘的横向距离应该接近相等,即s1≈s2,据此进一步可剔除噪声区域的MSER.具体操作步骤如下:1)在文中1.2提取颈椎椎体轮廓的基础上,对轮廓进行边缘跟踪并细化; 跟踪并细化的效果如图3(b)所示.2)同策略二,对经策略一处理后的所有矩形区域,按中心点位置的纵坐标大小对矩形进行排序,得到排序后的k个矩形Ri,其对应的矩形中心点坐标为(xi,yi),i=1,2,…,k.3)求k个矩形中心点坐标(xi,yi)横向对应的椎体前缘位置上的点坐标,假设坐标为(x'i,y'i).4)计算椎体的前后缘横向平均距离,有:
d1=(∑ki=1|xi-x'i|)/k.(5)
5)结合椎体的前后缘宽度先验知识,以d1加上一个增量ζ1作为阈值,剔除不符合条件的MSER区域.根据大量实验,本文中设定ζ1为0.7 mm.
通过提取椎体的MSER区域并排除噪声区域后,提取所有目标区域MSER的矩形中心点,采用最小二乘法拟合出椎体后缘曲线,但由于排除噪声区域后的MSER主要集中在颈椎第3~7节椎体,而且MSER 特征形状不规则,可能会导致后缘曲线的提取存在误差.为此,本文中采用对后缘曲线进行二次拟合的方法,首先在拟合出的后缘曲线上获取曲线的顶点位置坐标,并参照临床中常用的Boden氏测量方法[19]提取第2节颈椎椎体后上缘和第7节颈椎后下缘位置上的两点坐标,根据这两点坐标对曲线进行二次拟合,得到更准确的后缘曲线.
本文中颈椎椎体分割算法(简称本文算法)的具体流程如下:
1)采用基于图像密集度分布的分割算法自动提取颈椎区域;
2)对提取的颈椎区域进行图像增强;
3)采用C-V和GAC相结合的水平集算法(CV-GAC)提取颈椎椎体前缘轮廓;
4)采用MSER算法提取颈椎区域的MSER;
5)对提取的前缘轮廓进行跟踪并细化处理;
6)结合椎体结构刚性特征,通过文中提出的3种策略剔除噪声区域的MSER;
7)采用最小二乘法对剔除噪声区域后的所有MSER进行拟合,得到椎体后缘曲线;
8)参照Boden氏测量方法,分别提取第2节颈椎椎体后上缘(椎体后缘中偏上)和第7节颈椎后下缘位置上的两点坐标;
9)结合人工提取的两点坐标,对后缘曲线采用最小二乘法进行二次拟合;
10)融合前缘轮廓和后缘曲线,提取颈椎椎体区域.
本文中的图像分割实验是在Windows系统下,利用MATLAB R2014b实现的.Lenovo Thinkpad 笔记本具体配置为 CPU Intel Core(TM)i5 2.60 GHz,4 GB RAM,并采用 SPSS 22.0 软件进行统计学分析.
实验数据来源:本文中实验的所有颈椎X射线图像均由龙岩市第二医院放射科同一组专业技师完成,摄片条件相同(机器所设置的参数),且拍摄头颅侧位图像时要求受试者颈部姿势统一.实验选取颈椎病临床诊断中常见的生理曲度正常、反弓、变直以及旋转等4种类型X射线颈椎侧位图像共计170例进行算法验证.
根据本文算法流程,首先要提取椎体的前缘轮廓,图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)所示.
本研究中实验结果的好坏主要取决于后缘曲线提取的准确性,因目前文献尚未查到有公开自动分割方法用于颈椎后缘曲线提取的先例,无法进行相关的对比实验.在颈椎诊治中,颈椎曲度值与后缘曲线密切相关,医学通常根据颈椎曲度值来判断病情的严重程度与类型.为了验证本文算法的有效性和鲁棒性,实验中成立专家组,由专家组医生手动提取后缘曲线并测量出颈椎曲度值,与本文算法提取的后缘曲线的曲度值作比较,以颈椎曲度值这个重要指标来衡量分割算法的效果.专家组由放射科3位专业技师组成,采用目前常用的Boden氏测量方法,分别对粗分割后的170例颈椎侧位图片手工提取其曲度值,最后用3位专家的平均值进行比较; 本文算法在提取后缘曲线后,借助医学影像信息系统(picture archiving and communication systems,PACS)分别获取170例图片的曲度值.实验中的结果采用SPSS 22.0进行统计分析,曲度值采用平均值±标准差((-overx)±s)表示,p为统计相关性.部分分割结果如图6所示,统计分析结果如表1所示.
图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 本文算法和手动分割测得的颈椎曲度
Tab.1 Cervical curvature measured with proposed algorithm and manual segmentation
颈椎侧位X射线图像是目前颈椎诊断中最常用的影像学检查方法,颈椎椎体的分割在X射线颈椎图像数据测量中起着关键的作用.但是因为测量方法受限,临床实际操作中普遍采用人工分割为主,人工分割主观性较强、效率低,对此本文中提出了基于水平集和MSER的颈椎椎体分割方法.从实验结果看,该方法能完整提取颈椎椎体,对正常、反弓及变直的颈椎椎体特别有效,尤其是后缘曲线的提取与专家手工提取的结果非常接近,可为颈椎曲度的测量、相邻椎体屈曲度的测量以及颈椎骨龄判定中的测量等实际的测量操作提供较准确且客观的数据.因此,本文算法在颈椎的辅助诊断中具有一定的临床实用价值.