CN103871097B - 基于牙齿预备体的数据柔性融合方法 - Google Patents

Abstract

一种基于牙齿预备体的数据柔性融合方法，属于计算机辅助设计领域。本方法包括数据采集与预处理，特征曲线自适应提取技术，数据预配准和数据动态配准四个步骤。本发明的方法只包含有简单的人机交互操作，自动化程度高，融合后的模型不仅具有高精度表面信息，而且具有完善的内部结构，真实反映物理模型的完整信息，可以提高后续设计制造的效率和精度，在计算机辅助口腔设计领域有着重要的应用价值。

Description

基于牙齿预备体的数据柔性融合方法

所属技术领域

本发明属于计算机辅助设计领域，涉及逆向工程和生物医学工程领域。

背景技术

随着各领域科学技术的飞速发展，信息采集应用范围不断扩大，待处理的数据量也越来越大，来自超级计算机、地球卫星、各类扫描仪、各类摄像系统以及核磁共振等新兴领域的数据量都需要进行实时迅速的处理。数据源的多样性导致了数据操作的复杂性，在实际应用中，很多情况下单一数据已难以满足需求。

医学领域中，随着现代医学中计算机应用的不断发展，已可用数字化技术对口腔牙体缺损患者进行信息采集与三维重建，主要包括计算机断层扫描(CT)、激光扫描、结构光测量、三维立体摄影等。

CT是一种先进的无损检测（Nondestructive Testing，NDT）技术，在工业及医学等方面得到广泛应用，目前已经发展到第五代CT机。结构光扫描测量利用光栅投射在被测物体和参考平面上来测量三维形状，精度较高一般可到0.02～0.05mm，但无法获取模型内部数据信息。而CT扫描为断面扫描，常选用横断面，可在1～10mm范围内选择合适的层厚与层间距进行扫描，具有无层面外组织干扰、结构清晰的优点，但由于系统原理以及受人体所能接受放射剂量限制，CT精度通常只在0.2～0.4mm，这就导致三维重建出的模型虽然可以同时获得物体内外结构数据，但表面精度相对较低，丢失很多细节特征。

医学要求数字模型需要具有真实牙齿的内外结构，传统的单一数据源建模已难以满足设计制造的需求，而通过多源数据融合的方法可以对获得的几何模型进行优化组合得到物理模型的完整信息。

发明内容

针对现有方法的不足，本发明提出一种利用计算机辅助设计技术进行牙齿预备体的数据柔性融合方法，目的在于提高数字模型的精度，反映其真实物理内外结构，提高医生设计预备体模型的效率。

其方法由以下步骤组成：

步骤1、数据采集与预处理：通过CT扫描得到患者牙颌二维切片数据，然后取模灌制上下颌石膏模型，用结构光学扫描仪对石膏模型进行数据采集得到三维点云数据；通过对采集的两类数据进行三维数据重建，构建出上下颌精确的三角网格模型；由于扫描中不可避免的遇到缺陷，因此同时对重建模型进行去噪音、三角化、孔洞修复在内的预处理；

步骤2、融合区域提取：通过结构光扫描的整个牙颌模型，对于设计过程并不全是有效区域，首先需要将模型进行区域划分，分为融合区域和非融合区域，对融合区域曲面进行操作处理，非融合区域则约束固定；首先求解结构光模型的特征信息得到牙齿模型特征点，其次采用基于启发式搜索策略的龈缘线提取方法，裁剪提取融合区域；

步骤3、数据预配准：由于两种数据源采集方式不同，采集时并不在同一坐标系下，首先需要对数据进行初始融合，将两个模型转换到相同坐标系下。配准技术分为初始配准和精确配准两步，初始配准采用基于三点定位原理的人工交互配准，精确配准采用迭代最近点算法实现，初始配准是为了缩小模型之间错位，提供给后续精确配准良好的初始值以提高精确配准的效率和趋向。

步骤4、数据动态融合：主要包括通过迭代变形的方法，将CT模型的表面数据变形调整为结构光学模型的高精度表面。由于配准仅是两者之间在空间位置上的刚性对齐，所以需要采用基于微分属性的网格变形技术，利用网格顶点的微分属性隐式的描述模型的变形过程，以结构光学扫描数据为参照，运用拉普拉斯变形方法进一步调整优化CT模型的表面形态。融合过程中，通过将显著特征点定义为约束点及设计权重系数的方法进行优化。

通过两种数据源间的柔性融合方法，建立具有牙齿高精度表面，牙沟、牙窝等形态特征明显的数字模型，同时模型具有完善的内部组织结构，如牙髓腔等结构得以保留。融合过程平滑稳定，无交错、体积收缩现象。融合算法简单高效，具有较强的实用性。在齿科修复过程中，临床要求的修复体精度范围为0.05～0.1mm，融合额精度能够满足修复体的临床设计要求。在后续的修复体设计过程中，可以直接选取融合模型进行设计工作，同时可以在设计过程中实时观察牙齿的内部结构，进一步满足了口腔修复学的设计需求。

附图说明

图1磨牙CT扫描三角网格模型；

图2下颌结构光学扫描三角网格模型；

图3龈缘特征线提取；

图4融合区域裁剪提取；

图5 CT模型特征点拾取；

图6结构光学扫描模型特征点拾取；

图7模型精确配准结果；

图8模型第一次融合结果；

图9模型第二次融合结果；

图10模型第三次融合结果；

图11本发明流程图；

图中标号名称：1.龈缘特征线；2.颌面特征；3.初始配准特征点。

具体实施方案

步骤1、数据采集与预处理：通过CT扫描得到患者牙颌二维切片数据，然后取模灌制上下颌石膏模型，用结构光学扫描仪对石膏模型进行数据采集得到三维点云数据。通过对采集的两类数据进行三维数据重建，构建出上下颌精确的三角网格模型。由于扫描中不可避免的遇到噪声、扫描盲区等缺陷，因此同时需要对重建模型进行预处理（去噪音、三角化、孔洞修复等），CT三维重建与结构光点云数据三维重建分别如图1、图2所示。

步骤2、融合区域提取：通过结构光扫描的整个牙颌模型，对于设计过程并不全是有效区域，首先需要将模型进行区域划分，分为融合区域和非融合区域，对融合区域曲面进行操作处理，非融合区域则约束固定。首先求解结构光模型的特征信息得到牙齿模型特征点（脊点和谷点），其次采用基于启发式搜索策略的龈缘线提取方法，裁剪提取融合区域，具体过程如下：

(1)通过极值系数ε反映主曲率沿主方向上的变化率，将ε＝0处的点定义为

曲面的特征点，ε_max称为最大极值系数。进一步计算ε在主方向上的导数，将特征点划分为脊点与谷点。通过局部三次拟合曲面计算牙齿网格模型的离散微分量，在网格边v_iv_j上利用公式插值出所有的特征点p_v，其中v_i,v_j分别为网格边的两个顶点，依次连接模型上具有相邻关系的特征点。

(2)在搜索目标节点过程中，采用基于启发式的搜索策略，利用估值函数f＝ω_af_dir1+ω_bf_D+ω_cf_dir2+ω_df_C来评估所有节点的重要性，其中f_dir1，f_dir2表示方向函数，f_dir1表示当前搜索方向与前一搜索方向间的夹角，f_dir2表示下一搜索方向与目标点间的夹角关系；f_D表示距离函数，反映了当前搜索点与下一搜索点间的距离关系；f_C表示曲率函数，反映了当前节点与下一搜索点间曲率的差异，同时定义函数权重系数ω_a、ω_b、ω_c、ω_d来分别控制代价函数对于搜索路径结果的影响。通过启发函数可以提高整体搜索效率，f越小表示当前路径点与最优路径点越接近，将最先搜索到的最优点作为新的起始搜索点，迭代搜索至目标点，将搜索到的所有最优点连接构成最优路径。提取得到的龈缘线如图3所示，裁剪提取得到的融合区域如图4所示。

步骤3、数据预配准：由于两种数据源采集方式不同，采集时并不在同一坐标系下，首先需要对数据进行初始融合，将两个模型转换到相同坐标系下。配准技术分为初始配准和精确配准两步：

(1)初始配准采用基于三点定位原理的人工交互配准。根据三点重定位原理，不在同一直线上的三点可以定义一个坐标关系。为了实现初始配准，需要在不同模型上拾取彼此相互对应的至少三对特征点，如图5、图6所示，根据基于三点的坐标变换算法计算初始配准变换矩阵R和T，将矩阵作用于模型进行初始位置变换，完成两模型坐标系的归一化。

(2)精确配准采用迭代最近点算法实现。在每次迭代的过程中,对CT模型中的每一点v_i(i＝1,2,…,n)，在结构光扫描模型中寻找欧式距离最近点v_i′作为对应点，通过计算第k次迭代的最优旋转矩阵R^k与平移向量T^k，使得CT模型经过空间变化后与扫描模型M的二乘逼近目标函数

达到最小。通过刚体变换(旋转和平移)将CT模型定位到结构光学扫描三角网格模型当中，配准结果如图7所示。

步骤4、数据动态融合：采用基于拉普拉斯迭代变形的方法，将CT模型的表面数据变形调整为结构光学模型的高精度表面，融合过程中，通过将显著特征点定义为约束点及设计权重系数的方法进行优化，具体如下：

(1)首先采用将网格显著特征与莫尔斯理论相结合的方法实现扫描模型上显著特征点c_j的提取，将特征点c_j用作变形约束点。其次通过拉普拉斯方法，最小化对应能量函数 $E\left(V^{\′}\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\left|\right|{\mathrm{\δ}}_i^{\′}-{\mathrm{\δ}}_i\left|\right|}^2+\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\left|\right|v_j^{\′}-c_j\left|\right|}^2$ 求解变形后位置顶点坐标，δ_i,δ_i′分别为CT模型变形前后拉普拉斯坐标，v_j′为CT模型变形后笛卡尔坐标，n为CT模型顶点个数，m为提取的特征点个数，V′表示模型变形后的顶点向量。

(2)融合过程中，为保证图4中模型的牙尖、牙窝等主要颌面特征在变形后不会丢失，设计了一种包含距离权重函数w_d与角度权重函w_a的新的权重系数w_i，其中w_i＝w_d×w_a用来量化和评估点对的接近程度。同时，将变形总量分解为t次执行，每次变形量b＝1/t，迭代变形可以保证变形后网格自然平缓，而且很大程度上避免了网格自交现象。设定迭代次数为3进行示例说明，迭代融合过程分别如图8、图9、图10所示。

Claims (4)

1.一种基于牙齿预备体的数据柔性融合方法，其特征在于包括以下步骤：

步骤1、数据采集与预处理：通过CT扫描得到患者牙颌二维切片数据，然后取模灌制上下颌石膏模型，用结构光学扫描仪对石膏模型进行数据采集得到三维点云数据；通过对采集的两类数据进行三维数据重建，构建出上下颌精确的三角网格模型；由于扫描中不可避免的遇到缺陷，因此同时对重建模型进行去噪音、三角化、孔洞修复在内的预处理；

步骤2、融合区域提取：通过结构光扫描的整个牙颌模型，对于设计过程并不全是有效区域，首先需要将模型进行区域划分，分为融合区域和非融合区域，对融合区域曲面进行操作处理，非融合区域则约束固定；首先求解结构光模型的特征信息得到牙齿模型特征点，其次采用基于启发式搜索策略的龈缘线提取方法，裁剪提取融合区域；

步骤3、数据预配准：由于两种数据源采集方式不同，采集时并不在同一坐标系下，首先需要对数据进行初始融合，将两个模型转换到相同坐标系下；配准分为初始配准和精确配准两步；初始配准用于缩小模型之间错位，采用基于三点定位原理的人工交互配准；精确配准采用迭代最近点算法实现；

步骤4、数据动态融合：主要包括通过迭代变形的方法，将CT模型的表面数据变形调整为结构光学模型的高精度表面；由于配准仅是两者之间在空间位置上的刚性对齐，所以需要采用基于微分属性的网格变形技术，利用网格顶点的微分属性隐式的描述模型的变形过程，以结构光学扫描数据为参照，运用拉普拉斯变形方法进一步调整优化CT模型的表面形态；融合过程中，通过将显著特征点定义为约束点及设计权重系数的方法进行优化。

2.根据权利要求1所述的基于牙齿预备体的数据柔性融合方法，其特征在于：所述步骤2具体过程如下：

步骤2-1、通过极值系数ε反映主曲率沿主方向上的变化率，将ε＝0处的点定义为曲面的特征点，ε_max称为最大极值系数；进一步计算ε在主方向上的导数，将特征点划分为脊点与谷点；通过局部三次拟合曲面计算牙齿网格模型的离散微分量，在网格边v_iv_j上利用公式插值出所有的特征点p_v，其中v_i,v_j分别为网格边的两个顶点，依次连接模型上具有相邻关系的特征点；

步骤2-2、在搜索目标节点过程中，采用基于启发式的搜索策略，利用估值函数f＝ω_af_dir1+ω_bf_D+ω_cf_dir2+ω_df_C来评估所有节点的重要性，其中f_dir1，f_dir2表示方向函数，f_dir1表示当前搜索方向与前一搜索方向间的夹角，f_dir2表示下一搜索方向与目标点间的夹角关系；f_D表示距离函数，反映了当前搜索点与下一搜索点间的距离关系；f_C表示曲率函数，反映了当前节点与下一搜索点间曲率的差异，同时定义函数权重系数ω_a、ω_b、ω_c、ω_d来分别控制代价函数对于搜索路径结果的影响；通过启发函数可以提高整体搜索效率，f越小表示当前路径点与最优路径点越接近，将最先搜索到的最优点作为新的起始搜索点，迭代搜索至目标点，将搜索到的所有最优点连接构成最优路径。

3.根据权利要求1所述的基于牙齿预备体的数据柔性融合方法，其特征在于：所述步骤3具体过程如下：

步骤3-1、初始配准采用基于三点定位原理的人工交互配准；根据三点重定位原理，不在同一直线上的三点可以定义一个坐标关系；为了实现初始配准，需要在不同模型上拾取彼此相互对应的至少三对特征点，根据基于三点的坐标变换算法计算初始配准变换矩阵R和T，将矩阵作用于模型进行初始位置变换，完成两模型坐标系的归一化；

步骤3-2、精确配准采用迭代最近点算法实现；在每次迭代的过程中,对CT模型中的每一点v_i，在结构光扫描模型中寻找欧式距离最近点v′_i作为对应点，其中i＝1，2，···，n；通过计算第k次迭代的最优旋转矩阵R^k与平移向量T^k，使得CT模型经过空间变化后与扫描模型M的二乘逼近目标函数达到最小；通过旋转和平移变换将CT模型定位到结构光学扫描三角网格模型当中。

4.根据权利要求1所述的基于牙齿预备体的数据柔性融合方法，其特征在于：所述步骤4具体过程如下：

步骤4-1、首先采用将网格显著特征与莫尔斯理论相结合的方法实现扫描模型上显著特征点c_j的提取，将特征点c_j用作变形约束点；其次通过拉普拉斯方法，最小化对应能量函数求解变形后位置顶点坐标，δ_i,δ′_i分别为CT模型变形前后拉普拉斯坐标，v′_j为CT模型变形后笛卡尔坐标，n为CT模型顶点个数，m为提取的特征点个数，V′表示模型变形后的顶点向量；

步骤4-2、融合过程中，为保证模型的牙尖、牙窝在内的主要颌面特征在变形后不会丢失，设计了一种包含距离权重函数w_d与角度权重函w_a的新的权重系数w_i，其中w_i＝w_d×w_a用来量化和评估点对的接近程度；同时，将变形总量分解为t次执行，每次变形量b＝1/t，迭代变形可以保证变形后网格自然平缓，而且很大程度上避免了网格自交现象。