发明内容
针对现有技术的不足,本发明是基于三维模型检索技术的牙齿建模方法。先收集大量的牙齿模型对于进行牙齿模型库,接着将不同类型、不同特征的牙冠数据进行三维模型检索,查找到已有的相似牙冠模型后,将相似牙冠的牙齿根部和待建牙冠通过三点平移变换法进行模型配准,通过AABB包围盒算法对待建模牙冠与相似牙冠模型进行裁剪,使用平面裁剪三角网格模型时,可能会出现不被希望出现的四边形或者狭小三角形,这时需要对其进行调整,如将四边形继续分为两个或者多个三角形,对于狭小三角形则采用合并处理,之后需要建模的牙冠数据与已有的牙根数据合并成为一个完整的牙齿模型。
在研究建立牙齿模型方法的同时,本文还提出了一种基于区域分割的牙齿类型识别方法。通过分析不同类型恒牙牙冠特征的差异,利用高斯曲率和平均曲率两种曲率描绘的曲面特征对牙冠表面进行区域分割。根据不同类型恒牙牙冠特征的差异和分割区域所描述的牙冠曲面特征对牙冠进行类型识别。该方法相较于现有的技术,可以将牙冠模型被分割成形状相对简单,特征较为明显的多个区域。这些区域共同组成了整个牙冠曲面,代表了原始牙冠曲面的形状特征。
本发明所采用的技术方案为:
一种基于三维模型检索的牙齿建模方法,其特征在于,包括如下步骤:
(1)建立牙齿模型库,通过对牙冠进行区域分割而提取四个特征值;
(2)读入待建模的牙冠数据模型,将待建模的牙冠数据模型进行区域分割,根据分割的区域计算该牙冠的描述子并且通过类型识别判断其牙齿类别;
(3)根据牙冠的描述子及牙齿类别在模型库中检索出相似度最高的牙齿模型;
(4)使用三点平移变换法将待建模的牙冠数据模型与检索的牙齿模型配准;
(5)对两个模型裁剪与拼接,得到完整的牙齿模型。
步骤(1)在建立牙齿模型库时,需要对每个恒牙牙冠进行区域分割,从而便于不同区域的特征值提取;对牙冠进行区域分割时,通过分析不同类型恒牙牙冠特征的差异,利用高斯曲率和平均曲率两种曲率描绘的曲面特征对牙冠表面进行区域分割。
所述的区域分割是通过根据曲面的微分几何特性,曲面的平均曲率H和高斯曲率K反映了曲面的形状特征;Kf>0且Hf>0代表峰的区域称作I类区域,描述的是对应牙冠牙合面的牙尖特征,Kf<0且Hf>0代表嵴的区域称作II类区域,描述的是对应牙合面上的各种嵴,Hf<0代表谷的区域称作III类区域,描述的是牙冠牙合面上的沟、窝等各种凹陷。
所述步骤(1)中的四个特征值,分别是区域类型、区域总曲率、区域相对面积和邻接区域的相对边长。
根据权利要求4所述的一种基于三维模型检索的牙齿建模方法,其特征在于,所述区域类型分为三类,分别反映牙冠曲面上的牙尖特征区域(凸特征区域)、切嵴等嵴状特征区域(双曲点特征区域)和沟、窝等谷状特征区域(凹特征区域)。
所述的区域总曲率对分割区域内的高斯曲率进行面积分,以得到分割区域的总曲率为:
对于KA由微分几何知识可知,其具有缩放不变性。
所述区域相对面积ΔSt作为特征值评价曲面相似度。区域相对面积指该区域面积占曲面总面积的百分比:
其中St为该区域面积,S为曲面总面积。
所述邻接区域的相对边长,记录了区域与相邻的三种类型区域公共边的长度与区域总边长的比值,若与某一种区域不相邻,则其相对边长为0。
所述步骤(2)中所述的描述子是由步骤(1)中的四个特征值组成的一个6维的广义向量。
所述步骤(3)的相似度比较是根据四个特征值依次比较的。
所述相似度比较中,对于区域类型进行比较,如果两个相比较区域类型相同,则认为两个区域相似,其相似度为1;反之则认为两区域间不相似,相似度为0,表示为:
所述相似度比较中,对于区域总曲率进行比较,则是令两个区域总曲率为KA1和KA2,则区域间总曲率相似度表示为:
所述相似度比较中,对于区域相对面积进行比较,是令两个区域的相对面积为ΔS1和ΔS2,则区域间相对面积相似度表示为:
所述相似度比较中,对于对邻接区域相对边长进行比较,由于邻接区域相对边长记录了区域与相邻的三种类型区域公共边的长度与区域总边长的比值,是一个三维向量;在做相似度比较时,使用与上述相似度类似的计算方法,并求平均值:
步骤(3)中的相似度比较时,两个区域间相似度表示为:
其中α1、α2和α3代表相应权值。
所述步骤(4)使用三点平移变换法前需要选取种子点和目标点进行匹配。
所述的选取种子点,需在待建模牙冠曲面上选取不共线的三点S1、S2和S3作为种子点与目标牙齿模型建立联系;以种子点si的高斯曲率和平均曲率为匹配特征搜索目标牙齿模型,找出所有满足条件的对应点,其中约束条件为:
pj为目标牙齿模型中一点;δk和δh分别为高斯曲率误差及平均曲率误差。
在进行匹配时,种子点与对应点之间难免会出现一对多的对应联系,从而降低匹配成功率,因此需要建立一定的约束条件排除不正确的对应联系,根据种子点的法矢及距离,建立三角约束条件,进一步筛选对应点,获得目标点;三角约束条件如下:
代表种子点S1、S2和S3法矢间的夹角;代表种子点S1、S2和S3间的距离;δα和δd分别为法矢间夹角的误差和种子点距离的误差。通过进一步的筛选,最终得到最优目标点m1、m2和m3。
所述步骤(4)中的三点平移变换法步骤如下:
(1)由种子点S1、S2和S3构建局部坐标系,令局部坐标系为
CoorS=(coorx(S),coory(S),coorz(S))
(2)以S1为坐标原点,令S1到S3的方向为x轴方向,则
(3)再以S1到S2的方向向量与coorx(S)的叉乘作为y轴方向,则
(4)最后利用coorx(S)与coory(S)的叉乘确定z轴方向
coorz(S)=coorx(S)×coory(S)
(5)类似的方法可求出目标点m1、m2和m3的局部坐标系Coorm,通过公式推导可得坐标系Coors经旋转平移变换到坐标系Coorm的旋转矩阵R和平移向量T:
R=Coorm(Coorx)T
T=(m1+m2+m3)/3-R(S1+S2+S3)/3
所述步骤(5)对两个模型裁剪与拼接,在裁剪时会用到AABB包围盒算法。
所述步骤(5)中的裁剪是使用平面裁剪三角网格,可能会出现非理想状态,即四边形和狭小三角形;处理出现的四边形,则是将其分为两个或者更多的小三角形;处理出现的狭小三角形,一般对其进行合并处理。
为了证明本发明的可行性,则进行了试验,通过试验结果表明,该方法不仅简单快捷,更能得到良好的建模效果。
具体实施方式
如图1所示,本发明的具体建模步骤如下:
步骤(1)建立牙齿模型库。
在建立牙齿模型库时,需要对每个恒牙牙冠进行区域分割,从而便于不同区域的特征值提取,这样的特征值提取更为精确。对牙冠进行区域分割时,通过分析不同类型恒牙牙冠特征的差异,利用高斯曲率和平均曲率两种曲率描绘的曲面特征对牙冠表面进行区域分割。根据不同类型恒牙牙冠特征的差异和分割区域所描述的牙冠曲面特征对牙冠进行类型识别。
步骤(1)中对牙冠的区域分割则是根据牙冠的凹凸性及形态特征,可将特征划分为三种:峰、嵴、谷。其中峰对应牙冠牙合面上牙尖的特征;嵴与牙合面上的各种嵴的概念相同;谷则代表牙冠牙合面上的沟、窝等各种凹陷。
根据曲面的微分几何特性,曲面的平均曲率和高斯曲率反映了曲面的形状特征。平均曲率H根据其取值的正负反映了曲面上点局部的凹凸性,H>0的点局部表现为凸,H<0的点局部表现为凹。高斯曲率K则表现了曲面上点的形状信息,当K>0时,曲面上点表现为椭圆点,K=0表现为抛物点,K<0则为双曲点。基于上述微分几何知识,结合恒牙牙冠特征,则可将牙冠表面上的区域划分为一下三个区域,。Kf>0且Hf>0代表峰的区域称作I类区域,描述的是对应牙冠牙合面的牙尖特征,Kf<0且Hf>0代表嵴的区域称作II类区域,描述的是对应牙合面上的各种嵴,Hf<0代表谷的区域称作III类区域,描述的是牙冠牙合面上的沟、窝等各种凹陷。
步骤(1)中的特征值提取则是用高斯曲率和平均曲率来描述曲面。由微分几何特性可知,曲面上一点主方向上两个主曲率如果用k1和k2表示,则高斯曲率K为k1和k2的乘积,即:
K=k1×k2
平均曲率H为k1和k2和的一半,即:
H=(k1+k2)/2
本发明使用的三维牙齿模型数据是三角网格模型,而三角网格模型是一种分段连续性模型,不存在连续曲率。故无法使用上述两式计算三角网格模型上点的高斯曲率和平均曲率。对于三角网格模型,通常只计算三角形面片顶点处的法矢和曲率。
图2表示了顶点vi周围的邻域情况,为顶点vi的法矢,为由vi、vj和vj+1三点组成的三角面片fk的法矢,αk为三角面片fk在顶点vi处的顶角,具体算法如下:
第一步需要计算,可以如下定义:
式中eij和ei,j+1分别表示由顶点vi指向顶点vj和vj+1的两条边矢量。计算顶点vi的法矢时,一般常用三角面片fk,k∈planes(vi)的面积进行加权平均。其中planes(vi)表示点vi所有邻接三角面片的集合。然而相同面积的两个三角形其形状可能相差很大,所以使用顶角与面积共同加权去对进行加权平均:
式中Ak为三角面片fk的面积,αk为三角面片fk在顶点vi处的顶角。从而,利用计算出的公式和古典微分几何中可得离散法曲率:
式中vj代表vi周围相邻顶点。
同时,由古典微分几何中的平均曲率公式:
最终可得离散平均曲率公式:
根据经典微分几何中的Gauss-Bonnet定理可直接计算离散高斯曲率:
式中,αk为vivj与vivj+1的夹角,对Gauss-Bonnet定理计算的离散高斯曲率中积分进行离散,可以得到离散高斯曲率公式:
步骤(1)的牙冠描述子一共有四个,分别是区域类型、区域总曲率、区域相对面积和邻接区域的相对边长。
其中区域类型则是将步骤(1)中的区域分割分为三类,分别反映牙冠曲面上的牙尖特征区域(凸特征区域)、切嵴等嵴状特征区域(双曲点特征区域)和沟、窝等谷状特征区域(凹特征区域)。将分割后区域所对应的类型标志作为特征值提取、保存,用于检索时相似性的比较。
高斯曲率作为描述曲面形状特征的重要变量,具有平移和旋转不变性,但不具有缩放不变性。然而在进行三维检索评价曲面相似度时,希望评价因素具有缩放不变性。故对分割区域内的高斯曲率进行面积分,以得到分割区域的总曲率:
由微分几何知识可知,式中的KA具有缩放不变性。
由于三角网格模型所表示的曲面不连续,故KA可近似为:
其中Kf为分割区域中三角面片上的高斯曲率,Af为分割区域中三角面片的面积。
与区域总曲率Kf类似。由于区域面积不具有缩放不变性,所以使用区域相对面积ΔSt作为特征值评价曲面相似度。区域相对面积指该区域面积占曲面总面积的百分比:
其中St为该区域面积,S为曲面总面积。
邻接区域的相对边长记录了各分割区域间的连接关系,体现了各分割区域之间的拓扑特征。之所以使用相对边长,是为了消除缩放对边长的影响,其定义与相对面积类似。该特征值记录了区域与相邻的三种类型区域公共边的长度与区域总边长的比值,若与某一种区域不相邻,则其相对边长为0。
综上所述,本发明将以上四种特征综合起来形成一个6维的广义向量作为描述牙冠曲面特征的描述子。然后将各牙齿模型及其牙冠的描述子按照牙齿类型分类保存,建立检索所需的牙齿模型库,模型库中包含牙齿模型的STL文件及其牙冠描述子的二进制文件。
步骤(2)读入待建模的牙冠数据模型,将待建模的牙冠数据模型进行区域分割,根据分割的区域计算该牙冠的描述子并且通过类型识别判断其牙齿类别。
以编号2-1的牙冠模型为例,按照步骤(1)中的描述子读入待建模的牙冠数据模型,如图3所示。
步骤(3)根据牙冠的描述子及牙齿类别在模型库中检索出相似度最高的牙齿模型。
检索相似度的步骤如下:先对区域类型进行比较。如果两个相比较区域类型相同,则认为两个区域相似,其相似度为1。反之则认为两区域间不相似,相似度为0,表示为
再对区域总曲率进行比较。令两个区域总曲率为KA1和KA2,则区域间总曲率相似度表示为
接着对区域相对面积进行比较。令两个区域的相对面积为ΔS1和ΔS2,则区域间相对面积相似度表示为
最后对邻接区域相对边长进行比较。由于邻接区域相对边长记录了区域与相邻的三种类型区域公共边的长度与区域总边长的比值,是一个三维向量。在做相似度比较时,使用与上述相似度类似的计算方法,并求平均值。
根据上述四中特征值相似度的比较,两个区域间相似度表示为:
其中α1、α2和α3代表相应权值,经试验三个权值分别取0.25、0.25和0.5相似性比较效果最好。
牙冠的相似度由各个区域相似度的总和构成。两个牙冠之间区域相似度比较则可以看做由两组区域节点组成的完全二分图的最优匹配的过程,如图4所示,设S1和S2分别代表两个牙冠,ui为S1的区域节点,vi为S2的区域节点。
将ui与vi依次配对,计算出每对区域节点的相似度simAij,可得到S1与S2之间的相似度矩阵。如果S1与S2之间区域节点数不同,则通过补0的方法使矩阵成为方阵,即:
为了得到牙冠间总体配对最优方案,采用图论中计算赋权二分图最优匹配的Kuhn-Munkres算法计算最优匹配方案。最终可得:
其中simAm(j)j表示S1第m(j)个区域与S2中第j个区域的相似度;m(j)表示M中第j列最优匹配的行数;ΔSm(j)与ΔSj分别表示S1第m(j)个区域与S2中第j个区域的相对面积。
根据该牙冠模型的描述子,与模型库中数据进行比较后,图4给出了对编号2-1磨牙牙冠在模型库中进行检索的前7个检索结果,由图可知,编号u01牙冠模型的相似度相比较其他牙冠模型最高。因此将编号u01作为目标模型。
步骤(4)使用三点平移变换法将待建模的牙冠数据模型与检索的牙齿模型配准。
如图6所示,图中蓝色模型为待建模牙冠2-1,灰色模型为拥有相似牙冠的牙齿模型u01。从图中可看出,两个模型之间存在位置上的偏差。
曲面匹配算法可以很好地解决这一问题。本发明所使用的牙冠曲面匹配算法简捷、快速。算法以曲率作为匹配特征,在待建模牙冠与目标牙齿模型之间建立满足角度和距离约束的对应关系。利用三点平移变换法生成旋转矩阵及平移向量,实现两者间的匹配。
进行匹配算法前,需要在待建模牙冠曲面上选取不共线的三点S1、S2和S3作为种子点与目标牙齿模型建立联系。以种子点Si的高斯曲率和平均曲率为匹配特征搜索目标牙齿模型,找出所有满足条件的对应点,其中约束条件为:
其中pj为目标牙齿模型中一点,δk和δh分别为高斯曲率误差及平均曲率误差。
由于模型上相似曲率点的存在,种子点与对应点之间难免会出现一对多的对应联系,从而降低匹配成功率。因此需要建立一定的约束条件排除不正确的对应联系。本文根据种子点的法矢及距离,建立三角约束条件,进一步筛选对应点,获得目标点。三角约束条件如下:
代表种子点S1、S2和S3法矢间的夹角;代表种子点S1、S2和S3间的距离;δα和δd分别为法矢间夹角的误差和种子点距离的误差。通过进一步的筛选,最终得到最优目标点m1、m2和m3。
步骤(4)中的三点平移变换法是为了实现待建模牙冠与目标牙齿模型的匹配。其步骤如下:
(1)由种子点S1、S2和S3构建局部坐标系,令局部坐标系为:
CoorS=(coorx(S),coory(S),coorz(S))
以S1为坐标原点,令S1到S3的方向为x轴方向,则
(2)再以S1到S2的方向向量与coorx(S)的叉乘作为y轴方向,则
(3)最后利用coorx(S)与coory(S)的叉乘确定z轴方向:
coorz(S)=coorx(S)×coory(S)
类似的方法求出目标点m1、m2和m3的局部坐标系Coorm,通过公式推导可得坐标系Coors经旋转平移变换到坐标系Coorm的旋转矩阵R和平移向量T:
R=Coorm(Coorx)T
T=(m1+m2+m3)/3-R(S1+S2+S3)/3
最终将旋转矩阵R和平移向量T作为变换矩阵叠加到原始模型数据上,图6中所示的两个模型经匹配后,结果如图7所示。
步骤(5)对两个模型裁剪与拼接,得到完整的牙齿模型。
经过匹配算法处理的待建模牙冠模型在空间位置上已经大致与相似牙齿模型的牙冠位置相同,如图7所示。这时只需将相似牙齿模型的牙根部分按待建模牙冠模型的位置与大小进行裁剪,之后拼接到待建模牙冠模型上即可完成建模。
裁剪时需要确定裁剪平面的位置。将牙冠模型看做由空间中一个与其紧密贴合的长方体包裹,根据长方体的六个面确定裁剪平面的位置,此方法也称为包围盒算法。本发明中为了减少运算、提高速度,建立了一个模型自身的坐标系,如图8所示。
该坐标系由三个向量和一个点坐标组成,其中三个向量两两垂直,作为坐标系的x轴、y轴和z轴;点坐标作为坐标系的原点。通过人机交互的方式对建立的坐标系进行调整,使该坐标系的z轴方向垂直穿过牙冠牙合面;x轴方向与牙冠的唇舌方向一致,并指向唇面方向;y轴方向与牙冠的近远中方向一致,并指向远中方向。
构造旋转矩阵R和平移矩阵T使牙冠数据模型的自身坐标系通过变换矩阵作用后与世界坐标系重合。将该旋转矩阵R和平移矩阵T同时作用于待建模牙冠模型与相似牙齿模型,这时使用AABB包围盒算法,利用平行于坐标轴的最小长方体包裹待建模牙冠模型,该长方体z轴负半轴方向上平面的位置就是所需裁剪平面的位置,如图9所示。
在使用平面裁剪三角网格模型时,理想状态如图9中三角形ΔAHI。但也可能出现图9所示的四边形ΔADJI以及狭小的三角形ΔEML。
对于上述裁剪后的不理想状态,本文采用不同的方法对其进行处理。针对产生的四边形ΔADJI对其进行图10所示处理。
图10中(a)方法是利用四边形的对角线将四边形划分为两个三角形;(b)中的方法是将被平面裁剪产生边的中点与其对边的两个端点连接,使原四边形分割成三个三角形;(c)方法则较为繁琐,它将裁剪产生边的两个端点与其对边的中点相连,同时将其对边的中点与相邻三角形顶点相连。三种方法各有优劣。
然而无论是图10中所示的哪种方法,在实际处理四边形时,都不可避免的会出现裁剪后的第三种情况:产生狭小三角形。例如图9中四边形ENOM,无论使用哪种方法都会有新的狭小三角形出现。
处理狭小三角形一般对其进行合并处理,为了避免合并狭小三角形时产生的繁杂运算,本文在进行裁剪之前,对裁剪平面附近的三角网格顶点坐标做一定调整。使距离裁剪平面较近的三角网格的顶点移动到平面上,例如对图9中的三角网格处理后,可得图13。
相较于图9的裁剪结果,图13的结果不仅解决了狭小三角形的问题,同时也使裁剪过程中产生的新点数明显减少,简化了裁剪过程。
裁剪后,在待建模牙冠模型与裁剪得到的牙根模型之间建立三角网格将两者拼接完成建模。在将两者使用三角网格进行拼接的过程中,往往会出现待建模牙冠模型与相似牙齿模型尺寸大小不同的情况。这时需要再次利用确定裁剪平面时生成的AABB包围盒,根据包围盒中分别与x轴和y轴平行的四个面确定待建模牙冠和裁剪得到牙根的长宽利用两个模型的长宽数据对裁剪得到的牙根模型进行缩放,最后将待建模的牙冠模型与裁剪好的牙根模型间以三角面片相连完成拼接,如图14所示。
经过建模实验发现,基于三维检索的牙齿建模方法在对待建模牙冠数据进行建模时,当模型库中存在与待建模牙冠数据相似度较高的牙齿模型时,可以通过本发明所述的牙齿建模方法得到所需要的牙齿模型,经过计算得出,当相似度大于0.8时,该建模方法能得到较好的建模结果。然而当模型库中的牙齿模型与待建模牙冠数据相似度较低时,建模结果往往不太理想。所以在建立模型库时应该尽可能多的收集牙齿模型,从而使建模结果更准确。