一种基于混合水平集的三维牙齿建模方法
一种基于混合水平集的三维牙齿建模方法
【技术领域】
[0001]本发明属于图像处理领域，涉及一种水平集方法在三维牙齿建模中的应用。
【背景技术】
[0002] 近年来，随着三维数字化成像技术在口腔医学领域的飞速发展，计算机辅助诊疗 修复越来越多的应用在口腔修复当中，并逐渐成为这一领域的发展趋势。在计算机口腔修 复系统中，首先重要的就是要获取数字化的三维牙齿模型，而模型的精度和完整性直接关 系到后续的排牙、种植、正畸、以及生物力学分析的结果。
[0003] 目前牙齿建模方法最常用的方法就是利用图像处理技术分割牙齿CT图像序列，即 首先从每一层CT切片中分割出牙齿轮廓，然后利用这些层间轮廓重建出牙齿三维模型。由 于该类方法能够获得整个牙齿形状结构，为患者口腔病变提供完整的诊断依据，因此，基于 CT图像的牙齿建模方法越来越受到研究学者的广泛关注。
[0004] 由于口腔CT图像中牙齿和颂骨的密度和距离都较为接近，采用传统的图像分割方 法很难精确地提取出每颗牙齿的组织轮廓。CT图像牙齿分割一直是一个充满挑战性的课 题。王黎等（王黎，崔进，韩清凯等，基于CT图像的牙齿3维实体模型建立。中国图象图形学 报，2005,10(10): 1289-1292。）利用二值化和边界提取筛选出每层切片牙齿轮廓关键点，然 后利用3D-Delaunay四面体化算法得到整颗牙齿的实体模型；但二值化操作容易产生过分 割或欠分割的问题。Wu 等（X.Wu，H.Gao，H.Heo，etal.ImprovedB-splinecontour fitting using genetic algorithm for the segmentation of dental computerized tomography image sequences. The Journal of imaging science and technology, 2007，51 (4): 328-336.)利用基于遗传算法的B样条曲线拟合来提取每层切片的牙齿轮廓， 但B样条曲线无法处理牙齿拓扑结构变化的问题。Gao等(H.Gao，0.Chae · Individual tooth segmentation from CT image s using level set method with shape and intensity prior .Pattern Recognit ion，2010,43:2406-2417)采用基于水平集的活动轮廓模型进行 牙齿分割，并对不同牙层切片采用不同的分割模型，能够将每颗牙齿的牙冠和牙根都分割 出来;该模型主要依靠图像边缘梯度和先验形状的概率分布来指导水平集的演化，但由于 牙齿密度不均匀，且周围容易受到牙槽骨等结构的干扰，因此依靠先验形状周围区域的灰 度概率来控制水平集轮廓的收缩和扩张容易产生边界泄露的问题。
【发明内容】
[0005] 为克服上述现有技术的不足，提高牙齿分割的效率和精度，本发明综合考虑口腔 CT图像牙齿形状变化特点，提出一种基于混合水平集的三维牙齿建模方法。
[0006] 本发明提供:一种基于混合水平集的三维牙齿建模方法，包括如下步骤：
[0007] (1)从口腔CT图像序列中选取一张作为起始切片，并在该图像上勾画出每颗牙齿 轮廓以初始化水平集函数；
[0008] (2)起始切片以下的牙根层切片利用先验形状约束能量、基于Flux模型的边缘能 量、基于先验灰度的局部区域能量相结合构造的单相混合水平集模型分割牙齿轮廓；
[0009] (3)起始切片以上的牙冠层切片利用结合区域竞争约束的双相混合水平集模型分 割牙齿轮廓；
[0010] (4)将所有切片分割后的牙齿轮廓像素点转化为三维坐标，利用Delaunay三角剖 分方法进行重建以得到每颗牙齿的三角网格模型。
[0011] 步骤（1)按如下步骤进行:在牙颈部位的切片图像中选取一张所有牙都出现且牙 槽骨较少出现的切片作为起始切片，并在该切片图像上勾画出每颗牙齿的大致轮廓C 1G = l，2，...n，n为牙齿个数)作为水平集的初始轮廓，然后利用C1初始化η个水平集函数ΦΚ? = 1，2，. . .η)，Φ i的初始化通过计算图像上每个点到匕的带符号的距离来完成，即：
[0013]其中，d[(x)，Ci]表示像素点X与曲线Ci之间的欧式距离。
[0014] 步骤(2)所述的分割牙根层切片的单相混合水平集模型的能量泛函定义为先验形 状约束能量、基于Flux模型的边缘能量、基于先验灰度的局部区域能量等的加权和：
[0015] E3fg( Φ ) 一yEint( Φ )+ 7" Elength( Φ )+ClEprior ( Φ )+vEedge( Φ )+^EregionC Φ )
[0016]其中，μ，γ，α，ν，λ为各个能量项的权系数；
[0017] (3a)符号距离保持能量Eint(O)，用来保证水平集演化过程中的稳定性，定义为：
[0019] (3b)曲线弧长平滑能量Elength(O)，用来平滑水平集轮廓，定义为：
[0021] (3c)先验形状约束能量ΕΡ"α(Φ)，用来控制水平集的形状，将每次分割后的牙齿 轮廓映射到相邻切片图像，作为当前水平集函数演化的先验形状加以约束，其能量泛函定 义为：
[0023] 其中Φ为当前切片的水平集函数，ΦΡ为上一张切片分割完成后先验形状对应的 水平集函数，H(X)为Heaviside函数；
[0024] (3d)基于Flux模型的边缘能量Ee3dge3(O),用来探测牙齿的外边界轮廓，将图像梯 度方向和水平函数梯度方向之间的角度信息嵌入到传统Flux模型当中，其能量泛函定义 为：
[0026]其中Δ为Laplacian算子，I。代表高斯平滑后的图像，
[0028] 为梯度方向检测函数，▽为梯度算子，?代表点积；
[0029] (3e)基于先验灰度的局部区域能量Ereglcin(O),用来克服图像灰度不均匀问题，将 先验灰度信息嵌入到区域模型当中，其能量泛函定义为：
[0031] 其中fre3f_in(x)和frrfjtU)分别定义为先验形状上参考点的r邻域在先验形状曲 线内、外的灰度均值，参考点为先验形状上距离当前图像目标像素点X最近的点；
[0032] (3f)综合以上各能量项，将牙根层切片混合水平集模型的最小化能量泛函表示 为：
[0034]对上述能量泛函进行最小化，得到水平集曲线的演化方程：
[0036] 其中，δ(Φ)为Dirac 函数。
[0037] 步骤(3)所述中的分割牙冠层切片的双相混合水平集模型按如下步骤建立:对起 始切片利用所述的单相混合水平集模型进行分割，将分割后的所有水平集函数按照每隔一 个牙齿的规则将其组合成两个耦合的双相混合水平集函数，并加入区域竞争约束能量以克 服两个水平集产生重叠，其能量泛函定义为：
[0038] Eas( Φ 1 , Φ 2 ) = Φ I ) +Eafg( Φ 2 ) +PErepulse
[0039] 其中，Erepulse = jQH( Φ?)Η( C>2)dx
[0040] 为区域竞争约束能量，β用来控制两水平集函数区域重叠的程度；
[0041]对上述能量泛函进行最小化，得到双相水平集函数? 的演化方程分别为：
[0044] 其中Φρ?，ΦΡ2分别代表Φι，φ>2的先验值;ξ?(Χ)、ξ2(Χ)分别代表Φι，φ>2演化时的 梯度方向检测函数：
[0047]步骤(4)所述的三维牙齿重建按如下步骤进行：利用窄带法和半隐式差分方案求 解上 述水平集函数的演化方程，提取更新完成后的水平集函数在φ=0的像素点，得到当前 切片图像的牙齿轮廓，将所有切片图像分割后的牙齿轮廓像素点转化为三维坐标，利用基 于区域增长的Delaunay三角剖分方法进行重建，以生成每颗牙齿的三角网格模型。
[0048]本发明的有益效果:本发明提出的混合水平集模型结合了图像边缘、局部区域、先 验知识等多方面的信息，能有效克服传统模型边缘定位不准确以及无法处理图像灰度不均 匀等问题。本发明方法具有较少的人工干预，而且具有较好的分割效果和较高的准确率，重 建出的三维牙齿模型能正确反映牙齿的解剖形态，从而为下一步制定口腔修复规划、生物 力学分析等奠定坚实的基础，对于提高口腔修复的精度和效率具有重要意义。
【附图说明】
[0049]图1为本发明的牙齿建模技术流程图。
[0050] 图2为起始切片的水平集函数初始化示意图。
[0051] 图3(a)为牙齿与周围干扰结构的梯度方向示意图。
[0052]图3(b)为水平集函数的梯度方向示意图。
[0053]图4为牙冠层切片图像的双相水平集函数示意图。
[0054]图5为双相水平集分割区域示意图。
[0055] 图6为本发明与现有方法对起始切片图像的分割效果图。
[0056] 图7为本发明与现有方法对牙根切片磨牙部位的分割效果图。
[0057]图8为本发明磨牙部位若干切片分割轮廓和三维重建效果图。
【具体实施方式】
[0058]下面结合附图对本发明实施例作进一步说明：
[0059]如图1所示，本发明的具体实现步骤如下：
[0060]步骤1、起始切片的选取和水平集初始化:首先读取牙齿CT图像序列，然后在牙颈 部位的切片图像中选择一张所有牙都出现且牙槽骨较少出现的切片作为起始切片。在起始 切片图像上，勾画出每颗牙齿的大致轮廓(： 1(1 = 1，2，...11，11为牙齿个数)作为水平集的初 始轮廓，如图2所示，然后利用Ci(i = l，2, . . .η)初始化η个水平集函数C>i(i = l，2, . . .η)。 Φ:的初始化通过计算图像上每个点到匕的带符号的距离来完成，即：
[0062]其中，d[(x)，Ci]表示点X与曲线Ci之间的欧式距离。
[0063] 步骤2、牙根层切片分割:在牙根层图像中，牙齿易受到牙槽骨和牙髓腔等结构的 较多干扰。因此，提出一种新的混合水平集牙齿分割模型。该模型的能量泛函主要由先验形 状约束能量、基于Flux模型的边缘能量、基于先验灰度的局部区域能量等部分组成。
[0064] (1)在先验形状约束能量项中，由于相邻切片之间，对应牙齿的形状和位置都比较 接近，因此将上一张切片已完成分割的牙齿轮廓映射到当前切片图像，作为水平集演化的 先验形状加以约束。令Φ为当前切片的水平集函数，Φ Ρ为上一张切片分割完成后先验形状 对应的水平集函数，则先验形状约束项定义为：
[0066] 式中，H(x)为Heaviside 函数。
[0067] (2)在基于Flux模型的边缘能量中，利用基于梯度向量场的Flux模型来定位目标 边界，同时利用图像梯度方向和水平函数梯度方向之间的角度来克服周围牙槽骨和牙髓腔 的干扰问题。
[0068]由于图像在边界处梯度与边缘法向量同向，因此Flux模型通过计算图像梯度与曲 线法向量内积的边缘积分来定位图像边缘，该方法能有效解决弱边界目标的分割问题，且 演化速度快。Flux模型的最大化能量泛函表示为：
[0069] Efiux( Φ) = ΙωΔΙ( 1-Η( Φ )dx
[0070] 式中，Δ为Laplacian算子，I。代表高斯平滑后的图像。
[0071]然而，由于牙齿内部和外部的干扰结构较多，直接利用上述模型进行分割，水平集 会同时捕捉到牙齿内部的牙髓腔以及外部的邻牙或牙槽骨等边界上，因此本发明将梯度方 向约束加入到上述的Flux模型当中。
[0072]在牙齿CT图像中，牙齿为高亮度区域，背景为低亮度区域，因此可将牙齿与周围干 扰结构的图像梯度方向表示为图3a所示。由于本文定义的水平集函数Φ在水平集内部取 正、外部取负，则Φ的梯度方向可表示为图3b所示。可以看出，如果要使演化曲线只捕捉到 牙齿的外边界轮廓，则图像的梯度方向应该与水平集函数的梯度方向一致。因此，将本发明 的边缘能量的最小化泛函定义为：
[0075] 为梯度方向检测函数，▽为梯度算子，?代表点积。
[0076] (3)在结合先验灰度的局部区域能量中，将先验形状的局部邻域灰度均值取代全 局均值，以克服传统模型无法利用先验灰度信息和无法处理图像灰度不均匀的问题。
[0077] 由于相邻切片之间图像的灰度分布非常相同，且牙齿轮廓位置也接近，因此本发 明将先验灰度信息嵌入到能量模型中，以提高比较精度，其能量泛函定义为：
[0079] 其中fre3f_in(x)和frrfjtU)分别定义为先验形状上参考点的r邻域在先验形状曲 线内、外的灰度均值。参考点按照如下方法确定：当上一张切片分割完成后，计算先验形状 上每个点的r邻域在先验形状曲线内、外的灰度均值，然后在当前切片图像中，将先验形状 上距离当前图像目标像素点X最近的点作为参考点。r邻域半径的选取需根据图像的分辨率 而定，一般不应太大，以免包含过多的邻牙区域。
[0080] (4)为保证水平集在整个演化过程中的形状和稳定性，再加入符号距离保持能量：
[0082] 和曲线弧长平滑项：
σ
[0084]这样，综合以上各个能量信息，可将分割牙根层切片的单相混合水平集模型的最 小化能量泛函表示为：
[0086]对能量方程进行最小化，可得到牙根层切片的水平集函数Φ的演化方程为：
[0088] 其中 δ(Φ)为Dirac 函数。
[0089] 步骤3、牙冠层切片分割模型:在牙冠层切片图像中，牙齿逐渐紧连在一起，缝隙较 小，为有效提取相邻两牙之间的边界，对起始切片利用步骤2所述的单相混合水平集模型进 行牙齿分割，将分割后的水平集函数按照每隔一个牙齿的规则将其组合成两个耦合的双相 混合水平集函数，如图4所示。
[0090] 由多相水平集分割原理可知，两个水平集函数可以将图像划分为四个区 域，如图5所示。对于{ΦχΧ^ΦζΧ)}这个相，代表相邻牙齿的重叠区域。为克服 演化过程中产生重叠，在能量泛函中引入区域竞争准则。设A^A 2分别为水平集函数 对应的水平集曲线围成的内部区域面积。根据Heaviside函数H(X)的性质，可将A1,A2分别表 示为：
[0091] Ai = area ( Φι>0) = ?ΩΗ(Φι)?χ
[0092] A2 = area( C>2 2 0) = JqH( C>2)dx
[0093]要使这两个区域不产生重叠，则应满足：
[0095]因此可将区域惩罚项定义为：
[0096] Erepulse = jQH( Φ?)Η( ?2)dx
[0097] 对上式进行最小化相当于避免两个牙齿产生重叠。
[0098] 将区域惩罚项加入到水平集能量泛函当中，可得到牙冠层的双相混合水平集模型 的能量方程为：
[0099] Eas( Φ 1 , Φ 2 ) = Φ I ) +Egfg( Φ 2 ) +PErepulse
[0100] 其中β用来控制两水平集函数区域重叠的程度，对上述能量泛函进行最小化，可得 到双相水平集函数φ:，φ2的演化方程分别为：
[0104] 其中Φρ?，ΦΡ2分别代表Φι，φ>2的先验值;ξ?(Χ)、ξ2(Χ)分别代表Φι ，φ>2演化时的 梯度方向检测函数，BP:
[0107] 步骤4、牙齿三维模型重建:利用窄带法和半隐式差分方案求解上述水平集函数的 演化方程，提取更新完成后的水平集函数在Φ=〇的像素点，将所有切片图像分割后的这些 牙齿轮廓像素点转化为三维坐标，利用Delaunay三角剖分方法进行重建，即可得到每颗牙 齿的三角网格模型。
[0108] 本发明的有效性可以通过下面的实验进一步说明：
[0109] 以一位下颂牙齿完整的患者口腔CT图像序列作为例。图6为利用本发明和Li模型、 CV模型对该CT序列的起始切片图像进行牙齿分割结果对比图。图6(a)为在起始切片上画的 各个牙齿初始轮廓多边形。可以看出，基于边缘的Li模型由于过于依赖初始曲线附近的图 像梯度信息，但缺乏梯度方向约束，所以图6(b)中存在很多牙齿边缘定位不准确的现象。基 于区域的C-V模型依赖全局灰度信息，因此在图6(c)中距离较近的牙齿之间以及牙齿与颂 骨之间的弱边界处产生融合。图6(d)为本发明的分割结果，由于利用基于梯度向量场的 Flux模型，能够有效探测弱边界，并通过加入梯度方向和先验形状约束，从而实现准确定位 牙齿边界轮廓。
[0110] 图7(a)和图7(b)分别本发明和Gao等人的方法对该患者第二磨牙的牙根部位切片 图像的分割结果。可以看出，本发明具有较准确的分割结果，但Gao等人的方法依靠曲线内 外灰度概率来探测轮廓，由于牙髓腔的灰度较暗，会在牙齿弱边界处使探测的曲线偏向牙 髓腔。
[0111] 图8显示了在右侧第一磨牙部位的一些切片的分割结果。可以看出，基本上所有切 片图像的分割结果都能很好地贴合牙齿目标轮廓。但是，在第160层的切片中，水平集没有 捕捉到牙齿内部沟槽轮廓，这主要是受到水平集演化方程中窄带大小的影响。由于牙齿表 面沟槽的存在，会造成一些切片的轮廓出现极度的凹陷。如果此时窄带宽度取得较小，可能 就会限制水平集演化到目标轮廓上。但这并不影响最后的分割结果，由于目标轮廓线包含 在水平集曲线的内部，可通过灰度阈值提取出最终的轮廓线。图8右边显示了对分割后的每 层切片牙齿轮廓进行三维重建后的结果，可以看出本发明重建的三维模型能正确反映牙齿 的解剖形态，从而为下一步的计算机辅助诊断提供有效数据，使医生能够针对病人的实际 情况制定最佳的手术方案。
[0112]以上所述仅为本发明的较佳实例而已，并不用以限制本发明，凡在本发明的精神 和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。
【主权项】
1. 一种基于混合水平集的三维牙齿建模方法，其特征在于，包括如下步骤： (1) 从口腔CT图像序列中选取一张作为起始切片，并在该图像上勾画出每颗牙齿轮廓 以初始化水平集函数； (2) 起始切片以下的牙根层切片利用先验形状约束能量、基于Flux模型的边缘能量、基 于先验灰度的局部区域能量相结合构造的单相混合水平集模型分割牙齿轮廓； (3) 起始切片以上的牙冠层切片利用结合区域竞争约束的双相混合水平集模型分割牙 齿轮廓； (4) 将所有切片分割后的牙齿轮廓像素点转化为三维坐标，利用Delaunay三角剖分方 法进行重建以得到每颗牙齿的三角网格模型。2. 根据权利要求1所述的一种基于混合水平集的三维牙齿建模方法，其特征在于，步骤 (1) 按如下步骤进行:在牙颈部位的切片图像中选取一张所有牙都出现且牙槽骨较少出现 的切片作为起始切片，并在该切片图像上勾画出每颗牙齿的大致轮廓C 1Q = I, 2，···η，η为 牙齿个数)作为水平集的初始轮廓，然后利用C1初始化η个水平集函数O1Q = Id, ···!〇 ,CD1 的初始化通过计算图像上每个点到匕的带符号的距离来完成，即：其中，d [(X )，Ci ]表示像素点X与曲线Ci之间的欧式距离。3. 根据权利要求1所述的一种基于混合水平集的三维牙齿建模方法，其特征在于，步骤 (2) 所述的分割牙根层切片的单相混合水平集模型的能量泛函定义为先验形状约束能量、 基于Flux模型的边缘能量、基于先验灰度的局部区域能量等的加权和：I其中，μ，γ，α， ν，λ为各个能量项的权系数； (3a)符号距离保持能量Ειη?(Φ)，用来保证水平集演化过程中的稳定性，定义为：(3b)曲线弧长平滑能量Ele3ngth( Φ )，用来平滑水平集轮廓，定义为：(3c)先验形状约束能量ΕΡ"α(Φ)，用来控制水平集的形状，将每次分割后的牙齿轮廓 映射到相邻切片图像，作为当前水平集函数演化的先验形状加以约束，其能量泛函定义为：其中Φ为当前切片的水平集函数，Φρ为上一张切片分割完成后先验形状对应的水平集 函数，H(X)为Heaviside函数； (3d)基于Flux模型的边缘能量心以〇)，用来探测牙齿的外边界轮廓，将图像梯度方向 和水平函数梯度方向之间的角度信息嵌入到传统Flux模型当中，其能量泛函定义为： Eedge( Φ ) =-?ΩξΔΙσ(?-Η( Φ ) dx 其中Λ为Laplacian算子，I。代表高斯平滑后的图像，为梯度方向检测函数，▽:为梯度算子，?代表点积； (3e)基于先验灰度的局部区域能量Ereglcin( Φ )，用来克服图像灰度不均匀问题，将先验 灰度信息嵌入到区域模型当中，其能量泛函定义为：其中fre3f_ln(X)和frrf+cmtU)分别定义为先验形状上参考点的r邻域在先验形状曲线内、 外的灰度均值，参考点为先验形状上距离当前图像目标像素点X最近的点； (3f)综合以上各能量项，将牙根层切片混合水平集模型的最小化能量泛函表示为：对上述能量泛函进行最小化，得到水平集曲线的演化方程：其中，δ( Φ )为Dirac函数。4.根据权利要求3所述的一种基于混合水平集的三维牙齿建模方法，其特征在于，步骤 (3)所述中的分割牙冠层切片的双相混合水平集模型按如下步骤建立:对起始切片利用所 述的单相混合水平集模型进行分割，将分割后的所有水平集函数按照每隔一个牙齿的规则 将其组合成两个耦合的双相混合水平集函数，并加入区域竞争约束能量以克服两个水平集 产生重叠，其能量泛函定义为：为区域竞争约束能量，β用来控制两水平集函数区域重叠的程度； 对上述能量泛函进行最小化，得到双相水平集函数φ:，φ2的演化方程分别为：其中Φρ1，ΦΡ2分别代表Φ?，Φ2的先验值;ξ?(Χ)、ξ2(Χ)分别代表Φ?，Φ2演化时的梯度方 向检测函数：05.根据权利要求1所述的一种基于混合水平集的三维牙齿建模方法，其特征在于，步骤 (4)所述的三维牙齿重建按如下步骤进行:利用窄带法和半隐式差分方案求解上述水平集 函数的演化方程，提取更新完成后的水平集函数在Φ =O的像素点，得到当前切片图像的牙 齿轮廓，将所有切片图像分割后的牙齿轮廓像素点转化为三维坐标，利用基于区域增长的 Delaunay三角剖分方法进行重建，以生成每颗牙齿的三角网格模型。
【专利摘要】本发明提供一种基于混合水平集的三维牙齿建模方法，包括以下几个步骤：1)起始切片的选取和水平集初始化；2)利用先验形状约束能量、基于Flux模型的边缘能量、基于先验灰度的局部区域能量相结合构造的单相混合水平集模型分割牙根层切片图像；3)结合区域竞争约束的双相混合水平集模型分割牙冠层切片图像；4)基于分割轮廓的牙齿三维模型重建。本发明的混合水平集模型能有效克服传统分割模型边缘定位不准确以及无法处理图像灰度不均匀等问题，快速准确地构建出每颗牙齿独立的三维模型，从而为制定口腔修复规划、生物力学分析等奠定坚实的基础。
【IPC分类】G06T17/30
【公开号】CN105488849
【申请号】CN201510821682
【发明人】吴婷, 张礼兵
【申请人】嘉兴学院
【公开日】2016年4月13日
【申请日】2015年11月24日
【技术领域】
[0001]本发明属于图像处理领域，涉及一种水平集方法在三维牙齿建模中的应用。
【背景技术】
[0002] 近年来，随着三维数字化成像技术在口腔医学领域的飞速发展，计算机辅助诊疗 修复越来越多的应用在口腔修复当中，并逐渐成为这一领域的发展趋势。在计算机口腔修 复系统中，首先重要的就是要获取数字化的三维牙齿模型，而模型的精度和完整性直接关 系到后续的排牙、种植、正畸、以及生物力学分析的结果。
[0003] 目前牙齿建模方法最常用的方法就是利用图像处理技术分割牙齿CT图像序列，即 首先从每一层CT切片中分割出牙齿轮廓，然后利用这些层间轮廓重建出牙齿三维模型。由 于该类方法能够获得整个牙齿形状结构，为患者口腔病变提供完整的诊断依据，因此，基于 CT图像的牙齿建模方法越来越受到研究学者的广泛关注。
[0004] 由于口腔CT图像中牙齿和颂骨的密度和距离都较为接近，采用传统的图像分割方 法很难精确地提取出每颗牙齿的组织轮廓。CT图像牙齿分割一直是一个充满挑战性的课 题。王黎等（王黎，崔进，韩清凯等，基于CT图像的牙齿3维实体模型建立。中国图象图形学 报，2005,10(10): 1289-1292。）利用二值化和边界提取筛选出每层切片牙齿轮廓关键点，然 后利用3D-Delaunay四面体化算法得到整颗牙齿的实体模型；但二值化操作容易产生过分 割或欠分割的问题。Wu 等（X.Wu，H.Gao，H.Heo，etal.ImprovedB-splinecontour fitting using genetic algorithm for the segmentation of dental computerized tomography image sequences. The Journal of imaging science and technology, 2007，51 (4): 328-336.)利用基于遗传算法的B样条曲线拟合来提取每层切片的牙齿轮廓， 但B样条曲线无法处理牙齿拓扑结构变化的问题。Gao等(H.Gao，0.Chae · Individual tooth segmentation from CT image s using level set method with shape and intensity prior .Pattern Recognit ion，2010,43:2406-2417)采用基于水平集的活动轮廓模型进行 牙齿分割，并对不同牙层切片采用不同的分割模型，能够将每颗牙齿的牙冠和牙根都分割 出来;该模型主要依靠图像边缘梯度和先验形状的概率分布来指导水平集的演化，但由于 牙齿密度不均匀，且周围容易受到牙槽骨等结构的干扰，因此依靠先验形状周围区域的灰 度概率来控制水平集轮廓的收缩和扩张容易产生边界泄露的问题。
【发明内容】
[0005] 为克服上述现有技术的不足，提高牙齿分割的效率和精度，本发明综合考虑口腔 CT图像牙齿形状变化特点，提出一种基于混合水平集的三维牙齿建模方法。
[0006] 本发明提供:一种基于混合水平集的三维牙齿建模方法，包括如下步骤：
[0007] (1)从口腔CT图像序列中选取一张作为起始切片，并在该图像上勾画出每颗牙齿 轮廓以初始化水平集函数；
[0008] (2)起始切片以下的牙根层切片利用先验形状约束能量、基于Flux模型的边缘能 量、基于先验灰度的局部区域能量相结合构造的单相混合水平集模型分割牙齿轮廓；
[0009] (3)起始切片以上的牙冠层切片利用结合区域竞争约束的双相混合水平集模型分 割牙齿轮廓；
[0010] (4)将所有切片分割后的牙齿轮廓像素点转化为三维坐标，利用Delaunay三角剖 分方法进行重建以得到每颗牙齿的三角网格模型。
[0011] 步骤（1)按如下步骤进行:在牙颈部位的切片图像中选取一张所有牙都出现且牙 槽骨较少出现的切片作为起始切片，并在该切片图像上勾画出每颗牙齿的大致轮廓C 1G = l，2，...n，n为牙齿个数)作为水平集的初始轮廓，然后利用C1初始化η个水平集函数ΦΚ? = 1，2，. . .η)，Φ i的初始化通过计算图像上每个点到匕的带符号的距离来完成，即：
[0013]其中，d[(x)，Ci]表示像素点X与曲线Ci之间的欧式距离。
[0014] 步骤(2)所述的分割牙根层切片的单相混合水平集模型的能量泛函定义为先验形 状约束能量、基于Flux模型的边缘能量、基于先验灰度的局部区域能量等的加权和：
[0015] E3fg( Φ ) 一yEint( Φ )+ 7" Elength( Φ )+ClEprior ( Φ )+vEedge( Φ )+^EregionC Φ )
[0016]其中，μ，γ，α，ν，λ为各个能量项的权系数；
[0017] (3a)符号距离保持能量Eint(O)，用来保证水平集演化过程中的稳定性，定义为：
[0019] (3b)曲线弧长平滑能量Elength(O)，用来平滑水平集轮廓，定义为：
[0021] (3c)先验形状约束能量ΕΡ"α(Φ)，用来控制水平集的形状，将每次分割后的牙齿 轮廓映射到相邻切片图像，作为当前水平集函数演化的先验形状加以约束，其能量泛函定 义为：
[0023] 其中Φ为当前切片的水平集函数，ΦΡ为上一张切片分割完成后先验形状对应的 水平集函数，H(X)为Heaviside函数；
[0024] (3d)基于Flux模型的边缘能量Ee3dge3(O),用来探测牙齿的外边界轮廓，将图像梯 度方向和水平函数梯度方向之间的角度信息嵌入到传统Flux模型当中，其能量泛函定义 为：
[0026]其中Δ为Laplacian算子，I。代表高斯平滑后的图像，
[0028] 为梯度方向检测函数，▽为梯度算子，?代表点积；
[0029] (3e)基于先验灰度的局部区域能量Ereglcin(O),用来克服图像灰度不均匀问题，将 先验灰度信息嵌入到区域模型当中，其能量泛函定义为：
[0031] 其中fre3f_in(x)和frrfjtU)分别定义为先验形状上参考点的r邻域在先验形状曲 线内、外的灰度均值，参考点为先验形状上距离当前图像目标像素点X最近的点；
[0032] (3f)综合以上各能量项，将牙根层切片混合水平集模型的最小化能量泛函表示 为：
[0034]对上述能量泛函进行最小化，得到水平集曲线的演化方程：
[0036] 其中，δ(Φ)为Dirac 函数。
[0037] 步骤(3)所述中的分割牙冠层切片的双相混合水平集模型按如下步骤建立:对起 始切片利用所述的单相混合水平集模型进行分割，将分割后的所有水平集函数按照每隔一 个牙齿的规则将其组合成两个耦合的双相混合水平集函数，并加入区域竞争约束能量以克 服两个水平集产生重叠，其能量泛函定义为：
[0038] Eas( Φ 1 , Φ 2 ) = Φ I ) +Eafg( Φ 2 ) +PErepulse
[0039] 其中，Erepulse = jQH( Φ?)Η( C>2)dx
[0040] 为区域竞争约束能量，β用来控制两水平集函数区域重叠的程度；
[0041]对上述能量泛函进行最小化，得到双相水平集函数? 的演化方程分别为：
[0044] 其中Φρ?，ΦΡ2分别代表Φι，φ>2的先验值;ξ?(Χ)、ξ2(Χ)分别代表Φι，φ>2演化时的 梯度方向检测函数：
[0047]步骤(4)所述的三维牙齿重建按如下步骤进行：利用窄带法和半隐式差分方案求 解上 述水平集函数的演化方程，提取更新完成后的水平集函数在φ=0的像素点，得到当前 切片图像的牙齿轮廓，将所有切片图像分割后的牙齿轮廓像素点转化为三维坐标，利用基 于区域增长的Delaunay三角剖分方法进行重建，以生成每颗牙齿的三角网格模型。
[0048]本发明的有益效果:本发明提出的混合水平集模型结合了图像边缘、局部区域、先 验知识等多方面的信息，能有效克服传统模型边缘定位不准确以及无法处理图像灰度不均 匀等问题。本发明方法具有较少的人工干预，而且具有较好的分割效果和较高的准确率，重 建出的三维牙齿模型能正确反映牙齿的解剖形态，从而为下一步制定口腔修复规划、生物 力学分析等奠定坚实的基础，对于提高口腔修复的精度和效率具有重要意义。
【附图说明】
[0049]图1为本发明的牙齿建模技术流程图。
[0050] 图2为起始切片的水平集函数初始化示意图。
[0051] 图3(a)为牙齿与周围干扰结构的梯度方向示意图。
[0052]图3(b)为水平集函数的梯度方向示意图。
[0053]图4为牙冠层切片图像的双相水平集函数示意图。
[0054]图5为双相水平集分割区域示意图。
[0055] 图6为本发明与现有方法对起始切片图像的分割效果图。
[0056] 图7为本发明与现有方法对牙根切片磨牙部位的分割效果图。
[0057]图8为本发明磨牙部位若干切片分割轮廓和三维重建效果图。
【具体实施方式】
[0058]下面结合附图对本发明实施例作进一步说明：
[0059]如图1所示，本发明的具体实现步骤如下：
[0060]步骤1、起始切片的选取和水平集初始化:首先读取牙齿CT图像序列，然后在牙颈 部位的切片图像中选择一张所有牙都出现且牙槽骨较少出现的切片作为起始切片。在起始 切片图像上，勾画出每颗牙齿的大致轮廓(： 1(1 = 1，2，...11，11为牙齿个数)作为水平集的初 始轮廓，如图2所示，然后利用Ci(i = l，2, . . .η)初始化η个水平集函数C>i(i = l，2, . . .η)。 Φ:的初始化通过计算图像上每个点到匕的带符号的距离来完成，即：
[0062]其中，d[(x)，Ci]表示点X与曲线Ci之间的欧式距离。
[0063] 步骤2、牙根层切片分割:在牙根层图像中，牙齿易受到牙槽骨和牙髓腔等结构的 较多干扰。因此，提出一种新的混合水平集牙齿分割模型。该模型的能量泛函主要由先验形 状约束能量、基于Flux模型的边缘能量、基于先验灰度的局部区域能量等部分组成。
[0064] (1)在先验形状约束能量项中，由于相邻切片之间，对应牙齿的形状和位置都比较 接近，因此将上一张切片已完成分割的牙齿轮廓映射到当前切片图像，作为水平集演化的 先验形状加以约束。令Φ为当前切片的水平集函数，Φ Ρ为上一张切片分割完成后先验形状 对应的水平集函数，则先验形状约束项定义为：
[0066] 式中，H(x)为Heaviside 函数。
[0067] (2)在基于Flux模型的边缘能量中，利用基于梯度向量场的Flux模型来定位目标 边界，同时利用图像梯度方向和水平函数梯度方向之间的角度来克服周围牙槽骨和牙髓腔 的干扰问题。
[0068]由于图像在边界处梯度与边缘法向量同向，因此Flux模型通过计算图像梯度与曲 线法向量内积的边缘积分来定位图像边缘，该方法能有效解决弱边界目标的分割问题，且 演化速度快。Flux模型的最大化能量泛函表示为：
[0069] Efiux( Φ) = ΙωΔΙ( 1-Η( Φ )dx
[0070] 式中，Δ为Laplacian算子，I。代表高斯平滑后的图像。
[0071]然而，由于牙齿内部和外部的干扰结构较多，直接利用上述模型进行分割，水平集 会同时捕捉到牙齿内部的牙髓腔以及外部的邻牙或牙槽骨等边界上，因此本发明将梯度方 向约束加入到上述的Flux模型当中。
[0072]在牙齿CT图像中，牙齿为高亮度区域，背景为低亮度区域，因此可将牙齿与周围干 扰结构的图像梯度方向表示为图3a所示。由于本文定义的水平集函数Φ在水平集内部取 正、外部取负，则Φ的梯度方向可表示为图3b所示。可以看出，如果要使演化曲线只捕捉到 牙齿的外边界轮廓，则图像的梯度方向应该与水平集函数的梯度方向一致。因此，将本发明 的边缘能量的最小化泛函定义为：
[0075] 为梯度方向检测函数，▽为梯度算子，?代表点积。
[0076] (3)在结合先验灰度的局部区域能量中，将先验形状的局部邻域灰度均值取代全 局均值，以克服传统模型无法利用先验灰度信息和无法处理图像灰度不均匀的问题。
[0077] 由于相邻切片之间图像的灰度分布非常相同，且牙齿轮廓位置也接近，因此本发 明将先验灰度信息嵌入到能量模型中，以提高比较精度，其能量泛函定义为：
[0079] 其中fre3f_in(x)和frrfjtU)分别定义为先验形状上参考点的r邻域在先验形状曲 线内、外的灰度均值。参考点按照如下方法确定：当上一张切片分割完成后，计算先验形状 上每个点的r邻域在先验形状曲线内、外的灰度均值，然后在当前切片图像中，将先验形状 上距离当前图像目标像素点X最近的点作为参考点。r邻域半径的选取需根据图像的分辨率 而定，一般不应太大，以免包含过多的邻牙区域。
[0080] (4)为保证水平集在整个演化过程中的形状和稳定性，再加入符号距离保持能量：
[0082] 和曲线弧长平滑项：
σ
[0084]这样，综合以上各个能量信息，可将分割牙根层切片的单相混合水平集模型的最 小化能量泛函表示为：
[0086]对能量方程进行最小化，可得到牙根层切片的水平集函数Φ的演化方程为：
[0088] 其中 δ(Φ)为Dirac 函数。
[0089] 步骤3、牙冠层切片分割模型:在牙冠层切片图像中，牙齿逐渐紧连在一起，缝隙较 小，为有效提取相邻两牙之间的边界，对起始切片利用步骤2所述的单相混合水平集模型进 行牙齿分割，将分割后的水平集函数按照每隔一个牙齿的规则将其组合成两个耦合的双相 混合水平集函数，如图4所示。
[0090] 由多相水平集分割原理可知，两个水平集函数可以将图像划分为四个区 域，如图5所示。对于{ΦχΧ^ΦζΧ)}这个相，代表相邻牙齿的重叠区域。为克服 演化过程中产生重叠，在能量泛函中引入区域竞争准则。设A^A 2分别为水平集函数 对应的水平集曲线围成的内部区域面积。根据Heaviside函数H(X)的性质，可将A1,A2分别表 示为：
[0091] Ai = area ( Φι>0) = ?ΩΗ(Φι)?χ
[0092] A2 = area( C>2 2 0) = JqH( C>2)dx
[0093]要使这两个区域不产生重叠，则应满足：
[0095]因此可将区域惩罚项定义为：
[0096] Erepulse = jQH( Φ?)Η( ?2)dx
[0097] 对上式进行最小化相当于避免两个牙齿产生重叠。
[0098] 将区域惩罚项加入到水平集能量泛函当中，可得到牙冠层的双相混合水平集模型 的能量方程为：
[0099] Eas( Φ 1 , Φ 2 ) = Φ I ) +Egfg( Φ 2 ) +PErepulse
[0100] 其中β用来控制两水平集函数区域重叠的程度，对上述能量泛函进行最小化，可得 到双相水平集函数φ:，φ2的演化方程分别为：
[0104] 其中Φρ?，ΦΡ2分别代表Φι，φ>2的先验值;ξ?(Χ)、ξ2(Χ)分别代表Φι ，φ>2演化时的 梯度方向检测函数，BP:
[0107] 步骤4、牙齿三维模型重建:利用窄带法和半隐式差分方案求解上述水平集函数的 演化方程，提取更新完成后的水平集函数在Φ=〇的像素点，将所有切片图像分割后的这些 牙齿轮廓像素点转化为三维坐标，利用Delaunay三角剖分方法进行重建，即可得到每颗牙 齿的三角网格模型。
[0108] 本发明的有效性可以通过下面的实验进一步说明：
[0109] 以一位下颂牙齿完整的患者口腔CT图像序列作为例。图6为利用本发明和Li模型、 CV模型对该CT序列的起始切片图像进行牙齿分割结果对比图。图6(a)为在起始切片上画的 各个牙齿初始轮廓多边形。可以看出，基于边缘的Li模型由于过于依赖初始曲线附近的图 像梯度信息，但缺乏梯度方向约束，所以图6(b)中存在很多牙齿边缘定位不准确的现象。基 于区域的C-V模型依赖全局灰度信息，因此在图6(c)中距离较近的牙齿之间以及牙齿与颂 骨之间的弱边界处产生融合。图6(d)为本发明的分割结果，由于利用基于梯度向量场的 Flux模型，能够有效探测弱边界，并通过加入梯度方向和先验形状约束，从而实现准确定位 牙齿边界轮廓。
[0110] 图7(a)和图7(b)分别本发明和Gao等人的方法对该患者第二磨牙的牙根部位切片 图像的分割结果。可以看出，本发明具有较准确的分割结果，但Gao等人的方法依靠曲线内 外灰度概率来探测轮廓，由于牙髓腔的灰度较暗，会在牙齿弱边界处使探测的曲线偏向牙 髓腔。
[0111] 图8显示了在右侧第一磨牙部位的一些切片的分割结果。可以看出，基本上所有切 片图像的分割结果都能很好地贴合牙齿目标轮廓。但是，在第160层的切片中，水平集没有 捕捉到牙齿内部沟槽轮廓，这主要是受到水平集演化方程中窄带大小的影响。由于牙齿表 面沟槽的存在，会造成一些切片的轮廓出现极度的凹陷。如果此时窄带宽度取得较小，可能 就会限制水平集演化到目标轮廓上。但这并不影响最后的分割结果，由于目标轮廓线包含 在水平集曲线的内部，可通过灰度阈值提取出最终的轮廓线。图8右边显示了对分割后的每 层切片牙齿轮廓进行三维重建后的结果，可以看出本发明重建的三维模型能正确反映牙齿 的解剖形态，从而为下一步的计算机辅助诊断提供有效数据，使医生能够针对病人的实际 情况制定最佳的手术方案。
[0112]以上所述仅为本发明的较佳实例而已，并不用以限制本发明，凡在本发明的精神 和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。
【主权项】
1. 一种基于混合水平集的三维牙齿建模方法，其特征在于，包括如下步骤： (1) 从口腔CT图像序列中选取一张作为起始切片，并在该图像上勾画出每颗牙齿轮廓 以初始化水平集函数； (2) 起始切片以下的牙根层切片利用先验形状约束能量、基于Flux模型的边缘能量、基 于先验灰度的局部区域能量相结合构造的单相混合水平集模型分割牙齿轮廓； (3) 起始切片以上的牙冠层切片利用结合区域竞争约束的双相混合水平集模型分割牙 齿轮廓； (4) 将所有切片分割后的牙齿轮廓像素点转化为三维坐标，利用Delaunay三角剖分方 法进行重建以得到每颗牙齿的三角网格模型。2. 根据权利要求1所述的一种基于混合水平集的三维牙齿建模方法，其特征在于，步骤 (1) 按如下步骤进行:在牙颈部位的切片图像中选取一张所有牙都出现且牙槽骨较少出现 的切片作为起始切片，并在该切片图像上勾画出每颗牙齿的大致轮廓C 1Q = I, 2，···η，η为 牙齿个数)作为水平集的初始轮廓，然后利用C1初始化η个水平集函数O1Q = Id, ···!〇 ,CD1 的初始化通过计算图像上每个点到匕的带符号的距离来完成，即：其中，d [(X )，Ci ]表示像素点X与曲线Ci之间的欧式距离。3. 根据权利要求1所述的一种基于混合水平集的三维牙齿建模方法，其特征在于，步骤 (2) 所述的分割牙根层切片的单相混合水平集模型的能量泛函定义为先验形状约束能量、 基于Flux模型的边缘能量、基于先验灰度的局部区域能量等的加权和：I其中，μ，γ，α， ν，λ为各个能量项的权系数； (3a)符号距离保持能量Ειη?(Φ)，用来保证水平集演化过程中的稳定性，定义为：(3b)曲线弧长平滑能量Ele3ngth( Φ )，用来平滑水平集轮廓，定义为：(3c)先验形状约束能量ΕΡ"α(Φ)，用来控制水平集的形状，将每次分割后的牙齿轮廓 映射到相邻切片图像，作为当前水平集函数演化的先验形状加以约束，其能量泛函定义为：其中Φ为当前切片的水平集函数，Φρ为上一张切片分割完成后先验形状对应的水平集 函数，H(X)为Heaviside函数； (3d)基于Flux模型的边缘能量心以〇)，用来探测牙齿的外边界轮廓，将图像梯度方向 和水平函数梯度方向之间的角度信息嵌入到传统Flux模型当中，其能量泛函定义为： Eedge( Φ ) =-?ΩξΔΙσ(?-Η( Φ ) dx 其中Λ为Laplacian算子，I。代表高斯平滑后的图像，为梯度方向检测函数，▽:为梯度算子，?代表点积； (3e)基于先验灰度的局部区域能量Ereglcin( Φ )，用来克服图像灰度不均匀问题，将先验 灰度信息嵌入到区域模型当中，其能量泛函定义为：其中fre3f_ln(X)和frrf+cmtU)分别定义为先验形状上参考点的r邻域在先验形状曲线内、 外的灰度均值，参考点为先验形状上距离当前图像目标像素点X最近的点； (3f)综合以上各能量项，将牙根层切片混合水平集模型的最小化能量泛函表示为：对上述能量泛函进行最小化，得到水平集曲线的演化方程：其中，δ( Φ )为Dirac函数。4.根据权利要求3所述的一种基于混合水平集的三维牙齿建模方法，其特征在于，步骤 (3)所述中的分割牙冠层切片的双相混合水平集模型按如下步骤建立:对起始切片利用所 述的单相混合水平集模型进行分割，将分割后的所有水平集函数按照每隔一个牙齿的规则 将其组合成两个耦合的双相混合水平集函数，并加入区域竞争约束能量以克服两个水平集 产生重叠，其能量泛函定义为：为区域竞争约束能量，β用来控制两水平集函数区域重叠的程度； 对上述能量泛函进行最小化，得到双相水平集函数φ:，φ2的演化方程分别为：其中Φρ1，ΦΡ2分别代表Φ?，Φ2的先验值;ξ?(Χ)、ξ2(Χ)分别代表Φ?，Φ2演化时的梯度方 向检测函数：05.根据权利要求1所述的一种基于混合水平集的三维牙齿建模方法，其特征在于，步骤 (4)所述的三维牙齿重建按如下步骤进行:利用窄带法和半隐式差分方案求解上述水平集 函数的演化方程，提取更新完成后的水平集函数在Φ =O的像素点，得到当前切片图像的牙 齿轮廓，将所有切片图像分割后的牙齿轮廓像素点转化为三维坐标，利用基于区域增长的 Delaunay三角剖分方法进行重建，以生成每颗牙齿的三角网格模型。
【专利摘要】本发明提供一种基于混合水平集的三维牙齿建模方法，包括以下几个步骤：1)起始切片的选取和水平集初始化；2)利用先验形状约束能量、基于Flux模型的边缘能量、基于先验灰度的局部区域能量相结合构造的单相混合水平集模型分割牙根层切片图像；3)结合区域竞争约束的双相混合水平集模型分割牙冠层切片图像；4)基于分割轮廓的牙齿三维模型重建。本发明的混合水平集模型能有效克服传统分割模型边缘定位不准确以及无法处理图像灰度不均匀等问题，快速准确地构建出每颗牙齿独立的三维模型，从而为制定口腔修复规划、生物力学分析等奠定坚实的基础。
【IPC分类】G06T17/30
【公开号】CN105488849
【申请号】CN201510821682
【发明人】吴婷, 张礼兵
【申请人】嘉兴学院
【公开日】2016年4月13日
【申请日】2015年11月24日
最新回复(0)