帶斷層約束的Delaunay三角剖分混合算法＊

張群會，解子毅

（西安科技大學計算機科學與技術學院，陜西 西安 710054）

0 引言

目前比較流行的構建數(shù)據(jù)高程模型(DEM)的方法主要有2種:一種是基于不規(guī)則三角形(Triangular Irregular Networks，TIN)網(wǎng)格法，另一種是四邊形(GRID)網(wǎng)格法。和GRID相比，TIN能夠精確地表達網(wǎng)格邊界和斷層，是一種比較理想的三維層面建模方法，在地質(zhì)層面可視化方面得到了廣泛應用。Delaunay三角剖分法是建立TIN模型最常用的方法。由于Delaunay三角化滿足最小角最大準則和外接圓不包含其他點準則，總是能盡可能避免狹長的三角形，自動向等邊三角形逼近，具有網(wǎng)格形態(tài)優(yōu)美等特點。因此利用Delaunay三角剖分來建立帶斷層約束的地質(zhì)模型具有十分重要的價值。

現(xiàn)階段的 CDT生成算法較多，主要以2步法［1－2］為主，其基本思想是第一步構建 DT網(wǎng)格；第二步將約束線段強行嵌入到DT網(wǎng)格中進行整體剖分。傳統(tǒng)CDT兩步法在約束線段嵌入方面的方法很多，比較流行的有3種:一種［1］是利用生長法和分治法［3］的思想，首先刪除與約束線段相交的線段，形成影響域。然后以約束線段為基邊把影響域分成左右2部分分別進行剖分。由于約束線段的長度沒有進行細分，構網(wǎng)時在斷層邊界處可能出現(xiàn)狹長三角形。解決上述問題可以采用盧朝陽的對半劃分增量型附加點插入算法［4］，通過對約束線段上插入數(shù)據(jù)點來優(yōu)化三角網(wǎng)格。第一種算法簡單但效率不高，網(wǎng)格形態(tài)也不優(yōu)美。第二種是李立新優(yōu)化的多對角線交換算法［5］。該算法的思想是對影響區(qū)域的三角形組成的凸四邊形進行對角線交換，通過多次交換對角線可實現(xiàn)在任意影響域中強行嵌入約束線段。但該算法執(zhí)行效率較低、編程難度大，并且也會出現(xiàn)狹長三角形的情況。第三種是基于復連通區(qū)域的三角剖分算法［6］。該算法首先要搜索到與約束線段相交或靠近約束線段的三角形，并把這些三角形合并成一個單連通區(qū)域，即找到包含約束線段的最小外環(huán)。然后求取該區(qū)域的外邊界，將最小外環(huán)和其外邊界組成復連通區(qū)域，最后將復連通區(qū)域三角剖分轉(zhuǎn)化成了簡單多邊形的三角剖分。該方法雖然能夠有效的解決約束線段的嵌入問題，但是由于其步驟較多，構網(wǎng)效率不高。而文中的三角剖分混合算法可以較好的解決上述問題。

1 DT的基本性質(zhì)

1.1 空外接圓性質(zhì)

在由點集P所形成的DT網(wǎng)格中，其每個三角形的外接圓均不包含P中的其他任意點。在二維空間中，給定點集P，若P中任意4個點不共圓，則Delaunay三角化是唯一的，否則Delaunay三角化不唯一。對于Delaunay三角化，局部優(yōu)化可以保證全局最優(yōu)，如圖1(a)。

1.2 最小角度最大性質(zhì)

在由點集P所能形成的三角網(wǎng)中，DT網(wǎng)格中三角形的最小角度是最大的。即在三角形網(wǎng)格中，對任意相鄰的、構成嚴格凸四邊形的2個三角形來說，其6個內(nèi)角中的最小值大于交換四邊形的對角線所形成的另外2個三角形的6個內(nèi)角中的最小值，如圖1(b)。

圖1 Delaunay三角化的特點Fig.1 Characteristics of the Delaunay triangulation

2 數(shù)據(jù)結構

在插入點及約束邊的嵌入過程中，需要由三角網(wǎng)拓撲信息快速搜索到插入點和約束線段所影響的三角形，并不斷的更新三角網(wǎng)的拓撲信息。因此，保證三角網(wǎng)拓撲信息的正確性和有效性對于建立CDT是至關重要的。文中將三角網(wǎng)內(nèi)查詢的拓撲信息主要定義為點、邊、三角形3個方面，其各自的結構如下:

3 帶斷層約束的三角剖分混合算法

傳統(tǒng)CDT剖分算法所耗費的構網(wǎng)時間主要用在斷層數(shù)據(jù)插入和約束線段嵌入這2方面上。該算法由于要不斷的判斷插入點是否包含在三角形的外接圓內(nèi)，因此需要耗費大量的時間來計算三角形的外接圓。而插入點所影響的三角形一般都是它周圍鄰近的三角形，沒有必要對所有的三角形進行處理。針對上述問題，文中提出的混合算法在插點方面對傳統(tǒng)算法進行了優(yōu)化:首先以插入點作為中心，然后在以搜索半徑為圓的范圍內(nèi)尋找受插入點影響的三角形(只要搜索半徑設置合理，就一定能夠找到所有受插入點影響的三角形)，并只對搜索到的三角形進行處理，從而減少了處理過程，大大提高了構網(wǎng)效率。而在約束線段嵌入方面，文中算法首先利用均分法對約束線段進行細分，解決了約束線段過長而產(chǎn)生狹長三角形的問題，并且還能夠使約束線段盡量包含在DT網(wǎng)格中，減少后續(xù)過程對約束線段的嵌入處理。

文中分3步來實現(xiàn)帶斷層約束的三角剖分混合算法。第一步首先建立一個包含所有數(shù)據(jù)點的矩形框，用逐點插入法［7－8］將數(shù)據(jù)點逐一插入到矩形框中進行剖分，然后利用重心布點法［9］對剖分好的網(wǎng)格進行加密，最后通過局部優(yōu)化法(LOP)［10］對剖分好的三角形進行優(yōu)化，形成無約束的DT網(wǎng)格；第二步首先對斷層數(shù)據(jù)進行細分，并以細分好的數(shù)據(jù)作為中心，在一定的范圍內(nèi)對DT網(wǎng)格進行局部剖分，直到所有斷層點嵌入到DT網(wǎng)格中。然后連接相鄰的斷層點，形成斷層線段，判斷其是否包含在網(wǎng)格中，如果不包含，采用三角網(wǎng)生長法和分治法將斷層線段嵌入到網(wǎng)格中，直到所有斷層線段處理完畢，最后用局部優(yōu)化法優(yōu)化網(wǎng)格；第三步將由斷層線段圍成的斷層區(qū)域中的三角形刪除，形成帶斷層約束的CDT網(wǎng)格。

3.1 構造無約束的DT網(wǎng)格

1)設邊界點數(shù)組 P=(p1，p2，…，pn)，n 為邊界數(shù)據(jù)點個數(shù)，三角形數(shù)組 T=(t1，t2，tm)，m 為三角形個數(shù)。構造初始矩形框，用逐點插入法在矩形框中插入數(shù)據(jù)點pj，找到T中外接圓含pj的三角形ti，并將ti從T刪除，將ti的相鄰邊刪除形成多邊形空腔，把pj與空腔多邊形各頂點相連，生成新三角形存入到T中。

2)用重心布點法對網(wǎng)格進行加密，并用局部優(yōu)化法(LOP)進行優(yōu)化。重心布點法的思想是首先設置三角形的最大高度為H，并通過在三角形重心位置布點，使所生成的三角網(wǎng)相對均勻，且三角形高度均不超過H.

3.2 構造帶斷層約束的CDT網(wǎng)格

1)利用細分算法把由斷層數(shù)據(jù)構成的斷層邊進行細分，將細分結果有序的存儲在數(shù)組F中。具體步驟如下

步驟1 從P中取出相鄰的2點pi和pi+1，i=(1，2，…，n －1)，計算出 pi，pi+1的間距 di，并設置一定的步長Step，如果Step＜d轉(zhuǎn)到步驟2，否則繼續(xù)步驟1；

步驟2 設D=d/Step，將線段pipi+1均分為D+1段，將分段點按順序存放在F中；

2)將細分后的斷層數(shù)據(jù)點作為新點逐一插入到DT網(wǎng)格中進局部剖分，使斷層點嵌入到DT網(wǎng)格中。具體步驟如下

步驟1 根據(jù)三角形的最大高度H來設置搜索半徑r；

步驟2 從F中任取一點fi∈F，計算fi到ti∈T 3個頂點的距離，取其中最短距離作為di，如果di＜r，將ti存入到臨時數(shù)組V中。

步驟3 判斷fi是否在V中三角形的外接圓內(nèi)，如果在，將三角形進行標記；

步驟4 刪除標記的三角形鄰邊，形成多邊形空腔，將fi與多邊形各頂點進行V連接形成新的三角形存入到T中，清空V.迭代2，3，4步驟，直到所有斷層數(shù)據(jù)點處理完畢。

圖2為剖分過程圖(以f1∈F為例):從圖2(a)可知f1C為搜索半徑r，C為圓上的任意一點，f1到 t6，t7，t8，t9的距離為 d1=f1A，f1到 t1，t2，t3，t4，t5的距離為 d2=f1B，且 d1＜ r，d2＜ r因此 t1，t2，t3，t4，t5，t6，t7，t8∈V 為 f1鄰近三角形；從圖 2(b)可知，f1在 t4，t5，t6的外接圓內(nèi)，將 t4，t5，t6從 T 中刪除；從圖2(c)可知，刪除 t4，t5，t6的相鄰邊所形成多邊形空腔為Q={ABCDE}，將Q的頂點與f1進行連接，形成新的三角形 t10，t11，t12，t13，t14存入到 T中，完成該點局部剖分。

圖2 空外接圓法剖分過程圖Fig.2 Empty circum circle subdivision process diagram

3)將F中的第一個元素存入F，實現(xiàn)首尾相連，并采用生長法和分治法將斷層線段嵌入到DT網(wǎng)格中。具體步驟如下:

步驟1 取F中相鄰的兩點形成約束線段fifi+1，i≤n －1，判斷 fi，fi+1是否包含在 DT 網(wǎng)格中，如果包含繼續(xù)步驟1，否則轉(zhuǎn)到步驟2；

步驟2 找到與fifi+1的相交邊，將相交邊刪除形成多邊形空腔，并以fifi+1為基邊將多邊形空腔分成左右2個區(qū)域；

步驟3 利用三角網(wǎng)生長法對QL和QR2個區(qū)域分別進行剖分(以QL為例)。

1)建立一堆棧S，并把fifi+1邊壓入S；

2)在S中彈出一邊為fifi+1作為基邊，在QL頂點集合中找出距fifi+1最近的點A進行連接得到fi+1A，將三角形fifi+1A存入T中；

3)如果fiA或fi+1A是QL的一條邊界邊，則不入棧，反之把fiA或fi+1A壓入堆棧S中；

4)如果S不空，則返2，否則退出；

5)對QL和QR剖分后形成的三角形進行局部優(yōu)化(LOP)處理；

步驟4 重復上述步驟直到所有約束線段處理完畢。刪除F中的最后一個元素；

從圖3(a)中可以看出 fifi+1與 AD，BD，BC相交，刪除線段 AD，BD，BC形成多邊形空腔 Q={fiABfi+1CD}.圖3(b)中約束線段fifi+1將多邊形空腔分成QL={fiABfi+1}和QR={fifi+1CD}2個區(qū)域。圖3(c)中QL和QR分別進行了剖分，最終實現(xiàn)了約束線段的嵌入。

4)將斷層區(qū)域中的三角形刪除，最終形成帶斷層約束的CDT網(wǎng)格。具體步驟如下

圖3 影響域的形成和剖分過程Fig.3 Formation of influence domain and subdivision process

步驟1 首先利用F中的斷層數(shù)據(jù)點建立斷層多邊形 Q=(f1，f2，…，fn)；

步驟2 計算T中三角形ti的重心gi，判斷gi是否在Q內(nèi)，如果在就從T中刪除ti，重復步驟2，直到所有三角形處理完畢。

4 實驗結果與分析

文中算法在Microsoft Visual C++2010平臺上編寫測試通過，實現(xiàn)了帶斷層約束的Delaunay三角剖分。對Release版本進行了時間性能測試，測試環(huán)境為:操作系統(tǒng)Windows 7旗艦版(64位)；處理器:AMD A8－4500 M with Radeon(tm)HD Graphics 44 CPU 1.9 GHz；內(nèi)存:4 GB.圖4為文中算法構造出的(a)不帶斷層的DT網(wǎng)格和(b)帶斷層約束的CDT網(wǎng)格，其時間性能見表1.

表1 算法執(zhí)行效率測試結果Tab.1 Test results of algorithm performance

圖4 混合算法剖分效果圖Fig.4 Mixed algorithm subdivision renderings

從圖4可以看出，混合CDT剖分算法所構造的DT網(wǎng)格和CDT網(wǎng)格形態(tài)都比較均勻優(yōu)美，沒有出現(xiàn)狹長三角形的情況，并且斷層邊能夠完好的嵌入到網(wǎng)格中。圖4(b)中的黑點為斷層數(shù)據(jù)點。通過分析表1可得，在剖分相同數(shù)量的三角形情況下，文中算法在構建DT和CDT網(wǎng)格的時間效率上都明顯好于傳統(tǒng)算法。

5 結語

文中通過對幾種傳統(tǒng)算法進行分析和研究，逐一找出了各算法在構網(wǎng)過程中耗時最多的算法步驟，并在深入了解其核心思想的基礎上，通過實驗對比，判斷出幾種傳統(tǒng)算法的優(yōu)缺點，并根據(jù)研究結果將各算法的優(yōu)點進行結合和優(yōu)化，提出了一種帶斷層約束的Delaunay三角剖分混合算法。該算法通過對斷層數(shù)據(jù)的插入和約束線段的嵌入過程進行優(yōu)化，快速實現(xiàn)了帶斷層約束的CDT網(wǎng)格。最后通過實驗結果的分析和對比，說明文中算法在構網(wǎng)效率和構網(wǎng)質(zhì)量上都優(yōu)于傳統(tǒng)算法，在數(shù)據(jù)建模方面具有較好的應用性和實用價值。

References

［1］劉學軍，龔健雅.約束數(shù)據(jù)域的Delaunay三角剖分與修改算法［J］.測繪學報，2001，30(1):82 －88.LIU Xue-jun，GONG Jian-ya.Constrained Data Delaunay triangulation and modification algorithms［J］.Acta Geodaetica et Cartographica Sinica，2001，30(1):82 －88.

［2］劉少華，程朋根，史文中.約束Delaunay三角網(wǎng)生成算法研究［J］.測繪通報，2004，26(1):4 －7.LIU Shao-hua，CHENG Peng-gen，SHI Wen-zhong.Constrained delaunay triangulation algorithm research［J］.Bulletin of Surveying and Mapping，2004，26(1):4 －7.

［3］謝增廣.平面點集帶Delaunay三角剖分的分治算法［J］.計算機工程與設計，2012，33(7):70 －73.XIE Zeng-guang.Partition of planar point set with delaunay triangulation subdivision algorithm［J］.Computer Engineering and Design.2012，33(7):70 －73.

［4］盧朝陽，吳成柯，周幸妮.滿足全局Delaunay特性的帶特征約束的散亂數(shù)據(jù)最優(yōu)三角剖分［J］.計算機學報，1997，20(2):118 －124.LU Chao-yang，WU Cheng-ke，ZHOU Xing-ni.Meet the global feature of constraint delaunay properties of approximate optimal triangle subdivision［J］.Chinese Journal of Computers，1997，20(2):118 －124.

［5］李立新，譚建榮.約束Delaunay三角剖分中強行嵌入約束邊的多對角線交換算法［J］.計算機學報，1999，22(10):1114－1118.LI Li-xin，TAN Jian-rong.Constraint delaunay triangulation subdivision forcibly embedded constraints in multiple diagonal exchange algorithm［J］.Chinese Journal of Computers，1999，22(10):1114 －1118.

［6］劉少兵.帶斷層地震數(shù)據(jù)的Delaunay三角剖分算法［J］.煤田地質(zhì)與勘探，2008，36(6):71 －73.LIU Shao-bing.With seismic data of fault delaunay triangular subdivision algorithm［J］.Coal Geology and Exploration，2008，36(6):71 －73.

［7］Watson D F.Computing the n-dimension delaunay tessellation with application to voronoi polylopes［J］.Computer Journal，1981，24(2):167 －172.

［8］周雪梅，黎應飛.基于Bowyer-Watson三角網(wǎng)生成算法的研究［J］.計算機工程與應用，2013，49(6):198－201.ZHOU Xue-mei，LI Ying-fei.Based bowyer-watson triangulation generation algorithm research［J］.Computer Engineering and Applications，2013，49(6):198 －201.

［9］吳 芬.平面域Delaunay三角剖分新加密算法［J］.計算機與現(xiàn)代化，2007(7):19－22.WU Fen.Plane domain delaunay triangular subdivision new encryption algorithm［J］.Computer and Modernization，2007(7):19 －22.

［10］Lawson C L.Software for C1surface interpolation mathematical software Ⅲ［M］.NewYork:Academic Press，1977.