分治法重建数字地形的子网凸包合并算法
郑辑涛1, 2, 何红红2, 张涛2, 朱纪洪1     
1. 清华大学 计算机科学与技术系, 北京 100084;
2. 北京航空气象研究所, 北京 100085
摘要:该文深入研究了以往基于分治策略的数字地形重建方法,在实践基础上分析了参与合并的两个子网凸包的各种可能情况,针对传统合并算法的局限性和弊端,给出了一个子网凸包合并算法。该算法根据公共支撑线的性质,通过判断凸包顶点投影位置的关系确定公共支撑线和支撑点,然后在确保合理的前提下,在两个凸包之间交替生成新三角形完成两子网凸包的合并。实验结果表明: 该算法稳定可靠,能够实现各种复杂情况下两子网凸包的成功合并。
关键词子网    凸包    合并    数字地形    重建    
Subnet convex hull merging algorithm for reconstructing digital terrain models
ZHENG Jitao1, 2, HE Honghong2, ZHANG Tao2, ZHU Jihong1     
1. Department of Computer Science and Technology, Tsinghua University, Beijing 100084, China;
2. Beijing Aviation Meteorological Institute, Beijing 100085, China
Abstract: There are many digital terrain reconstructing methods based on divide and conquer. This study analyzes all possible cases of the merging of two subnet convex hulls. A subnet merging algorithm is then developed to overcome the limitations and drawbacks of traditional algorithms. The public support line and support points are found from the projection position of the vertices of the convex hulls. Then, the convex hull vertices between the support points are linked to generate new triangles so that no two triangles are intersect. Tests show that the algorithm is stable and reliable and that any two subnet convex hulls can be successfully merged.
Key words: subnet    convex hull    merge    digital terrain model    Reconstruction    

离散点云重建数字地形三角网的算法有: 插入法[1, 2]、 生长法[3, 4]和分治法[5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]。三维激光扫描技术的进步使得获取海量地形表面点云更加容易。分治法由于其建网效率高而受到研究者的普遍关注。分治法通常需要将完整的地形点云分割成若干子集,分别建立子网再合并。子网的建立可以采用逐点插入法[5, 6]、 三角网生长法[7]和扫描线法[8],通常子网建立后其边界为凸多边形,因此本文也称其为子网凸包。子网凸包合并是分治法快速重建数字地形的关键。Shamos和Hoey[9]最早提出了基于分治法生成维诺图的算法,随后Lewis和Robinson[10]加以改进用于重建Delaunay三角网。其他学者也分别对分治法数字地形建模进行了研究和改进。例如,Lee和Schachter[11]在假设凸包边不存在3点以上共线的情况下提出的合并算法在凸包边多点共线情况下并不适用; 文[12]中的算法从凸包中一个给定初始点开始搜索符合条件的顶点,通过斜率和最大张角来判断识别另一个凸包中的支撑点,由于斜率的变化存在正负突变,有的情况下会由于找不到符合条件的支撑点而使算法失效; 文[13, 14]通过最高最低点来判断支撑点和公共支撑线。

本文在海量数据数字地形建模的实践中发现已有的合并算法由于存在缺陷常常不能正确合并。因此,本文在前人的研究基础上,深入研究了海量数据分治法构建数字地形三角网中的凸包合并问题,对分割后参与合并的凸包所有可能出现的情况进行了仔细的分析,提出了一个新的凸包合并算法,通过2个子网交替外推搜索二者公共支撑线(公切线)和支撑点(公切点),然后再对2个凸包内侧公共支撑点之间的部分重剖分,使得新生成的三角面片没有重叠和交叉,完成2个子网的合并。

1 基本概念

分治法重构数字地形,首先要把完整地形点云分割成若干数量较小的点云子集,然后分别对分割后的点云子集建立子三角网,子三角网的外边界构成的多边形是凸多边形,这个凸多边形就被称为子三角网的凸包。例如,图1中2个粗实线构成的凸多边形分别是2个子网的凸包。通常可以用构成凸包的顶点集合的有序序列(不失一般性本文以逆时针序为例)来表示一个凸包。因此,图1中的2个凸包可以分别表示为2个集合A={a1,a2,…,a7}和B={b1,b2,…,b8}。凸包的每一条边的方向也是指向凸包的逆时针方向。如图1中的合并前凸包边矢量$\overrightarrow {{a_7}{a_1}} $、 $\overrightarrow {{b_6}{b_7}} $等,合并后凸包边矢量$\overrightarrow {{a_6}{b_6}} $。对于凸包中任一顶点ai,以ai为尾的边矢量称为ai的前边,以ai为头的边矢量称为ai的后边。例如,凸包A中$\overrightarrow {{a_6}{b_7}} $为顶点a6的前边,$\overrightarrow {{a_5}{a_6}} $为a6后的边。

图 1 子网凸包与凸包合并示意图

如果一个凸包不存在任意3个以上的顶点共线,则被称为标准凸包,否则就是一般凸包。

凸包的一条支撑线是指通过凸包的某个顶点的一条直线,且满足凸包完全在该直线的一侧。不相交或不重叠的两凸包存在共同的支撑线,称之为公共支撑线(或称公切线)。对应的顶点称为支撑点(或称公切点)。图1中$\overrightarrow {{a_6}{b_6}} $和$\overrightarrow {{b_2}{a_2}} $即为AB两凸包的公共支撑线,其中a2a6b2b6称为两凸包的公共支撑点。

不重叠的两凸包合并后生成一新凸包,其中2条公共支撑线成为合并后新的凸包边。图1中凸包AB合并后的凸包可以用集合H={b1,b2,a2,a3,a4,a5,a6,b6,b7,b8}表示。处于任意位置的2个不重叠的凸包总可以通过旋转变换成左右位置关系。不失一般性,为了描述问题方便,本文仅考虑左右凸包的合并问题,结果很容易推广到任意位置关系。

凸包合并需要正确找到参与合并的2个凸包的公共支撑线和对应的支撑点。

2 凸包合并

合并合并的步骤包括: 1) 找到参与合并的2个凸包的公共支撑线和支撑点;2) 正确连接2个凸包内侧位于公共支撑点之间的顶点,生成合理三角形;3) 更新合并后的凸包边。

2.1 凸包间相对位置关系

经过对大量数据的分析和分割试验,参与合并的2个子网凸包相对公共支撑线可总结为6种情况,如图2所示。下面以下公共支撑线为例说明。

图 2 2个凸包位置关系

1) 2个凸包中y坐标最小的顶点即是下支撑点,不需调整,如图2a所示。

2) 2个凸包中至少有一个凸包存在多点共线且与公共支撑线重合。此时若多点共线出现在左侧凸包中应取共线顶点中最右侧顶点,若多点共线出现在右侧凸包中应取共线顶点中最左侧顶点,此种情况如图2b所示。传统的凸包合并算法此时往往会失效。

3) 2个凸包中,其中一个凸包的最低点是支撑点,另一个凸包的最低点不是支撑点,通常需要通过判断调整,此情况如图2c和2d所示。

4) 2个凸包尺寸大小悬殊,有时导致在小凸包中根本不存在满足另一凸包完全在其一侧的边。例如图2e中左侧较小凸包中,假设每个凸包边沿逆时针方向,没有任何一条边满足另一凸包完全在其左侧。

5) 参与合并的2个凸包的最低点都是非支撑点,通常需要在两凸包中交替调整多次,才能最终确定公共支撑点和公共支撑线,如图2f所示。

2.2 公共支撑线性质

凸包中x坐标值最小和最大的点分别称为最左和最右点,当具有相同的最小或最大x值的顶点存在多个时,取其中y值最小的点。图3中用alblarbr分别表示2个凸包的最左和最右点,从最左点到最右点的逆时针序点集称为凸包的下半链,如图3中两个半链{al,al+1,…,ar-1,ar}和{bl,bl+1,…,br-1,br}。同理可以定义上半链。

图 3 支撑点与半链之间的关系

定理1 逆时针序表示的凸包,合并后的凸包总是位于公共支撑线的左侧。

由公共支撑线的概念显然得到本性质。

定理2 顶点逆时针序排列的凸包,上支撑点总是出现在上半链中,下支撑点总是出现在下半链中。

证明 以图3左侧凸包下半链为例。假设左侧凸包的下支撑点出现在下半链之外的任一顶点,不妨设为ai,连接alar做射线ls,如果过ai的公共支撑线如luls上方,显然部分顶点将位于其右侧,如果过ai的公共支撑线如ldls下方,由于其永远不能到达极限位置lc,因此至少顶点al会位于其右侧,这与下支撑线定义矛盾,假设不成立。同理可以证明上支撑线情况。

2.3 支撑点检测

不失一般性,仅以下支撑点的检测为例,很容易推广至上支撑点。为方便描述,用Prj$\overrightarrow u $ai表示点ai在矢量$\overrightarrow u $上的投影位置。

2点在同一矢量上投影,沿矢量方向靠前的投影位置大于靠后的投影位置; 如2个投影位置重合,则定义为相等。便如图4中,可以表示为Prj$\overrightarrow u $a>Prj$\overrightarrow u $b

图 4 矢量上的投影位置关系

根据定理2可以缩小检测范围,只需在凸包的下半链中查找。首先选择参与合并的两凸包下半链中的最右侧点为初始点,连接此2点作为下公共支撑线,然后交替检测左右凸包是否满足都在此支撑线左侧。具体检测方法和过程如图5所示。

图 5 支撑线检测示意图

选择初始点支撑点a6b8,连接a6b8,把矢量$\overrightarrow {{a_6}{b_8}} $右转90°,生成一个与其垂直的矢量$/overrightarrow {{u_1}}$,参与合并的2个凸包下半链各顶点向$/overrightarrow {{u_1}}$投影,分别找到各自投影位置最大的顶点,即满足

$\begin{array}{*{20}{c}} {\Pr j_u^{{a_5}} = \max \Pr j_u^{{a_i}},i = 3,4,5,6;}\\ {\Pr j_u^{{a_6}} = \max \Pr j_u^{{b_j}},i = 3,4,5,6,7,8.} \end{array}$

判断下半链中在当前ui→方向上投影位置最大的顶点是否与当前选定的支撑点相同。若有不同则更新为当前投影位置最大的顶点。重复上述检测过程,直至某次检测后两个支撑点即为最大投影位置点,则当前支撑点即为最终支撑点。同理可进行上支撑点检测。

显然对于不相交的2个凸包,肯定存在上下2条公共支撑线和对应的4个支撑点,因此总能经过有限步检测后确定公共支撑线和支撑点。

2.4 子网合并

检测到两个子网凸包的4个支撑点,需要把二者合并成一个子网,同时保证边缘仍然为1个凸包。

顺序排列两凸包支撑点之间内侧的顶点序列,如图6所示,分别为{ad,ad+1,…,au}和{bd,bd-1,…,bu},然后从任意公共支撑线开始,逐渐生成新三角面,直至另一条公共支撑线。每生成一个新的三角面需要判断其合理性,即应保证三角形之间不能重叠或相交,通常可以根据新生成的三角形内部是否包含其他顶点来判断其合理性。

图 6 子网合并示意图 Fig.6

图6中,设当前边为$\overline {{a_d}{b_d}} $,可以取左侧凸包的相邻点ad+1构成△adbdad+1,此时需要判断右侧凸包的相邻点bd-1是否在该三角形内; 也可以取右侧凸包的相邻点bd-1构成△adbdbd-1,同样需要判断左侧凸包的相邻点ad+1是否在其内部。判断后取合理者,如果都合理任取其一即可。实践中在都合理的情况下可以左右交替取相邻点生成新三角形,可以尽量避免狭长三角形的出现。

合理性判断只需对相邻的一个点进行判断,这是由凸包的性质保证的。例如图6中取右边凸包的相邻点bd-1构成△adbdbd-1,只需对左侧凸包中的相邻点ad+1 判断。这里分别考察左右侧凸包AB,因为$\overline {{b_{d - 1}}{b_d}} $是B的凸包边并且adB的外侧,△adbdbd-1中不可能包含B中顶点; A中若有顶点在△adbdbd-1内,则ad+1一定在其内部,因为若ad+1不在△adbdbd-1内而其他点在,则必然会出现某顶点在凸包边$\overline {{a_d}{a_{d + 1}}} $另一侧的情况,使得A不满足凸包性而产生矛盾。合并后更新凸包边界。

2.5 子网凸包合并算法

输入为2个子网凸包A={a1,…,an}和B={b1,…,bm};

输出为合并后子网及其凸包C={c1,c2,…,cq}。

算法步骤为:

步骤1 分别计算两凸包上下半链AdBdAuBu

步骤2 选择初始下支撑点adbd;

步骤3 当前矢量$\overrightarrow {{a_d}{b_d}} $右旋90°生成新矢量$\overrightarrow u $;

步骤4 2个凸包下半链顶点在$\overrightarrow u $方向投影排序,找到各自投影位置最大的顶点,如果与当前两下支撑点不同,则替换当前下支撑点,转步骤3;

步骤5 选择初始上支撑点aubu;

步骤6 当前矢量$\overrightarrow {{a_u}{b_u}} $右旋90°生成新矢量$\overrightarrow \upsilon $;

步骤7 两凸包上半链顶点在$\overrightarrow \upsilon $投影排序,找到各自投影位置最大的顶点,如果与当前两上支撑点不同,则替换当前上支撑点,转步骤5;

步骤8 从下支撑线$\overrightarrow {{a_d}{b_d}} $到上支撑线$\overrightarrow {{a_u}{b_u}} $依次合并生成新三角面片;

步骤9 合并后更新凸包边。

3 复杂度分析与实验

假设两个子网凸包的顶点个数分别为mn,则算法中涉及的主要操作有上下两个半链的查找,比较次数最多(m+n)次; 检索上下公共支撑线过程中投影位置的计算,最多不会超过(m+n)次; 合并时内部三角面的生成,最坏情况(m+n)个; 更新凸包边,删除增加新的顶点序列,最坏不超过(m+n)。因此整个算法的时间复杂度为O(m+n)。

采用真实地形数据实验,数据分割后采用四叉树组织,其中子网的建立采用文[8]中的扫描线算法。每个节点的合并首先是上下2个子集分别进行左右合并,然后再进行一次上下合并,最终合并成一个完整的子网,如图7所示。

图 7 四叉树组织的地形节点合并过程

图8为某地形的建模结果局部放大图,点云数量1 469 801个,生成三角面2 836 234个。

图 8 某地形数字地形重建结果

采用图8中地形进行实验。计算机配置为: CPU为Intel酷睿i3,主频2.6 GHz,内存2 GB,采用不同阈值分割点云数据,测试结果如表1所示。

表 1 不同阈值分割后的合并次数与耗时
阈值 合并次数 耗时/s
10 108 332 12.208
20 55 296 11.312
100 11 474 12.263
500 2 262 18.984
1 000 1 155 25.771

实验结果表明,本文的分割算法稳定性好效率高,并且分割阈值取lbN时效最高。阈值过大导致子网规模过大,建立子网的时间较长; 阈值过小,合并次数太多同样使得时间开销增大。取阈值为lbN时对不同的地形数据进行试验,结果如表2所示。可以看出,采用分治法建立地形三角网效率高,阈值取lbN时建网的时间复杂度近似O(n),生成三角面数量为点云数量的2倍。

表 2 不同地形分治法测试结果
点数 阈值 合并次数 三角面数 耗时/s
67 102 16 3 209 129 380 0.515
347 683 18 15 566 694 686 3.652
1 469 801 20 55 296 2 836 234 11.312

可见本文的合并算法稳定性好,大数据量情况下能有效地完成子网合并最终建立数字地形三角网。

4 结 论

本文在试验基础上分析列举了数字地形重建中分割后参与合并的子网凸包的各种情况,针对传统合并算法的局限性和弊端,证明了公共支撑线的性质,根据此性质提出了基于投影位置关系的支撑点检测规则,给出了子网合并中新三角面生成的原则与合理性判据,在此基础上给出了一个子网凸包合并新算法。本文从理论上证明了该算法的正确性,通过实验证实了其稳定可靠,能够成功用于海量数据三维数字地形快速重建。

参考文献
[1] 袁正午, 侯林, 彭军还. 点集收集分配的 Delaunay 三角网快速生成算法及实现 [J]. 测绘科学, 2011, 36(5): 223-225.YUAN Zhengwu, HOU Lin, PENG Junhuan. A fast algorithm of Delaunay triangulation generation based on collecting and distributing points [J]. Science of Surveying and Mapping, 2011, 36(5): 223-225. (in Chinese)
[2] 王龙浩, 王解先. 基于逐点插入法的 Delaunay 三角网快速生成算法 [J]. 工程勘察, 2013, 10: 75-79.WANG Longhao, WANG Jiexian. Fast Delaunay triangulation generation algorithm based on incremental insertion method [J]. Geotechnical Investigation and Surveying, 2013, 10: 75-79. (in Chinese)
[3] 詹曦, 张建生. 点云边界提取及三角网格生成的集成算法研究 [J]. 计算机仿真, 2013, 30(11): 272-275.ZHAN Xi, ZHANG Jiansheng. Research on algorithm of integrating boundary extraction and triangle mesh generation of point cloud [J]. Computer Simulation, 2013, 30(11): 272-275. (in Chinese)
[4] 刘永和, 王燕平, 齐永安. 一种快速生成平面Delaunay三角网的横向扩张法 [J]. 地球信息科学, 2008, 10(1): 20-25.LIU Yonghe, WANG Yanping, QI Yongan. Horizontal expanding method: A quick algorithm for generating Delaunay triangulation from points in the plane [J]. Geo-Information Science, 2008, 10(1): 20-25. (in Chinese)
[5] 武晓波, 王世新, 肖春生. 一种生成Delaunay三角网的合成算法 [J]. 遥感学报, 2000, 4(1): 32-35.WU Xiaobo, WANG Shixin, XIAO Chunsheng. A hybridized method for building Delaunay triangulation [J]. Journal of Remote Sensing, 2000, 4(1): 32-35.(in Chinese)
[6] 吴宇晓, 张登荣. 生成Delaunay三角网的快速合成算法 [J]. 浙江大学学报: 理学版, 2004, 31(3): 343-348.WU Yuxiao, ZHANG Dengrong. Improved algorithm for building Delaunay triangulation [J]. Journal of Zhejiang University: Science Edition, 2004, 31(3): 343-348. (in Chinese)
[7] 徐青, 常歌, 杨力. 基于自适应分块的TIN三角网建立算法 [J]. 中国图像图形学报, 2000, 5A(6): 461-465. (in Chinese)XU Qing, CHANG Ge, YANG Li. The algorithm of TIN generation based on self-adapt clump organization [J], Journal of Image and Graphics, 2000, 5A(6): 461-465. (in Chinese)
[8] 芮一康, 王结臣. Delaunay三角形构网的分治扫描线算法 [J]. 测绘学报, 2007, 36(3): 358-362.RUI Yikang, WANG Jiechen. A new study of compound algorithm based on sweepline and divide-and-conquer algorithms for constructing Delaunay triangulation [J]. Acta Geodaetica et Cartographica Sinica, 2007, 36(3): 358-362. (in Chinese)
[9] Shamos M I, Hoey D. Closest-point problems [C]// Proceedings of the 16th IEEE Symposium on Foundations of Computer Science. Washinghon, USA: IEEE, 1975: 151-162.
[10] Lewis B A, Robinson J S. Triangulation of planar regions with application [J]. The Computer Journal, 1978, 21(4): 324-332.
[11] Lee D T, Schachter B J. Two algorithms for constructing a Delaunay triangulation [J]. International Journal of Computer and Information Sciences, 1980, 9(3): 219-242.
[12] 谢传节, 万洪涛. 基于四叉树结构的数字地表模型快速生成算法设计 [J]. 中国图像图形学报, 2002, 7A(4): 394-399.XIE Chuanjie, WAN Hongtao. Algorithm for rapidly generating digital terrain model based on quad-tree [J]. Journal of Image and Graphics, 2002, 7A(4): 394-399. (in Chinese)
[13] 胡金星, 潘懋, 马照亭, 等. 高效构建Delaunay三角网数字地形模型算法研究 [J].北京大学学报, 2003, 39(5): 736-741. HU Jinxing, PAN Mao, MA Zhaoting, et al. Study on faster algorithm for constructing Delaunay triangulations DTM [J]. Scicentiarum Naturalum Universitis Pekinesis, 2003, 39(5): 736-741. (in Chinese)
[14] 罗胜, 王 鑫, 孙玉平. 机载LiDAR点云的Delaunay三角网快速生成算法 [J]. 海洋测绘, 2014, 34(2): 18-21.LUO Sheng, WANG Xin, SUN Yuping. An algorithm for quick generation of Delaunay triangular net for airborne LiDAR point cloud [J]. Hydrographic Surveying and Charting, 2014, 34(2): 18-21. (in Chinese)
[15] 刘合辉, 罗勇军. 基于等高线的三角网快速构建与处理 [J]. 测绘工程, 2009, 18(2): 55-58.LIU Hehui, LUO Yongjun. Speedy constructing and processing TIN based on contours [J]. Engineering of Surveying and Mapping, 2009, 18(2): 55-58. (in Chinese)
[16] 刘永和, 王燕平, 齐永安. 一种简单快速的Delaunay三角网逐块生成算法 [J]. 测绘科学, 2008, 33(6): 133-135.LIU Yonghe, WANG Yanping, QI Yongan. A simple and quick block-by-block generating method for generating delaunay triangulation from points in the plane [J]. Science of Surveying and Mapping, 2008, 33(6): 133-135. (in Chinese)