CN110599506B - 一种复杂异形曲面机器人三维测量的点云分割方法 - Google Patents

一种复杂异形曲面机器人三维测量的点云分割方法 Download PDF

Info

Publication number
CN110599506B
CN110599506B CN201910908581.1A CN201910908581A CN110599506B CN 110599506 B CN110599506 B CN 110599506B CN 201910908581 A CN201910908581 A CN 201910908581A CN 110599506 B CN110599506 B CN 110599506B
Authority
CN
China
Prior art keywords
point
points
point cloud
curved surface
curvature
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201910908581.1A
Other languages
English (en)
Other versions
CN110599506A (zh
Inventor
王耀南
史雅兰
毛建旭
朱青
吴昊天
彭伟星
苏学叁
田吉委
黄亨斌
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hunan University
Original Assignee
Hunan University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Hunan University filed Critical Hunan University
Priority to CN201910908581.1A priority Critical patent/CN110599506B/zh
Publication of CN110599506A publication Critical patent/CN110599506A/zh
Application granted granted Critical
Publication of CN110599506B publication Critical patent/CN110599506B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/23—Clustering techniques
    • G06F18/232—Non-hierarchical techniques
    • G06F18/2321—Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions
    • G06F18/23211—Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions with adaptive number of clusters
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/11—Region-based segmentation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/136—Segmentation; Edge detection involving thresholding
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10028—Range image; Depth image; 3D point clouds
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20024—Filtering details

Abstract

本发明公开了一种复杂异形曲面机器人三维测量的点云分割方法，包括以下步骤：S100、输入以地面、桌面为背景的叶片点云X，利用体素滤波滤除背景点，得到目标叶片点云Y；S200、利用PCA算法计算出Y中点的法向量、平面轮廓度，并剔除离群点，将关联点集合记为一致集CS；S300、利用法向量、平面轮廓度偏差建立成对连接，确定聚类中心后进行搜索，并搜索所有与它相连的点，生成聚类C；S400、利用Delaunay三角剖分方法对聚类C进行曲面拟合；S500、对于每一个拟合的曲面切片，计算其曲率，并设定曲率偏差阈值，若相邻两个曲面切片之间的曲率偏差小于阈值，则合并；否则不合并，从而得到与背景点云分割开的完整叶片点云Y。具有分割准确、输入参数少、鲁棒性强的优点。

Description

一种复杂异形曲面机器人三维测量的点云分割方法
技术领域
本发明涉及视觉测量技术领域，特别是涉及一种复杂异形曲面机器人三维测量的点云分割方法。
背景技术
复杂异形曲面是一类形状非结构化、曲面复合弯曲并在空间任意堆叠、曲率变化大且边缘检测较难的曲面。以航空发动机涡轮叶片为例，复杂异形曲面在航空航天、汽车制造等高端智能制造领域应用十分普遍，人们对于如何获得其三维信息并进行精密测量的需求日益提高。然而，三维视觉测量获得的复杂异形曲面点云是非结构化的，且噪声点多。点云分割作为对三维点云进行加工时的重要预处理过程，是三维测量的关键步骤，点云分割的效果好坏与否直接影响着后续的三维测量乃至最后的加工、打磨的精确度。
对于复杂异形曲面的点云分割，主要存在的难点有难分割、误分割、分割慢。首先，由于复杂异形曲面结构具有形状非结构化、曲面复合弯曲并在空间任意堆叠等特殊性，因此，难以找到特征点的同时也难以对其进行拟合，从而造成分割难；其次，由于其曲率变化大且边缘检测较难，这会导致背景点以及边缘噪声点被误分割入前景点云，从而无法达到精密测量所需的精度；最后，复杂异形曲面点云以及背景点云数据量大，处理过程慢，不适用于效率要求高的工业制造领域。以上的难点会使得分割精度与效率降低，影响后续测量。
针对这一问题，目前采取的大多都是RANSAC(随机采样一致)算法、区域生长算法等传统分割算法。RANSAC算法通过创建一个分割器，设置目标几何形状，设置误差容忍范围，输入并分割点云。但其缺点在于，随着局外点的增多，迭代次数会变多，耗时也会变长。如果设置了迭代上限，可能会导致抽样不充分从而计算出错误的模型；如果不控制迭代次数，便会造成算法效率低下难以收敛。此外，只能在特定的数据集中时估计出一个模型，对于数据集中存在两个或者两个以上的模型，RANSAC不能识别。区域生长算法首先人工或随机选择种子点，将种子周围的点和种子相比，如果同时满足法线方向足够近、曲率是否足够小，则该点可用做种子。如果只满足法线方向足够近，则归类而不做种子。从某个种子出发，其“子种子”不再出现，则一类聚集完成。该方法简单易实现，但是太依赖于种子点的选择或者区域生长的规则，通常只能用于初步的语义分割，精细化处理阶段需与其他方法混合使用，导致其精度不高。
发明内容
有鉴于此，本发明提出一种复杂异形曲面机器人三维测量的点云分割方法，在体素滤波与成对连接聚类算法的基础上，筛除了背景点云及误分割点，从而提高了分割的精度和效率，大大节省了所需的迭代时间及运算量。
一方面，本发明提供了一种复杂异形曲面机器人三维测量的点云分割方法，包括以下步骤：
S100、输入以地面、桌面为背景的叶片点云X，利用体素滤波滤除背景点，得到目标叶片点云Y；
S200、利用PCA算法计算出Y中点的法向量、平面轮廓度，并剔除离群点，将关联点集合记为一致集CS；
S300、利用法向量、平面轮廓度偏差建立成对连接，确定聚类中心后进行搜索，并搜索所有与它相连的点，生成聚类C；
S400、利用Delaunay三角剖分方法对步骤S300产生的聚类C进行曲面拟合；
S500、对于每一个拟合的曲面切片，计算其曲率，并设定曲率偏差阈值，若相邻两个曲面切片之间的曲率偏差小于阈值，则合并；否则不合并，从而得到与背景点云分割开的完整叶片点云Y。
进一步地，所述步骤S100具体包括如下步骤:
S101、将整个点云场景在XY平面上垂直划分为一系列局部点云块Bi，i＝1,2,...,n，n为划分后生成的局部点云块的总数；
S102、根据Octree索引，对于每个局部点云块Bi进行划分，按一定的宽度r将其划分为空间上一系列连续的体素Vj，j＝1,2,...,m，m为划分后生成的体素的总数；
S103、将与体素Vj相邻接并在空间上位于Vj上方的9个邻域体素定义为向上邻域(L1,L2,...,L9)，并随机规定一个的向上邻域为向上生长起点,并按同样的方式，分别沿着Vj的向上邻域继续向上生长，直到所有已执行向上生长操作的体素无向上邻域可生长，向上生长过程结束；
S104、对于已执行向上生长操作的体素，计算出其整体高度值、局部高度值，并选出具有最大局部高度值的体素Vh，并将其定为该生长区域的最高点；
S105、计算地面高度，定义Hg为地面高度的临界值：如果Vh与Hg的高度差小于第一阈值hv1，则Vh被标记为地面体素，同时滤除包含Vh的局部点云块中的所有点；如果Vh与Hg的高度差大于第一阈值hv1，则Vh被标记为非地面体素，同时保留包含Vh的局部点云块中的所有点，为非地面点；
S106、计算桌面高度，定义Hd为桌面高度的临界值：如果Vh与Hd的高度差小于第二阈值hv2，则Vh被标记为桌面体素，同时滤除包含Vh的局部点云块中的所有点；如果Vh与Hd的高度差大于第二阈值hv1，则Vh被标记为非桌面体素，同时保留包含Vh的局部点云块中的所有点，为非桌面点；
S107、非桌面点和非地面点的交集则为目标叶片点云Y。
进一步地，所述步骤S200具体包括如下步骤：
S201、在目标叶片点云Y中随机选取一个点，记为Pi，构建k-d树，搜索出K个与Pi距离最近的临近点，记为KNN，将KNN归为同一个邻域内，此邻域记为Ii；
S202、将Ii中的点按与Pi之间的距离升序排列，选出前
Figure BDA0002214029330000031
个点，利用式(1)计算出协方差矩阵Σ：
Figure BDA0002214029330000041
其中，
Figure BDA0002214029330000042
为平均向量，T表示向量的转置；
S203、利用奇异值分解(SVD)产生标准特征值方程：
λV＝∑V (2)
式中，λ为特征值矩阵，V为特征向量矩阵；
S204、取前三个最大的特征值λ0、λ1、λ2及其对应的特征向量v0、v1、v2，且λ2>λ1>λ0,令λ0表示Pi的平面轮廓度，v0表示Pi所在平面的法向量；
S205、利用如下公式查找Pi的关联点、离群点：
Figure BDA0002214029330000043
式中，Rz表示判断一个点是否为点Pi关联点的权值得分，
Figure BDA0002214029330000044
是Pi的KNN与估计平面的正交距离,NOD是/>
Figure BDA0002214029330000045
的集合，median(NOD)为NOD的中位数，其中：
Figure BDA0002214029330000046
式中，a为设定常数，取值为1.4826；
当Rz<2.5时，认为该点为Pi的关联点，将其归到一致集CS中；反之为离群点，将其剔除。
进一步地，所述步骤S300具体包括如下步骤：
S301、在Pi的CS中搜索，找出CS中λ比λ(Pi)小的邻域点，并从这些邻域中选择n与n(Pi)偏差最小的邻域点，作为CNP(Pi)，则CNP(Pi)为与点Pi法向量、平面轮廓度偏差最小的点，并建立Pi与CNP(Pi)的成对连接；
S302、若CS中已经不存在λ比λ(Pi)小的邻域点了，则如下公式计算阈值：
Figure BDA0002214029330000047
式中，thλ表示阈值，
Figure BDA0002214029330000048
是所有N个数据点平面轮廓度的平均值，其中：
Figure BDA0002214029330000049
其中，N表示数据点的数量，
Figure BDA0002214029330000051
是所有N个数据点平面轮廓度的平均值；
S303、如果λ(Pi)<thλ，则将Pi记为聚类中心，并归入Ccenter；否则，不记为聚类中心；
S304、找到Ccenter，搜索所有与其间接、直接建立了成对连接的点，并将这些点生成聚类C。
进一步地，所述步骤S400具体包括如下步骤：
S401、针对聚类C中的点，按x坐标从小到大进行排序，根据离散点的最大分布确定包含了点集中所有点的超级三角形；
S402、初始化边缓存数组edge buffer，并遍历temp triangles中的每一个三角形，计算该三角形的外接圆的圆心和半径；
S403、遍历超级三角形中的每一个数据点，若该点在外接圆的右侧，则该三角形为Delaunay三角形，保存至triangles中，并在temp triangles中去掉；若该点在外接圆的内侧，则该三角形不为Delaunay三角形，保存至edge buffer中，并在temp triangles中去掉；若该点在外接圆的左侧，则该三角形为不确定；
S404、对edge buffer数组进行去重，将此时edgebuffer中保留的边与当前的点相结合，形成若干三角形，保存至temp triangles中，并将triangles和temp triangles中的三角形合并，组成最终的超级三角形；
S405、除去与超级三角形三个点有关的所有三角形，对于聚类C进行的曲面拟合完成，将拟合好的曲面记为曲面切片Sp。
进一步地，步骤S400中Delaunay三角剖分方法中三角网格模型通过如下线性表表示：
M＝(V,F) (7)
式中，V＝{vi；1≤i≤nv}表示顶点集，F＝{fk；1≤k≤nk}表示三角片集。
进一步地，所述步骤S500中曲面切片Sp的曲率值通过如下步骤计算：
(1)通过如下公式计算曲面切片Sp中每个三角面片fk的法向量Nfk：
Figure BDA0002214029330000061
式中，vi和vj+1分别表示曲面切片Sp上任意选取的两点，
Figure BDA0002214029330000062
表示曲面切片Sp点vi到vj+1的向量，/>
Figure BDA0002214029330000063
表示曲面切片Sp点vj+1到vj的向量。
(4)通过如下公式计算任意点vi的法向量
Figure BDA0002214029330000064
Figure BDA0002214029330000065
式中，|Ni|表示除顶点vi外其它顶点组成集合的顶点个数，Fi为包含点vi三角面片fk的集合；
(5)根据如下公式计算点vi的法曲率Kij：
Figure BDA0002214029330000066
式中，
Figure BDA0002214029330000067
表示曲面切片Sp点vj到vi的向量；
(4)根据如下公式：
Figure BDA0002214029330000068
式中，K(SP)为高斯曲率，用于反映曲面总的弯曲程度，H(SP)为平均曲
率，用于反映空间曲面上某一点任意两个相互垂直的正交曲率的平均值，K1为最大法曲率的极大值，K2为垂直于极大曲率面曲率的极小值。
进一步地，所述最大法曲率通过如下方式选取：遍历三角拟合经过点vi的所有条曲线，计算并选出具有最大法曲率的曲线，此时的法曲率为最大法曲率。
进一步地，所述步骤S500中相邻两个曲面切片之间的曲率偏差包括高斯曲率偏差和平均曲率偏差，通过如下公式表达：
Figure BDA0002214029330000069
式中，Sp1、Sp2为两个相邻曲面切片，K(SP1)、K(SP2)分别为曲面切片Sp1、Sp2的高斯曲率，H(SP1)、H(SP2)分别为曲面切片Sp1、Sp2的平均曲率，ΔK(Sp)为高斯曲率偏差，ΔH(Sp)为平均曲率偏差。
进一步地，所述步骤S100中叶片点云X通过3D位移传感器扫描获取。
本发明提供的复杂异形曲面机器人三维测量的点云分割方法，相比于传统的点云分割过程，使用体素滤波直接滤除背景点云，有效减小输入数据规模和运算量，极大的提高了分割速度和鲁棒性；并采用改进的成对连接聚类算法，实现的精切快速的点云主目标分割，且利用了Delaunay三角剖分进行曲面拟合，该方法可以很好地运用于非结构化点云，对于曲率变化快、高次非常规函数的复杂异形曲面点云，不受限于分割物体的形状、结构，能达到较好的分割与曲面拟合效果；进一步地，本发明还通过曲率偏差的曲面合并方法，有效抑制可能的过度分割情况，为后续三维测量提供整体性良好的目标点云。
附图说明
构成本发明的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：
图1为本发明实施例提供的复杂异形曲面机器人三维测量的点云分割方法的流程图；
图2为超级三角形示意图；
图3-a至图3-d为超级三角形内点的遍历过程图；
图4为将triangles和temp triangles合并后得到的数组图；
图5为除去与超级三角形中点有关的三角形后的结果图；
图6为利用Delaunay三角剖分方法将聚类C拟合成的曲面切片Sp的过程图；
图7为切片Sp的示意图。
具体实施方式
需要说明的是，在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。下面将参考附图并结合实施例来详细说明本发明。
本发明主要针对的是复杂异形曲面，其曲面复合弯曲并在空间任意积叠，边缘曲率复杂。传统的分割方法会将背景点误分割为目标点，从而影响后续三维目标测量的精度，且直接进行分割会导致运算量大、分割速度慢，除此之外，传统分割方法在分割完成进行形状拟合时，只能拟合平面，或者是能用常规表达式拟合、曲率变化较小的曲面，而不能运用于非结构化点云，拟合效果较差。
为此，本发明提供了一种复杂异形曲面机器人三维测量的点云分割方法包括以下步骤：
S100、输入以地面、桌面为背景的叶片点云X，利用体素滤波滤除背景点，得到目标叶片点云Y；优选地，叶片点云X通过3D位移传感器扫描获取，3D位移传感器可采用扫描精度达0.01mm的Codex相机ds1101。
具体地，步骤S100中目标叶片点云Y通过如下步骤得到：
S101、将整个点云场景在XY平面上垂直划分为一系列局部点云块Bi，i＝1,2,...,n，n为划分后生成的局部点云块的总数；
S102、根据Octree(八叉树)索引，对于每个局部点云块Bi进行划分，按一定的宽度r将其划分为空间上一系列连续的体素Vj，j＝1,2,...,m，m为划分后生成的体素的总数，优选地，此处取r＝0.01mm；
S103、将与体素Vj相邻接并在空间上位于Vj上方的9个邻域体素定义为向上邻域(L1,L2,...,L9)，并随机规定一个的向上邻域为向上生长起点,并按同样的方式，分别沿着Vj的向上邻域继续向上生长，直到所有已执行向上生长操作的体素无向上邻域可生长，向上生长过程结束；
S104、对于已执行向上生长操作的体素，计算出其整体高度值、局部高度值，并选出具有最大局部高度值的体素Vh，并将其定为该生长区域的最高点；
S105、计算地面高度，定义Hg为地面高度的临界值：如果Vh与Hg的高度差小于第一阈值hv1，则Vh被标记为地面体素，同时滤除包含Vh的局部点云块中的所有点；如果Vh与Hg的高度差大于第一阈值hv1，则Vh被标记为非地面体素，同时保留包含Vh的局部点云块中的所有点，为非地面点；
S106、计算桌面高度，定义Hd为桌面高度的临界值：如果Vh与Hd的高度差小于第二阈值hv2，则Vh被标记为桌面体素，同时滤除包含Vh的局部点云块中的所有点；如果Vh与Hd的高度差大于第二阈值hv1，则Vh被标记为非桌面体素，同时保留包含Vh的局部点云块中的所有点，为非桌面点；
S107、非桌面点和非地面点的交集则为目标叶片点云Y。
S200、利用PCA(Principal ComponentAnalysis，主成分分析)算法计算出Y中点的法向量、平面轮廓度，并剔除离群点，将关联点集合记为一致集CS；
具体地，步骤S200可分解为如下步骤：
S201、在目标叶片点云Y中随机选取一个点，记为Pi，构建k-d树，搜索出K个与Pi距离最近的临近点，记为KNN(k-NearestNeighbor，邻近算法)，将KNN归为同一个邻域内，此邻域记为Ii；
S202、将Ii中的点按与Pi之间的距离升序排列，选出前
Figure BDA0002214029330000091
个点，利用式(1)计算出协方差矩阵Σ：
Figure BDA0002214029330000092
其中，
Figure BDA0002214029330000093
为平均向量，T表示向量的转置；
S203、利用奇异值分解(SVD)产生标准特征值方程：
λV＝∑V (2)
式中，λ为特征值矩阵，V为特征向量矩阵；
S204、取前三个最大的特征值λ0、λ1、λ2及其对应的特征向量v0、v1、v2，且λ2>λ1>λ0,令λ0表示Pi的平面轮廓度，v0表示Pi所在平面的法向量；
S205、利用如下公式查找Pi的关联点、离群点：
Figure BDA0002214029330000101
式中，Rz表示判断一个点是否为点Pi关联点的权值得分，
Figure BDA0002214029330000102
是Pi的KNN与估计平面的正交距离,NOD是/>
Figure BDA0002214029330000103
的集合，median(NOD)为NOD的中位数，其中：
Figure BDA0002214029330000104
式中，a为设定常数，取值为1.4826；
当Rz<2.5时，认为该点为Pi的关联点，将其归到一致集CS中；反之为离群点，将其剔除。
S300、利用法向量、平面轮廓度偏差建立成对连接，确定聚类中心后进行搜索，并搜索所有与它相连的点，生成聚类C。
优选地，步骤S300具体包括如下步骤：
S301、在Pi的CS中搜索，找出CS中λ比λ(Pi)小的邻域点，并从这些邻域中选择n与n(Pi)偏差最小的邻域点，作为CNP(Pi)，则CNP(Pi)为与点Pi法向量、平面轮廓度偏差最小的点，并建立Pi与CNP(Pi)的成对连接；
S302、若CS中已经不存在λ比λ(Pi)小的邻域点了，则如下公式计算阈值：
Figure BDA0002214029330000105
式中，thλ表示阈值，
Figure BDA0002214029330000106
是所有N个数据点平面轮廓度的平均值，其中：
Figure BDA0002214029330000107
其中，N表示数据点的数量，
Figure BDA0002214029330000108
是所有N个数据点平面轮廓度的平均值；
S303、如果λ(Pi)<thλ，则将Pi记为聚类中心，并归入Ccenter；否则，不记为聚类中心；
S304、找到Ccenter，搜索所有与其间接、直接建立了成对连接的点，并将这些点生成聚类C。
S400、利用Delaunay三角剖分方法对步骤S300产生的聚类C进行曲面拟合。
具体地，步骤S400包括如下步骤：
S401、针对聚类C中的点，按x坐标从小到大进行排序，根据离散点的最大分布确定包含了点集中所有点的超级三角形。具体方法是根据相似三角形定理，以图2三个点为例的超级三角形示意图为例，先求出连接三个点的矩形的一半，并划分成两个对角三角形，扩大一倍后，则扩大后的直角三角形斜边经过点(Xmax,Ymax)，但是为了将所有的点包含在超级三角形内，对三角形的右下角顶点进行底和高的扩展，并要保证这个扩展三角形底大于高，才能实现包含点集中所有点。此方法求解过程简单，运算量小，且求得的超级三角形大小适中。
S402、初始化边缓存数组edge buffer，并遍历temp triangles(未确定三角形列表)中的每一个三角形，计算该三角形的外接圆的圆心和半径。
S403、遍历超级三角形中的每一个数据点，若该点在外接圆的右侧，则该三角形为Delaunay三角形，保存至triangles(已确定的三角形列表)中，并在temp triangles中去掉；若该点在外接圆的内侧，则该三角形不为Delaunay三角形，保存至edge buffer中，并在temp triangles中去掉；若该点在外接圆的左侧，则该三角形为不确定。以图3-a至3-d为例，对遍历过程进行详细阐述：
①对temp triangles中的三角形遍历画外接圆，先对左边的第一个点进行判断，其在圆内，所以该三角形不为Delaunay三角形，将其三边保存至edge buffer中，temptriangles中删除该三角形；
②将该点与edge buffer中的每一个边相连，组成三个三角形，加入到temptriangles中,如图3-b所示；
③遍历temp triangles中的三个三角形1、2和3，并画外接圆，如图3-c所示，并使用按坐标排序的第二个点来进行判断，有以下三种情况：
a.该点在三角形1外接圆右侧，则表示左侧三角形为Delaunay三角形，将该三角形保存至triangles中；
b.该点在三角形2外接圆外侧，为不确定三角形，跳过，但并不在temp triangles中删除；
c.该点在三角形3外接圆内侧，则向清空后的edge buffer加入该三角形的三条边，并用该点写edge buffer中的三角边进行组合，组合成了三个三角形并加入到temptriangles中；
④将第二个点与edge buffer中的每一条边相连接，生成三个三角形。此时该数组里则含有四个三角形，如图3-d所示，一个是上次检查跳过的三角形和新根据第二个点生成的三个三角形，再次对temp triangles进行遍历，使用按坐标排序的第四个点判断，有以下四种情况：
a.该点在三角形1外接圆右侧，则该三角形为Delaunay三角形，保存至triangles中，并在temp triangles中删除；
b.该点在三角形2外接圆外侧，跳过；
c.该点在三角形3外接圆内侧，将该三边保存至temp buffer(保存三角形的临时数组)中，并在temp triangles中删除；
d.该点在三角形4外接圆内侧，将该三边保存至temp buffer中，并在temptriangles中删除。
这时，temp buffer中有六条边，triangles中有两个三角形，temp triangles中有1个三角形，对temp buffer中的六条边进行去重，得到五条边，将该点与这五条边组合成五个三角形并加入到temp triangles中，这时temp triangles中有6个三角形。三个点已经遍历结束，则不再对第三个点形成的三角形做外接圆，这时则将triangles与temp trianlges合并，合并后的数组表示包含已经确定的Delaunay三角形和剩下的三角形。图4即为将triangles和temp triangles合并后得到的数组。除去合并后数组中的和超级三角形三个点有关的所有三角形，即进行数组坐标的限定，则得到了最后的结果，具体如图5所示。
S404、对edge buffer数组进行去重，将此时edge buffer中保留的边与当前的点相结合，形成若干三角形，保存至temp triangles中，并将triangles和temp triangles中的三角形合并，组成最终的超级三角形；
S405、除去与超级三角形三个点有关的所有三角形，对于聚类C进行的曲面拟合完成，将拟合好的曲面记为曲面切片Sp。
进一步地，步骤S400中Delaunay三角剖分方法中三角网格模型通过如下线性表表示：
M＝(V,F) (7)
式中，V＝{vi；1≤i≤nv}表示顶点集，F＝{fk；1≤k≤nk}表示三角片集。
利用Delaunay三角剖分方法将聚类C拟合成的曲面切片Sp的过程图具体见图6。
S500、对于每一个拟合的曲面切片，计算其曲率，并设定曲率偏差阈值，若相邻两个曲面切片之间的曲率偏差小于阈值，则合并；否则不合并，从而得到与背景点云分割开的完整叶片点云Y。具体地该步骤通过以下过程展开：
(1)通过如下公式计算曲面切片Sp中每个三角面片fk的法向量Nfk：
Figure BDA0002214029330000131
式中，vi和vj+1分别表示曲面切片Sp上任意选取的两点，
Figure BDA0002214029330000132
表示曲面切片Sp点vi到vj+1的向量，/>
Figure BDA0002214029330000133
表示曲面切片Sp点vj+1到vj的向量。
(2)假设1-环邻域是与点vi相邻的三角形集合，对于离散三角网格模型M＝(V,F)，任意一点vi的法向量可定义为1-环三角形法向量的平均值，故此，可通过如下公式计算任意点vi的法向量
Figure BDA0002214029330000134
Figure BDA0002214029330000135
其中，在vi的1-环邻域中，除顶点vi外其它顶点组成的集合记为Vi。如果顶点vj属于Vi,则vj是vi的相邻点。Vi中顶点的个数称为其顶点的度，记为|Ni|。Fi为包含点vi三角面片fk的集合，vi为切片上任意选取的一点，fk为点vi的1-环邻域中的其中一个三角面片。
曲面切片Sp的示意图具体如图7所示。
(3)根据如下公式计算点vi的法曲率Kij：
Figure BDA0002214029330000141
式中，
Figure BDA0002214029330000142
表示曲面切片Sp点vj到vi的向量。
(4)根据如下公式计算曲面切片Sp的高斯曲率K(SP)和平均曲率H(SP)：
Figure BDA0002214029330000143
式中，K(SP)为高斯曲率，用于反映曲面总的弯曲程度，H(SP)为平均曲率，用于反映空间曲面上某一点任意两个相互垂直的正交曲率的平均值，K1为最大法曲率的极大值，K2为垂直于极大曲率面曲率的极小值。
此处需要说明的是，最大法曲率通过如下方式选取：遍历三角拟合经过点vi的所有条曲线，计算并选出具有最大法曲率的曲线，此时的法曲率即为最大法曲率。
在进一步的技术方案中，步骤S500中相邻两个曲面切片之间的曲率偏差包括高斯曲率偏差和平均曲率偏差，通过如下公式表达：
Figure BDA0002214029330000144
式中，Sp1、Sp2为两个相邻曲面切片，K(SP1)、K(SP2)分别为曲面切片Sp1、Sp2的高斯曲率(总曲率)，H(SP1)、H(SP2)分别为曲面切片Sp1、Sp2的平均曲率，ΔK(Sp)为高斯曲率偏差，ΔH(Sp)为平均曲率偏差。
设定曲率偏差阈值为δ1和δ2，其中δ1表示高斯曲率偏差阈值，δ2为平均曲率偏差阈值，则有：
Figure BDA0002214029330000145
需要说明的是，本步骤中若相邻切片Sp1与Sp2之间的曲率偏差小于阈值，说明两片切片弯曲程度相近，高斯曲率偏差阈值δ1和平均曲率偏差阈值δ2的大小可以通过实验仿真求得一个最好效果。
相比于传统的点云分割过程，本发明使用体素滤波直接滤除背景点云，有效减小输入数据规模和运算量，极大的提高了分割速度和鲁棒性；并采用改进的成对连接聚类算法，实现的精切快速的点云主目标分割，且利用了Delaunay三角剖分进行曲面拟合，该方法可以很好地运用于非结构化点云，对于曲率变化快、高次非常规函数的复杂异形曲面点云，不受限于分割物体的形状、结构，能达到较好的分割与曲面拟合效果；进一步地，本发明还通过曲率偏差的曲面合并方法，有效抑制可能的过度分割情况，为后续三维测量提供整体性良好的目标点云。
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (9)

1.一种复杂异形曲面机器人三维测量的点云分割方法，其特征在于，包括以下步骤：
S100、输入以地面、桌面为背景的叶片点云X，利用体素滤波滤除背景点，得到目标叶片点云Y；
所述步骤S100具体包括如下步骤:
S101、将整个点云场景在XY平面上垂直划分为一系列局部点云块Bi，i＝1,2,...,n，n为划分后生成的局部点云块的总数；
S102、根据Octree索引，对于每个局部点云块Bi进行划分，按一定的宽度r将其划分为空间上一系列连续的体素Vj，j＝1,2,...,m，m为划分后生成的体素的总数；
S103、将与体素Vj相邻接并在空间上位于Vj上方的9个邻域体素定义为向上邻域(L1,L2,...,L9)，并随机规定一个的向上邻域为向上生长起点,并按同样的方式，分别沿着Vj的向上邻域继续向上生长，直到所有已执行向上生长操作的体素无向上邻域可生长，向上生长过程结束；
S104、对于已执行向上生长操作的体素，计算出其整体高度值、局部高度值，并选出具有最大局部高度值的体素Vh，并将其定为该生长区域的最高点；
S105、计算地面高度，定义Hg为地面高度的临界值：如果Vh与Hg的高度差小于第一阈值hv1，则Vh被标记为地面体素，同时滤除包含Vh的局部点云块中的所有点；如果Vh与Hg的高度差大于第一阈值hv1，则Vh被标记为非地面体素，同时保留包含Vh的局部点云块中的所有点，为非地面点；
S106、计算桌面高度，定义Hd为桌面高度的临界值：如果Vh与Hd的高度差小于第二阈值hv2，则Vh被标记为桌面体素，同时滤除包含Vh的局部点云块中的所有点；如果Vh与Hd的高度差大于第二阈值hv1，则Vh被标记为非桌面体素，同时保留包含Vh的局部点云块中的所有点，为非桌面点；
S107、非桌面点和非地面点的交集则为目标叶片点云Y；
S200、利用PCA算法计算出Y中点的法向量、平面轮廓度，并剔除离群点，将关联点集合记为一致集CS；
S300、利用法向量、平面轮廓度偏差建立成对连接，确定聚类中心后进行搜索，并搜索所有与它相连的点，生成聚类C；
S400、利用Delaunay三角剖分方法对步骤S300产生的聚类C进行曲面拟合；
S500、对于每一个拟合的曲面切片，计算其曲率，并设定曲率偏差阈值，若相邻两个曲面切片之间的曲率偏差小于阈值，则合并；否则不合并，从而得到与背景点云分割开的完整叶片点云Y。
2.根据权利要求1所述的复杂异形曲面机器人三维测量的点云分割方法，其特征在于，所述步骤S200具体包括如下步骤：
S201、在目标叶片点云Y中随机选取一个点，记为Pi，构建k-d树，搜索出K个与Pi距离最近的临近点，记为KNN，将KNN归为同一个邻域内，此邻域记为Ii；
S202、将Ii中的点按与Pi之间的距离升序排列，选出前
Figure FDA0004006590950000021
个点，利用式(1)计算出协方差矩阵Σ：
Figure FDA0004006590950000022
其中，
Figure FDA0004006590950000023
为平均向量，T表示向量的转置；
S203、利用奇异值分解(SVD)产生标准特征值方程：
λV＝∑V (2)
式中，λ为特征值矩阵，V为特征向量矩阵；
S204、取前三个最大的特征值λ0、λ1、λ2及其对应的特征向量v0、v1、v2，且λ2>λ1>λ0,令λ0表示Pi的平面轮廓度，v0表示Pi所在平面的法向量；
S205、利用如下公式查找Pi的关联点、离群点：
Figure FDA0004006590950000031
式中，Rz表示判断一个点是否为点Pi关联点的权值得分，
Figure FDA0004006590950000032
是Pi的KNN与估计平面的正交距离,NOD是
Figure FDA0004006590950000033
的集合，median(NOD)为NOD的中位数，其中：
Figure FDA0004006590950000034
式中，a为设定常数，取值为1.4826；
当Rz<2.5时，认为该点为Pi的关联点，将其归到一致集CS中；反之为离群点，将其剔除。
3.根据权利要求2所述的复杂异形曲面机器人三维测量的点云分割方法，其特征在于，所述步骤S300具体包括如下步骤：
S301、在Pi的CS中搜索，找出CS中λ比λ(Pi)小的邻域点，并从这些邻域中选择n与n(Pi)偏差最小的邻域点，作为CNP(Pi)，则CNP(Pi)为与点Pi法向量、平面轮廓度偏差最小的点，并建立Pi与CNP(Pi)的成对连接；
S302、若CS中已经不存在λ比λ(Pi)小的邻域点了，则如下公式计算阈值：
Figure FDA0004006590950000035
式中，thλ表示阈值，
Figure FDA0004006590950000036
是所有N个数据点平面轮廓度的平均值，其中：
Figure FDA0004006590950000037
其中，N表示数据点的数量，
Figure FDA0004006590950000038
是所有N个数据点平面轮廓度的平均值；
S303、如果λ(Pi)<thλ，则将Pi记为聚类中心，并归入Ccenter；否则，不记为聚类中心；
S304、找到Ccenter，搜索所有与其间接、直接建立了成对连接的点，并将这些点生成聚类C。
4.根据权利要求3所述的复杂异形曲面机器人三维测量的点云分割方法，其特征在于，所述步骤S400具体包括如下步骤：
S401、针对聚类C中的点，按x坐标从小到大进行排序，根据离散点的最大分布确定包含了点集中所有点的超级三角形；
S402、初始化边缓存数组edge buffer，并遍历temp triangles中的每一个三角形，计算该三角形的外接圆的圆心和半径；
S403、遍历超级三角形中的每一个数据点，若该点在外接圆的右侧，则该三角形为Delaunay三角形，保存至triangles中，并在temp triangles中去掉；若该点在外接圆的内侧，则该三角形不为Delaunay三角形，保存至edge buffer中，并在temp triangles中去掉；若该点在外接圆的左侧，则该三角形为不确定；
S404、对edge buffer数组进行去重，将此时edgebuffer中保留的边与当前的点相结合，形成若干三角形，保存至temp triangles中，并将triangles和temp triangles中的三角形合并，组成最终的超级三角形；
S405、除去与超级三角形三个点有关的所有三角形，对于聚类C进行的曲面拟合完成，将拟合好的曲面记为曲面切片Sp。
5.根据权利要求4所述的复杂异形曲面机器人三维测量的点云分割方法，其特征在于，步骤S400中Delaunay三角剖分方法中三角网格模型通过如下线性表示：
M＝(V,F) (7)
式中，V＝{vi；1≤i≤nv}表示顶点集，F＝{fk；1≤k≤nk}表示三角片集。
6.根据权利要求5所述的复杂异形曲面机器人三维测量的点云分割方法，其特征在于，所述步骤S500中曲面切片Sp的曲率值通过如下步骤计算：
(1)通过如下公式计算曲面切片Sp中每个三角面片fk的法向量Nfk：
Figure FDA0004006590950000041
式中，vi和vj+1分别表示曲面切片Sp上任意选取的两点，
Figure FDA0004006590950000042
表示曲面切片Sp点vi到vj+1的向量，
Figure FDA0004006590950000043
表示曲面切片Sp点vj+1到vj的向量；
(2)通过如下公式计算任意点vi的法向量
Figure FDA0004006590950000044
Figure FDA0004006590950000051
式中，|Ni|表示除顶点vi外其它顶点组成集合的顶点个数，Fi为包含点vi三角面片fk的集合；
(3)根据如下公式计算点vi的法曲率Kij：
Figure FDA0004006590950000052
式中，
Figure FDA0004006590950000053
表示曲面切片Sp点vj到vi的向量；
(4)根据如下公式：
Figure FDA0004006590950000054
式中，K(SP)为高斯曲率，用于反映曲面总的弯曲程度，H(SP)为平均曲率，用于反映空间曲面上某一点任意两个相互垂直的正交曲率的平均值，K1为最大法曲率的极大值，K2为垂直于极大曲率面曲率的极小值。
7.根据权利要求6所述的复杂异形曲面机器人三维测量的点云分割方法，其特征在于，所述最大法曲率通过如下方式选取：遍历三角拟合经过点vi的所有条曲线，计算并选出具有最大法曲率的曲线，此时的法曲率即为最大法曲率。
8.根据权利要求7所述的复杂异形曲面机器人三维测量的点云分割方法，其特征在于，所述步骤S500中相邻两个曲面切片之间的曲率偏差包括高斯曲率偏差和平均曲率偏差，通过如下公式表达：
Figure FDA0004006590950000055
式中，Sp1、Sp2为两个相邻曲面切片，K(SP1)、K(SP2)分别为曲面切片Sp1、Sp2的高斯曲率，H(SP1)、H(SP2)分别为曲面切片Sp1、Sp2的平均曲率，ΔK(Sp)为高斯曲率偏差，ΔH(Sp)为平均曲率偏差。
9.根据权利要求1至8中任一项所述的复杂异形曲面机器人三维测量的点云分割方法，其特征在于，所述步骤S100中叶片点云X通过3D位移传感器扫描获取。
CN201910908581.1A 2019-10-16 2019-10-16 一种复杂异形曲面机器人三维测量的点云分割方法 Active CN110599506B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910908581.1A CN110599506B (zh) 2019-10-16 2019-10-16 一种复杂异形曲面机器人三维测量的点云分割方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910908581.1A CN110599506B (zh) 2019-10-16 2019-10-16 一种复杂异形曲面机器人三维测量的点云分割方法

Publications (2)

Publication Number Publication Date
CN110599506A CN110599506A (zh) 2019-12-20
CN110599506B true CN110599506B (zh) 2023-03-24

Family

ID=68863182

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910908581.1A Active CN110599506B (zh) 2019-10-16 2019-10-16 一种复杂异形曲面机器人三维测量的点云分割方法

Country Status (1)

Country Link
CN (1) CN110599506B (zh)

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111227444A (zh) * 2020-01-17 2020-06-05 泉州装备制造研究所 一种基于k最近邻的3D鞋底喷胶路径规划方法
CN111496789B (zh) * 2020-04-23 2021-09-28 佛山科学技术学院 一种离线复杂曲面喷涂轨迹规划系统及控制方法
CN112446952B (zh) * 2020-11-06 2024-01-26 杭州易现先进科技有限公司 三维点云法向量的生成方法、装置、电子设备及存储介质
CN112308761A (zh) * 2020-11-13 2021-02-02 济南浪潮高新科技投资发展有限公司 基于切片与插值的机器人深度相机点云降采样滤波方法
CN113239500B (zh) * 2021-07-12 2021-09-21 四川大学 基于协方差矩阵的参考点邻域特征匹配方法
CN116197910B (zh) * 2023-03-16 2024-01-23 江苏集萃清联智控科技有限公司 风电叶片轮式移动打磨机器人的环境感知方法和装置
CN116912312B (zh) * 2023-09-15 2023-12-01 湖南大学 一种面向复杂曲面构件的三维孔类定位方法
CN117593515B (zh) * 2024-01-17 2024-03-29 中数智科(杭州)科技有限公司 一种轨道车辆用螺栓松动检测系统、方法及存储介质

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105957076A (zh) * 2016-04-27 2016-09-21 武汉大学 一种基于聚类的点云分割方法及系统
CN107369161A (zh) * 2017-07-19 2017-11-21 无锡信捷电气股份有限公司 一种基于改进欧式聚类的散乱工件点云分割方法
WO2019174236A1 (zh) * 2018-03-14 2019-09-19 浙江大学 一种基于ViBe的三维声纳点云图像分割方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105957076A (zh) * 2016-04-27 2016-09-21 武汉大学 一种基于聚类的点云分割方法及系统
CN107369161A (zh) * 2017-07-19 2017-11-21 无锡信捷电气股份有限公司 一种基于改进欧式聚类的散乱工件点云分割方法
WO2019174236A1 (zh) * 2018-03-14 2019-09-19 浙江大学 一种基于ViBe的三维声纳点云图像分割方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于特征融合的林下环境点云分割;樊丽等;《北京林业大学学报》;20160531;第38卷(第05期);第133-138页 *

Also Published As

Publication number Publication date
CN110599506A (zh) 2019-12-20

Similar Documents

Publication Publication Date Title
CN110599506B (zh) 一种复杂异形曲面机器人三维测量的点云分割方法
CN108010116B (zh) 点云特征点检测方法和点云特征提取方法
CN109887015B (zh) 一种基于局部曲面特征直方图的点云自动配准方法
CN111080684B (zh) 一种点邻域尺度差异描述的点云配准方法
Shi et al. Adaptive simplification of point cloud using k-means clustering
Ji et al. A novel simplification method for 3D geometric point cloud based on the importance of point
CN100559398C (zh) 自动的深度图像配准方法
CN109325993B (zh) 一种基于类八叉树索引的显著性特征强化采样方法
US20070036434A1 (en) Topology-Based Method of Partition, Analysis, and Simplification of Dynamical Images and its Applications
CN114926699B (zh) 一种室内三维点云语义分类方法、装置、介质及终端
CN110222642A (zh) 一种基于全局图聚类的平面建筑构件点云轮廓提取方法
CN107680168B (zh) 三维重建中基于平面拟合的网格简化方法
CN109523582B (zh) 一种顾及法向量和多尺度稀疏特征的点云粗配准方法
CN111275724A (zh) 一种基于八叉树和边界优化的机载点云屋顶平面分割方法
CN111652855A (zh) 一种基于存活概率的点云精简方法
CN109919955A (zh) 地基式激光雷达点云的隧道轴线提取和分割方法
CN108595631B (zh) 基于图谱理论的三维cad模型双层检索方法
CN116402866A (zh) 基于点云的零件数字孪生几何建模与误差评定方法及系统
CN114663373A (zh) 一种用于零件表面质量检测的点云配准方法及装置
CN106355178B (zh) 基于分层聚类和拓扑连接模型的海量点云自适应简化方法
CN112241676A (zh) 一种地形杂物自动识别的方法
Lu et al. Automatic point cloud registration algorithm based on the feature histogram of local surface
CN108629315B (zh) 一种针对三维点云的多平面识别方法
CN112734934B (zh) 一种基于相交边映射的stl模型3d打印切片方法
CN114943130A (zh) 分割表示机械部件的3d建模对象

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant