矢量多邊形柵格化算法快速并行化方法研究

陳振杰，周琛，李飛雪，李滿春，任沂斌

(南京大學(xué) 江蘇省地理信息技術(shù)重點(diǎn)實(shí)驗(yàn)室，南京 210023)

1 引 言

矢量數(shù)據(jù)和柵格數(shù)據(jù)是地理信息系統(tǒng)(GIS)兩種基本的地理數(shù)據(jù)模型[1]。柵格數(shù)據(jù)更適合進(jìn)行空間分析和空間模擬，在遙感地學(xué)分析、空間決策分析等領(lǐng)域被廣泛應(yīng)用[2-3]。在大區(qū)域地理計(jì)算中，經(jīng)常需要將矢量多邊形數(shù)據(jù)轉(zhuǎn)換為柵格數(shù)據(jù)。目前，已有不少多邊形柵格化方法的研究，形成了內(nèi)部點(diǎn)擴(kuò)散法、復(fù)數(shù)積分法、包含檢驗(yàn)法、掃描線算法和邊界代數(shù)法[4-6]等算法。然而，現(xiàn)有算法大都是串行算法，柵格化效率提升面臨瓶頸，難以滿足海量矢量多邊形快速柵格化的需求[7-8]。

隨著并行計(jì)算技術(shù)的快速發(fā)展及并行硬件的門檻不斷降低，研發(fā)柵格化并行算法成為解決海量數(shù)據(jù)快速柵格化的有效途徑，一些學(xué)者開始研究單個(gè)矢量多邊形柵格化算法的并行實(shí)現(xiàn)。Healey等人設(shè)計(jì)了一種矢量數(shù)據(jù)向柵格數(shù)據(jù)轉(zhuǎn)換的并行算法，并給出矢量柵格化并行算法的選擇標(biāo)準(zhǔn)[9]；Wang設(shè)計(jì)并實(shí)現(xiàn)了一種基于掃描線法的矢量柵格化并行算法[10]。但將逐一改造現(xiàn)有串行算法不僅工作量大，且效率不高。

本文將在分析典型多邊形柵格化算法的基礎(chǔ)上，研究串行算法并行化思路，提出一種多邊形柵格化算法并行框架，以期解決矢量多邊形柵格化算法快速并行化的問題。

2 柵格化算法特征分析

多邊形柵格化一般過程可歸納為：逐個(gè)遍歷所有多邊形，判定每個(gè)多邊形內(nèi)部及邊界上的柵格單元，將多邊形屬性值賦給這些柵格單元。典型的多邊形柵格化算法有內(nèi)部點(diǎn)擴(kuò)散法、復(fù)數(shù)積分法、包含檢驗(yàn)法、掃描線算法和邊界代數(shù)法等。內(nèi)部點(diǎn)擴(kuò)散法通過重復(fù)設(shè)定種子點(diǎn)，填充位于多邊形內(nèi)部及邊界上的種子點(diǎn)柵格，直至多邊形內(nèi)部區(qū)域被填滿(圖1(a))；復(fù)數(shù)積分法與包含檢驗(yàn)法類似，都對(duì)每個(gè)柵格單元，逐個(gè)判定其是否包含在多邊形之內(nèi)，并將多邊形內(nèi)部的柵格單元進(jìn)行填充(圖1(b)、圖1(c))；掃描線算法通過逐行掃描，識(shí)別多邊形內(nèi)部柵格像元條帶，并用多邊形的屬性值將其填充(圖1(d))；邊界代數(shù)法是一種基于積分思想的柵格化方法，通過簡(jiǎn)單的加減代數(shù)運(yùn)算將多邊形屬性值賦給多邊形內(nèi)部及邊界上的柵格單元，實(shí)現(xiàn)多邊形的柵格化(圖1(e))。

上述柵格化算法雖然實(shí)現(xiàn)原理不同，但具有一些共性：①判定多邊形內(nèi)部及邊界上的柵格單元是多邊形柵格化的關(guān)鍵步驟；②多邊形遍歷不依賴于具體算法；③單個(gè)多邊形柵格化的處理都在多邊形最小外包矩形內(nèi)。

圖1 矢量數(shù)據(jù)柵格化算法基本原理示意圖

3 柵格化串行算法并行化方法

3.1 串行算法并行化思路

將串行算法并行化，需要考慮采用何種并行模式、如何將問題域分解。此外，為了能便于將多種柵格化串行算法并行化，需要考慮如何將串行算法的核心功能(或多邊形柵格化基本算子)嵌入進(jìn)來。在分析矢量多邊形數(shù)據(jù)特點(diǎn)和多邊形柵格化算法特征的基礎(chǔ)上，設(shè)計(jì)了一種通用的并行框架。該并行框架主要包括：雙層并行模式、矢量多邊形劃分方法、多邊形柵格化基本算子調(diào)用接口。其中，雙層并行模式是實(shí)現(xiàn)算法并行的基礎(chǔ)和核心，通過充分發(fā)揮進(jìn)程級(jí)并行與線程級(jí)并行的各自優(yōu)勢(shì)，確保盡可能提升算法的并行效率；矢量多邊形劃分方法是問題域劃分的重點(diǎn)，通過對(duì)多邊形數(shù)據(jù)進(jìn)行合理劃分，實(shí)現(xiàn)各進(jìn)程之間的任務(wù)負(fù)載均衡；算法調(diào)用接口確保柵格化串行算法能直接嵌入并行框架中，以實(shí)現(xiàn)串行算法的快速并行化。

3.2 雙層并行模式

在多邊形柵格化過程中，利用進(jìn)程級(jí)并行可實(shí)現(xiàn)算法的全局并行，但在各進(jìn)程內(nèi)部仍是串行處理各個(gè)多邊形。可見，柵格化填充部分在利用進(jìn)程級(jí)并行后仍具有較大并行潛力。線程級(jí)并行由于其共享存儲(chǔ)等特點(diǎn)，可以彌補(bǔ)該不足，從而實(shí)現(xiàn)對(duì)進(jìn)程內(nèi)處理過程的加速[11-12]。

本文將MPI和OpenMP相結(jié)合，設(shè)計(jì)了一種進(jìn)程級(jí)和線程級(jí)雙層并行模式(圖2)。在該并行模式中，MPI進(jìn)程負(fù)責(zé)讀寫數(shù)據(jù)，并對(duì)源矢量多邊形數(shù)據(jù)進(jìn)行劃分；OpenMP線程在各進(jìn)程內(nèi)部發(fā)揮作用，通過循環(huán)獲取各多邊形及其外包矩形區(qū)域的柵格單元，并對(duì)各多邊形調(diào)用多邊形柵格化基本算子，完成對(duì)多邊形區(qū)域的柵格化。具體過程如圖2所示，由2個(gè)計(jì)算節(jié)點(diǎn)中的4個(gè)MPI進(jìn)程分別讀取相應(yīng)的矢量多邊形，進(jìn)而在各進(jìn)程內(nèi)部分別派生出2個(gè)OpenMP線程，并將多邊形分發(fā)給多線程并行處理；每個(gè)OpenMP線程同時(shí)處理1個(gè)多邊形，通過調(diào)用多邊形柵格化基本算子完成柵格化處理；當(dāng)多邊形處理結(jié)束后，由各MPI進(jìn)程將柵格化后的結(jié)果寫入目標(biāo)柵格文件。

圖2 多邊形柵格化雙層并行模式示意圖

3.3 顧及負(fù)載均衡的矢量多邊形劃分方法

針對(duì)柵格化雙層并行模式，本文設(shè)計(jì)了一種顧及負(fù)載均衡的劃分方法，利用多邊形包含頂點(diǎn)數(shù)評(píng)估多邊形計(jì)算復(fù)雜度，并根據(jù)進(jìn)程數(shù)目進(jìn)行均勻劃分(圖3)。具體過程如下：

①多邊形復(fù)雜度排序。復(fù)雜度指標(biāo)的選擇應(yīng)兼顧代表性和計(jì)算簡(jiǎn)便。在柵格化過程中，頂點(diǎn)數(shù)影響柵格化過程中判定與計(jì)算的時(shí)間。因此，本文將頂點(diǎn)數(shù)作為衡量多邊形復(fù)雜度的指標(biāo)，頂點(diǎn)數(shù)越多，則認(rèn)為多邊形復(fù)雜度越高。通過統(tǒng)計(jì)各多邊形的頂點(diǎn)數(shù)，將多邊形按復(fù)雜度由大到小排序，形成待處理多邊形隊(duì)列。

②進(jìn)程間多邊形劃分。根據(jù)多邊形復(fù)雜度、進(jìn)程數(shù)目求得各進(jìn)程處理頂點(diǎn)個(gè)數(shù)Vprocess，順序從待處理多邊形隊(duì)列中取出多邊形歸入各進(jìn)程數(shù)據(jù)隊(duì)列中，使得各進(jìn)程數(shù)據(jù)隊(duì)列中的多邊形頂點(diǎn)數(shù)之和不超過Vprocess，使各進(jìn)程待處理的多邊形復(fù)雜度基本相當(dāng)。其中，Vprocess的計(jì)算方法如下：

Vprocess=∑Vi/nProcess

(1)

其中，i為多邊形編號(hào)，Vi為第i個(gè)多邊形的頂點(diǎn)數(shù)，nProcess為進(jìn)程數(shù)目。

③線程間多邊形動(dòng)態(tài)分配。在各MPI進(jìn)程內(nèi)部，動(dòng)態(tài)地將待處理的多邊形分配各OpenMP線程。各進(jìn)程數(shù)據(jù)每次只分配給各線程1個(gè)待處理多邊形，并標(biāo)識(shí)已分配的多邊形；當(dāng)某個(gè)線程處理完后，再從進(jìn)程數(shù)據(jù)隊(duì)列取出1個(gè)多邊形分配給該線程，直至所有待處理多邊形被處理。這樣，即使進(jìn)程數(shù)據(jù)隊(duì)列中的多邊形復(fù)雜度有較大差異，也可確保線程間的負(fù)載均衡。

圖3 矢量多邊形劃分方法示意圖

3.4 多邊形柵格化基本算子調(diào)用接口

多邊形柵格化基本算子用于實(shí)現(xiàn)單個(gè)多邊形的柵格化，在多邊形柵格化過程中被反復(fù)調(diào)用，是多邊形柵格化的核心和基本功能單元。為確保不同的多邊形柵格化基本算子能直接嵌入到并行框架中，必須規(guī)范多邊形柵格化基本算子的調(diào)用接口。多邊形柵格化基本算子內(nèi)部以串行方式執(zhí)行，通過被并行框架調(diào)用，實(shí)現(xiàn)串行算法的并行化。因此，多邊形柵格化基本算子可直接從現(xiàn)有的多邊形柵格化串行算法中提取。

通過對(duì)常見多邊形柵格化算法參數(shù)的分析，設(shè)計(jì)了多邊形柵格化填充函數(shù)接口參數(shù)。參數(shù)包括：①多邊形頂點(diǎn)信息：將頂點(diǎn)X坐標(biāo)值與Y坐標(biāo)值分開存儲(chǔ)，分別用double*padfX、double*padfY表示。②多邊形屬性值：用double*pValue表示。③目標(biāo)柵格數(shù)據(jù)：用char*pCBData表示，用于臨時(shí)存儲(chǔ)多邊形外包矩形內(nèi)的目標(biāo)柵格數(shù)據(jù)，柵格化后函數(shù)返回更新后的目標(biāo)柵格數(shù)據(jù)pCBData。任何多邊形柵格化填充算法，只要滿足上述接口規(guī)范的要求，且程序開發(fā)語言相同，就可直接嵌入到并行框架中，實(shí)現(xiàn)算法的并行化。

4 測(cè)試結(jié)果與分析

4.1 測(cè)試環(huán)境與測(cè)試數(shù)據(jù)

測(cè)試環(huán)境為浪潮Linux高性能集群。集群包含4個(gè)節(jié)點(diǎn)，每個(gè)節(jié)點(diǎn)包含2顆四核八線程Intel Xeon CPU、12GB內(nèi)存、292GB硬盤、千兆以太網(wǎng)卡。軟件配置為Centos Linux 6.0操作系統(tǒng)、Openmpi 1.4.1、OpenMP 2.0、PostGIS 2.0.3+PostgreSQL 9.0.5數(shù)據(jù)庫。

測(cè)試數(shù)據(jù)為某省土地現(xiàn)狀地類圖斑數(shù)據(jù)。數(shù)據(jù)格式為PostGIS格式，數(shù)據(jù)空間參考系為1980西安坐標(biāo)系，圖斑總數(shù)為1186萬個(gè)，數(shù)據(jù)量為5.5GB。數(shù)據(jù)共有4個(gè)類別，分別以分類1、分類2、分類3和分類4表示。原始試驗(yàn)區(qū)的矢量多邊形如圖4所示。

圖4 矢量多邊形測(cè)試數(shù)據(jù)示意圖

本文選取掃描線算法和邊界代數(shù)法兩種常用多邊形柵格化算法作為算法用例，利用本文形成的并行框架將掃描線與邊界代數(shù)串行算法快速并行化。試驗(yàn)設(shè)定輸出柵格單元大小為20m×20m，輸出的柵格數(shù)據(jù)大小為25880×24225。試驗(yàn)結(jié)果表明，經(jīng)并行化后的掃描線并行算法與邊界代數(shù)并行算法均運(yùn)行穩(wěn)定，且結(jié)果正確。

4.2 并行效率分析

運(yùn)行時(shí)間和加速比是度量算法并行效率的重要指標(biāo)。其中，算法的運(yùn)行時(shí)間是從算法啟動(dòng)直至最后1個(gè)進(jìn)程執(zhí)行完所花費(fèi)的時(shí)間；加速比是同一個(gè)任務(wù)在串行及并行環(huán)境下運(yùn)行時(shí)間的比值[13]。通過測(cè)試本文實(shí)現(xiàn)的不同并行算法的運(yùn)行時(shí)間和加速比，比較并行算法較于串行算法的性能提升和不同算法并行效率的差異。其中，并行算法計(jì)算單元數(shù)目用NP×NT數(shù)表示，NP表示并行算法使用的MPI進(jìn)程個(gè)數(shù)，NT表示OpenMP線程個(gè)數(shù)。對(duì)于同一NP×NT數(shù)，常有不同的進(jìn)程與線程組合，本文測(cè)試采用取得最佳效率的組合來表示。測(cè)試結(jié)果如圖5所示。圖5(a)、圖5(b)為不同并行算法與串行算法的對(duì)比結(jié)果。從中可看出，利用本文的并行化方法改造后的多邊形柵格化算法大大減少了柵格化時(shí)間，取得了良好的加速比，從而驗(yàn)證了本文設(shè)計(jì)的并行化方法的有效性。圖5(c)、圖5(d)為兩種并行算法的效率對(duì)比結(jié)果。結(jié)果表明，邊界代數(shù)法較掃描線算法更為快速地實(shí)現(xiàn)柵格化轉(zhuǎn)換，取得更好的并行效率。原因在于邊界代數(shù)法以整個(gè)柵格單元為基本操作單位，從而實(shí)現(xiàn)對(duì)柵格矩陣的直接操作；掃描線算法則需要進(jìn)行掃描線與多邊形交點(diǎn)的多次判別，因而效率低于邊界代數(shù)法。

4.3 不同并行模式對(duì)比

本文并行算法采用雙層并行模式實(shí)現(xiàn)，為比較其與常規(guī)方法的效率差異，本文實(shí)現(xiàn)了基于MPI并行模式的柵格化算法，并將兩種并行模式進(jìn)行了效率對(duì)比，測(cè)試結(jié)果如圖6所示。從圖中可明顯地看出，對(duì)于每一種算法，雙層并行模式效率均優(yōu)于MPI并行模式。對(duì)于MPI并行模式，當(dāng)進(jìn)程數(shù)超過32時(shí)運(yùn)行時(shí)間緩慢回升；但對(duì)于雙層并行模式，隨著NP×NT數(shù)的增大，運(yùn)行時(shí)間一直降低，直至達(dá)到64。原因在于本文實(shí)驗(yàn)集群共包含32個(gè)使用超線程技術(shù)的計(jì)算核心，對(duì)于MPI并行模式，可用計(jì)算單元數(shù)為32；對(duì)于雙層并行模式，可用計(jì)算單元數(shù)為64。可見本研究所選的雙層并行模式可以更好地利用系統(tǒng)資源，能獲得更好的并行效率。

圖5 并行算法效率測(cè)試結(jié)果

圖6 算法并行模式效率對(duì)比結(jié)果

5 結(jié)束語

本文在分析典型多邊形柵格化算法的基礎(chǔ)上，提出多邊形柵格化算法并行框架，從而形成了一種多邊形柵格化串行算法快速并行化的方法，為解決地理計(jì)算串行算法快速并行化問題提供有

益的借鑒。試驗(yàn)表明，本文設(shè)計(jì)并實(shí)現(xiàn)的并行化方法可實(shí)現(xiàn)柵格化串行算法的快速并行化，利用此方法改造后的算法具有良好的穩(wěn)定性和加速比。在后續(xù)的研究中，將進(jìn)一步探索該并行框架在其他矢量數(shù)據(jù)處理中的適用性，并拓展和深化相關(guān)研究與應(yīng)用。

參考文獻(xiàn)：

[1] MAGUIRE D J，GOODCHILD M F，WRHIND D.Geographical information systems：Principles and applications[M].England：Longman Scientific and Technical，1991.45-54.

[2] GOODCHILD M F.Scale in GIS：An overview[J].Geomorphology，2011，130(1-2)：5-9.

[3] CONGALTON R G.Exploring and evaluating the consequences of vector-to-raster and raster-to-vector conversion[J].Photogrammetric Engineering and Remote Sensing，1997，63(4)：425-434.

[4] 龔健雅.地理信息系統(tǒng)基礎(chǔ)[M].北京：科學(xué)出版社，2001：167-169.

[5] 吳立新，史文中.地理信息系統(tǒng)原理與算法[M].北京：科學(xué)出版社，2003：202-212.

[6] JIMéNEZ J J，F(xiàn)EITO F R，SEGURA R J.A new hierarchical triangle-based point-in-polygon data structure[J].Computers and Geosciences，2009，35(9)：1843-1853.

[7] WANG S，ARMSTRONG M P.A quadtree approach to domain decomposition for spatial interpolation in grid computing environments[J].Parallel Computing，2003，29(10)：1481-1504.

[8] GUAN Q，CLARKE K C.A general-purpose parallel raster processing programming library test application using a geographic cellular automata model[J].International Journal of Geographical Information Science，2010，24(5)：695-722.

[9] HEALEY R，DOWERS S，GITTINGS B，et al.Parallel processing algorithms for GIS[M].London：Taylor and Francis，1997.

[10] WANG Y，CHEN Z，CHENG L，et al.Parallel scanline algorithm for rapid rasterization of vector geographic data[J].Computers and Geosciences，2013，59：31-40.

[11] BOVA S W，BRESHEARS C P，CUICCHI C E，et al.Dual-level parallel analysis of harbor wave response using MPI and OpenMP[J].International Journal of High Performance Computing Applications，2000，14(1)：49-64.

[12] JIN H，JESPERSEN D，MEHROTRA P，et al.High performance computing using MPI and OpenMP on multi-core parallel systems[J].Parallel Computing，2011，37(9)：562-575.

[13] XIE J.Implementation and performance optimization of a parallel contour line generation algorithm[J].Computers and Geosciences，2012，49：21-28.