顧及障礙物的一般圖形Voronoi圖及其加權圖的ArcGIS柵格實現

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

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

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

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

1 基本定義

1.1 一般圖形加權Voronoi圖定義

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

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

1.2 一般圖形障礙加權Voronoi圖定義

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

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

2 顧及障礙物的一般圖形Voronoi圖及其加權圖結晶算法

2.1 算法的基本思想

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

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

2.2 算法的實現步驟

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

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

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

2.3 關于算法的幾點說明

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

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

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

2.3.3 結晶速度確定 算法結晶速度確定如下：

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

根據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 應用實例

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

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

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

圖6 一般圖形障礙加權Voronoi圖Fig.6 Weighted Voronoi diagram with obstacles for general figures

4 結論

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

［1］ 范熙偉.加權Voronoi圖矢量生成算法研究及其實現［D］.西安：西北大學，2011.

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

［3］ 陳軍.Voronoi動態空間數據模型［M］.北京：測繪出版社，2000.

［4］ 張有會，李秀麗，楊立平，等.Voronoi圖畫法的改進與實現［J］.計算機科學，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圖的空間關系計算研究［D］.長沙：中南大學，2002.

［7］ 王新生，劉紀遠，莊大方，等.Voronoi圖用于確定城市經濟影響區域的空間組織［J］.華中師范大學學報（自然科學版），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］ 覃瑜，師學義.利用Voronoi圖的城鄉居民點布局優化研究［J］.測繪科學，2012，37（1）：136－138.

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

［11］ 張有會，淺野哲夫，小保方幸次，等.關于一般圖形Voronoi圖的近似構造法的研究［J］.數值計算與計算機應用，2002（3）：216－225.

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

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

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

［15］ 董雪.障礙Voronoi圖性質及其應用研究［D］.哈爾濱：哈爾濱理工大學，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］ 章毓晉.圖像處理和分析基礎［M］.北京：高等教育出版社，2002.