一種Delaunay三角剖分的改進(jìn)算法

余代俊，蒲朝旭，朱逍賢

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

一、引 言

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

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

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

二、算法設(shè)計(jì)

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

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

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

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

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

{

double X;∥X坐標(biāo)

double Y;∥Y坐標(biāo)

double Z;∥Z坐標(biāo)

}

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

{

Point3d StartPt;∥起點(diǎn)

Point3d EndPt;∥終點(diǎn)

Triangle LeftTri;∥左三角形

Triangle RightTri;∥右三角形

}

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

{

Edge edge1;∥邊1

Edge edge2;∥邊2

Edge edge3;∥邊3

}

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

2. 算法實(shí)現(xiàn)

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

(1) 生成凸包

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

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

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

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

3) 排序點(diǎn)：對(duì)進(jìn)行清理之后的點(diǎn)集按照X坐標(biāo)升序?yàn)橹?、Y坐標(biāo)升序?yàn)檩o的方式進(jìn)行排序，以便用于后續(xù)凸包的生成。對(duì)于點(diǎn)集的排序可以采用LINQ方式進(jìn)行，該方法是所有排序算法中最高效的。

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

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

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

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

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

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

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

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

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

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

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

(3) LOP優(yōu)化

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

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

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

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

圖3中，LP2P4=51.71，LP1P3=54.02，如果按照文獻(xiàn)[15]的說法，取最短的邊作為該四邊形的分割線，則應(yīng)該取圖3中的虛線(對(duì)角線)，而實(shí)際上應(yīng)該取圖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°相減是因?yàn)镈-TIN需要滿足最大最小角特性，即三角形盡量接近等邊，等邊三角形中各角均為60°，同時(shí)也很好地滿足空?qǐng)A特性。

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

基于上述計(jì)算和檢驗(yàn)筆者提出了角度判別對(duì)角線的方法，該方法能夠很好地處理LOP優(yōu)化，同時(shí)該方法的時(shí)間復(fù)雜度僅為O(n)，能夠很好地完成三角網(wǎng)的局部優(yōu)化。

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

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

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

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

三、算法運(yùn)行效果

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

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

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

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

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

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

四、結(jié)束語

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

參考文獻(xiàn)：

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

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

[3] 凌海濱,吳兵.改進(jìn)的自連接Delaunay三角網(wǎng)生成算法[J].計(jì)算機(jī)應(yīng)用，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].測(cè)繪學(xué)報(bào),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].測(cè)繪學(xué)報(bào),2007，36(3):358-362.

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

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

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

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

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