局部信息约束的形状导向水平集算法对下颌管的分割

杨 玲 侯小叶 王中科 饶妮妮

（电子科技大学生命科学与技术学院，成都 610054）2（成都信息工程学院电子工程学院，成都 610225）3（成都信息工程学院网络工程学院，成都 610225）

引言

牙种植手术是口腔医学领域中应用最广泛的手术之一，如何确保在理想的位置植入种植体，最大限度地利用骨量，并确保邻近重要生理组织结构的安全，仍是临床医学界的一大难题［1］。下牙槽神经穿过整个下颌神经管，直接决定着种植体的植入位置，若植入不当会引起下牙槽神经束损伤，造成下唇麻木甚至丧失感觉。由于下牙槽神经的CT值与包含在下颌管中的许多组织非常相似，区分异常困难，临床中主要通过下颌管位置来确定种植体的植入位置，以避免对下牙槽神经的识别。下颌管是下颌骨内一条主要的血管神经通道，其结构不易被界定，且常与相邻区域粘连，因而其识别也成为了口腔界一直关注的问题，并影响着牙种植技术的发展［2］。例如，1987年 Tamas等对下颌管的位置进行了X线观察测量［3］，1993年 Obradovic等通过活体测量得到下颌管的位置走向［4］，李宁毅、浦铁民、卜令学等分别通过解剖成人的下颌骨样本总结测出了下颌管在下颌骨中的位置［5－7］，王吉昌等随机选取60例青年女性，对其下颌骨的CT断层图像进行测量总结，获得了下颌管在下颌骨中位置的应用解剖资料［8］。但因为解剖学个体差异较大及医生经验不同，均属经验性总结，对临床指导意义有限。近期，比利时Materialise公司专为口腔种植研发的交互式计算机辅助种植手术软件 Simplant系统，能够通过CT图像重建三维模型，医生可以根据下颌管端口的位置大体画出下颌管在下颌骨中的走向；但此项功能主要通过头颅CT切片数据进行三维重建，并需根据解剖结构对下颌管进行手动标记，对医生的经验要求较高。为进一步推进下颌管的自动提取，Yau等在重建下颌管时根据下颌管区域的相似性，提出了一种自适应区域增长的方法来分割下颌管［9］，但这种方法对种子点的选取比较敏感，同时对较模糊的成像情况其区域增长条件较难控制，容易造成分割失败。Bresson等提出基于Mumford-Shah模型和全局图像信息的形状导向算法模型，在图像分割中起到了很好的效果［10］，但是由于下颌管相对于周围组织对比度较低，在选取宽泛的初始轮廓时，活动轮廓不能在下颌管区域边缘停止演化。针对这种情况，本研究利用PCA获取下颌管先验轮廓的主分量，改进形状导向水平集算法中的能量函数模型，引入下颌管的局部信息，将局部、边界、区域和形状信息结合，共同约束水平集的演化过程，使其有效地分割出下颌管。

1 材料与方法

1.1 数据选择

所用的CT图像数据来自华西口腔医院，选取对象为2组（女）健康下颌要求牙种植手术的患者，身体健康，下颌骨左右基本对称，无畸形。CT机为日本森田制造所的3DX MULTI-Image Micro CT，扫描参数为电源AC100 V，最大耗电量为2.0 kVA，X线管焦点尺寸为 0.5 mm×0.5 mm，旋转角度为360°，扫描时间为18 s，扫描范围为直径40 mm×高30 mm的圆柱空间，扫描从颅顶到下颌范围共480层，立体像素尺寸为 0.25 mm×0.25 mm×0.25 mm，选取的图像是头颅下半部分的260张CT图像，原始图像的尺寸为683像素 ×683像素，像素深度为16位。

1.2 方法

1.2.1 下颌侧断面的重建

CT图像是层面图像，常用的是轴位图像，为分割下颌管，需要重建出下颌的侧断面，可以直观地看到下颌管的解剖结构。应用文献［9］重建下颌侧断面的方法来重建出下颌的侧断面，在CT图像的轴位图面（见图1（a））上沿着牙弓手动选取一系列点，用 B样条曲线拟合出牙弓线（见图1（a）中曲线），然后根据式（1）做牙弓线上任意一点的法线（见图1（a）中的直线）。

式中，Sx（p），Sy（p）分别表示 B样条曲线S（p）在点p的x和y方向上的向量。

取260张CT图像同样法线位置的数据，重建下颌骨侧断面结果（如图1（b）所示），其中下颌管的侧断面是闭合圆形。

1.2.2 先验形状模型构建

为准确分割下颌管，需用先验形状进行引导。为了将形状信息结合到图像分割过程中，基于一组训练集，应用文献［11］中提出的基于符号距离函数的方式来描述先验形状，并采用主成分分析法（principle component analysis，PCA）构建出下颌管的统计模型。

为了捕获形状信息，首先利用水平集方法表示目标形状，在符号距离函数的集合{Φ1，Φ2，…，Φn}之上建立一个形状模型，其平均水平集函数为

然后，确定一组最佳正交基 {U1，U2，…，Un}来表示n个训练集 {Φ1，Φ2，…，Φn}的最小二乘拟合，矢量Ui是由协方差矩阵 Σ =（1/n）MMT经过SVD分解得到的，其中M是一个由n个对齐训练集的符号距离函数构成的矩阵。于是，训练集的每一个形状矢量可以表示为

式中，Uk是一个由U的前k列组成的矩阵，α是由k个主成分特征值构成的k维向量。

在由α描述的形状服从高斯分布的假设下，可以计算一个给定形状的概率，即

式中，A（n×n）是由矩阵的特征值构成的一个对角矩阵，Ak是由A的前k行和列组成的矩阵。

1.2.3 局部信息约束的形状导向水平集分割模型

形状导向水平集分割模型是由Bresson等提出［10］，其主要思想是在水平集演化时增加目标物的形状参数，以此控制水平集演化时不受目标物遮挡等的影响，其能量泛函定义为

式中，Fs、Fb和Fr分别表示形状信息、边界信息和区域信息，可计算为

式中，C是活动轮廓，其作用是将图像 I（x，y）划分为若干近似同质区域；Xt是参数向量XT给出的几何变换；Ωi和Ωo是Φ*（由式（4）给定）零水平集的内部和外部区域；g是边缘提取函数，即 g：[0，+∞）→R+，g（0）=1，g（x）→0（x→∞），单调下降；βs、βb、βr是各个能量项的权重系数。根据下颌管的特征，设其先验形状为圆形几何形状。

上述模型的能量函数仅仅利用了图像的区域、边缘梯度和形状信息，缺少曲线所在处的灰度值的局部信息，这对于处理同质区域内强度分布不均匀（即灰度渐变）的图像，必然在阈值选取上会忽略掉某些局部信息，特别是在目标边界微弱变化的区域，当闭合曲线到达的时候，很容易受到由整体信息决定的阈值影响，造成分割的最终结果不尽人意。进一步，由于下颌管相对于周围组织对比度较低，且结构偏小，对初始轮廓要求较高，因此采用上述普通的形状导向水平集分割方法难以得到理想的分割效果。结合下颌管与周围组织的模糊特性，在形状导向水平集分割方法的基础上引入局部信息［12］，即在形状导向能量泛函的基础上，加入曲线到达图像局部强度的分布信息，使得算法对于强度分布不均匀的图像或区域仍然有着很好的处理能力。此时目标分割模型的能量泛函定义为式中，βl是其权重系数，Fl表示增加的局部信息分量，可计算为

此时，新的能量函数应用水平集方程表示为

式中，δ（x）=H'（x）。

若在式（13）中令λ1=λ2=1，则在构造高斯核函数时，其模板的大小（n×n）选取必须满足n＞4σ。σ越小，处理的效果越好。选取σ=3，n=13 。

2 结果

根据图2 OsinX软件全景片的定位和口腔医师的指导，为缩小分析范围，对图1（b）下颌骨侧断面的重建结果截取以下颌管为中心、尺寸为101像素×101像素的图像，结果如图3（a）所示。然后，提取二值化的下颌管样本，如图3（b）所示。接着，应用先验形状模型，取20个训练样本轮廓曲线（浅色线），得到平均轮廓（深色线），如图3（c）所示，其轮廓放大图像如图3（d）所示。选取较宽松的初始轮廓区域为 SDF（35：65，35：65）=1，如图3（a）中的矩形框所示。形状能量项的几何变换参数包括：先验形状放大的幅度值 mu0=1/1.1，平移分量 tx0=－2、ty0=－1。最终利用局部信息约束的形状导向水平集分割模型，得到两组患者的分割结果，分别如图4和图5所示，其中图（a）～图（e）分别对应图2中穿过曲线所示下颌管位置的1～5标签处的下颌骨侧断面，图（f）～图（j）和图（k）～图（o）均为对应图（a）～图（e）的下颌管区域的放大图像。在图（f）～图（j）中，实线表示的活动轮廓曲线基本在虚线表示的下颌管边缘停止演化，勾勒出了下颌管的形状。其中，参数设置为 βs=0.1，βb=0.2，βr=4，βl=0.01，μ =50，Δt=0.4，共迭代 7 次。

图2 OsinX软件全景片Fig.2 Panoramic image in OsinX software

图3 下颌管平均轮廓的提取。（a）下颌骨侧断面（矩形区域为初始轮廓）；（b）二值化的下颌管样本；（c）下颌管的平均轮廓（深色线）以及训练样本（浅色线）；（d）为（c）的中心区域放大效果Fig.3 Extraction fo the averagecontour of the mandibular canal.（a）The mandibular cross sectional image（Rectangular area is the initial contour）；（b）Binary sampleofthemandibularcanal；（c）The average contour of the mandibular canal（dark line）and training shapes（light line）；（d）The enlarged image of figure（c）'s center region

3 比较与讨论

为说明局部信息约束的形状导向水平集算法的优越性，在此通过展示不加局部信息的分割结果与上述分割结果进行对比。图4和图5中的（k）～（o）分别为（f）～（j）的不加局部信息（即常规形状导向水平集算法）的分割结果，其中分割参数及初始轮廓保持不变。

图4 第1组患者的下颌管分割中活动轮廓（实线）和先验形状（虚线）。（a）～（e）分别为图2中位置1～5处的下颌骨侧断面；（f）～（j）分别对应（a）～（e）的局部放大图像；（k）～（o）分别对应（f）～（j）的不加局部信息的分割结果Fig.4 The first patient’s segmentation contour（solid line）and shape prior（dotted line）.（a）～（e）The mandibular cross-sections at the position 1～5 in Fig.2，respectively；（f）～（j）The local enlarged images corresponding to（a）～（e），respectively；（k）～（o）The segmentation results without local information corresponding to（f）～（j），respectively

图5 第2组患者的下颌管分割中活动轮廓（实线）和先验形状（虚线）。（a）～（e）分别为图2中位置1～5处的下颌骨侧断面；（f）～（j）分别对应（a）～（e）的局部放大图像；（k）～（o）分别对应（f）～（j）的不加局部信息的分割结果Fig.5 The second patient’s segmentation results（contour（solid line）and shape prior（dotted line））.（a）～（e）The mandibular cross-sections at the position 1～5 in Fig.2，respectively；（f）～（j）The local enlarged images corresponding to（a）～（e），respectively；（k）～（o）The segmentation results without local information corresponding to（f）～（j），respectively

这里以图4中的（f）和（k）为例进行对比说明。从图4（f）中可以看出，由于加入了局部信息，实线表示的活动轮廓从演化到接近虚线表示的下颌管边缘时，克服了较模糊的周围组织低对比度的影响，最终在下颌管边缘处停止了演化，勾勒出下颌管的边缘。在图4（k）中，由于未加局部信息，实线表示的活动轮廓在虚线表示的下颌管边缘外侧便停止了演化，未能有效逼近真实的下颌管边缘，分割效果较差。图5中对应图像的分割效果与上述分析结果一致。为进一步对比分割效果，按文献［14］方法，将算法分割结果与参考图像（口腔医师对图4和图5中（a）～（e）的10幅图像手动分割下颌管结果）的相对变化率作为分割精度Mf进行评价。Mf=Rf－Sf，其中Rf表示从参考图像获得的特征值，Sf为从算法分割图像获得的特征值，Rf为Rf的度量值，Rf－Sf表示参考图像特征和算法分割图像特征差异的度量值。一般Rf、Sf是感兴趣目标的像素个数，Rf－Sf是两种分割结果中不同像素的个数。Mf值越小，说明分割的质量越高，表1是利用本算法和形状导向的水平集算法对感兴趣区域的分割精度Mf的对比结果。

从表1中可以看出，本算法的平均分割精度为1.82％，而形状导向水平集算法的平均分割精度为4.19％，由此可见，本算法由于其速度函数融合了图像的局部信息、区域信息、边缘梯度信息和先验形状信息，因此它适应待分割模糊边界图像的能力大大增强，且对初始轮廓要求不高，分割结果较理想。

表1 图像分割精度对比结果Tab.1 The accuracy comparison results of image segmentation

4 结论

在下颌管分割时，常规的形状导向水平集算法由于受到待分割区域模糊边缘的影响，不能有效地提高分割准确性。针对上述问题，本研究引入下颌管图像的局部信息，将局部、区域、边界和形状信息进行融合，重新定义水平集函数，使其能够克服形状导向水平集算法的缺陷，以解决目标边界区域微弱变化的下颌管图像分割效果较差的问题。实验结果证明，所提出的局部信息约束的形状导向水平集算法在下颌管的分割中取得了很好的效果，为临床上下颌管的定位提供了一种行之有效的新方法。

［1］Pogrel MA.Damage to the inferior alveolar nerve as the result of root canal therapy［J］.JADA，2007，138（1）：65 －69.

［2］Gyehyun K，Jeongjin L，Ho L，et al.Automatic extraction of inferior alveolar nerve canal using feature-enhancing panoramic volume rendering［J］.IEEE Trans on Technol Biomed，2011，2（58）：253－246.

［3］Tamas F. Position ofthe mandibularcanal［J］. JOral Maxillofac Surg，1987，16：65－69.

［4］Obradovic O，Todorovic L.Morphometric analysis of mandibular canal：clinic aspects［J］.Bull Group Int Rech Sci Stomatol Odontol，1993，36：109 －113.

［5］李宁毅，马彦博，谷芳，等.数字化可视下颌骨的初步研究［J］.华西口腔医学杂志，2007，25（1）：83－86.

［6］Pu Tiemin， Zhu Xinyi， LiJianfeng， etal. Multi-section measurements of mandibular canal of adult，ex vivo mandibles［J］.中国组织工程研究与临床康复，2009，13（48）：9592－9596.

［7］Bu Lingxue，Wang Ke，Chen Xin，et al.Anatomic structure of the mandibular canal［J］.中国组织工程研究与临床康复，2011，15（2）：377－380.

［8］王吉昌，归来，张智勇，等.汉族青年女性下颌神经管的三维CT定位测量［J］.中华整形外科杂志，2007，23（3）：212－214.

［9］Yau HT，Lin YK，Tsou LS，et al.An adaptive region growing method to segment inferior alveolar nerve canal from 3D medical images for dental implant surgery［J］.Computer-Aided Design and Applications，2008，5（5）：743－752.

［10］Bresson X，Vandergheynst P，Thiran JP.A variational model for object segmentation using boundary information and shape prior driven by the mumford-shah functional［J］.International Journal of Computer Vision，2006，68（2）：145－162.

［11］Leventon M E，Grimson WEL，Faugers Q.Statistical shape influence in geodesic active contour［C］//IEEE Conference of Computer Vision and Patten Recognition.Hilton Head Island：IEEE，2000：1：316 －323.

［12］Li Chunming，Kao CY，Gore JC，et al.Implicit active contours driven by local binary fitting energy［C］ //IEEE Conference on Computer Vision and Pattern Recognition.Minneapolis：IEEE，2007：1－7.

［13］Wang Li，Li Chunming，Sun Quansen，et al.Brain MR image segmentation using localand globalintensity fitting active contours/surfaces［C］//11th InternationalConference on Medical Image Computing and Computer-Assisted Intervention.New York：Springer，2008：5241：384－392.

［14］Pardo XM，Carreira MJ，Mosquera A，et al.A snake for CT image segmentation integrating region and edge information［J］.Image and Vision Computing，2001，19（7）：461－475.