CN105184760B - 牙体影像的接合方法 - Google Patents

Abstract

一种牙体影像的接合方法，先由一探头取得一第一影像与一第二影像，再由一影像处理装置接收并存储该两影像以进行影像接合动作，其中该影像处理装置首先计算该两影像之间的重复区域以及该第二影像的一第一位移量，接着获取该第二影像中的齿面特征点以及非齿面特征区域的样本点Bsi，接着计算Bsi在该第一影像中的对应点Asi；最后根据一旋转矩阵与一第二位移量调整该第二影像的位置以贴合第一影像，避免牙齿轮廓不连续的情况。通过本发明根据第一位移量以初步移动该第二影像，再根据该旋转矩阵与第二位移量以进行第一影像与第二影像彼此相对位置的微调，通过此微调手段可将第一影像与第二影像对齐，避免牙齿轮廓不连续的情况。

Description

牙体影像的接合方法

技术领域

本发明是关于一种牙体影像的接合方法，特别是指可根据多个牙体影像的重复区域彼此叠合而接合影像的方法。

背景技术

在建立口腔内牙体的影像时，牙科医生会持一探头伸入患者的口腔中以进行牙体的拍摄。以下颚为例，若要取得整排牙齿的舌侧面影像，则牙科医生需要对下颚的牙齿进行拍摄一系列的影像，例如沿着患者左边的大臼齿往右边的大臼齿方向的轨迹拍摄。于拍摄完成后，由与该探头连线的电脑装置接收这些影像，该电脑装置可分辨出这些影像彼此间的关联性，将该多个影像接合而得到完整的一下颚牙齿舌侧面影像。

然而，由于牙科医生进行拍照时，手持探头的拍照角度不尽相同，因此当该多个影像彼此接合后，牙齿轮廓有可能产生不连续的状况，例如牙缝的线条产生错位，造成不精确的牙齿影像。

发明内容

本发明的目的是提供一种牙体影像的接合方法，用以解决现有技术中接合影像不连续的问题。

本发明的技术方案是提供一种牙体影像的接合方法包含有：

由一探头取得一第一影像与一第二影像；

由一影像处理装置接收并存储这些影像，以计算该第一影像与第二影像之间的重复区域以及该第二影像的一第一位移量，以利用该第一位移量移动该第二影像；

获取该第二影像中的齿面特征点以及于非齿面特征区域的样本点Bsi；

计算该第二影像齿面特征点与样本点Bsi在第一影像中的对应点Asi；

根据Asi与Bsi建立一旋转矩阵与该第二影像的一第二位移量，以根据该旋转矩阵与第二位移量移动该第二影像与第一影像接合。

根据本发明的方法，根据第一位移量以初步移动该第二影像，再根据该旋转矩阵与第二位移量以进行第一影像与第二影像彼此相对位置的微调，通过此微调手段可将第一影像与第二影像对齐，避免牙齿轮廓不连续的情况。

附图说明

图1：本发明的系统方块示意图。

图2：本发明较佳实施例的流程示意图。

图3：本发明中第一影像的示意图。

图4：本发明中第二影像的示意图。

图5：本发明中第一影像的重复区域示意图。

图6：本发明中第二影像的重复区域示意图。

图7：本发明中第一影像与第二影像根据第一位移量接合的示意图。

图8：本发明中第一影像与第二影像进一步由第二位移量接合的示意图。

主要元件符号说明

10 探头

11 姿态感测器

20 影像处理装置

31 第一影像

310 重复区域

311 缺损部分

320 重复区域

32 第二影像

40 牙齿

41 牙缝处

42 牙缝处

具体实施方式

请参考图1所示，执行本发明方法的系统包含有一探头10与一影像处理装置20，该影像处理装置20可为电脑，连线到该探头10以进行数据传输。图2为本发明方法的流程图。

首先该探头10取得多个牙体的影像，各影像的三维数据以空间坐标参数所表示(步骤101)。该探头10可为手持式探头或固定式探头，以手持式探头为例，牙科医师将探头10伸入患者的口腔内以进行牙齿的拍摄作业，其中该探头10内设有一姿态感测器11，该姿态感测器11可为陀螺仪，用以产生该影像当下的空间坐标参数，在直角坐标系统下，空间坐标参数包含有x轴旋转索引值(yaw)、y轴旋转索引值(pitch)、z轴旋转索引值(roll)、x坐标(dx)、y坐标(dy)与z坐标(dz)，是以，每张拍摄的影像都有对应的空间坐标参数。以下颚牙齿为例，该探头10可沿着患者左侧的大臼齿、小臼齿、犬齿至门牙方向的轨迹对不同牙齿的舌侧面依序拍摄，进而得到多个具有连续关系的牙体影像。

以探头10拍摄得到多个影像后，由该影像处理装置20接收并存储这些影像与对应的空间数据参数(步骤102)，于取得这些影像后，该影像处理装置20负责将这些影像彼此接合。请参考图3与图4所示，本发明以两个接续拍摄的第一影像31与第二影像32的接合动作为例说明，该两影像31、32各自拍到相同的牙齿40，其中第一影像31与第二影像32为点云数据(point cloud)。

该影像处理装置20可在直角坐标系统进行影像接合的运算。在进行接合之前，请参考图5与图6，首先该影像处理装置20计算该第一影像31与第二影像32之间的重复区域310、320以及第二影像32的一第一位移量(步骤103)。所述重复区域310、320是指该两影像31、32所拍到相同物体(即牙齿40)的分布区域。在直角坐标系统中，该第一影像31与第二影像32主要是排列在X轴上，因此若要接合该两影像31、32，则是令该第一影像31朝第二影像32沿着X轴平移，或令第二影像32朝第一影像31沿着X轴平移，直到第一影像31与第二影像32的重复区域310、320彼此重叠，本较佳实施例是令第二影像32朝第一影像31沿着X轴平移。

假设平移前该第二影像32包含一参考点，该参考点的坐标为Bi，则平移后该参考点新坐标为Bnewi，则Bnewi＝Bi+Dx，Dx为第二影像32的第一位移量，其中i＝1，2，3，…，NB，NB为第二影像32点云数据的总点数。Dx是在[Dini-d，Dini+d]范围内搜寻取得，[Dini-d，Dini+d]的范围是由使用者操作该影像处理装置20而设定，Dini为初始位移量，d为搜寻范围。

关于第一位移量Dx的建立，本发明是根据三维物体(即第一影像31与第二影像32)的深度图(depth image)，深度图可以D(x，y)表示之，并计算深度图的梯度大小(gradientmagnitude)，本发明是将梯度沿着Y轴方向投影为例，进而分别计算该第一影像31与第二影像32于X坐标的梯度累计值。是以，该两影像31、32共产生两个梯度累计值。其中请参考图5，拍摄第一影像31时，因拍摄死角或唾液反光而导致第一影像31有缺损部分311，缺损部分311不列入运算，其梯度累计设为0。

本发明深度图的梯度大小表示为|D(x,y)|，通过滤波器计算如下，本发明是采用Sobel滤波器：

其中

假设在X坐标的梯度累计值表示为：

则该第一影像31与第二影像32的梯度累计值则可分别表示为AccA(x)与AccB(x)，该第一影像31与第二影像32的梯度累计值相关性可表示如下：

该第一位移量为

请参考图7所示，即是将该第二影像32根据第一位移量平移后，与第一影像31初步接合的示意图，由图5～图7中可见，该第一影像31与第二影像32的重复区域310、320已大致重叠，惟牙缝处41(由较密集的点所构成的区域代表牙缝)产生不连续的情形。

为了微调第一影像31与第二影像32的相对位置以达到对齐目的时，首先获取该第二影像32中的齿面特征点(步骤104)，接着获取该第二影像32在非齿面特征区域的样本点(步骤105)，如此一来，因为是利用获取出来的特征点与样本点进行后续计算，不是影像中的每一点数据进行计算，除了可以降低运算数据量，还可以提升较具代表性特征点的比例，以提高运算速度并增进比对的精确度。

在第104步骤中，齿面特征点为齿面的三维几何特征点，即具有明显几何变化的特征点。本发明可利用曲率计算三维几何变化。以S(p)(shape index)表示几何特征，考虑点云数据中的任一点p，则

S(p)介于0与1之间，其中k₁(p)与k₂(p)分别为该p点的两个主曲率(principalcurvature)，且k₁(p)≥k₂(p)；S(p)为0时代表p的周边形状是球杯状(spherical cup)，S(p)为1时代表p的周边形状是球冠状(spherical cap)，至于中间值0.5则表示p点为马鞍点(saddle point)，至于其他的值则介于这些形状之间。由于牙齿上包含有凹陷处、凸起处与马鞍处，故S(p)可分别表示牙齿的凹陷处、凸起处与马鞍处，其中S(p)为1～0.9时代表p的周边形状是球冠状，S(p)为0～0.1时代表p的周边形状是球杯状，S(p)为0.45～0.55时代表p点为马鞍点，以上S(p)的数值范围仅是较佳实施例而已，不用以限制本发明。因此，对该第二影像32的点云数据进行S(p)的计算，可取得该第二影像32特征点。

于齿面特征点取完后，进行第105步骤以在非齿面特征区域获取样本点，本发明可利用均匀取样手段(uniform sampling)对该第二影像32的非齿面特征区域获取样本点，该非齿面特征区域指该第二影像32中齿面特征点以外的区域。本较佳实施例中，定义三维坐标系的Z轴为深度，因此沿着X轴与Y轴，并于每间隔距离△x与△y进行取样，即可完成样本点的获取动作。

根据第104步骤与第105步骤分别取得齿面特征点与样本点后，接着计算第一影像31与第二影像32的对应关系。以第二影像32为例，取其特征点与样本点的联集后，其联集内包含有一点Bsi，i＝1，2，3…NBS，NBS为第二影像32齿面特征点数量与样本点数量的总和，若第一影像31包含有一点Asi，i＝1，2，3…NAS，NAS＝NBS，Asi为Bsi的对应点(corresponding point)，则(Asi，Bsi)形成两影像的对应关系。因为探头10在进行口内扫描拍照时，口腔内被遮蔽的地方无法被扫描到，例如从颊侧或舌侧拍摄时凹陷处易被遮蔽，导致第二影像32的一点Bsi无法与第一影像31对应而产生误差，因此本发明以内插手段修补第一影像31与第二影像32因遮蔽所产生的空缺处，以解决数据缺失造成的对应困难。本发明的内插手段可以是线性内插(linear interpolation)、三次内插(cubicinterpolation)或其他内插方式，所述内插计算为公知常识，再此不赘述。

计算第一影像31与第二影像32的对应关系时，以第二影像32中任一Bsi点为例，将该Bsi与第一影像31的所有点云数据进行Ak-d tree(Approximated k-d tree，Ak-d tree)演算(可参考Michael Greenspan等人所著的”Approximate K-D tree Search forEfficient ICP”，如附件)，计算Bsi与第一影像31所有点云数据的距离，进而决定该第一影像31中距离该Bsi最近的点作为对应的Asi，演算时限制两对应点Asi、Bsi间的距离在一门槛距离内，借以忽略距离过远的对应点，避免影响后续的估算，并节省运算量(步骤106)。前述实施例是先计算第二影像32的特征点与样本点后，再计算出与Bsi对应的Asi；于另一实施例中，亦可先计算第一影像31的特征点与样本点后，再计算出与Asi对应的Bsi。本发明是以Ak-d tree演算法计算最接近的Asi与Bsi，但不以此为限。

当Asi与Bsi的对应关系建立后，该影像处理装置20接着进行刚体转换(rigidtransformation)以微调第一影像31与第二影像32的位置(步骤107)，所述刚体转换包含旋转与位移的动作。各第一影像31与第二影像32通过合适的旋转与位移可彼此叠合，而旋转与位移通过三维点的互协方差矩阵(cross covariance matrix)M_AB进行计算，互协方差矩阵如下：

其中

为第一影像31齿面特征点坐标与样本点坐标的平均值，为第二影像32齿面特征点坐标与样本点坐标的平均值，前述互协方差矩阵建立后，产生如下的4×4矩阵：

其中tr(M_AB)代表矩阵M_AB主对角线的和，E(M_AB)最大特征值(eigenvalue)对应的特征向量(eigenvector)为使得均方差(mean square error)为最小的旋转。若E(M_AB)p＝λp的关系成立，p与λ分别为E(M_AB的特征向量与特征值，可由特征多项式(characteristicpolynomial)的根计算特征值，再代入E(M_AB)p＝λp，解线性方程组求得特征向量p。本发明以电脑执行数值分析方法，以Householder法将E(M_AB)转换为三对角线矩阵(tridiagonalmatrix)，再以QL演算法取得特征值与特征向量，可参考Numerical Recipes in C(ISBN:0521 43108 5)书中所提到的函式tred2()与tqli()计算特征向量与特征值(第108～109页、第113～115页、第469～480页)，以下说明Householder演算流程。

要找到特征值与特征向量，需先化简矩阵，若原始矩阵是对称矩阵，则可化简为三对角线矩阵(tridiagonal matrix)，也就是除了主对角线与其上下相邻的两对角线之外，元素均为零的矩阵，再以叠代的方式求解。我们使用的化简方法为Householder法，此法将一个n乘n的对称矩阵化简为tridiagonal matrix，首先考虑一Householder矩阵P，此矩阵的形式为P＝1-2w·wT，w为一个长度为一的向量，也就是|w|2＝1，由于P2＝1-4w·wT+4w·(wT·w)·wT＝1，所以P-1＝P；此外，P为对称矩阵，PT＝P，可得P-1＝PT，所以P为正交矩阵(orthogonal matrix)。

把P改写为

其中

假设向量x为要化简矩阵A的第一行(column)，令u＝x±|x|e₁，e₁＝[1,0,0...,0]^T为单位向量，则

其中X₁为X的第一个元素。这样的结果显示P矩阵把矩阵A的第一行(column)化简为第一元素不为零且其余元素为零的向量。如果略过矩阵A第一行的第一元素，而以其第二元素之后的部分形成向量x来建构Householder矩阵若以此Householder矩阵建立一n×n矩阵：

则

则也有类似的效果，此时P₁·A第一行第二元素之后均为零，为P₁·A扣除第一列与第一行所得的矩阵；若P₁·A·P₁ ^T(P₁ ^T＝P₁)，则可以让第一列第二元素之后变成零，而得到

此时A1的第一列与第一行的第二元素之后为零，为P₁·A·P₁扣除第一列与第一行所得的矩阵，由于P₁为正交矩阵，所以直接将上式的P₁ ^T写成P₁。接着，重复上述步骤，也就是取矩阵A第二行第三元素之后的部分形成向量来建构Householder矩阵若以此Householder矩阵建立一n×n矩阵：

则A₂＝P₂·A·P₂可得前两行与前两列除三个主对角线之外均为零的矩阵，重复上述步骤，可建立tridiagonal matrix。

而为方便运算，可重新整理流程，A·P可化简如下：

从A的最后一行开始建Householder矩阵，演算流程如下：

第m次叠代，m＝1，2，3，…，n-2，进行下列步骤：

1、建立for i＝n-m+1＝n，n-1，…，3

其中且前的正负号由a_i,i-1的正负号决定，应设定成与a_i,i-1的正负号一样。此外，如果σ很小，则略过本次的叠代。

2、取得向量其中常数

3、建立向量q＝p-Ku，其中常数

4、有了向量q，可进行矩阵化简A'＝A-q·u^T-u·q^T。

完成上述步骤之后，就可以将主对角三个对角线以外的元素变成零而得到tridiagonal matrix，而能够对得到的tridiagonal matrix求取特征向量，取得的特征向量再加以反转换，也就是乘上P₁·P₂·P₃…·P_n-2即可得矩阵A的特征向量。以本发明而言，上述的矩阵A以E(M_AB)代入执行。

由于E(M_AB)为4×4矩阵，所以特征向量长度为4，表示为[q1q2q3q4]，可转换表示为旋转矩阵R：

于取得旋转矩阵R后，利用第一影像31与旋转后的第二影像32计算一第二位移量T，是以，该第二影像32的新坐标Bs_i'为Bs_i'＝RBs_i+T，并评估邻近Asi与Bs_i'的误差是否够落入一由使用者设定的容许范围，若误差落在容许范围，表示该第一影像31与第二影像32可完整对齐接合，如图8所示，牙缝处42(由较密集的点所构成的区域代表牙缝)已经没有错位的情形，牙缝处42为平滑的线条，则该影像处理装置20存储第一影像31与第二影像32的最佳接合位置的计算结果。若邻近Asi与Bs_i'的误差落超出容许范围，回复执行第106步骤以再次计算旋转矩阵R与第二位移量T，直到误差落到容许范围。或者，若重复执行第106步骤的次数达到一门槛值，则以产生最小的误差的该次第二位移量T为最后的第二位移量T。

若欲完成半颚牙齿影像，则可于第一影像31与第二影像32完成接合后，可根据前述步骤进行第二影像32与一连续的第三影像的进行接合，并依此类推，直到完成单颚所有牙齿的接合，得到一精确的半颚牙齿影像。

Claims (6)

1.一种牙体影像的接合方法，其特征在于，所述牙体影像的接合方法包含：

由一探头取得一第一影像与一第二影像，所述第一影像与第二影像为点云数据；

由一影像处理装置接收并存储所述影像，以计算该第一影像与第二影像之间的重复区域以及该第二影像的一第一位移量，以利用该第一位移量移动该第二影像；

获取该第二影像中的齿面特征点以及于非齿面特征区域的样本点Bsi；

计算该第二影像的齿面特征点与样本点Bsi在该第一影像中的对应点Asi；

根据Asi与Bsi建立一旋转矩阵与该第二影像的一第二位移量，以根据该旋转矩阵与第二位移量移动该第二影像与第一影像接合。

2.如权利要求1所述牙体影像的接合方法，其特征在于，取得该第一位移量的步骤包含：

根据深度图分别计算该第一影像与第二影像的梯度大小；

根据该第一影像与第二影像的梯度大小计算该第一影像与第二影像的梯度累计值；

根据该第一影像与第二影像的梯度累计值计算该第一位移量。

3.如权利要求2所述牙体影像的接合方法，其特征在于，在取得齿面特征点的步骤中，齿面特征点为齿面的三维几何特征点，根据几何特征S(p)取得所述齿面特征点：

p为该第二影像中的任一点，S(p)介于0与1之间，k₁(p)与k₂(p)分别为该p点的两个主曲率，且k₁(p)≥k₂(p)；齿面特征点指S(p)为1～0.9、0～0.1或0.45～0.55的点。

4.如权利要求3所述牙体影像的接合方法，其特征在于，在取得样本点的步骤中，利用均匀取样手段对该第二影像的非齿面特征区域获取样本点。

5.如权利要求1至4中任一项所述牙体影像的接合方法，其特征在于，产生该旋转矩阵的步骤包含：

产生一矩阵M_AB，其中 为该第一影像特征点坐标与样本点坐标的平均值，为该第二影像齿面特征点坐标与样本点坐标的平均值；

产生4×4矩阵E(M_AB)，

对E(M_AB)通过Householder与QL演算法进行特征值与特征向量的计算以得出该旋转矩阵。

6.如权利要求5所述牙体影像的接合方法，其特征在于，该第二位移量为其中R为该旋转矩阵，该第二影像的特征点或样本点与该第一影像接合后的位置表示为Bs_i'＝RBs_i+T。