具体实施方式

本发明提出的正畸矫治过程中的牙齿矫治状态的检测方法，包括：

获取牙齿矫治状态信息；

对数字化口腔网格模型进行预处理，得到数字化口腔网格模型；

使用基于投影的二维-三维配准方法，对每颗牙齿的三维网格模型进行刚性变换,使得每颗牙齿的三维网格模型投影到该照片的二维平面上，将每颗牙齿的二维投影牙齿轮廓与患者牙齿轮廓进行配准，得到该照片坐标系下每颗牙齿的真实刚性变换矩阵；

根据每张照片坐标系下牙齿的真实刚性变换矩阵，获取两张照片坐标系的关联矩阵；

根据坐标系关联矩阵和刚性变换矩阵，获得整副牙齿模型，实现正畸矫治过程中的牙齿矫治状态的检测。

以下结合附图，详细介绍本发明的内容：

本发明提出的正畸矫治过程中的牙齿矫治状态的检测方法，其流程框图如图1所示，包括以下步骤：

(1)获取牙齿矫治状态信息，具体过程如下：

(1-1)在隐形正畸治疗前，使用三维扫描枪对患者的牙齿进行扫描，获取患者的数字化口腔网格模型，其中数字化网格模型使用行业常用的obj(3D模型文件格式)文件格式存储，数字化口腔网格模型中存储网格模型的三维空间点、线和面的位置信息；

(1-2)从不同角度拍摄带有标签示意图的患者牙齿照片，口腔照片使用手机拍摄，患者牙齿照片为彩色的M*N二维图像，标签示意图是由黑色和白色的类似于二维码的正方形图像，M为二维图像的宽度，N为二维图像的高度；

(2)对步骤(1)的数字化口腔网格模型进行预处理，具体过程如下：

(2-1)对患者在隐形正畸前的口内扫描网格模型进行分牙操作。分牙操作为对步骤(1-1)的数字化口腔网格模型进行分割，得到每颗牙齿的三维网格模型，并根据正畸学对现实中牙齿的编号，并对每颗牙齿的三维网格模型进行编号，得到编号为j的牙冠网格模型A_j；

获取拍摄照片的相机内参数(相机内参数是与相机自身特性相关的参数,比如相机的焦距、像素大小等)。

(2-2)对步骤(1-2)的患者牙齿照片进行去畸变处理，得到去畸变处理后的患者口腔M*N二维图像；使用图像处理工具OpenCV(一个基于BSD许可发行的跨平台计算机视觉和机器学习软件库)，对拍摄的每一张患者照片进行去畸变处理。

(2-3)获取步骤(1-2)中拍摄照片时的相机内参数和相机外参数；(相机外参数是在世界坐标系中的参数，比如相机的位置、旋转方向等)。在对患者牙齿进行拍照时，用定制的标签(1cm*1cm，类似于二维码)的标签贴在患者嘴唇上，使得标签能出现在照片中。根据照片中标签的位置，来获取相机的外参数，该方法为公知公用方法。

(2-4)从步骤(2-2)的去畸变处理后的患者口腔M*N二维图像中提取患者牙齿轮廓，具体过程如下：

(2-4-1)采用Canny边缘检测方法，对步骤(2-2)患者口腔M*N二维图像进行轮廓边缘检测，得到带有边缘信息的M*N二维图像；

本发明提出的低成本牙齿矫治状态检测方法中，主要使用患者拍摄的牙齿照片中的轮廓作为二维信息，来进行二维三维配准。因此需要对患者提交的牙齿照片进行轮廓提取。

关于照片中牙齿轮廓的处理，本发明使用多种图片轮廓提取方法提取照片中牙齿的轮廓，然后将这些方法的结果进行融合，获得牙齿的轮廓(牙齿轮廓指牙冠处的轮廓，不包括牙齿和牙龈交界线处的轮廓

口腔内部光线相对不足，而且患者的拍摄设备一般为手机等设备，拍摄的牙齿照片容易出现很多噪点，因为Canny算法不容易受到照片中噪点的影响，能够准确地提取边缘，所以本发明使用Canny算法对照片中的轮廓边缘进行提取。

(2-4-2)对步骤(2-2)患者口腔M*N二维图像进行Garbor滤波处理，得到带有轮廓信息的M*N二维图像；

牙齿轮廓形状是特殊的(大部分由横着或者竖着的弧线组成)，而且因为患者在拍摄口内牙齿照片时，光线的变化非常大(例如有无开闪光灯、室内或室外拍摄、光源是否有遮挡等)，因此提取患者的牙齿轮廓需要考虑照片拍摄时的光照。因为Garbor滤波对于光照不敏感，而且能够从频域提取空域难以提取的信息，能够选择提取轮廓的方向。因此，本发明使用Garbor滤波从0度、45度、90度、135度四种方向对患者提交的照片进行边缘轮廓提取。

(2-4-3)利用生成式边缘学习方法(Boost Edge Learning(简称BEL))对步骤(2-2)的患者口腔M*N二维图像进行处理，得到只有牙齿轮廓信息的M*N二维图像；

因为Canny算法和Garbor滤波是对照片中全局的边缘轮廓进行提取，而本发明需要提取的是患者的牙齿轮廓(不包括牙龈和牙齿的边缘)。所以需要有一种边缘轮廓提取算法，来提取患者的牙齿轮廓。因为本发明的牙齿轮廓训练数据数据量相对较少，所以本发明使用BEL方法，用概率提升树，通过少量的训练样本，来学习患者的牙齿轮廓。

BEL算法在具体实现时，对图像中的每个点独立地判断其是否是轮廓点。在判断时，以此点为中心，使用一个非常大的孔径(本发明中为49*49)作为判断的背景。在学习阶段，BEL算法选择并组合不同尺度的大量特征，以便使用概率提升树分类算法的扩展版本来学习辨别模型。

(2-4-4)采用轮廓融合提取方法，将上述步骤(2-4-1)、步骤(2-4-2)和步骤(2-4-3)得到的三张M*N二维图像进行融合，得到带有患者牙齿轮廓的M*N二维图像，包括以下步骤：

Canny算法具有良好的对噪声不敏感性和能进行准确边缘提取，但是会提取全局边缘轮廓；Garbor滤波对患者拍摄牙齿时的光照不敏感，对于轮廓的方向有很好的选择，但是也会提取边缘轮廓；BEL算法只提取患者牙齿的轮廓，但是准确性有待提高，且对噪声相对敏感。通过三种方法的融合，扬长避短，获得各个方法的优势，就可以获得较为准确的牙齿轮廓。

(2-4-4-1)设定图像中像素值的阈值，对步骤(2-4-2)的Garbor滤波后得到的M*N二维图像进行二值化处理，将该二值化处理后的二维图像中像素值高于阈值的边缘判定为边缘轮廓；

(2-4-4-2)将步骤(2-4-4-1)、步骤(2-4-1)和步骤(2-4-3)的M*N二维图像进行“交”运算，即步骤(2-4-4-1)、步骤(2-4-1)和步骤(2-4-3)同时判定为边缘轮廓的点才确定为牙齿轮廓的边缘点；

(2-5)对步骤(2-4)得到的带有患者牙齿轮廓的M*N二维图像中的每颗牙齿轮廓进行编号，即得到B_ij，B_ij表示第i张牙齿照片中编号为j牙齿的轮廓；

因为每个人的牙齿形状都不尽相同，且从不同角度拍摄牙齿得到的牙齿轮廓也不尽相同，所以通过机器学习来对得到的牙齿轮廓进行编号是一件相对困难的事情。而且前文提到的轮廓提取融合算法也无法做到100％的准确率，所以需要借助人工对前文轮廓提取融合算法得到的牙齿轮廓进行简单修饰并依据牙位表对牙齿轮廓进行编号。

为了确保下一步二维三维配准的精确度，在人工对牙齿轮廓进行编号时，可以舍去轮廓模糊的牙齿，而保留轮廓清晰的牙齿轮廓。虽然本步骤仍需要进行人工交互，但是和人工进行轮廓提取相比，本发明的方法已经大大简化了人工交互的复杂度，人工只需要在轮廓提取融合算法的基础上进行编号和简单修饰即可。

(3)将步骤(2-5)的患者牙齿轮廓B_ij与步骤(2-1)的牙齿数字化三维模型A_j进行配准，(简称二维三维配准)，即使用基于投影的二维-三维配准方法，对步骤(2-1)的每颗牙齿的三维网格模型进行刚性变换,使得每颗牙齿的三维网格模型投影到二维平面上，投影得到每颗牙齿的M*N二维投影牙齿轮廓，(该牙齿轮廓不包含牙齿牙龈交接线)，将每颗牙齿的M*N二维投影牙齿轮廓与步骤(2-4)的患者牙齿轮廓进行配准，得到每颗牙齿的真实刚性变换矩阵其中，i表示步骤(2-4)中第i张照片，j表示步骤(2-4)中牙齿的编号，获取第i张照片中编号为j的牙齿的真实刚性变换矩阵包括以下步骤：

(3-1)分别对每一张照片的牙齿轮廓，以上颌面牙列和下颌面牙列为整体进行二维三维配准，即根据步骤(2-5)得到第i张照片的每颗牙齿轮廓B_ij，按照牙齿所在的上颌面或下颌面，分成两个颌面的牙列轮廓，将上颌面牙列轮廓或下颌面牙列轮廓分别与步骤(2-1)得到的牙齿数字化三维模型A_j进行二维三维配准，得到第i张牙齿照片中相应颌面牙齿的刚性变换矩阵T′_i，T_i表示第i张牙齿照片中颌面牙齿的刚性变换矩阵理论值，T′_i为T_i的计算值；其中，上颌面牙列处理步骤和下颌面牙列处理步骤完全一致。整副牙列二维三维配准的步骤用以下公式表示：

使用第i张牙齿照片中整个颌面牙齿的刚性变换矩阵T_i，对步骤(2-1)得到的牙齿数字化三维模型A_j进行刚性变化，得到刚性变化后牙齿三维模型其中，其中N_i为牙齿照片i中上颌面或下颌面上被标出的牙齿轮廓数量；刚性变换矩阵T由六个参数构成，分别是t_x、t_y、t_z、α、β、γ，其中，t_x、t_y、t_z表示空间内平移量，α、β、γ表示空间内旋转角度，其中T和t_x、t_y、t_z、α、β、γ的关系为：

整副牙列的牙齿配准过程如下：

(3-1-1)使用函数对得到的刚性变化后牙齿三维模型从三维网格模型投影到二维平面，得到二维投影图形其中，函数由步骤(2-3)获取的相机内参数和相机外参数决定；本发明中使用开放式图形库OpenGL进行投影，通过相机内参数和相机外参数，将三维网格模型投影至二维成像平面，此方法为本技术领域的公知技术。

(3-1-2)使用函数C,对得到的二维投影图形进行牙齿轮廓提取，得到投影牙齿轮廓本发明中，将步骤(3-1-1)投影得到的二维图形二值化(即背景为黑色，牙齿的投影为白色)，再使用图形处理软件OpenCV进行投影牙齿轮廓提取，此轮廓提取方法也为本技术领域的公知技术。

(3-1-3)使用在二维层面两个轮廓的距离度量函数F^2D，对投影牙齿轮廓和步骤(2-5)得到的牙齿轮廓B_ij进行距离度量，获得度量距离

在本发明实现中，采用最近点对的平均距离(欧几里得距离)作为两个轮廓的距离。用{x^project}表示某颗牙齿投影生成的轮廓点点集，K^project投影轮廓点的数量；用{x^real}表示对应牙齿在照片中真实轮廓的点集，K^real表示真实轮廓点的数量；则距离度量函数F^2D可以用下式表示：

其中，是真实轮廓点在投影轮廓点集中最近的对应点。

(3-1-4)遍历一副牙列上的所有牙齿，对每一颗牙齿重复上述步骤(3-1-1)、(3-1-2)和(3-1-3)，得到所有每颗牙齿的度量距离，求所有每颗牙齿的度量距离的平均值，得到整副牙列的距离度量均值其中N_i为牙齿照片i中上颌面或下颌面上被标出的牙齿轮廓数量；

(3-1-5)使用最优化方法，不断调整牙列的刚性变换矩阵T_i，使得牙齿数字化三维模型A_j的投影轮廓与步骤(2-5)得到的整副牙列轮廓的度量距离最小，当度量距离最小时，得到此时的刚性变换矩阵T′_i；

(3-2)根据步骤(3-1)的整副牙列刚性变换矩阵T′_i，对每一颗牙齿进行二维三维配准。即根据步骤(3-1)的整副牙列刚性变换矩阵T′_i，将步骤(2-5)得到第i张照片的每颗牙齿轮廓B_ij和步骤(2-1)的牙齿数字化三维模型A_j,进行进一步刚性变化，得到进一步刚性变换的刚性变换矩阵T′_ij，单颗牙齿二维三维配准的步骤如以下公式所示：

其中，T_ij表示照片i中编号为j的牙齿，在整副牙列进行刚性变换T′_i的基础上进行的刚性变换。当单颗牙齿投影和真实轮廓距离最小时，得到刚性变化矩阵T′_ij；具体步骤如下：

(3-2-1)根据步骤(3-1)获得的整副牙列刚性变换矩阵T′_i和第i张照片编号j的牙齿数字化三维模型A_j，得到进一步刚性变换后的牙齿三维模型

(3-2-2)使用函数对得到的刚性变化后牙齿三维模型从三维网格模型投影到二维平面，得到二维投影图形其中，函数由步骤(2-3)获取的相机内参数和相机外参数决定；

(3-2-3)使用函数C,对得到的二维投影图形进行牙齿轮廓提取，得到投影牙齿轮廓

(3-2-4)使用在二维层面两个轮廓的评价函数(用于度量距离)F^2D，对投影牙齿轮廓和步骤(2-5)得到的牙齿轮廓B_ij进行距离度量，获得度量距离

(3-2-5)使用最优化方法，不断调整牙齿的刚性变换矩阵T_ij，使得牙齿数字化三维模型A_j的投影轮廓与步骤(2-5)得到的整副牙列轮廓的度量距离最小。当度量距离最小时，得到照片i上编号为j的牙齿的刚性变换矩阵T′_ij；

(3-2-6)对每颗牙齿均进行上述(3-2-1)、(3-2-2)、(3-2-3)、(3-2-4)、(3-2-5)操作，得到每颗牙齿在步骤(3-1)进行整副牙列刚性变换T′_i基础上进一步刚性变换的刚性变换矩阵T′_ij；

(3-3)根据步骤(3-1)获得的整副牙列刚性变换矩阵T′_i和步骤(3-2)获得的每颗牙齿进一步刚性变换的刚性变换矩阵T′_ij，获取第i张照片中编号为j的牙齿的真实刚性变换矩阵

(4)根据步骤(3)得到的获取第i张照片中编号为j的牙齿的真实刚性变换矩阵获取拥有n颗相同编号牙齿的照片i₁和照片i₂的坐标系关联矩阵包括以下步骤：

(4-1)获取照片i₁和照片i₂上具有相同编号的n颗牙齿，记为n颗共同牙齿；牙齿编号为口腔学通用的编号；

(4-2)在n颗共同编号牙齿中，以编号k牙齿为例，k不等于i,k不等于j，求出以编号为k牙齿为基准，在坐标系关联矩阵关联下，得到照片i₁和照片i₂的n颗共同牙齿的投影度量距离D_k，其中，表示从i₁变换到i₂坐标系的坐标系关联矩阵， 表示的逆矩阵，D_k的计算公式如下：

以该编号k牙齿为基准，在照片i₁的坐标系下，求出在照片i₂中的n颗共同牙齿的投影到照片i₁成像平台的投影度量距离获取投影度量距离的步骤如下：

(4-2-1)根据步骤(3-3)，获取照片i₂中编号为j的牙齿的真实刚性变换矩阵以及经过刚性变换后的三维牙齿网格模型

(4-2-2)根据照片i₁和照片i₂的坐标系关系矩阵从步骤(4-1-1)得到的照片i₂中编号为j的牙齿的三维网格模型得到转移到照片i₁的牙齿的三维网格模型

(4-2-3)使用函数对步骤(4-1-2)的刚性变化后牙齿的三维网络模型从三维网格模型投影到二维平面，得到二维投影图形

(4-2-4)使用函数C,对得到的二维投影图形进行牙齿轮廓提取，得到投影牙齿轮廓

(4-2-5)使用在二维层面两个轮廓的距离度量函数F^2D，对投影牙齿轮廓和步骤(2-5)得到的牙齿轮廓B_ij进行距离度量，获得度量距离

(4-2-6)以该编号k牙齿为基准，在照片i₁的坐标系下，求出在照片i₂中的n颗共同牙齿的投影到照片i₁成像平台的投影平均度量距离即

(4-2-7)重复步骤(4-1-2)、(4-1-3)、(4-1-4)、(4-1-5)、(4-1-6)，求出以该编号k牙齿为基准，在照片i₂的坐标系下，求出在照片i₁中的n颗共同编号牙齿的投影到照片i₁成像平台的投影平均度量距离即

(4-2-8)投影平均度量距离和投影平均度量距离相加，得到以该编号k牙齿为基准，照片i₁和照片i₂的N颗共同牙齿的投影度量距离D_k；

(4-2-9)重复步骤(4-2-1)至(4-2-8)，获得n颗共同牙齿中，以每一颗牙齿为基准的投影度量距离。

(4-3)在n颗共同牙齿中，根据步骤(4-2)得到以每一颗牙齿为基准的投影度量距离。从所有牙齿的度量距离中找到度量距离最小的牙齿编号，记为牙齿编号k，使用其(编号为k的牙齿)坐标系关联矩阵作为照片i₁和照片i₂的坐标系关联矩阵

(5)根据步骤(4)的相同编号的牙齿j的照片i₁和照片i₂的坐标系关联矩阵和步骤(3)获得的在照片i上编号j的牙齿的刚性变换矩阵获得整副牙齿模型，实现正畸矫治过程中的牙齿矫治状态的检测。

通过上述方法生成牙齿模型后，可以使用数字化正畸系统中的一些计量方式进行牙齿矫治状态的检测。牙齿矫治状态的具体测量根据不同正畸公司的要求有所不同，各个公司可以根据需求对生成的模型用PAR指数、相邻牙齿间距等状态测量指标来进行牙齿矫治状态检测。