论文推荐|朱杰：多约束的平面点集形状重构方法

《测绘学报》

构建与学术的桥梁        拉近与权威的距离

多约束的平面点集形状重构方法

朱杰^1,2, 孙毅中^1,2     

1. 南京师范大学虚拟地理环境教育部重点实验室, 江苏 南京 210023; 
2. 江苏省地理信息资源开发与利用协同创新中心, 江苏 南京 210023

收稿日期：2016-03-30； 修回日期：2016-10-27

基金项目：国家自然科学基金（41671392）；公安部科技强警基础工作专项（2015GABJC39）

第一作者简介： 朱杰(1989-),男,博士生,研究方向为城市空间数据表达与挖掘。E-mail： 

通信作者： 孙毅中

E-mail:

摘要：针对平面点集空间分布的复杂性，本文提出了一种基于Delaunay三角网的平面点集形状重构方法。首先采用一种简单且实用的数据结构以表达Delaunay三角网中嵌入的几何信息和拓扑信息，然后由外向内迭代过滤Delaunay三角网得到一个大概边界，最后进一步考虑边界的凹凸信息和空洞现象，获取最终的精细边界。试验结果表明与其他典型的Delaunay三角网重构方法相比，本文提出的算法能更好地适用于平面点集空间分布的复杂性，通过所构建的数学模型实现了凸凹多边形内外边界提取。

关键词： 平面点集     形状重构     Delaunay三角网     多约束     GIS    

An Efficient Approach to Shape Reconstruction from Planar Point Set Based on Multi-constraints

ZHU Jie^1,2, SUN Yizhong^1,2     

Abstract: An efficient algorithm to boundary representation from a planar point set in order to adapt the complexity of spatial distribution was presented in this paper. At first, an appropriate and practical data structure was designed to express geometric information and topological information, which provides an easy access to links embedded in DT serving as a basis for the filtering procedures; then the algorithm generates rough boundary based on an iterative removal of Delaunay triangulation. Furthermore, a mathematic formulation for cavities and holes was given and a statistical method to detect them was designed. Finally, a series of experiments including both simulated and real data sets to validate the effectiveness and practicability of our algorithm was conducted.

Key words: planar point set     shape reconstruction     Delaunay triangulation     multi-constraints     GIS    

平面点集形状重构在GIS相关应用领域如地图综合^[-]、建筑物轮廓线提取^[-]、地理范围确定^[-]以及地理信息检索(GIR)^[]中是一项重要而基础的工作，旨在从一堆无序的点集(仅有坐标信息)中提取出平面点集的分布范围，近似地表达真实的轮廓信息。如何考虑点集空间分布的复杂性是许多平面点集形状重构方法面临的主要问题。通过对国内外学者的研究分析，将这种复杂性归纳为以下4个方面：①点集的空间关系表达，包括几何信息和拓扑信息；②空间密度分布不均；③局部凹凸信息的精确描述；④内部空洞问题。此外，为更加精确地描述点集形状轮廓信息，重构过程中还需要顾及重构结果的唯一性和规则性。

现有平面点集形状重构方法可以分为3类：①基于凸壳的方法；②基于Delaunay三角网的方法；③基于曲线重构的方法。基于凸壳的方法^[-]大都采用凸壳內缩的策略获取边界的凹凸特征，在一定情况下能够处理空间密度分布不均的情况，但由于该算法缺少点集的拓扑信息，难以保证边界的规则性，也无法识别空洞现象。基于Delaunay三角网的重构过程可以分为两步：①构建Delaunay三角网；②基于图形结构特征提取一些几何单元来近似地表达真实边界情况。最早的这类重构方法是Sculpture方法^[]和α-shapes方法^[]。Delaunay Sculpture是一种剥离三角网算法，由外向内不断删除不规则的三角形直到产生一个符合条件的点集形状。基于这种思想衍生了许多算法^[-]，这些算法能够提取任意形状的点集，剥离规则也保证了最终边界的规则性，但这些算法不够稳健，难以适应复杂的凹部和空洞现象，如χ-shape算法^[]通过迭代判断Delaunay三角网边界边长，将满足阈值的边长保留作为最终的边界提取结果，这种全局参数对密度分布不均的样本数据重构效果不是很稳定，该算法也没有考虑空洞问题；为克服全局度量受密度变化影响较大的缺陷，∂RGG算法^[]采用了局部度量方式对三角网进行过滤，避免了算法的参数化，针对三角网中特定的结构特征能够有效识别空洞现象，但该算法对含有规则的凹边界和狭长空洞是失效的。α-shapes算法类似Delaunay三角网，参数α控制了轮廓信息的含量,由于这个原因，需要不断变化参数α值以得到理想的轮廓信息。但是，α-shapes不能准确地表达点集密度分布不均的情况。为了顾及局部密度信息,后期出现了若干α-shapes改进算法，如r-shape^[]、A-shapes^[]、LDA-α-complex^[]等，这些算法采用了局部度量和全局度量策略，能够满足平面点集分布不均、空洞等复杂情况，但是输入参数过多。基于曲线重构的方法采用曲线拟合方法如泊松法^[]、最小二乘法^[]、径向基函数法^[]估算点集的曲率，将满足曲率阈值条件的点判定为边界特征候选点，最后将这些边界点连接即达到了边界提取的目的。曲线重构方法可以针对无规则的点云数据提取出精度较高的边界点，但这种方法需要计算每一个数据点的曲率值，其计算过程非常复杂。

综上分析，基于Delaunay三角网的平面点集重构方法原理较简单，其结构能够很好地弥补不规则点集的几何信息和拓扑信息，较其他方法相对系统和稳健。针对该方法存在的问题，本文基于Delaunay三角网提出了一种有效的平面点集重构方法，首先采用一种局部度量方式对整个Delaunay三角网(凸壳)进行过滤得到初始边界，然后通过构建有效的数学模型实现凸凹多边形内外边界提取。

1 理论基础1.1 数据结构

Delaunay三角网中包含了3种实体要素：点、边和三角形，本文采用了一种简单且实用的数据结构，即在Delaunay三角网中每个三角形记录了其相应的三条边，每条边存储了两个端点信息，包括每个端点的坐标信息。此外，为了更有效地对边界进行操作，每条边还记录了邻近三角形信息，每个点存储了邻域边信息。采用以上数据组织方式来表达Delaunay三角网中3种要素之间的拓扑信息对边界搜索、提取以及拓扑检查有着重要的作用。

1.2 相关定义

设平面点集为S，DTS为平面点集S所构建的Delaunay三角网，DTS中的三角形、边、点分别表示为t、e、v。

(1) 边界边：如果DTS中一条边e属于边界边，则这条边只属于一个三角形，即e不是公共边。基于Delaunay三角网的平面点集形状重构方法最终结果是由一系列边界边构成。

(2) 边界三角形：如果DTS中一个三角形属于边界三角形，则这个三角形至少包含一条边界边。

(3) 边界点：如果DTS中一个点v属于边界点，则这个点必定在边界边上，其余的点都是内点。

(4) 边界规则性：在对DTS过滤过程中要保持边界的规则性，必须排除以下几种情况：

a. 悬挂边、链式边：它们不属于任何一个三角形()。

图 1 边界规则性示意图Fig. 1 Illustration of regular and non-regular boundary

b. 交叉点：每个边界点由两条边界边交汇而成，交叉点是两条以上边界边的交点()。

当一个边界边从DTS中删除后，需要对新生成的边界进行拓扑检验，如果边界出现以上几种情况，则边界边予以保留，其他情况下则可以删除。

(5) 凹部^[]：给定一个覆盖点集P的多边形Q，Conv(P)表示P的凸壳,C表示Conv(P)与Q之间的开放区域，每个开放区域用一条封闭曲线C来表示，这些封闭曲线构成的区域称之为凹部。在DTS中，凹部内的边称之为凹边，每个凹部由一条可见边界边和若干凹边构成。

(6) 空洞：从Delaunay三角网的结构而言，空洞是一个封闭的凹部区域，即每个空洞由若干凹边构成而不存在可见边界边，从某种意义上而言，凹部和空洞具有相似的结构特征。

2 算法过程

本文提出的算法是一种由外向内有序过滤Delaunay三角网的过程。首先按照设定的数据结构采用生长法^[]对平面点集构建Delaunay三角网，将点群围在一个凸壳内，找出初始边界；然后根据一种有序的筛选策略对初始边界向内过滤，得到一个大概边界，过滤的过程中要保证边界的规则性；最后进一步考虑边界的凹凸信息和空洞现象，得到最终的准确边界。算法可分为两部分：粗边界提取和精细边界生成。

2.1 粗边界提取

Delaunay三角网提供了3个参数：周长、角度和边长。周长和边长阈值都属于全局度量方式，容易受密度变化的影响，角度作为一种局部度量方式要优于周长和边长。根据格式塔邻近性原则^[]可知平面点集的边界由点群外围且距离相近的点连接而成，换言之，在一个边界三角形中，如果边界边所对应的角度越大，则两个边界点之间的距离也就越大，两个边界点之间不符合邻近性原则，应当删除它们之间的连线(即边界边)。上述性质能够适应空间实体密度的变化，为此本文采用角度作为重构参数。给定一个平面点集，设边界边为Be,边界三角形中边界边所对应的角度表示为Angle(Be)，角度阈值为R，其粗边界提取的具体流程如所示。

图 2 粗边界提取流程Fig. 2 Rough boundary extraction flow chart

粗边界提取过程采用了排序-过滤的策略，这是因为边界处的几何特征与内部相比具有明显的“不规则性”，排序的效果是为了过滤掉最不规则的元素，保证边界提取结果的唯一性；另一方面粗边界提取算法虽然在一定程度上能够反映形状轮廓的凹凸信息，但是无法识别复杂的凹部和空洞现象，需要进一步施加约束条件，生成精细边界。

2.2 精细边界生成

粗边界提取过程采用了基于角度的过滤策略对以下两种结构是失效的。

(1) 含有规则三角形的凹部。如果凹部存在如所示的规则三角形时，其作为边界三角形边界边所对应的角度很小，为删除此边界边必须降低阈值R，但阈值过小，边界会过度收缩，难以得到理想的重构效果。如果无法删除此边界边，势必会阻碍凹部的进一步探测。

图 3 含有规则三角形的凹部Fig. 3 An example of cavity with regular triangle

(2) 空洞。粗边界提取是由外向内过滤边界边的內缩过程，由边界边的定义可知，空洞内的边()不属于边界边，因此，空洞现象无法被识别。

图 4 空洞现象Fig. 4 An example of a hole

在1.2节中指出凹部和空洞可以视为同一种结构，它们通常位于较周边密度稀疏的区域，采用Voronoi图来表示这种密度变化，从和中可以看出凹部和空洞区域密度较小，其周边区域的密度相对较高。选用密度变化来定义这两种结构的原因在于Delaunay三角网是一个带有密度性质的空间邻近性模型，不需要任何参数设置，参数信息完全包含在三角网中。从和中可以看出，位于这两种结构的边界点其邻域内连线长短不一，梯度变化较大。基于此，采用一个相对参数来量化这种特征，即通过判断点的松散度来探测凹部和空洞，具体表达如下。

点的松散度：给定一个平面点集D，由点集生成的Delaunay三角网表示为DT。针对DT内任一点实体P，F(P)表示点实体P的松散度，表示为

 (1)

 (2)

 (3)

式中，Local_Mean_Length(P)表示点P的邻域内所有边长的均值；Local_SD(P)表示点P邻域内所有边长的标准差；d(P)表示点P直接邻近对象的个数；e_i表示点P邻域内的边长。

F(P)值越小，说明该对象邻域内边长梯度变化小。真实轮廓内的点由于邻域内边长变化较小，其F(P)值较小，相反，空洞和凹部处的边界点F(P)值较大(和)。因此，针对DT内任一点实体P，凹部和空洞的松散度可以表达为

 (4)

式中，Sets表示凹部和空洞的松散度集合；T表示松散度阈值。通过集合Sets可以找到空洞和凹部的边界点，这些边界点邻域内连线长短不一，如何筛选这些边长是关键问题。从空间聚类的角度出发，松散度较大的主要原因在于长边的存在，为删除这些不合理的长边，本文采用了AUTOCLUST算法^[]，其根据密度信息来探测长边，具体表达如下。

长边约束：给定一个平面点集D，由点集生成的Delaunay三角网表示为DT。针对DT内任一点实体P，Longedge表示与P连接的所有边的长边约束条件，表示为

 (5)

 (6)

式中，Local_Mean_Length(P)表示P邻域内所有边长的均值；Local_SD(P_i)表示点P_i邻域内所有边长的标准差；Global_SD表示DT内所有对象Local_SD(P_i)的平均值；N表示平面点集数目。

针对网内任一点实体P，若与其直接相连接的边的长度大于Longedge，则删除该边。对于保留下来的短边可以调用粗边界提取进行约束。基于此，精细边界提取的具体过程如下。

(1) 在粗边界提取的结果基础上，计算整个点集的松散度并获取集合Sets。

(2) 将Sets分为两部分：边界点的松散度和内点松散度，如果点集同时存在凹部和空洞，凹部识别优先于空洞，即边界点遍历优先于内点。

精细边界的提取过程如和所示。

图 5 凹部识别流程Fig. 5 Concave recognition flow chart

图 6 空洞识别流程Fig. 6 Hole recognition flow chart

精细边界提取过程中顾及了两个顺序，一是优先级，凹部和空洞识别的优先级，邻域内边长和边界边的遍历优先级，优先级能够确保过滤掉最不规则的元素；另一个顺序是遍历顺序，松散度遍历和邻域边长遍历都采用了降序排列，这确保了边界提取的唯一性。

2.3 松散度阈值T的确定

T值是识别凹部和空洞边界的重要参数，为了保证阈值T能够适应不同点集的需求，设计了一种自适应的阈值确定方法。由点的松散度内涵可知，边界处点的松散度要大于真实轮廓内点的松散度，本文将阈值选择看作是一种二类分割问题，通过目标函数对点集进行判别以寻找到边界处与真实轮廓之间过渡的点。具体过程为首先计算整个点集的松散度F(P)并按照升序进行排列，记为集合M=F(P₁)，F(P₂)，…，F(P_N)；然后通过目标函数对集合M进行顺序判别，如果出现不满足目标函数的点，判别结束则此点的松散度即为阈值T。自适应阈值确定方法的目标函数表示为

 (7)

式中，x_i表示FP_i值；表示集合M的均值；σ表示集合M的标准差；k是调节参数,k越大，则松散度阈值T越大，网内较多的边将会被保留。

目标函数通过探测过渡点将松散度划分为两部分，为达到二类分割的最佳效果即类内方差最小或类间方差最大，引入PBM^[]指数来求解参数k。PBM指数作为一个相对评价指标能够满足紧密性与分离性的聚类准则，其具体表示为

 (8)

 (9)

式中，N_c表示簇的数量；N_i表示簇C_i中实体数量；v_i表示簇i的质心；E_i表示簇C_i的内部距离(簇内所有空间对象到其质心的距离之和)；E₁表示数据集只分为一类的聚类内部距离；E_Nc表示数据集分为N_c个簇的聚类内部距离；D_Nc表示空间簇间的分离度，其随着N_c增大而增大，最大值为数据集中最远两个簇的质心距离。PBM指数越大，则表示得到的聚类结果越可靠。

结合粗边界和精细边界提取过程，分别对含有规则三角形的凹部和空洞机进行识别。首先调用粗边界提取流程对含有凹部或者空洞的Delaunay三角网(参考和)进行过滤，其结果如和所示；在粗边界提取的基础上调用精细边界提取流程，通过PBM指数确定松散度阈值T(参考—和—)以此找到位于凹部或者空洞处的边界点(和)，筛选边界点的邻域边从而得到最终的边界(和)。

图 7 含有规则三角形的凹部识别过程Fig. 7 Different stages of the proposed algorithm for cavity detection

图 8 空洞的识别过程Fig. 8 Different stages of the proposed algorithm for hole detection

3 试验分析3.1 参数R的选择

粗边界和精细边界的生成过程中，角度阈值R是边界收敛的重要约束条件。平面点集存在两种形式：一是含有空洞和凹部的点集，另一个是不含有以上两种现象的点集。针对不含有空洞和凹部的点集，Delaunay三角网内的三角形近似等边三角形，此类点集的凸壳边界即为理想的边界收敛结果。对于含有空洞和凹部的点集，需要深入讨论角度阈值R对此类点集边界收敛的作用过程，为此，本文设计了几组模拟试验来评价参数R对重构结果的精度影响，以获取合适的参数取值范围为用户在参数R选择方面提供参考。模拟试验采用既定的形状和点集分布类型，L2误差范数^[]作为重构结果的精度评价方法，其具体表达为

 (10)

式中,O代表原始多边形；S代表重构结果。L2误差范数通过原始多边形与重构结果的面积比例大小评价重构效果，L2误差范数等于0，则表示重构结果与原始多边形边界完全吻合。

模拟实验采用了5组形状数据：云南、内蒙古、F字(宋体)、K字(宋体)和A字(宋体)，5组形状数据边界具有明显的不规则性，包含了大量的凹部和空洞；考虑到点集分布类型(如随机分布)有可能使点集内部产生空洞现象，因此点集分布类型选择均匀分布；每组形状数据内部生成的点集数目确保能够覆盖到形状边界，模拟实验采用的点集数目分别是70、100、130和160，代表了不同的点集密度。由角度参数的内涵可知阈值R的取值范围在[0, 180](试验取值间隔为10)，显示了基于参数R每组形状数据边界收敛的L2误差范数变化曲线，5组数据的边界收敛过程大体相似这主要是因为空洞和凹部是等价的结构；另外，每组数据的边界收敛过程并没有随点集密度变化而发生改变。当R取值180时粗边界中Delaunay三角网没有边被删除(如所示)，网内存在大量的“不规则”长边，此时精细边界中的松散度阈值T也较高(如所示)，只能删除部分不规则长边，无法彻底识别凹部或者空洞，边界收敛效果不理想，表现为L2误差范数值较高(如所示)；当R取值0时，网内的边本该全部删除，但是由于受到边界规则性的约束边界边不能存在悬挂边、链式边和交叉点，保留下来的边界呈“锯齿状”，如和所示，边界收敛效果同样不理想，L2误差范数值较高，如所示，因此要想获得理想的边界收敛效果R的取值范围应该是0＜R＜180。在180到0这个区间内随着参数R取值逐渐减少粗边界不断向内过滤，如所示，精细边界中的松散度阈值T也相应变小，如所示，阈值R在一定区间内凹部或者空洞内的边得以彻底识别并删除，边界收敛效果较理想，如中R取[60,110]的结果，表现为L2误差范数值较低，如所示。点集中的凹部或者空洞彻底识别后，Delaunay三角网内松散度差异会很小，二类分割的结果有可能是松散度同属一类，那么粗边界结果与精细边界结果是等价的，如和中R取[40,80]的结果，如果继续降低阈值R向内过滤，边界收敛效果会变差，表现为L2误差范数值不断提高，如所示。

图 9 粗边界提取过程Fig. 9 Rough boundary extraction process

图 10 精细边界提取过程Fig. 10 Refined boundary extraction process

图 11 边界收敛的精度变化Fig. 11 Precision variation of boundary convergence

结合模拟试验的边界收敛过程和中不同点集密度边界收敛的L2误差范数变化过程可以得出，角度阈值R取值范围在[80,120]内就含有空洞和凹部的点集而言可以取得理想的边界收敛结果。

3.2 试验结果对比

为了验证本文方法的有效性与优势性，将提出的方法与其他两种基于Delaunay三角网的典型方法χ-shape和∂RGG算法进行比较分析，试验内容主要包含以下几个方面：

3.2.1 非均匀分布

分别显示了在平面点集分布不均情况下3种方法的边界提取结果。χ-shape算法由于采用边长阈值作为度量参数，当点集分布不均时边界提取效果不是很稳定，会产生大块的锯齿现象；∂RGG算法和本文方法都采用了局部度量参数(三角形外接圆心和角度)能够获取合理的重构结果；另外，平面点集形状重构本质上属于不确定性问题，即合理的重构结果并不唯一。虽然∂RGG是一种自动构建形状轮廓的方法，不需要任何参数，但从不确定性角度而言，本文提出的方法能够产生若干合理的重构结果，提供了更多的重构细节。

图 12 针对点集分布不均χ-shape(绿色)、∂RGG(蓝色)和本文方法(白色)边界提取结果的定性比较Fig. 12 Qualitative comparison of χ-shape (in green)) and ∂RGG(in cyan)) with the proposed algorithm (in black boundary) for non-uniform point set

3.2.2 含有规则三角形的凹部

分别显示了3种方法对含有规则三角形的凹部处理结果。当凹部内点集分布均匀时()，∂RGG算法采用的局部度量方式只考虑到了狭长的边界三角形，无法探测角度较小的三角形，因此凹部探测不彻底；本文方法和χ-shape算法都能够准确地识别这一类凹部。当凹部内点集分布不均时()，与χ-shape算法相比，本文针对凹部给出了合理的数学定义以及识别方法，具有较强的理论保证，从而能够获取准确的重构结果。

图 13 针对含有规则三角形的凹部χ-shape(绿色)、∂RGG(蓝色)和本文方法(白色)边界提取结果的定性比较Fig. 13 Qualitative comparison of χ-shape (in green)) and ∂RGG(in cyan)) with the proposed algorithm (in black boundary) for Cavities with regular triangles

3.2.3 空洞现象

分别显示了3种方法对空洞现象的处理结果。χ-shape算法无法识别空洞；∂RGG算法针对三角网的结构可以对空洞现象构造边界，但它无法识别狭长的空洞现象；本文方法能够识别多样的空洞现象，包括狭长空洞()和球状空洞()。

图 14 针对空洞现象χ-shape(绿色)、∂RGG(蓝色)和本文方法(白色)边界提取结果的定性比较Fig. 14 Qualitative comparison of χ-shape (in green)) and ∂RGG(in cyan)) with the proposed algorithm (in black boundary) for thin hole detection

3.2.4 噪声和非连通区域

3种方法的重构结果都包含了平面点集中所有点包括噪声，若要解决噪声问题可以采用空间聚类的方法如AMOEA算法^[]对平面点集进行预处理，剔除噪声；在一些实例中，有些点集是非连通区域如“i”或者“=”,考虑到边界的规则性，三种方法同样不能生成非连通区域。

3.3 应用实例

本节将引入两个应用实例来验证本文方法的实用性，以发现平面点集重构方法在不同领域内的应用潜力。城市轮廓是城市空间形态分析中一个重要的组成部分，本文基于POI数据利用提出的方法提取城市轮廓线。显示了洛阳市地名地址POI数据(灰色部分)，由于城市轮廓线表达的是一个城市的大概分布范围，因此这里采用粗边界提取方法和凹部识别算法，不需要构造城市内部的边界。给出了洛阳市城市轮廓提取结果(黑色曲线)，将提取的城市轮廓线与洛阳市遥感影像图()叠加比较，结果如所示，红色曲线代表洛阳市遥感影像边界线，黑色曲线代表轮廓提取结果，二者边界线有很高的吻合度。

图 15 洛阳城市轮廓提取结果及对比Fig. 15 The urban contour of Luoyang generated by our algorithm and comparison result

另一个实例是基于LIDAR数据的建筑轮廓线提取，数据来源于文献[]给出的马来西亚吉隆坡城市中心区各类不同几何形状的建筑，基于本文方法提取的各种形状建筑轮廓线结果(如所示)与文献[]的提取结果相似，这表明提出的方法可适用于凸凹多边形内外边界提取。

图 16 各种形状建筑轮廓线提取结果Fig. 16 The extracting results of different building boundary with our algorithm

4 结 论

本文提出了一种多约束的平面点集形状重构方法，相比已有方法，该方法能够更好地适应平面点集空间分布的复杂性。通过模拟试验和应用实例可以发现：①本文方法可以适应空间密度分布不均、凹凸信息精确描述、内部空洞等复杂情况下的平面点集形状重构；②需要的输入参数只有一个且在试验部分中给出了指导信息，方便用户实际使用；③重构过程中对平面点集中存在的凹部和空洞给出了一种有效的数学描述以及识别方法，具有较强的理论保证。进一步改进的方面有：实验分析中缺少对算法自身的定量分析包括点集密度变化和点集分布类型的适应性分析；在应用领域中还局限在二维平面中，后期进一步考虑拓展到三维空间中。

【引文格式】朱杰，孙毅中。 多约束的平面点集形状重构方法[J]. 测绘学报，2017，46(2)：253-264. DOI: 10.11947/j.AGCS.2017.20160122

更多精彩内容：

权威 | 专业 | 学术 | 前沿

微信投稿邮箱 | [e-mail]

微信公众号中搜索「测绘学报」，关注我们，长按上图二维码，关注学术前沿动态。

欢迎加入《测绘学报》作者QQ群： 297834524

进群请备注：姓名+单位+稿件编号