壽華好， 袁子薇， 繆永偉， 王麗萍

（1. 浙江工業(yè)大學(xué)理學(xué)院，浙江 杭州 310023；2. 浙江工業(yè)大學(xué)計算機科學(xué)與技術(shù)學(xué)院，浙江 杭州 310023；3. 浙江工業(yè)大學(xué)經(jīng)貿(mào)管理學(xué)院，浙江 杭州 310023）

一種平面點集Voronoi圖的細分算法

壽華好1， 袁子薇1， 繆永偉2， 王麗萍3

（1. 浙江工業(yè)大學(xué)理學(xué)院，浙江 杭州 310023；2. 浙江工業(yè)大學(xué)計算機科學(xué)與技術(shù)學(xué)院，浙江 杭州 310023；3. 浙江工業(yè)大學(xué)經(jīng)貿(mào)管理學(xué)院，浙江 杭州 310023）

Voronoi圖是計算幾何中的重要概念之一，在計算機圖形學(xué)、計算幾何、計算機輔助幾何設(shè)計、有限元網(wǎng)格劃分、機器人軌跡控制、模式識別、氣象學(xué)和地質(zhì)學(xué)研究中得到廣泛應(yīng)用。借助于四叉樹和區(qū)間算術(shù)，提出了一種新的構(gòu)造平面點集Voronoi圖的細分算法, 并且和經(jīng)典的增量算法、柵格擴張法進行了比較, 結(jié)果顯示新細分算法更為有效。最重要的是細分算法原理簡單，很容易編程實現(xiàn)。

Voronoi圖；細分算法；增量算法；柵格擴張法；區(qū)間算術(shù)

Voronoi圖的應(yīng)用領(lǐng)域非常廣泛，在Voronoi圖明確定義以前，1644年法國數(shù)學(xué)家 Descartes（笛卡爾）發(fā)表的太陽系及其周邊天體的發(fā)布圖，就是一種 Voronoi圖，而且還是一種加權(quán)Voronoi圖[1]。1932年B.N.Delone給出了“Voronoi區(qū)域”的概念，標(biāo)志著Voronoi圖作為計算幾何的一個研究分支的正式誕生[1]。其后在生物學(xué)、化學(xué)、流體力學(xué)、醫(yī)學(xué)等諸多領(lǐng)域得到廣泛應(yīng)用。作為當(dāng)前計算幾何學(xué)科的一個研究熱點，就其重要性來說，Voronoi圖是僅次于凸殼的一個重要的幾何結(jié)構(gòu)[2]。在不同的領(lǐng)域，Voronoi圖有時也被稱為 Thiessen多邊形[1]、Dirichrit網(wǎng)格、或Wigner-Seitz域等[3]。Voronoi圖的基本定義和算法可見于許多計算幾何教科書。它是關(guān)于空間鄰近關(guān)系一種基礎(chǔ)數(shù)據(jù)結(jié)構(gòu)。Voronoi圖有二維和三維、狹義和廣義、一階和高階之分[2]。其中最基本、應(yīng)用最廣泛和研究最深入的還是二維歐氏空間平面點集 Voronoi圖，平面線集和面集Voronoi圖可以通過平面點集Voronoi圖處理近似獲得[4]。平面點集 Voronoi圖常用構(gòu)造算法主要包括矢量方法[3]和柵格方法[5-6]，基于矢量的方法有對偶生成法、增量構(gòu)造算法、分治法、減量算法、平面掃描算法[3]；基于柵格的方法有鄰域柵格擴張法和柵格鄰近歸屬法[5-6]。在矢量方法中，增量構(gòu)造算法的時間復(fù)雜度在最壞的情況下為O(n2)，其中n為平面點集中點的個數(shù)，在一般情況下其運行時間為，在改進的 Voronoi圖增量構(gòu)造方法中，時間復(fù)雜度為 O( n logn)[7]，分治法、減量算法和平面掃描算法的時間復(fù)雜度為O( n log2n)，間接法的時間復(fù)雜度取決于其對偶Delaunay三角網(wǎng)獲取的時間復(fù)雜度[68]。國內(nèi)外研究者根據(jù)其應(yīng)用對傳統(tǒng)的 Voronoi劃分或Delaunay三角剖分算法進行了一系列的改進[9-11]。柵格法與矢量法相比,思想較為簡單，易于向三維空間擴展，但是它的精度受到柵格單元大小的限制，一般耗時長，精度不高[12-13]，柵格法中相對比較新的方法為一種經(jīng)過改進的柵格擴張法[13]，其時間復(fù)雜度為 O( m2- n2)，其中m為柵格的總數(shù)量。我們在本文中提出的平面點集Voronoi圖的細分算法是一種全新的柵格法，由于我們在細分算法中使用了四叉樹和區(qū)間算術(shù)，當(dāng)判斷出某一平面區(qū)域中沒有Voronoi圖的時候可以整個地把這個平面區(qū)域拋掉，不需要對每一個柵格一個一個地處理，從而極大地提高了柵格法的計算速度，這一點在我們所運行的實例中可以看得很清楚。此外，更重要的是細分算法原理更為簡單，非常容易編程實現(xiàn)。細分算法在最壞的情況下的時間復(fù)雜度為 O (mn)，實際計算的時候由于細分算法可以成片地拋掉不包含 Voronoi圖的平面區(qū)域，從而實際計算量遠遠小于這個數(shù)。

在Voronoi圖中，被用來劃分空間的各個基本圖形元素一般被稱為站點。最基本的 Voronoi圖是以平面點集 p={pi, i =1,2,…,n}為站點的Voronoi圖，它將平面劃分成凸多邊形形狀的Voronoi 區(qū)域，p中的每個站點 pi，對應(yīng)一個這樣區(qū)域 vi，使得 vi內(nèi)的任何點距離 pi比距離其它站點近。本文的基本思想是在給定的矩形區(qū)域中找出所有到站點 pi的最短距離至少在兩處或以上達到的像素（也就是柵格）集合，此像素集合即是平面點集p的Voronoi圖，在求解過程中借助于四叉樹數(shù)據(jù)結(jié)構(gòu)和區(qū)間運算技術(shù)進行加速[14]。

1 細分算法

假設(shè) p1( x1, y1)， p2(x2, y2)… pn(xn,yn)是給定的二維平面上的n個點，所考慮的平面矩形區(qū)域是，像素點大?。聪袼攸c的長度和寬度中較大的那個）為ε。計算該平面矩形區(qū)域內(nèi)這n個點的Voronoi圖，相當(dāng)于計算該平面矩形區(qū)域內(nèi)到這些點的最短距離至少有兩處或兩處以上達到的像素點全體。

根據(jù)以上細分算法的基本原理, 我們給出細分算法的具體步驟如下：

1） 輸入n個平面點 p1( x1, y1)，p2( x2, y2)… pn( xn, yn)以及所在的平面矩形區(qū)域和像素的大小ε。

3） 畫出 p1( x1, y1)， p2(x2,y2)… pn(xn,yn)以及Voronoi圖M，算法結(jié)束。

以上細分算法充分利用了Voronoi區(qū)域的連貫性，當(dāng)算法判斷出某矩形區(qū)域內(nèi)不可能包含Voronoi圖的時候，可以把整個矩形區(qū)域丟掉，從而運行效率比較高。

由于到目前為止，雖然構(gòu)造平面點集Voronoi圖的算法有很多，但是經(jīng)典的增量算法是最常用的構(gòu)造Voronoi圖的方法[7]。柵格法也出現(xiàn)了很多新的算法，本文對比的就是一種新的柵格擴張方法。我們考慮將這里新提出的細分算法與增量算法、柵格擴張法作一個比較，以便通過比較看一看細分算法的表現(xiàn)如何。 經(jīng)典增量算法的基本思想是：在已構(gòu)造的輸入了站點的Voronoi圖的基礎(chǔ)上, 逐步加入新的站點, 再利用局部特性, 通過局部修改已有的Voronoi圖來生成新的Voronoi圖，即對于平面離散的站點集合{p1,p2,… ,pn}，在站點集合{p1,p2,… ,pi}(i ＜ n)的 Voronoi圖的基礎(chǔ)上, 再加入新的站點 pi+1, 然后通過局部修改來構(gòu)造{p1,p2,… ,pi,pi+1}的Voronoi圖。如此不斷加入新的站點, 最后即可得到點集{p1,p2,…,pn}的Voronoi圖[7]。增量算法的具體步驟如下[2]：

1） 產(chǎn)生 p1, p2，做的中垂線，輸出Voronoi圖為中垂線。

2） 產(chǎn)生 p3，連接 p1, p2,p3成三角形，做三邊的中垂線，交點為Voronoi點，從該點引出的三條中垂線構(gòu)成Voronoi圖。

3） 產(chǎn)生 pi，判斷 pi落入那個Voronoi多邊形域內(nèi)，修改該Voronoi多邊形及相應(yīng)Voronoi多邊形的邊與頂點。

4） 直至產(chǎn)生點的工作終止[2]。

其中第 3）步，判斷 pi落入那個 Voronoi多邊形域內(nèi)比較耗費時間。

柵格擴張法的基本思想是：在平面的柵格空間中，兩個柵格 p1( x1, y1)和 p2(x2, y2)間的歐氏距離定義為d(p1.p2)=，計算每一個柵格與其臨近的幾個發(fā)生元柵格之間的歐氏距離，以距離最近的發(fā)生元柵格代碼作為該柵格的隸屬代碼，直至定義區(qū)域中所有柵格單元的歸屬都被檢索完為止，柵格擴張法的具體步驟如下[15]：

1） 構(gòu)建柵格。

2） 生成點集{p1, p2,… ,pn}。

3） 查找每個發(fā)生元的臨近發(fā)生元柵格，即每個點的臨近發(fā)生元柵格。

4） 計算某個柵格與其臨近的發(fā)生元柵格之間的歐氏距離，以距離最近的發(fā)生元柵格的代碼作為該柵格的隸屬代碼。

5） 合并集合。

6） 作圖[15]。

其中第4）步計算距離耗時巨大。其計算機時隨著柵格變小而增加, 但隨著發(fā)生元所占柵格數(shù)量的增加而減小。

2 細分算法的計算復(fù)雜度分析

假設(shè)所考慮的平面區(qū)域中總的柵格數(shù)量為m=k2，而 k=2l，假設(shè)平面區(qū)域中站點的個數(shù)為n。那么在最壞的情況下，也就是說在細分算法的執(zhí)行過程中一個區(qū)域都拋不掉的情況下，總的計算復(fù)雜度為 n +4 n +42n +… +4ln =

l+1n=O( 22ln)=O (k2n)=O (mn)，實際計算的時候由于細分算法可以成片地拋掉不包含Voronoi圖的平面區(qū)域，從而實際計算量遠遠小于這個數(shù)。從這個結(jié)果可以看出，細分算法的計算復(fù)雜度不但跟站點的個數(shù)n有關(guān)，而且跟柵格的總數(shù)量m也有關(guān)。注意到增量算法的時間復(fù)雜度在最壞的情況下為 O( n2)，也就是說增量算法的時間復(fù)雜度只跟站點的個數(shù)n有關(guān)，跟柵格的總數(shù)量m無關(guān)。那么在柵格的總數(shù)量m固定的情況下，隨著站點個數(shù)n的增加，細分算法的計算時間只是線性地增長，而增量算法的計算時間是以平方的速度增長。從這里就可以看出細分算法的優(yōu)勢所在。由于柵格擴張法的時間復(fù)雜度為O( m2- n2)，而站點數(shù) n一般遠遠小于柵格總數(shù)m，所以細分算法的計算速度顯然比柵格擴張法要快，這一點在以下的實例中也可以看得很清楚。

3 實 例

我們用Mathematica8.0編程實現(xiàn)了上面提出的細分算法，經(jīng)典的增量算法和柵格擴張法，并在中央處理器為 Intel(R) Core(TM) CPU i5-2410M @ 2.30 GHz的微機系統(tǒng)里運行程序進行了一些實例的計算和比較，下面給出幾個例子。

例1 輸入兩點 p1(0.24, 0.36), p2(0.52, 0.83),圖1為細分算法的結(jié)果，圖2為增量算法的結(jié)果，圖3為柵格擴張法的結(jié)果。

圖1 2個點的Voronoi圖（細分算法）CPU Time Used:0.156 Second

圖2 2個點的Voronoi圖（增量算法）CPU Time Used:0.374 Second

圖3 2個點的Voronoi圖（柵格擴張法）CPU Time Used:105281 Second

圖4 5個點的Voronoi圖（細分算法）CPU Time Used:1.279 Second

例2 輸入5點：p1={0 .35,0}, p2={0.25, 0.45}, p3={0 .5,1},p4={0 .75,0.5},p5={0 .85,0.1},圖4為細分算法的結(jié)果，圖5為增量算法的結(jié)果，圖6為柵格擴張法的結(jié)果。

圖5 5個點的Voronoi圖（增量算法）CPU Time Used:1.423 Second

圖6 5個點的Voronoi圖（柵格擴張法）CPU Time Used:104441 Second

圖8 10個點的Voronoi圖（增量算法）CPU Time Used:6.115 Second

例3 輸入10點：p1={0 .35,0},p2={0 .25,0.45}, p3={0.5,1},p4={0 .75,0.5},p5={0 .85,0.1}, p6={0 .15,0.3}, p7={0 .95,0.65},p8={0.4,0.2}, p9={0.6,0.7},p10={1 ,0.8},圖 7為增量算法的結(jié)果，圖8為細分算法的結(jié)果，圖9為柵格擴張法的結(jié)果。

圖9 10個點的Voronoi圖（柵格擴張法）CPU Time Used:85742.2 Second

圖10 100個點的Voronoi圖（細分算法）CPU Time Used:144.394 Second

為了說明細分算法的有效性我們還特別增加了一個多點Voronoi圖的例子。

例4 隨機產(chǎn)生的100個點，圖10為細分算法結(jié)果。

4 結(jié) 論

從以上實例可以看出，本文提出的細分算法可以有效的計算出平面點集的Voronoi圖，并且細分算法的運算速度比柵格擴張法要快很多，比經(jīng)典的增量算法也要快，從算法復(fù)雜度分析可以看出增量算法的計算時間隨著點數(shù)的增加以平方的速度增加，而細分算法的計算時間只是線性地增加，從而多點的時候細分算法的優(yōu)勢就更明顯。最重要的是細分算法思想十分簡單，非常容易編程實現(xiàn)，而增量算法和柵格擴張法的原理相對比較復(fù)雜，實現(xiàn)起來要麻煩多了。此外需要指出的是細分算法得到的是一個包含Voronoi圖的像素點的集合，由于區(qū)間算術(shù)的保守性，部分鄰近Voronoi圖但是不應(yīng)包含于Voronoi圖的像素點有可能會因為無法排除而被保留了下來，這使得用細分算法得到的Voronoi圖的圖像在有些地方比實際要粗一些，如何細化是一個值得進一步研究的問題。

[1] 蔡 強. 限定Voronoi網(wǎng)格剖分的理論及應(yīng)用研究[M]. 北京:北京郵電大學(xué)出版社, 2010:5-6.

[2] 周培德. 計算幾何—算法分析與設(shè)計[M]. 北京:清華大學(xué)出版社, 2000:146-216.

[3] 代曉巍, 李樹軍, 劉曉紅. Voronoi圖增點構(gòu)造算法研究[J]. 測繪工程, 2007, 16(1):19-22.

[4] 普雷帕拉塔 F P, 沙莫斯 M I. 計算幾何導(dǎo)論[M].北京:科學(xué)出版社, 1990:226-277.

[5] 王新生, 劉紀(jì)遠, 莊大方, 等. 一種新的構(gòu)建Voronoi圖的柵格方法[J]. 中國礦業(yè)大學(xué)學(xué)報, 2003, 32(3):293-296.

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

[7] 孟 雷, 張俊偉, 王筱婷, 等. 一種改進的Voronoi圖增量構(gòu)造算法[J]. 中國圖像圖形學(xué)報, 2010, 15(6):978-984.

[8] 趙仁亮. 基于Voronoi圖的GIS空間關(guān)系計算[M]. 北京:測繪出版社, 2006:44-45.

[9] 程 丹, 楊 欽, 李吉剛. 二維黎曼流形的Voronoi圖生成算法[J]. 軟件學(xué)報, 2009, 20(9):2407-2416.

[10] Evazi M, Mahani H. Generation of Voronoi grid based on vorticity for coarse-scale modeling of flow in heterogeneous formations [J]. Transport in Porous Media, 2010, 83(3):541-572.

[11] ?alik B. An efficient sweep-line Delaunay triangulation algorithm [J]. Computer Aided Design, 2005, 37(10):1027-1038.

[12] Okabe A, Boots B, Sugihara K, et al. Spatial tessellations:concepts and applications of voronoi diagrams (second edition) [M]. New York:JohnWiley & Sons Ltd, 2000:287-290.

[13] Li Cheng ming, Chen Jun. Raster methods of the generation of Voronoi diagrams for spatial entities [J]. International Journal of Geographical Information Seience, 1999, 13:209-225.

[14] Martin R, Shou H, Voiculescu I, et al. Comparison of interval methods for plotting algebraic curves [J]. Computer Aided Geometric Design, 2002, 19(7):553-587.

[15] 王新生, 劉紀(jì)遠, 莊大方, 等. 一種新的構(gòu)建Voronoi圖的柵格方法[J]. 中國礦業(yè)大學(xué)學(xué)報, 2003, 32(5):293-296.

A Subdivision Algorithm for Voronoi Diagram of Planar Point Set

Shou Huahao1, Yuan Ziwei1, Miao Yongwei2, Wang Liping3
( 1. College of Science, Zhejiang University of Technology, Hangzhou Zhejiang 310023, China; 2. College of Computer Science and Technology, Zhejiang University of Technology, Hangzhou Zhejiang 310023, China; 3. College of Business and Administration, Zhejiang University of Technology, Hangzhou Zhejiang 310023, China )

Voronoi diagram is one of the most important concepts in computational geometry, It is applied widely in computer graphics, computational geometry, computer aided geometric design, finite element grid partition, robot trajectory control, pattern recognition, meteorology and geology. Based on quadtree data structure and interval arithmetic technique, a new subdivision algorithm for Voronoi diagram of a planar point set is proposed. A comparison of this subdivision algorithm with the well known incremental algorithm and grid expansion method is conducted. Test results show that the subdivision algorithm is more efficient. The most important is that the idea of subdivision algorithm is very simple and therefore it is easy to implement.

Voronoi diagram; subdivision algorithm; incremental algorithm; grid expansion method; interval arithmetic

TP 391.72

A

2095-302X (2013)02-0001-06

2012-09-02；定稿日期：2012-01-25

國家自然科學(xué)基金資助項目（61272309，61070135）

壽華好（1964-），男，浙江諸暨人，教授，博士，主要研究方向為計算機輔助幾何設(shè)計與圖形學(xué)。E-mail：shh@zjut.edu.cn