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

目前，关于牙齿3维重建方法主要分为两种类型：基于模型的方法和基于CT图像的方法。第1种类型首先利用光学扫描仪获取口腔结构的三角网格模型，然后再将牙齿一个个分割开来^[1-5]；这种方法获取的牙齿模型具有较高的精度，但由于缺少牙根信息，往往需要借助X片或统计模型进行牙根估计^[6]。第2种类型利用图像分割技术从每一层牙齿CT切片中分割出牙齿轮廓，然后利用这些层间轮廓重建出牙齿3维模型；由于该类方法能够获得整个牙齿形状结构，为患者口腔病变提供完整的诊断依据，而且还可以和光学扫描模型进行数据融合获得更精确的结果^[7-8]。因此，基于CT图像的牙齿建模方法越来越受到研究学者的广泛关注。

现今，基于CT图像的牙齿建模方法主要有：1)区域分割法。王黎等人^[9]利用二值化和边界提取筛选出每层切片牙齿轮廓关键点，然后利用3D-Delaunay四面体化算法得到整颗牙齿的实体模型；但二值化操作容易产生过分割或欠分割的问题。张飞等人^[10]利用参考切片的牙包围盒作为当前切片的操作区间，然后在此区间内通过区域生长算法提取牙齿轮廓；但包围盒有时不能完全覆盖当前切片的牙齿区域，还有可能将其他干涉结构也包围进来，因此这种方法在实际应用中并不可行。2)参数曲线法。韩国庆熙大学的研究学者^[11-12]利用基于遗传算法的B样条曲线拟合来提取每层切片的牙齿轮廓，但参数曲线无法处理牙齿拓扑结构变化的问题。3)水平集分割法。其核心思想是将演化的分割曲线嵌入到一个更高维水平集函数的零水平集，从而将曲线演化问题转化为求解一系列隐式表达的偏微分方程。该类方法最具有代表性的是Gao等人^[13]采用的基于统计水平集模型的方法，其主要依靠图像边缘梯度和先验形状的概率分布来指导水平集的演化，但因为牙齿密度不均匀且牙槽骨和牙齿的灰度接近，这种通过水平集曲线内外像素灰度直方图的概率比值来控制水平集轮廓收缩和扩张的方法容易产生不稳定的结果^[14]。

由于基于水平集的图像分割方法具有强大的拓扑变化处理能力，因此在Gao等人^[13]模型的基础上，根据口腔CT图像序列中牙齿形态变化和排列特点，提出一种新的混合水平集牙齿分割模型。该模型融合梯度向量场、先验图像的局部灰度和形状等信息，以提高牙齿分割的效率和精度。最后，利用分割后的层间轮廓重建出每颗牙齿独立的三角网格模型。

口腔CT图像除了对比度低、边界轮廓模糊，还具有如下特点(图 1)：1)牙齿容易发生拓扑结构变化(例如：在牙根处分裂为多根分支)；2)牙齿内部密度不均匀，且存在牙髓腔等孔洞结构；3)在牙根部位切片图像中，牙根被牙槽骨包围，二者密度相近，边缘不清晰；4)在牙冠部位切片图像，相邻牙齿粘连在一起，缝隙较小甚至消失。

从上述分析可以看出，单一水平集方法很难将每颗牙齿独立分割开来，因此提出一种混合水平集模型方法，其步骤如下：

1)首先从口腔CT图像序列中选取一张切片标记为起始切片，从起始切片开始利用水平集模型分割牙齿轮廓，每次分割后的轮廓映射到相邻切片作为初始分割和先验形状；

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

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

4)将每层分割后的水平集轮廓转换为3维坐标，利用Delaunay方法重建牙齿的三角网格模型。

整个牙齿重建路线如图 2所示。

水平集方法对图像质量和初始轮廓较为敏感，所以选取一幅好的图像和初始轮廓能够获得准确的分割结果。由于口腔CT图像在牙冠部位相邻牙齿之间距离太近，在牙根部位牙齿又被牙槽骨包围，而在牙颈部位的切片中牙齿之间相互独立，且没有牙槽骨干扰。因此，在牙颈部位附近的切片中选取一张所有牙齿都出现且牙槽骨较少的切片标记为起始切片。

在起始切片图像上，勾画出每颗牙齿的大致轮廓C_i(i=1,2,…,n,n为牙齿个数)作为水平集的初始轮廓，如图 3所示，然后利用C_i(i=1,2,…,n)初始化n个水平集函数Φ_i(i=1,2,…,n)。Φ_i的初始化通过计算图像上每个点到C_i的带符号距离来完成，即

${{\Phi }_{i}}\left( x \right)=\left\{ \begin{matrix} d\left[ \left( x \right),{{C}_{i}} \right] & x\in \text{inside}\left( {{C}_{i}} \right) \\ 0 & x\in {{C}_{i}} \\ -d\left[ \left( x \right),{{C}_{i}} \right] & x\in \text{outside}\left( {{C}_{i}} \right) \\ \end{matrix} \right.$

(1)

式中，d[(x),C_i]表示点x与曲线C_i之间的欧氏距离，inside(C_i)和outside(C_i)分别代表曲线C_i的内部和外部区域。

在牙根层图像中，牙齿易受到牙槽骨和牙髓腔等结构的较多干扰。因此，提出一种混合水平集牙齿分割模型，该模型的能量泛函主要由先验形状约束能量、基于Flux模型的边缘能量、基于先验灰度的局部区域能量等部分组成。

由于相邻切片之间，对应牙齿的形状和位置都比较接近，因此可以将上一张切片已完成分割的牙齿轮廓映射到当前切片图像，作为水平集演化的先验形状加以约束。令Φ为当前切片的水平集函数，Φ_p为上一张切片分割完成后先验形状对应的水平集函数，则先验形状约束项定义为

${{E}_{\text{prior}}}\left( \Phi \right)=\int_{\Omega }{{{\left( H\left( \Phi \right)-H\left( {{\Phi }_{p}} \right) \right)}^{2}}\text{d}x}$

(2)

式中，H(x)为Heaviside函数。

由于图像在边界处的梯度与边缘法向量同向，因此Flux模型^[15]通过计算图像梯度与曲线法向量内积的边缘积分来定位图像边缘，该方法能有效解决弱边界目标的分割问题，且演化速度快。Flux模型通过最大化能量函数得到水平集演化方程，即

$E=\underset{C}{\mathop{\oint}}\,\nabla I\cdot N\text{d}s$

(3)

式中，I为定义在平面域Ω上的待分割图像，N为水平集函数的零水平集曲线的法向量，C为Ω上的一条轮廓曲线，▽为梯度算子，·代表点积。引入变分水平集方法，可将Flux模型的最小化能量泛函表示为

${{E}_{\text{Flux}}}\left( \Phi \right)=\int_{\Omega }{H\left( \Phi \right)\Delta I\text{d}x}$

(4)

式中，Δ为Laplacian 算子。

然而实际上牙齿内部和外部的干扰结构较多，若直接利用上述模型进行分割，水平集会同时捕捉到牙齿内部的牙髓腔以及外部的邻牙或牙槽骨等边界上，因此将图像梯度方向和水平函数梯度方向之间的角度信息嵌入到上述的Flux模型当中。

在牙齿CT图像中，牙齿为高亮度区域，背景为低亮度区域，因此可将牙齿与周围干扰结构的图像梯度方向表示为图 4(a)所示。由于定义的水平集函数Φ在水平集内部取正、外部取负，则Φ的梯度方向可表示为图 4(b)所示。可以看出，如果要使演化曲线只捕捉到牙齿的外边界轮廓，则图像的梯度方向应该与水平集函数的梯度方向一致。因此，将式(4)修改为

${{E}_{\text{edge}}}\left( \Phi \right)=\int_{\Omega }{\xi H\left( \Phi \right)\Delta I\text{d}x}$

(5)

式中

$\xi \left( x \right)=\left\{ \begin{matrix} 1 & \nabla \Phi \cdot \nabla {{I}_{\sigma }}>0 \\ 0 & \nabla \Phi \cdot \nabla {{I}_{\sigma }}\le 0 \\ \end{matrix} \right.$

(6)

为梯度方向检测函数，I_σ代表高斯平滑后的图像。

传统经典的基于区域的水平集模型如CV模型^[16]利用曲线内外区域的灰度均值来进行水平集演化，但对于灰度不均匀的牙齿而言，这种全局方法不能很好地实现这种异质图像的分割。虽然Li等人提出了局部二值拟合LBF模型^[17]有效解决了这个问题，但该模型计算量大；而且以上两种方法都无法有效地利用先验区域的灰度信息，对于牙齿这类周围干扰结构较多的图像并不适用。

由于相邻切片之间图像的灰度分布非常相同，且牙齿位置也接近，因此将先验灰度信息嵌入到能量模型中，以提高比较精度。对于给定平面域Ω上的当前图像I，定义其局部区域能量为

$\begin{align} & {{E}_{\text{region}}}\left( \Phi \right)=\iint\limits_{\Omega }{{{\left( I-{{f}_{\text{ref }\!\!\_\!\!\text{ in}}} \right)}^{2}}H\left( \Phi \right)}\text{d}x+ \\ & \iint\limits_{\Omega }{{{\left( I-{{f}_{\text{ref }\!\!\_\!\!\text{ out}}} \right)}^{2}}\left( 1-H\left( \Phi \right) \right)}\text{d}x \\ \end{align}$

(7)

式中，f_{ref_out}(x)和f_{ref_out} (x)分别定义为先验形状上参考点的r邻域(如图 5所示)在先验形状曲线内、外的灰度均值。参考点按照如下方法确定：当上一张切片的先验轮廓分割完成后，计算先验形状上每个点的r邻域在先验形状曲线内、外的灰度均值，然后将先验形状曲线上距离当前切片图像上目标像素点x最近的点作为参考点。

为保证水平集在整个演化过程中的形状和稳定性，加入曲线弧长平滑项^[16]

${{E}_{\text{length}}}=\int_{\Omega }{\left| \nabla H\left( \Phi \right) \right|\text{d}x}$

(8)

和符号距离保持项^[18]

${{E}_{\text{int}}}\left( \Phi \right)=\int_{\Omega }{\frac{1}{2}{{\left( \left| \nabla H \right|-1 \right)}^{\text{2}}}\text{d}x}$

(9)

综合以上各个能量信息，可将牙根层切片水平集模型的最小化能量泛函表示为

$\begin{align} & {{E}_{\text{root}}}\left( \Phi \right)=\underbrace{\mu \int_{\Omega }{\frac{1}{2}{{\left( \left| \nabla \Phi \right|-1 \right)}^{2}}\text{d}x+}}_{符号距离保持项} \\ & \underbrace{\gamma \int_{\Omega }{\left| \nabla H\left( \Phi \right) \right|\text{d}x+}}_{弧长平滑项} \\ & \underbrace{\alpha \int_{\Omega }{{{\left( H\left( \Phi \right)-H\left( {{\Phi }_{p}} \right) \right)}^{2}}\text{d}x+}}_{先验形状约束项} \\ & \underbrace{\upsilon \int_{\Omega }{\xi \left( x \right)H\left( \Phi \right)\Delta I\text{d}x+}}_{边缘能量项} \\ & \underbrace{\lambda \left[ \begin{align} & \iint\limits_{\Omega }{{{\left( I-{{f}_{\text{ref}\_\text{in}}} \right)}^{2}}H\left( \Phi \right)\text{d}x+} \\ & \iint\limits_{\Omega }{{{\left( I-{{f}_{\text{ref}\_\text{out}}} \right)}^{2}}\left( 1-H\left( \Phi \right) \right)\text{d}x} \\ \end{align} \right]}_{局部区域能量项} \\ \end{align}$

(10)

式中，μ,γ,α,ν,λ为各个能量项的权值参数。

另外，尽管相邻切片之间牙齿轮廓形状相似，但是由于牙齿表面沟槽的存在，可能会造成一些切片图像的牙齿外轮廓存在极度的凹陷，从而较大地偏离先验形状。因此，通过计算曲线的曲率来解决。

水平集曲线上每个点的曲率为

$\kappa =-\text{div}\left( \frac{\nabla \Phi }{\left| \nabla \Phi \right|} \right)$

(11)

式中，div为散度算子。为使水平集曲线更好地捕捉到凹陷部位，在曲线的凹陷部位即负曲率时，降低形状约束能量项的参数α，同时增大边缘和区域能量项的参数ν, λ,即

$\left\{ \begin{matrix} \alpha *=\omega \alpha ,v*=\frac{1}{\omega }v,\lambda *=\frac{1}{3\omega }\lambda & \kappa <0 \\ \alpha *=\alpha ,v*=v & \kappa \ge 0 \\ \end{matrix} \right.$

(12)

式中，0<ω<1。这样，将式(12)代入能量方程式(10)，并对能量方程进行最小化，可得到牙根层切片的水平集函数Φ的演化方程为

$\begin{align} & \frac{\partial \Phi }{\partial t}=\mu \left[ \nabla \Phi -\text{div}\left( \frac{\nabla \Phi }{\left| \nabla \Phi \right|} \right) \right]+ \\ & \delta \left( \Phi \right)\left[ \begin{matrix} +\gamma \text{div}\left( \frac{\nabla \Phi }{\left| \nabla \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}_{\text{ref}\_\text{in}}} \right)}^{2}}-{{\left( I-{{f}_{\text{ref}\_\text{out}}} \right)}^{2}} \right] \\ \end{matrix} \right] \\ \end{align}$

(13)

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

在牙冠层切片图像中，牙齿逐渐紧连在一起，缝隙较小，为有效提取相邻两牙之间的边界，将起始切片分割后的n个水平集的内部区域按照每隔一个牙齿的规则将其分为两组，然后分别初始化为两个耦合的双相水平集函数，如图 6所示。

由多相水平集分割原理可知，两个水平集函数Φ₁、Φ₂可以将图像划分为4个区域^[19]，如图 7所示。其中，{Φ₁>0，Φ₂>0}这个区域，代表相邻牙齿的重叠区域。为克服Φ₁、Φ₂在演化过程中产生重叠，在能量泛函中引入区域竞争约束^[20]以对牙冠层切片进行有效分割。设A₁、A₂分别为水平集函数Φ₁、Φ₂对应的水平集曲线围成的内部区域面积。根据H(x) 函数的性质，可将A₁，A₂分别表示为

$\begin{align} & {{A}_{1}}=area\left( {{\Phi }_{1}}\ge 0 \right)=\int_{\Omega }{H\left( {{\Phi }_{1}} \right)}\text{d}x \\ & {{A}_{2}}=area\left( {{\Phi }_{2}}\ge 0 \right)=\int_{\Omega }{H\left( {{\Phi }_{2}} \right)}\text{d}x \\ \end{align}$

(14)

要使{Φ₁>0，Φ₂>0}这个重叠区域面积最小，可将区域惩罚项定义为

${{E}_{\text{repulse}}}=\int_{\Omega }{H\left( {{\Phi }_{1}} \right)H\left( {{\Phi }_{2}} \right)\text{d}x}$

(15)

对上式进行最小化相当于避免两个相邻牙齿产生重叠。

将区域惩罚项加入到水平集能量泛函当中，可得到牙冠层的双相水平集混合模型的能量方程为

$\begin{align} & {{E}_{\text{crown}}}\left( {{\Phi }_{1}},{{\Phi }_{2}} \right)={{E}_{\text{root}}}\left( {{\Phi }_{1}} \right)+ \\ & {{E}_{\text{root}}}\left( {{\Phi }_{2}} \right)+\beta {{E}_{\text{repulse}}} \\ \end{align}$

(16)

式中，β用来控制两水平集函数区域重叠的程度。对能量泛函(16)进行最小化，可得到双相水平集函数Φ₁，Φ₂的演化方程分别为

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

(17)

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

(18)

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

${{\xi }_{1}}\left( x \right)=\left\{ \begin{matrix} 1 & \nabla {{\Phi }_{1}}\cdot \nabla {{I}_{1\sigma }}>0 \\ 0 & \nabla {{\Phi }_{1}}\cdot \nabla {{I}_{1\sigma }}\le 0 \\ \end{matrix} \right.$

(19)

${{\xi }_{2}}\left( x \right)=\left\{ \begin{matrix} 1 & \nabla {{\Phi }_{2}}\cdot \nabla {{I}_{2\sigma }}>0 \\ 0 & \nabla {{\Phi }_{2}}\cdot \nabla {{I}_{2\sigma }}\le 0 \\ \end{matrix} \right.$

(20)

为验证本文方法的有效性，从口腔医院的CT数据库中选取若干数据进行测试。实验在CPU为2.3 GHz、8 GB内存的PC机上运行，操作系统为64位Windows 8，软件开发平台为Matlab R2014b。所有图像的像素间距为0.3 mm×0.3 mm，层间距为0.3 mm。

由于标准CT数据一般以DICOM格式存储，灰度级分辨率为16 bit。为加快计算速度，首先将CT数据统一转换为8 bit的灰度图像，然后采用半隐式差分方案和窄带法求解水平集函数的演化方程，窄带宽度设为8个像素。在整个分割实验当中，取参数μ=0.01,γ=0.005×255×255,ν=50,α=100,λ=0.1, ω=0.25, β=1 000。

实验1 利用本文模型和CV模型^[16]、Li模型^[18]对某CT图像序列的起始切片图像进行牙齿分割，其分割结果如图 8所示。图 8(a)为在起始切片上画的各个牙齿的初始轮廓多边形。可以看出，基于区域的CV模型依赖全局灰度信息，因此在距离较近的牙齿之间以及牙齿与颌骨之间的弱边界处产生了融合。基于边缘的Li模型由于过于依赖初始曲线附近的图像梯度信息，但缺乏梯度方向约束，所以存在边缘定位不准确的现象，有些定位到牙齿内边缘处。本文所提出的混合模型方法，由于利用基于梯度向量场的Flux模型，能够有效探测弱边界，并通过加入梯度方向和先验形状约束，从而实现准确定位牙齿边界轮廓。

实验2 利用本文模型和传统的单CV模型^[16]、双CV模型^[19]对牙冠层切片进行分割实验，其结果如图 9所示。可以发现，由于牙冠层牙齿之间粘连在一起，采用单CV模型会将几颗牙齿分割成一个整体，即使采用耦合的双CV模型进行分割，由于没有区域竞争惩罚项，两个水平集相互交错在一起，无法将两两相邻的牙齿轮廓分割开来。本文的双水平集混合模型由于加入了区域竞争惩罚项，能够很好地避免两个水平集轮廓之间产生融合。

实验3 由于牙根部位的切片图像容易受到牙髓腔和牙槽骨的较多干扰，因此利用本文模型和Gao等人^[13]的模型对10幅不同牙根部位的切片图像进行分割实验，并利用相似系数S_agr和不相似系数S_dis^[12]进行定量分析比较，这两个系数都是用来衡量算法分割和手动分割区域一致和不一致的程度，其计算公式分别为

$\begin{align} & {{S}_{\arg }}=2\frac{\left| {{A}_{\text{man}}}\cap {{A}_{\text{alg}}} \right|}{{{\left| A \right.}_{\text{man}}}+{{A}_{\text{alg}}}} \\ & {{S}_{\text{dis}}}=2\frac{\left| {{A}_{\text{man}}}\cup {{A}_{\text{alg}}}-{{A}_{\text{man}}}\cap {{A}_{\text{alg}}} \right|}{\left| {{A}_{\text{man}}}+{{A}_{\text{alg}}} \right|} \\ \end{align}$

式中，A_man和A_alg分别表手动分割和算法分割的区域，|·|表示区域面积。

图 10为一张典型牙根切片的分割结果对比，可以看出Gao模型产生的分割轮廓很容易在牙齿内部较暗区域产生错误的内凹，这是因为该方法是主要依靠先验形状曲线窄带内外直方图分布的概率比值来控制水平集的演化，即如果当前图像目标像素点的灰度值在先验形状内部区域的灰度概率大于在先验形状外部区域的灰度概率，则曲线向内收缩，否则向外扩张；但由于牙齿密度不均匀，这种方法很容易将牙齿内部较暗部分误判为牙齿外部区域。本文模型通过利用先验形状局部区域的灰度值，能够较好处理牙齿灰度不均问题，所以分割结果较Gao模型有较大改善。从最后表 1的分割参数结果也可看出，本文方法得到的相似系数S_agr比Gao方法提高了将近3%，而不相似系数S_dis下降了5%，所以本文方法较Gao方法在分割精度上有一定的改善。

本文的分割方法在牙根层的能量方程中一共需要6个权值参数μ、γ、ν、α、λ、ω，牙冠层能量方程还需要增加一个区域竞争参数β，这些参数用来控制每个能量项的影响。在起始分割中，边缘能量项起主要作用，而在后续的分割中先验形状和先验灰度区域能量项则逐渐发挥作用。通过大量的测试发现，实验所使用的参数值适用于大部分相似分辨率的CT图像。

由于本文分割策略基于层间映射机制，因此相邻切片之间的形状差异对分割结果起着重要的影响。尤其在牙根底部和牙冠上部，由于拓扑结构变化，可能造成相邻切片之间的形状差异较大，如果采用固定的权值参数可能会产生不理想的分割结果，因此自适应权值参数ω的调整十分重要。图 11为利用自适应权值参数和常值参数方法对3张连续切片图像进行分割实验的结果。可以看出，在编号#1和编号#2的切片图像中，该牙齿形状在中间部位差别比较大，即该牙齿在切片#2图像中突然分裂为两部分。如果采用常值的权值参数，由于先验形状以及窄带的约束，会限制水平集演化到牙齿轮廓，即使到了切片#3也不能产生正确的结果，如图 11(a)所示。但本文可以通过水平集曲率的大小自动调整权值参数(本实验设ω=0.3)，来降低形状约束能量项的参数α，并增大边缘和区域能量项的参数ν,λ(详见式(12))。从图 11(b)的结果可以看出，虽然在切片#2中水平集没有完全贴合牙齿(主要是由于计算过程中窄带宽度的限制)，但演化到切片#3时，水平集能够正确捕捉牙齿轮廓，从而较好地改善相邻切片之间牙齿形状突变问题。参数ω的值越小，则水平集函数受曲率的影响越大，但在分割过程中ω不宜过小，以免改变水平集的牙齿形状。通过测试发现，ω不易小于0.1为合适。

对一名下颌牙齿完整的患者口腔CT图像序列进行牙齿3维重建实例验证，图 12显示了若干切片的分割结果。可以看出，大部分切片图像的分割结果较为理想，但是，在第100层的牙冠上部切片图像中，水平集只捕捉到该牙齿的外轮廓，内部的极限凹槽结构没有捕捉到，其原因即为上面所描述的牙齿不规则凹陷或拓扑变化结构会导致相邻切片的形状差异较大，而先验形状和窄带大小有可能限制水平集的演化，但通过自适应参数的调整，可以慢慢演化到正确的轮廓。

将所有切片图像分割后的牙齿轮廓像素点转化为3维坐标，利用Delaunay三角剖分方法进行重建，即可生成每颗牙齿的三角网格模型，如图 13所示。从最后的结果可以看出，本文方法能够将每颗牙齿从牙槽骨里分割出来，并正确反映牙齿的解剖形态，这样每颗牙齿都可以单独的操作，从而可以为下一步的计算机辅助诊断提供有效的口腔结构数据，使医生能够为病人制定最佳的手术方案。

针对牙齿CT图像序列的形状和排列特点，提出一种基于混合水平集模型的CT图像牙齿分割方法。该模型融合梯度向量场、先验图像的局部灰度和形状等信息，能有效克服牙髓腔和周围颌骨的干扰以及图像灰度不均的问题；并针对牙冠层切片中相邻牙相互粘连现象，采用双相水平集模型和区域竞争约束以提取相邻牙齿之间的共同边界。与现有的方法相比，本文方法具有较好的分割效果和较高的准确率，而且具有较少的人工干预。但对于牙齿存在的极度拓扑变化结构，本文方法可能在这些位置的切片图像中无法获得精确的分割结果，这些将是我们后续工作的重点。