一种基于路径规划的三角网格牙齿分割方法
技术领域
本发明属于计算机图形学领域,涉及一种从三角网格牙颌模型上自动分割所有牙齿的基于路径规划的三角网格牙齿分割方法。
背景技术
随着三维测量、计算机图形学技术的不断发展,计算机辅助诊断修复已广泛的应用于口腔正畸当中。在正畸系统中,为了模拟牙齿牙龈运动、制定矫治方案,需要首先从数字化的三维牙颌模型上将每颗牙齿分割出来。因此,分割的精度和完整性直接关系到最终的矫治效果和设计成本。
目前三维牙齿分割方法大部分是利用图形图像处理技术提取牙颌模型的牙齿形态特征来分离单颗独立的牙齿,然而由于三维扫描精度的有限性,三角网格牙颌模型上牙齿和牙龈的边界线以及相邻牙齿之间的边界线通常融合在一起,再加上牙齿形态各异,因此传统的分割方法很难获得良好的分割结果,尤其当患者具有畸形牙时,这一问题尤为突出。Wu等人(K.Wu,L.Chen,J.Li,et al.Tooth segmentation on dental meshes usingmorphologic skeleton,Computers&Graphics.2014,38:199–211)首先基于骨架线法提取牙齿和牙龈的边界线,然后通过匹配边界线的角点来获得相邻牙齿之间的分界线,但是在匹配过程中容易受到骨架线的毛刺以及中心牙窝线的干扰。Zou等人(B.Zou,S.Liu,S.Liao,et al.Interactive tooth partition of dental mesh base on tooth-targetharmonic field,Computers in Biology and Medicine,2015,56:132–144)提出一种基于调和场的牙齿分割方法,但需要大量的用户交互来设置牙齿分界的约束点。
发明内容
为克服上述现有技术的不足,提高牙齿分割的效率和精度,提出一种新颖的基于路径规划的三角网格牙齿分割方法。
本发明的技术方案是:一种基于路径规划的三角网格牙齿分割方法,其步骤如下:
1)确定咬合平面,通过在三维牙颌模型上交互式选择牙齿解剖特征点,所述特征点包括:第一磨牙的近中舌尖点、中切牙的切缘点,并利用主成分分析法拟合所述特征点以确定咬合平面,并进行坐标变换以使三维牙颌模型的XY平面对齐到咬合平面,Y轴穿过两中切牙的牙缝;
2)提取特征区域F,其中,第一步,对三维牙颌模型进行离散曲率估算,并对得到的曲率值进行基于直方图均衡化的拉伸变换以突出凹凸特征,第二步,对曲率提取的顶点集合进行连通性过滤以剔除噪音,以获取牙齿和牙齿以及牙齿和牙龈分界的特征区域F;
3)拟合牙弓线,将2)中获得特征区域投影到1)中的咬合平面上以构建特征二值图像,对二值图像进行形态学操作以获取牙弓曲线,并利用牙弓曲线将特征区域进行分区;
4)搜索牙龈线,将特征区域F映射为无向连通图,基于二次优路径规划算法搜索牙龈线:第一次搜索时,无向连通图中边的权值定义为两顶点之间的欧式距离;第二次搜索在第一次搜索结果的路径邻域附近的特征区域F上再次搜索,并将无向连通图中边的权值定义为两顶点欧式距离和曲率的乘积,以实现牙龈线的精确定位;
5)搜索牙缝线,第一步,根据牙龈线和牙弓曲线的特征信息找到牙龈线在牙缝处弯曲的位置,从而将整个三维牙颌模型的牙龈线分割为每颗牙齿的牙龈线,然后搜索每颗牙齿牙龈线两侧的最优路径,获得牙缝线;
6)分割牙齿,将每颗牙齿的牙龈线和牙缝线组合,以构成每颗牙齿封闭的分割线,然后利用区域生长算法将每颗牙齿从三维牙颌模型上分割出来。
步骤2)中包括以下步骤:
(1)利用局部三次曲面拟合法计算每个网格顶点的最大主曲率值,并将所有网格顶点的曲率值都归一化到[0,1],然后将归一化后的曲率值集合{ki,i=1,2,…,N,N为网格顶点总数}平均划分为L个区间:计算曲率在每个区间的累计分布概率密度:
其中,表示曲率值满足的网格顶点个数;
然后,每个网格顶点的曲率ki按照如下公式进行变换:
(2)利用曲率值提取三维牙颌模型的初始特征区域F0:
F0={pi|k(pi)≥H,(i=1,2,...,N)}
其中,k(pi)为网格顶点pi变换后曲率值,H为曲率阈值;
将初始特征区域F0根据网格顶点之间的连接关系映射为无向连通图G0,G0中边的权值edge(vi,vj)按下式定义:
其中,E为三维牙颌模型所有三角边的集合,计算无向连通图G0的连通分量,并统计每个连通分量中的网格顶点个数,保留具有最多网格顶点个数的连通分量,即可删除噪音分量,得到牙齿和牙龈以及相邻牙齿之间分界的特征区域F。
步骤3)中包括以下步骤:
(1)将步骤2)提取的特征区域F投影到咬合平面上,并构建特征二值图像;
(2)对特征二值图像进行孔洞填充以及细化操作获取牙弓骨架线;
(3)将牙弓骨架线利用最小二乘法线拟合为四次多项式曲线:y(x)=ax4+bx3+cx2+dx+e,从而得到牙弓曲线arch;
(4)搜索牙弓曲线与特征二值图像外轮廓的交点,并映射回到三维特征区域当中,得到两个特征点m、n,并以这两个特征点为界,定义内侧特征区域为舌侧区域,外侧特征区域为颊侧区域。
步骤4)中包括以下步骤:
(1)先将特征区域F根据网格顶点之间的连接关系映射为无向连通图G1,G1中边的权值edge(vi,vj)定义为三维牙颌模型两网格顶点vi,vj之间的欧式距离:
利用Dijkstra算法搜索图G1从点m到点n的最优路径,即为舌侧牙龈线;然后将搜索到的舌侧牙龈线邻域附近的边权值赋为无穷大,再次搜索从特征点m到特征点n的最优路径,即为颊侧牙龈线,E为三维牙颌模型所有三角边的集合;
(2)找到舌侧和颊侧牙龈线邻域附近的特征区域F’,并根据网格顶点连接关系映射为无向连通图G2,G2中边的权值定义为两顶点欧式距离和曲率的乘积:
其中,k(vi)和k(vj)分别表示网格顶点vi和vj的曲率值,||vi-vj||表示两网格顶点之间的欧式距离,然后重新利用Dijkstra算法搜索舌侧牙龈线和颊侧牙龈线。
步骤5)包括以下步骤:
(1)将步骤4)中得到的舌侧牙龈线和颊侧牙龈线利用能量法分别拟合为B样条曲线:flingual和fbuccul,并计算B样条曲线的曲率:ρlingual和ρbuccal以及B样条曲线的二阶导数:和然后结合牙弓曲线的特征信息来确定舌侧牙龈线和颊侧牙龈线在牙缝弯曲处的位置:
其中,T为牙龈线的曲率阈值,archi为牙弓曲线上距离待测牙龈线的最近点,为牙弓曲线上点archi的法向量,将牙龈线在牙缝处的拐点删除,并将剩余的舌侧牙龈线和颊侧牙龈线上的点按顺序进行组合,以获得每颗牙齿的牙龈分割线;
(2)利用步骤4)所构建的无向连通图G2,基于Dijkstra算法分别搜索每颗牙齿的舌侧牙龈线起点到颊侧牙龈线起点的最优路径,以及舌侧牙龈线终点到颊侧牙龈线终点的最优路径,由此获得每颗牙齿两侧的牙缝线。
步骤6)中将每颗牙齿的牙龈线和牙缝线组合构成分割线后,计算该分割线的质心,并将三维牙颌模型上与该质心的x、y坐标最相近的点作为种子点,不断搜索种子点的邻域,直到达到当前牙齿的分割线为止。
本发明具有较少的人工干预和参数调整,而且由于该方法在搜索牙齿分割线时,将牙龈线和牙缝线分开搜索,能够有效避免牙冠网格特征的复杂计算,因此该方法适用于各种畸形牙和牙弓拥挤问题的牙颌模型,对于提高口腔正畸治疗效果具有重要意义。
附图说明
图1为本发明的牙齿分割技术流程图。
图2为牙齿特征区域示意图。
图3a、图3b、图3c、图3d为拟合牙弓曲线示意图。
图4为牙龈线搜索结果示意图。
图5为牙龈线分割结果示意图。
图6为牙缝线搜索结果示意图。
图7为牙齿分割结果示意图。
具体实施方式
下面针对附图对本发明的实施例作进一步说明:
如图1所示,本发明的主要包括六个步骤:1、选择牙齿解剖特征点以确定咬合平面;2、利用牙颌模型的曲率信息提取牙齿特征区域;3、投影特征区域以拟合牙弓线;4、基于二次路径规划算法搜索牙龈线;5、分割牙龈线并搜索牙缝线;6、组合牙龈线和牙缝线构成牙齿分割线,并基于区域生长算法分割牙齿。
步骤1、确定咬合平面:在三维牙颌模型上交互式的选择牙齿解剖特征点:第一磨牙的近中舌尖点、中切牙的切缘点,然后利用主成分分析法拟合这些点以确定咬合平面,并进行坐标变换以使模型的XY平面对齐到咬合平面,Y轴穿过两中切牙牙缝。
步骤2、提取特征区域F:首先对三维牙颌模型进行离散曲率估算,并对曲率值进行基于直方图均衡化的拉伸变换以突出凹凸特征,然后对曲率提取的顶点集合进行连通性过滤以剔除噪音,从而获取牙齿和牙齿以及牙齿和牙龈分界的特征区域F,具体为:
(1)利用局部三次曲面拟合法计算每个网格顶点的最大主曲率值,并将所有顶点的曲率值都归一化到[0,1],然后将归一化后的曲率值集合{ki,i=1,2,…,N(N为网格顶点总数)}平均划分为L个区间:计算曲率在每个区间的累计分布概率密度:
其中,表示曲率值满足的网格顶点个数。然后,每个网格顶点的曲率ki按照如下公式进行变换:
(2)利用曲率值提取三维牙颌模型的初始特征区域F0:
F0={pi|k(pi)≥H,(i=1,2,...,N)}
其中,k(pi)为网格顶点pi变换后曲率值,H为曲率阈值。将初始特征区域F0根据网格顶点之间的连接关系映射为无向连通图G0,G0中边edge(vi,vj)的权值按如下规则进行定义:
其中,E为三维牙颌模型上所有三角边的集合。计算图G0的联通分量,并统计每个联通分量中的顶点个数,保留具有最多网格顶点个数的连通分量,即可删除噪音分量,得到牙齿和牙龈以及相邻牙齿之间分界的特征区域F,如图2所示。
步骤3、拟合牙弓线:将特征区域F投影到咬合平面上以构建特征二值图像,对特征二值图像进行形态学操作以获取牙弓曲线,并利用牙弓曲线将特征区域进行分区。具体为:
(1)将特征区域F投影到咬合平面上,并构建特征二值图像I(图3(a)所示)。对特征二值图像I进行孔洞填充(图3(b)所示)和细化操作,获取牙弓骨架线(图3(c)所示)。将牙弓骨架线利用最小二乘法线拟合为四次多项式曲线:y(x)=ax4+bx3+cx2+dx+e,从而得到牙弓曲线arch。
(2)计算特征二值图像I的外轮廓,并搜索外轮廓与牙弓曲线的交点,然后映射回到三维特征区域F当中,得到两个特征点m、n。以这两个特征点为界,内侧特征区域为舌侧区域,外侧为颊侧区域,如图3(d)所示。
步骤4、搜索牙龈线:为将特征区域F薄化为具有一个顶点宽度的特征线,将F映射为无向连通图,并基于二次路径规划算法搜索牙龈线,具体为:
(1)首先将特征区域F根据顶点之间的连接关系映射为无向连通图G1,G1中边的权值edge(vi,vj)定义为三维牙颌模型上两顶点vi,vj之间的欧式距离,即:
利用Dijkstra算法搜索无向连通图G1从点m到点n的最优路径,即为舌侧牙龈线。然后将搜索到的舌侧牙龈线邻域附近的边权值赋为无穷大,再次搜索从点m到点n的最优路径,即为颊侧牙龈线。
(2)为使牙龈线准确定位到牙齿和牙龈分界的最凹处,对上一步搜索的牙龈线进行修正。首先找到舌侧和颊侧牙龈线邻域附近的特征区域F’,并根据顶点连接关系映射为无向连通图G2,无向连通图G2中边的权值定义为两点欧式距离和曲率的乘积,即:
其中,k(vi)和k(vj)分别表示顶点vi和vj的曲率值,||vi-vj||表示两网格顶点之间的欧式距离,然后重新利用Dijkstra算法搜索舌侧牙龈线和颊侧牙龈线,其结果如图4所示。
步骤5、搜索牙缝线:首先根据牙龈线和牙弓曲线的特征信息找到牙龈线在牙缝处弯曲的位置,从而将整个三维牙颌模型的牙龈线分割为每颗牙齿的牙龈线,然后搜索每颗牙齿牙龈线两侧的最优路径,获得牙缝线。具体为:
(1)将步骤4得到的舌侧牙龈线和颊侧牙龈线利用能量法分别拟合为B样条曲线:flingual和fbuccul,并计算B样条曲线的曲率:ρlingual和ρbuccal,以及B样条曲线的二阶导数:和然后结合牙弓曲线的特征信息来确定舌侧牙龈线和颊侧牙龈线在牙缝弯曲处的位置:
其中,T为牙龈线的曲率阈值,archi为牙弓曲线上距离待测牙龈线的最近点,为牙弓曲线上点archi的法向量。将牙龈曲线在牙缝处的拐点删除,并将剩余的舌侧牙龈线和颊侧牙龈线上的点按顺序进行组合,以获得每颗牙齿的牙龈分割线,如图5所示。
(2)利用步骤4所构建的无向连通图G2,基于Dijkstra算法分别搜索每颗牙齿的舌侧牙龈线起点到颊侧牙龈线起点的最优路径,以及舌侧牙龈线终点到颊侧牙龈线终点的最优路径,以获得每颗牙齿两侧的牙缝线,其结果如图6所示。
步骤6、分割牙齿:将每颗牙齿的牙龈线和牙缝线组合,以构成每颗牙齿封闭的分割线,然后利用区域生长算法将每个牙齿从三维牙颌模型上分割出来。具体为:计算每颗牙齿分割线的质心,将三维牙颌模型上与该质心的x、y坐标最相近的点作为种子点,不断搜索种子点的邻域,直到达到当前牙齿的分割线为止,最终分割后结果如图7所示。
实施例不应视为对发明的限制,但任何基于本发明的精神所作的改进,都应在本发明的保护范围之内。