CN113763563A - 一种基于平面识别的三维点云几何网格结构生成方法 - Google Patents

一种基于平面识别的三维点云几何网格结构生成方法 Download PDF

Info

Publication number
CN113763563A
CN113763563A CN202111044729.5A CN202111044729A CN113763563A CN 113763563 A CN113763563 A CN 113763563A CN 202111044729 A CN202111044729 A CN 202111044729A CN 113763563 A CN113763563 A CN 113763563A
Authority
CN
China
Prior art keywords
points
plane
boundary
point cloud
triangle
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.)
Pending
Application number
CN202111044729.5A
Other languages
English (en)
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.)
Daiwu Intelligent Technology Shanghai Co ltd
Original Assignee
Daiwu Intelligent Technology Shanghai Co ltd
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 Daiwu Intelligent Technology Shanghai Co ltd filed Critical Daiwu Intelligent Technology Shanghai Co ltd
Priority to CN202111044729.5A priority Critical patent/CN113763563A/zh
Publication of CN113763563A publication Critical patent/CN113763563A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation
    • G06T17/205—Re-meshing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
    • G06F18/213—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods
    • G06F18/2135—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods based on approximation criteria, e.g. principal component analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/25—Fusion techniques
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q50/00—Systems or methods specially adapted for specific business sectors, e.g. utilities or tourism
    • G06Q50/08—Construction
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/005—Tree description, e.g. octree, quadtree
    • 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/13—Edge detection
    • 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

Abstract

本发明公开了一种基于平面识别的三维点云几何网格结构生成方法，所述方法包括以下步骤：步骤S1：数据预处理，输入数据为通用的三维空间点云数据，对所述输入数据进行降采样，使用主成分分析来计算点云的法向量值；步骤S2：平面识别，构建能量函数，得到局部一致的平面提取结果，并对结果进行优化，用主成分分析计算其所述内点在空间中的分布，用优化步骤得到的法向量与经过点即可表示平面；步骤S3：求解平面自身轮廓边界；步骤S4：求解平面相交边界；步骤S5：边界融合，得到一个平面的闭合边界；步骤S6：网格生成。本发明可以对平面识别算法的改进与优化；可以实现平面多边形轮廓的双方法生成与融合。

Description

一种基于平面识别的三维点云几何网格结构生成方法
技术领域
本发明涉及计算机三维重建领域，尤其涉及一种基于平面识别的三维点云几何网格结构生成方法。
背景技术
三维重建的过程中，在扫描得到点云后，下一步通常是要三角化得到网格化的模型。原因是点云不管多密，在三维空间中实际上还是一种稀疏的表达方式，即点与点之间会有距离的存在，而网格化的表示则是非稀疏的，能够闭合、完整地表达一个物体的表面。同时，对于一个巨大的长方形平面，如果用点来表达，则需要大量的点，但用三角形的话只需要两个，所以也可以省下大量的存储空间。网格(mesh)：指模型的网格，3D模型是由多边形拼接而成，而多边形实际上又是由多个三角形拼接而成的。即一个3D模型的表面其实是由多个彼此相连的三角面构成。
目前常用的建立三角网格的方法是泊松重建或是截断有向距离函数(TSDF)，它们都是通过一些方法计算栅格化的三维空间中各提速点到网格正面的有向距离，并通过移动立方体法(Marching Cubes)来计算三维离散体素中的等值面，得到网格。
现有方法都是把三维空间栅格化后计算其中的三维面，那么计算过程中的内存消耗、结果的精度都严重依赖于栅格化过程中八叉树的精度。
一个经过多个体素的面会由多个小三角形组成，这虽然能够表示具有弧度的面，但是面对一个大型的平直面，例如一整面墙时，并不需要大量三角形来变现凹凸，同时两个大三角形能在完整表现墙面的同时节省大量的空间。
尤其时泊松重建，由于其中平滑函数的存在，最后的结果中会有较多的圆角，而非建筑结构中应有的那些90度的边缘。
因此，本领域的技术人员致力于开发一种基于平面识别的三维点云几何网格结构生成方法，本方法并未将三维空间栅格化，因而不会出现类似的精度限制与内存消耗；同时在平面所有方向都与其它平面相交的情况下，边缘完全由交线构成，能形成简洁的网格；而相交带来的交线也保证了边角的锐度。
发明内容
有鉴于现有技术的上述缺陷，本发明所要解决的技术问题涉及计算机视觉中的三维重建领域的网格化部分，具体地，在用三维传感器扫描点得到环境场景点云后，如何三角化环境点云，并得到场景的网格化表示。针对该问题，本发明在对三维点云进行面识别后，基于每个面自身的边界以及面与面之间的相交关系，建立出精简可靠的网格模型。
为实现上述目的，本发明提供了一种基于平面识别的三维点云几何网格结构生成方法，所述方法包括以下步骤：
步骤S1：数据预处理，输入数据为通用的三维空间点云数据，对所述输入数据进行降采样，使用主成分分析来计算点云的法向量值；
步骤S2：平面识别，构建能量函数，得到局部一致的平面提取结果，并对结果进行优化，用主成分分析计算其所述内点在空间中的分布，用优化步骤得到的法向量与经过点即可表示平面；
步骤S3：求解平面自身轮廓边界；
步骤S4：求解平面相交边界，对每个面的所有点建立k维树，双重便利所有的面，对任意两个面都计算彼此距离另一个面的点集的距离小于0.05m的点的数量，如果点数不到10则说明两者不相交；如果两个面的点的数量都大于等于10个就算这两个面法向量夹角，夹角小于20度就不相交，否则判定为相交；
步骤S5：边界融合，得到一个平面的闭合边界；
步骤S6：网格生成，对每一个平面的多边形边界，使用耳切法计算其网格，然后把所有平面的网格合并，则可以得到最终的环境网格数据。
进一步地，所述降采样的方法是体素滤波，其通过体素化三维空间，来把所有的点归到三维中的不同小方块中，并用小方块的重心来近似地表达其内的所有点，得到过滤后的三维点云。
进一步地，所述主成分分析来计算点云的法向量值，其主要步骤为：
步骤S11：对数据集X＝{x1,x2,…,xn}进行去中心化，即每个特征维度减去各自的平均值，其中x是向量，X是向量组成的矩阵；
步骤S12：计算协方差矩阵
Figure BDA0003250745980000021
步骤S13：用特征值分解方法求协方差矩阵
Figure BDA0003250745980000022
的特征值与特征向量，并按特征值从大到小排序得到特征值λ1,λ2,λ3与特征向量
Figure BDA0003250745980000023
其中
Figure BDA0003250745980000024
即为法向量。
进一步地，所述步骤S2中结果优化包含以下步骤：
若区域内点在法向量的发向上的长度大于0.1m，则判断其可能存在双层墙，调用基于密度的聚类算法对其进行聚类，其在空间中使用k维树加速，具体步骤为：
步骤S21：设定参数eps＝0.05m，min_samples根据疏密调整，其中eps和min_samples是算法中的两个参数；
步骤S22：随机选取一个点，找到所有到这个点距离≤eps的点，若点个数小于min_samples，则该点被标记为噪声，否则被标记为核心，并被分配一个新的簇标签；
步骤S23：若被标记为核心，则访问eps范围内的所有邻居，如果它们还没有被分配一个簇，则将刚创建的新簇标签分配给它们，如果它们是核心样本，那么就依次访问其邻居，以此类推，直到在簇的eps距离内没有更多的核心样本为止；
步骤S24：回到步骤S22，直到所有点都被访问。
进一步地，所述步骤S2中计算所述内点在空间中的分布，还包括以下步骤：
用主成分分析计算其所述内点在空间中的分布，得到对应特征值最大的特征向量
Figure BDA0003250745980000031
结合该区域法向量
Figure BDA0003250745980000032
则可得到子区域坐标系到全局坐标系的旋转矩阵Rwn与从全局坐标系到子区域坐标系的旋转矩阵Rnw：
Figure BDA0003250745980000033
Figure BDA0003250745980000034
用Rnw将区域内点旋转到子区域坐标系内，排序得到z值为中位数z＝median(zi)的点，将其用Rwn旋转回全局坐标系，则得到平面在全局坐标系中经过的点pmedian。
进一步地，所述步骤S3中所述计算平面自身轮廓边界，还包括以下子步骤：
先用Rnw把这个平面的点投到局部坐标系并仅保留x、y坐标，用阿尔法形状边缘轮廓算法做出粗略边界，需要设定半径值alpha，子步骤如下：
S31：根据点集建立德洛内三角网；
S32：若三角形中某条边的长度大于2*alpha，则删除该三角形；
S33：对三角形的每条边进行判断：若过某条边的两点且半径为alpha的圆包含其他点，则删除该三角形；
在第所述步骤S32、S33步中提到的不符合阿尔法形状边缘轮廓算法要求的三角形后所得到的三角网上求出三角网的边缘，该边缘即为点集的边缘线。
进一步地，所述步骤S4还包括：
如果两个面存在交线，假设两个平面的法向量为
Figure BDA0003250745980000035
和
Figure BDA0003250745980000036
则定义旋转矩阵：
Figure BDA0003250745980000037
Figure BDA0003250745980000038
用R′nw把两个平面旋转到新的局部坐标系中，在新坐标系里求出交线经过的点的坐标，然后再用R′wn旋转回原坐标系。
进一步地，所述步骤S4中还包括：
把两边的领域点投到直线上，取两边方向的最大最小值作为端点；同时每个面都要记录所有由它相交得到的线段，具体是：
步骤S41：每个面用Rnw法向量投到二维，然后在二维平面上把这平面上的所有线两两相交取交点，比较交点和端点的距离，如果接近就把端点拉过去，这样每一根线就会在每个产生它的平面上具有一个矫正过的投影，再用Rwn旋转回原坐标系；
步骤S42：在全局三维坐标系中，每根线的端点取它们在每个平面上矫正值的端点的平均值。
进一步地，所述步骤S5还包括以下子步骤：
步骤S51：对于每一个面所属的所有交线边界，查找自身边界上距离其两个端点最近的两个点，若距离小于阈值，且自身边界这两个端点之间超过80％的点距离该线段的距离都小于同样的阈值，则用这段交线替换这段粗略边界；
步骤S52：替换完成后，对于剩余的自身边界，遍历并找出相连的斜率相近的短线段，若该线段集合两边的端点都属于所述交线边界，则直接将其作为一根线段相连；反之则用该线段集合的端点来拟合直线，从而得到新线段方向，把两端属于自身边界的点投影到直线上，得到新的线段。
进一步地，所述步骤S6中的耳切法流程为：
步骤S61：建立所有顶点组成的双向链表；
步骤S62：循环双向链表，对于顶点
Figure BDA0003250745980000041
检查其是否为凸点，以及链表内是否由其它点在三角形
Figure BDA0003250745980000042
内，如为凸点且无链表内其它点，则将三角形
Figure BDA0003250745980000043
添加到网格中，并从链表中删除vi；
步骤S63：循环直到链表中只剩下最后一个三角形，加入到网格中；
其中，判断是否为凸点的方法是计算
Figure BDA0003250745980000044
若为正则凸，否则为凹；判断其它点是否在三角形内则使用重心法，具体地，对任一点
Figure BDA0003250745980000045
建立
Figure BDA0003250745980000046
可解得
Figure BDA0003250745980000047
若λ2＝D1/D、λ2＝D2/D、1-λ2-λ3都在0与1之间，则该点在三角形内。
其中，
Figure BDA0003250745980000048
是平面上的顶点，
Figure BDA0003250745980000049
是三个点组成的三角形，
Figure BDA00032507459800000410
和
Figure BDA00032507459800000411
是两个点组成的向量，D、D1、D2、λ2、λ3都是普通的标量。
本发明可以对平面识别算法的改进与优化；可以实现平面多边形轮廓的双方法生成与融合，较现有技术均有实质性的进步。对方向做出非强制性限制，使其能够适应一些非曼哈顿的结构，例如“人”字型的建筑楼层，同时又不会散乱；生成的边界基于自身轮廓与相交关系，保证了每个面的闭合以及交线处的锐度；耳切法生成的网格是基于提供的多边形下最精简的，保证了结果的简洁性。
以下将结合附图对本发明的构思、具体结构及产生的技术效果作进一步说明，以充分地了解本发明的目的、特征和效果。
附图说明
图1是本发明的一个较佳实施例的方法流程示意图；
图2是本发明的一个较佳实施例的融合过程示意图；
图3是本发明的一个较佳实施例的平面边界轮廓生成网格示意图。
具体实施方式
以下参考说明书附图介绍本发明的多个优选实施例，使其技术内容更加清楚和便于理解。本发明可以通过许多不同形式的实施例来得以体现，本发明的保护范围并非仅限于文中提到的实施例。
在附图中，结构相同的部件以相同数字标号表示，各处结构或功能相似的组件以相似数字标号表示。附图所示的每一组件的尺寸和厚度是任意示出的，本发明并没有限定每个组件的尺寸和厚度。为了使图示更清晰，附图中有些地方适当夸大了部件的厚度。
本发明的目的是提供一种三维重建过程中栅格模型的生成方法，其能够通过扫描得到的三维点云，基于平面自身信息与相交关系得到精简可靠的栅格模型。
本发明的方法流程如图1所示，包含数据预处理、平面识别、确定平面自身轮廓边界、确定平面相交边界、边界融合和网格生成步骤。
步骤S1：数据预处理
本发明所需的输入数据为通用的三维空间点云数据，其一般基于三维传感器以及相应slam算法扫描融合得到，原始数据中并不需要具备法向量(Normal)值，因为首先需要对数据进行降采样，来避免点云的疏密影响平面的计算，并提高平面检测的计算效率。
本发明中使用的降采样方法是体素滤波，其通过体素化三维空间，来把所有的点归到三维中的不同小方块中，并用小方块的重心来近似地表达其内的所有点，得到过滤后的三维点云。在本发明中，为了在加速平面识别的同时尽可能保持袭击，本发明选择体素立方体的长宽高大小为0.03m。
在降采样后，本发明使用目前被广泛地使用的主成分分析来计算点云的法向量值，其主要步骤为：
①对数据集X＝{x1,x2,…,xn}进行去中心化，即每个特征维度减去各自的平均值；
②计算协方差矩阵
Figure BDA0003250745980000051
③用特征值分解方法求协方差矩阵
Figure BDA0003250745980000052
的特征值与特征向量，并按特征值从大到小排序得到特征值λ1,λ2,λ3与特征向量
Figure BDA0003250745980000053
其中
Figure BDA0003250745980000054
即为法向量。
其中，x是向量，比如一个点的(x,y,z)，大写的X是向量组成的矩阵，就是整个点集，右上角有T的是转置，λ1,λ2,λ3是特征值，
Figure BDA0003250745980000055
是特征向量。
而法向量能用这种统计中的方法计算的原因则是该平面上的所有向量都应该垂直于法向量的方向，即：
Figure BDA0003250745980000061
步骤S2：平面识别
本发明使用了一种已有的平面识别算法，其首次被提出时在2020年于“Fastregularity-constrained plane fitting”一文中。其核心思想在于构建了一个兼顾了全局与局部的能量函数，使其能够在限制了全局方向总数的前提下得到局部一致的平面提取结果。
该方法中作为核心的能量函数可以被表示为：
Figure BDA0003250745980000062
其中N是点的数量，Ni是第i个点的邻域点集合，V是重建后的法向量集合，zi是第i个点的法向量在V中的序号，I是预处理阶段计算得到的法向量集合，Ii则是第i个点的原始法向量。
在该式中，||Vzi-Ii||2和
Figure BDA0003250745980000063
分别是为了引导输出的法向量能够尽可能的与输入一致，以及限制空间中相邻的点最后被分配到同样的法向量，该式前的λ是限制在能量函数中所占的比重，通常根据数据情况在1到10之间调整。为了便于优化，该问题被拆分为了两个子问题：
Figure BDA0003250745980000064
Figure BDA0003250745980000065
即第一个子问题着重于优化法向量方向，第二个子问题则在尽可能保证法向量准确的情况下，把法向量方向的数量限制在m个以内。
在该方法中，第一个式子通过合并相邻平面域(构成一个平面的点集)的方法进行优化，其计算合并两个子域时能量函数的变化：
Figure BDA0003250745980000066
其中cij为两个子域中相连单元的个数，而wi则为域i的点集数量。如果该式结果为正，则能量函数会减小，即合并两子域，否则不合并。合并后新子域的法向量更新为：
Figure BDA0003250745980000067
而为了优化第二个式子，即挑选子区域，通过构建满足：
Figure BDA0003250745980000071
的次模函数：
Figure BDA0003250745980000072
其中
Figure BDA0003250745980000073
是不相连区域的集合，
Figure BDA0003250745980000074
是第i个子区域的法向量，
Figure BDA0003250745980000075
是包含
Figure BDA0003250745980000076
的法向量集合。然后按照以下步骤进行优化：
①初始化
Figure BDA0003250745980000077
②遍历区域内点数量大于阈值的所有区域，对第i个区域
Figure BDA0003250745980000078
若
Figure BDA0003250745980000079
则X←X∪{Si}，
Figure BDA00032507459800000710
反之
Figure BDA00032507459800000711
最后得到的V即为结果。显然，这是为了挑选增加的贡献比减少的贡献更大的子区域
该方法的结果相比传统的区域增长存在优势的同时也具有缺陷。比如由于限制了法向量的数量，双层墙在扫描的结果中会被归到一起。本发明为了解决这一问题，若区域内点在法向量的发向上的长度大于0.1m，则判断其可能存在双层墙，调用DBSCAN对其进行聚类，其在空间中使用kdtree加速，大致流程如下：
①设定参数eps＝0.05m，min_samples根据疏密调整；
②随机选取一个点，找到所有到这个点距离≤eps的点，若点个数小于min_samples，则该点被标记为噪声，否则被标记为核心，并被分配一个新的簇标签；
③若被标记为核心，则访问eps范围内的所有邻居，如果它们还没有被分配一个簇，则将刚创建的新簇标签分配给它们，如果它们是核心样本，那么就依次访问其邻居，以此类推，直到在簇的eps距离内没有更多的核心样本为止。
④回到步骤2，直到所有点都被访问。
其中，kdtree是k维树的意思，简单说就是把2叉树的概念扩展到了k维空间中。其所有非叶子节点都可以看作是用一个超平面把该处空间里的点分为两部分。
其中，DBSCAN是一种基于密度的聚类算法，eps和min_samples是算法中的两个参数名。
该方法相比传统的区域增长却有两个极大的优点：由于其具有全局的限制，平面的法向量不会散乱，而局部的限制则驱使其平面交线附近不会丢失点。由于平面交线部分的法向量通过PCA计算通常是斜的，所以传统区域增长的结果不会包括这些点，而该方法则会包括。
平面的表达除了法向量之外，还有经过的点。对此，本发明用主成分分析(PCA)计算其内点在空间中的分布，得到对应特征值最大的特征向量
Figure BDA00032507459800000712
结合该区域法向量
Figure BDA00032507459800000713
则可得到子区域坐标系到全局坐标系的旋转矩阵Rwn与从全局坐标系到子区域坐标系的旋转矩阵Rnw：
Figure BDA0003250745980000081
Figure BDA0003250745980000082
用Rnw将区域内点旋转到子区域坐标系内，排序得到z值为中位数z＝median(zi)的点，将其用Rwn旋转回全局坐标系，则得到平面在全局坐标系中经过的点pmedian。
至此，用优化步骤得到的法向量与经过点即可表示平面。
步骤S3：平面自身轮廓边界
先用Rnw把这个平面的点投到局部坐标系并仅保留x、y坐标，用alpha shape做出粗略边界，其是一种既有的二维边界运算方法，需要设定半径值alpha(一般设定为0.3m)，简要步骤如下：
①根据点集建立Delaunay三角网；
②若三角形中某条边的长度大于2*alpha，则删除该三角形；
③对三角形的每条边进行判断：若过某条边的两点且半径为alpha的圆包含其他点，则删除该三角形；
在第2-3步提到的不符合Alpha shape要求的三角形后所得到的三角网上求出三角网的边缘，该边缘即为点集的边缘线。
其中，alpha shape(阿尔法形状)是一种求边缘轮廓的算法；Delaunay(德洛内三角网)，其是指把一个点集三角剖分为只包含一种叫Delaunay边的三角化方式。Delaunay边指经过一条边两个端点的圆中不含有点集中的任何其它点。
步骤S4：平面相交边界
对每个面的所有点建立kdtree，然后双重便利所有的面，对任意两个面都计算彼此距离另一个面的点集的距离小于0.05m的点的数量，如果点数不到10则说明两者不相交；如果两个面的点的数量都大于等于10个就算这两个面法向量夹角，夹角小于20度就不相交，否则判定为相交。
如果两个面存在交线，假设两个平面的法向量为
Figure BDA0003250745980000083
和
Figure BDA0003250745980000084
则定义旋转矩阵：
Figure BDA0003250745980000085
Figure BDA0003250745980000086
用R′nw把两个平面旋转到新的局部坐标系中，在新坐标系里求出交线经过的点的坐标，然后再用R′wn旋转回原坐标系。这样做的原因是在原坐标系里想把两个平面方程联立求交点的话要把所有未知数前面的系数是否为0判断一遍，否则直接除可能会分母为0，即需要考虑2^6/2种情况，然后新坐标系里就直接转成了一个平面类似y＝b的形式，代到另一个里就结束了。
求出直线方程后把两边的领域点投到这根直线上，取两边方向的最大最小值作为端点。同时每个面都要记录所有由它相交得到的线段。
由于旋转等运算会产生计算误差，而且端点本身就由于稀疏性具有误差，接下来的一步是把因为误差而可能分开的线段端点融合，具体是：
①每个面用Rnw法向量投到二维，然后在二维平面上把这平面上的所有线两两相交取交点，比较交点和端点的距离，如果接近就把端点拉过去，这样每一根线就会在每个产生它的平面上具有一个矫正过的投影，再用Rwn旋转回原坐标系。
②然后在全局三维坐标系中，每根线的端点取它们在每个平面上矫正值的端点的平均值。
步骤S5：边界融合
对于每一个面所属的所有交线边界，查找自身边界上距离其两个端点最近的两个点，若距离小于阈值(可取0.5m)，且自身边界这两个端点之间超过80％的点距离该线段的距离都小于同样的阈值，则用这段交线替换这段粗略边界。这段比较主要是为了防止非边界的交线，比如楼层天花板和楼层中心一个房间的某一面墙。
替换完成后，对于剩余的自身边界，遍历并找出相连的斜率相近的短线段，若该线段集合两边的端点都属于交线边界，则直接将其作为一根线段相连；反之则用该线段集合的端点来拟合直线，从而得到新线段方向，把两端属于自身边界的点投影到直线上，得到新的线段。
如此，就得到了一个平面的闭合边界。图2即为融合过程的简单示意。
步骤S6：网格生成
对每一个平面的多边形边界，使用耳切法(Ear Clipping)计算其网格，其流程为：
①建立所有顶点组成的双向链表；
②循环双向链表，对于顶点
Figure BDA0003250745980000091
检查其是否为凸点，以及链表内是否由其它点在三角形
Figure BDA0003250745980000092
内，如为凸点且无链表内其它点，则将三角形
Figure BDA0003250745980000093
添加到网格中，并从链表中删除vi；
③循环直到链表中只剩下最后一个三角形，加入到网格中。
其中，判断是否为凸点的方法是计算
Figure BDA0003250745980000094
若为正则凸，否则为凹；判断其它点是否在三角形内则使用重心法，具体地，对任一点
Figure BDA0003250745980000095
建立
Figure BDA0003250745980000096
可解得：
Figure BDA0003250745980000097
若λ2＝D1/D、λ2＝D2/D、1-λ2-λ3都在0与1之间，则该点在三角形内。
其中，
Figure BDA0003250745980000098
是平面上的顶点，
Figure BDA0003250745980000099
是三个点组成的三角形，
Figure BDA00032507459800000910
和
Figure BDA00032507459800000911
是两个点组成的向量，D、D1、D2、λ2、λ3都是普通的标量。
图3即为一个平面边界轮廓生成网格的示意图。
由于各平面之间相交部分的点都经过矫正，是共有的顶点，所以直接把所有平面的网格合并，则可以得到最终的环境网格数据。
在用常规方法建立网格后，也可合并相邻的法向量相近的面。
以上详细描述了本发明的较佳具体实施例。应当理解，本领域的普通技术无需创造性劳动就可以根据本发明的构思作出诸多修改和变化。因此，凡本技术领域中技术人员依本发明的构思在现有技术的基础上通过逻辑分析、推理或者有限的实验可以得到的技术方案，皆应在由权利要求书所确定的保护范围内。

Claims (10)

1.一种基于平面识别的三维点云几何网格结构生成方法，其特征在于，所述方法包括以下步骤：
步骤S1：数据预处理，输入数据为通用的三维空间点云数据，对所述输入数据进行降采样，使用主成分分析来计算点云的法向量值；
步骤S2：平面识别，构建能量函数，得到局部一致的平面提取结果，并对结果进行优化，用主成分分析计算其所述内点在空间中的分布，用优化步骤得到的法向量与经过点即可表示平面；
步骤S3：求解平面自身轮廓边界；
步骤S4：求解平面相交边界，对每个面的所有点建立k维树，双重便利所有的面，对任意两个面都计算彼此距离另一个面的点集的距离小于0.05m的点的数量，如果点数不到10则说明两者不相交；如果两个面的点的数量都大于等于10个就算这两个面法向量夹角，夹角小于20度就不相交，否则判定为相交；
步骤S5：边界融合，得到一个平面的闭合边界；
步骤S6：网格生成，对每一个平面的多边形边界，使用耳切法计算其网格，然后把所有平面的网格合并，则可以得到最终的环境网格数据。
2.如权利要求1所述的基于平面识别的三维点云几何网格结构生成方法，其特征在于，所述降采样的方法是体素滤波，其通过体素化三维空间，来把所有的点归到三维中的不同小方块中，并用小方块的重心来近似地表达其内的所有点，得到过滤后的三维点云。
3.如权利要求1所述的基于平面识别的三维点云几何网格结构生成方法，其特征在于，所述主成分分析来计算点云的法向量值，其主要步骤为：
步骤S11：对数据集X＝{x1,x2,…,xn}进行去中心化，即每个特征维度减去各自的平均值，其中x是向量，X是向量组成的矩阵；
步骤S12：计算协方差矩阵
Figure FDA0003250745970000011
步骤S13：用特征值分解方法求协方差矩阵
Figure FDA0003250745970000012
的特征值与特征向量，并按特征值从大到小排序得到特征值λ1,λ2,λ3与特征向量
Figure FDA0003250745970000013
其中
Figure FDA0003250745970000014
即为法向量。
4.如权利要求1所述的基于平面识别的三维点云几何网格结构生成方法，其特征在于，所述步骤S2中结果优化包含以下步骤：
若区域内点在法向量的发向上的长度大于0.1m，则判断其可能存在双层墙，调用基于密度的聚类算法对其进行聚类，其在空间中使用k维树加速，具体步骤为：
步骤S21：设定参数eps＝0.05m，min_samples根据疏密调整，其中eps和min_samples是算法中的两个参数；
步骤S22：随机选取一个点，找到所有到这个点距离≤eps的点，若点个数小于min_samples，则该点被标记为噪声，否则被标记为核心，并被分配一个新的簇标签；
步骤S23：若被标记为核心，则访问eps范围内的所有邻居，如果它们还没有被分配一个簇，则将刚创建的新簇标签分配给它们，如果它们是核心样本，那么就依次访问其邻居，以此类推，直到在簇的eps距离内没有更多的核心样本为止；
步骤S24：回到步骤S22，直到所有点都被访问。
5.如权利要求1所述的基于平面识别的三维点云几何网格结构生成方法，其特征在于，所述步骤S2中计算所述内点在空间中的分布，还包括以下步骤：
用主成分分析计算其所述内点在空间中的分布，得到对应特征值最大的特征向量
Figure FDA0003250745970000021
结合该区域法向量
Figure FDA0003250745970000022
则可得到子区域坐标系到全局坐标系的旋转矩阵Rwn与从全局坐标系到子区域坐标系的旋转矩阵Rnw：
Figure FDA0003250745970000023
Figure FDA0003250745970000024
用Rnw将区域内点旋转到子区域坐标系内，排序得到z值为中位数z＝median(zi)的点，将其用Rwn旋转回全局坐标系，则得到平面在全局坐标系中经过的点pmedian。
6.如权利要求1所述的基于平面识别的三维点云几何网格结构生成方法，其特征在于，所述步骤S3中所述计算平面自身轮廓边界，还包括以下子步骤：
先用Rnw把这个平面的点投到局部坐标系并仅保留x、y坐标，用阿尔法形状边缘轮廓算法做出粗略边界，需要设定半径值alpha，子步骤如下：
S31：根据点集建立德洛内三角网；
S32：若三角形中某条边的长度大于2*alpha，则删除该三角形；
S33：对三角形的每条边进行判断：若过某条边的两点且半径为alpha的圆包含其他点，则删除该三角形；
在第所述步骤S32、S33步中提到的不符合阿尔法形状边缘轮廓算法要求的三角形后所得到的三角网上求出三角网的边缘，该边缘即为点集的边缘线。
7.如权利要求1所述的基于平面识别的三维点云几何网格结构生成方法，其特征在于，所述步骤S4还包括：
如果两个面存在交线，假设两个平面的法向量为
Figure FDA0003250745970000025
和
Figure FDA0003250745970000026
则定义旋转矩阵：
Figure FDA0003250745970000031
Figure FDA0003250745970000032
用R′nw把两个平面旋转到新的局部坐标系中，在新坐标系里求出交线经过的点的坐标，然后再用R′wn旋转回原坐标系。
8.如权利要求1所述的基于平面识别的三维点云几何网格结构生成方法，其特征在于，所述步骤S4中还包括：
把两边的领域点投到直线上，取两边方向的最大最小值作为端点；同时每个面都要记录所有由它相交得到的线段，具体是：
步骤S41：每个面用Rnw法向量投到二维，然后在二维平面上把这平面上的所有线两两相交取交点，比较交点和端点的距离，如果接近就把端点拉过去，这样每一根线就会在每个产生它的平面上具有一个矫正过的投影，再用Rwn旋转回原坐标系；
步骤S42：在全局三维坐标系中，每根线的端点取它们在每个平面上矫正值的端点的平均值。
9.如权利要求1所述的基于平面识别的三维点云几何网格结构生成方法，其特征在于，所述步骤S5还包括以下子步骤：
步骤S51：对于每一个面所属的所有交线边界，查找自身边界上距离其两个端点最近的两个点，若距离小于阈值，且自身边界这两个端点之间超过80％的点距离该线段的距离都小于同样的阈值，则用这段交线替换这段粗略边界；
步骤S52：替换完成后，对于剩余的自身边界，遍历并找出相连的斜率相近的短线段，若该线段集合两边的端点都属于所述交线边界，则直接将其作为一根线段相连；反之则用该线段集合的端点来拟合直线，从而得到新线段方向，把两端属于自身边界的点投影到直线上，得到新的线段。
10.如权利要求1所述的基于平面识别的三维点云几何网格结构生成方法，其特征在于，所述步骤S6中的耳切法流程为：
步骤S61：建立所有顶点组成的双向链表；
步骤S62：循环双向链表，对于顶点
Figure FDA0003250745970000033
检查其是否为凸点，以及链表内是否由其它点在三角形
Figure FDA0003250745970000034
内，如为凸点且无链表内其它点，则将三角形
Figure FDA0003250745970000035
添加到网格中，并从链表中删除vi；
步骤S63：循环直到链表中只剩下最后一个三角形，加入到网格中；
其中，判断是否为凸点的方法是计算
Figure FDA0003250745970000036
若为正则凸，否则为凹；判断其它点是否在三角形内则使用重心法，具体地，对任一点
Figure FDA0003250745970000037
建立
Figure FDA0003250745970000041
可解得
Figure FDA0003250745970000042
若λ2＝D1/D、λ2＝D2/D、1-λ2-λ3都在0与1之间，则该点在三角形内；
其中，
Figure FDA0003250745970000043
是平面上的顶点，
Figure FDA0003250745970000044
是三个点组成的三角形，
Figure FDA0003250745970000045
和
Figure FDA0003250745970000046
是两个点组成的向量，D、D1、D2、λ2、λ3都是普通的标量。
CN202111044729.5A 2021-09-07 2021-09-07 一种基于平面识别的三维点云几何网格结构生成方法 Pending CN113763563A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111044729.5A CN113763563A (zh) 2021-09-07 2021-09-07 一种基于平面识别的三维点云几何网格结构生成方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111044729.5A CN113763563A (zh) 2021-09-07 2021-09-07 一种基于平面识别的三维点云几何网格结构生成方法

Publications (1)

Publication Number Publication Date
CN113763563A true CN113763563A (zh) 2021-12-07

Family

ID=78793551

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111044729.5A Pending CN113763563A (zh) 2021-09-07 2021-09-07 一种基于平面识别的三维点云几何网格结构生成方法

Country Status (1)

Country Link
CN (1) CN113763563A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116935231A (zh) * 2023-09-14 2023-10-24 湖北工业大学 一种隧道围岩结构面信息提取及关键块体识别方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103985155A (zh) * 2014-05-14 2014-08-13 北京理工大学 基于映射法的散乱点云Delaunay三角剖分曲面重构方法
CN105046688A (zh) * 2015-06-23 2015-11-11 北京工业大学 一种三维点云中的多平面自动识别方法
CN109993783A (zh) * 2019-03-25 2019-07-09 北京航空航天大学 一种面向复杂三维建筑物点云的屋顶及侧面优化重建方法
CN111563965A (zh) * 2019-02-14 2020-08-21 贝壳技术有限公司 一种通过优化深度图生成全景图的方法及装置
CN111709981A (zh) * 2020-06-22 2020-09-25 高小翎 特征线融合的激光点云与模拟图像的配准方法
CN112070870A (zh) * 2020-07-31 2020-12-11 广州景骐科技有限公司 点云地图评估方法、装置、计算机设备和存储介质

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103985155A (zh) * 2014-05-14 2014-08-13 北京理工大学 基于映射法的散乱点云Delaunay三角剖分曲面重构方法
CN105046688A (zh) * 2015-06-23 2015-11-11 北京工业大学 一种三维点云中的多平面自动识别方法
CN111563965A (zh) * 2019-02-14 2020-08-21 贝壳技术有限公司 一种通过优化深度图生成全景图的方法及装置
CN109993783A (zh) * 2019-03-25 2019-07-09 北京航空航天大学 一种面向复杂三维建筑物点云的屋顶及侧面优化重建方法
CN111709981A (zh) * 2020-06-22 2020-09-25 高小翎 特征线融合的激光点云与模拟图像的配准方法
CN112070870A (zh) * 2020-07-31 2020-12-11 广州景骐科技有限公司 点云地图评估方法、装置、计算机设备和存储介质

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
丁鸽 等: "基于离群点探测准则和主成分分析的点云平面拟合效果研究", 测绘地理信息, vol. 42, no. 3, 28 April 2017 (2017-04-28), pages 25 - 28 *
马学磊 等: "基于改进区域生长法的羊体点云分割及体尺参数测量", 中国农业大学学报, vol. 25, no. 3, 15 March 2020 (2020-03-15), pages 99 - 105 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116935231A (zh) * 2023-09-14 2023-10-24 湖北工业大学 一种隧道围岩结构面信息提取及关键块体识别方法
CN116935231B (zh) * 2023-09-14 2023-11-28 湖北工业大学 一种隧道围岩结构面信息提取及关键块体识别方法

Similar Documents

Publication Publication Date Title
CN111932688B (zh) 一种基于三维点云的室内平面要素提取方法、系统及设备
Bauchet et al. Kinetic shape reconstruction
Zhang et al. Online structure analysis for real-time indoor scene reconstruction
Vo et al. Octree-based region growing for point cloud segmentation
Shi et al. Adaptive simplification of point cloud using k-means clustering
CN105993034B (zh) 用于增强表面重构的轮廓求全
Vosselman et al. Recognising structure in laser scanner point clouds
Li et al. Feature-preserving 3D mesh simplification for urban buildings
Xu et al. Reconstruction of scaffolds from a photogrammetric point cloud of construction sites using a novel 3D local feature descriptor
Liu et al. Automatic segmentation of unorganized noisy point clouds based on the Gaussian map
Li et al. Modelling of buildings from aerial LiDAR point clouds using TINs and label maps
Lee et al. Perceptual organization of 3D surface points
CN115661374B (zh) 一种基于空间划分和模型体素化的快速检索方法
Zhang et al. Resolving topology ambiguity for multiple-material domains
Magillo et al. Algorithms for parallel terrain modelling and visualisation
Bauchet et al. City reconstruction from airborne LiDAR: A computational geometry approach
Qiu et al. An adaptive down-sampling method of laser scan data for scan-to-BIM
JP2978786B2 (ja) 物体の構造性立体幾何表示を作る方法及び装置並びに物体の構造性立体幾何表示を局所的に最適化する方法
US20220261512A1 (en) Segmenting a 3d modeled object representing a mechanical part
US20230281350A1 (en) A Computer Implemented Method of Generating a Parametric Structural Design Model
CN113763563A (zh) 一种基于平面识别的三维点云几何网格结构生成方法
Thiemann et al. 3D-symbolization using adaptive templates
US20180113958A1 (en) Method for Immediate Boolean Operations Using Geometric Facets
CN113673005A (zh) 在拓扑优化期间分配材料时在实体模型中保持形状的方法
KR100340080B1 (ko) 정렬되지 않은 3차원 거리 데이터로부터 캐드모델 생성 방법

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