王結(jié)臣， 蒲英霞， 崔 璨， 陳 剛， 馬勁松

（江蘇省地理信息技術(shù)重點實驗室，江蘇 南京 210093；南京大學(xué)地理信息科學(xué)系，江蘇 南京 210093）

Voronoi圖是計算幾何中的一個經(jīng)典數(shù)學(xué)問題，它在地理信息系統(tǒng)、氣象學(xué)、生態(tài)學(xué)、分子生物學(xué)等學(xué)科有著較成熟的應(yīng)用，自動構(gòu)建Voronoi圖也因此得到較多的研究[1-4]。針對點線面體等多維空間實體的Voronoi圖生成[5-7]都有著廣泛的研究，其中平面點集的Voronoi圖應(yīng)用最為廣泛，其生成算法主要有增量算法、掃描線算法、分治算法等幾種。

增量算法[8-9]是較早建立的一類Voronoi圖生成算法，它通過逐個加入點，并局部修改Voronoi圖結(jié)構(gòu)信息來構(gòu)造新的Voronoi圖。增量算法設(shè)計思路簡單，易于實現(xiàn)，且Voronoi圖的動態(tài)更新較為容易。但其時間復(fù)雜度為 O(n2)，在數(shù)據(jù)量較大的情況下，難以滿足實際應(yīng)用。

掃描線算法[10]是Fortune于上世紀80年代提出，它給定一條垂直掃描線對平面點集從左至右進行掃描，每掃過一個點就為該點構(gòu)造一條以掃描線為基線的拋物線（稱之為點事件），離掃描線最近的一些拋物線在相交和連接后可構(gòu)造出一條岸線（Beach Line）。隨著掃描線的推進，岸線形狀不斷改變。當岸線上某交點到掃描線距離與到該交點對應(yīng)的兩平面點距離相等時，獲得Voronoi圖多邊形的頂點（稱之為圓事件）。按此思路，將掃描線推進到平面點集最右端，Voronoi圖構(gòu)建完成。

分治算法最早由 Shamos和 Hoey提出[11]，該算法預(yù)先將平面點分為若干組，為每個組生成各自的Voronoi圖后，再進行合并。具體實施步驟是將所有點按照x值進行排序，使用一條垂線將點分為兩組，每組點數(shù)大致相等，分別將兩組點生成 Voronoi圖，最后在垂線處實現(xiàn) Voronoi圖的合并。該算法一般通過遞歸方式分組，即在子組內(nèi)構(gòu)造垂直線完成進一步分組。垂線兩側(cè)的圖形合并采用構(gòu)造分割鏈的方法實現(xiàn)，將分割鏈兩側(cè)的完整邊和非完整邊剔除，并把分割鏈補充給新的Voronoi圖，實現(xiàn)合并。Oishi等在此基礎(chǔ)上將Voronoi圖中元素的拓撲方向引入計算，以提升分治算法速度，取得了一定的效果[12]。

一些學(xué)者從其他角度入手研究Voronoi圖的生成及其重要性質(zhì)，如凸包距離函數(shù)法[13]、遞歸算法[14]、離散Voronoi圖算法等[15-16]，這些研究對Voronoi圖生成也取得了較好的效果。

隨著計算機數(shù)據(jù)處理規(guī)模的不斷增大，現(xiàn)有方法越來越難以滿足應(yīng)用需求，針對海量空間數(shù)據(jù)的高效計算方法成為新的研究重點。當代計算機在多核、分布式技術(shù)等方面的發(fā)展，使Voronoi生成算法在并行環(huán)境下的實現(xiàn)成為可能。在將并行計算技術(shù)引入Voronoi圖的生成方面目前也有一些研究成果，如 Cole等在分治算法的基礎(chǔ)上進行并行計算[17]，Ulrith等采用凸包距離函數(shù)實現(xiàn)了一種隨機并行算法[16]。然而，受當時軟硬件條件限制，現(xiàn)有的并行算法多著重于理論分析，算法測試數(shù)據(jù)量較小，一些研究由于受計算環(huán)境的限制難以推廣應(yīng)用。此外，已有的并行算法通常以分治法為基礎(chǔ)，此時分組形狀不易控制，易形成狹長的條帶，在進行分組子網(wǎng)合并時會出現(xiàn)大量邊界點，導(dǎo)致子網(wǎng)間出現(xiàn)過多未封閉多邊形，在一定程度上影響了算法的整體效率。

為此，本文在分析分治算法思想的基礎(chǔ)上，考慮將平面點集進行自適應(yīng)網(wǎng)格劃分，使得每個分組中點的數(shù)量近似相等，且盡量避免每個矩形分組出現(xiàn)狹長條帶。同時，考慮到掃描線算法作為一種較經(jīng)典的Voronoi圖生成算法，其時間復(fù)雜度為 O(nlgn)，算法性能和穩(wěn)定性都比較好。本文采用其作為各點集分組的子網(wǎng)生成算法，并采用并行計算方式完成這一過程。在相鄰網(wǎng)格之間，依據(jù)掃描線算法中的岸線理論提取邊界點，將邊界點集生成Voronoi圖并更新原圖所在的多邊形邊以完成多邊形的合并。

1 平面點集自適應(yīng)分組

1.1 自適應(yīng)分組

為加快平面點集Voronoi圖生成的效率，較有效的策略是對點集進行空間分組，而后進行分組間Voronoi圖的合并。這種“各個擊破”的策略在諸多空間分析算法中有應(yīng)用范例，如分組構(gòu)建凸包、分組生成TIN等。

圖1 基于格網(wǎng)的二維平面點集分組

如圖1所示，點集被分配到6行×9列的空間網(wǎng)格中，若將每個分組單獨生成Voronoi圖，則后續(xù)處理只需完成組與組之間的Voronoi子圖合并即可。這種分組策略與常見的分治算法有一定的相似之處，它也體現(xiàn)了一種分而治之的思想。

規(guī)則格網(wǎng)分組在一定程度上可提高原算法的效率，但當點的空間分布不均勻時，會導(dǎo)致各網(wǎng)格內(nèi)點數(shù)差異太大，影響算法執(zhí)行的整體效率。尤其是進行算法并行化實現(xiàn)時，這種影響尤為顯著，此時由于各組耗時相差太大，導(dǎo)致處理器之間難以均衡而降低計算效能。

圖2 平面點集自適應(yīng)分組

為避免現(xiàn)有分治算法中采用沿掃描線單方向進行分組易形成狹長條帶等問題，同時借鑒規(guī)則網(wǎng)格分組的思想，本文設(shè)計了一種基于二叉樹的點集自適應(yīng)分組方案。如圖2所示，假設(shè)分組后最小點集中點數(shù)不少于20，首先為二叉樹設(shè)定根結(jié)點（圖2中I(1)），該結(jié)點覆蓋原始點集；然后以一條水平分割線將該初始網(wǎng)格單元一分為二，兩個子單元中點數(shù)近似相等；再后分別針對上述每個子網(wǎng)格單元分別采用一條垂向分割線進行劃分，依此類推，直至網(wǎng)格單元中的點數(shù)小于預(yù)設(shè)值2倍時終止分割，得到最終的葉子結(jié)點。劃分過程中采用水平還是垂直分割線，主要取決于網(wǎng)格單元的長寬。該過程很容易使用遞歸方式實現(xiàn)，圖3為圖2的二叉樹分裂示意圖。

圖3 平面點集分組的二叉樹表達

1.2 數(shù)據(jù)結(jié)構(gòu)

筆者稱上述方法為“基于二叉樹的點集自適應(yīng)分組”，在數(shù)據(jù)組織時，按算法中需要表達的空間要素從小到大可分為4個層次：平面點集、Voronoi邊、Voronoi多邊形和二叉樹矩形網(wǎng)格（限于篇幅，本文省略了Fortune掃描線算法中的幾個中間過程結(jié)構(gòu)體），其中二叉樹網(wǎng)格只有在其葉子結(jié)點上才賦值，其他單元賦值為空。這些結(jié)構(gòu)體用C語言可描述為下述形式（見1.2節(jié)）。

//平面點結(jié)構(gòu)體

//點的平面坐標

//該點在網(wǎng)格單元中的編號

//該點在整個平面點集中的編號

//Voronoi邊結(jié)構(gòu)體

//構(gòu)成該邊的直線方程系數(shù)

//Voronoi邊構(gòu)成的線段的兩端點坐標

//Voronoi邊兩側(cè)所對應(yīng)的兩個原始平面點

//標記該邊所在的線段是否與矩形網(wǎng)格單元邊界線相交

//Voronoi多邊形結(jié)構(gòu)體

//構(gòu)成該Voronoi多邊形的邊集合

//Voronoi多邊形對應(yīng)的平面點

//網(wǎng)格單元結(jié)構(gòu)體

//網(wǎng)格單元內(nèi)平面點的個數(shù)

//網(wǎng)格單元內(nèi)平面點集合

//網(wǎng)格單元的邊界最下角點和最上角點//網(wǎng)格內(nèi)Voronoi多邊形邊個數(shù)

//網(wǎng)格內(nèi)Voronoi邊集合

//二叉樹存儲的網(wǎng)格單元集合

//網(wǎng)格單元

//網(wǎng)格單元的左右(上下)兩側(cè)子網(wǎng)格單元

在算法實現(xiàn)時，可以置每個非葉子結(jié)點中的BlockSite指針為空，因為實際的Voronoi子圖計算都是在葉子結(jié)點內(nèi)進行的，根結(jié)點只是為查詢檢索需要，不必實例化其BlockSite指針。

2 邊界點提取與子圖合并

2.1 邊界點提取

在完成自適應(yīng)分組后，每個子集均采用掃描線算法生成子Voronoi圖，如何合并子圖成為該算法的關(guān)鍵環(huán)節(jié)。在實施子圖合并時，位于子圖四周的部分點對應(yīng)的多邊形結(jié)構(gòu)會發(fā)生局部改變，包括多邊形構(gòu)成邊的更新與多邊形閉合，姑且將這類點稱為邊界點，邊界點的界定與提取方法將直接影響算法的穩(wěn)定性。

在生成各組點集的Voronoi圖時，無論網(wǎng)格內(nèi)的點集是按照掃描線從左至右、從右至左、從上至下或從下至上得到，其結(jié)果Voronoi圖是唯一的。設(shè)想將掃描線從無限遠處拉回到網(wǎng)格邊界線上（此時網(wǎng)格邊界線也可看作掃描線），則每個網(wǎng)格對應(yīng)著4條掃描線，這些掃描線處于向外推進的狀態(tài)。

如圖4所示，圖中的鉛垂線為兩相鄰網(wǎng)格的公共邊界，若將此邊界線看作掃描線，向兩側(cè)網(wǎng)格內(nèi)分別作岸線，左側(cè)岸線由點P1、P2、P3產(chǎn)生的拋物線構(gòu)成，右側(cè)岸線由點 Q1、Q2、Q3產(chǎn)生的拋物線構(gòu)成。掃描線兩側(cè)的圖形欲做合并，若固定左側(cè)，將掃描線向右推進，左側(cè)岸線的左后方不受影響；同樣地，若固定右側(cè)，將掃描線向左推進，右側(cè)岸線的右后方也不受影響。因為Voronoi圖中每個頂點正是共享此頂點的多邊形所在平面點的外接圓的圓心，岸線后方所在多邊形頂點及其對應(yīng)的各邊都不會因為Voronoi圖的合并發(fā)生改變，其原始平面點也不必再提取作為邊界點。若上述結(jié)論不成立，例如在圖中掃描線左側(cè)插入新的一點P'，其不為邊界點，但是P3、Q3和P'能夠構(gòu)成空外接圓，而按照掃描線算法理論，滿足這一條件的前提是該點P'必須是組成左側(cè)岸線的點，因為岸線一側(cè)的點與掃描線另一側(cè)的點發(fā)生圓事件時該點必須為構(gòu)成岸線的點，故已經(jīng)排除出構(gòu)造岸線的平面點不能再次發(fā)生圓事件。

圖4 網(wǎng)格邊界掃描線及其兩側(cè)岸線

在圖5中，加粗豎線為兩網(wǎng)格單元的分界線，加粗折線為兩網(wǎng)格單元合并時新加入的 Voronoi多邊形邊，這些新增加的折線段實現(xiàn)了原網(wǎng)格單元邊界處Voronoi多邊形的封閉，每條折線段正是邊界線兩側(cè)某兩點之間的垂直平分線。而網(wǎng)格內(nèi)位于邊界附近的一些Voronoi圖頂點將因子圖合并而舍棄，如圖中的頂點O，它是邊界點A、B和C的外接圓圓心，在進行子圖合并時掃描線右側(cè)的某點與A、B、C構(gòu)建新的3條垂直平分線，這3條平分線將A、B、C間的原平分線進行了分割，從而頂點O不復(fù)存在（它不再是新垂直平分線的交點）。

圖5 網(wǎng)格邊界Voronoi多邊形合并

在算法實現(xiàn)時，針對每條與邊界線相交的邊，將記錄該邊是與網(wǎng)格單元4條邊界線的哪條邊界線相交，根據(jù)該信息可以較方便地提取各邊界線附近的初始邊界點?；诔跏歼吔琰c，按照岸線理論可建立初始岸線，此岸線并不一定與以邊界線為掃描線的最終岸線重合。在網(wǎng)格單元內(nèi)形成的Voronoi圖中，某些平面點生成的多邊形雖然閉合，但是該點與邊界線構(gòu)成拋物線可能與初始岸線仍然有相交情形出現(xiàn)，通常出現(xiàn)在與初始邊界點共Voronoi邊的點中。在檢索其他邊界點時，可采用以下步驟進行：

1）構(gòu)造一個隊列，將初始邊界點逐個加入隊列，成為邊界點隊列，標記每個點已經(jīng)被判定過；

2）逐個遍歷隊列中的邊界點，提取與該邊界點共邊的相鄰點；

3）遍歷這些相鄰點，若該相鄰點已經(jīng)判定過，則跳過；否則為該點構(gòu)造拋物線，計算拋物線與岸線的相交關(guān)系，若相交則更新岸線，并將該點加入隊列，標記該點已被判定過；

4）當隊列中的邊界點都被遍歷過，則終止判定，隊列中的點集就是該條邊界線對應(yīng)的邊界點。

上述步驟的關(guān)鍵點在第3步，當某邊界點相鄰的點都沒有能夠繼續(xù)更新岸線，則終止尋找這些相鄰點的“二階”相鄰點，即停止對該邊界點的外延。因為邊界點所在的未閉合多邊形組合正好是位于邊界點某側(cè)的“第1層”Voronoi多邊形，此時需要以這些多邊形為基礎(chǔ)進行外延，提取后續(xù)滿足能夠更新岸線的新的多邊形。

2.2 Voronoi子圖合并

子網(wǎng)合并是要找到并插入如圖5所示的一條折線，它由兩側(cè)邊界點的垂直平分線所組成，可通過以下方法實現(xiàn)：在為網(wǎng)格單元生成 Voronoi圖的同時，提取每個網(wǎng)格單元的邊界點，將邊界點作為一組點集并生成Voronoi圖；那么，任一網(wǎng)格單元內(nèi)的邊界點將具有兩個 Voronoi多邊形，一個是網(wǎng)格單元自身生成的，一個則是邊界點集生成的。邊界點對應(yīng)的兩個多邊形的交集即是該點合并后的Voronoi多邊形。

盡管邊界點在兩種子圖中的多邊形形狀不盡相同，但是網(wǎng)格內(nèi)部邊界點與邊界點之間的Voronoi邊關(guān)系基本保持不變。因此，同一邊界點對應(yīng)的兩個Voronoi多邊形必定有幾條邊其構(gòu)成的直線方程相同，或者邊的長度位置都未曾發(fā)生變化。事實上，合并后的邊界點對應(yīng)的Voronoi多邊形正是上述兩個多邊形區(qū)域的交集。如圖6所示，圖中的水平線為兩個點集分組的邊界線，兩條連續(xù)的拋物線鏈是該邊界線對應(yīng)的兩側(cè)岸線。兩條岸線之間所夾的一條連續(xù)折線正是在進行合并時需要計算得到的。這條折線在為邊界點集生成的Voronoi圖中可直接得到，如圖6中邊界點 O在分組內(nèi)生成的 Voronoi多邊形包含F(xiàn)ABC，其中FA為射線，AB為線段，BC為射線；而點O在邊界點生成的Voronoi圖中多邊形包含BCDEF。針對這兩個未封閉的多邊形合并方法較為簡單，主要是求取多邊形的交集部分并構(gòu)造新的封閉多邊形，對于已經(jīng)形成的線段可以不用修改，僅僅對射線間進行求交運算即可實現(xiàn)封閉。

圖6 分組間Voronoi多邊形合并

3 分組并行計算

長期以來，并行計算理論在數(shù)值計算及非數(shù)值計算領(lǐng)域有了較快的發(fā)展，先后提出了多種并行計算模型，如PRAM、BSP、LogP等?；谶@些計算模型的常用并行算法主要有兩類：一類是SIMD（單指令流多數(shù)據(jù)流），另一類是MIMD（多指令流多數(shù)據(jù)流）。前一類并行算法中，各處理機同時處理相同的指令，只是所處理的數(shù)據(jù)不同，這主要體現(xiàn)在數(shù)據(jù)處理上的并行性；后一類并行算法中，各處理機處理不同的指令，并可以對不同的數(shù)據(jù)進行獨立操作，處理器之間互不干擾，這是對算法流程自身的一種并行處理。

前人對二維歐式空間平面點集Voronoi圖生成的并行算法也進行了相關(guān)嘗試，如文獻[17]為平面上有n個點的點集提供n個處理器，該算法屬于一種遞歸方法：給定一條垂直分界線每次將平面上的點集分為左右兩組，各組點個數(shù)大約為n/2，并為每個組分配n/2個處理器，然后計算穿越垂直分界線兩側(cè)的“等值線”。算法在每執(zhí)行一次遞歸時就計算分界線兩側(cè)的“等值線”，當遞歸結(jié)束時，合并“等值線”。文獻[18]基于對稱凸距離函數(shù)，提出一種隨機并行算法。該方法采用隨機采樣數(shù)據(jù)構(gòu)建Voronoi圖，在算法中使用固定數(shù)目的處理器?，F(xiàn)有基于并行技術(shù)的Voronoi圖構(gòu)建算法主要分為兩類：一類在構(gòu)建分組后為每個分組分配一個處理器，另一類為固定處理器個數(shù)，分組數(shù)量不定。本文采用后一種方式實現(xiàn)并行運算，在現(xiàn)有多核處理器的基礎(chǔ)上，為平面點集構(gòu)建分組。

本文算法中，采用點集自適應(yīng)分組構(gòu)建并行計算環(huán)境，是一種典型的SIMD算法，每個處理機的工作內(nèi)容相同，只是各自處理不同的分組數(shù)據(jù)。由于算法中還涉及到點集分組、邊界點提取構(gòu)建Voronoi圖、子Voronoi圖合并等環(huán)節(jié)，這3個環(huán)節(jié)是前后相繼處理的。為此，可以在處理機中確定一臺主機，其他處理機作為從機，主機負責點集自適應(yīng)分組、分組數(shù)據(jù)向從機分發(fā)以及接收處理結(jié)果，并在主機上完成子Voronoi圖合并，從機只需對接收到的點集構(gòu)建Voronoi子圖并提取組內(nèi)邊界點。具體處理過程如下：

1）主機按照給定的網(wǎng)格參數(shù)對平面點集進行自適應(yīng)分組，假定分為N組；

2）針對M個處理機，主機逐個為每個處理機分配一個分組點集，并即時監(jiān)聽各從機的運行結(jié)果；

3）主機監(jiān)聽到某從機運算結(jié)束，從機將運算結(jié)果發(fā)送到主機，若主機上還有點集分組未處理完，則繼續(xù)為該從機分配一個點集分組；

4）按上述步驟，直到將主機上的點集分組處理完畢；

5）主機根據(jù)各個點集分組生成的 Voronoi子圖和收集到的邊界點，為邊界點集生成Voronoi子圖；

6）針對上述N+1組Voronoi子圖，完成子圖合并。

4 測試與分析

本文算法首先對二維平面點集進行自適應(yīng)分組，然后在各分組內(nèi)采用掃描線算法構(gòu)建Voronoi圖，最后進行各網(wǎng)格單元子圖的合并。算法耗時主要由點集自適應(yīng)分組、各組點集的Voronoi圖構(gòu)建、邊界點提取與Voronoi圖構(gòu)建、子Voronoi圖合并4部分構(gòu)成。為評估各環(huán)節(jié)的耗時情況，我們進行了一些測試，測試結(jié)果如表1所示。

測試結(jié)果表明：

1）算法耗時主要集中在Voronoi子圖構(gòu)建環(huán)節(jié)，與總耗時相比，點集分組、子圖合并等所占時間甚少；

2）由于 Fortune掃描線算法的時間復(fù)雜度為O(nlogn)，即使不采用并行計算技術(shù)，在單處理機環(huán)境下采用分組計算策略亦能一定程度上提高算法效率；

3）本文算法針對最耗時的環(huán)節(jié)（構(gòu)建Voronoi子圖）進行并行化處理，可顯著提高算法效率。但需補充說明的是，分組數(shù)量并非越多越好，因為它會影響到所提取的邊界點數(shù)量，最極端的情形下（例如每個點作為 1組），原點集中的點均成為邊界點，此時分組策略將失去意義。

表1 算法效率測試結(jié)果表

5 結(jié) 束 語

基于二維歐式空間平面點集的Voronoi圖是Voronoi圖研究領(lǐng)域里面最為成熟也是應(yīng)用最廣的一類Voronoi圖。前人針對這類Voronoi圖的生成給出了數(shù)十種算法實現(xiàn)，主要可以歸納為矢量和柵格方法兩種，柵格方法由于難以保證其數(shù)據(jù)精度，限制了Voronoi圖的應(yīng)用場合，而矢量方法在算法時間復(fù)雜度上難以突破 nlogn，在遇到數(shù)據(jù)量較大的情況下，其運算耗時較長，難以滿足實際應(yīng)用。本文借鑒矢量和柵格方法各自的優(yōu)缺點，提出采用基于二叉樹的自適應(yīng)分組方法，將平面空間點集預(yù)先進行點集分組，并采用并行計算技術(shù)提升算法效率。本算法以 Fortune的掃描線算法為基礎(chǔ)，引入了網(wǎng)格自適應(yīng)分組方法提升算法運算效率，因為分組間的Voronoi圖生成比較獨立，它便于在并行計算環(huán)境下實現(xiàn)，在數(shù)據(jù)量較大時，并行環(huán)境能有效降低單機環(huán)境下算法對內(nèi)存及磁盤空間的限制。

[1] Ryu J,Park R,Kim D S. Molecular surfaces on proteins via beta shapes [J]. Computer Aided Design,2007,39(12): 1042-1057.

[2] Okabe A,Boots B,Sugihara K,et al. Spatial tessellations: concepts and applications of Voronoi diagrams [M]. Wiley,Chichester,2000: 1-42.

[3] Qian Bo,Zhang Lichao,Shi Yusheng,et al. Voronoi approach to recursive generation of tool path for SLS [J].Computer Aided Drafting,Design and Manufacturing,2008,18(2): 13-20.

[4] Li Yang,Feng Wei,Zhou Huamin. Un-thorough delaunay approach for tetrahedral mesh generation [J].Computer Aided Drafting,Design and Manufacturing,2007,17(2): 82-90.

[5] Li Jin,Donguk K,Lisen M,et al. A sweepline algorithm for Euclidean Voronoi diagrams of circles [J]. Computer Aided Design,2006,38(3):260-272.

[6] Yap C K. An O(n log n) algorithm for the Voronoi diagram of a set of simple curve segments [J]. Discrete& Computational Geometry,1987,(2): 365-393.

[7] Kim D S,Cho Y,Kim D. Euclidean Voronoi diagram of 3D balls and its computation via tracing edges [J].Computer Aided Design,2005,37(13): 1412-1424.

[8] Ohya T. Improvements of the incremental method for the Voronoi diagram with computational comparision of various algorithms [J]. Journal of Operation&Research Sociaty,1984,27(4): 306-336.

[9] Held M,Huber S. Topology-oriented incremental computation of Voronoi diagrams of circular arcs and straight line segmengts [J]. Computer Aided Design,2009,41(5): 327-338.

[10] Fortune S. A sweepline algorithm for Voronoi diagrams [J]. Algorithmica,1987,(2): 153-174.

[11] Shamos M I,Hoey D. Closest-point problems[C]//Proceedings of the 16th Annual Symposium on the Fundations of Computer Science,1975: 151-162.

[12] Oishi Y,Sugihara K. Topology-oriented divideand-conquer algorithm for Voronoi diagrams [J].Graphic Models and Image Procedings,1995,57(4):303-314.

[13] Chew L P,Drysedale R L. Voronoi diagrams based on convex distance functions[C]//Proceedings of the Symposium on Computational Geometry,Baltimore,1985: 235-244.

[14] Mark D M. Recursive algorithm for determination of proximal ( thiessen ) polygons in any metric space [J].Geogr. Anal. 1987,19(3): 264-272.

[15] Schwarzkopf O. Parallel computation of discrete Voronoi diagrams[C]// Tech. Rep. Fachbereich Mathematik,Freie Universitat,Berlin,1988: 193-204.

[16] Schueller A. A nearest neighbor sweep circle algorithm for computing discrect voronoi tessellations [J]. Journal of Mathmetical Anaysis and Applications,2007,336(2): 1018-1025.

[17] Cole R,Goodrich M T,et al. A nearly optimal deterministic paranell Voronoi diagram algorithm [J].Algorithmica,1996,16: 569-617.

[18] Ulrich K. A randomized parallel algorithm for Voronoi diagrams based on symmetric convex distance functions [J]. Discrete Applied Mathematics,2001,109(1-2): 177-196.