一種Delaunay三角剖分的改進算法

余代俊，蒲朝旭，朱逍賢

(1. 成都理工大學 現代工程測量技術及應用研究所，四川 成都 610059; 2. 四川科技職業學院 土木與建筑工程學院，四川 成都 611745)

一、引 言

不規則三角網(TIN)是二維平面空間對任意離散點實行GIS數據表達、管理、可視化，以及地學分析、DEM分析、計算機視覺等方面的一項不可或缺的應用技術與手段[1]。Delaunay三角網即D-TIN，它能夠很好地滿足TIN的3點基本要求，即唯一性、最大最小角特性、空圓特性，能夠對給定區域點集進行最佳的三角剖分[2]。

正因為TIN模型具有眾多優點，同時其數據結構比較簡單，因此它在許多領域得到了廣泛應用。但隨著計算機技術和地理信息技術的發展，如何快速高效地處理大量數據和對數據進行高精度模擬與仿真就成了該方面研究的主要內容[3]。

本文就如何改進和優化Delaunay三角網的構建算法來滿足海量數據的處理需求，以及如何將構建算法利用較為簡潔的程序進行實現，進行了深入研究和詳細探討。

二、算法設計

筆者結合已有三角網的生成算法[4-12]，取長補短，提出了基于凸包實現的Delaunay三角網剖分算法。

1. 數據存儲結構

一個好的數據結構直接影響著D-TIN的生成效率，D-TIN的數據結構主要包括點、線、面之間的組織關系和拓撲關系。

本文所使用的數據結構如下

struct Point3d∥點結構

{

double X;∥X坐標

double Y;∥Y坐標

double Z;∥Z坐標

}

struct Edge∥邊結構

{

Point3d StartPt;∥起點

Point3d EndPt;∥終點

Triangle LeftTri;∥左三角形

Triangle RightTri;∥右三角形

}

struct Triangle∥三角形結構

{

Edge edge1;∥邊1

Edge edge2;∥邊2

Edge edge3;∥邊3

}

基于上述結構的設計，三角形的各個頂點和邊均采用逆時針的方式進行存儲，在程序實現過程中采用鏈表和字典(字典是采用鍵-值的存儲結構，能夠通過鍵直接取出值)等數據存儲結構，能夠快速實現數據的存取及快速定位。

2. 算法實現

該算法主要分3步完成：首先，對輸入數據進行預處理，生成凸包；其次，在凸包的基礎之上采用逐點插入法生成三角網；最后對三角網采用LOP優化算法生成最終的三角網。

(1) 生成凸包

凸包的生成是基于Akl-Toussaint啟發式函數[13]和斜率判別法來建立的，如圖1所示。

1) 構建初始凸包：將輸入點集中maxx、maxy、minx和miny的4個點或是與該點集最小外圍邊框的西南角和東北角最近的4個點作為初始凸包，如圖1中的P11、P14、P6和P2點組成的虛多邊形，并且按照逆時針的方式將此4點存入凸包鏈表中，如果存在多個minx最小的點，則取其中y坐標最小的點，其余類似情況依此處理。

2) 清理點：建立初始凸包之后，該初始凸包的4個點肯定位于最終凸包點集中，而該4點所組成的四邊形中所包含的點則肯定不會出現在最終凸包點集中，為此使用Akl-Toussaint啟發式函數將位于初始凸包中的所有點全部刪除，將剩余的點用于凸包的構建計算。該Akl-Toussaint啟發式函數的性能為O(n)，對于大量的隨機點而言，該法能夠排除掉幾乎一半的點，能夠明顯提高后續建立凸包的性能。

圖1中的初始凸包由P11、P14、P6和P2組成，將該4點組成的四邊形內所包含的點全部剔除，即剔除6個點(P8、P9、P5、P4、P10、P13)，只留下剩余的4個點(P1、P3、P7、P12)來判斷是否為凸包點，明顯提高了程序的運行效率。

3) 排序點：對進行清理之后的點集按照X坐標升序為主、Y坐標升序為輔的方式進行排序，以便用于后續凸包的生成。對于點集的排序可以采用LINQ方式進行，該方法是所有排序算法中最高效的。

4) 修改凸包：取出初始凸包鏈表中相鄰的兩個點，檢查位于其右側斜率最小的點，該點亦即是其右側距離該兩點連線距離最大的點(可以采用三角形面積法判斷)。若該點存在，則將該點插入到上述兩點之間；否則不插入，進行下一條邊的判斷。循環執行上述過程，直至所有邊均沒有插入點則完成凸包的搜索。在進行點搜索的時候，每添加一個點，也可以添加Akl-Toussaint啟發式函數的判斷，將位于新加入點與該邊組成的三角形中的點剔除掉。

5) 生成凸包：將修改之后的凸包代替初始凸包則完成了點集凸包的搜索。

在構建凸包過程中，其時間復雜度主要是修改凸包所花費的時間，其時間復雜度和Akl-Toussaint啟發式函數的性能一致，均為O(n)。

(2) 建立三角網

1) 預處理：將所有凸包點從原始點集中剔除，所有凸包點已經是點集的最外圍點，不再用于插入構網。

2) 構建初始三角網：從剔除凸包點之后的點集中任意取出一點P9，如圖2所示。將該插入點與凸包的所有頂點連接，將生成的三角形存入三角形鏈表中，將各邊存入邊鏈表中。

圖2 Delaunay三角網建立過程示意圖

3) 修改三角網：將剩余的點依次插入到該初始三角網中，每插入一個點需要判斷該點位于哪一個三角形內，同時對該三角形進行重構，而不用重構整個三角網，能夠大大提高構網效率。判斷點位于哪一個三角形時需要遍歷三角形鏈表，在遍歷比較時可以進行排斥試驗，如果該點在該三角形的最小外圍邊框(BoundingBox，矩形邊與WCS的X、Y、Z軸平行)之外，則肯定不在該三角形內；如果在該矩形框內，則再使用內角和法、同向法及重心法等方法進行判斷，也可以使用射線法進行判斷等，但進行排斥試驗之后能夠節約很多程序運行時間。

在建立初始三角網之后，任意取一點，此處以點P13作為示例(如圖2所示)，插入到初始三角網中，經過計算可知點P13位于△P9P12P14中，將點P13與該三角形的3個頂點分別連接，形成新的3個三角形，并加入到三角形存儲鏈表中，新生成的3條邊存入邊鏈表中，并且更新原始邊鏈表中△P9P12P14的3條邊的屬性，如圖2中的虛線所示，同時在三角形鏈表中刪除原始的△P9P12P14[14]。

4) 生成三角網：所有點插入完成之后，將最終的三角形鏈表和邊鏈表代替初始三角網鏈表，則得到最終的三角網。

該部分程序的算法在當數據點排列比較規則、分散比較均勻的情況下，其時間復雜度為O(n)，在最壞的情況下，其時間復雜度為O(n2)，其時間復雜度平均保持在O(nlgn)。

(3) LOP優化

LOP優化即局部優化方法，是Lawson于1977年提出的，它的核心思想是通過交換凸四邊形的對角線來獲得等角性更好的TIN。

1) 預處理：將凸包邊從邊鏈表中剔除，凸包邊是最外圍的邊界，只包含在一個三角形中，故不再用于優化處理。

2) 改進LOP優化方法：LOP優化的目的是滿足空圓特性和最大最小角特性，即保證最鄰近的點構成三角形且三角形盡量接近等邊。文獻[15]中提出了通過比較凸四邊形中對角線的長短來進行三角形優化處理的方法，該方法選取最短的一條對角線作為該凸四邊形的劃分線，但是該法在某些情況下無法處理，如圖3所示。

圖3 改進LOP優化算法示意圖

圖3中，LP2P4=51.71，LP1P3=54.02，如果按照文獻[15]的說法，取最短的邊作為該四邊形的分割線，則應該取圖3中的虛線(對角線)，而實際上應該取圖3中邊P1P3作為該四邊形的分界線。從圖中可以看出，∠P1P4P3為55.96°，∠P1P2P3為102.93°，∠P2P3P4為147.44°，∠P2P1P4為53.67°，如果將∠P1P4P3與∠P1P2P3相加，將∠P2P3P4與∠P2P1P4相加，并將它們各自的和分別與120°相減可以得出：∠P1P4P3+∠P1P2P3-120°<∠P2P3P4+∠P2P1P4-120°。

之所以要與120°相減是因為D-TIN需要滿足最大最小角特性，即三角形盡量接近等邊，等邊三角形中各角均為60°，同時也很好地滿足空圓特性。

從上述討論可以看出，與120°相差最小的那一對角所對的對角線恰好為所選擇的對角線，如圖3中的邊P1P3，該方法同樣適用于其他凸多邊形。

基于上述計算和檢驗筆者提出了角度判別對角線的方法，該方法能夠很好地處理LOP優化，同時該方法的時間復雜度僅為O(n)，能夠很好地完成三角網的局部優化。

3) LOP優化：基于上面所提出的角度判別對角線的方法，使用改進后的LOP優化方法進行三角網優化的處理過程如下：

從邊鏈表中任意取出一條邊，判斷該邊的左三角形和右三角形組成的四邊形是否為凸四邊形，如果是凹四邊形則不進行后續操作而繼續取下一條邊；若是凸四邊形則需要按照介紹的方法進行LOP優化處理，如果交換了對角線，則需要將交換前的兩個三角形從三角形鏈表中刪除，同時加入新生成的三角形，還需要將交換前的邊從邊鏈表中刪除，同時加入新生成的邊到邊鏈表中。

由于設計邊鏈表時存儲了邊的左右三角形，因此進行LOP優化時可以直接取出該兩個三角形進行優化，省去了對左右三角形的搜索時間，同時三角形與三角形之間通過邊進行鏈接，各個三角形之間形成有向鏈接，能夠很方便地進行等值線生成等后續處理。

依次循環檢查各條邊，當邊鏈表循環完畢則LOP優化處理結束，此時輸出的三角形鏈表和邊鏈表為最終Delaunay三角網構網結果。

三、算法運行效果

將本文介紹的方法采用C#語言實現，對某一測區10 235個數據點進行構網，先利用Akl-Toussaint啟發式函數建立凸包，然后再將凸包以外的數據點進行插入，最后利用角度判別對角線法進行LOP優化以生成Delaunay三角網。圖4是原始數據點集，圖5是生成的Delaunay三角網的效果圖。

圖4 原始離散點數據集

圖5 三角網生成后的效果

筆者利用筆記本電腦，處理器為i5，采用本文所提出的方法生成上述三角網的時間(包括數據的讀取時間)為2.67 s，與分治算法的2.70 s基本齊平。在綜合對比了其他數量級的數據生成算法時間之后，得出該算法比較穩定，且唯一性較好，算法的效率大部分情況下與同樣唯一性較好的分治算法相當，偶爾比分治算法高效。

在常用的幾種三角剖分算法中[16-17]，逐點插入法的時間復雜度平均為O(n2)，三角網生長算法的時間復雜度平均也為O(n2)，分治算法的時間復雜度雖然也可以降到O(n)，但是其子塊的遞歸算法的時間復雜度仍然達到了O(n2)。

本算法從凸包的生成、點的插入、三角網的優化等方面入手，采用了Akl-Toussaint啟發式函數、邊的有向性存儲和角度判別對角線法進行LOP優化等，大大降低了程序的時間復雜度，使時間復雜度平均保持在O(nlgn)，在最好的情況下能夠達到O(n)，在最壞情況下也小于O(n2)，明顯小于常用三角剖分的時間復雜度。

四、結束語

本文將凸包法與逐點插入法相結合，提出了對于三角剖分算法的改進方法。此法從凸包生成、三角網存儲結構、三角形的快速定位，以及LOP優化等方面對算法進行優化，應用Akl-Toussaint啟發式函數、三角形邊的有向性存儲，以及角度判別對角線法等策略，大幅度提高算法的運行效率。從算法分析和程序運行效果可以看出，該算法不僅能夠完全滿足大數據量的Delaunay三角剖分對高效和高精度的需求，而且其算法比分治算法更加容易理解、容易實現。

參考文獻：

[1] 胡鵬,楊傳勇,吳艷蘭，等.新數字高程模型:理論、方法、標準和應用[M].北京：測繪出版社，2007.

[2] 張宏,溫永寧,劉愛利，等.地理信息系統算法基礎[M].北京：科學出版社，2006.

[3] 凌海濱,吳兵.改進的自連接Delaunay三角網生成算法[J].計算機應用，1999，19(12):10-12.

[4] LAWSON C L. A Software for C Surface Interpolation[J].Mathematical Software，1977(3): 161-194.

[5] LEE D T, SCHACHTER B J. Two Algorithms for Constructing a Delaunay Triangulation[J]. International Journal of Computer and Information Science, 1980,9(3):219-242.

[6] BOWYER A. Computing Dirichlet Tesselations[J]. Computer Journal,1981(24): 162-166.

[7] MACEDONIA G, PARESCHI M T. An Algorithm for the Triangulation of Arbitrarily Distributed Points: Applications to Volume Estimate and Terrain Fitting[J]. Computers & Geosciences,1991(17): 859-874.

[8] FLORIANI L D， PUPPO E. An On-line Algorithm for Constrained Delaunay Triangulation[J].CVGIP: Graphical Models and Image Processing, 1992, 54(3): 290-300.

[9] 武曉波,王世新,肖春生.Delaunay 三角網的生成算法研究[J].測繪學報,1999，28(1):28-35.

[10] MCCULLAGH M J.Delaunay Triangulation of a Random Data Set for Lirarithmic Mapping [J].The Cartographic Journal,1980(17):93-99.

[11] MIRANTE A,WEINGARTEN N.The Radical Sweep Algorithm for Constructing Triangulated Irregular Networks[J]. IEEE Computer Graphics and Application,1982，2(3):11-21.

[12] 芮一康,王結臣.Delaunay三角形構網的分治掃描線算法[J].測繪學報,2007，36(3):358-362.

[13] HEINEMAN G T, POLLICE G，SELKOW S.算法技術手冊[M].楊晨,李明，譯.北京：機械工業出版社，2010.

[14] 徐道柱,劉海硯.Delaunay三角網建立的改進算法[J].測繪與空間地理信息,2007，30(1)：38-41.

[15] 魏向輝,夏春林,魯慶.一種基于凸包的Delaunay三角網算法設計[J].測繪科學,2010，35(5):152-153.

[16] 邵春麗,胡鵬,黃承義,等.Delaunay三角網的算法詳述及其應用發展前景[J].測繪科學，2004，29(6):68-71.

[17] 余杰,呂品,鄭昌文.Delaunay三角網構建方法比較研究[J].中國圖象圖形學報,2010，15(8):1158-1167.