CN105678747A - 一种基于主曲率的牙齿网格模型自动分割方法 - Google Patents

Abstract

本发明属于三维图像处理技术领域，具体涉及一种基于主曲率的牙齿网格模型自动分割方法，将牙齿三角网格模型视为曲面计算网格的主曲率，然后利用主曲率与区域生长法分割出单个牙齿的牙齿曲面。本发明的基于主曲率的牙齿网格模型自动分割方法，能够避开复杂的模糊聚类，将牙齿三角网格模型视为曲面计算网格的主曲率，利用主曲率计算方法与区域生长法分割单个牙齿。能简化医生操作、让单个牙齿模型清晰地展现，帮助进行牙齿的全面检查和矫正。

Description

一种基于主曲率的牙齿网格模型自动分割方法

技术领域

本发明属于三维图像处理技术领域，具体涉及一种基于主曲率的牙齿网格模型自动分割方法。

背景技术

三维网格模型分割是根据一定的几何及拓扑特征将由网格组成的多面体或可定向的二维流形分解为一组数目有限、各自具有简单形状意义、各自连通的子部位的过程。多数三维网格分割算法是在二维图像分割的基础上，推广到三维网格空间，如：分水岭算法、K-means、Mean-shift、区域增长法等。在分割中引入几何与拓扑信息会提高分割效果，如：基于高斯曲率、基于体元、基于Reeb图、基于骨架提取和拓扑结构扫描等使用三维网格曲面形状特征的算法。

加拿大的西蒙弗雷泽大学计算机科学学院提出了一种基于谱聚类的3D网格分割方法。给出一组由三角网格构造的人脸数据集，为属于同组的每一对人脸进行似然编码以此定义相似度矩阵；利用距离矩阵避免脸部属于不同部位的凹凸区域被误分入同一组；随后计算相似度矩阵的前k个特征向量，使用这些特征向量把人脸网格映射到一个有k维的球上，这样就可以直接对球上面的点做k-means聚类，最终根据聚类结果进行三维分割。同时提出了一种启发式算法为聚类选择一个合适的起点，以此提高收敛速度并且避免局部最小值。此方法需要前期训练，计算复杂度较高，有时无法获得很平坦的分割边界或匹配过程较耗时，在精确度上有一定提升空间。

江南大学、浙江理工大学、北京航空航天大学共同提出一种基于凸凹信号的网格分割算法。将三维图形的某顶点分为平坦点、凸点、凹点和特征点四类，从平坦点扩展平坦区域后，从凸凹点出发进行凸凹区域扩展,最后对未缝合的区域进行基于边界点和边界边光滑度的区域竞争扩展。对于简化程度高的部位，很难寻找到平坦点或凹凸点，因此直接计算面的凹凸信号进行最后分割。但此方法会出现凹凸区域在边界处不完全缝合的状况，因此该方法最后对于不贴合部位进行区域合并，并且区域合并没有明确的标准，对不同的物体、甚至是同物体的不同区域需要做不同调整。此方法比较适合部位分割明显的物体，比如汽车、茶壶等，容易把同一颗牙齿上的凹凸分割成几部分。

上海交通大学的研究人员认为，牙齿的分割过程实质上是模型中三角网格上每个顶点在内力和外力的作用下不断运动的过程。首先，通过图像灰度信息得到抑制图像中过亮或过暗的体素点所产生影响的灰度值，然后通过观察牙齿截面图中牙齿的形状及角度，估计出初始模型的中心位置在x，y，z三个方向上的半径及旋转角度，并利用这些数据建立初始模型；最后，将初始模型放到牙齿的内部，由内力和外力同时控制模型向牙齿边缘迭代运动，使模型在演化过程中始终保持表面平滑及表面三角片的均匀分布。他们还定义了外力的权重系数，此系数不但使顶点的每步运动距离小于相邻顶点间的平均距离，而且要保证外力能够反映该处图像灰度信息，并不能影响模型的平滑性。但此方法的难点在于选择权重系数，不同的牙齿几何形状的多边性会导致三角网格模型在合外力的作用下收敛不够准确快速。

目前分割时容易出现过分割、边缘不够精确、需要大量人工干预、参数难以选择、或算法复杂度高等问题。如今很难有一种分割方法适用所有应用，大多数算法都是针对具体问题提出的。

发明内容

本发明提出一种基于主曲率的牙齿网格模型自动分割方法，本方法将牙齿三角网格模型视为曲面计算网格的主曲率，然后利用主曲率与区域生长法分割出单个牙齿的牙齿曲面。

为了实现上述目的，本发明采用如下技术方案：一种基于主曲率的牙齿网格模型自动分割方法，包括如下步骤：

(1)计算牙齿网格模型在每个顶点处的主曲率

(2)根据主曲率对牙齿网格模型进行初步分割；

(3)通过区域生长法初步分割出牙齿齿面；

(4)根据联通包围性找回牙齿咬合面的孤立区域，形成完整牙齿齿面。

其中，所述步骤(2)中根据主曲率对牙齿网格进行初步分割的步骤为：

(21)针对牙齿网格模型任意顶点P，计算P的主曲率绝对值；

(22)计算P的主曲率绝对值的最小值

(23)若则P点为平坦点，加入平坦点集Ps，否则P点为弯曲点，加入弯曲点集Pc，其中t为设定的阈值。

其中，所述步骤(3)中通过区域生长法初步分割出牙齿齿面的步骤为：

(31)人工选取牙齿齿面上任意顶点N，N属于平坦点集Ps，然后将放入待处理顶点集Pf；

(32)从待处理顶点集Pf取出一个顶点M，将M放入牙齿顶点集Pt；

(33)获得顶点M的1-ring邻域顶点Mn；

(34)若Mn属于平坦点集Ps，则将Mn放入待处理顶点集Pf；

(35)若待处理顶点集Pf不为空，则跳到步骤(32)，否则退出，得到牙齿顶点集Pt。

本发明的基于主曲率的牙齿网格模型自动分割方法，能够避开复杂的模糊聚类，将牙齿三角网格模型视为曲面计算网格的主曲率，利用主曲率计算方法与区域生长法分割单个牙齿。能简化医生操作、让单个牙齿模型清晰地展现，帮助进行牙齿的全面检查和矫正。

附图说明

图1为本发明牙齿网格模型自动分割处理流程图；

图2为本发明主曲率示意图；

图2a为峰点(谷点)处的主曲率；

图2b为鞍点处的主曲率；

图2c为平坦点处的主曲率；

图3为本发明的牙齿咬合面补齐示意图；

图3a为区域生长初始分割结果示意图；

图3b为牙齿咬合面补齐结果示意图。

具体实施方式

下面结合具体实施方式对本发明进行详细的描述。

本发明的基于主曲率的牙齿网格模型自动分割方法是针对牙齿模型的三维扫描图像进行分割，从而将牙齿模型三维扫描图像上的单颗牙齿的牙齿曲面分割下来。由于计算机不能直接表达曲面，因此计算机表达的牙齿模型的三维扫描图像是三角网格状的牙齿模型，也就是牙齿网格模型，该牙齿网格模型是由若干三角面片拼接而成的。

一、计算牙齿网格模型在每个顶点处的主曲率

本发明的方法需要先求取牙齿网格模型在每个顶点处的主曲率，即三角面片的顶点处的主曲率。下面对网格主曲率的理论基础和计算方法进行介绍：

1.主曲率

对某曲面S上的一点p，设T为点p上的单位切向量，取一对正交单位向量{T₁,T₂}，则p上沿T方向的法曲率为：

${\stackrel{\‾}{k}}_p\left(T\right)={\left(\begin{array}{c}{t}_1\\ t_2\end{array}\right)}^{\′}\left(\begin{array}{cc}{k}_p^{11}& k_p^{12}\\ k_p^{21}& k_p^{22}\end{array}\right)\left(\begin{array}{c}{t}_1\\ t_2\end{array}\right)---\left(1\right)$

当时，{T₁,T₂}为曲面S在点p的主方向，相应的为主曲率。在正交基{T₁,T₂}中加入法向量N得：{N,T₁,T₂}，因此，法曲率也可表示为：

${\stackrel{\‾}{k}}_p\left(T\right)={\left(\begin{array}{c}n\\ t_1\\ t_2\end{array}\right)}^{\′}\left(\begin{array}{ccc}0& 0& 0\\ 0& k_p^1& 0\\ 0& 0& k_p^2\end{array}\right)\left(\begin{array}{c}n\\ t_1\\ t_2\end{array}\right)---\left(2\right)$

其中，T＝nN+t₁T₁+t₂T₂是任意向量。

引入-π≤θ≤+π，则单位切向量可表示为T_θ＝T₁cosθ+T₂sinθ，法曲率可表示为：

$k_p\left(T_{\mathrm{\θ}}\right)=k_p^1{\cos }^2\mathrm{\θ}+k_p^2{\sin }^2\mathrm{\θ}---\left(3\right)$

其中，是{T₁,T₂}方向上的主曲率，且θ是与主方向T₁,T₂的夹角。

定义一对称矩阵：

$M_P=\frac{1}{2\mathrm{\π}}{\∫}_{-\mathrm{\π}}^{\mathrm{\π}}{k}_p\left(T_{\mathrm{\θ}}\right)T_{\mathrm{\θ}}{T}_{\mathrm{\θ}}^{1}d\mathrm{\θ}---\left(4\right)$

由于是秩为1的矩阵，因此法向量N是该矩阵特征值为0时对应的特征向量。然后将M_P分解：

$M_P=\left(\begin{array}{cc}{T}_1& T_2\end{array}\right)\left(\begin{array}{cc}\frac{1}{2\mathrm{\π}}{\∫}_{-\mathrm{\π}}^{\mathrm{\π}}{k}_p\left(T_{\mathrm{\θ}}\right){\cos }^2\mathrm{\θ}d\mathrm{\θ}& \frac{1}{2\mathrm{\π}}{\∫}_{-\mathrm{\π}}^{\mathrm{\π}}{k}_p\left(T_{\mathrm{\θ}}\right)\cos \mathrm{\θ}\sin \mathrm{\θ}d\mathrm{\θ}\\ \frac{1}{2\mathrm{\π}}{\∫}_{-\mathrm{\π}}^{\mathrm{\π}}{k}_p\left(T_{\mathrm{\θ}}\right)\sin \mathrm{\θ}\cos \mathrm{\θ}d\mathrm{\θ}& \frac{1}{2\mathrm{\π}}{\∫}_{-\mathrm{\π}}^{\mathrm{\π}}{k}_p\left(T_{\mathrm{\θ}}\right){\sin }^2\mathrm{\θ}d\mathrm{\θ}\end{array}\right)\left(\begin{array}{c}{T}_1\\ T_2\end{array}\right)---\left(5\right)$

令由列向量和组成的矩阵，则可以重写上述M_P矩阵：

$M_P=T_{12}\left(\begin{array}{cc}{m}_p^{11}& m_p^{12}\\ m_p^{21}& m_p^{22}\end{array}\right){T_{12}}^{\′}---\left(6\right)$

其中，。

除了特征向量N之外，M_P的另外两个特征向量是主方向和，但其相应的特征值不是主曲率。特征值与主对角线元素的关系如下：

$m_p^{11}=\frac{3}{8}{k}_p^1+\frac{1}{8}{k}_p^2---\left(7\right)$

$m_p^{22}=\frac{1}{8}{k}_p^1+\frac{3}{8}{k}_p^2---\left(8\right)$

因此可得主曲率是对称矩阵M_P的非零特征值函数：

$k_p^1=3m_p^{11}-m_p^{22}---\left(9\right)$

$k_p^2=3m_p^{22}-m_p^{11}---\left(10\right)$

于是，只要确定M_P中的便可求得主曲率。可由M_P的特征向量和特征值求解得主曲率和主方向。对于单位切向量T的方向采用如下公式计算

${\stackrel{\&OverBar;}{k}}_p\left(T\right)=\underset{s\&RightArrow;0}{\lim }\frac{2<N,\mathrm{\&alpha;}\left(s\right)-p>}{\left|\right|\mathrm{\&alpha;}\left(s\right)-p|{|}^2}---\left(11\right)$

对于三角网络上的顶点p_i，顶点p_j是顶点p_i的1-ring邻域上的点，沿着切方向，最终采用此公式估计顶点p_i在方向T_ij上的法曲率：

$k_{ij}=\frac{2<N_{p_i},p_j-p_i>}{\left|\right|p_j-p_i|{|}^2}---\left(12\right)$

2.主曲率离散计算步骤

(1)估计三角网络上顶点的法向量。

由于每个三角形面片f_k都定义的单位长度的法向量N_fk。曲面是有方向的，对与当前点v_i的法向量Taubin给出的面积加权和的定义：

$N_{v_i}=\frac{{\mathrm{\&Sigma;}}_{f_k\&Element;F^{\&prime;}}\left|f_k\right|N_{f_k}}{\left|\right|{\mathrm{\&Sigma;}}_{f_k\&Element;F^{\&prime;}}\left|f_k\right|N_{f_k}\left|\right|}---\left(13\right)$

注意：原始的Taubin算法才面积加权来求顶点的法线，我们也可以才采用角度加权和来求顶点的法线。

(2)估计矩阵

估计矩阵定义为顶点v_i的1-ring邻点面积和沿着方向T_ij曲率的乘积的加权和，具体计算积分矩阵M_p如下：

$M_{v_i}={\mathrm{\&Sigma;}}_{v_j\&Element;v^{\&prime;}}{w}_{ij}{k}_{ij}{T}_{ij}{T}_{ij}^{\&prime;}---\left(14\right)$

对于顶点v_i的每个1-ring邻点，定义单位向量T_ij是向量在切平面(平面法线为)上的投影，即

$T_{ij}=\frac{(I-N_{v_i}{N_{v_i}}^{\&prime;})\stackrel{\&OverBar;}{v_i-v_j}}{\left|\right|(I-N_{v_i}{N_{v_i}}^{\&prime;})\stackrel{\&OverBar;}{v_i-v_j}\left|\right|}---\left(15\right)$

沿着方向T_ij曲率估计方法为公式(12)。权值w_ij为面积加权。

(3)矩阵的Householder变换

令 $E=\left(\begin{array}{c}1\\ 0\\ 0\end{array}\right),W_{v_i}=\frac{E\&PlusMinus;N_{v_i}}{\left|\right|E\&PlusMinus;N_{v_i}\left|\right|}.$

如果成立时，符号取负，反之符号取正。

由此，得到Householder矩阵为：

$Q_{v_i}=I-2W_{v_i}{W}_{v_i}^{\&prime;}---\left(16\right)$

该矩阵是正交的，并且第一列为或(取决于前面的符号选择)，其余两列决定了切空间上一对标准正交基。

设两个向量和由于是矩阵的特征值为0的特征向量，则有：

$Q_{v_i}^{\&prime;}{M}_{v_i}{Q}_{v_i}=\left(\begin{array}{ccc}0& 0& 0\\ 0& {\stackrel{\&OverBar;}{m}}_{v_i}^{11}& {\stackrel{\&OverBar;}{m}}_{v_i}^{12}\\ 0& {\stackrel{\&OverBar;}{m}}_{v_i}^{21}& {\stackrel{\&OverBar;}{m}}_{v_i}^{22}\end{array}\right)---\left(17\right)$

其中， ${\stackrel{\&OverBar;}{m}}_{v_i}^{12}={\stackrel{\&OverBar;}{m}}_{v_i}^{21}.$

(4)Givens变换将矩阵对角化，

我们采用Givens变换，将矩阵对角化，获得角θ和该矩阵的特征值。可以求出cosθ和sinθ，然后可以求出主曲率的方向分别为：

$T_1=cos\mathrm{\&theta;}\stackrel{\&OverBar;}{T_1}-sin\mathrm{\&theta;}\stackrel{\&OverBar;}{T_2}---\left(18\right)$

$T_2=sin\mathrm{\&theta;}\stackrel{\&OverBar;}{T_1}+cos\mathrm{\&theta;}\stackrel{\&OverBar;}{T_2}---\left(19\right)$

向量是矩阵的特征向量，即为曲面上点v_i处的主方向，则得到主曲率的计算公式(9)和(10)。

3.顶点主曲率计算步骤

(1)计算法曲率

对于三角网络上的顶点p_i，顶点p_j是顶点p_i的1-ring邻域上的点，沿着切方向，最终采用此公式估计顶点p_i在方向T_ij上的法曲率：

为每个mesh三角网格f_k都定义单位长度的法向量N_fk。曲面是有方向的，对某点p_i的法向量N_pi，Taubin面积加权和的定义是：

$N_{p_i}=\frac{{\mathrm{\&Sigma;}}_{f_k\&Element;F^{\&prime;}}\left|f_k\right|N_{f_k}}{\left|\right|{\mathrm{\&Sigma;}}_{f_k\&Element;F^{\&prime;}}\left|f_k\right|N_{f_k}\left|\right|}$

(2)估计矩阵

估计矩阵定义为顶点p_i的1-ring邻点面积和沿T_ij方向曲率的乘积的加权和，具体计算积分矩阵方法： $M_{v_i}={\mathrm{\&Sigma;}}_{p_j\&Element;v^{\&prime;}}{w}_{ij}{k}_{ij}{T}_{ij}{T}_{ij}^{\&prime;}$

对于顶点p_i的每个1-ring邻点，定义单位向量T_ij是向量在切平面(平面法线为)上的投影，即：沿方向T_ij曲率估计方法为： $k_{\mathrm{ij}}=\frac{2<N_{p_i},p_j-p_i>}{{\left|\right|p_j-p_i\left|\right|}^2};$ 权值w_ij为面积加权。

(3)求曲率主方向和主曲率

令 $E=\left[\begin{array}{c}1\\ 0\\ 0\end{array}\right],W_{p_i}=\frac{E\&PlusMinus;N_{p_i}}{\left|\right|E\&PlusMinus;N_{p_i}\left|\right|}.$ 若 $\left|\right|E-N_{p_i}\left|\right|>\left|\right|E+N_{p_i}\left|\right|$ 成立时，取负号，反之取正。由此得Householder矩阵：该矩阵正交，且第一列为其余两列决定了切空间上一对标准正交基。

设两个向量和由于是矩阵的特征值为0的特征向量，则： $Q_{p_i}^{\&prime;}{M}_{p_i}{Q}_{p_i}=\left[\begin{array}{ccc}0& 0& 0\\ 0& {\stackrel{\&OverBar;}{m}}_{p_i}^{11}& {\stackrel{\&OverBar;}{m}}_{p_i}^{12}\\ 0& {\stackrel{\&OverBar;}{m}}_{p_i}^{21}& {\stackrel{\&OverBar;}{m}}_{p_i}^{22}\end{array}\right]$ 其中， ${\stackrel{\&OverBar;}{m}}_{p_i}^{12}={\stackrel{\&OverBar;}{m}}_{p_i}^{21}.$

随后用Givens变换将矩阵对角化，获得角θ和该矩阵的特征值求出cosθ和sin，随后求出主曲率的方向：向量是矩阵的特征向量，即为曲面上点p_i处的主方向，则得到主曲率的计算公式：

$\left\{\begin{array}{c}{k}_p^1=3m_p^{11}-m_p^{22}\\ k_p^2=3m_p^{22}-m_p^{11}\end{array}\right.$

二、根据主曲率对牙齿网格模型进行初步分割；

如图2所示，采用上述方式，牙齿网格模型的每个顶点处都计算得到两个主曲率这两个曲率分别为两个正交方向(T1和T2)的曲率。对于峰点处的顶点，这两个正交方向的曲率绝对值都比较大，如图2a所示，对于谷点也是一样。但对于鞍形点，存在一个方向的主曲率绝对值比较大，一个方向的主曲率绝对值可能比较小，通常可能快接近0，如图2b所示。对于鞍形点，采用平均曲率就无法刻画该处的真实弯曲程度。因此，我们采用主曲率绝对值的最小值作为分割网格的依据，根据主曲率最小值和曲率预先设定的阈值对网格顶点进行分割。

在具体分割中，针对于网格模型任意顶点p，根据主曲率阈值t，将该顶点初步分割为平坦点和弯曲点，具体步骤如下：

1.针对牙齿网格模型任意顶点P，计算P的主曲率绝对值；

2.计算P的主曲率绝对值的最小值

3.若则P点为平坦点，加入平坦点集Ps，否则P点为弯曲点，加入弯曲点集Pc，其中t为设定的阈值。

采用这种方式，将整个网格模型初步分割为平坦点集Ps和弯曲点集Pc。一般来说，在牙齿和牙龈交界处的顶点曲率比较大，这些顶点被分割为弯曲点集Pc，但在牙齿咬合面上也存在曲率比较大的弯曲点，故需要进一步处理去除这些顶点的干扰。

三、通过区域生长法初步分割出牙齿齿面；

为了能够自动分割出牙齿，我们需要在牙齿网格模型上采用人工交互的方式，拾取牙齿齿面上任意平坦点，然后采用区域生长(RegionGrowing)自动分割出指定牙齿的网格模型顶点。区域生长是指将成组的像素或区域发展成更大区域的过程。从种子点的集合开始，从这些点的区域增长是通过将与每个种子点有相似属性像强度、灰度级、纹理颜色等的相邻像素合并到此区域。在本发明中，采用的是牙齿网格模型上顶点是属于平坦点集Ps还是弯曲点集Pc进行区域生长自动分割的。具体牙齿预分割方法如下：

1.人工选取牙齿齿面上任意顶点N，N属于平坦点集Ps，然后将放入待处理顶点集Pf；

2.从待处理顶点集Pf取出一个顶点M，将M放入牙齿顶点集Pt；

3.获得顶点M的1-ring邻域顶点Mn；

4.若Mn属于平坦点集Ps，则将Mn放入待处理顶点集Pf；

5.若待处理顶点集Pf不为空，则跳到步骤2，否则退出，得到牙齿顶点集Pt。

通过上述区域生长迭代，我们得到了指定牙齿的预分割顶点集合。但是，由于牙齿咬合面存在曲率比较大的顶点，因此还需要后续方法将这些曲率比较大的顶点加入到牙齿顶点集合，形成完整的牙齿模型。此时，牙齿牙龈模型被分割三种类型区域，也就是牙齿顶点集Pt、平坦点集Ps和弯曲点集Pc。

四、根据联通性找回牙齿咬合面的孤立区域，形成完整牙齿齿面。

牙齿咬合面处的顶点曲率比较大，故无法采用曲率和区域增长法自动将这些顶点分割为牙齿顶点。也就是，在牙齿顶点集Pt的内部可能存在弯曲点集Pc的子集。如图3a所示，牙齿顶点集Pt内部存在弯曲点区域R1、R2和R3，这些区域通常为牙齿咬合面顶点。同时，在剩余的平坦区域也存在孤立的弯曲点区域R4、R5、R6和R7。

但我们发现，此时牙齿咬合面顶点R1、R2和R3被初步分割后的牙齿顶点集Pt区域包围。通过分析顶点网格的这种联通包围性，就可以将这些牙齿咬合面顶点R1、R2和R3补齐到牙齿顶点集Pt中，从而完成牙齿的准确分割，即得到新的牙齿顶点集合，如图3b所示。

Claims (3)

1.一种基于主曲率的牙齿网格模型自动分割方法，其特征在于，包括如下步骤：

(1)计算牙齿网格模型在每个顶点处的主曲率

(2)根据主曲率对牙齿网格模型进行初步分割；

(3)通过区域生长法初步分割出牙齿齿面；

(4)根据联通包围性找回牙齿咬合面的孤立区域，形成完整牙齿齿面。

2.根据权利要求1所述的基于主曲率的牙齿网格模型自动分割方法，其特征在于，所述步骤

(2)中根据主曲率对牙齿网格进行初步分割的步骤为：

(21)针对牙齿网格模型任意顶点P，计算P的主曲率绝对值；

(22)计算P的主曲率绝对值的最小值

(23)若则P点为平坦点，加入平坦点集Ps，否则P点为弯曲点，加入弯曲点集Pc，其中t为设定的阈值。

3.根据权利要求1所述的基于主曲率的牙齿网格模型自动分割方法，其特征在于，所述步骤

(3)中通过区域生长法初步分割出牙齿齿面的步骤为：

(31)人工选取牙齿齿面上任意顶点N，N属于平坦点集Ps，然后将放入待处理顶点集Pf；

(32)从待处理顶点集Pf取出一个顶点M，将M放入牙齿顶点集Pt；

(33)获得顶点M的1-ring邻域顶点Mn；

(34)若Mn属于平坦点集Ps，则将Mn放入待处理顶点集Pf；

(35)若待处理顶点集Pf不为空，则跳到步骤(32)，否则退出，得到牙齿顶点集Pt。