屠文森，汪佳佳

（南京理工大學(xué) 計算機科學(xué)與工程學(xué)院，江蘇 南京210094）

Voronoi圖柵格生成算法GPU并行實現(xiàn)

屠文森，汪佳佳

（南京理工大學(xué) 計算機科學(xué)與工程學(xué)院，江蘇 南京210094）

針對矢量法生成Voronoi圖計算與存儲復(fù)雜的缺點，重點分析研究了Voronoi圖的柵格生成方法。對不同的柵格生成算法的復(fù)雜性和效率進行了比較分析，并針對以往方法速度較慢的問題，提出一種CUDA平臺下GPU并行柵格掃描的方法。該方法利用GPU的多線程特性，將各個柵格的計算分散到不同的線程中并行處理。相比其他柵格生成方法，該方法不需要考慮柵格的規(guī)模，能夠以幾乎線性的時間完成Voronoi圖的生成，極大地提高了生成速度。

Voronoi圖；柵格法；GPU；CUDA

Voronoi圖是一種空間分割算法。其是對空間中的n個離散點而言的，它將平面分割為n個區(qū)域，每個區(qū)域包括一個點，此區(qū)域是到該點距離最近的點的集合。由于Voronoi圖具有最鄰近性，鄰接性等眾多性質(zhì)和完善的理論體系，其被廣泛的應(yīng)用在地理學(xué)、氣象學(xué)、結(jié)晶學(xué)、航天、機器人等領(lǐng)域。

Voronoi圖的生成主要有矢量方法和柵格方法[1?3]。矢量法中，典型的方法有增量法、分治法和間接法[4?6]。分治法是一種遞歸方法，算法思路簡單，但是很難在應(yīng)用過程中實現(xiàn)動態(tài)更新。間接法則是根據(jù)其對偶圖Delaunay三角網(wǎng)來構(gòu)造Voronoi圖，因此其性能的高低由所采用的Delaunay三角網(wǎng)的構(gòu)造算法所決定。增量法通過不斷向已生成的Voronoi圖中增加點來動態(tài)構(gòu)建Voronoi圖。相對于前兩種方法，增量法構(gòu)造簡單并且容易實現(xiàn)動態(tài)化，所以被廣泛應(yīng)用[7]。矢量方法的優(yōu)勢是生成Voronoi圖精度高，但是存在存儲復(fù)雜，生長元只能是點和線，以及難以向三維及高維空間擴展等問題[8]。因此本文重點研究了Voronoi圖的柵格生成方法，首先比較了常見的柵格方法生成Voronoi圖的優(yōu)缺點，然后結(jié)合CUDA的出現(xiàn)，提出一種基于GPU的Voronoi圖并行柵格生成算法。

1 柵格法簡介

柵格方法生成Voronoi圖主要是將二值圖像轉(zhuǎn)化為柵格圖像，然后確定各個空白柵格歸屬。主要方法有兩類，一類以空白柵格為中心，計算每個空白柵格到生長目標的距離，以確定其歸屬，常見的方法有代數(shù)距離變換法，逐個空白柵格確定法等；另一類以生長目標為中心，不斷擴張生長目標的距離半徑，填充其中的空白柵格，直到將整個圖像填充完成，主要有圓擴張法，數(shù)學(xué)形態(tài)學(xué)距離變換法等。代數(shù)距離變換法對距離圖像進行上行掃描（從上到下，從左到右）和下行掃描（從下向上，從右到左）兩次掃描，計算出每個空白柵格最鄰近的生長目標，以此生長目標作為其歸屬。此方法中柵格距離的定義直接影響了空白柵格的歸屬和Voronoi圖的生成精度，通常使用的柵格距離定義有街區(qū)距離、八角形距離、棋盤距離等[8]。距離變換的柵格生成方法精度

低、耗時長[2，7]，所需要花費的時間和柵格的數(shù)量成正比，當柵格為n×n大小時，其時間復(fù)雜度為O（n×n）。圓檢測法[9]以生長目標為圓心，以一定的步長為初始半徑，所有生長目標同時對其構(gòu)成的圓內(nèi)的空白柵格進行覆蓋。通過不斷擴大生長目標的半徑，將會有越來越多的空白柵格被各個圓所覆蓋，直到最終覆蓋完整個圖像。數(shù)學(xué)形態(tài)學(xué)距離變換法與圓檢測法類似，其思想來源于數(shù)學(xué)形態(tài)學(xué)中膨脹操作，膨脹操作起到了擴大圖像的效果，通過不斷的對生長目標進行膨脹操作，最終擴張到所有的空白柵格。這兩種方法有個共同的缺點，在每次擴張后，都需要判斷整個柵格圖像是否已完成擴張，而這需要遍歷柵格圖像，十分耗時。

2 GPU下的柵格生成方法

2.1 CUDA編程模型與GPU

CUDA是一個并行編程模型和一個軟件編程環(huán)境，其采用了C語言作為編程語言，提供了大量的高性能計算指令開發(fā)能力，使開發(fā)者能夠在GPU的強大計算能力上建立起一種更加高效的密集數(shù)據(jù)計算解決方案[10]。

CUDA將CPU作為主機端，GPU作為設(shè)備端，一個主機端可以有多個設(shè)備端。其采用CPU和GPU協(xié)同工作的方式，CPU主要負責(zé)程序中的串行計算的部分，GPU主要負責(zé)程序中的并行計算的部分。GPU上運行的代碼被稱為內(nèi)核函數(shù)，其能夠被GPU上內(nèi)置的多個線程并行執(zhí)行。一個完整的任務(wù)處理程序由CPU端串行處理代碼和GPU端并行內(nèi)核函數(shù)共同構(gòu)成。當CPU中執(zhí)行到GPU代碼時，其首先將相關(guān)數(shù)據(jù)復(fù)制到GPU中，然后調(diào)用GPU的內(nèi)核函數(shù)，GPU中多個線程并行執(zhí)行此內(nèi)核函數(shù)，當完成計算后，GPU端再把計算的結(jié)果返回給CPU，程序繼續(xù)執(zhí)行。通過將程序中耗時的且便于并行處理的計算轉(zhuǎn)移到GPU中使用GPU并行處理，以提高整個程序的運行速度。CUDA是以線程網(wǎng)格（Grid），線程塊（Block），線程（Thread）為三層的組織架構(gòu)[11]，每一個網(wǎng)格由多個線程塊構(gòu)成，而一個線程塊又由多個線程構(gòu)成，如圖1所示。在GPU中，線程是并行運行的最小單元，由此可見，當存在大量的線程時，程序的并行程度將會十分高。目前的GPU上一個網(wǎng)格最多包含65 535×65 535個線程塊，而一個線程塊通常有512個或1 024個線程，所以理論上可以對65 535× 65 535×512個柵格同時進行計算。

2.2 并行Voronoi圖柵格生成算法

傳統(tǒng)的柵格生成算法中，不論是采用以空白柵格為中心確定其歸屬的方法，還是以生長目標為中心通過不斷增長生長目標半徑對空白柵格進行覆蓋的方法，他們在計算每個空白柵格距離時，只能通過遍歷柵格，逐一處理。而柵格處理過程中的一個重要特點是，各個柵格的計算并不依賴于其他柵格的計算結(jié)果。即各個柵格的計算是相互獨立的，而由于CPU的串行性，導(dǎo)致了各個柵格只能順序處理，降低了處理速度。

圖1 GPU組織架構(gòu)

由于GPU下的多個線程都是硬件實現(xiàn)的，各個線程的處理都是并行的，因此將柵格距離的計算分散到GPU端各個線程，必然能夠提高其生成速度。為了并行處理柵格化圖像，可以采用如下的想法，將每一個柵格點對應(yīng)于一個線程，此線程計算此柵格到所有的生長目標的距離，取最小距離的生長目標作為其歸屬。即采用一個線程用來確定一個空白柵格歸屬的方法。

確定方法后，就需要對GPU端內(nèi)核函數(shù)進行設(shè)計，由于內(nèi)核函數(shù)是并行處理的執(zhí)行單元，其設(shè)計方式直接決定了GPU端的程序運行效率。因此如何設(shè)計良好的內(nèi)核函數(shù)是提高并行速度的關(guān)鍵。本文采用如下方式進行內(nèi)核函數(shù)的設(shè)計，假設(shè)共分配了K個并行處理線程，柵格規(guī)模為M×N，設(shè)A[i]為第i個線程處理的柵格編號。當K

由于顯卡上的內(nèi)存是動態(tài)隨機存儲（DRAM），因此最有效率的存取方式，是以連續(xù)的方式存取。當采用第一種方式時，看似是一種連續(xù)的存取方式，實際上恰好是非連續(xù)的，當?shù)趇個線程處理第i個柵格時，由于處理需要一定的時間，此時GPU自動將下個一線程i+1需要的內(nèi)存數(shù)據(jù)取出給其使用，此時下一個線程的內(nèi)存數(shù)據(jù)

卻是在i+C處，內(nèi)存變成了間斷存取。而在使用第二種方式進行處理時，恰好是一種連續(xù)的存取方式，由于第i個線程正在處理第i個柵格數(shù)據(jù)，此時GPU為第i+1個線程準備數(shù)據(jù)，而此時的數(shù)據(jù)正好為第i+1內(nèi)存處。滿足了內(nèi)存的連續(xù)存取特性。因此本文采用第二種方式，內(nèi)核函數(shù)的設(shè)計偽代碼如下：

具體步驟如下：（這里假設(shè)柵格的規(guī)模為M×N）：

Step1：根據(jù)柵格圖像的規(guī)模，確定GPU端線程塊和線程的分配方式和分配數(shù)量，初始化GPU端的參數(shù)。

Step2：程序調(diào)用GPU端內(nèi)核函數(shù)，同時將待處理柵格圖像數(shù)據(jù)傳入GPU中。數(shù)據(jù)主要是圖像的柵格距離，一般是二維數(shù)組，0表示空白柵格，其他各生長目標可由1，2等不同的數(shù)字定義。

Step3：GPU分配M×N個thread對柵格進行處理，當M×N大于所有的thread的總數(shù)時，可以將M×N個柵格分塊處理，即將其分成A行×B列×C塊，其中A×B小于thread的總數(shù)。對于分成了C塊的柵格來說，每個線程只需要處理C個柵格。

Step4：當生長目標數(shù)目不多時，每一個線程計算其對應(yīng)的柵格到所有的生長目標點的距離，取距離最小的生長目標，為此線程對應(yīng)的空白柵格的歸屬，轉(zhuǎn)Step6。當生長目標過多時，則轉(zhuǎn)Step5。

Step5：當生長目標較多時，為了減少遍歷生長目標的時間，通過借鑒王新生[13]的算法，不計算柵格點到每一個生長目標的距離，通過對空白柵格不斷的進行鄰域擴張，直到遇到目標生長點的方法確定此柵格的歸屬。

Step6將生成后的數(shù)據(jù)返回CPU端，CPU端完成柵格圖像的顯示與后處理。

3 實驗與總結(jié)

在CPU參數(shù)為Intel?Xeon?CPU E5?2609，2.4 GHz，2處理器8核心，GPU參數(shù)為TeslaC2075，448CUDA核心，顯存5.25 GB的試驗平臺下，做了不同方法在不同柵格規(guī)模下生成Voronoi圖的對比試驗，試驗中生長目標的個數(shù)定義為100個。由于不同的方法都采用了相同的距離定義，因此各種方法的Voronoi圖生成結(jié)果都是相同的，即他們之間的生成精度是相同的，所以這里重點比較了不同方法的生成耗時。表1列出了不同方法生成Voronoi圖的用時，圖2為表1的折線圖，從圖2中可以明顯看出，當柵格數(shù)量較少時，GPU并行技術(shù)的使用并不能提升生成速度，但是當柵格點數(shù)量增加時，逐點法和距離變換法用時明顯增加，但GPU并行算法用時幾乎不變。

表1 不同方法生成Voronoi圖用時表ms

圖2 Voronoi圖生成時間對比圖

4 結(jié)語

通過實驗結(jié)果可以看出，采用GPU對Voronoi圖的生成進行并行加速，能夠很好的提高生成速度。其生成Voronoi圖所需時間與只與生長目標的數(shù)量有關(guān)，而與柵格規(guī)模沒有關(guān)系，當生長目標數(shù)量為n時，其時間復(fù)雜度近似于O（n），為線性的生成時間。相對于前面的幾種CPU下串行算法，尤其是在柵格規(guī)模過大的情況下，能夠很好的提高Voronoi圖的生成速度。

[1]LI Cheng?ming，CHEN Jun.Raster methods of the generation of Voronoi diagrams for spatial entities[J].International Jour?nal of Geographical Information Science，1999，13：209?225.

[2]陳軍.Voronoi動態(tài)空間數(shù)據(jù)模型[M].北京：測繪出版社，2002.

[3]KOKICHI Sugihara.VORONOI2：a fortran program for eon?structing the Voronoi Diagram[J].Geographic Systems，1994（1）：347?349.[4]OKABE A，BOOTS B，SUGIHARA K，et al.Spatial tessel?lations：concepts and applications of Voronoi diagrams[M]. 2nd ed.New York：John Wiley and Sons，2000.

[5]HAMMER Auren.Voronoi diagrams：a survey of a fundamen?tal geometric data Strueture[J].ACM Computing Surveys，1991（23）：345?405.

[6]李成名，陳軍.Voronoi圖生成的柵格算法[J].武漢測繪科技大學(xué)學(xué)報，1998，23（3）：208?210.

[7]李成名.基于Voronoi圖的空間關(guān)系描述、表達與推斷[D].武漢：武漢測繪科技大學(xué)，1998.

[8]陳軍，趙仁亮.Voronoi動態(tài)空間數(shù)據(jù)模型[M].北京：測繪出版社，2002.

[9]羅以寧，王開升.散平面Voronoi圖的光柵圖形算法[J].四川大學(xué)學(xué)報，2003（3）：596?599.

[10]李波，趙華成，張敏芳.CUDA高性能計算并行編程[J].微型電腦應(yīng)用，2009，25（9）：55?57.

[11]吳焰斌.CUDA編程模型[J].科技風(fēng)，2009（3）：63?64.

[12]王新生，劉紀遠，莊大方，等.一種新的構(gòu)建Voronoi圖的柵格方法[J].中國礦業(yè)大學(xué)學(xué)報，2003（3）：293?296.

Raster?based method for Voronoi diagram using GPU parallel technology

TU Wen?sen，WANG Jia?jia
（School of Computer Science and Engineering，Nanjing University of Science and Technology，Nanjing 210094，China）

Aimed at the complexity of calculation and storage in Vector?based method for Voronoi diagram，the raster?based method is researched emphatically.different methods’complexity and efficiency of generating the Voronoi diagram are an?alyzed.A raster?based method for Voronoi diagram generating with GPU parallel technology is raised to resolve the problem of low speed.Compared with other methods，grid size was not took into account in this method.It improves the generation speed ob?viously

Voronoi diagram；Raster?based method；GPU；CUDA

TN919?34

A

1004?373X（2015）04?0066?03

屠文森（1989—），男，碩士研究生。研究方向為計算機可視化。

汪佳佳（1990—），女，碩士研究生。研究方向為計算機可視化。

2014?08?15

國家重大科學(xué)儀器設(shè)備開發(fā)專項（2012YQ05025004）