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

余代俊，蒲朝旭，朱逍賢

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

一、引 言

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

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

本文就如何改進和優(yōu)化Delaunay三角網(wǎng)的構(gòu)建算法來滿足海量數(shù)據(jù)的處理需求，以及如何將構(gòu)建算法利用較為簡潔的程序進行實現(xiàn)，進行了深入研究和詳細探討。

二、算法設計

筆者結(jié)合已有三角網(wǎng)的生成算法[4-12]，取長補短，提出了基于凸包實現(xiàn)的Delaunay三角網(wǎng)剖分算法。

1. 數(shù)據(jù)存儲結(jié)構(gòu)

一個好的數(shù)據(jù)結(jié)構(gòu)直接影響著D-TIN的生成效率，D-TIN的數(shù)據(jù)結(jié)構(gòu)主要包括點、線、面之間的組織關系和拓撲關系。

本文所使用的數(shù)據(jù)結(jié)構(gòu)如下

struct Point3d∥點結(jié)構(gòu)

{

double X;∥X坐標

double Y;∥Y坐標

double Z;∥Z坐標

}

struct Edge∥邊結(jié)構(gòu)

{

Point3d StartPt;∥起點

Point3d EndPt;∥終點

Triangle LeftTri;∥左三角形

Triangle RightTri;∥右三角形

}

struct Triangle∥三角形結(jié)構(gòu)

{

Edge edge1;∥邊1

Edge edge2;∥邊2

Edge edge3;∥邊3

}

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

2. 算法實現(xiàn)

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

(1) 生成凸包

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

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

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

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

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

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

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

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

(2) 建立三角網(wǎng)

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

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

圖2 Delaunay三角網(wǎng)建立過程示意圖

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

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

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

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

(3) LOP優(yōu)化

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

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

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

圖3 改進LOP優(yōu)化算法示意圖

圖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優(yōu)化，同時該方法的時間復雜度僅為O(n)，能夠很好地完成三角網(wǎng)的局部優(yōu)化。

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

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

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

依次循環(huán)檢查各條邊，當邊鏈表循環(huán)完畢則LOP優(yōu)化處理結(jié)束，此時輸出的三角形鏈表和邊鏈表為最終Delaunay三角網(wǎng)構(gòu)網(wǎng)結(jié)果。

三、算法運行效果

將本文介紹的方法采用C#語言實現(xiàn)，對某一測區(qū)10 235個數(shù)據(jù)點進行構(gòu)網(wǎng)，先利用Akl-Toussaint啟發(fā)式函數(shù)建立凸包，然后再將凸包以外的數(shù)據(jù)點進行插入，最后利用角度判別對角線法進行LOP優(yōu)化以生成Delaunay三角網(wǎng)。圖4是原始數(shù)據(jù)點集，圖5是生成的Delaunay三角網(wǎng)的效果圖。

圖4 原始離散點數(shù)據(jù)集

圖5 三角網(wǎng)生成后的效果

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

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

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

四、結(jié)束語

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

參考文獻：

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

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

[3] 凌海濱,吳兵.改進的自連接Delaunay三角網(wǎng)生成算法[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 三角網(wǎng)的生成算法研究[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] 芮一康,王結(jié)臣.Delaunay三角形構(gòu)網(wǎng)的分治掃描線算法[J].測繪學報,2007，36(3):358-362.

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

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

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

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

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