基于三角網的等值線填充算法研究

李 強，李 超，甘建紅

1.成都信息工程學院 計算機學院，成都 6102252.成都信息工程學院 發展規劃處，成都 6102253.成都信息工程學院 軟件學院，成都 610225

基于三角網的等值線填充算法研究

李 強1，李 超2，甘建紅3

1.成都信息工程學院 計算機學院，成都 610225
2.成都信息工程學院 發展規劃處，成都 610225
3.成都信息工程學院 軟件學院，成都 610225

1 引言

等值線能夠以圖形的方式表示離散數據，而色斑圖則是由等值線輔以分層著色，以進一步提高等值線圖的直觀效果。由于色斑圖（如溫度，降水色斑圖等）看起來形象直觀，它在GIS、氣象學、地質，石油勘探等許多學科中有著廣泛的應用。色斑圖繪制過程一般分為離散數據三角化及插值[1]，等值線追蹤[2-3]，等值線平滑[4-5]，等值線填充[6-11]及裁剪幾個步驟。

等值線填充就是根據等值線的拓撲關系，對等值線區域進行快速、準確的填充。等值線填充需要解決的兩個關鍵問題是：（1）如何確定兩條相鄰等值線之間的區域；（2）確定這個區域應該用何種顏色填充。在圖形學中有許多不同的算法能對等值線進行填充，如格網法[7]、掃描線填充[8]、Voronoi圖環評[9]等。格網法需要制定基于單個三角形的填充，來拼湊任意形狀的填充，如果多邊形太復雜勢必影響效率。掃描線填充，對填充顏色的確定比較復雜。而Voronoi圖環評用Voronoi圖來構建拓撲二叉樹，雖然算法新穎，但算法不易編程實現。而且這些算法對等值線的特點都沒有進行充分利用，增加了確定填充區域的難度。本文的目的就是從等值線的特點著手，無需判斷等值線首尾點在邊界上的走向[11]，使對等值線拓撲關系及填充顏色的判定達到最簡化，減小傳統算法的時間空間復雜度并結合AS3.0語言的特點提出一種可靠性高，時間復雜度小，等值線拓撲關系及填充顏色判定簡單且適合在web環境中運行的算法。

2 等值線填充算法的基本思路

本文以氣象實時業務系統中基于WEBGIS繪制色斑圖為研究背景，以中國氣象科學數據共享服務網提供的194個氣象基準站的經緯度、溫度以及全國行政邊界的shp文件為數據源。氣象實時業務的汛期（4～9月）數據訪問集中，分析業務集中而頻繁，希望色斑圖繪制分析等業務不負載于服務器，而是使用良好的終端資源，正是基于這樣的需求，結合Flash實現氣象實時業務分析系統，其核心是要實現色斑圖的繪制。本文結合Flash AS3.0語言優勢直接在客戶端讀取shp文件，對氣象站點基于Delaunay快速構建三角網，插值，追蹤，生成等值線，然后對等值線進行平滑及填充，最終生成等值線色斑圖。

3 等值線精度分析

影響等值線的精度的因素有：離散數據網格化，網格點插值以及等值線平滑，而等值線填充不會影響等值線的精度。離散數據網格化可分為矩形網法和三角網法。相對于矩形網格，三角網法的特點[12]有：（1）三角網法不需要將離散點經插值轉換為矩形網格點，減少了計算時間且三角形網格點數據就是離散點數據本身，最大程度上保留了原始數據，其精度顯然高于經插值得到的矩形網格點精度。（2）能通過任意形狀的凸包區域生成等值線，而矩形網法的邊界則只能是矩形。（3）三角網法則對特征點部位任意小的等值線都能繪出，而矩形網格法可能會丟掉很小的閉合等值線。很明顯三角網法在精度上優于矩形網法，所以本文采用已經非常成熟的Delaunay三角形來對氣象站點進行構網，但是由于站點分布不均勻，可能會導致Delaunay三角網大小差別較大，比如產生狹長三角形或面積較小的三角形，最終繪制出來的等值線精度低，可視化效果差。這時可以對初始delaunay網格進行加密[13]，使其成為相對均勻分布的網格，并對新插入的點采用克里金（kriging）或者距離加權反比法（Inverse Distance Weighting，IDW）插值獲得數據值，最后繪制等值線，這樣繪制出來的等值線可靠性更高。注意由于側重點不一樣本文沒有對初始delaunay進行加密，具體算法可以參考文獻[13]，本文主要討論等值線平滑對精度的損失以及等值線追蹤填充算法。

繪制出來的等值線是一條一條的折線，為了使等值線連續平滑，還需要對其進行平滑處理。由于氣象數據是高度離散的，如果偏差一點，可能就會造成很大的精度損失，為了確保空間分析的精度，平滑過后的曲線應該盡量通過原始離散點。常用的等值線平滑算法有二次、三次貝塞爾曲線和二次、三次B樣條曲線等，其中以二次B樣條曲線最為常用。但這些算法平滑過后的曲線都沒有通過原始離散點，對等值線精度造成了一定的損失。三次貝塞爾曲線有兩個固定點（起點和終點），另加兩個決定曲線走向的控制點（CP），本文采用Maxim Shemanarev提出的一種改進的三次貝塞爾曲線（cubic Bezier）算法[14]（以下簡稱MS Bezier算法），使平滑過后的等值線仍然穿過離散點，如圖1所示。其大致思路就是先計算出相鄰原始點的中點，再把相鄰中點連接起來的線段平移到對應的原始點，以平移后的線段所在直線上的點作為控制點，并可以使用一個與控制點和頂點距離相關的系數K，用來沿直線移動控制點，控制點離頂點越近，拐點看起來就越趨于銳利。然后以相鄰原始點為起始點畫三次貝塞爾曲線，這樣就得到了一條穿過所有原始點的光滑曲線，減小了等值線平滑造成的精度損失。其最終實例效果如4.3節圖5所示。

圖1 MS Bezier與二次B樣條平滑效果示意圖

4 三角網構建與等值線追蹤

為了便于闡述，整個算法實現過程的數據結構見圖2。

圖2 算法數據結構說明

4.1 構建三角網及插值

本文采用Delaunay三角形對氣象站點構建三角網。其算法基本思路是由地圖邊界盒先初始化兩個三角形，ΔC1C4C2和ΔC2C4C3如圖3所示。然后不斷向三角網中添加新的點，在添加點的同時，采用Lawson提出的LOP方法對其進行優化，使優化后的三角形滿足空外接圓，最小角最大化的屬性[1-2]。在站點三角化之后，在三角形的邊上進行插值，求得等值點的坐標。在所有三角形邊中找出凸包邊（此處即為C1C2，C2C3，C3C4，C4C1），其判斷標識為該邊只被一個三角形使用。如三角形一條邊為E，則該邊為凸包邊的標識為：E.tris.length==1。并按逆時針（或順時針）依次記錄凸包點（凸包邊端點）及凸包邊上的等值點（插值點），存入數組convexHull中，使其連接形成一個首尾閉合的凸包多邊形，如圖3所示。

圖3 站點三角網及凸包示意圖

4.2 等值線追蹤

等值線追蹤就是在已經插值過的三角形邊上找出值相同的插值點，即等值點，然后連接成等值線。如圖4虛線部分為由氣象站點構建的三角網，從邊AB開始到EF結束的曲線為一條凸包邊封閉等值線（Convex Hull Closed Isoline，CHCL），而圖中較粗的環狀曲線則為一條自封閉等值線（Self Closed Isoline，SCL），為了方便后面等值線的填充，把CHCL和SCL區分開，首先追蹤CHCL，然后再追蹤SCL。只要從凸包邊上某點開始追蹤，那么等值線的尾點必然落在一條凸包邊上，得到的必然是一條CHCL。則遍歷convexHull數組順著凸包邊依次追蹤每一條CHCL，等值線追蹤算法如下：

（1）如圖4從等值點a開始追蹤，首先判斷a點是否已經被等值線穿過（a.lineThrough?=true），已穿過則跳過，否則記錄a點為等值線（a到i，記為a-i）的第一個點，即a-i. coords[0]=a。將a的lineThrough屬性標識為true，并標識currTri（當前三角形）=T1(ΔABC)。

（2）判斷三角形T1(ΔABC)的其他兩條邊AB，AC，即lineThrough=false的兩條邊上是否有與點a值相同的等值點。

（3）通過判斷發現邊BC上有點b與點a值相同，記錄點b為等值線a-i的第二個點，即a-i.coords[1]=b。并將b 的lineThrough屬性標識為true。

（4）判斷點BC是否為凸包邊即BC.tris.length是否等于1，是則終止，否則繼續。此時使用邊BC的有兩個三角形T1，T2，則currTri=[（currTri==BC.tris[0]）?1∶0]，很明顯此時currTri=T2則繼續判斷三角形T2(ΔBCD)的其他兩條邊BD，CD是否有與點b相同的等值點。

（5）重復上面三步，追蹤到凸包邊EF結束（EF.tris. length==1凸包邊一定只被一個三角形使用），這樣按順序記錄了所有等值點的一條CHCL就追蹤完成了。同時記錄下等值線a-i的isolineID為存放等值線的數組的當前索引值，此處為第一條即為0。

（6）然后繼續以逆時針方向開始追蹤凸包邊BE上的下一個等值點j，直到所有CHCL追蹤完畢，這樣就得到一個按凸包邊上等值點依次存放的CHCL數組：chclArr。

（7）最后遍歷非凸包邊，用同樣的算法追蹤SCL，只是SCL追蹤的結束標識為：首尾點為同一個點，即當lineThrough=true時，即可結束追蹤。此時又得到一個存放SCL的數組：sclArr。

圖4 CHCL和SCL等值線示意圖

4.3 等值線繪制及平滑

將追蹤得到的每條等值線繪制出來其實是由逐個點連接起來的一條折線，為了使后面填充的等值區域看起來連續平滑，本文采用前面提到的MS Bezier算法來對等值線進行平滑，其效果如圖5所示。其中底圖為經過MS Bezier平滑的等值線，兩個小方塊中折線為未平滑的等值線，平滑曲線分別為MS Bezier和二次B樣條平滑過后的等值線。

圖5 MS Bezier與二次B樣條平滑效果實例對比

5 等值線填充

等值線填充算法主要解決如何確定兩條等值線之間的區域，以及確定用何種顏色填充的問題。但是由于任意兩條等值線之間區域的形狀是極不規則的，因此單純從確定兩條等值線間區域的思路出發將會使問題變得極為復雜。先來看看等值線的特征與性質：

（1）同一條等值線上的數值一定相等。

（2）相鄰的兩條等值線數值相差一個等值距，只有在凸包邊封閉等值線包含自封閉等值線時有可能相等。

（3）等值線可能自封閉的，也可能是凸包邊封閉的，而不封閉等值線的端點必然落在凸包邊上。

（4）等值線一定不會相互交錯。

（5）對于任意連通區域來說，圍成它的所有等值線值最多為2個。

（6）“大于大的，小于小的”原則：就是在兩條等值線a，b之間有一條子封閉等值線c。若c的數值與a，b中較小者相等，則c內的數值小于c的數值。若c的數值與a，b中較大者相等，則c內的數值大于c的數值。通過對等值線的特征分析，本文提出了一種快速且容易編程實現的算法。

5.1 凸包邊封閉等值線（CHCL）填充

（1）如圖6所示，若以a點為開始點的等值線L為chclArr中的第一條等值線，則L=chclArr[a.isolineID]= chclArr[0]，從L的線頭（L.coords[0]）開始，順著邊界逆時針比較凸包邊上點（凸包點及等值點）的X，Y值是否等于L. coords[n－1]（n=L.coords.length）的X，Y值，直到相等時結束，在convexHull中標識L的線頭和線尾為已追蹤（convexHull[0]. isTraced=convexHull[n].isTraced=true），并標記L的線尾為第一條填充等值線的線尾（convexHull[n].isFirstEndpoint= true），得到了第一個封閉區域。

（2）而該封閉區域的顏色則由逆時針順序的下一個點的值來確定，此處即為＜20，16＞，如果圖例從0開始以4為間距，則該封閉區域應該填充圖例中16～20區間的顏色。

（3）繼續從下一個點開始追蹤，首先判斷點的isTraced屬性，若為true直接跳過，反之則繼續。此時繼續追蹤左邊值為16的等值線，得到第二個封閉區域，此處需要與凸包邊端點b進行比較即＜16，14＞，所以給該封閉區域填充圖例中16～12區間的顏色。

（4）直到第一條填充的等值線L的結束點變成了開始點，即點的isFirstEndpoint屬性為true，此時等值線（L）的下半部分的所有CHCL等值區域都已填充好，則應順著邊界逆時針比較凸包邊上點的X，Y值是否等于L.coords[0]的X，Y值。只對此等值線進行特殊處理，其他等值線都按原來的方法填充。

圖6 CHCL填充區域追蹤示意圖

具體的填充過程，如圖7。

5.2 自封閉等值線（SCL）填充

對于SCL需要根據等值線之間的包含關系以最外層的SCL為根節點構造為一棵棵的拓撲樹，即把所有SCL轉換為數據結構中的森林。首先填充外層SCL，然后依次遞歸填充其子SCL，其填充顏色由其子SCL的值來確定。具體算法實現如下：

圖7 CHCL的填充過程圖

（1）如圖6右下角兩條SCL，12到16之間就應該填充圖例中12～16區間的顏色值。而值為16的SCL內部則應該填充圖例中16～20區間的顏色值。

（2）而對于圖4中左下角的那條SCL，它既沒有父節點又沒有子節點。雖然直觀從等值線圖上以“大于大的，小于小的”原則很容易判定該SCL內部應該填充什么顏色，但是讓程序去判斷包圍該SCL的等值線的拓撲結構是難以辦到的。通過觀察發現SCL內部必然有1個或者多個三角形頂點，可根據其值判斷，如圖4所示。因此判斷標準如下：若E1（EC）和E2（AC）為該SCL穿過的兩條邊，P1，P2分別為其端點，有：

三角形頂點值=E1.P1.Z

三角形頂點值=E1.P2.Z

（3）現在需要解決的問題是如何構建這棵拓撲樹。根據等值線性質4，等值線不能相互交錯，則判斷兩條自封閉等值線的包含關系即可簡化為，判斷某條等值線是否包含另一條等值線的任意一個點，即判斷點是否在多邊形內。判斷點在多邊形內的算法很多，如：水平交叉點數判別法[15]。其算法為O(N)，N為多邊形邊的數目。

（4）建立拓撲樹的偽代碼如下：

假如自封閉等值線SCL的條數為M，由于排除了不符合等值線間距規則的等值線，所以兩次遍歷closedIsos的時間T應該小于等于 M×M，再加上insidePolygon的時間復雜度為 O(n)，則整個算法的時間復雜度小于等于O(M×M×N)，N為多邊形邊的數目，即該等值線插值點個數。以全國194站的溫度為例，填充好的色斑圖如圖8所示。

圖8 全國194站溫度色斑圖

5.3 色斑圖裁剪

等值線圖繪制好后，還需要對超出地區邊界區域的部分進行裁剪。本文使用Actionscript3.0內置了蒙板遮罩的函數，其原理就是，遮罩層決定看到的形狀，被遮罩層決定顯示的內容。那么可以將底圖顯示對象復制為兩份，一份用于畫等值線作為被遮罩層，一份用于決定顯示形狀作為遮罩層。用全國行政邊界的數據和內置函數得到裁剪后的色斑圖如圖9所示。

圖9 溫度色斑圖裁剪前后的對比圖

6 算法效率測試

取全國194站的溫度，用AS3.0語言在CPU為Core 2.1 GHz，內存1 GB的PC上實現了全部算法。表1中列出了本算法及近年來發表的等值線填充算法執行效率比較，衡量標準為站點數比上運行時間（單位：s）。

7 結論

當前，隨著氣象探測系統的進一步發展與完善，往往一個省的觀測站網將達到2 000～3 000個。另一方面，氣象實時分析業務所采用的硬件平臺（PC）比本文實驗中采用的硬件平臺（見表1）性能要好得多，因此，利用該算法在終端處理氣象實時業務的色斑圖分析完全是可行的，從而保證實時業務的穩定性（特別是在每年4～9月汛期）。實踐表明，該算法充分結合AS3.0語言的特點及優勢，易于編程實現，等值線可信度高，適合在web中運行，且運算速度能滿足氣象業務需求。

表1 幾種等值線填充算法效率對比

[1]Bourke P.Efficient triangulation algorithm suitable for terrain modelling[C]//Pan Pacific Computer Conference.Beijing：[s.n.]，1989.

[2]郭新奇，嚴建鋼，楊士鋒，等.基于VC++的等值線追蹤與填充算法[J].兵工自動化，2011，30（4）：81-84.

[3]張順謙.基于Delaunay三角網的區域等值線繪制關鍵算法[J].計算機與網絡，2005，91（25）：39-41.

[4]蔣瑜，杜斌，盧軍，等.基于Delaunay三角網的等值線繪制算法[J].計算機應用研究，2010，27（1）：102-103.

[5]Wen Yihong，Liu Yongjiang.An isoline generating algorithm based on Delaunay[C]//International Conference on Computer Engineering and Technology.Chengdu：[s.n.]，2010：173-175.

[6]黃本宇，王家華，王湘波.復雜地質構造的等值線填充及實現[J].阜陽師范學院，2006，23（4）：51-52.

[7]楊清.C#實現基于三角網的等值線追蹤及填充算法[EB/OL]. （2010-11-29）.http：//blog.sciencenet.cn/u/moustudio.

[8]鄧飛，王美平.基于掃描線轉換的快速等值線填充算法[J].電子技術應用，2006（3）：39-40.

[9]吳培寧，譚建榮.基于Voronoi圖的環評等值線快速拓撲填充[J].浙江大學學報：工學版，2009，43（2）：322-327.

[10]吳自銀，高金耀.一種基于格網的快速等值線充填算法[J].測繪學報，1999，28（4）：351-354.

[11]湯子東，鄭明璽，王思群.一種基于三角網的等值線自動填充算法[J].中國圖象圖形學報，2009，14（12）：2578-2581.

[12]巫昌海.提高三角網法繪制等值線精度的算法[J].武測科技，1993（2）：10-12.

[13]胡金虎.基于不規則三角網的高精度等值線生成方法[J].工程勘測，2011（2）：64-68.

[14]Interpolation with Bezier a very simple method of smoothing polygons[EB/OL].（2006-08-16）.http：//www.antigrain.com/research/bezier_interpolation/index.html.

[15]Bourke P.Determining if a point lies on the interior of a polygon[EB/OL].（2010-12-05）.http：//paulbourke.net/geometry/ insidepoly/.

LI Qiang1,LI Chao2,GAN Jianhong3

1.College of Computer Science and Technology,Chengdu University of Information and Technology,Chengdu 610225,China
2.Development and Planning Department,Chengdu University of Information and Technology,Chengdu 610225,China
3.College of Software,Chengdu University of Information and Technology,Chengdu 610225,China

With deep research on isoline filling algorithm,a new cover algorithm of isoline filling is proposed in this paper.This algorithm classifies isolines as Convex Hull Closed Isoline（CHCL）and Self Closed Isoline（SCL）.Taking full advantage of isoline's characteristic,the judgement of topological relation and filling color is simplified mostly.For CHCL isoline,tracing convex hull isoline-points and convex points are employed to determine filling area,and a multi-branches tree which is established by topological relation is adopted to determine the order of filling for SCL isoline.By repeating cover filling isoline-area orderly,all the isoline-area filling is accomplished.For smoothing isoline,a new method which can decrease the precision loss of isoline smoothing is adopted.Algorithm model is realized using client side programming language-ActionScript3.0,and it shows that this algorithm is easy to implement and the time cost can satisfy the business demands.

Delaunay triangle mesh;isoline tracing;isoline filling;contour plot

通過對現有等值線填充算法的深入研究，提出了一種覆蓋填充等值線的算法。該算法把等值線分類為凸包邊封閉等值線（CHCL）和自封閉等值線（SCL），充分利用等值線的特點，使對等值線拓撲關系及填充顏色的判定達到最簡化。對于CHCL采用對凸包邊等值點及凸包點追蹤來確定填充區域，對于SCL則根據拓撲關系以最外層的SCL為根節點構建一棵多叉樹以確定填充順序.通過對等值區域依次反復覆蓋填充，最終完成所有等值線的填充。采用了一種穿過原始離散點平滑等值線的算法，減小了等值線平滑造成的精度損失。并運用客戶端語言AS3.0（ActionScript3.0）實現了算法模型。實驗結果表明，該算法簡單易于實現，而且運算速度能滿足業務需要。

Delaunay三角網；等值線追蹤；等值線填充；色斑圖

A

TP391.7

10.3778/j.issn.1002-8331.1203-0485

LI Qiang,LI Chao,GAN Jianhong.Study on algorithm of isoline filling based on triangle mesh.Computer Engineering and Applications,2013,49（5）：185-189.

成都信息工程學院人才引進項目（No.KYTZ201040）。

李強（1986—），男，碩士，主要研究方向為WEBGIS空間數據分析及WEB3D；李超（1964—），男，教授，主要研究方向為基于網絡的計算機應用技術；甘建紅（1978—），男，博士，主要研究方向為數字圖像處理。E-mail：lqlunwen@126.com

2012-03-21

2012-07-03

1002-8331（2013）05-0185-05