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

關(guān)鍵詞： Voronoi圖； 柵格法； GPU； CUDA

中圖分類號(hào)： TN919?34 文獻(xiàn)標(biāo)識(shí)碼： A 文章編號(hào)： 1004?373X（2015）04?0066?03

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)用過程中實(shí)現(xiàn)動(dòng)態(tài)更新。間接法則是根據(jù)其對(duì)偶圖Delaunay三角網(wǎng)來構(gòu)造Voronoi圖，因此其性能的高低由所采用的 Delaunay三角網(wǎng)的構(gòu)造算法所決定。增量法通過不斷向已生成的Voronoi圖中增加點(diǎn)來動(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ò)展等問題[8]。因此本文重點(diǎn)研究了Voronoi圖的柵格生成方法，首先比較了常見的柵格方法生成Voronoi圖的優(yōu)缺點(diǎn)，然后結(jié)合CUDA的出現(xiàn)，提出一種基于GPU的Voronoi圖并行柵格生成算法。

1 柵格法簡(jiǎn)介

柵格方法生成Voronoi圖主要是將二值圖像轉(zhuǎn)化為柵格圖像，然后確定各個(gè)空白柵格歸屬。主要方法有兩類，一類以空白柵格為中心，計(jì)算每個(gè)空白柵格到生長(zhǎng)目標(biāo)的距離，以確定其歸屬，常見的方法有代數(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)行覆蓋。通過不斷擴(kuò)大生長(zhǎng)目標(biāo)的半徑，將會(huì)有越來越多的空白柵格被各個(gè)圓所覆蓋，直到最終覆蓋完整個(gè)圖像。數(shù)學(xué)形態(tài)學(xué)距離變換法與圓檢測(cè)法類似，其思想來源于數(shù)學(xué)形態(tài)學(xué)中膨脹操作，膨脹操作起到了擴(kuò)大圖像的效果，通過不斷的對(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語言作為編程語言，提供了大量的高性能計(jì)算指令開發(fā)能力，使開發(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í)行。通過將程序中耗時(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)行的最小單元，由此可見，當(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)為中心通過不斷增長(zhǎng)生長(zhǎng)目標(biāo)半徑對(duì)空白柵格進(jìn)行覆蓋的方法，他們?cè)谟?jì)算每個(gè)空白柵格距離時(shí)，只能通過遍歷柵格，逐一處理。而柵格處理過程中的一個(gè)重要特點(diǎn)是，各個(gè)柵格的計(jì)算并不依賴于其他柵格的計(jì)算結(jié)果。即各個(gè)柵格的計(jì)算是相互獨(dú)立的，而由于CPU的串行性，導(dǎo)致了各個(gè)柵格只能順序處理，降低了處理速度。

由于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è)線程用來確定一個(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ì)偽代碼如下：

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塊的柵格來說，每個(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)過多時(shí)，則轉(zhuǎn)Step5。

Step5：當(dāng)生長(zhǎng)目標(biāo)較多時(shí)，為了減少遍歷生長(zhǎng)目標(biāo)的時(shí)間，通過借鑒王新生[13]的算法，不計(jì)算柵格點(diǎn)到每一個(gè)生長(zhǎng)目標(biāo)的距離，通過對(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í)幾乎不變。

4 結(jié) 語

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

參考文獻(xiàn)

[1] LI Cheng?ming， CHEN Jun. Raster methods of the generation of Voronoi diagrams for spatial entities [J]. International Journal 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 eonstructing 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 fundamental 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] 羅以寧，王開升.散平面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.