专利摘要

本发明提供一种基于混合水平集的三维牙齿建模方法，包括以下几个步骤：1)起始切片的选取和水平集初始化；2)利用先验形状约束能量、基于Flux模型的边缘能量、基于先验灰度的局部区域能量相结合构造的单相混合水平集模型分割牙根层切片图像；3)结合区域竞争约束的双相混合水平集模型分割牙冠层切片图像；4)基于分割轮廓的牙齿三维模型重建。本发明的混合水平集模型能有效克服传统分割模型边缘定位不准确以及无法处理图像灰度不均匀等问题，快速准确地构建出每颗牙齿独立的三维模型，从而为制定口腔修复规划、生物力学分析等奠定坚实的基础。

权利要求

1.一种基于混合水平集的三维牙齿建模方法，其特征在于，包括如下步骤：

(1)从口腔CT图像序列中选取一张作为起始切片，并在该图像上勾画出每颗牙齿轮廓以初始化水平集函数；

(2)起始切片以下的牙根层切片利用先验形状约束能量、基于Flux模型的边缘能量、基于先验灰度的局部区域能量相结合构造的单相混合水平集模型分割牙齿轮廓；

(3)起始切片以上的牙冠层切片利用结合区域竞争约束的双相混合水平集模型分割牙齿轮廓；

(4)将所有切片分割后的牙齿轮廓像素点转化为三维坐标，利用Delaunay三角剖分方法进行重建以得到每颗牙齿的三角网格模型。

2.根据权利要求1所述的一种基于混合水平集的三维牙齿建模方法，其特征在于，步骤(1)按如下步骤进行：在牙颈部位的切片图像中选取一张所有牙都出现且牙槽骨较少出现的切片作为起始切片，并在该切片图像上勾画出每颗牙齿的大致轮廓C_i(i＝1，2，…n，n为牙齿个数)作为水平集的初始轮廓，然后利用C_i初始化n个水平集函数Ф_i(i＝1，2，…n)，Ф_i的初始化通过计算图像上每个点到C_i的带符号的距离来完成，即：

${\Phi }_i\left(x\right)=\left\{\begin{array}{cc}d\left[\left(x\right),C_i\right],& x\in inside\left(C_i\right)\\ 0,& x\in C_i\\ -d\left[\left(x\right),C_i\right],& x\in outside\left(C_i\right)\end{array}\right.$

其中，d[(x)，C_i]表示像素点x与曲线C_i之间的欧式距离。

3.根据权利要求1所述的一种基于混合水平集的三维牙齿建模方法，其特征在于，步骤(2)所述的分割牙根层切片的单相混合水平集模型的能量泛函定义为先验形状约束能量、基于Flux模型的边缘能量、基于先验灰度的局部区域能量等的加权和：

E_牙根(Ф)＝μE_int(Ф)+γE_length(Ф)+αE_prior(Ф)+νE_edge(Ф)+λE_region(Ф)其中，μ，γ，α，ν，λ为各个能量项的权系数；

(3a)符号距离保持能量E_int(Ф)，用来保证水平集演化过程中的稳定性，定义为：

$E_{\mathrm{int}}\left(\Phi \right)={\int }_{\Omega }\frac{1}{2}{\left(\right|▿\Phi |-1)}^{2}dx$

(3b)曲线弧长平滑能量E_length(Ф)，用来平滑水平集轮廓，定义为：

$E_{length}={\int }_{\Omega }\left|▿H\left(\Phi \right)\right|dx$

(3c)先验形状约束能量E_prior(Ф)，用来控制水平集的形状，将每次分割后的牙齿轮廓映射到相邻切片图像，作为当前水平集函数演化的先验形状加以约束，其能量泛函定义为：

$E_{prior}\left(\Phi \right)=\underset{\Omega }{\int }{\left(H\right(\Phi )-H({\Phi }_p\left)\right)}^{2}dx$

其中Ф为当前切片的水平集函数，Ф_p为上一张切片分割完成后先验形状对应的水平集函数，H(x)为Heaviside函数；

(3d)基于Flux模型的边缘能量E_edge(Ф)，用来探测牙齿的外边界轮廓，将图像梯度方向和水平函数梯度方向之间的角度信息嵌入到传统Flux模型当中，其能量泛函定义为：

E_edge(Ф)＝-∫_Ωξ△I_σ(1-H(Ф)dx

其中△为Laplacian算子，I_σ代表高斯平滑后的图像，

$\xi \left(x\right)=\left\{\begin{array}{ccc}1,& if& ▿\Phi ·▿I_{\sigma }>0\\ 0,& if& ▿\Phi ·▿I_{\sigma }\le 0\end{array}\right.$

为梯度方向检测函数， 为梯度算子，·代表点积；

(3e)基于先验灰度的局部区域能量E_region(Ф)，用来克服图像灰度不均匀问题，将先验灰度信息嵌入到区域模型当中，其能量泛函定义为：

$E_{region}\left(\Phi \right)=\underset{\Omega }{\int \int }{\left(I-f_{ref\_in}\right)}^{2}H\left(\Phi \right)dx+\underset{\Omega }{\int \int }{\left(I-f_{ref\_out}\right)}^2\left(1-H\left(\Phi \right)\right)dx$

其中f_{ref_in}(x)和f_{ref_out}(x)分别定义为先验形状上参考点的r邻域在先验形状曲线内、外的灰度均值，参考点为先验形状上距离当前图像目标像素点x最近的点；

(3f)综合以上各能量项，将牙根层切片混合水平集模型的最小化能量泛函表示为：

对上述能量泛函进行最小化，得到水平集曲线的演化方程：

$\begin{array}{c}\frac{\partial \Phi }{\partial t}=\mu \left[\Delta \Phi -div\left(\frac{▿\Phi }{\left|▿\Phi \right|}\right)\right]\\ +\delta \left(\Phi \right)\left[\begin{array}{c}+\gamma div\left(\frac{▿\Phi }{\left|▿\Phi \right|}\right)\\ -2\alpha \left(H\left(\Phi \right)-H\left({\Phi }_p\right)\right)\\ -v\xi \left(x\right)\Delta I\\ -\lambda \left[{\left(I-f_{ref\_in}\right)}^2-{\left(I-f_{ref\_iout}\right)}^2\right]\end{array}\right]\end{array}$

其中，δ(Ф)为Dirac函数。

4.根据权利要求3所述的一种基于混合水平集的三维牙齿建模方法，其特征在于，步骤(3)所述中的分割牙冠层切片的双相混合水平集模型按如下步骤建立：对起始切片利用所述的单相混合水平集模型进行分割，将分割后的所有水平集函数按照每隔一个牙齿的规则将其组合成两个耦合的双相混合水平集函数，并加入区域竞争约束能量以克服两个水平集产生重叠，其能量泛函定义为：

E_牙冠(Ф₁，Ф₂)＝E_牙根(Ф₁)+E_牙根(Ф₂)+βE_repulse

其中，

E_repulse＝∫_ΩE(Ф₁)H(Ф₂)dx

为区域竞争约束能量，β用来控制两水平集函数区域重叠的程度；

对上述能量泛函进行最小化，得到双相水平集函数Ф₁，Ф₂的演化方程分别为：

$\frac{\partial {\Phi }_1}{\partial t}=\mu [{\mathrm{\Delta \Phi }}_1-div\left(\frac{▿{\Phi }_1}{|▿{\Phi }_1|}\right)]$

和$+\delta \left({\Phi }_1\right)\left[\begin{array}{c}+\gamma div\left(\frac{▿{\Phi }_1}{\left|▿{\Phi }_1\right|}\right)\\ -2{\alpha }^{*}\left(H\left({\Phi }_1\right)-H\left({\Phi }_{p1}\right)\right)\\ -v^{*}{\xi }_1\Delta I\\ -\lambda \left[{\left(I-f_{ref\_in}\right)}^2-{\left(I-f_{ref\_iout}\right)}^2\right]\\ -\beta H\left({\Phi }_2\right)\end{array}\right]$

和

$\begin{array}{c}\frac{\partial {\Phi }_2}{\partial t}=\mu \left[{\mathrm{\Delta \Phi }}_2-div\left(\frac{▿{\Phi }_2}{\left|▿{\Phi }_2\right|}\right)\right]\\ +\delta \left({\Phi }_2\right)\left[\begin{array}{c}+\gamma div\left(\frac{▿{\Phi }_2}{\left|▿{\Phi }_2\right|}\right)\\ -2{\alpha }^{*}\left(H\left({\Phi }_2\right)-H\left({\Phi }_{p2}\right)\right)\\ -v^{*}{\xi }_2\Delta I\\ -\lambda \left[{\left(I-f_{ref\_in}\right)}^2-{\left(I-f_{ref\_iout}\right)}^2\right]\\ -\beta H\left({\Phi }_1\right)\end{array}\right]\end{array}$

其中Φ_p1，Φ_p2分别代表Ф₁，Ф₂的先验值；ξ₁(x)、ξ₂(x)分别代表Ф₁，Ф₂演化时的梯度方向检测函数：

${\xi }_1\left(x\right)=\left\{\begin{array}{ccc}1,& if& ▿{\Phi }_1·▿I_{1\sigma }>0\\ 0,& if& ▿{\Phi }_1·▿I_{1\sigma }\le 0\end{array}\right.$

${\xi }_2\left(x\right)=\left\{\begin{array}{ccc}1,& if& ▿{\Phi }_2·▿I_{2\sigma }>0\\ 0,& if& ▿{\Phi }_2·▿I_{2\sigma }\le 0\end{array}\right..$

5.根据权利要求1所述的一种基于混合水平集的三维牙齿建模方法，其特征在于，步骤(4)所述的三维牙齿重建按如下步骤进行：利用窄带法和半隐式差分方案求解上述水平集函数的演化方程，提取更新完成后的水平集函数在Ф＝0的像素点，得到当前切片图像的牙齿轮廓，将所有切片图像分割后的牙齿轮廓像素点转化为三维坐标，利用基于区域增长的Delaunay三角剖分方法进行重建，以生成每颗牙齿的三角网格模型。

说明书

技术领域

本发明属于图像处理领域，涉及一种水平集方法在三维牙齿建模中的应用。

背景技术

近年来，随着三维数字化成像技术在口腔医学领域的飞速发展，计算机辅助诊疗修复越来越多的应用在口腔修复当中，并逐渐成为这一领域的发展趋势。在计算机口腔修复系统中，首先重要的就是要获取数字化的三维牙齿模型，而模型的精度和完整性直接关系到后续的排牙、种植、正畸、以及生物力学分析的结果。

目前牙齿建模方法最常用的方法就是利用图像处理技术分割牙齿CT图像序列，即首先从每一层CT切片中分割出牙齿轮廓，然后利用这些层间轮廓重建出牙齿三维模型。由于该类方法能够获得整个牙齿形状结构，为患者口腔病变提供完整的诊断依据，因此，基于CT图像的牙齿建模方法越来越受到研究学者的广泛关注。

由于口腔CT图像中牙齿和颌骨的密度和距离都较为接近，采用传统的图像分割方法很难精确地提取出每颗牙齿的组织轮廓。CT图像牙齿分割一直是一个充满挑战性的课题。王黎等(王黎，崔进，韩清凯等，基于CT图像的牙齿3维实体模型建立。中国图象图形学报，2005，10(10)：1289-1292。)利用二值化和边界提取筛选出每层切片牙齿轮廓关键点，然后利用3D-Delaunay四面体化算法得到整颗牙齿的实体模型；但二值化操作容易产生过分割或欠分割的问题。Wu等(X.Wu，H.Gao，H.Heo，etal.ImprovedB-splinecontourfittingusinggeneticalgorithmforthesegmentationofdentalcomputerizedtomographyimagesequences.TheJournalofimagingscienceandtechnology，2007，51(4)：328-336.)利用基于遗传算法的B样条曲线拟合来提取每层切片的牙齿轮廓，但B样条曲线无法处理牙齿拓扑结构变化的问题。Gao等(H.Gao，O.Chae.IndividualtoothsegmentationfromCTimagesusinglevelsetmethodwithshapeandintensityprior.PatternRecognition，2010，43：2406-2417)采用基于水平集的活动轮廓模型进行牙齿分割，并对不同牙层切片采用不同的分割模型，能够将每颗牙齿的牙冠和牙根都分割出来；该模型主要依靠图像边缘梯度和先验形状的概率分布来指导水平集的演化，但由于牙齿密度不均匀，且周围容易受到牙槽骨等结构的干扰，因此依靠先验形状周围区域的灰度概率来控制水平集轮廓的收缩和扩张容易产生边界泄露的问题。

发明内容

为克服上述现有技术的不足，提高牙齿分割的效率和精度，本发明综合考虑口腔CT图像牙齿形状变化特点，提出一种基于混合水平集的三维牙齿建模方法。

本发明提供：一种基于混合水平集的三维牙齿建模方法，包括如下步骤：

(1)从口腔CT图像序列中选取一张作为起始切片，并在该图像上勾画出每颗牙齿轮廓以初始化水平集函数；

(2)起始切片以下的牙根层切片利用先验形状约束能量、基于Flux模型的边缘能量、基于先验灰度的局部区域能量相结合构造的单相混合水平集模型分割牙齿轮廓；

(3)起始切片以上的牙冠层切片利用结合区域竞争约束的双相混合水平集模型分割牙齿轮廓；

(4)将所有切片分割后的牙齿轮廓像素点转化为三维坐标，利用Delaunay三角剖分方法进行重建以得到每颗牙齿的三角网格模型。

步骤(1)按如下步骤进行：在牙颈部位的切片图像中选取一张所有牙都出现且牙槽骨较少出现的切片作为起始切片，并在该切片图像上勾画出每颗牙齿的大致轮廓C_i(i＝1，2，...n，n为牙齿个数)作为水平集的初始轮廓，然后利用C_i初始化n个水平集函数Φ_i(i＝1，2，...n)，Φ_i的初始化通过计算图像上每个点到C_i的带符号的距离来完成，即：

Φi(x)=d[(x),Ci],x∈inside(Ci)0,x∈Ci-d[(x),Ci],x∈outside(Ci)]]>

其中，d[(x)，C_i]表示像素点x与曲线C_i之间的欧式距离。

步骤(2)所述的分割牙根层切片的单相混合水平集模型的能量泛函定义为先验形状约束能量、基于Flux模型的边缘能量、基于先验灰度的局部区域能量等的加权和：

E_牙根(Φ)＝μE_int(Φ)+γE_length(Φ)+αE_prior(Φ)+vE_edge(Φ)+λE_region(Φ)

其中，μ，γ，α，ν，λ为各个能量项的权系数；

(3a)符号距离保持能量E_int(Φ)，用来保证水平集演化过程中的稳定性，定义为：

Eint(Φ)=∫Ω12(|▿Φ|-1)2dx]]>

(3b)曲线弧长平滑能量E_length(Φ)，用来平滑水平集轮廓，定义为：

Elength=∫Ω|▿H(Φ)|dx]]>

(3c)先验形状约束能量E_prior(Φ)，用来控制水平集的形状，将每次分割后的牙齿轮廓映射到相邻切片图像，作为当前水平集函数演化的先验形状加以约束，其能量泛函定义为：

Eprior(Φ)=∫Ω(H(Φ)-H(Φp))2dx]]>

其中Φ为当前切片的水平集函数，Φ_p为上一张切片分割完成后先验形状对应的水平集函数，H(x)为Heaviside函数；

(3d)基于Flux模型的边缘能量E_edge(Φ)，用来探测牙齿的外边界轮廓，将图像梯度方向和水平函数梯度方向之间的角度信息嵌入到传统Flux模型当中，其能量泛函定义为：

E_edge(Φ)＝-∫_ΩξΔI_σ(1-H(Φ)dx

其中Δ为Laplacian算子，I_σ代表高斯平滑后的图像，

ξ(x)=1,if▿Φ·▿Iσ>00,if▿Φ·▿Iσ≤0]]>

为梯度方向检测函数， 为梯度算子，·代表点积；

(3e)基于先验灰度的局部区域能量E_region(Φ)，用来克服图像灰度不均匀问题，将先验灰度信息嵌入到区域模型当中，其能量泛函定义为：

Eregion(Φ)=∫∫Ω(I-fref_in)2H(Φ)dx+∫∫Ω(I-fref_out)2(1-H(Φ))dx]]>

其中f_{ref_in}(x)和f_{ref_out}(x)分别定义为先验形状上参考点的r邻域在先验形状曲线内、外的灰度均值，参考点为先验形状上距离当前图像目标像素点x最近的点；

(3f)综合以上各能量项，将牙根层切片混合水平集模型的最小化能量泛函表示为：

对上述能量泛函进行最小化，得到水平集曲线的演化方程：

∂Φ∂t=μ[ΔΦ-div(▿Φ|▿Φ|)]+δ(Φ)+γdiv(▿Φ|▿Φ|)-2α(H(Φ)-H(Φp))-vξ(x)ΔI-λ[(I-fref_in)2-(I-fref_iout)2]]]>

其中，δ(Φ)为Dirac函数。

步骤(3)所述中的分割牙冠层切片的双相混合水平集模型按如下步骤建立：对起始切片利用所述的单相混合水平集模型进行分割，将分割后的所有水平集函数按照每隔一个牙齿的规则将其组合成两个耦合的双相混合水平集函数，并加入区域竞争约束能量以克服两个水平集产生重叠，其能量泛函定义为：

E_牙冠(Φ₁，Φ₂)＝E_牙根(Φ₁)+E_牙根(Φ₂)+βE_repulse

其中，E_repulse＝∫_ΩH(Φ₁)H(Φ₂)dx

为区域竞争约束能量，β用来控制两水平集函数区域重叠的程度；

对上述能量泛函进行最小化，得到双相水平集函数Φ₁，Φ₂的演化方程分别为：∂Φ1∂t=μ[ΔΦ1-div(▿Φ1|▿Φ1|)]]]>

和+δ(Φ1)+γdiv(▿Φ1|▿Φ1|)-2α*(H(Φ1)-H(Φp1))-v*ξ1ΔI-λ[(I-fref_in)2-(I-fref_iout)2]-βH(Φ2)]]>

和∂Φ2∂t=μ[ΔΦ2-div(▿Φ2|▿Φ2|)]+δ(Φ2)+γdiv(▿Φ2|▿Φ2|)-2α*(H(Φ2)-H(Φp2))-v*ξ2ΔI-λ[(I-fref_in)2-(I-fref_iout)2]-βH(Φ1)]]>

其中Φ_p1，Φ_p2分别代表Φ₁，Φ₂的先验值；ξ₁(x)、ξ₂(x)分别代表Φ₁，Φ₂演化时的梯度方向检测函数：

ξ1(x)=1,if▿Φ1·▿I1σ>00,if▿Φ1·▿I1σ≤0]]>

ξ2(x)=1,if▿Φ2·▿I2σ>00,if▿Φ2·▿I2σ≤0.]]>

步骤(4)所述的三维牙齿重建按如下步骤进行：利用窄带法和半隐式差分方案求解上述水平集函数的演化方程，提取更新完成后的水平集函数在Φ＝0的像素点，得到当前切片图像的牙齿轮廓，将所有切片图像分割后的牙齿轮廓像素点转化为三维坐标，利用基于区域增长的Delaunay三角剖分方法进行重建，以生成每颗牙齿的三角网格模型。

本发明的有益效果：本发明提出的混合水平集模型结合了图像边缘、局部区域、先验知识等多方面的信息，能有效克服传统模型边缘定位不准确以及无法处理图像灰度不均匀等问题。本发明方法具有较少的人工干预，而且具有较好的分割效果和较高的准确率，重建出的三维牙齿模型能正确反映牙齿的解剖形态，从而为下一步制定口腔修复规划、生物力学分析等奠定坚实的基础，对于提高口腔修复的精度和效率具有重要意义。

附图说明

图1为本发明的牙齿建模技术流程图。

图2为起始切片的水平集函数初始化示意图。

图3(a)为牙齿与周围干扰结构的梯度方向示意图。

图3(b)为水平集函数的梯度方向示意图。

图4为牙冠层切片图像的双相水平集函数示意图。

图5为双相水平集分割区域示意图。

图6为本发明与现有方法对起始切片图像的分割效果图。

图7为本发明与现有方法对牙根切片磨牙部位的分割效果图。

图8为本发明磨牙部位若干切片分割轮廓和三维重建效果图。

具体实施方式

下面结合附图对本发明实施例作进一步说明：

如图1所示，本发明的具体实现步骤如下：

步骤1、起始切片的选取和水平集初始化：首先读取牙齿CT图像序列，然后在牙颈部位的切片图像中选择一张所有牙都出现且牙槽骨较少出现的切片作为起始切片。在起始切片图像上，勾画出每颗牙齿的大致轮廓C_i(i＝1，2，...n，n为牙齿个数)作为水平集的初始轮廓，如图2所示，然后利用C_i(i＝1，2，...n)初始化n个水平集函数Φ_i(i＝1，2，...n)。Φ_i的初始化通过计算图像上每个点到C_i的带符号的距离来完成，即：

Φi(x)=

一种基于混合水平集的三维牙齿建模方法专利购买费用说明

Q：办理专利转让的流程及所需资料

A：专利权人变更需要办理著录项目变更手续，有代理机构的，变更手续应当由代理机构办理。

1：专利变更应当使用专利局统一制作的“著录项目变更申报书”提出。

2：按规定缴纳著录项目变更手续费。

3：同时提交相关证明文件原件。

4：专利权转移的，变更后的专利权人委托新专利代理机构的，应当提交变更后的全体专利申请人签字或者盖章的委托书。

Q:专利著录项目变更费用如何缴交

A：（1）直接到国家知识产权局受理大厅收费窗口缴纳,（2）通过代办处缴纳，（3）通过邮局或者银行汇款,更多缴纳方式

Q：专利转让变更，多久能出结果

A：著录项目变更请求书递交后，一般1-2个月左右就会收到通知，国家知识产权局会下达《转让手续合格通知书》。