ACTA Scientiarum Naturalium Universitatis Pekinensis
A Plane Sweep Based Arc Splitting and Polygon Auto-construction Algorithm
LIU Yuefeng, SUN Ying, ZHANG Kai, et al
Institude of Remote Sensing and Geographic Information System, Peking University, Beijing 100871; † E-mail: yuefengliu@pku.edu.cn
Abstract Aiming at the deficiency of traditional polygon auto-construction algorithm in automation and time efficiency, which leads to the insufficiency of commercial GIS softwares’ data processing and spatial-temporal analysis ability, an arc splitting and polygon auto-construction algorithm based on plane sweep idea is proposed. Our algorithm contains three features as follows. First, it is a complete process from intersection testing until polygon construction. Next, it takes full advantage of useful information during intersection testing to realize arc splitting and polygon auto-construction at the cost of little algorithm complexity and computing resources. Finally, it avoids the calculation of nested relation and handles degenerate cases of bridge and dangling edge. The result of experiments proves the proposed algorithm improves efficiency significantly in comparison with traditional algorithms. Key words plane sweep algorithm; arc splitting; polygon auto-construction
在地理信息系统(geographical information system, GIS)中, 弧段分割及多边形自动构建在矢量数据处理和空间分析等方面被广泛应用。这个问题的描述如下: 给定平面内一个弧段集合(允许存在自相交和互相交), 输出平面内由这些弧段划分的所有闭合连通子集(即多边形集合)。
在 GIS拓扑关系研究的早期, 研究者们关注从地图数字化到拓扑关系建立的过程[1‒2], 算法的效
率不是研究的重点。随着数据的来源越来越充分,有许多研究探讨如何从既有点线拓扑关系数据或多边形数据中建立多边形拓扑关系[3‒8], 其中包含使用网络大数据WEB-GIS作为数据来源的思想。另外一些研究试图从图论的角度建立多边形拓扑关系, 但图的构建也需要计算方位角度并排序, 依然是一项非常耗费时间的任务[9‒10]。部分研究借鉴栅格数据的处理手段来生成拓扑关系[11‒12], 但是栅格
数据的应用范围有限。也有研究采用分治策略来生成多边形的拓扑信息[13‒14], 但此类算法的实现过程比较复杂, 拼接的时候会耗费很多计算资源, 并且不一定能获得效率的提升。目前, 大多数关于多边形生成的算法都是基于方位角比较的左转多边形追踪算法[15‒23],虽然研究者们构建了很多策略进行优化搜索方向角的选择, 但这一部分的计算仍然耗费大量的时间, 并且这些研究没有解决如何进行弧段相交断链的问题, 也没有关注一些特殊情况(比如悬边、岛等)。
归纳已有的研究成果可知, 多边形拓扑关系的自动构建一般分为4个阶段: 1) 求弧段交点; 2) 对弧段进行分割; 3) 建立结点与弧段的拓扑关系; 4)搜索多边形并建立与弧段的拓扑关系。已有的研究集中在第3和第4阶段, 不足之处在于割断了各阶段之间的有机联系, 不能充分地利用前面阶段的信息, 使得后面阶段算法的复杂度和计算量很大: 在第3阶段, 需要进行结点匹配(将由于计算误差导致的不完全重合的弧段端点合并为一个结点)以及结点和弧段连接关系的全局搜索; 在第4 阶段, 不仅要对构成多边形的弧段进行全局搜索, 还会形成无效多边形(如所有子区域的边界多边形), 同时要分析和计算任意两个多边形之间的包含关系, 以便实现带洞多边形的正确构建。另外, 悬边和桥也是一个不容易处理的问题。
针对上述不足, 本文提出一种基于扫描思想的弧段分割与多边形自动生成算法。该算法有以下特点: 1) 从求交开始至生成多边形结束, 实现多边形自动生成的完整算法; 2) 以 Benttley & Ottmann线段
[24]求交算法 为基本框架进行扩展, 在求交的同时,利用相关信息, 同步实现弧段分割和多边形自动生成, 并建立拓扑关系, 结构紧凑; 3) 基于扫描思想,避免了传统方法中的多边形嵌套关系的复杂计算,并能有效地处理桥和悬边问题。
1 数据结构的扩充
本算法依然使用Bentley & Ottmann算法中表征扫描线状态的数据结构事件点队列Q和状态结构T (第 2节将详细介绍), 并使用平衡二叉树进行组织。Bentley & Ottmann 算法中的线段不包含任何属性, 为了方便弧段分割与多边形生成, 本文对普通线段添加表征方向和拓扑关系的字段, 将其扩充成一个新的数据结构——边(表1)。初始化时, 边的数量等于所有弧段的边数之和。在扫描过程中, 有的边被一分为二(在交点处被分割)。
为了将边有效地组织起来, 表达边与弧段的逻辑关系, 并在弧段分割过程中实现快速更新维护,本文设计一个数据结构: 活化边链表E。对每个弧段的每个边建立一个边对象, 并加入E中, 通过设置它们的 First, Last, Previous 和 Next属性, 建立它们之间的逻辑关系。在扫描过程中不断地更新E,但对E只有末尾插入操作, 而没有删除、中间插入和改变顺序等操作。边与弧段之间逻辑关系的更新仅仅通过对部分边的属性的更新实现, 这样可以避免更新维护过程中边的频繁移动。
图1是一个活化边链表E数据结构应用的简单例子: 两个弧段在交点N处被分割成4个弧段。分
割前的两个弧段分别为Arc(e1-e2-e3)和 Arc(e4-e5-e6),分割后的 4个弧段为 Arc(e1-e2)、 Arc(e7-e3)、Arc (e4-e8)和 Arc(e5-e6)。而边的数量由分割前的6个变成分割后的8个。活化边链表E的数据结构变化见图2。
在实际问题中, 经常出现带洞多边形。为了描述这种现象, 我们将环作为多边形的基本组成部分,一个无洞的多边形由一个环组成, 一个带洞多边形则由一个外环和一个或多个内环组成。
在扫描过程中, 为了将所有可能形成多边形的边记录下来, 并方便提取出来作为输出结果的最终多边形, 本文设计一个数据结构: 活化多边形。一个活化多边形储存的是当前状态下可能形成此多边形的所有边的集合。随着扫描的进行, 活化多边形将不断地被更新, 包括增加边以及修改边的属性等,直至形成一个完整有效的多边形或者被删除。在扫描开始前, 建立一个空的活化多边形集合S, 在扫描过程中, 会有新的活化多边形加入S, 或者已有的活化多边形被删除。由于存在删除操作, 活化多边形集合S采用一棵平衡二叉树进行组织。
2 基于扫描思想的弧段分割与多边形生成算法2.1 整体描述
Bentley & Ottmann算法针对的问题是: 输入一个线段集合, 输出所有交点。此算法有两个重要概念: 事件点和扫描线状态。事件点包括线段端点以及扫描过程中逐渐发现的交点; 扫描线状态指与当前扫描线相交的所有线段构成的集合, 相对于扫描线, 这些线段是自左向右排序的。Bentley & Ottmann算法的总体思路是: 假想有一条水平线自上而下地扫过整个平面。在扫描过程中, 每遇到一个事件点,则对状态结构T进行维护, 并检测新的交点。本文算法是在 Bentley & Ottmann算法的框架下进行扩展, 重点是在处理事件点时增加弧段分割和活化多边形生成, 并在算法末尾增加最终多边形提取。本文算法流程如图3所示, 其中下画线部分为本文算法增加的内容。算法描述如下。
在算法开始时, 初始化一个事件队列Q, 并将所有边的端点插入Q中, 以事件点为上端点的边也要和事件点一起存储。此外, 初始化一个空的状态结构T。
初始化一个空的活化边链表E, 将所有弧段的所有边插入E中。在扫描过程中进行弧段分割时,不断地更新E。
初始化一个活化多边形集合S, 在扫描之前S是空的, 随着扫描的进行, 不断有新的活化多边形加入, 也会有活化多边形被删除。
在完成上述初始化任务后, 逐一处理每个事件点。在处理事件点时, 除维护状态结构T和发现新
的交点外, 还要进行弧段分割和活化多边形生成的操作。
在算法的最后阶段, 从活化多边形集合S中提取并输出所有多边形。本文算法的核心部分为弧段分割、活化多边形生成以及最终多边形提取, 下面分别加以描述。
2.2 弧段分割
设当前处理的事件点为p, 令 U(p)为所有以 p为上端点(如是水平边则为左端点)的边构成的集合,这些边与事件点一起存储, 且接下来将被插入状态结构T中。令L(p)为 T中以p为下端点的边的集合, C(p)为 T中内部包含p点的边的集合。下面根据事件点的不同情况进行描述。
1) |U(p)∪ L(p)∪ C(p)| = 1 (即 3个集合中的边的数量之和为1, 下同), 此时p必定为某弧段的端点(end point)(图 4(a)), 不需要对p点所在的弧段进行分割。
2) |U(p)∪ L(p)∪ C(p)| = 2, 且|C(p)|=0。此时存在以下 3 种情况, 均无需进行弧段分割操作。
① 这两条边属于同一弧段, 并且彼此为相邻边(图 4(b)), 无需进行弧段分割。
② 这两条边属于同一条弧段, 并且一个为起始边, 一个为终止边(图 4(c)), 意味着当前弧段为一条闭合弧段, 也无需进行弧段分割。
③ 这两条边属于不同的弧段, 并且, 或均为起始边, 或均为终止边, 或一个为起始边, 另一个为终止边(图 4(d)), 此时与p点相关的弧段虽然有两条, 但是p点作为两条弧段的端点, 仍无需进行弧段分割。
3) |U(p) ∪ L(p) ∪ C(p)| = 2, 且 |C(p)|>0 。如图4(e)所示, 此时必有两条弧段经过p点, 对这两条弧段进行分割操作。
4) |U(p)∪l(p)∪c(p)| ≥ 3。如图4(f)所示, 此时p点必为弧段交点, 对经过p点的所有弧段进行分割操作。
进行弧段分割时, 只需要对活化边链表进行维护, 算法相对简单, 不详述。
2.3 活化多边形生成
弧段分割完成后, 就可以进行活化多边形生成。此时 T中已经没有包含 p的边存在了, 即|C(p)|=0。如果弧段分割操作前C(p)中有边, 则此时这些边已经被分割并分别插入U(p)和 L(p)中了。下面根据当前事件点分割操作后U(p)和 L(p)的不同情况, 阐述活化多边形的生成算法(为了叙述的简洁, 依然用“多边形”指代“活化多边形”这一概念)。
1) |U(p)| > 1, 且无论L(p)中是否有边。如图5(a)所示, 此时若L(p)中有边, 则这些边在之前均已被处理, 因此无论L(p)中是否有边, 只需考虑U(p)中的边。此时U(p)中的N条边可能构成N−1个多边形(每两个相邻边之间形成1个多边形), 因此新建N−1个多边形, 然后设置N条边的左右多边形指针,并将其加入相应的左右多边形中。其中, 最左侧边(图 5(a)中的e1)的左侧多边形设置为该边左侧紧邻边(图 5(a)中的 el)的右侧多边形, 如果不存在紧邻边, 则设置为空。最右侧边(图 5(a)中的 e3)的右侧多边形的设置方法同理, 下文中不再说明。
2) |U(p)| = 1, 且无论 L(p)中是否有边。如图5(b)所示, 设置 e1左右的多边形, 并将 e1加入其左右多边形中。
3) |U(p)| = 0, 且|L(p)| = 1。如图 5(c)所示, 分以
下几种情况。① 若 e1左右侧均为空, 则不做任何操作。② 若 e1左右仅有一侧为空, 则从活化多边形集合S中删除另一侧多边形, 并将被删除多边形中的所有边指向该多边形一侧的指针调整为空。
③ 若 e1左右两侧均存在多边形且为同一多边形, 则不做任何操作。
④ 若 e1左右两侧均存在多边形且不为同一多边形, 则将左右两侧多边形合并, 并调整相应边的多边形指针(实现方法不详述)。
4) |U(p)| = 0, 且|L(p)| > 1。如图 5(d)所示, 分以下几种情况。
① 若最左侧边(图 5(d) 中的 e1)的左侧和最右侧边(图 5(d)中的 e3)的右侧均为空, 则不做任何操作。
② 若最左侧边的左侧为空, 最右侧边的右侧不为空, 则从活化多边形集合S中删除最右侧边的右侧多边形, 并将被删除多边形中的所有边指向该多边形一侧的指针调整为空; 反之亦然。③ 若最左侧边的左侧和最右侧边右侧均存在多边形且为同一多边形, 则不做任何操作。④ 若最左侧边的左侧和最右侧边右侧均存在多边形且不为同一多边形, 则将最左侧边的左侧多边形和最右侧边的右侧多边形合并, 并调整相应边的多边形指针。
图 6是一个活化多边形生成操作的简单示例。
2.4 从活化多边形中提取最终多边形
扫描结束后, 活化多边形的生成便全部完成。下一步, 需要从活化多边形集合中的每个活化多边形中提取出作为最终输出的多边形。对每个多边形的提取分为两步: 第一步, 通过删除活化多边形中左侧和右侧多边形相同的边, 剔除可能存在的悬链和桥; 第二步, 提取所有环, 这一步需要注意以下问题。
1) 确定追踪法则。由于可能存在内环, 因此不能简单地采用顺时针或逆时针追踪法则, 而应采用右手螺旋法则, 即多边形的内部区域总是处于追踪
方向的左侧。
2) 可能存在多边形自身相切的情况, 即在弧段端点处可能连接了 3 个或 3 个以上的弧段。遇到这种情况, 采用最左转原则。
3) 从活化多边形的边集合中第一条边所在的弧段开始搜索, 根据扫描的原理, 第一个提取的环必定为外环, 之后提取的环必定为内环。
限于篇幅, 提取算法不详述。
3 算法复杂度分析与数据实验3.1 复杂度分析
设输入弧段数为a , 边数为e , 输出的交点数为 I , 弧段数为ao , 边数为eo , 多边形数为no 。
时间代价方面, 线段求交与Bentley & Ottmann算法完全相同, 时间代价为 O(e log e + I log e )。维护活化边链表时仅需顺序插入, 分割弧段时会新建弧段, 二者耗时总和为O(eo)+o(ao-a)。活化多边形数目不超过边数eo, 其生成、维护、合并和删除均为线性操作, 时间代价为O(eo)。因此, 算法总的时间代价为 O(e log e + I log e )+ O(ao) + O(no) + O(eo)。
空间代价方面, 状态结构和事件队列的空间规模正如 Bentley & Ottman 算法所阐述的, 分别为O(eo)和 O(E+I), 活化的边链表的空间规模为O(eo),活化多边形集合与多边形集合的空间规模皆为O(no)。因此, 算法总的空间代价为O( e + I )+ O ( e o) + O ( n o)。
3.2 数据实验
为了验证本算法的有效性, 我们利用模拟数据和实际数据, 通过与左转算法和ARCGIS中的多边形自动生成功能(feature to polygon)进行多次对比测试, 并取平均值。实验采用的计算机配置为2.8 GHZ Intel Core i5 CPU和4 GB RAM, 代码用C#语言实现。
模拟数据如图7所示, 实验结果如图8和表2所示。实际数据包含全国县级行政区域边界数据、全国铁路网数据和全国主要路网数据(图 9, 包括高速公路、一级公路、二级公路和三级公路), 实验结果如表3所示。综合模拟数据和实际数据的实验结果可知, 与传统算法和ARCGIS相比, 本文算法在数据处理效率上明显提升。
4 结论与展望
针对传统算法的不足, 本文提出一种基于扫描思想的弧段分割和多边形自动构建算法。根据理论分析和实验验证, 得到如下结论: 1) 除求交的复杂度与传统算法相同外, 本文算法中的弧段分割以及多边形自动生成均为线性复杂度, 远远低于传统算法; 2) 除求交外, 本文算法通过设计合理的数据结构, 使得弧段分割以及多边形自动生成基本上只有比较、调整指针等运算, 很少有几何计算, 大大减少了计算量, 从而极大地提高了算法的效率; 3) 针对多边形嵌套的情况, 本文算法的效率更稳定(传统左转算法的效率对于嵌套多边形的数量非常敏感, 算法的效率不稳定); 4) 算法更优美, 结构更紧凑, 与传统算法相比, 需要处理的特殊情况极少,
对于桥或者悬链的处理更稳定、更有效; 5) 基于模拟数据和实际数据的实验结果表明, 与传统算法和ARCGIS相比, 本文算法在数据处理效率上有明显的提升。
本文算法的另一个特点是易于扩展, 这是因为算法是以事件点的处理为核心。显然, 事件点(包括弧段端点、中间顶点和交点)是数据的骨架, 围绕数据的多种运算都可以围绕事件点展开。例如,除实现弧段分割和多边形构建外, 通过建立合理的数据结构, 容易实现结点和弧段、多边形和弧段等拓扑关系的输出; 通过适当的扩展, 还容易实现多边形集合的布尔运算等。另外, 对于空间数据, 点是最简单的, 也最容易实现并行计算和分布式计算,因此基于扫描的思想, 有望发展出更多的可以实现并行和分布式计算的空间数据处理算法, 从而有效地应对大数据时代对空间数据处理能力的挑战。
参考文献
[1] 周心铁. 地理信息系统中由线段数字化记录重建多边形的方法. 遥感学报, 1988, 3(1): 57‒64 [2] 陆守一, 韩兆云. 地理信息系统中自动获取多边形几何信息和拓扑信息算法的研究. 北京林业大学学
报, 1996, 18(增刊 2): 44‒47 [3] 杜清运. 地图数据库中多边形数据的自动组织. 测绘学报, 1989, 18(3): 204‒212 [4] 张军, 李玉祥. 地理信息系统中建立多边形拓扑关系的算法研究. 测绘科学, 1995, 20(2): 24‒28 [5] 陈春, 张树文. GIS中多边形图拓扑信息生成的数学基础. 测绘学报, 1996, 25(4): 266‒271 [6] 程双伟. GIS中拓扑关系的建立与更新[D]. 郑州:中国人民解放军信息工程大学, 2002 [7] Yang B, Purves R, Weibel R. Efficient transmission of vector data over the internet. International Journal of Geographical Information Science, 2007, 21(2): 215‒ 237 [8] 刘勖, 张镠亮, 宣国富. 基于有向弧的多边形拓扑关系生成算法. 计算机与信息技术, 2009, 17(5): 65‒67 [9] 周立新, 严静, 潘云鹤. 一个基于图的多边形拓扑关系生成算法. 计算机应用, 1999, 19(10): 37‒39 [10] 陈占龙, 吴亮, 刘焕焕. 基于图模型的多边形自动构建算法. 微电子学与计算机, 2011, 28(8): 143‒145 [11] 王杰臣. 多边形拓扑关系构建的栅格算法. 测绘学报, 2002, 31(3): 249‒254 [12] Mineter M J. A software framework to create vectortopology in parallel GIS operations. International Journal of Geographical Information Science, 2003, 17(3): 203‒222 [13] Žalik B. A topology construction from line drawings using a uniform plane subdivision technique. Computer-aided Design, 1999, 31(5): 335‒348 [14] Teng J, Huang W, Lou X. A new algorithm of construction of polygon topological relationships in geographic information systems // IEEE International Symposium on Geoscience and Remote Sensing. Denver, 2006: 2822‒2824 [15] 闫浩文, 杨维芳, 陈全功, 等. 基于方位角计算的拓扑多边形自动构建快速算法. 中国图象图形学报, 2000, 5(7): 563‒567 [16] Liang X W, Liu Z Q. An algorithm of polygon autoconstruction based on angle changing tendence. Journal of Image & Graphics, 2005, 6(6): 785‒789 [17] 李大军, 刘波, 赵宝贵, 等. 拓扑多边形自动构建的一种改进算法. 计算机工程与应用, 2005, 41(16): 80‒82 [18] 付胜博, 戴冠中. GIS中实时建立区域拓扑关系的快速算法. 信息安全与通信保密, 2006, 28(11): 66‒68 [19] Ling Y, Ouyang Y, Zhang M. Polygon auto-construction algorithm based on vector external product of virtual arc // International Symposium on Intelligent Information Technology Application Workshops. Shanghai, 2008: 24‒27 [20] Krivograd S, Trlep M, Žalik B. A rapid algorithm for topology construction from a set of line segments // Wseas International Conference on Telecommunications and Informatics. Istanbul, 2006: 133‒138 [21] 卢浩, 钟耳顺, 王天宝, 等. 一种基于拓扑信息的多边形数据自动生成算法. 地理与地理信息科学, 2012, 28(4): 38‒41 [22] Vinh N N, Le B. Constructing and modeling parcel boundaries from a set of lines for querying adjacent spatial relationships // International Conference on Computational Science and Its Applications. Barcelona, 2013: 540‒549 [23] Jia Q, Che D, Xiu C. A fast algorithm of polygon auto-construction based on left-turn algorithm in coal mine geological mapping. Science of Surveying & Mapping, 2016, 41(12): 70‒74 [24] Bentley J L, Ottmann T A. Algorithms for reporting and counting geometric intersections. IEEE Transactions on computers, 1979, 28(9): 643‒647