基于GPU的任意多邊形相交面積計算方法

高 藝,羅健欣,裘杭萍,唐 斌,吳 波

(1.解放軍理工大學 指揮信息系統(tǒng)學院，江蘇 南京 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.解放軍理工大學 指揮信息系統(tǒng)學院，江蘇 南京 210007； 2.解放軍61175部隊，江蘇 南京 210049)

一直以來，任意多邊形相交面積的高效計算都是地理信息系統(tǒng)中空間分析算法研究的重點。文中提出了一種基于GPU的柵格化多邊形相交面積算法GPURAS，在此基礎上，分別采用蒙特卡羅方法和遮擋查詢技術進一步提出GPURASMC算法和GPURASQ算法，并證明了上述算法的正確性。實驗對簡單多邊形、任意復雜多邊形及大數(shù)據(jù)量多邊形進行了測試對比，結果表明：GPURAS算法精度高，通用性較好但效率受CPU與GPU通信延遲的影響；GPURASMC算法效率較高但犧牲了部分精度；GPURASQ算法精度高、效率高但局限于特定運行環(huán)境。與基于CPU的傳統(tǒng)算法相比，文中所提3種算法效率更高，在處理包含大量頂點的多邊形時，效率提升尤為明顯。

多邊形; 相交面積計算; GPU; 柵格化; 蒙特卡羅

多邊形指同一平面上不同的點首尾順次連結，任一個頂點都在邊上，且任意不相鄰的兩條邊不相交。求任意多邊形的相交面積是計算幾何學、流體力學、計算機圖形學的基本問題。在全國國土資源“一張圖”基礎數(shù)據(jù)服務系統(tǒng)中，經(jīng)常涉及大量的多邊形套合相交計算。因此，對這一問題的研究具有重要的理論意義和應用價值。

多邊形相交面積的計算離不開多邊形相交集的判定，這是計算相交面積的前提。對這一問題的討論在多邊形布爾運算中已有大量的研究成果，如處理凸多邊形的經(jīng)典算法有Shamos算法[1]和O’Rourke算法[2]；支持任意多邊形裁剪的經(jīng)典算法有Weiler算法[3]和Vatti算法[4]。朱雅音等[5]提出將多邊形的邊分為奇偶邊的思想，根據(jù)邊的拓撲關系來確定相交；劉永奎等[6]利用單鏈表結構減少交點計算的時間；候寶明等[7]對Weiler算法進行了改進，采用接縫技術解決了有空洞多邊形的求交問題；Murta[8]實現(xiàn)了對Vatti算法的改進,增強了算法的穩(wěn)定性。最近，齊東洲等[9]提出了一種新的改進算法，利用循環(huán)單鏈表避開復雜的出入點以計算提高算法的執(zhí)行效率。

上述算法在一定程度上提高了計算效率，降低了時間復雜度，但大多原理復雜，編程實現(xiàn)困難。許多學者注意到柵格數(shù)據(jù)結構比矢量數(shù)據(jù)結構簡單且模型復雜度對數(shù)據(jù)量影響較小，提出了數(shù)據(jù)柵格化方法。范俊甫等[10]提出了一種基于柵格化思想的多邊形裁剪算法；Zhou等[11]提出了多核CPU柵格化方法。然而由于CPU的串行性，各個柵格只能順序處理，降低了處理速度。

隨著圖形處理硬件的飛速發(fā)展，GPU依靠其強大的計算能力越來越多應用于計算機圖形學[12]、碰撞檢測[13]等領域。與基于CPU的串行處理方式不同，GPU的并行特性減少了程序運算過程的時間代價，極大提升了程序的執(zhí)行效率。趙斯思等[14]將GPU用于多邊形疊加分析，有效地實現(xiàn)了算法的加速。

本文將GPU用于多邊形的交集運算，提出了基于GPU的任意多邊形判交及相交面積計算方法。算法將頂點表示的多邊形轉換成柵格表示的圖像，在GPU中生成采樣柵格，根據(jù)柵格位置標識符確定對象間是否相交，利用相交柵格數(shù)計算相交面積。算法不需對多邊形的凹凸性做任何假設，能工作在任意可光柵化的幾何圖元上。

1 多邊形相交面積的CPU計算方法

現(xiàn)有的面積計算方法分為兩類，一類是確定性算法，計算結果是絕對精確值；另一類是數(shù)值計算方法，計算結果是近似值，Monte Carlo方法是其中一項非常重要的數(shù)值計算方法。

問題描述：算法的輸入是任意兩個多邊形P{P1P2…PK}和Q{Q1Q2…QL},頂點坐標分別為Pi(xi,yi)和Qi(xi,yi),頂點按逆時針方向排列，輸出P和Q的相交面積為S。

1.1 確定性算法

算法步驟：

1)計算兩個多邊形每條邊之間的交點。判斷線段之間的位置關系并計算交點。若兩個多邊形分別有n條邊和m條邊，則兩個多邊形需要計算邊數(shù)的時間復雜性為O(nm)。

2)判斷包含在多邊形內(nèi)部的點，即P(Q)的頂點是否在Q(P)的內(nèi)部。判斷點是否在多邊形內(nèi)部的方法有很多，如射線法、環(huán)繞數(shù)法及改進的射線掃描法，時間復雜度為O(nlogn)，n為多邊形的頂點數(shù)。

3)將交點和多邊形內(nèi)部的點，按逆時針(或順時針)排序，得出最終的點集。

4)將多邊形分割成多個三角形，求三角形面積代數(shù)和。三角形的面積由3個頂點構成的兩個平面向量的外積求得。

1.2MonteCarlo算法

算法步驟：

1)獲取多邊形P和多邊形Q的頂點數(shù)據(jù)及多邊形的邊的數(shù)量，并根據(jù)坐標值求得這兩個多邊形分別在x方向最大和最小值right和left及y方向最大最小值top和bottom。

2)生成m個隨機頂點，通過隨機頂點的坐標值判斷其與兩個多邊形的包含關系。

3)統(tǒng)計隨機點落在兩個多邊形相交部分的個數(shù)count。

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

上述兩種算法各有優(yōu)缺點，確定性算法計算結果精確，但是大量的跨立實驗及點包含測試實現(xiàn)復雜，當多邊形數(shù)量較多時計算量較大；Monte Carlo算法原理簡單易實現(xiàn)，但需對大量的隨機點做點包含測試，算法的整體性能受限于多邊形的數(shù)量。

2 多邊形相交面積的GPU計算方法

2.1 基于GPU的柵格化算法GPURAS

多邊形柵格化首先要根據(jù)多邊形集合中x和y方向的坐標極值來確定柵格化圖幅的范圍。圖幅寬度TextureWidth是頂點坐標x方向的坐標極值[Vxmin，Vxmax]，圖幅高度TextureHeight是頂點坐標y方向的坐標極值[Vymin，Vymax]，由此確定柵格圖幅面積STexArea。

GPURAS算法按照從P到Q的順序柵格化多邊形，采用兩路模板測試來判交。一路測試初始化緩沖區(qū)模板值為0.0，測試對象標記為“完全可見”，此時P在相應模板位置的模板值將變?yōu)?.0，其余位置保持初始值0.0不變；執(zhí)行二路測試時將測試對象標記為“大于等于”，具體而言，若Q與P相交，則其在緩沖區(qū)模板位置的模板值將變?yōu)?.0，當且僅當大于等于2.0的模板值才能通過二路模板測試，此時對象間處于相交狀態(tài)。使用glReadPixels()對緩沖區(qū)的矩形柵格區(qū)域進行讀取，記錄大于等于2.0的模板值數(shù)量count。

形式化描述如下：分別用I(P)和I(Q)來表示P和Q柵格化后在柵格圖幅中的柵格內(nèi)域。一路測試初始化模板值F=0.0，柵格化P，I(P)的模板值F=1.0；二路測試柵格化Q，若P和Q相交，則相交部分模板值F=2.0。用數(shù)學方式表達為：若?P∩Q?FI(P)∩I(Q)=2.0。

相交面積S通過像素比與柵格場面積求得。相交面積計算公式為

(1)

其中，RES是柵格分辨率。從上式可知算法的復雜度只與柵格分辨率RES有關，與多邊形的頂點數(shù)或邊數(shù)無關。

上述討論了兩個多邊形相交的一般情況，但兩個多邊形相交還有幾種特殊情況，如圖1所示。對于多邊形的特殊位置關系，算法采用設置相交面積閾值[0.0,0.1]的方式來處理交集為空集的情況,當相交面積值在該區(qū)間則為空。

1)兩個完全分離的多邊形，如圖1(a)所示，交集為空集。

2)兩個多邊形分離，但部分邊重疊,如圖1(b)、圖1(c)、圖1(d)所示，交集為空集。

3)兩個多邊形分離，但共點，如圖1(e)所示，交集為空集。

4)兩個多邊形相互包含，其中一個多邊形在另一個多邊形內(nèi)部，如圖1(f)所示，交集為內(nèi)層多邊形。

圖1 兩個多邊形之間的特殊位置關系

以上討論了兩個多邊形相交的判斷，下面對有n個多邊形的相交判斷進行說明。

若存在n個多邊形的相交判斷，那么在最壞情況下需要(n-1)+(n-2)+…+1=n(n-1)/2=O(n2)次相互測試，時間復雜度為O(n2)。如果能減少相互測試所產(chǎn)生的消耗，將會線性降低測試運行時間。因此，本文將大量多邊形相交測試的過程分為兩個階段：粗粒度階段和細粒度階段。

給出粗粒度階段的快速排除策略：若任意一點R對應的柵格標示符F為0.0，且它的8鄰域柵格的F全為0.0，則該點一定不在多邊形內(nèi)部。

如圖2(a)所示，對于10個多邊形對象，傳統(tǒng)的算法需要45次單獨的對組合測試，而經(jīng)過本文粗粒度階段后，可以分成5個獨立的子集。如圖2(b)(虛線框區(qū)域)所示，隨后的細粒度階段只需在5個子集分組中執(zhí)行上述算法，即只需10次獨立的多邊形對兩兩測試，極大地減少了計算工作量。

圖2 粗粒度階段確定分離子集

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

GPURAS算法需用glReadPixels()將存儲在GPU內(nèi)存緩沖區(qū)中的數(shù)據(jù)讀回CPU，這將造成CPU和GPU之間的通信延遲。為解決這一問題，提出了一些改進算法，比如直方圖查詢，雙線性過濾等，目的都是為了減少緩沖區(qū)數(shù)據(jù)的讀取量，從而降低緩沖區(qū)的讀取次數(shù)。

本節(jié)在GPURAS算法的基礎上，進行Monte Carlo取樣，提出基于GPU的Monte Carlo柵格化算法(GPURASMC)，用少量隨機柵格子塊模擬整個柵格場以減少緩沖區(qū)的數(shù)據(jù)讀取量。

按照柵格圖幅分辨率的大小，生成行、列隨機數(shù)m，對每一個隨機選取的圖像子塊讀取大于等于2.0的模板值，統(tǒng)計數(shù)量count'，點位于兩相交部分的概率為S/STexArea。給出相交面積公式如下：

(2)

證明：不失一般性，設STexArea為柵格場面積，S為兩多邊形相交部分的面積。

設事件A：投擲一次，并投在相交區(qū)域內(nèi)。

因為投擲點服從二維均勻分布，所以p(A)=S/STexArea，設count′是m次投擲中投在相交區(qū)域內(nèi)的次數(shù)，即柵格位置標識符大于等于2.0的柵格數(shù)，ε>0為任意正數(shù)，根據(jù)伯努利大數(shù)定律

復雜性：算法并不存儲多邊形的頂點，存儲的僅是多邊形所離散化的柵格范圍，從式(1)和式(2)可知，其復雜度只與分辨率有關O(RES)。

2.3 遮擋查詢的GPU柵格化算法GPURASQ

GPURASMC算法通過減少緩沖區(qū)的讀取量降低了CPU和GPU之間的傳輸時間，但該查詢工作仍是以讀取緩沖區(qū)的方式實現(xiàn)。采用遮擋查詢(occlusionquery)技術可以完全避免讀取GPU緩沖區(qū)。

基于GPU遮擋查詢的基本思想是判斷一個物體的可見性，將物體的包圍盒送入GPU柵格化后，其像素與預先設定的模板閾值比較，再由GPU返回可見像素的數(shù)量?；谶@一思想本文提出利用遮擋查詢的GPU柵格化算法(GPURASQ)。

OPENGL1.5及后續(xù)版本以及OpenGL ES 3.0都支持這一特性。ARB_occlusion_query擴展執(zhí)行GPU遮擋查詢的命令，其查詢過程是由GPU來確定最終在屏幕上可見像素的數(shù)量。

在兩遍模板測試前開啟遮擋查詢功能，通過glGetQueryObjectuiv()對通過二路渲染的圖元中的采樣數(shù)量進行計數(shù)，自動返回采樣數(shù)量。

3 實驗與分析

實驗使用C++與GLSL語言，在Microsoft Windows 7操作系統(tǒng)，OpenGL 4.4.0上實現(xiàn)本文算法，實驗硬件為Inter?Core(TM)i5-3337U處理器，4G內(nèi)存，顯卡為NVIDIA GeForce GT 620M。實驗結果及分析如下：

3.1 算法誤差分析

柵格化處理本身是一個有損處理過程，不可避免地存在誤差[15]，因此必須對柵格化處理后的面積誤差進行分析。相對面積誤差E的定義如下：

%.

(3)

式中:S為分別采用本文3種算法得到的相交面積值；SCPUD為采用基于CPU的確定性算法得到的相交面積值。

實驗數(shù)據(jù)采用隨機產(chǎn)生的兩個簡單多邊形，確定性算法得到的相交面積為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%。本文算法的誤差隨著分辨率的提高呈現(xiàn)下降的趨勢，這樣的舍入誤差在很多工程和一些大型軟件的可接受范圍內(nèi)。

表1 本文算法與確定性算法相交面積誤差統(tǒng)計結果

3.2 算法效率分析

本文3種算法與基于CPU的Monte Carlo算法同為數(shù)值解，因此,在相同的硬件實驗條件下進行效率對比。實驗數(shù)據(jù)采用兩個有“孔洞”的復雜多邊形，實驗結果如表2所示。基于CPU的Monte Carlo算法在表2中簡稱CPUMC,在隨機點數(shù)～分辨率基本相同的數(shù)量級下進行效率對比。

256×256～104時, CPUMC算法的計算時間大約為0.074 6 s，本文3種算法的計算時間分別為0.004 s，0.003 s，0.001 s；數(shù)量級增加3倍(1 024×1 024～106)時，CPUMC為5.147 s，本文3種算法的計算時間分別為0.032 s， 0.024 s，0.002 s。

表2 本文算法與CPUMC算法運行效率統(tǒng)計結果 s

將本文算法的分辨率分別固定在256×256，CPUMC取104個隨機數(shù)和分辨率固定在1 024×1 024，CPUMC取106個隨機數(shù)，比較具有相同形狀但不同頂點數(shù)量的多邊形間的求交計算時間，對比兩種算法的效率差異。如表3和表4所示，其中N為兩個多邊形的頂點數(shù)之和。

表3 本文算法256×256分辨率與CPUMC算法104個隨機數(shù)的運行效率統(tǒng)計結果

表4 本文算法1 024×1 024分辨率與CPUMC算法106個隨機數(shù)的運行效率統(tǒng)計結果

當多邊形所包含的頂點數(shù)量較少時, CPUMC算法具有較高的計算效率,隨著多邊形頂點數(shù)量的增加, CPUMC算法的計算效率快速降低,上萬個點的點包含測試造成了性能的急劇下降，106個隨機數(shù)時，經(jīng)過十多分鐘仍未得到結果。本文算法在處理包含大量頂點的多邊形時更有優(yōu)勢。

4 結束語

本文提出了基于GPU柵格化的多邊形相交面積計算算法。通過實驗，對比和分析了本文算法與基于CPU確定值算法的相對面積誤差率，以及與基于CPU的Monte Carlo算法的執(zhí)行效率。實驗結果表明，本文算法的計算效率不受多邊形頂點數(shù)量的影響，僅對分辨率敏感，隨著分辨率的增加，誤差下降但時間開銷有所增加?；贑PU的確定值算法和Monte Carlo算法在處理小數(shù)據(jù)時較為高效；而本文算法是一種以精度的適當犧牲換取優(yōu)異綜合性能的算法，在處理復雜多邊形及包含大量頂點的多邊形數(shù)據(jù)時具有優(yōu)勢。

[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].計算機研究與發(fā)展，2003,40(4): 576-583.

[6] 劉永奎，高云，黃有群.一個有效的多邊形裁剪算法[J].軟件學報，2003,14(4): 845-856.

[7] 候寶明，劉雪娜.任意多邊形區(qū)域交的有效算法[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-),女,工程師,博士研究生.