基于改进边界代数法的相交多边形提取方法

摘要

本发明涉及一种基于改进边界代数法的相交多边形提取方法，包括以下步骤：对所有图层中的多边形顺序进行编号；计算包含所有图层的MBR，数组

说明书

技术领域

本发明涉及一种基于改进边界代数法的相交多边形提取方法，属于图像处理技术领域。

背景技术

多边形相交计算是地理信息系统(Geographical Information System,GIS)中的重要计算类型，被广泛地应用于多边形拓扑检查、缓冲区生成、叠置分析等GIS空间分析类型，具有典型的计算复杂、密集的算法特征。该类型空间计算的一般步骤包括：(1)遍历所有多边形数据集并搜索、确定相交多边形；(2)提取相交多边形并完成对相交部分结果的求解。目前已有大量学者针对该类型算法的处理效率和精度展开研究，根据其使用方法的不同，可分为基于矢量的方法和基于栅格的方法两种类型。

基于矢量的方法主要面向空间对象描述，各多边形对象之间的空间关系判断需要基于计算几何算法进行定位、分析和检索。其基本处理过程包括：利用多边形的最小外接矩形(Minimal Bounding Rectangle,简称MBR)表征多边形的空间形状，首先过滤出可能相交的候选多边形，进而对候选多边形进行精炼确定相交多边形组，最后对相交多边形组进行裁剪。尽管不同的多边形空间计算算法原理不同，但具有相同的算法特征，即本质操作均是过滤出具有相交关系的多边形，并对多边形点集集合进行裁剪。因此，计算密集型多边形空间计算的核心问题在于多边形之间的裁剪及相交多边形的确定。现有成熟的多边形裁剪算法主要是针对两个多边形进行；针对海量多边形数据，通用的方法是采用“暴力”算法逐个遍历源数据判断当前多边形与其他多边形的空间关系，并反复调用裁剪算法实现相交多边形的裁剪，不但计算复杂度高、效率低下，而且难以适用于规模化的多边形数据集的相交计算。但利用空间索引只能减少后续多边形的判断次数，其形成的多边形组中仍然包含大量不相交多边形，因而该方法仍然会浪费大量的时间在不相交多边形的计算上，从而降低了计算效率。

栅格数据模型是与矢量数据模型是同等重要的空间对象表示方法。相对于矢量地理数据格式而言，栅格数据分布规则、独立性较高。因此采用基于栅格的计算方法可极大地提高空间计算效率。基于栅格的计算方法是在简单的栅格图像中通过栅格单元之间属性值的简单代数运算实现空间分析。该方法的一般计算过程包括三个步骤：(1)将矢量多边形图层分别进行栅格化填充，获得两个处理图层的栅格化结果；(2)对两个图层的栅格化图像根据其属性值进行栅格计算，以获得指定空间计算结果；(3)对计算结果进行栅格矢量化，提取多边形边界及其所属岛洞，并构建多边形拓扑关系，从而完成计算。目前，已有较多学者采用该方法实现对多边形数据的空间计算，但现有的方法往往只能处理两个或两个以上图层中的多边形的相交、合并等栅格计算，无法处理同一图层内多边形相互重叠的情形时的相交计算。

发明内容

本发明要解决技术问题是：克服上述技术的缺点，提供一种计算复杂度低、适用于规模化的多边形数据集的基于改进边界代数法的相交多边形提取方法。

为了解决上述技术问题，本发明提出的技术方案是：一种基于改进边界代数法的相交多边形提取方法，包括以下步骤：

步骤一、对所有图层中的多边形顺序进行编号，每个多边形对应唯一的多边形ID；

步骤二、计算包含所有图层的MBR，根据栅格尺寸，生成初始值为0的数组hDstDS以存放栅格单元的属性值，其大小为nYSize行×nXSize列，同时生成与数组hDstDS相同大小、初始值为-1的数组pIDArray以存放多边形ID；

生成一个新的数组RLEGroup以存放游程；其中，一个栅格行中栅格单元的属性值为pNum的游程定义为RLE{StartX，EndX，LocateY，(pID₀，pID₁，…，pID_pNum-1)，pNum}；其中，StartX和EndX为该栅格行中属性值为pNum的起始栅格单元与结束栅格单元的横坐标，LocateY为该栅格行的纵坐标，(pID₀，pID₁，…，pID_pNum-1)为存放相交多边形ID的数组pGroup

步骤三、对步骤一中的所有多边形使用边界代数算法依次进行栅格化，在栅格化过程中赋予各多边形的属性值均为1，并更新数组hDstDS的值为栅格单元的属性值；

在数组hDstDS中获取当前多边形MBR包含的栅格单元，并逐行读取，对于纵坐标为LocateY的栅格行中的任一栅格单元(StartX,LocateY)，若该栅格单元位于该多边形内部，则计算其在hDstDS中的存放位置locate＝LocateY×nXSize+StartX，并获取其属性值value＝hDstDS[locate]；

1)若该栅格单元的属性值value<2，则该栅格单元仅位于一个多边形内部，不予处理；

2)若该栅格单元的属性值value＝2，则该栅格单元位于两个相交多边形内部；在该栅格行中以当前栅格单元为起点继续扫描后续栅格单元，以形成从栅格单元(StartX，LocateY)至栅格单元(EndX，LocateY)的栅格序列，其中StartX≤EndX；上述栅格序列中的任一栅格单元(LocateX_i，LocateY)，StartX≤LocateX_i≤EndX，须位于该多边形内部，并满足以下条件：

这样，该栅格序列中的所有栅格单元均位于两个相交多边形内部；同时，这两个相交多边形ID分别为pIDArray[locate]和pCurrentID；创建一个新的游程以存放上述栅格序列，该游程可表示为{StartX,EndX,LocateY,(pIDArray[locate],pCurrentID),2}，并将该游程存放至数组RLEGroup中；

3)若该栅格单元的属性值value>2，则该栅格单元位于三个或三个以上相交多边形内部；在该栅格行中以当前栅格单元为起点继续扫描后续栅格单元，以形成从栅格单元(StartX,LocateY)至栅格单元(EndX,LocateY)的栅格序列，其中StartX≤EndX；上述栅格序列中的任一栅格单元(LocateX_i,LocateY)，其中StartX≤LocateX_i≤EndX，须位于多边形内部，并满足以下条件：

这样，该栅格序列中的所有栅格单元均位于value个相交多边形内部，其中一个多边形ID为pCurrentID，另外value-1个多边形ID存放于已有的游程中；

从栅格单元(max(StartX_j,StartX),LocateY)至栅格单元(min(EndX_j,EndX),LocateY)，搜索包含value-1个多边形ID的游程RLE_j{StartX_j,EndX_j,LocateY,(pID₀,pID₁,…,pID_value-2),value-1}，首先将其分解为两个新的游程，分别为{StartX_j,max(StartX_j,StartX)-1,LocateY,(pID₀,pID₁,…,pID_value-2),value-1}和{min(EndX_j,EndX)+1,EndX_j,LocateY,(pID₀,pID₁,…,pID_value-2),value-1}，其中min(EndX_j,EndX)≤EndX_j，StartX_j≤min(StartX_j,StartX)，并将两个新的游程放入数组RLEGroup中，然后将原游程RLE_j从数组RLEGroup中删除；同时生成一个新的游程{max(StartX_j,StartX),min(EndX_j,EndX),LocateY,(RLE_j.pGroup,pCurrentID)，value}，并将其放入数组RLEGroup中；

步骤五、更新数组pIDArray，将当前多边形ID值赋给数组pIDArray中位于当前多边形内部的栅格单元；

步骤六、从数组RLEGroup存储的游程中提取相应的相交多边形组，即每个游程中的数组pGroup即对应一个相交多边形组。

相较于矢量多边形数据，栅格数据具有规整的数据结构，栅格数据更适合用来进行GIS空间分析，并较多边形数据具有更高的计算效率。因此，将矢量多边形转换为栅格数据以快速搜索相交多边形可提高GIS空间分析效率。

边界代数法是一种基于积分思想的多边形栅格化方法，通过简单的加减代数运算将多边形属性值赋给多边形内部及边界上的栅格单元，从而实现多边形栅格化。矢量多边形通常由一个外环和若干个内环组成，边界代数法对任一多边形的处理过程包括以下步骤，如图1所示：(1)将覆盖该多边形的面域进行整体栅格化，并对栅格阵列的栅格单元值初始化为零；(2)由其边界上任意一点开始，对多边形外环逆时针方向搜索其边界，对内环顺时针方向搜索其边界；(3)当边界线段方向为上行时，对该线段左侧具有相同行坐标的所有栅格全部加上属性值a；当边界线段方向为下行时，对该线段左侧具有相同行坐标的所有栅格全部减去属性值a；当边界线段平行于栅格行时，则不做运算；(4)重复步骤(3)直至多边形所有边界处理完毕。在完成上述计算后，属性值为a的栅格单元即为本次多边形栅格化的计算结果。相较于其他多边形栅格化算法，边界代数法的优势在于其计算效率更高；更重要的是该算法可快速识别位于多边形相交部分的栅格单元。

本发明利用边界代数法的思想，赋予栅格单元的属性值固定为1，则在边界代数法处理过程中，若栅格单元属性值为1，则该栅格单元只位于一个多边形内部；若栅格单元属性值为n(n>1)，则该栅格单元位于n个多边形相交部分内，如图2所示。本发明利用该特性，实现了对多边形相交部分的栅格单元的快速定位。

附图说明

下面结合附图对本发明作进一步说明。

图1是现有技术中边界代数法的示意图。

图2是本发明的原理示意图。

图3是本发明实施例的流程示意图。

具体实施方式

实施例

如图1所示，本实施例的基于改进边界代数法的相交多边形提取方法，包括以下步骤：

步骤一、对所有图层中的多边形顺序进行编号，每个多边形对应唯一的多边形ID；图1所示共有三个多边形，三个多边形的ID分别是0，1和2。

步骤二、计算包含所有图层的MBR(即最小外接矩形)，根据栅格尺寸，生成初始值为0的数组hDstDS以存放栅格单元的属性值，其大小为nYSize行×nXSize列，同时生成与数组hDstDS相同大小、初始值为-1的数组pIDArray以存放多边形ID。本例中数组hDstDS和数组pIDArray的大小均为6×6。

生成一个新的数组RLEGroup以存放游程；其中，一个栅格行中栅格单元的属性值为pNum的游程定义为RLE{StartX，EndX，LocateY，(pID₀，pID₁，…，pID_pNum-1)，pNum}；其中，StartX和EndX为该栅格行中属性值为pNum的起始栅格单元与结束栅格单元的横坐标，LocateY为该栅格行的纵坐标，(pID₀，pID₁，…，pID_pNum-1)为存放相交多边形ID的数组pGroup。

步骤三、对步骤一中的所有多边形使用边界代数算法依次进行栅格化，在栅格化过程中赋予各多边形的属性值均为1，并更新数组hDstDS的值为栅格单元的属性值；

在数组hDstDS中获取当前多边形MBR包含的栅格单元，并逐行读取，对于纵坐标为LocateY的栅格行中的任一栅格单元(StartX,LocateY)，若该栅格单元位于该多边形内部，则计算其在hDstDS中的存放位置locate＝LocateY×nXSize+StartX，并获取其属性值value＝hDstDS[locate]；

1)若该栅格单元的属性值value<2，则该栅格单元仅位于一个多边形内部，不予处理；

2)若该栅格单元的属性值value＝2，则该栅格单元位于两个相交多边形内部；在该栅格行中以当前栅格单元为起点继续扫描后续栅格单元，以形成从栅格单元(StartX，LocateY)至栅格单元(EndX，LocateY)的栅格序列，其中StartX≤EndX；上述栅格序列中的任一栅格单元(LocateX_i，LocateY)，StartX≤LocateX_i≤EndX，须位于该多边形内部，并满足以下条件：

这样，该栅格序列中的所有栅格单元均位于两个相交多边形内部；同时，这两个相交多边形ID分别为pIDArray[locate]和pCurrentID；创建一个新的游程以存放上述栅格序列，该游程可表示为{StartX,EndX,LocateY,(pIDArray[locate],pCurrentID),2}，并将该游程存放至数组RLEGroup中；

在对一个多边形的栅格化过程中，其MBR中被赋值的栅格单元即位于该多边形内部，即，判断栅格单元位于多边形内部的过程是与栅格化过程同时进行，不再需要额外的计算过程。3)若该栅格单元的属性值value>2，则该栅格单元位于三个或三个以上相交多边形内部；在该栅格行中以当前栅格单元为起点继续扫描后续栅格单元，以形成从栅格单元(StartX,LocateY)至栅格单元(EndX,LocateY)的栅格序列，其中StartX≤EndX；上述栅格序列中的任一栅格单元(LocateX_i,LocateY)，其中StartX≤LocateX_i≤EndX，须位于多边形内部，并满足以下条件：

这样，该栅格序列中的所有栅格单元均位于value个相交多边形内部，其中一个多边形ID为pCurrentID，另外value-1个多边形ID存放于已有的游程中；

从栅格单元(max(StartX_j,StartX),LocateY)至栅格单元(min(EndX_j,EndX),LocateY)，搜索包含value-1个多边形ID的游程RLE_j{StartX_j,EndX_j,LocateY,(pID₀,pID₁,…,pID_value-2),value-1}，首先将其分解为两个新的游程，分别为{StartX_j,max(StartX_j,StartX)-1,LocateY,(pID₀,pID₁,…,pID_value-2),value-1}和{min(EndX_j,EndX)+1,EndX_j,LocateY,(pID₀,pID₁,…,pID_value-2),value-1}，其中min(EndX_j,EndX)≤EndX_j，StartX_j≤min(StartX_j,StartX)，并将两个新的游程放入数组RLEGroup中，然后将原游程RLE_j从数组RLEGroup中删除；同时生成一个新的游程{max(StartX_j,StartX),min(EndX_j,EndX),LocateY,(RLE_j.pGroup,pCurrentID)，value}，并将其放入数组RLEGroup中；

图3(a)所示为栅格化多边形0，此时尚无相交多边形。图3(b)所示为栅格化多边形1，此时在多边形1的MBR中有两个栅格行中有多边形0和多边形1相交的栅格序列，通过游程RLE₀和RLE₁进行记录。图3(c)所示为栅格化多边形2，此时在多边形2的MBR中即存在有相交多边形的栅格序列，通过游程RLE₂、sub-RLE₁、RLE₃、RLE₄和RLE₅进行记录。

考虑到不同的游程可能存放重复的多边形组，因此在相交多边形组提取完成后需要删除该类型的多边形组，以避免后续重复计算。当两个多边形组存储个数相同、ID相同的多边形时，其中一个多边形组应当予以删除；当一个多边形组存储的多边形为另一个多边形组的子集时，该多边形组也应当予以删除。在形成的多边形组中，若一个多边形组包含两个多边形，则这两个多边形相交；若一个多边形组包含两个以上的多边形，则每两个多边形之间均相交。

步骤五、更新数组pIDArray，将当前多边形ID值赋给数组pIDArray中位于当前多边形内部的栅格单元。

步骤六、从数组RLEGroup存储的游程中提取相应的相交多边形组，即每个游程中的数组pGroup即对应一个相交多边形组，如图3(d)所示。当相互独立的多边形组提取完成后，即可对各多边形组实现相交结果的计算，以获得相交结果多边形。

为了验证本实施例的提取方法，发明人采用包含2个真实土地利用数据和1个模拟数据进行提取相交多形边组后进行相交运算，并与传统的传统空间索引方法进行相交计算作了对比。其中，数据1为上海市2009年土地利用现状数据，其数据量为1.04GB、多边形个数为1,371,765；数据2为上海市县级行政区划数据，其数据量为376KB、多边形个数为17。数据1和数据2图层内均无相交多边形。数据3为模拟数据，主要通过移动、旋转数据1中一定数量的多边形，从而形成数据3中的相交多边形。数据3数据量为1.04GB，多边形个数为1,371,765，包含163,934个多边形组。

验证时采用高性能IBM计算节点作为计算工具，其硬件配置为：CPU 2颗，规格为Intel(R)Xeon(R)CPU E5-2620(主频2.00GHz，六核十二线程)；内存为16GB(4根4GB内存条，规格为DDR3 RDIMM1600MHz)；硬盘为2TB(2个1TB硬盘，规格为6Gbps 2.5”7.2Krpm NLSAS)，网络为集成的双口千兆以太网。软件配置：操作系统为Centos Linux 6.3，文件系统为Lustre系统。矢量与栅格数据的读写采用开源地理空间数据转换库GDAL(Geospatial DataAbstraction Library)实现；多边形组空间几何计算通过开源算法库GEOS(GeometryEngine Open Source)实现。

结果表明，本实施例的方法较传统空间索引方法的运行时间减少50％左右，即本实施例的方法较传统方法计算效率更高。

本发明不局限于上述实施例所述的具体技术方案，除上述实施例外，本发明还可以有其他实施方式。对于本领域的技术人员来说，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等形成的技术方案，均应包含在本发明的保护范围之内。