CN1928921A - 三维扫描系统中特征点云带的自动搜索方法 - Google Patents

Abstract

三维扫描系统中特征点云带的自动搜索方法主要是涉及到逆向工程中，用三维扫描系统对产品建立点云模型过程中离散点特征点云带自动计算的一种方法。属于图像三维信息重构的技术领域。该方法是先对点云进行空间划分，将整个点云划分成一个个的立方块，每个点分属于且仅属于一个立方块，在此基础上求取每个点的邻近点，然后由这些邻近点用最小二乘法求取此点法向量，再由法向量求得此点曲率，由曲率的大小就可以分割整个点云数据，从而得到特征点云带，该方法既可靠又避免了必须判断法向量方向的缺点，仿真表明该算法效果较好。

Description

三维扫描系统中特征点云带的自动搜索方法

技术领域

本发明主要是涉及到逆向工程中，用三维扫描系统对产品建立点云模型过程中离散点特征点云带自动计算的一种方法。属于图像三维信息重构的技术领域。

背景技术

逆向工程中，很重要的一个研究方向是寻找物体的特征区域。从海量散乱点云数据中，寻找物体的特征区域，不仅是做特征曲线的直接前提，而且在滤除杂点，拼合点云，由点生成曲面(根据特征区域，将整个自由曲面(整个点云)分割成不同组件(各个点云)，对各个点云生成多个小片曲面，再根据需要拼合)，都有十分重要的意义。

目前，虽然许多学者对点云数据处理做了大量研究，但对如何由海量点云数据得到特征区域这一方向，研究得较少。一些文献提出的利用相邻点的法向量差值的大小为原则来判断原曲面在该处的凹凸程度，这样即使对简单曲面也需要判断其法向量的正反方向。这种方法一方面耗时，另一方面当曲面凹凸变化剧烈时，法向量的正反方向并不存在较好的法则。另一些文献提出了将曲面点云按照特征区域划分成不同多个组件的算法，但是这些算法为了寻找特征区域要求先对点云三角网格化，而三角网格化是十分耗时的。还有一些文献提出的方法是建立在曲面方程上的，这在只给出原始数据而不附加任何其它信息的条件下是很难做到的。

发明内容

技术问题：本发明的目的是提供一种三维扫描系统中特征点云带的自动搜索方法，该方法既可靠又避免了必须判断法向量方向的缺点，仿真表明该算法效果较好。

技术方案：在对产品进行数字化的过程中，基于海量散乱点云的特征区域搜索在逆向工程中有着重要应用，本发明给出了一种以物体表面上不附加任何几何信息(包括测点法向量、曲面边界信息)的海量数据为处理对象，自动搜索物体表面特征区域的方法。

对海量点云的预处理：工程实际中的点云一般是数以十万计，甚至上百万。为了计算每个点的曲率，我们必须找出每个点的邻域，但在点云中，点的存储是无序的，所以理论上必须遍历整个点云，显然这是非常耗时的。本发明的做法是先对点云数据进行空间划分，将整个点云划分成一个个的立方块(每个点分属于且仅属于一个立方块)，在此基础上求取每个点的邻近点，然后由这些邻近点用最小二乘法求取此点法向量，再由法向量根据本发明提供的算法求得此点曲率。我们知道，在特征区域，由于形状变化剧烈，曲率比较大。因此由曲率的大小就可以分割整个点云数据，从而得到特征点云带。

本发明的三维扫描系统中特征点云带的自动搜索方法是先对点云进行空间划分，将整个点云划分成一个个的立方块，每个点分属于且仅属于一个立方块，在此基础上求取每个点的邻近点，然后由这些邻近点用最小二乘法求取此点法向量，再由法向量求得此点曲率，由曲率的大小就可以分割整个点云数据，从而得到特征点云带，具体方法为：

1)在滤除了杂点后，求得了相邻点的规范距离γ：即通常情况下，求得相邻点的最近距离，然后设定25倍的γ为小立方体的边长，对点云中的每一个点而言，它一定属于并且仅属于一个小立方块；当考虑到边界情况，这个点邻近小立方块的顶点或者面时，该点的邻域就不仅包含这个小立方块，而且包含该小立方块周围的多个小立方块，这样的小立方块共有27个，在这27个小立方块中，搜索出距离该点最近的若干个点，即构成了该点的邻域点，并将该点的所有邻域点依次存储在一个生成的临时的数据结构CCloseDot链表中；

2)求出法向量：根据各个点的邻域点，近似求出各点的切平面，根据该点的邻域点到待求平面的距离为最小的优化原则，用最小二乘法线性拟合出该点所在的切平面；

3)求取空间点的曲率：构造曲面M在点p处的离散数据表示，点p及其n个邻近点组成的邻域κ(n)，不含点p本身，π(p)为这些点的最小二乘拟合平面，将其作为构造曲面M在点p处的近似切平面， p是κ(n)的形心，称为点p的中心点， $\stackrel{\‾}{p}=\frac{1}{(n+1)}(\underset{p_i\∈\mathrm{\κ}\left(n\right)}{\mathrm{\Σ}}{p}_i+p);$ 各个邻近点到π(p)的距离为d_i，各个邻近点到点p的距离为τ_i。对点p存在函数f_i(p)，对点p的每一个邻近点可以唯一的确定一个值：f_i(p)＝d_i/τ_i则称 $f\left(p\right)=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{f}_i\left(p\right)$ 为点p的曲率函数，曲面M上的每一点的曲率函数就构成了曲面M的曲率场函数，这样可以用点p的曲率场函数近似表示曲面M在p点的凹凸程度，根据场函数值的大小来确定点云的特征区域，指定一个阀值，则场函数值大于阀值所对应的点位于特征区域，以不同颜色区分出来，

该方法的具体步骤如下：

步骤1：算法初始化，令γ＝相邻点的规范距离，

步骤2：从存储邻域点链表CCloseDot取出一个点，初始化整数num(邻域点链表中的邻近点数)为零，初始化f(p)为零，

步骤3：从该点的CNode(在这27个邻近小立方块中，遍历每一个点(该点本身除外)，将与该点最近的30个点用数据结构中的链表记录于链表CNode)中获得它的一个邻近点及其到该点的距离d_i，若d_i大于5倍的规范距离γ，计算f_i(p)＝d_i/τ_i，f(p)增加f_i(p)，整数num加1。若遍历完该点的领域点，转步骤4，否则，重复执行这一步，

步骤4：f(p)＝f(p)/num，将f(p)作为该点的曲率值记录该点中，转步骤2，

步骤5.计算得到特征点云带：指定一个阀值，当一点的由式 $f\left(p\right)=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{f}_i\left(p\right)$ 计算的曲率函数值大于该阀值时，就认为这个点所在区域为特征区域，而阀值的取值是各点曲率平均值的3倍左右。

有益效果：在对产品进行数字化的过程中，基于海量散乱点云的特征区域搜索在逆向工程中有着重要应用，本发明给出了一种以物体表面上不附加任何几何信息(包括测点法向量、曲面边界信息)的海量数据为处理对象，自动搜索物体表面特征区域的新的算法。具体的优点如下：

(1)避免了必须判断法向量方向的缺点，不需要判断其法向量的正反方向。

(2)本发明对海量散乱点云进行分块处理，减小了算法的时间复杂度。

(3)本发明提出的方法不是建立在曲面方程上的。

(4)仿真表明该算法效果较好。

(5)本发明以PC上的Windows NT为开发平台，以Visual C++6.0为开发环境，利用集成在Visual C++中的三维图形软件接口的新一代标准OpenGL进行编程。用Visual C++编的程序相对于其它编程工具有运行效率高的优势，而且功能强大，而三维图形处理运算量大，对效率要求高，因此选用Visual C++作为开发平台。而OpenGL始终贯彻了独立于软件平台的精神，它可以在Windows系列、Unix、Linux等不同系统上实现与Visual C++的完美结合，使得OpenGL编程得心应手。

附图说明

图1是逆向工程通用流程框图。

图2是光栅式三维扫描系统组成示意图。

图3是邻域点的数学模型示意图。

具体实施方式

下面结合附图示例对本发明的具体实施方式作进一步描述。本发明以PC上的Windows NT为开发平台，以Visual C++6.0为开发环境，利用集成在VisualC++中的三维图形软件接口的新一代标准OpenGL进行编程。主要包括以下步骤：

1)在滤除了杂点后，求得了相邻点的规范距离γ(即通常情况下，相邻点的最近距离)，然后设定25倍的γ为小立方体的边长。这样，当扫描仪的扫描精度发生变化(这时规范距离可以相差几倍)时，该算法依然适用，即使将来扫描精度大大提高了，该算法也不需作改变。在划分了空间区域后，邻域点的求取效率将得到极大的提高。对点云中的每一个点而言，它一定属于并且仅属于一个小立方块。但是当考虑到边界情况，比如这个点邻近小立方块的顶点或者面时，该点的邻域就不仅包含这个小立方块，而且包含该小立方块周围的多个小立方块，这样的小立方块共有27个。在这27个小立方块中，搜索出距离该点最近的若干个点，即构成了该点的邻域点。分块算法的主要步骤：

步骤1：算法初始化，令

Point Dis tan ce＝预先指定的相邻点的规范距离

InternalX，InternalY，InternalZ＝预先指定的点云在x轴，y轴，z轴上的最大间距

步骤2：计算NumX＝(int)InternalX/(25×Point Dis tan ce)+1，

           NumY＝(int)InternalY/(25×Point Dis tan ce)+1，

           NumZ＝(int)InternalZ/(25×Point Dis tan ce)+1，分别得到小立方块在x轴，y轴，z轴上的个数

步骤3：从点云中取出一个点的地址，构成CNode数据结构，记录于CCloseDot中。计算IndexX＝(int)InternalX/(25×Po int Dis tan ce)，

IndexY＝(int)InternalY/(25×Po int Dis tan ce)，

IndexZ＝(int)InternalZ/(25×Point Dis tan ce)，就可得到该点所在小立方块在x轴，y轴，z轴上的索引m＝IndexZ×NumX×NumY+IndexY×NumX+IndexX，考虑到其它26个小立方块，用一个3×3×3的三维数组array[i][j][k]＝IndexZ×(i-1)×NumX×NumY+IndexY×(j-1)×NumX+IndexX×(k-1)+m，(此处i，j，k＝0，1，2)，应注意的是对x轴，y轴，z轴上的边界情况要特殊考虑

步骤4：在这27个邻近小立方块中，遍历每一个点(该点本身除外)，将与该点最近的30个点用数据结构中的链表记录于CNode中，并将这些点距离该点的距离记录在该点中。遍历完这27个小立方块，转步骤3。

2)法向量是否精确求出关系到特征区域能否准确求出。根据各个点的邻域点，近似求出该点的切平面，当然越准确越好。这里作者根据一定的优化原则(该点的领域点到待求平面的距离为最小)，用最小二乘法线性拟合出该切平面。过程如下：

在测量点所在空间中，给定N个点

$P_i=(x_1^{\left(i\right)},x_2^{\left(i\right)},x_3^i),$ i＝1，2，...，N                           (1)

其中，i表示欲求的第i组领域点的法式，x₁ ⁽ⁱ⁾，x₂ ⁽ⁱ⁾，x₃ ⁱ表示第i组点的3维坐标，欲求一平面H：K^Tx＝k₀，(P₁为待求法向量的点，其余P_i为P₁的邻近点)。其中K＝(k₁，k₂，k₃)^T，表示欲求平面的法矢。x＝(x₁，x₂，x₃)^T，使得各点P_i到H的距离平方和为最小。

这样求得的H为点组{P₁，P₂，...，P_N}的最小二乘线性拟合。P_i到H的距离的平方为： $d_i^2=\frac{{(K^T{P}_i-k_0)}^2}{K^{T}K},$ i＝1，2，3，...，N。

从而诸点P_i到平面H的距离平方和

$\mathrm{\σ}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{d}_i^2=\frac{1}{K^{T}K}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{(K^T{P}_i-k_0)}^2---\left(2\right)$

很显然，平面H

K^Tx＝k₁x₁+k₂x₂+k₃x₃＝k₀                            (3)

方程的系数k₁，k₂，k₃构成的向量K为H的法向量。欲使距离平方和达到极小，

式(2)的一个必要条件是： $\frac{\∂\mathrm{\σ}}{\∂k_0}=-\frac{2}{K^{T}K}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(K^T{P}_i-k_0)=0$

从而得到： $k_0=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{K}^T{P}_i=K^T\left(\frac{1}{K}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{P}_i\right),$ 令

$\stackrel{\‾}{P}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{P}_i---\left(4\right)$

它为给定点的几何中心。从而有平面H

K^T(x- P)＝0                                          (5)

即最佳超平面H一定过中心点 P，由式(2)、(4)得到

$\mathrm{\σ}=\frac{1}{K^{T}K}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{(K^T{P}_i-k_0)}^2=\frac{1}{K^{T}K}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\left[(K^T{P}_i-\stackrel{\‾}{P})\right]}^2---\left(6\right)$

若令Y_i＝P_i- P∈R^m， $\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{Y}_i{Y}_i^T=A\∈R_m^m,$ 将矩阵A称为协变矩阵。根据 $K^T{Y}_i=Y_i^{T}K,$ 再联立式(6)，可得

$\mathrm{\σ}=\frac{1}{K^{T}K}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\left[K^T{Y}_i\right]}^2=\frac{1}{K^{T}K}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{K}^T{Y}_i{Y}_i^{T}K=\frac{1}{K^{T}K}{K}^T\mathrm{AK}---\left(7\right)$

又令Z＝K/‖K‖( $\left|\right|K\left|\right|=\sqrt{K^{T}K},$ ‖Z‖＝1)，则σ＝Z^TAZ。引入Lagrange乘子μ，即求函数φ(Z，μ)＝Z^TAZ+μ(Z^TZ-1)的极小点，该函数分别对Z、μ求导得

$\frac{\∂\mathrm{\φ}}{\∂Z}=2\mathrm{AZ}+2\mathrm{\μZ}=0---\left(8\right)$

$\frac{\∂\mathrm{\φ}}{\∂\mathrm{\μ}}=Z^{T}Z-1=0---\left(9\right)$

由式(8)、(9)得到：AZ＝-μZ，‖Z‖＝1，令μ＝-λ，得到

       AZ＝λZ                                           (10)

则λ为矩阵A的特征值，所以得：σ＝-μZ^TZ＝-μ‖Z‖²＝-μ＝λ

σ＝λ                                                   (11)

由式(11)易知σ的最小值为矩阵A的最小特征值，而矩阵A的最小特征值对应的单位特征向量Z即为待求平面H的法向量K的单位化结果，现在求平面H的法向量的问题就转化为求矩阵A的最小特征值对应的单位特征向量。对该特征向量的求取有两种方法：或者运用反乘幂法直接求取矩阵A的最小特征值对应的单位特征向量，或者求出矩阵A的逆矩阵，再运用乘幂法直接求取该逆矩阵的最大特征值对应的单位特征向量，由于乘幂法过程简单、可靠性高，故这里采用后一种方法。若矩阵A不可逆，则认为该点位于特征区域，由于发生这种情况的概率非常小(实用中这种情况大约为万分之一)，故它不会对结果有实质影响。

求法向量算法的主要步骤：

Step 1.在链表CCloseDot中取出一个点，从它的链表CNode中获得它的30个邻近点，生成这31个点的中心点，由此生成协变矩阵Matrix。

Step 2.求出该矩阵的行列式，若该行列式小于0.01，该点就属于特征区域，转Step 1，否则求其逆矩阵ReverseMatrix，用乘幂法求出ReverseMatrix的相应于最大特征值的特征向量。

Step 3.将Step 2求得的特征向量单位化，将该单位向量记录于该点中，转Step 1。

3)求取空间点的曲率

构造曲面M在点p处的离散数据表示如图3所示：

设点p及其n个邻近点组成的邻域κ(n)(不含点p本身)，π(p)为这些点的最小二乘拟合平面，将其作为构造曲面M在点p处的近似切平面，p是κ(n)的形心，称为点p的中心点

$\stackrel{\‾}{p}=\frac{1}{(n+1)}(\underset{p_i\∈\mathrm{\κ}\left(n\right)}{\mathrm{\Σ}}{p}_i+p)---\left(12\right)$

各个邻近点到π(p)的距离为d_i，各个邻近点到点p的距离为τ_i。相关定义：(假设已对点云进行了预处理)对点p存在函数f_i(p)，对点p的每一个邻近点可以唯一的确定一个值：f_i(p)＝d_i/τ_i(对不同邻近点，f_i(p)的值是不同的)，则称

$f\left(p\right)=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{f}_i\left(p\right)---\left(13\right)$

为点p的曲率函数，曲面M上的每一点的曲率函数就构成了曲面M的曲率场函数。这样可以用点p的曲率场函数近似表示曲面M在p点的凸凹程度，根据场函数值的大小来确定点云的特征区域。指定一个阀值，则场函数值大于阀值所对应的点位于特征区域，应以不同颜色区分出来。

假设κ(n)真实地包含了待重建曲面M在p点附近的几何形状信息，则M在p点处的曲率越大，p点的邻近点到p点的距离就会越大，其在点的曲率场函数就越大，由此可得出结论：待重建曲面M的曲率场函数表达了曲面M在该点的凸凹程度。从场函数公式(13)分析，邻近点的个数n越小越能逼近M在该点该点处的真实曲率，而邻近点的个数n越大越能全面的表达曲面M在p点附近的几何形状信息，从以上分析可以看出，邻近点的个数n的选取是比较重要的，并且，n值的选取应保证曲面M在κ(n))范围内是单凸或者单凹的，而不能是振荡的；若n值过小，当数据点在各向分布不太均匀时，会使π(p)不能代表曲面M在点p处的切平面；n过大，则邻近点生成的f_i(P)＝d_i/τ_i可能会由于d_i大而τ_i小，而导致场函数偏大。因而，n的取值要根据曲面M的凹凸程度和测点的分布情况而定，实践中一般将n取为20-30，效果较好。

在由式(4)确定了测点的中心点 P和式(10)确定了测点所在切平面π(p)的法向量后，这个最佳拟合平面H也就确定了。运用式(13)计算并求取每一测点的曲率函数值即可。值得关注的是，尽管我们已经对海量点云进行了预处理并剔除了绝大部分的杂点(非常接近甚至相重叠的点及背景点)，但杂点或不期望的点总是存在的。为此，在计算邻近点的函数值f_i(p)＝d_i/τ_i时，若d_i过大则将该点舍弃。由于该点距离点p很远，故它并不能反映点p处的几何信息，测试表明舍弃该f_i(p)后对整体效果的影响可以忽略。实践中，当d_i大于5-7倍的γ(γ为相邻点的规范距离)时，可以舍弃该邻近点。

算法的主要步骤：

Step 1.算法初始化，令γ＝相邻点的规范距离

Step 2.从链表CCloseDot中取取出一个点，初始化整数num为零，初始化f(p)为零

Step 3.从该点的链表CNode中获得它的一个邻近点及其到该点的距离d_i，若d_i大于5倍的规范距离γ，计算f_i(p)＝d_i/τ_i，f(p)增加f_i(p)，整数num加1。若遍历完该点的领域点，转Step 4，否则，重复执行这一步

Step 4.f(p)＝f(p)/num，将f(p)作为该点的曲率值记录该点中，转Step 2。

4)计算得到特征点云带

指定一个阀值，当一点的由式(13)计算的曲率函数值大于该阀值时，就认为这个点所在区域为特征区域，而阀值的取值是各点曲率平均值的3倍左右。根据本算法求得的场函数值不会过大或过小，数据相对比较整齐。应用中为节时，可以抽取若干个点的场函数值求它们的平均值，这样不仅省时，而且当点的个数适当时，并不影响精度。例如，当抽取的点数达到100左右时，效果已经较好。

对海量测点集进行数据简化，对于后续的模型重建有重要意义。本发明首先提出了对海量数据进行空间划分的算法，提高了算法的处理效率。建立了测点集的邻近域，并且不需要判断法式的正反方向，最后，给出了相应的算法与优化准则。应用本发明算法，我们对一系列测量数据寻找它的特征区域。应用实例表明本发明提出的算法具有较高的效率，达到了预期的效果。

Claims (2)

1.一种三维扫描系统中特征点云带的自动搜索方法，其特征在于该方法是先对点云进行空间划分，将整个点云划分成一个个的立方块，每个点分属于且仅属于一个立方块，在此基础上求取每个点的邻近点，然后由这些邻近点用最小二乘法求取此点法向量，再由法向量求得此点曲率，由曲率的大小就可以分割整个点云数据，从而得到特征点云带，具体方法为：

1.)在滤除了杂点后，求得了相邻点的规范距离γ：即通常情况下，求得相邻点的最近距离，然后设定25倍的γ为小立方体的边长，对点云中的每一个点而言，它一定属于并且仅属于一个小立方块；当考虑到边界情况，这个点邻近小立方块的顶点或者面时，该点的邻域就不仅包含这个小立方块，而且包含该小立方块周围的多个小立方块，这样的小立方块共有27个，在这27个小立方块中，搜索出距离该点最近的若干个点，即构成了该点的邻域点，并将该点的邻域点依次存储在一个生成的临时的数据结构CCloseDot链表中；

2.)求出法向量：根据各个点的邻域点，近似求出各点的切平面，根据该点的领域点到待求平面的距离为最小的优化原则，用最小二乘法线性拟合出该点所在的切平面；

3)求取空间点的曲率：构造曲面M在点p处的离散数据表示，点p及其n个邻近点组成的邻域κ(n)，不含点p本身，π(p)为这些点的最小二乘拟合平面，将其作为构造曲面M在点p处的近似切平面， p是κ(n)的形心，称为点p的中心点， $\stackrel{\‾}{p}=\frac{1}{(n+1)}(\underset{p_i\∈\mathrm{\κ}\left(n\right)}{\mathrm{\Σ}}{p}_i+p);$ 对点p存在函数f_i(p)，对点p的每一个邻近点可以唯一的确定一个值：f_i(p)＝d_i/τ_i则称 $f\left(p\right)=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{f}_i\left(p\right)$ 为点p的曲率函数，曲面M上的每一点的曲率函数就构成了曲面M的曲率场函数，这样可以用点p的曲率场函数近似表示曲面M在p点的凸凹程度，根据场函数值的大小来确定点云的特征区域，指定一个阀值，则场函数值大于阀值所对应的点位于特征区域，以不同颜色区分出来。

2.根据权利要求1所述的三维扫描系统中特征点云带的自动搜索方法，其特征在于该方法的具体步骤如下：

步骤1：算法初始化，令γ＝相邻点的规范距离，

步骤2：从存储邻域点链表CCloseDot取出一个点，初始化整数num(邻域点链表中的邻近点数)为零，初始化f(p)为零，

步骤3：从该点的链表CNode中获得它的一个邻近点及其到该点的距离d_i，若d_i大于5倍的规范距离γ，计算f_i(p)＝d_i/τ_i，f(p)增加f_i(p)，整数num加1。若遍历完该点的领域点，转步骤4，否则，重复执行这一步，

步骤4：f(p)＝f(p)/num，将f(p)作为该点的曲率值记录该点中，转步骤2，

步骤5.计算得到特征点云带：指定一个阀值，当一点的由式 $f\left(p\right)=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{f}_i\left(p\right)$ 计算的曲率函数值大于该阀值时，就认为这个点所在区域为特征区域，而阀值的取值是各点曲率平均值的3倍左右。

Images