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

屠文森，汪佳佳

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

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

屠文森，汪佳佳

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

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

Voronoi圖；柵格法；GPU；CUDA

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

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

1 柵格法簡(jiǎn)介

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

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

2 GPU下的柵格生成方法

2.1 CUDA編程模型與GPU

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

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

2.2 并行Voronoi圖柵格生成算法

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

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

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

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

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

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

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

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

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

Step3：GPU分配M×N個(gè)thread對(duì)柵格進(jìn)行處理，當(dāng)M×N大于所有的thread的總數(shù)時(shí)，可以將M×N個(gè)柵格分塊處理，即將其分成A行×B列×C塊，其中A×B小于thread的總數(shù)。對(duì)于分成了C塊的柵格來(lái)說(shuō)，每個(gè)線程只需要處理C個(gè)柵格。

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

Step5：當(dāng)生長(zhǎng)目標(biāo)較多時(shí)，為了減少遍歷生長(zhǎng)目標(biāo)的時(shí)間，通過(guò)借鑒王新生[13]的算法，不計(jì)算柵格點(diǎn)到每一個(gè)生長(zhǎng)目標(biāo)的距離，通過(guò)對(duì)空白柵格不斷的進(jìn)行鄰域擴(kuò)張，直到遇到目標(biāo)生長(zhǎng)點(diǎn)的方法確定此柵格的歸屬。

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

3 實(shí)驗(yàn)與總結(jié)

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

表1 不同方法生成Voronoi圖用時(shí)表ms

圖2 Voronoi圖生成時(shí)間對(duì)比圖

4 結(jié)語(yǔ)

通過(guò)實(shí)驗(yàn)結(jié)果可以看出，采用GPU對(duì)Voronoi圖的生成進(jìn)行并行加速，能夠很好的提高生成速度。其生成Voronoi圖所需時(shí)間與只與生長(zhǎng)目標(biāo)的數(shù)量有關(guān)，而與柵格規(guī)模沒(méi)有關(guān)系，當(dāng)生長(zhǎng)目標(biāo)數(shù)量為n時(shí)，其時(shí)間復(fù)雜度近似于O（n），為線性的生成時(shí)間。相對(duì)于前面的幾種CPU下串行算法，尤其是在柵格規(guī)模過(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動(dòng)態(tài)空間數(shù)據(jù)模型[M].北京：測(cè)繪出版社，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].武漢測(cè)繪科技大學(xué)學(xué)報(bào)，1998，23（3）：208?210.

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

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

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

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

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

[12]王新生，劉紀(jì)遠(yuǎn)，莊大方，等.一種新的構(gòu)建Voronoi圖的柵格方法[J].中國(guó)礦業(yè)大學(xué)學(xué)報(bào)，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—），男，碩士研究生。研究方向?yàn)橛?jì)算機(jī)可視化。

汪佳佳（1990—），女，碩士研究生。研究方向?yàn)橛?jì)算機(jī)可視化。

2014?08?15

國(guó)家重大科學(xué)儀器設(shè)備開(kāi)發(fā)專項(xiàng)（2012YQ05025004）