田松，崔希民，孫云華，崔偉宏，李文杰

（1.中國礦業(yè)大學(xué)（北京）地球科學(xué)與測繪工程學(xué)院，北京 100083；2.中國科學(xué)院遙感與數(shù)字地球研究所，北京 100101）

Voronoi圖是一種對空間無縫不重疊的最短距離分割方法［1］，具有勢力范圍、側(cè)向鄰近、局域動態(tài)等特性［2，3］，成功解決了找最近點、求最大空圓、求n個點的凸包、求最小樹等問題［4］，是計算幾何的重要分支，廣泛應(yīng)用于物理、化學(xué)、分子生物、醫(yī)學(xué)、地學(xué)等領(lǐng)域。在地學(xué)領(lǐng)域，Voronoi圖可用于空間插值［5］、空間關(guān)系表達［6］、城市影響范圍分析［7］、最優(yōu)路徑選?。?］、城市規(guī)劃［9］、設(shè)施選址分析［10］等方面。

一般圖形Voronoi圖由普通Voronoi圖生成元擴展為點、線、面而成，是對Voronoi圖理論和應(yīng)用的擴充，當(dāng)前，該領(lǐng)域研究較為成熟。例如，張有會等［11］研究了一般圖形Voronoi圖的近似構(gòu)造法；趙曄等［12］以距離表方式離散生成一般圖形Voronoi圖；曹清潔、安志宏、董雪等［13－15］先后研究了障礙物Voronoi圖的離散生成及應(yīng)用；Gong等［16］利用矢量方式生成了一般圖形加權(quán)Voronoi圖，并實現(xiàn)了點、線、面實體的插入和刪除。關(guān)于一般圖形Voronoi圖與 GIS結(jié)合的應(yīng)用研究，Dong［17］和范熙偉［1］分別以柵格方式和矢量方式結(jié)合GIS開發(fā)組件生成了點、線、面的Voronoi圖及加權(quán)圖，但二者算法均只能單獨處理點、線、面生成元，限制了其應(yīng)用范圍，且未考慮障礙物情況。本文利用ArcEngine，以柵格結(jié)晶方式生成了顧及障礙物的一般圖形Voronoi圖及其加權(quán)圖，該方法不需要考慮生成元形狀，避免了計算Voronoi邊界形狀的大量計算，實現(xiàn)簡單。同時與ArcGIS軟件緊密結(jié)合，為研究和發(fā)展GIS空間數(shù)據(jù)結(jié)構(gòu)和模型提供了重要方法。

1 基本定義

1.1 一般圖形加權(quán)Voronoi圖定義

定義1［11］：設(shè)p為平面上的點，G為平面上的幾何圖形，定義p和G之間的距離為：

定義2：設(shè)Gi（i＝1，2，…，n）為平面上互不相交的幾何圖形，λi（i＝1，2，…，n）是給定的n個正實數(shù)，稱V（Gi）是由生成元Gi確定的權(quán)重為λi的一般圖形加權(quán)Voronoi圖，當(dāng)λ1＝λ2＝…＝λn時，V（Gi）為一般圖形Voronoi圖。

1.2 一般圖形障礙加權(quán)Voronoi圖定義

定義3［18］：設(shè)G為平面上互不相交的幾何圖形，O為平面上的線狀或面狀障礙物，p1、p2為平面上的兩個點。若從任一目標(biāo)點開始，到達另一目標(biāo)點的路徑不與障礙物集合O相交或僅與邊界相交，則稱此路徑經(jīng)過的最短歐氏距離為p1與p2的障礙距離，記為d0（p1，p2）。定義p和G 的障礙距離為：

定義4：設(shè)Gi（i＝1，2，…，n）為平面上互不相交的幾何圖形，λi（i＝1，2，…，n）是給定的n個正實數(shù)，Oj（j＝1，2，3，…，m）為平面上的線狀、面狀障礙物，稱V（Gi）是由生成元Gi確定的權(quán)重為λi的一般圖形障礙加權(quán)Voronoi圖，當(dāng)λ1＝λ2＝…＝λn時，V（Gi）為一般圖形障礙Voronoi圖。

2 顧及障礙物的一般圖形Voronoi圖及其加權(quán)圖結(jié)晶算法

2.1 算法的基本思想

所有被障礙物覆蓋的柵格（記為障礙物柵格）賦予統(tǒng)一顏色值。所有被生成元覆蓋的柵格（記為生成元柵格）指定相應(yīng)顏色，保證不同生成元柵格之間顏色不同，同時與障礙物柵格顏色不同，賦予生成元柵格相應(yīng)權(quán)重值。對于每類生成元柵格，選取邊界柵格作為擴散中心，以指定規(guī)則向外結(jié)晶，并將結(jié)晶區(qū)域賦予同中心柵格相同的顏色。結(jié)晶過程中需檢驗：如果結(jié)晶區(qū)域已著色（包括已賦予其他類別生成元柵格顏色或障礙物柵格顏色），則越過；否則就以中心柵格對應(yīng)的顏色著色結(jié)晶區(qū)域，當(dāng)屏幕所有像素著色完成時，算法結(jié)束。這時，不同顏色區(qū)域間的邊界即為Voronoi圖的邊界近似曲線。

圖1 算法流程Fig.1 Algorithm flow chart

2.2 算法的實現(xiàn)步驟

輸入：生成元集合Gi（i＝1，2，…，n），為平面上互不相交的任意幾何圖形，包含屬性PointID：生成元Gi的唯一標(biāo)識；Color：Gi顏色，Weight：Gi權(quán)重，Gi.Weight＝λi，障礙物集合Oj（j＝1，2，…，m）。

輸出：一般圖形Voronoi柵格圖或加權(quán)柵格圖VR（Gi）（i＝1，2，…，n）；一般圖形 Voronoi矢量圖或加權(quán)矢量圖Vv（Gi）（i＝1，2，…，n）。

算法具體實現(xiàn)步驟如下：1）數(shù)據(jù)獲取。從GIS矢量數(shù)據(jù)集中獲取生成元集合Gi（i＝1，2，…，n），障礙物集合Oj（j＝1，2，…，m），保證任意生成元不重合，任意生成元與障礙物不重合；結(jié)晶速度集合Si（i＝1，2，…，n）；速度閾值thresh。2）生成Gi的最小外包矩形 MBR。3）柵格劃分。設(shè)定柵格大小，將MBR劃分為一系列正方形柵格集合Rk（k＝1，2，…，z），包括屬性GridID：柵格Rk的唯一標(biāo)識；PointID：生成元標(biāo)識，擁有相同PointID的柵格為一類生成元柵格；Color：柵格顏色；CurrentSpeed：記錄算法執(zhí)行過程中的實時結(jié)晶速度；Original－Speed：記錄原始結(jié)晶速度，該值恒定，由生成元權(quán)重確定；State：當(dāng)前柵格狀態(tài)，由Before、Now和Next 3種狀態(tài)組成，分別表示柵格“已參與完操作”、“正參與操作”和“等待參與操作”。4）障礙物區(qū)域賦值。遞歸循環(huán)R，如果Rk∩Oj≠Φ，則令Rk.Corlor＝ζ，ζ為恒定正整數(shù)，表示該柵格為障礙物區(qū)域。5）生成元區(qū)域賦值。遞歸循環(huán)R，如果Rk∩Gi≠Φ，則令Rk.PointID＝Gi.PointID，Rk.Color＝Gi.Color，Rk.OriginalSpeed＝Si，Rk.CurrentSpeed＝Si，Rk.State＝Now。如果Rk∩Gi＝Φ，則令Rk.PointID＝－1，Rk.Color＝－1，Rk.OriginalSpeed＝－1，Rk.CurrentSpeed＝－1，Rk.State＝Next。6）區(qū)域結(jié)晶。遞歸循環(huán)R，如果Rk.State＝Now，同時Rk.PointID＞0，表示Rk是“正參與操作”的生成元柵格，可以參與結(jié)晶過程。如果選擇一般圖形障礙加權(quán)Voronoi結(jié)晶，則計算Rk結(jié)晶速度speed，speed＝Rk.Original－Speed＋Rk.CurrentSpeed，若speed＞thresh，按相應(yīng)結(jié)晶規(guī)則獲取Rk周圍柵格Rq（q＝1，2，…，x），如果Rq.PointID＞0或Rq.Color＝ζ，說明該位置已經(jīng)被填充，繼續(xù)下一次遞歸；否則令Rq.PointID＝Rk.Point－ID，Rq.Color＝ Rk.Color，Rq.OriginalSpeed＝Rk.OriginalSpeed，Rq.CurrentSpeed＝speed－thresh，Rk.State＝Before。如果speed＜thresh，則令Rq.CurrentSpeed＝speed，繼續(xù)下一次遞歸。如果選擇一般圖形障礙Voronoi圖結(jié)晶，則按相應(yīng)結(jié)晶規(guī)則獲取Rk周圍柵格Rq（q＝1，2，…，x），如果柵格Rq.PointID＞0或Rq.Color＝ζ，說明該位置已經(jīng)被填充，繼續(xù)下一次遞歸，否則令Rq.PointID＝Rk.PointID，Rq.Color＝ Rk.Color，Rk.State＝Before。7）狀態(tài)更新。遞歸循環(huán)R，如果Rk.PointID＞0，同時Rk.State＝Next，令Rk.State＝Now。如果全部PointID＞0，算法結(jié)束，否則執(zhí)行步驟6。8）數(shù)據(jù)輸出。輸出文件為一般圖形Voronoi柵格圖或加權(quán)柵格圖VR，根據(jù)實際需要可以將其轉(zhuǎn)化為矢量圖或加權(quán)矢量圖Vv，一并輸出。

2.3 關(guān)于算法的幾點說明

2.3.1 結(jié)晶規(guī)則 Voronoi圖離散生成方式包括距離表法和結(jié)晶法等。由于結(jié)晶法模型眾多，實現(xiàn)靈活，固選作本文算法的主要實現(xiàn)方式，可根據(jù)實際需要，選擇相應(yīng)結(jié)晶方式。常見結(jié)晶規(guī)則如下［13］：4－模板（圖2a），生成元以菱形方式生長；8－模板（圖2b），生成元以正方形方式生長；4－模板、8－模板交替（圖2c），生成元以近似圓形方式生長。3種結(jié)晶規(guī)則分別對應(yīng)數(shù)字圖像中3種不同的距離度量方法［19］：

圖2 圖像結(jié)晶規(guī)則及距離度量方法［13］Fig.2 Image crystallization rule and distance measurement method

2.3.2 顏色的確定及邊界提取 Gi的顏色屬性可以是灰度值或RGB值，相應(yīng)生成VR（Gi）灰度圖或RGB圖。算法執(zhí)行前設(shè)定顏色方式，Color值為隨機生成。在VR（Gi）生成后，通常需要標(biāo)定不同類型生成元柵格間的Voronoi邊，稱為Voronoi邊的近似抽出，其實質(zhì)是柵矢轉(zhuǎn)換的過程。在ArcEngine中，提供相關(guān)接口實現(xiàn)柵矢圖轉(zhuǎn)換。

2.3.3 結(jié)晶速度確定 算法結(jié)晶速度確定如下：

由式（8）和式（9）可知：Si的最大值為1。

根據(jù)2.2步驟（6）給出的公式：Rk.OriginalSpeed＝Si，Rk.CurrentSpeed＝Si，speed＝Rk.OriginalSpeed＋Rk.CurrentSpeed。令thresh＝1，判斷柵格Rk擴張的條件是如果speed＞1，則擴張，Rk.CurrentSpeed＝speed－1；否則Rk.CurrentSpeed＝speed。

3 應(yīng)用實例

利用本文算法，針對不同生成元和障礙物（圖3）可生成不同類型Voronoi圖：點Voronoi圖及加權(quán)圖、線Voronoi圖及加權(quán)圖、面Voronoi圖及加權(quán)圖和一般圖形Voronoi圖及加權(quán)圖等，同時以上情況還可分為顧及障礙物及未顧及障礙物兩種情況。下面給出利用本算法得到的3個實例：1）一般圖形Voronoi圖（圖4）；2）一般圖形加權(quán) Voronoi圖（圖5）；3）顧及障礙物的一般圖形加權(quán)Voronoi圖（圖6）。圖4、圖5的生成元位置、形狀一致，由于其權(quán)重不同，造成Voronoi圖形狀的差異。

圖3 生成元和障礙物Fig.3 Generators and obstacles

圖4 一般圖形Voronoi圖Fig.4 Voronoi diagram for general figures

圖6 一般圖形障礙加權(quán)Voronoi圖Fig.6 Weighted Voronoi diagram with obstacles for general figures

4 結(jié)論

本文借鑒了一般圖形Voronoi圖、加權(quán)圖及障礙物Voronoi圖等研究，利用C＃和ArcEngine，以柵格結(jié)晶方式生成了顧及障礙物的一般圖形Voronoi圖及其加權(quán)圖，取得較好的實驗效果。相比以往研究，該算法有如下優(yōu)點：1）與文獻［1］、［16］等矢量實現(xiàn)方式相比，該算法中不需要考慮生成元的形狀，點、線、面均可作為生成元，避免了矢量方法中計算Voronoi邊界形狀的大量計算。同時，考慮了線狀、面狀障礙物的情況，擴展了算法的應(yīng)用范圍。2）與文獻［12－15］、［17］等柵格實現(xiàn)方式相比，該算法利用GIS開發(fā)組件編程實現(xiàn)，可移植性強，功能靈活，如可選擇不同的結(jié)晶方式（棋盤距離、城區(qū)距離和歐式距離等）、不同的生成元輸出方式（點生成元輸出、線生成元輸出、面生成元輸出或三者的任意組合輸出）、不同的文件輸出格式（灰度圖像、RGB圖像、GIS矢量圖）及是否加權(quán)、是否考慮障礙物等；在輸出GIS矢量文件時，可選擇是否包含生成元標(biāo)識、顏色、權(quán)重等屬性信息。該算法為研究和發(fā)展GIS空間數(shù)據(jù)結(jié)構(gòu)和模型提供了重要方法，擴展了Voronoi圖在地學(xué)領(lǐng)域的應(yīng)用。

［1］ 范熙偉.加權(quán)Voronoi圖矢量生成算法研究及其實現(xiàn)［D］.西安：西北大學(xué)，2011.

［2］ GOLD C M.The meaning of neighbor［J］.Lecture Notes in Computing Science，1992（39）：220－235.

［3］ 陳軍.Voronoi動態(tài)空間數(shù)據(jù)模型［M］.北京：測繪出版社，2000.

［4］ 張有會，李秀麗，楊立平，等.Voronoi圖畫法的改進與實現(xiàn)［J］.計算機科學(xué)，1999，26（11）：86－87.

［5］ LOWELL K，GOLD C.Using a fuzzy surface＿based cartographic representation to decrease digitizing efforts for natural phenomena［J］.Cartography and GIS，1995，22（3）：222－231.

［6］ 趙仁亮.基于Voronoi圖的空間關(guān)系計算研究［D］.長沙：中南大學(xué)，2002.

［7］ 王新生，劉紀(jì)遠，莊大方，等.Voronoi圖用于確定城市經(jīng)濟影響區(qū)域的空間組織［J］.華中師范大學(xué)學(xué)報（自然科學(xué)版），2003，37（2）：256－260.

［8］ SHARIFZADEH M，SHAHABI C.Processing optimal sequenced route queries using Voronoi diagrams［J］.GeoInformatica，2008，12（4）：411－433.

［9］ 覃瑜，師學(xué)義.利用Voronoi圖的城鄉(xiāng)居民點布局優(yōu)化研究［J］.測繪科學(xué)，2012，37（1）：136－138.

［10］ 張偉松.基于Voronoi圖的數(shù)字電視地面廣播站選址分析［D］.北京：中國測繪科學(xué)研究院，2011.

［11］ 張有會，淺野哲夫，小保方幸次，等.關(guān)于一般圖形Voronoi圖的近似構(gòu)造法的研究［J］.數(shù)值計算與計算機應(yīng)用，2002（3）：216－225.

［12］ 趙曄，張有會，趙志輝，等.關(guān)于一般圖形Voronoi圖的離散構(gòu)造法的研究［J］.計算機應(yīng)用與軟件，2004，21（6）：76－78.

［13］ 曹清潔.障礙Voronoi圖的結(jié)晶生成［D］.石家莊：河北師范大學(xué)，2004.

［14］ 安志宏.線段障礙城市Voronoi圖的結(jié)晶生成［D］.石家莊：河北師范大學(xué)，2007.

［15］ 董雪.障礙Voronoi圖性質(zhì)及其應(yīng)用研究［D］.哈爾濱：哈爾濱理工大學(xué)，2011.

［16］ GONG Y X，LI G C，TIAN Y，et al.A vector－based algorithm to generate and update multiplicatively weighted Voronoi diagrams for points，polylines，and polygons［J］.Computer ＆Geoscience，2012，34：118－125.

［17］ DONG P L.Generating and updating multiplicatively weighted Voronoi diagrams for point，line and polygon features in GIS［J］.Computer ＆ Geoscience，2008，34：411－421.

［18］ BERTIN E，CHASSERY J M.3－D Voronoi diagram：Application to segmentation［A］.International Symposium on Voronoi Diagrams in Science and Engineering［C］.1992.197－200.

［19］ 章毓晉.圖像處理和分析基礎(chǔ)［M］.北京：高等教育出版社，2002.