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

張群會，解子毅

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

0 引言

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

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

1 DT的基本性質

1.1 空外接圓性質

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

1.2 最小角度最大性質

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

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

2 數據結構

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

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

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

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

3.1 構造無約束的DT網格

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

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

3.2 構造帶斷層約束的CDT網格

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

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

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

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

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

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

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

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

圖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的外接圓內，將 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，實現首尾相連，并采用生長法和分治法將斷層線段嵌入到DT網格中。具體步驟如下:

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

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

步驟3 利用三角網生長法對QL和QR2個區域分別進行剖分(以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剖分后形成的三角形進行局部優化(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個區域。圖3(c)中QL和QR分別進行了剖分，最終實現了約束線段的嵌入。

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

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

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

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

4 實驗結果與分析

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

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

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

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

5 結語

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

References

［1］劉學軍，龔健雅.約束數據域的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三角網生成算法研究［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特性的帶特征約束的散亂數據最優三角剖分［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］劉少兵.帶斷層地震數據的Delaunay三角剖分算法［J］.煤田地質與勘探，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三角網生成算法的研究［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］.計算機與現代化，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.