CN113705582A - 一种建筑立面边缘特征关键点提取方法 - Google Patents

Abstract

本发明公开了一种建筑立面边缘特征关键点提取方法，属于三维建模领域。本发明首先通过对局部邻域内建立的协方差矩阵的特征值来计算三维点云中每个点的置信度；然后定义三维空间中每个点的梯度，再将三维点云的梯度信息编码为结构张量；结构张量的特征值可以表示局部点云的梯度分布情况，从而将建筑物点云的立面特征提取问题转换为分析结构张量的三个特征值的问题；将每个点的置信度和结构张量作为双阈值法的输入，判断当前被处理点是否为关键点；最后采用边缘细化算法对边缘特征点进行细化处理。本发明所提出的算法不仅精度高于现有直接从三维中提取边缘特征的方法，也优于基于Canny算子到二维图像进行边缘提取和细化的方法。

Description

一种建筑立面边缘特征关键点提取方法

技术领域

本发明属于三维建模领域，更具体地说，涉及一种建筑立面边缘特征关键点提取方法。

背景技术

边缘特征是物体几何特征的最显著的组成部分，作为最底层且最重要的特征边缘特征也是实现对地物场景理解、语义标注和三维抽象感知的基础。在计算机视觉的发展初阶阶段，绝大多数地物边缘检测算法均基于二维图像，毕竟在保证分辨率的基础上，图像的“线状”特征更加精确。另外在二维图像中地物的边缘特征具有明确的定义：是图像中特性(如像素灰度、纹理等)分布的不连续处，图像周围特性有阶跃变化或屋脊状变化的那些像素集合，也是图像信息最集中的地方。地物边缘是图像分割、图像理解及图像识别的重要特征。然后由于基于图像的边缘探测不能够充分反应复杂地物的三维几何，在基于三维的场景认知、几何表达等方面限制其进一步的发展。

近二十年来随着激光雷达传感器技术的提升，点云密度和精度显著提高，激光雷达点云已经能够充刻画物体的全局结构特征和局部细节几何，这使得从三维点云中精细提取出物体的边缘特征成为可能。即便如此，由于地物场景千变万化且扫描获取的点云中仍然会存在大量的噪声、离值点和由于地物遮挡和自遮挡带来的大量的点云数据的缺失，使得对地物边界特征的智能提取带来了巨大的挑战。另外本申请的发明人也注意到在三维点云中，对地物边缘的界定亦尚不明确，缺乏统一的定义。相比传统的图像特征提取，基于点云的地物特征的提取仍亟待完善。结合现有文献中对于三维点云中数据边缘的讨论，以及图像处理中对于边缘的定义，本发明确将三维点云边缘定义为以下三类点的集合：角点：三个相互不平行平面的交点；边界线：两个平面相交形成的交线，如建筑物立面与立面、建筑物立面与屋顶形成的边界线；边缘点：点云物体的外边界、物体内孔的边界(如：建筑物内窗户的窗框)或者是因为数据缺失引起的边界。

激光雷达点云的边缘特征提取是许多应用的基础，如三维重建、配准、目标检测、数据简化等。根据是否要借助于二维图像中的边缘特征提取方法，可以将现有的三维激光雷达点云中的边缘特征提取方法分为“基于图像的轮廓提取”和“激光雷达点云的轮廓提取”两类。

基于图像的轮廓提取：二维图像上的边缘检测算法比三维上的边缘检测算法成熟，因此很多方法都是将三维点云数据转换到二维图像上，利用二维图像上的边缘特征检测与提取方法，来对三维边缘特征进行处理。有的方法是将三维点云数据转换为二维灰度图像，如论文名称为：《从激光雷达数据和高分辨率图像中精确提取建筑物边界的新方法》，论文作者为：李慧、钟诚等，发表在《传感器检查》2013年3月，第33卷第2期，在该论文中首先根据三维点云生成对应的二维灰度图，将在二维灰度图中探测的边缘与三维点云配准以确定边缘点，该方法只针对建筑物屋顶轮廓提取；也有的方法是将三维点云转换为二值图，如论文名称为：点云数据自动建模框架，发表在《模式分析与机器智能学报》2013年3月，第35卷第11期，在该论文中首先将激光雷达三维点云中的建筑物的屋顶点转换为二值图像，然后用二值图像中的边缘检测方法去提取屋顶轮廓；还有的方法是利用三维点云对应的二维灰度图和距离图像，如论文名称：《融合三维点云和二维图像数据的边缘提取》，发表在《智能世界新兴技术国际会议暨博览会(CEWIT)》2013年10月，它是在三维点云对应的二维灰度图和距离图像中检测边缘，合并多组边缘点得到最终的三维点云边缘数据，该方法可以提取建筑物的边界并规则化。“基于图像的轮廓提取”中还存在着无需将三维点云向二维图像转化的方法，但需要将图像与三维点云数据进行配准以确定精确的特征对应关系。例如，论文名称：《基于光学影像与激光雷达数据融合的改进建筑物边界提取算法》，发表在《光学》2013年11月，第124卷第22期，在该论文中首先在激光雷达点云中用过滤、建筑物检测、墙体点去除和屋顶斑块检测四个步骤来检测出屋顶粗糙边缘，然后把这些粗糙边缘投射到图像空间中进行精细边缘的提取，该方法在复杂环境下具有很强的鲁棒性；论文名称：《提出了一种基于激光雷达数据和图像的建筑物边界重建方法》，发表在《激光传感与成像及其应用》2013年9月，第8905卷第11期，在该论文中将单景的三维点云数据与单幅图像进行匹配，在对应图像中进行边缘检测，将它与三维点云相对应来获取三维点云的边缘，该方法计算复杂度低。

激光雷达点云的轮廓提取：三维点云数据包含物体结构信息与拓扑关系，因此也可以利用三维点云数据固有的曲率、法向量、梯度等信息来提取点云边缘特征。由于边缘点的曲率高于非边缘区域点的曲率，论文名称：《利用曲线自动配准密集地面激光扫描点云》，发表在《摄影测量与遥感学报》2014年9月，第95卷第10期，在该论文中使用局部曲率作为边缘指标，选择曲率高于给定阈值的点作为边缘,该方法适用于高质量点云中的小型物体的边缘提取。三维点云包含的信息中也可以获取除曲率之外的特征，如论文名称：《基于正态估计和图论的点云闭合锐边检测》，发表在《计算机辅助设计》2007年4月，第39卷第4期，在该论文中首先基于法向量估计来对边缘特征点进行提取，并将提取出的特征点处理为一个图结构来恢复尖锐的特征线，该算法速度快，且在工业应用中获得了较好的实际效果。在三维点云中，大多数边缘也可以看作是表面的交点。论文名称：《基于视觉的定位，使用从3D激光距离数据提取的边缘地图》，发表在《2010年机器人与自动化国际会议》2010年5月，在该论文中首先将点云分割成平面，然后通过计算相邻平面之间的相交线来定位边缘，但这种方法并不适用于提取小型物体的边缘。为了解决不能提取小型物体边缘的问题，论文名称：《基于邻域几何特性分析的三维点云边缘检测与特征线跟踪》，发表在《遥感》2016年9月，第8卷第9期，在该论文中将每个点的邻域点投影到局部平面上，然后根据邻域点在平面上的角度间隙来提取边缘点,这个方法能在小场景的高质量点云中取得良好的效果，但是对于大型复杂场景，涉及的参数调整将会比较麻烦。

参考上述对“基于图像的轮廓提取”和“激光雷达点云的轮廓提取”的描述可以发现，“基于图像的轮廓提取”是将三维点云数据转换到二维图像上，利用了较为成熟的二维图像边缘检测算法来提取各类边缘特征，但“基于图像的轮廓提取”忽略了三维点云的空间信息与拓扑关系，且在三维与二维转换过程中产生了很多的信息丢失，会导致结果精度偏低；而“激光雷达点云的轮廓提取”虽然利用了三维物体的空间信息与拓扑关系，但该类方法大多数只适用于特定物体(如：建筑物、道路等)的边缘检测，也要对原始的激光雷达点云数据进行分割、目标识别等步骤且很难控制所提取出的边缘的粗细程度。

发明内容

针对现有基于图像的轮廓提取法会导致信息丢失，结果精度偏低以及激光雷达点云的轮廓提取法只适用于特定物体的边缘检测且很难控制所提取出的边缘的粗细程度等问题，本发明提供一种建筑立面边缘特征关键点提取方法，能够从三维点云数据中精确提取建筑外立面的边缘特征关键点，并且能够自由控制边缘关键点细化的粒度，以增强建筑立面抽象表达的灵活性和紧凑型，同时克服当前激光雷达点云关键点提取大多需要实施“面片分割”、“部件识别”等操作所带来的效率低下等缺陷。

为解决上述问题，本发明采用如下的技术方案。

一种建筑立面边缘特征关键点提取方法，包括如下步骤：

步骤一：获取建筑立面的三维点云数据；

步骤二：由三维点云数据中的点p_i和其邻域球N内的点构建协方差矩阵，并根据协方差矩阵的特征值来逐点计算每个点的置信度Cⁱ，邻域球N的球心为p_i；

步骤三：计算三维点云数据中每个点的梯度G(x,y,z)，并将梯度G(x,y,z)沿空间坐标系中三个坐标轴的方向分解为三个分量；

步骤四：根据梯度G(x,y,z)的三个分量构建梯度G(x,y,z)的结构张量M，并对结构张量M采用高斯权函数进行平滑处理；

步骤五：结合每个点的置信度Cⁱ、梯度G(x,y,z)以及平滑处理后结构张量M，采用双阈值法判别当前点是否为初始关键点，并将所得到的所有初始关键点提取出来从而得到初始关键点集；

步骤六：采用边缘细化算法对步骤五中得到的初始关键点集进行细化。

进一步的，步骤二中，每个点p_i的置信度Cⁱ通过拟合质量C_f和均匀分布性C_s来表示：

Cⁱ∈[0，1]，其中，n表示每个点置信度的静态尺度的数目，w表示每个点置信度的静态尺度的计数数目；

拟合质量C_f表示局部切平面在点p_i处的拟合质量，计算公式为：

均匀分布性C_s为局部均匀采样性，计算公式为：

其中，

和

表示协方差矩阵的三个特征值；p_i的邻域点p_j的坐标与邻域球N的几何中心的坐标相减得到3×1列向量，3×1列向量与其转置后的1×3行向量相乘得到3×3矩阵，将所有矩阵相加得到协方差矩阵；i、j分别表示点p_i和邻域点p_j的下标。

进一步的，邻域球N的半径为三维点云数据的平均密度的1.0倍、1.5倍或2.0倍。

进一步的，步骤三中，点p_i的梯度G(x，y，z)的计算如下：

定义函数f(x，y，z)为可微函数，三维空间中的方向向量u＝cos(α)A+cos(β)B+cos(γ)C，其中，A、B、C分别表示沿x、y、z坐标轴方向的单位向量，α，β，γ分别表示方向向量u与x、y、z坐标轴正方向的夹角，函数f(x，y，z)沿方向向量u的方向导数

的计算公式为：

其中t表示当前点p_i和邻域点p_j的距离；

在离散点云中，可微函数f(x，y，z)表示点p_i处的置信度，点p_i的梯度G(x，y，z)与方向导数

的关系公式如下：

其中，Cⁱ和C^j分别表示当前点p_i和邻域点p_j的置信度，d_ij表示点p_i和点p_j的欧几里得距离；

点p_i处的梯度G(x，y，z)沿x、y、z坐标轴方向的分解为：

其中，

分别表示点p_i和点p_j在x、y、z坐标轴方向上的坐标差，g_x、g_y和g_z分别表示当前点的梯度在x、y、z坐标轴三个方向上的分解值。

进一步的，步骤四中，结构张量M的具体构建步骤为：

在三维点云空间中，点p_i在邻域球N内的梯度变化E(Δx，Δy，Δz)通过向不同方向(Δx，Δy，Δz)的移动位移来计算，梯度变化E(Δx，Δy，Δz)的公式为：

上述公式的泰勒一阶展开式为：E(Δx，Δy，Δz)＝∑_{p(x，y，z)∈N}[Δx·g_x+Δy·g_y+Δz·g_z]²＝(Δx，Δy，Δz)M(Δx，Δy，Δz)^T；

其中，Δx表示在x轴方向上的坐标移动位移，Δy表示在y轴方向上的坐标移动位移，Δz表示在z轴方向上的坐标移动位移；M为以点p_i为中心的邻域球N的结构张量；

M为半正定对称矩阵，其三个特征值

和

两两相互正交。

进一步的，对结构张量M使用高斯权函数进行平滑处理，平滑处理如下公式：

其中，符号

代表点p_i生成的结构张量M_i和邻域点p_j生成的结构张量M_j进行卷积操作，σ是点p_i和其所有邻域点p_j之间距离的标准差，

为高斯平滑的结果。

进一步的，步骤五中，采用结构张量M的三个特征值

和

险测建筑立面中的初始关键点，

初始关键点的响应函数公式为：R_M＝det(M)-k[trace(M)]³，参数k取值为1；

若R_M≥T_M且G(x，y，z)≥T_G，则保留该初始关键点，否则删除该初始关键点；其中，梯度G(x，y，z)为正态分布，G(x，y，z)～N(ψ，ω)，ψ和ω分别是三维点云数据中所有点云的梯度G(x，y，z)的平均值和标准差；T_G为双阈值法中预定义的梯度阈值，T_G∈[0.6ψ，ψ]；T_M为双阈值法中预定义的阈值，T_M∈[5，30]。

进一步的，步骤六中，边缘细化算法采用网格简化法。

进一步的，步骤六中，边缘细化算法采用层次简化法。

进一步的，步骤六中，边缘细化算法采用加权局部最优投影算法。

相比于现有技术，本发明的有益效果为：

(1)本发明的一种建筑立面边缘特征关键点提取方法，能够自由控制所提建筑立面关键点的粒度，增强了建筑立面抽象表达的灵活性和紧凑度。

(2)本发明的一种建筑立面边缘特征关键点提取方法，不仅精度高于现有直接从三维中提取边缘特征的方法，也优于基于Canny算子到二维图像进行边缘提取和细化的方法。

(3)本发明的一种建筑立面边缘特征关键点提取方法，建筑立面关键点提取算法对点云噪声不敏感、能够处理各类泛在点云，包括：机载、车载、地面、手持式SLAM点云和无人机获取的倾斜摄影点云。

(4)本发明的一种建筑立面边缘特征关键点提取方法，建筑立面关键点提取算法无需要对建筑面片进行分割和部件探测等操作，高效较高。

附图说明

图1为本发明的流程图。

具体实施方式

下面结合具体实施例和附图对本发明进一步进行描述。

实施例1

本实施提供了一种建筑立面边缘特征关键点提取方法，算法流程图如图1所示，先输入建筑立面的三维数据点云，为了生成建筑物立面边缘特征关键点，具体包括如下步骤：

S1、逐点计算三维点云的置信度：通过对局部邻域内建立的协方差矩阵的特征值分析来估计三维点云中每个点的置信度；

置信度是指以单个点为中心、以固定半径形成的邻域球内的点云局部质量。由于本发明的目的是从建筑物立面中提取边缘特征点，置信度的定义应该与这些关键点的特征直接相关。基于这个原则，本申请的发明人采用两个几何性质的组合来估计每个点的置信度：拟合质量C_f和均匀分布性C_s。拟合质量C_f表示局部切平面在点p_i处的拟合质量，计算方法为：

其中

点p_i及其邻域球N内所有的邻域点p_j构成协方差矩阵，

和

表示协方差矩阵的三个特征值。协方差矩阵的构成：p_i的邻域点p_j的坐标与邻域球N的几何中心的坐标相减得到3×1列向量，3×1列向量与其转置后的1×3行向量相乘得到3×3矩阵，协方差矩阵由所有邻域点计算出的矩阵相加得到。如果点p_i和它的邻域点p_j能够完全拟合局部切平面，则点p_i的拟合质量

越接近0。相反，如果点p_i的邻域是不均匀分布的，它们很可能是分布在建筑物立面边缘点、立面角点等。在这种情况下，p_i点的拟合质量

会接近1。均匀分布性C_s为局部均匀采样性，它可以由如下公式被计算出来：

如果p_i点和其邻域点p_j呈线状分布，则p_i点的C_s值会接近0；如果均匀分布，则Cs值会接近1。因此，该方法是对建筑物立面的外边界和窗户边框的有效检测，在这些外边界中，由于激光束穿透窗玻璃，经常会出现点云密度变化大和数据丢失的情况。另外，为了使置信度的计算更加稳健，通过改变以当前点为中心的邻域球的半径大小，在多个尺度上计算上述两个几何属性。比如，通过设置三个邻域球N的尺度，半径分别为三维点云数据的平均密度的1.0、1.5和2.0倍。所以，点p_i完整的置信度Cⁱ∈[0，1]由公式(1)给出：

其中，参数n表示每个点置信度的静态尺度的数目，w表示每个点置信度的静态尺度的计数数目。很明显的是，如果每个点的置信度Cⁱ接近1，这意味着p_i点很有可能来自建筑物立面边界点、角点、窗户边框和其他形状变化很大的地方。相反，如果某个点的置信度接近0，p_i点具有较高的局部拟合质量，则这个点很有可能来自建筑物立面的内部点。

S2、定义数据空间中三维点云的梯度：计算三维空间中每个点的梯度G(x，y，z)，将梯度G(x，y，z)沿三个坐标轴方向分解为三个分量；

一旦定义了点云的置信度，就可以根据置信度来定义三维空间离散点云的梯度。首先分析方向导数，并由此可以推导出梯度。方向导数表示函数在任意给定方向上的变化率。梯度的大小是当前点所有方向导数中最大的值，梯度方向是与方向导数相对应的方向。让函数f(x，y，z)代表有三个自变量的可微函数，u＝cos(α)A+cos(β)B+cos(γ)C代表三维空间中的方向向量u，其中，A、B、C分别表示沿x、y、z坐标轴方向的单位向量。函数f(x，y，z)沿方向向量u的方向导数

计算如公式(2)：

其中α，β，γ分别代表方向u与x、y、z坐标轴正方向的夹角，t代表当前点p_i和邻域点p_j的距离。当梯度

方向和方向向量u的方向相同时，方向导数的值达到最大。

在离散点云中，函数f(x，y，z)代表p_i点处的置信度，点p_i处的梯度值G(x，y，z)可以通过最大化局部区域内的方向导数得到，它近似等于

在三维空间离散形式下，梯度与方向导数的关系如公式(3)所示：

点p_i处的梯度值G(x，y，z)可以沿三个坐标轴方向分解为如公式(4)：

其中，i、j分别表示当前点p_i和它的邻域点p_j的下标，

分别代表点p_i和邻域点p_j在x、y、z三个坐标轴方向上的坐标差，N是点p_i的邻域球，Cⁱ、C^j分别表示当前点p_i和其邻域点p_j的置信度，d_ij代表点p_i和邻域点p_j的欧几里得距离，g_x、g_y和g_z分别表示当前点的梯度在x、y、z坐标轴三个方向上的分解值。

S3、构建结构张量：将三维点云的梯度信息编码为一个3×3结构张量，结构张量的特征值可以反映每个点局部邻域的梯度分布；

通过在三维点云空间中向不同方向(Δx，Δy，Δz)少量移动位移来计算梯度变化。点p(x，y，z)在其邻域球N内的梯度变化E(Δx，Δy，Δz)可以被表达为公式(5)：

E(Δx，Δy，Δz)＝∑_{p(x，y，z)∈N}[G_{x+Δx，y+Δy，z+Δz}-G_x，y，z]² (5)

公式(5)的泰勒一阶展开为公式(6)：

E(Δx，Δy，Δz)＝∑_{p(x，y，z)∈N}[Δx·g_x+Δy·g_y+Δz·g_z]²≈(Δx，Δy，Δz)M(Δx，Δy，Δz)^T (6)

其中，Δx表示在x轴方向上的坐标移动位移，Δy表示在y轴方向上的坐标移动位移，Δz表示在z轴方向上的坐标移动位移；

，M表示以点p_i为中心的邻域球N的结构张量。M是一个半正定对称矩阵，它的三个特征值

和

是两两相互正交的。邻域球N内的梯度变化可以通过分析结构张量M的三个特征值来得到。

当前点p_i存在着如下四种情况：

角点：当前点p_i位于三个互不平行的平面交点处，此时，结构张量M的三个特征值都比较大；

边缘点：当前点p_i可能位于建筑物立面与立面或建筑物立面与屋顶相交形成的边缘点，在这种情况下，结构张量M的三个特征值中有两个特征值相对较大；

边界点：当前点p_i可能来源于建筑物的外边界、内孔边界(窗户边框)或者由数据缺失引起的“假”边界。在这种情况下，结构张量M的三个特征值中只有一个特征值较大。

连续点：当前点p_i的邻域内保持着近似恒定的梯度值，即以p_i为中心的三维体素窗口移动时，E的变化值很小。这种情况下，结构张量M的三个梯度值都较小。

在本发明中，关键点是指角点、边缘点和边界点，这三部分构成了建筑物立面的骨架。结构张量M作为梯度矩阵，不是连续表示，而是以离散形式计算，计算出的特征值是不稳定和不鲁棒的。也就是说，得到的结构张量M不能完全反映局部邻域的梯度分布。为解决这一缺陷，需要对结构张量M使用高斯权函数进行平滑处理，平滑处理如公式(7)：

其中，运算符号“·”代表数字与矩阵相乘，运算符号

代表点p_i生成的结构张量M_i和邻域点p_j生成的结构张量M_j进行卷积操作，σ是点p_i和其所有邻域点p_j之间距离的标准差，

为高斯平滑的结果。

S4、利用双阈值法判别初始关键点：结每个点p_i的置信度Cⁱ、梯度G以及结构张量M，采用双阈值法判别当前被处理点是否为初始关键点，并将所得到的所有初始关键点提取出来从而得到初始关键点集；

由于特征值的“大”或“小”很难去定义，于是利用结构张量M的三个特征值来定义响应函数，从而来检测建筑立面中的关键点。关键点响应函数定义如公式(8)：

R_M＝det(M)-k[trace(M)]³ (8)

R_M值可以表达当前点所形成的结构张量M的三个特征值组合情况，k是可以调整的参数，通常设为1。也就是说，如果R_M的值大于或等于双阈值法中预定义的阈值TM，则将其对应的点视为关键点。但是，设置不恰当的阈值T_M会导致建筑物立面的关键点分割不足或过分割。这很可能是因为R_M的计算值是一个局部变量，它是通过置信度、梯度张量生成、张量平滑等多个步骤确定的，这些步骤使用了不同的局部邻域大小的尺度。相比之下，阈值T_M是在整个建筑立面数据空间中定义的全局阈值，这导致固定的阈值不能完全代表在不同的局部邻域中估计的每个点的梯度变化。

通过实验本申请的发明人发现，如果只用R_M≥T_M这一个条件来判断当前点是否为关键点的话，所提取出的建筑物骨架较厚且屋顶上的许多杂点也会被错误提取出来。为了缓解这个问题，本申请的发明人通过引入另一个约束条件，即：G(x，y，z)≥T_G来实现双阈值提取关键点，其中TG是预定义的梯度阈值。也就是说，如果一个点被标记为关键点，则它同时满足两个条件：R_M≥T_M并且G(x，y，z)≥T_G，符合上述两个条件的初始关键点保留，不符合的删除，由所有的初始关键点构成初始关键点集。

实际上，梯度符合一个正态分布，即G(x，y，z)～N(ψ，ω)，ψ和ω分别是三维点云数据中所有点云的梯度G(x，y，z)的平均值和标准差；它们可以通过梯度分布的直方图得到。TG的阈值设置建议取小一点，大概在[0.6ψ，ψ]之间。阈值T_M控制了判断一个点是否为关键点的严格程度，选择小或大的值将导致欠分割或过分割。然而，阈值T_M对大范围的值不太敏感，建议设置在[5，30]这个范围内。与单阈值所得的结果相比，双阈值法得到的关键点结果相对较薄，更能合理地描述建筑物立面和屋顶的骨架。此外，屋顶上一些零散的异常点也被消除，结果更加干净简洁。

S5、细化初始关键点集：为了进一步提升初始关键点提取的精度，采用了边缘细化算法进一步细化初始关键点集，提升建筑立面关键点表达的紧凑度，增强建筑立面抽象表达的灵活性。

虽然也可以通过联合调优两个阈值，即T_M和T_G来获得一个最优的初始关键点集，以保证建筑物立面表示的几何精度和高紧致性，但需要注意的是，同时调优两个参数并不是一件简单的工作。它需要尝试和错误的实验，导致耗时高和劳动密集。从这个意义上说，本申请倾向于获得一个相对较厚的建筑物骨架，边缘上有较厚的点，并在后续的简化步骤中对结果进行细化。为此，本申请的发明人将建筑物骨架抽象问题转化为简化问题。也就是说，本申请的发明人通过简化的概念进一步细化双阈值法所得结果，以增强建筑物立面抽象的灵活性。

为此，本申请采用了现有的三种经典算法中的任意一种算法对初始关键点集进行细化，三种经典算法为网格简化法、层次简化法和加权局部最优投影算法(WeightedLocally Optimal Projection,WLOP)。

网格简化算法将点云空间划分为多个超体素，从中选取具有代表性的一个点云来表示对应的超体素。

层次简化法通过使用二进制空间划分递归分割点云，来提供初始关键点集的自适应简化结果。

本实施例中选用了加权局部最优投影算法，与其他两种简化算法相比，加权局部最优投影算法不仅简化了建筑物立面的关键点云，而且对其进行了正则化。这意味着得到的点云是新生成的，而不是简单地从关键点集中选择得到的。

实施例2

本实施例中提供了对实施例提出的一种建筑立面边缘特征关键点提取方法的具体实验数据。

第一个实验数据集是Semantic3D，它提供了一个大规模的超过40亿标记点云的不同户外场景。本申请的发明人选取了8个具有代表性的建筑来验证关键点提取算法。从Semantic3D中选择建筑的原因有三个：

第一、精确标注的点云：手动为Semantic3D中的每个点云指定一个特定的类标签，可以很容易地根据建筑标签选择不同的建筑来验证所本发明提出的算法。

第二、多样的建筑：所有发布的场景都是在中欧拍摄的，从中提供了多样的欧洲建筑。

第三、非均匀扫描点：该数据集具有特殊的挑战性，因为其扫描是使用测量级陆地激光扫描仪(Terrestrial Laser Scanning,TLS)获取的，测量范围较长，因此会造成极值点密度变化和遮挡。

正因为如此，包含的建筑非常适合测试本发明所提出的算法。

为了评估本发明提出的算法处理极稀疏点云下的建筑立面的能力，本申请的发明人从荷兰AHN3(Actueel Hoogtebestand Nederland)和都柏林城市数据集中选择了两组建筑。荷兰数据集在2014年和2019年由一架飞机获得。在该数据集的采集中，Riegl LMS-Q680i激光扫描传感器是最常用的，有时也采用Riegl VQ-780i。点云的平均密度约为16个点/平方米。AHN3数据集通过中央分布式平台PDOK向公众发布。都柏林数据集于2015年3月通过直升机来采集的。整个数据集由14条飞行路径组成，覆盖都柏林市中心主要区域约5.6平方公里。采用的传感器是TopEye system S/N 443。整个数据集由500×500的矩形贴图的形式发布。多重路径配准后的点云提供了250-348点/平方米的平均点密度。这两个ALS(Airborne Laser Scanning,ALS)数据集具有清晰的语义建筑标识，便于选择具有不同几何形状的建筑。因为从ALS数据集中选择，这两组数据中的建筑物立面点云是非常稀疏的，并且包含了由于建筑自遮挡而缺失的数据。此外，所选的两组数据集在不同的方向上有不同的尺寸和屋顶形状/结构。这两个数据集提供了一个测试所提出算法的抽象能力的机会，不仅关注建筑物立面，而且关注屋顶。

为了评估算法的可扩展性，从TerraMobilita/iQmulus 3D城市分析基准中选择了一个区域的MLS点云(Mobile Laser Scanning,MLS)。整个基准包含11个区域，在法国巴黎中心有3亿点云。该数据集于2013年1月由法国国家测绘局开发的MLS系统Stereopolis II获得。两个Riegl LMS-Q12Oi和一个HDL-64E Velodyne激光雷达传感器集成到StereopolisII中。选择的区域包含一个完全注释的200米长的街道部分，有1200万个点云。可以认为这部分来自巴黎密集城市环境的MLS点云可以很好地大规模测试算法的可扩展性。

利用无人机摄影测量转换所得点云和手持式扫描点云验证了算法的鲁棒性。利用DJI Phantom 4RTK无人机(Unmanned Aerial Vehicle,UAV)获取了南京林业大学A楼和B楼影像数据。DJI Phantom 4RTK飞行高度为100米，地面采样距离(Ground SamplingDistance,GSD)为2.74cm，航向和旁向的重叠率分别设置为80％和70％。将采集到的倾斜图像输入Bentley ContextCapture商用软件中生成高密度点云。由于密度过于高，本申请将原始点云降采样到3cm。本发明还采用手持式激光雷达扫描仪(Hand-holding LaserScanning,HLS)GeoSLAM ZEB-HORIZON Scanner扫描南京林业大学C楼，基于SLAM(Simultaneous Localization and Mapping,SLAM)技术实时实现三维点云配准和拼接。这两种类型的点云具有点精度低、异常值和噪声高的特点，因此对关键点的提取造成了极大的挑战。

算法参数设置：

算法在实施过程中，针对各类实验数据，算法所用参数设置如表1所示：

表1针对不同实验数据算法参数的具体设置

实验结果-紧凑性评价：

紧凑性是评价所提取的关键点集是否轻量化的关键指标。也就是说，如果结果包含更少数目的关键点，那么它代表着建筑物立面抽象更加紧凑。高度紧凑的关键点集便于网络传输、存储保存和大规模建筑点云渲染。本申请用生成关键点的数量除以原始建筑点云的数量计算出的压缩比C_a来表示紧凑程度。此外，本申请还使用另一种压缩比C_b，它是将创建的关键点的数目除以参考结果点云的数目。参考结果点云是通过使用CloudCompare开源工具进行大量手工制作生成的。

不同数据集选取的建筑定量评价结果如表2所示，可以看到，对于8个Semantic3D建筑，双阈值结果的平均紧凑度接近36.96％。此时关键点的数目与参考点的数目比较一致，如双阈值列C_b所示。这一比例可以根据不同的抽象程度，通过三种细化算法进一步压缩到20％以下，甚至10％以下。值得注意的是，高的紧凑性(例如网格细化后的语义建筑8的紧凑性为0.92％)并不意味着结果更合理。虽然具有很高的紧凑度，但在某些情况下，细节形状特征可能会被忽略，从而影响抽象精度。对于其他两个ALS和MLS数据集，比值的紧凑度C_a比TLS数据集要低得多(例如，经过层次细化后的都柏林数据集只有0.72)。这是因为ALS点云由于点密度稀疏，没有足够的能力来描述建筑形状的细节。由于高密度的点云使得UAV和HLS数据集包含了足够的形状细节，当进行较大的抽象(经过网格或层次细化后小于1％)时，可以保证提取精度。综上所述，紧凑性与场景复杂度、点云密度、抽象程度密切相关。为了保证关键点的抽象既紧凑又准确，必须在两者之间取得平衡。

表2不同数据类所提关键点型紧凑度指标

实验结果-精度评价：

本申请的发明人在Semantic3D数据集中为每个建筑提供了关键点云的参考点集。该参考点集是通过使用CloudCompare开源工具手动标注生成的。在此基础上，本申请的发明人评估了双阈值法和网格简化法、分层简化法和WLOP算法所产生的结果的质量。定量评价统计如表3所示。建筑物3和4的最大值大于3米，这是因为激光束穿透建筑物立面窗户玻璃产生了离群点，这些离群点被错误提取为关键点了。另一个原因是，结果中包含了一些伪边缘/边界，这些伪边缘/边界是出现在建筑物立面的由于遮挡而造成数据缺失的区域。事实上，建筑物3、4的关键点参考点云中并没有包含这些伪边缘，这也会造成最大值偏大。本申请的发明人还观察到，双阈值的最大值全都大于网格简化法和分层简化法的结果，因为这两个简化方法的结果是从初始双阈值结果中进行降采样得到的。虽然WLOP也是由初始的双阈值结果生成的，但是建筑物7的WLOP简化法后的最大值为1.3606m，大于对应的双阈值结果1.3306m。这可能是因为WLOP简化法不是简单地从初始的双阈值关键点集中进行降采样，而是通过优化技术重新生成点集。除了建筑物2和4的分层简化法均值为0.0271m和0.0442m外，双阈值法结果的均值都小于三种简化法结果的均值。这意味着通过简化的概念进行细化会削弱几何精度。均方根误差(RMSE)度量产生的关键点集和参考点集之间的差异。为了评估这种差异的重要性，并使这种差异独立于单元，本申请使用另一个相关度量相对均方根误差，即将均方根误差的值除以相应建筑的边界盒的对角线长度。从而可以看到，相对均方根误差的最大值为建筑物6的0.0088，这意味着本申请的结果和参考结果之间的最大差距只有建筑对角线的0.88％。

表3 TLS数据中8幢建筑的关键点提取精度统计

与现有方法的对比评价：

本申请的发明人将本发明提出的算法与其他两种经典算法(Canny边缘检测器和Xia的方法)进行了比较。由于Canny边缘检测器需要在二维图像上工作，而不是在三维点云上，需要根据建筑物立面的法向量和位置来将三维点云投影到二维灰度图像上。每个像素的值是当前像素内所有点云的梯度平均值。像素值被归一化为0到255的范围。与Canny边缘检测器在图像上的应用相比，Xia的方法和本申请的方法是直接在三维点云上进行处理的。对比两种方法的结果可以发现，本申请的方法可以很好地从形状高度复杂的立面点中提取关键点。对于复杂的形状，本申请的方法可以得到一个合理的建筑物骨架，但是在Xia的结果中缺少了一些边缘点。由于本申请的方法直接在三维点云上工作，就像Xia的方法一样，这意味着这两种方法可以很好地保持建筑物立面的几何形状而不需要人工干预。但对于Canny方法来说，在从三维点云到二维灰度图像的转换过程中不可避免地会带来失真。另外，需要注意的是，Canny检测器经常生成“双边”效应，即由于边缘的厚度/宽度，一条边被两条等高线表示。还应该注意到，本申请的方法对点云的密度也不敏感。

除了与Canny检测器和Xia的方法进行定性比较外，本申请中还利用Semantic3D数据集中的8个建筑物与Xia的方法进行了定量比较。首先根据关键点参考点云来计算每个建筑的平均值、标准差、均方根误差和紧凑度。接下来，对所有这些数据取平均值，从而得到整个Semantic3D数据集的整体值。统计数据如表4所示，本申请的方法在几何精度和紧致性方面都优于Xia的方法。

表4所提方法与Canny边缘检测器和Xia的方法定量对比结果

本发明所述实例仅仅是对本发明的优选实施方式进行描述，并非对本发明构思和范围进行限定，在不脱离本发明设计思想的前提下，本领域工程技术人员对本发明的技术方案作出的各种变形和改进，均应落入本发明的保护范围。

Claims (10)

1.一种建筑立面边缘特征关键点提取方法，其特征在于，包括如下步骤：

步骤一：获取建筑立面的三维点云数据；

步骤二：由三维点云数据中的点p_i和其邻域球N内的所有点构建协方差矩阵，并根据协方差矩阵的特征值来逐点计算每个点的置信度Cⁱ，邻域球N的球心为p_i；

步骤三：计算三维点云数据中每个点的梯度G(x，y，z)，并将梯度G(x，y，z)沿空间坐标系中三个坐标轴的方向分解为三个分量；

步骤四：根据梯度G(x，y，z)的三个分量构建结构张量M，并对结构张量M采用高斯权函数进行平滑处理；

步骤五：结合每个点的置信度Cⁱ、梯度G(x，y，z)以及平滑处理后结构张量M，采用双阈值法判别当前点是否为初始关键点，并将所得到的所有初始关键点提取出来从而得到初始关键点集；

步骤六：采用边缘细化算法对步骤五中得到的初始关键点集进行细化。

2.根据权利要求1所述的一种建筑立面边缘特征关键点提取方法，其特征在于：步骤二中，每个点p_i的置信度Cⁱ通过拟合质量C_f和均匀分布性C_s来表示：

其中，n表示每个点置信度的静态尺度的数目，w表示每个点置信度的静态尺度的计数数目；

拟合质量C_f表示局部切平面在点p_i处的拟合质量，计算公式为：

均匀分布性C_s为局部均匀采样性，计算公式为：

其中，

和

表示协方差矩阵的三个特征值；p_i的邻域点p_j的坐标与邻域球N的几何中心的坐标相减得到3×1列向量，3×1列向量与其转置后的1×3行向量相乘得到3×3矩阵，将所有矩阵相加得到协方差矩阵；i、j分别表示点p_i和邻域点p_j的下标。

3.根据权利要求2所述的一种建筑立面边缘特征关键点提取方法，其特征在于：邻域球N的半径为三维点云数据的平均密度的1.0倍、1.5倍或2.0倍。

4.根据权利要求2或3所述的一种建筑立面边缘特征关键点提取方法，其特征在于：步骤三中，点p_i的梯度G(x，y，z)的计算如下：

定义函数f(x，y，z)为可微函数，三维空间中的方向向量u＝cos(α)A+cos(β)B+cos(γ)C，其中，A、B、C分别表示沿x、y、z坐标轴方向的单位向量，α，β，γ分别表示方向向量u与x、y、z坐标轴正方向的夹角，函数f(x，y，z)沿方向向量u的方向导数

的计算公式为：

其中t表示当前点p_i和邻域点p_j的距离；

在离散点云中，可微函数f(x，y，z)表示点p_i处的置信度，点p_i的梯度G(x，y，z)与方向导数

的关系公式如下：

其中，Cⁱ和C^j分别表示当前点p_i和邻域点p_j的置信度，d_ij表示点p_i和点p_j的欧几里得距离；

点p_i处的梯度G(x，y，z)沿x、y、z坐标轴方向的分解为：

其中，

分别表示点p_i和点p_j在x、y、z坐标轴方向上的坐标差，g_x、g_y和g_z分别表示当前点的梯度在x、y、z坐标轴三个方向上的分解值。

5.根据权利要求4所述的一种建筑立面边缘特征关键点提取方法，其特征在于：步骤四中，结构张量M的具体构建步骤为：

在三维点云空间中，点p_i在邻域球N内的梯度变化E(Δx，Δy，Δz)通过向不同方向(Δx，Δy，Δz)的移动位移来计算，梯度变化E(Δx，Δy，Δz)的公式为：

上述公式的泰勒一阶展开式为：

其中，Δx表示在x轴方向上的坐标移动位移，Δy表示在y轴方向上的坐标移动位移，Δz表示在z轴方向上的坐标移动位移；M为以点p_i为中心的邻域球N的结构张量；

M为半正定对称矩阵，其三个特征值

和

两两相互正交。

6.根据权利要求5所述的一种建筑立面边缘特征关键点提取方法，其特征在于：对结构张量M使用高斯权函数进行平滑处理，平滑处理如下公式：

其中，符号

代表点p_i生成的结构张量M_i和邻域点p_j生成的结构张量M_j进行卷积操作，σ是点p_i和其所有邻域点p_j之间距离的标准差，

为高斯平滑的结果。

7.根据权利要求6所述的一种建筑立面边缘特征关键点提取方法，其特征在于：步骤五中，采用结构张量M的三个特征值

和

检测建筑立面中的初始关键点，

初始关键点的响应函数公式为：R_M＝det(M)-k[trace(M)]³，参数k取值为1；

若R_M≥T_M且G(x，y，z)≥T_G，则保留该初始关键点，否则删除该初始关键点；其中，梯度G(x，y，z)为正态分布，G(x，y，z)～N(ψ，ω)，ψ和ω分别是三维点云数据中所有点云的梯度G(x，y，z)的平均值和标准差；T_G为双阈值法中预定义的梯度阈值，T_G∈[0.6ψ，ψ]；T_M为双阈值法中预定义的阈值，T_M∈[5，30]。

8.根据权利要求7所述的一种建筑立面边缘特征关键点提取方法，其特征在于：步骤六中，边缘细化算法采用网格简化法。

9.根据权利要求7所述的一种建筑立面边缘特征关键点提取方法，其特征在于：步骤六中，边缘细化算法采用层次简化法。

10.根据权利要求7所述的一种建筑立面边缘特征关键点提取方法，其特征在于：步骤六中，边缘细化算法采用加权局部最优投影算法。

Images