基于GPU的任意多边形相交面积计算方法

高 艺,罗健欣,裘杭萍,唐 斌,吴 波

(1.解放军理工大学 指挥信息系统学院，江苏 南京 210007； 2.解放军61175部队，江苏 南京 210049)

引用著录：高艺,罗健欣,裘杭萍，等.基于GPU的任意多边形相交面积计算方法[J].测绘工程，2017，26(12)：55-59.

DOI:10.19349/j.cnki.issn1006-7949.2017.12.011

基于GPU的任意多边形相交面积计算方法

高 艺1,2,罗健欣1,裘杭萍1,唐 斌1,吴 波1

(1.解放军理工大学 指挥信息系统学院，江苏 南京 210007； 2.解放军61175部队，江苏 南京 210049)

一直以来，任意多边形相交面积的高效计算都是地理信息系统中空间分析算法研究的重点。文中提出了一种基于GPU的栅格化多边形相交面积算法GPURAS，在此基础上，分别采用蒙特卡罗方法和遮挡查询技术进一步提出GPURASMC算法和GPURASQ算法，并证明了上述算法的正确性。实验对简单多边形、任意复杂多边形及大数据量多边形进行了测试对比，结果表明：GPURAS算法精度高，通用性较好但效率受CPU与GPU通信延迟的影响；GPURASMC算法效率较高但牺牲了部分精度；GPURASQ算法精度高、效率高但局限于特定运行环境。与基于CPU的传统算法相比，文中所提3种算法效率更高，在处理包含大量顶点的多边形时，效率提升尤为明显。

多边形; 相交面积计算; GPU; 栅格化; 蒙特卡罗

多边形指同一平面上不同的点首尾顺次连结，任一个顶点都在边上，且任意不相邻的两条边不相交。求任意多边形的相交面积是计算几何学、流体力学、计算机图形学的基本问题。在全国国土资源“一张图”基础数据服务系统中，经常涉及大量的多边形套合相交计算。因此，对这一问题的研究具有重要的理论意义和应用价值。

多边形相交面积的计算离不开多边形相交集的判定，这是计算相交面积的前提。对这一问题的讨论在多边形布尔运算中已有大量的研究成果，如处理凸多边形的经典算法有Shamos算法[1]和O’Rourke算法[2]；支持任意多边形裁剪的经典算法有Weiler算法[3]和Vatti算法[4]。朱雅音等[5]提出将多边形的边分为奇偶边的思想，根据边的拓扑关系来确定相交；刘永奎等[6]利用单链表结构减少交点计算的时间；候宝明等[7]对Weiler算法进行了改进，采用接缝技术解决了有空洞多边形的求交问题；Murta[8]实现了对Vatti算法的改进,增强了算法的稳定性。最近，齐东洲等[9]提出了一种新的改进算法，利用循环单链表避开复杂的出入点以计算提高算法的执行效率。

上述算法在一定程度上提高了计算效率，降低了时间复杂度，但大多原理复杂，编程实现困难。许多学者注意到栅格数据结构比矢量数据结构简单且模型复杂度对数据量影响较小，提出了数据栅格化方法。范俊甫等[10]提出了一种基于栅格化思想的多边形裁剪算法；Zhou等[11]提出了多核CPU栅格化方法。然而由于CPU的串行性，各个栅格只能顺序处理，降低了处理速度。

随着图形处理硬件的飞速发展，GPU依靠其强大的计算能力越来越多应用于计算机图形学[12]、碰撞检测[13]等领域。与基于CPU的串行处理方式不同，GPU的并行特性减少了程序运算过程的时间代价，极大提升了程序的执行效率。赵斯思等[14]将GPU用于多边形叠加分析，有效地实现了算法的加速。

本文将GPU用于多边形的交集运算，提出了基于GPU的任意多边形判交及相交面积计算方法。算法将顶点表示的多边形转换成栅格表示的图像，在GPU中生成采样栅格，根据栅格位置标识符确定对象间是否相交，利用相交栅格数计算相交面积。算法不需对多边形的凹凸性做任何假设，能工作在任意可光栅化的几何图元上。

1 多边形相交面积的CPU计算方法

现有的面积计算方法分为两类，一类是确定性算法，计算结果是绝对精确值；另一类是数值计算方法，计算结果是近似值，Monte Carlo方法是其中一项非常重要的数值计算方法。

问题描述：算法的输入是任意两个多边形P{P1P2…PK}和Q{Q1Q2…QL},顶点坐标分别为Pi(xi,yi)和Qi(xi,yi),顶点按逆时针方向排列，输出P和Q的相交面积为S。

1.1 确定性算法

算法步骤：

1)计算两个多边形每条边之间的交点。判断线段之间的位置关系并计算交点。若两个多边形分别有n条边和m条边，则两个多边形需要计算边数的时间复杂性为O(nm)。

2)判断包含在多边形内部的点，即P(Q)的顶点是否在Q(P)的内部。判断点是否在多边形内部的方法有很多，如射线法、环绕数法及改进的射线扫描法，时间复杂度为O(nlogn)，n为多边形的顶点数。

3)将交点和多边形内部的点，按逆时针(或顺时针)排序，得出最终的点集。

4)将多边形分割成多个三角形，求三角形面积代数和。三角形的面积由3个顶点构成的两个平面向量的外积求得。

1.2MonteCarlo算法

算法步骤：

1)获取多边形P和多边形Q的顶点数据及多边形的边的数量，并根据坐标值求得这两个多边形分别在x方向最大和最小值right和left及y方向最大最小值top和bottom。

2)生成m个随机顶点，通过随机顶点的坐标值判断其与两个多边形的包含关系。

3)统计随机点落在两个多边形相交部分的个数count。

4)计算相交面积S=(count/m)×Smax，(Smax=(right-left)×(top-bottom))。

上述两种算法各有优缺点，确定性算法计算结果精确，但是大量的跨立实验及点包含测试实现复杂，当多边形数量较多时计算量较大；Monte Carlo算法原理简单易实现，但需对大量的随机点做点包含测试，算法的整体性能受限于多边形的数量。

2 多边形相交面积的GPU计算方法

2.1 基于GPU的栅格化算法GPURAS

多边形栅格化首先要根据多边形集合中x和y方向的坐标极值来确定栅格化图幅的范围。图幅宽度TextureWidth是顶点坐标x方向的坐标极值[Vxmin，Vxmax]，图幅高度TextureHeight是顶点坐标y方向的坐标极值[Vymin，Vymax]，由此确定栅格图幅面积STexArea。

GPURAS算法按照从P到Q的顺序栅格化多边形，采用两路模板测试来判交。一路测试初始化缓冲区模板值为0.0，测试对象标记为“完全可见”，此时P在相应模板位置的模板值将变为1.0，其余位置保持初始值0.0不变；执行二路测试时将测试对象标记为“大于等于”，具体而言，若Q与P相交，则其在缓冲区模板位置的模板值将变为2.0，当且仅当大于等于2.0的模板值才能通过二路模板测试，此时对象间处于相交状态。使用glReadPixels()对缓冲区的矩形栅格区域进行读取，记录大于等于2.0的模板值数量count。

形式化描述如下：分别用I(P)和I(Q)来表示P和Q栅格化后在栅格图幅中的栅格内域。一路测试初始化模板值F=0.0，栅格化P，I(P)的模板值F=1.0；二路测试栅格化Q，若P和Q相交，则相交部分模板值F=2.0。用数学方式表达为：若∃P∩Q⟺FI(P)∩I(Q)=2.0。

相交面积S通过像素比与栅格场面积求得。相交面积计算公式为

(1)

其中，RES是栅格分辨率。从上式可知算法的复杂度只与栅格分辨率RES有关，与多边形的顶点数或边数无关。

上述讨论了两个多边形相交的一般情况，但两个多边形相交还有几种特殊情况，如图1所示。对于多边形的特殊位置关系，算法采用设置相交面积阈值[0.0,0.1]的方式来处理交集为空集的情况,当相交面积值在该区间则为空。

1)两个完全分离的多边形，如图1(a)所示，交集为空集。

2)两个多边形分离，但部分边重叠,如图1(b)、图1(c)、图1(d)所示，交集为空集。

3)两个多边形分离，但共点，如图1(e)所示，交集为空集。

4)两个多边形相互包含，其中一个多边形在另一个多边形内部，如图1(f)所示，交集为内层多边形。

图1 两个多边形之间的特殊位置关系

以上讨论了两个多边形相交的判断，下面对有n个多边形的相交判断进行说明。

若存在n个多边形的相交判断，那么在最坏情况下需要(n-1)+(n-2)+…+1=n(n-1)/2=O(n2)次相互测试，时间复杂度为O(n2)。如果能减少相互测试所产生的消耗，将会线性降低测试运行时间。因此，本文将大量多边形相交测试的过程分为两个阶段：粗粒度阶段和细粒度阶段。

给出粗粒度阶段的快速排除策略：若任意一点R对应的栅格标示符F为0.0，且它的8邻域栅格的F全为0.0，则该点一定不在多边形内部。

如图2(a)所示，对于10个多边形对象，传统的算法需要45次单独的对组合测试，而经过本文粗粒度阶段后，可以分成5个独立的子集。如图2(b)(虚线框区域)所示，随后的细粒度阶段只需在5个子集分组中执行上述算法，即只需10次独立的多边形对两两测试，极大地减少了计算工作量。

图2 粗粒度阶段确定分离子集

2.2 基于GPU的Monte Carlo栅格化算法GPURASMC

GPURAS算法需用glReadPixels()将存储在GPU内存缓冲区中的数据读回CPU，这将造成CPU和GPU之间的通信延迟。为解决这一问题，提出了一些改进算法，比如直方图查询，双线性过滤等，目的都是为了减少缓冲区数据的读取量，从而降低缓冲区的读取次数。

本节在GPURAS算法的基础上，进行Monte Carlo取样，提出基于GPU的Monte Carlo栅格化算法(GPURASMC)，用少量随机栅格子块模拟整个栅格场以减少缓冲区的数据读取量。

按照栅格图幅分辨率的大小，生成行、列随机数m，对每一个随机选取的图像子块读取大于等于2.0的模板值，统计数量count'，点位于两相交部分的概率为S/STexArea。给出相交面积公式如下：

(2)

证明：不失一般性，设STexArea为栅格场面积，S为两多边形相交部分的面积。

设事件A：投掷一次，并投在相交区域内。

因为投掷点服从二维均匀分布，所以p(A)=S/STexArea，设count′是m次投掷中投在相交区域内的次数，即栅格位置标识符大于等于2.0的栅格数，ε>0为任意正数，根据伯努利大数定律

复杂性：算法并不存储多边形的顶点，存储的仅是多边形所离散化的栅格范围，从式(1)和式(2)可知，其复杂度只与分辨率有关O(RES)。

2.3 遮挡查询的GPU栅格化算法GPURASQ

GPURASMC算法通过减少缓冲区的读取量降低了CPU和GPU之间的传输时间，但该查询工作仍是以读取缓冲区的方式实现。采用遮挡查询(occlusionquery)技术可以完全避免读取GPU缓冲区。

基于GPU遮挡查询的基本思想是判断一个物体的可见性，将物体的包围盒送入GPU栅格化后，其像素与预先设定的模板阈值比较，再由GPU返回可见像素的数量。基于这一思想本文提出利用遮挡查询的GPU栅格化算法(GPURASQ)。

OPENGL1.5及后续版本以及OpenGL ES 3.0都支持这一特性。ARB_occlusion_query扩展执行GPU遮挡查询的命令，其查询过程是由GPU来确定最终在屏幕上可见像素的数量。

在两遍模板测试前开启遮挡查询功能，通过glGetQueryObjectuiv()对通过二路渲染的图元中的采样数量进行计数，自动返回采样数量。

3 实验与分析

实验使用C++与GLSL语言，在Microsoft Windows 7操作系统，OpenGL 4.4.0上实现本文算法，实验硬件为Inter®Core(TM)i5-3337U处理器，4G内存，显卡为NVIDIA GeForce GT 620M。实验结果及分析如下：

3.1 算法误差分析

栅格化处理本身是一个有损处理过程，不可避免地存在误差[15]，因此必须对栅格化处理后的面积误差进行分析。相对面积误差E的定义如下：

%.

(3)

式中:S为分别采用本文3种算法得到的相交面积值；SCPUD为采用基于CPU的确定性算法得到的相交面积值。

实验数据采用随机产生的两个简单多边形，确定性算法得到的相交面积为5.5 m2，在相同的硬件条件下，本文3种算法的实验结果如表1所示。在256×256分辨率下，相交面积依次为5.525 0 m2，5.656 0 m2, 5.525 0 m2，相对面积误差分别为0.45%，2.84%，0.45%；1 024×1 024分辨率下相交面积分别为5.502 9 m2，5,496 3 m2和5.502 9 m2，相对面积误差分别为0.05%，0.067%，0.05%。本文算法的误差随着分辨率的提高呈现下降的趋势，这样的舍入误差在很多工程和一些大型软件的可接受范围内。

表1 本文算法与确定性算法相交面积误差统计结果

3.2 算法效率分析

本文3种算法与基于CPU的Monte Carlo算法同为数值解，因此,在相同的硬件实验条件下进行效率对比。实验数据采用两个有“孔洞”的复杂多边形，实验结果如表2所示。基于CPU的Monte Carlo算法在表2中简称CPUMC,在随机点数～分辨率基本相同的数量级下进行效率对比。

256×256～104时, CPUMC算法的计算时间大约为0.074 6 s，本文3种算法的计算时间分别为0.004 s，0.003 s，0.001 s；数量级增加3倍(1 024×1 024～106)时，CPUMC为5.147 s，本文3种算法的计算时间分别为0.032 s， 0.024 s，0.002 s。

表2 本文算法与CPUMC算法运行效率统计结果 s

将本文算法的分辨率分别固定在256×256，CPUMC取104个随机数和分辨率固定在1 024×1 024，CPUMC取106个随机数，比较具有相同形状但不同顶点数量的多边形间的求交计算时间，对比两种算法的效率差异。如表3和表4所示，其中N为两个多边形的顶点数之和。

表3 本文算法256×256分辨率与CPUMC算法104个随机数的运行效率统计结果

表4 本文算法1 024×1 024分辨率与CPUMC算法106个随机数的运行效率统计结果

当多边形所包含的顶点数量较少时, CPUMC算法具有较高的计算效率,随着多边形顶点数量的增加, CPUMC算法的计算效率快速降低,上万个点的点包含测试造成了性能的急剧下降，106个随机数时，经过十多分钟仍未得到结果。本文算法在处理包含大量顶点的多边形时更有优势。

4 结束语

本文提出了基于GPU栅格化的多边形相交面积计算算法。通过实验，对比和分析了本文算法与基于CPU确定值算法的相对面积误差率，以及与基于CPU的Monte Carlo算法的执行效率。实验结果表明，本文算法的计算效率不受多边形顶点数量的影响，仅对分辨率敏感，随着分辨率的增加，误差下降但时间开销有所增加。基于CPU的确定值算法和Monte Carlo算法在处理小数据时较为高效；而本文算法是一种以精度的适当牺牲换取优异综合性能的算法，在处理复杂多边形及包含大量顶点的多边形数据时具有优势。

[1] PREPARATA F P, SHAMOS M I. Computational geometry: an introduction[M]. Berlin: Springer,1985.

[2] O’ROURKE J. Computational geometry in c [M]. New York: Cambridge University Press,1994.

[3] WEILER K, ATHERTON P. Hidden surface removal using polygon area sorting[C].Proceedings of the 4th Annual Conference on Computer Graphic and Interactive Techniques. New York: ACM Press,1977:214-222.

[4] VATTI B R. A generic solution to polygon clipping [J].Communications of the ACM,1992,35(7):56-63.

[5] 朱雅音，王化文，万丰，等.确定两个任意简单多边形交、并、差的算法[J].计算机研究与发展，2003,40(4): 576-583.

[6] 刘永奎，高云，黄有群.一个有效的多边形裁剪算法[J].软件学报，2003,14(4): 845-856.

[7] 候宝明，刘雪娜.任意多边形区域交的有效算法[J].计算机辅助工程，2009,18(2):73-76.

[8] MURTA A. A General Polygon Clipping Library[EB/OL].http://www.cs.man.ac.uk/-toby/alan/software/gpc.html.2013-11-01.

[9] 齐东洲，吴敏.高效的多边形布尔计算方法[J].计算机应用,2014,34(增2):78-82.

[10] 范俊甫，孔维华，马廷，等.RaPC：一种基于栅格化思想的多边形裁剪算法及其误差分析[J].测绘学报,2015，44(3):338-345.

[11] ZHOU Chen, CHEN Zhenjie, LIU Yongxue, et al. A strategy for parallelising polygon rasterisation algorithms using multi-core CPUs [J]. Journal of Spatial Science. 2016. 47-48.

[12] LAINE S, KARRAS T. Efficient sparse voxel octrees[C]. ACM SIGGRAPH Symposium on Interactive 3D Graphics and Games (I3D),2011.17(8),1048-1059.

[13] CHENTANEZ N, MÜLLER M,MACKLIN M. GPU accelerated grid-free surface tracking[J]. Computers & Graphics. 2016,57, 1-11.

[14] 赵斯思，周成虎.GPU加速的多边形叠加分析[J].地理科学进展.2013,32(1):114-120.

[15] BAI Yan,LIAO Shunbao,SUN Jiulin.Scale effect and methods for accuracy evaluation of attribute information Loss in rasterization[J].Journal of Geographical Sciences,2011,21(6):1089-1100.

[责任编辑：刘文霞]

Anefficientalgorithmofarbitrarypolygonsintersectionarea

GAO Yi1,2, LUO Jianxin1, QIU Hangping1, TANG Bin1，WU Bo1

(1.College of Command Information System, PLA University of Science and Technology, Nanjing 210007, China;2. PLA Troops 61175, Nanjing 210049, China)

For a long time, the fast processing the intersection area of arbitrary polygons has been an important research topic in spatial analysis algorithm of geographic information system. A GPU-based rasterized polygon intersection area algorithm is proposed as GPURAS. And its accelerated versions: GPURASMC and GPURASQ, which use Monte Carlo method and occlusion query technique respectively, are also presented and the correctness of these algorithms is proved. Experiments and comparisons are performed using simple, arbitrary complex, and large-scale polygons. The result shows the GPURAS algorithm has better universality but the efficiency is affected by CPU and GPU communication delay. The GPURASMC algorithm has higher efficiency but compromises a bit on accuracy. The GPURASQ algorithm has excellent accuracy and the highest efficiency, but its usage is limited by the running environment. The average efficiency of all the three proposed algorithms is higher than their CPU counterparts, especially in handling large-scale polygons.

polygon; intersection area; GPU; rasterization; Monte Carlo

P208

A

1006-7949(2017)12-0055-05

2016-11-01

国家863计划资助项目(2012AA01A509);江苏省青年科学基金资助项目(BK20150722)

高 艺(1982-),女,工程师,博士研究生.