随着3维数字化成像技术在口腔医学领域的飞速发展，计算机辅助诊疗修复越来越多的应用在口腔修复当中，并逐渐成为这一领域的发展趋势^[1]。在计算机口腔正畸系统当中，一个关键的步骤就是将数字化3维牙颌模型上的牙齿一个个单独分割出来，而分割的精度和效率直接关系到后续的虚拟牙排列和治疗规划的结果^[2]。

受3维测量的分辨率和网格重建精度的限制，三角网格牙颌模型上龈缘和牙齿缝隙往往存在着粘连问题，再加上牙齿表面凹凸不平，使得单颗牙齿的自动分割变得极为困难，尤其当遇到患者存在牙齿畸形的状况，这一问题变得更为严重。尽管目前出现了许多商业化的牙科软件系统，例如“3shape”，实现了一定程度的牙齿分割自动化，但是它们对手工干预的程度较高，分割效率较低。由于牙齿形状各异，并且不同人的牙齿排列也有很大的差异，因此单颗牙齿的分割即使在牙科计算机软件中也是一项较为耗时的工作。

近几年来，国内外许多研究人员对单颗牙齿的自动分割进行了研究。这些方法主要分为两种类型：基于2维图像的分割方法和基于3维模型的分割方法。第一种类型最为典型的代表是Kondo等人^[3]提出的深度图像分割法：首先从咬合平面方向上获取牙颌模型的深度图像，并通过计算该图像的牙弓曲线获得全景深度图像，然后结合两幅图像的边缘特征得到牙齿分割轮廓；然而，由于牙缝特征点是通过做牙弓曲线的垂线来获得，当患者存在畸形牙的时候，这种方法可能会产生较大的分割偏差。Grzegorzek等人^[4]提出一种基于多层深度图像分割方法，避免了仅依靠单一图像而产生精度差的问题，该方法首先探测每一层图像的牙齿外轮廓和牙齿之间的非凸顶点，然后连接这些非凸顶点并映射回3维空间当中获得牙齿分割边界；然而这些深度图像也仅仅是从咬合平面视角方向获得，当患者存在不整齐的牙齿时，这种方法也不能产生较好的分割结果。

由于基于2维图像的分割无法准确表达牙齿的真实形态，因此基于3维模型的分割方法越来越受到研究学者的关注。Kronfeld等人^[5]结合snake活动轮廓模型和边缘相切流的方法实现了牙齿的自动分割，然而由于相切流容易受到噪声的影响，因此分割精度较低。Yuan等人^[6]首先删除齿间粘连区域，然后采用曲面能量约束方式恢复缺失的开口区域以获得单颗牙齿完整形态，最后利用骨架线方法获得每颗牙齿分割边界，但该方法在识别粘连区域时需要大量的手动交互。Kumar等人^[7]首先利用曲率分析和漫水填充(flood fill)算法探测牙齿和牙龈的分割边界，然后探测此边界上的角点来获取牙缝分割边界；但漫水填充算法在噪音和干涉分支的干扰下容易产生欠分割或过分割的现象，该方法还需要复杂的路径搜索算法进行修复。Wu等人^[8]提出了一种改进的形态学骨架分割算法较为精确的获取牙齿分割边界，但骨架线算法容易产生较多的骨刺和干涉分支。Zou等人^[9]和Wang等人^[10]提出一种基于调和场的牙齿分割方法，不仅能够处理严重畸形的牙齿，还能够处理牙齿边缘模糊的情况，然而该方法需要手动设定每颗牙齿的分割约束点，因此人工干预较多。

针对上述问题，本文提出一种新颖的利用路径规划技术的单颗牙齿自动分割方法，以利用路径搜索强大的抗干扰抗断裂能力最大程度地减少人工干预。该方法同时结合牙颌模型的微分特征、图像形态学以及B样条拟合技术，以确保分割的精度和稳定性。

由于三角网格牙颌模型上牙齿和牙龈以及牙齿和牙齿之间的边界往往交织在一起，因此，为避免在探测分割边界时产生相互干扰，首先规划牙齿和牙龈之间的分割路径，然后对牙齿之间的牙缝分割路径进行规划，具体路线如下：

1) 牙龈分割路径规划(图 1)。首先利用模型的曲率特征提取牙齿和牙龈分界的特征区域，并基于一种双重路径规划法搜索牙龈分割路径，然后对得到的路径利用图像形态学和B样条延伸拟合方法构建牙弓曲线。

2) 牙缝分割路径规划(图 2)。综合牙龈路径和牙弓曲线的形态特征搜索牙龈路径上的凹点，以确定牙齿之间的缝隙边界点，然后基于路径规划法搜索牙缝分割路径，最后细化路径以得到每颗牙齿精确的分割轮廓。

牙齿和牙龈之间的分界特征区域一般位于牙颌模型表面曲率较大的地方，因此首先采用局部曲面拟合法^[11]对模型进行离散曲率分析，并基于最大主曲率原则获取牙颌模型上每个顶点$\mathit{\boldsymbol{p}}$的曲率值$\kappa \left( \mathit{\boldsymbol{p}} \right)$。为提高曲率可视化的效果，突出模型的凹凸特征，将曲率值映射到[0, 1]区间，然后对曲率值进行基于直方图均衡化的拉伸变换^[12]，变换后的曲率映射如图 3(a)所示。

对变换后的曲率值进行阈值操作{$\kappa \left( \mathit{\boldsymbol{p}} \right) > h$($h$为阈值)}，得到初始特征区域，图 3(b)所示。区域中的一些杂点和断点利用3维形态学开闭运算方法^[13]处理。剩余的一些干涉分支，可通过连通性进行过滤，即如果两个特征点$\mathit{\boldsymbol{p}}$, $\mathit{\boldsymbol{q}}$在一条网格边上或存在一条路径$\mathit{\boldsymbol{p}} \to {\mathit{\boldsymbol{v}}_1} \to {\mathit{\boldsymbol{v}}_2} \to \cdots \to {\mathit{\boldsymbol{v}}_i} \to \mathit{\boldsymbol{q}}$(${\mathit{\boldsymbol{v}}_i} \in \mathit{\boldsymbol{V}}$, $\mathit{\boldsymbol{V}}$为模型顶点集合)，则这两个点属于同一个连通分量。由于分支区域内点的数量较少，因此利用这个属性删除这些干涉连通分量，得到最终的牙齿和牙龈分界的特征区域$\mathit{\boldsymbol{F}}$，如图 3(c)所示。

路径规划方法是一种基于图论的节点路径搜索算法，其最大的特点是在获取最优分割线的同时，避免分割线断裂和干涉分支。但由于其搜索方向无法控制，直接进行路径搜索，容易使路径从牙缝处绕行，导致不正确的结果，如图 4所示。因此，采用一种基于双重规划的牙龈路径搜索方法，以确保分割结果的稳定性。

路径规划需要首先确定路径的起点和终点，为减少人工干预，利用坐标转换法进行自动搜索。将特征区域$\mathit{\boldsymbol{F}}$内的顶点由直角坐标转化为柱面坐标，然后对柱面坐标中的极角进行排序，则起始节点$\mathit{\boldsymbol{m}}$和终止节点$\mathit{\boldsymbol{n}}$分别定义为这些角度差别中最大的两个点，如图 5(a)所示。

路径搜索的基本方法是将数据构成一个带权的图，图中两个节点之间边的权值表示相应的代价，最优路径就是从起点到终点中所有边的代价之和最小的那条路径^[14]。根据特征区域$\mathit{\boldsymbol{F}}$的分布情况，首先将此区域数据构建一个无向连通图$\mathit{\boldsymbol{G}}$，为使路径搜索不从牙缝绕行，将图$\mathit{\boldsymbol{G}}$中边的权值代价函数定义为

$ {w_1}\left( {{\mathit{\boldsymbol{v}}_i}, {\mathit{\boldsymbol{v}}_j}} \right) = \left\{ \begin{array}{l} \left\| {{\mathit{\boldsymbol{v}}_i} - {\mathit{\boldsymbol{v}}_j}} \right\|\;\;\;\;\;\left( {{\mathit{\boldsymbol{v}}_i}, {\mathit{\boldsymbol{v}}_j}} \right) \in \mathit{\boldsymbol{E}}\\ \infty \;\;\;\;\;\;\;\;\;\;\;\;\;其他 \end{array} \right. $

(1)

式中，$\left\| {{\mathit{\boldsymbol{v}}_i} - {\mathit{\boldsymbol{v}}_j}} \right\|$代表两个顶点的欧氏距离，$\mathit{\boldsymbol{E}}$代表网格模型上所有边的集合。在此函数定义下，边的权值与两个邻接点之间的距离成正比。通过利用Dijkstra算法^[15]进行求解，可得到从点$\mathit{\boldsymbol{m}}$到点$\mathit{\boldsymbol{n}}$的最短路径，如图 5(a)所示的绿色线。可以看出，该方法搜索到了舌侧的牙龈分割路径${P_{{\rm{in}}}}$，并避免了路径从牙缝以及其他干涉分支处绕行的问题。为进一步搜索颊侧的牙龈分割路径，将搜索到的舌侧牙龈路径邻域边的权值赋为无穷大，即断开通过舌侧的路线，再次搜索从点$\mathit{\boldsymbol{m}}$到点$\mathit{\boldsymbol{n}}$的最优路径，即为颊侧牙龈线${P_{{\rm{out}}}}$，如图 5(b)所示的蓝色线。

值的注意的是，舌侧牙龈路径搜索后，需更新路径的起点和终点，如图 5(c)所示，因为利用坐标转换法得到的起点和终点可能会偏离所需的位置。

由于第1次搜索路径为避开牙缝，得到的路径仅仅是“最短的”，并没有很好地贴合在牙齿和牙龈分界最为凹陷的地方，因此利用模型的离散曲率信息对路径进行修复。首先收集上一步搜索得到的舌侧牙龈路径${P_{{\rm{in}}}}$的邻域集合，然后根据顶点连接关系构建一个无向连通图${\mathit{\boldsymbol{G}}_{{\rm{in}}}}$，图中边的权值函数定义为

$ {w_2}\left( {{\mathit{\boldsymbol{v}}_i}, {\mathit{\boldsymbol{v}}_j}} \right) = \left\{ \begin{array}{l} \frac{{\left\| {{\mathit{\boldsymbol{v}}_i} - {\mathit{\boldsymbol{v}}_j}} \right\|}}{{\kappa \left( {{\mathit{\boldsymbol{v}}_\mathit{i}}} \right)\kappa \left( {{\mathit{\boldsymbol{v}}_\mathit{j}}} \right)}}\;\;\;\;\;\left( {{\mathit{\boldsymbol{v}}_i}, {\mathit{\boldsymbol{v}}_j}} \right) \in \mathit{\boldsymbol{E}}\\ \infty \;\;\;\;\;\;\;\;\;\;\;\;\;其他 \end{array} \right. $

(2)

式中，${\kappa \left( {{\mathit{\boldsymbol{v}}_\mathit{i}}} \right)}$和${\kappa \left( {{\mathit{\boldsymbol{v}}_\mathit{j}}} \right)}$分别为顶点${{\mathit{\boldsymbol{v}}_i}}$和${{\mathit{\boldsymbol{v}}_j}}$的曲率值。该函数意义在于，当两个邻接点具有较高的曲率值以及较近的距离，其边的代价就越小，以吸引最优路径走到分界最为明显的节点上。求解图${\mathit{\boldsymbol{G}}_{{\rm{in}}}}$中点$\mathit{\boldsymbol{m}}$到点$\mathit{\boldsymbol{n}}$的最优路径，即可得到修正后的舌侧牙龈路径${{P'}_{{\rm{in}}}}$。用同样的方法，构造颊侧的连通图${\mathit{\boldsymbol{G}}_{{\rm{out}}}}$，可以得到修正后的颊侧牙龈路径${{P'}_{{\rm{out}}}}$。图 5(d)对比显示了牙龈路径修正前后的效果，可以看出，修正后的路径较好地贴合在牙齿和牙龈分界最为凹陷的位置。

牙弓曲线是医学上描述牙齿在牙槽骨上排列的形态。在口腔正畸中，牙弓曲线是评定牙齿咬合不正的一个重要参考标准。目前，用来描述牙弓曲线的几种数学形式，主要有多项式曲线，悬链曲线、beta函数曲线以及B样条曲线等，由于在实际的操作过程中，医生可能会对患者的牙弓曲线进行修改^[16]，因此利用B样条曲线的局部修改性质进行调整比较方便。但传统方法进行B样条曲线拟合时，往往需要人工交互选择特征型值点，本文提出一种基于图像形态学和B样条延伸拟合的方法对搜索到的牙龈分割路径自动创建牙弓曲线。

首先将搜索到的牙龈分割路径投影到咬合平面上，构建一个牙龈路径二值图像(图 6(a))，然后对该二值图像进行形态学填充(图 6(b))以及细化操作，以获取牙弓骨架线(图 6(c))。

将牙弓骨架线作为型值点数据$\left\{ {{\mathit{\boldsymbol{Q}}_i}, i = 1, 2, \cdots, m} \right\}$，拟合为B样条曲线$C\left( s \right)$。为使拟合后的曲线长度能够覆盖整个牙颌模型范围，提出一种延伸拟合的方法^[17]。假设拟合后的B样条曲线在型值点数据的起始点两端各自延长比例分别为${\alpha _1}$、${\alpha _m}$，$L$为型值点数据原始总弦长, 即

$ L = \sum\limits_{i = 2}^m {\left\| {{\mathit{\boldsymbol{Q}}_i} - {\mathit{\boldsymbol{Q}}_{i - 1}}} \right\|} $

(3)

则延伸后的总长可表示为

$ {L_{{\rm{extend}}}} = L\left( {1 + {\alpha _0} + {\alpha _m}} \right) $

(4)

因此，B样条曲线$C\left( s \right)$的参数$s$，可计算为

$ \left\{ \begin{array}{l} {s_1} = \frac{{{\alpha _0}}}{{1 + {\alpha _0} + {\alpha _m}}}\\ {s_i} = {s_{i - 1}} + \frac{{\left\| {{\mathit{\boldsymbol{Q}}_i} - {\mathit{\boldsymbol{Q}}_{i - 1}}} \right\|}}{{{L_{{\rm{extend}}}}}}\;\;\;\;\;\;i = 2, 3, \cdots, m \end{array} \right. $

(5)

通过最小化能量函数

$ E = \sum\limits_{i = 1}^m {\left\| {\mathit{\boldsymbol{C}}\left( {{s_i}} \right) - {\mathit{\boldsymbol{Q}}_i}^2} \right\|} + \lambda \sum\limits_{i = 1}^m {{{\left\| {\mathit{\boldsymbol{C''}}\left( {{s_i}} \right)} \right\|}^2}} \Delta s $

(6)

即可获得延伸拟合后的牙弓B样条曲线$C\left( s \right)$，如图 6(d)所示。

牙齿之间的缝隙位置一般位于牙龈路径的凹点处，如图 7所示，因此在牙缝路径规划前，需先进行凹点探测。由于牙龈路径波浪起伏不定，这里提出一种综合牙弓曲线和牙龈路径形态特征的自动探测方法。

曲线上的凹点和凸点一般都具有较大的曲率，因此首先将得到的舌侧牙龈路径${{P'}_{{\rm{in}}}}$和颊侧牙龈路径${{P'}_{{\rm{out}}}}$分别拟合为光滑曲线${\mathit{\boldsymbol{f}}_{{\rm{in}}}}\left( t \right)$和${\mathit{\boldsymbol{f}}_{{\rm{out}}}}\left( t \right)$，并通过计算曲线的曲率得到牙龈路径上各点的曲率，如图 8(a)(b)所示。

仔细观察牙龈曲线上的凹点特性，可以发现，凹点的法向量方向与凹点距离牙弓曲线上的最近点和该凹点连线的方向一致(图 8(c)(d))，而法向量可以利用二阶导数确定，因此可利用约束条件进行凹点探测，即

$ \begin{array}{l} \mathit{\boldsymbol{f}}_{{\rm{in}}}\left( t \right) = \left\{ \begin{array}{l} 凹点\;\;\;\;\;\;\;\;{{\mathit{\boldsymbol{f''}}}_{{\rm{in}}}}\left( t \right) \cdot \left( {{\mathit{\boldsymbol{f}}_{{\rm{in}}}}\left( t \right) - \mathit{\boldsymbol{C}}\left( s \right)} \right) > 0\;\& \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rho _{{\rm{in}}}}\left( t \right) \ge T\\ 非凹点\;\;\;\;\;\;\;\;\;\;其他 \end{array} \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\mathit{\boldsymbol{f}}_{{\rm{out}}}}\left( t \right) = \\ \left\{ \begin{array}{l} 凹点\;\;\;\;\;\;\;\;{{\mathit{\boldsymbol{f''}}}_{{\rm{out}}}}\left( t \right) \cdot \left( {{\mathit{\boldsymbol{f}}_{{\rm{out}}}}\left( t \right) - \mathit{\boldsymbol{C}}\left( s \right)} \right) > 0\;\& \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rho _{{\rm{out}}}}\left( t \right) \ge T\\ 非凹点\;\;\;\;\;\;\;\;\;\;其他 \end{array} \right. \end{array} $

(7)

式中，${\mathit{\boldsymbol{C}}\left( s \right)}$为牙弓曲线上距离待测牙龈曲线最近的参数点，${{\mathit{\boldsymbol{f''}}}_{{\rm{in}}}}\left( t \right)$和${{\mathit{\boldsymbol{f''}}}_{{\mathit{out}}}}$分别为舌侧牙龈曲线${\mathit{\boldsymbol{f}}_{{\rm{in}}}}\left( t \right)$和颊侧牙龈曲线${\mathit{\boldsymbol{f}}_{{\rm{out}}}}\left( t \right)$的二阶导数，“·”为两个向量的点积；${\rho _{{\rm{in}}}} $和${\rho _{{\rm{out}}}} $为相应牙龈曲线的曲率，$T$为曲率阈值.

将探测到的凹点(图 8(e))删除，并将剩余的舌侧和颊侧牙龈路径上的点按顺序进行组合，获得每颗牙齿的牙龈分割线，如图 8(f)所示。

由于牙缝位于模型表面相邻牙齿之间的凹状区域，因此也利用路径规划方法进行搜索。将特征区域$\mathit{\boldsymbol{F}}$构建成无向连通图${\mathit{\boldsymbol{G}}_{\rm{2}}}$，${\mathit{\boldsymbol{G}}_{\rm{2}}}$中边的权值按基于距离和曲率的式(2)定义。然后，利用Dijkstra算法分别搜索每颗牙的舌侧牙龈分割线的一个端点到对应颊侧牙龈分割线端点的最优路径，以及舌侧牙龈分割线的另一个端点到对应颊侧牙龈线另一个端点的最优路径，从而获得每颗牙齿两侧的牙缝路径，如图 9所示。

利用路径规划方法得到的路径点都只是网格模型上的顶点，其路径的光滑性取决于网格的密度。因此，本文系统提供给用户一个额外的选择：细化路径。目前，关于细化分割线的方法大都采用插值^[8]或网格细分^[18]等方法，然而网格细分容易导致分割轮廓偏离原始路径，而插值法需要将插值点投影到网格上才能获得正确的结果, 且投影求交繁琐。本文采用一种的简单而有效的策略^[19]以对获得的原始分割路径进行细化。

对于每颗牙齿封闭的路径，首先搜索其邻域网格，将路径邻域内的每一条网格特征边进行对分，使对分后的每段边长小于给定的值${d_{{\rm{min}}}}$。每个对分顶点的坐标${\mathit{\boldsymbol{v}}_x}$和曲率值$\kappa \left( {{\mathit{\boldsymbol{v}}_x}} \right)$根据其所在的网格边两端的原始顶点坐标及其曲率值插值确定，即

$ \begin{array}{l} \left\{ \begin{array}{l} {\mathit{\boldsymbol{v}}_x} = {\mathit{\boldsymbol{v}}_i} + x\frac{{\left\| {{\mathit{\boldsymbol{v}}_i} - {\mathit{\boldsymbol{v}}_j}} \right\|}}{{N + 1}}\\ \kappa \left( {{\mathit{\boldsymbol{v}}_x}} \right) = \kappa \left( {{\mathit{\boldsymbol{v}}_i}} \right) + x\frac{{\left| {\kappa \left( {{\mathit{\boldsymbol{v}}_i}} \right) - \kappa \left( {{\mathit{\boldsymbol{v}}_j}} \right)} \right|}}{{N + 1}} \end{array} \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;x = 1, 2, \cdots, N\\ \;\;\;\;\;\;\;\;\;\;\;N = floor\left( {\frac{{{\mathit{\boldsymbol{v}}_i} - {\mathit{\boldsymbol{v}}_j}}}{{{d_{{\rm{min}}}}}}} \right) \end{array} $

(8)

式中，$N$代表该网格边上对分顶点的个数，$floor\left( {} \right)$为向下取整函数，${{\mathit{\boldsymbol{v}}_i}}$和${{\mathit{\boldsymbol{v}}_j}}$为该特征边两端的顶点坐标，${\kappa \left( {{\mathit{\boldsymbol{v}}_i}} \right)} $和${\kappa \left( {{\mathit{\boldsymbol{v}}_j}} \right)}$为其相应的曲率值。

将细分后的所有顶点构建无向连通图${\mathit{\boldsymbol{G}}_{{\rm{refine}}}}$，以利用路径规划方法搜索精确的最终路径。为避免路径线偏离网格模型表面，在进行对分顶点时，需记录对分顶点所在的网格边，以方便判断图中两个点${{\mathit{\boldsymbol{v}}_x}}$, ${{\mathit{\boldsymbol{v}}_y}}$是否属同一三角片。因此，图${\mathit{\boldsymbol{G}}_{{\rm{refine}}}}$中两个顶点边的权值定义为

$ {w_3}\left( {{\mathit{\boldsymbol{v}}_x}, {\mathit{\boldsymbol{v}}_y}} \right) = \left\{ \begin{array}{l} \frac{{\left\| {{\mathit{\boldsymbol{v}}_x} - {\mathit{\boldsymbol{v}}_y}} \right\|}}{{\kappa \left( {{\mathit{\boldsymbol{v}}_x}} \right)\kappa \left( {{\mathit{\boldsymbol{v}}_y}} \right)}}\;\;\;{\mathit{\boldsymbol{v}}_x}, {\mathit{\boldsymbol{v}}_y}\ 属于同一三角片\\ \infty \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;其他 \end{array} \right. $

(9)

然后，在此连通图内，基于每颗牙齿的牙齿分割轮廓，再次利用Dijkstra算法搜索精确的分割路径。图 10显示路径细化前后的对比结果，可以看出，细化后的路径更加光滑，同时也很好地贴合在边界较明显的位置。

为验证本文方法的有效性，对不同年龄、不同程度畸形牙患者的牙颌模型进行测试(图 11，图 12)。患者模型利用高精度3维激光扫描仪进行数据采集，获取数字化三角网格牙颌模型，网格顶点间距设为0.15 mm。实验在CPU为2.3 GHz、内存为8 GB的微机上运行，软件开发平台为Matlab R2015b。

由于篇幅所限，图 11列出了实验中5种典型的牙颌模型的分割结果，并按行显示了具体的分割过程。这5个模型的牙齿畸形程度依次递增。可以看出，尽管在特征提取阶段，一些模型的特征区域包含很多牙沟牙窝以及噪音等干涉区域，但是在牙龈路径搜索过程中，仍然获得了较好的结果，这是因为路径规划具有强大的抗干扰能力，能够避免产生干涉分支；同时由于我们提出的双重路径规划方法，较好的避免了路径的绕行问题。另外，也可以发现，即使在像模型4和模型5这样牙齿畸形较严重的情况下，本文方法也能产生正确的结果。

在整个实验中，大部分过程基本都是自动的，需要调整的参数较少，除了一些模型需要修改默认的曲率阈值。但是，当存在牙齿分割边界较为模糊的情况下，例如像如图 12显示的模型，由于牙齿错位造成相邻牙齿分界位置较为不明显，因此在牙缝探测过程中可能会出现图 12中所示的遗漏现象，需要手工干预调整。但相比于3shape软件以及其他文献(如文献[9]设置每颗牙齿周围的约束点)的方法，本文的手工调整方便而简单，能够同时保证分割结果的效率和准确性。

为方便比较算法消耗时间，将整个系统的操作过程分为4个部分：“特征区域提取”、“牙龈分割路径规划”、“牙缝分割路径规划”和“分割路径细化”。各部分耗时如表 1所示。其中，“特征区域提取”包含离散曲率分析、阈值提取和连通性过滤，“牙龈分割路径规划”包含牙龈路径探测和牙弓曲线拟合，“牙缝分割路径规划”则包含凹点探测和牙缝路径搜索，最后的“分割路径细化”包含网格对分和牙齿切割。

从表 1可以看出，牙龈分割路径规划和牙缝分割路径规划都各自只需不到1 s的时间即可完成，相比于其他文献方法，大大降低了搜索时间。本文算法在特征区域提取和分割路径细化上耗时较多，其中特征区域提取平均需耗时约为5~6 s，并且其时间随网格模型顶点数量的增多而增加，该过程需要对网格的每个顶点进行局部曲面拟合以计算曲率，因此耗时较长；其次，分割路径细化耗时大约为4~5 s，这主要是因为细化路径需要花费较多的时间对原始分割路径周围的网格进行对分，并判断这些细分顶点之间是否属于同一个三角片，因此耗时较多。每个模型在整个分割过程中，加上打开模型和设置参数，基本控制在20 s以内。即使特殊情况加上手动干预，本系统也能在30 s之内完成一个牙齿模型的分割。而文献[8-9]都至少需要1~2 min的时间；著名的“3shape”软件则更是需要数分钟时间的手动交互才能完成分割。

本文针对3维牙颌解剖形态和排列特点，提出一种新颖的基于路径规划技术的单颗牙齿自动分割方法。该方法采用牙龈分割路径和牙缝分割路径分开搜索的策略，避免了牙龈和牙缝分割线的相互干扰；提出的双重路径规划法避免了传统分割路径容易产生的分割线绕行以及断裂问题，并有效解决牙沟和牙缝干涉分支这一难点；结合图像形态学和B样条拟合的单颗牙齿路径划分方法，避免了传统方法导致误分以及人工干预较多的问题。与现有的方法相比，本文方法简单快速，而且具有较少的人工干预和参数调整，适用于各种畸形牙患者的牙齿分割，对后续正畸治疗中的虚拟牙排列具有重要应用价值。但对于牙齿边界较为模糊的情况，本方法还需要一些手动调整，如何进一步地减少手动干预、提高分割的精度，将是后续工作的重点。