一種含斷層數據的變尺度加密三角剖分算法

史曉楠 王夢凡 王子童

(西安科技大學計算機科學與技術學院 陜西 西安 710054)

0 引 言

Delaunay三角網，是由一系列相連的但并不重合的三角形組成的一個集合，可以用來表示線性特征和迭加任意形狀的區域邊界，與不規則的地面特征很協調。近年來，鑒于Delaunay三角網滿足最小角度最大、任意四點不會共圓等特征，被大量使用在有限元分析和信息可視化領域。長期以來，在含斷層數據地形三角剖分方面有兩大主流思想，即先構建還是先加密。學者們提出了很多高價值的算法，其中最常用的三種剖分算法就是逐點插入法、生長法和分治算法[2]。許多學者對這三種算法進行了深入研究，并做出了相應的改進。袁小翠等[3]在傳統的生長法上提出了消去遞歸的生長法和基于微切平面的生長法，實驗證明三種生長法剖分效果近似，但是改進的兩種生長法的計算速度卻有著非常大的提高。在網格索引方面，姜志偉等[4]提出的基于格網索引和方向法搜索的算法很大程度上提高了構網的速率。在點定位算法上，國內外學者有著較多的研究和改進，其中最常用的算法有鐘正等[5]提出的關于在四面體體積坐標基礎上的點定位算法，李小秋等[6]提出的點定位算法，但是這種算法思路比較繁瑣，執行效果不佳。李剛等[7]通過設置角度和長度閾值以及結合特征點優化建立具有新特征約束的三角網，該方法在擬合方面較為逼真，但是閾值的盲目選擇會導致數據冗余和原始數據喪失。綜上所述，本文在逐點插入法上應用了網格索引，通過設置剖分尺度有效地解決了上述難以解決的斷層數據的合理嵌入問題。

本文主要研究的是包含斷層的數據如何較高效率且高質量的實現Delaunay三角剖分的問題。基于兩步法，本文首先使用逐點插值方法和可變尺度加密算法生成無約束DT(Delaunay Triangulation)網格，然后先用細分算法細分斷層邊，避免了狹長三角形的出現，再結合生長法和分治法思想構建帶斷層約束的CDT(Constrained Delaunay Triangulation)網格。對比其他常見三角剖分方法，本文方法可以較好地解決斷層線段的難以合理的插入的問題，有效提高了構網質量和效率。

1 算法整體流程

算法整體流程如圖1所示。算法首先融合所有數據(包括斷層約束數據)，在網格索引的基礎上，通過遞歸算法，讀取數據形成凸殼；其次測試出邊界多邊形邊長的約束值，用來構建更高精度的凸殼[8]；關于點位置的確立，這里采用面積法進行構建；遞歸確定最佳的加密尺度，然后完成整體網格的加密；最后將細化后的斷層數據用生長法嵌入網格，其中定位算法采用高效率的余弦算法。最終實現了含斷層線段合理嵌入的三角剖分。

圖1 算法流程圖

2 關鍵算法設計

2.1 無約束三角網的形成

(1)

(2)

(3)

S=(Xmax-Xmin)×(Ymax-Ymin)

(4)

將所有離散點按照上述規則投影到網格，如圖2所示。

圖2 虛擬網格

(2) 建立區域凸包。本文借鑒Larkin改進和完善的凸殼算法生成凸包，并在此基礎上“邊界收縮”，建立高精細度的凸包，更貼切斷層數據。在這里，建立一個存放所生成的凸包點的鏈表，且這些點按照逆時針順序存放，具體的數據結構是點的坐標(x，y)。具體步驟如下[10]：

Step1將具有x,y,x+y,x-y的極大值和極小值點連接，圍成最初凸包，如圖3所示。

圖3 凸殼初步建立

Step2選取相鄰的兩個點pi、pi+1，查詢pipi+1右側距離最大的點，如果該點在pipi+1之間，將其插入，重復上述過程，直到對鏈表中任意兩點構成的線段都判斷一次。

查詢出線段右側距離最大的點，可以結合網格和三角形面積判斷。已知三點坐標，由面積公式得面積最大對應的點(x3,y3)：

(5)

Step3當鏈表中頂點個數等于每次沒有插入點的次數，停止操作，初始凸包成功建立。

(3) 高邊界精度凸殼的生成。

Step1設置最大凸殼邊界約束值。

Step2在已知凸殼中遞歸搜尋比最大邊長長的線段，若存在，得到其集合M。

Step3依次取出集合M中的線段AB，在線段附近搜尋距離線段最近的點J。

Step4若J點存在，且若AJ、BJ符合邊界約束值，繼續取出下一條線段，否則繼續搜尋次之距離最近的點，重復Step 3。最終結果如圖4所示，精度提高。

圖4 “高精細度”凸殼

(4) 初始三角剖分。任取凸殼內一點與凸殼點進行連接，構成三角網的雛形。隨機選取凸殼內任意一點，依次連接該點和圖中所有的點進行圖5所示的初始三角剖分。

圖5 初始三角剖分

(5) 逐點插入完成Delaunay三角網。主要思想是在初始三角網的基礎上，將余下的點逐一插入，最后使用LOP算法來確定這些點組成Delaunay三角網。

Step1隨機選取最初三角形，開始檢索。

Step2將剩余的點依次插入，插入過程中需要判斷該點的位置，找出包含它的三角形T，刪除T的公共邊。

這里給定點p(x0,y0)，判斷三角形T中是否包含點P，設三頂點坐標分別為A(x1,y1)、B(x2,y2)、C(x3,y3)，S1為三角形ABC，S2為三角形PBC，S3為三角形APC，為三角形ABP，使用面積法判斷：如在式(6)中，如果形成的三個小三角形的面積的絕對值之和等于初始大三角形面積的絕對值，則意味著點p在內部或者邊上。

(6)

反之，若出現如下所示的不等關系，表示點在外部。

abs(S1)

(7)

如圖6所示，S1為0，表示P在BC上，S2為0，表示P在AC上，S3為0，表示P在AB上。

圖6 面積法判斷點的位置

Step3找出T之后，形成新三角形T1、T2、T3。如圖7所示。

圖7 三個新的三角形的形成

Step4LOP算法優化[11]。

Step5重復進行Step 1，當所有點都完成插入時，三角網構建完成。

圖8是逐點插入法的流程框架。

圖8 逐點插入法框架圖

2.2 變尺度加密算法

由于斷層數據雜亂無章且分布較為稠密，為了使三角剖分得到的三角網更加規整，即方便貼合實際地形進行處理，需對整個網格進行更高的加密。本文提出了變尺度加密方案。設置一定尺度，并逐步縮小，一步一步規格化網格。

搜索形成的三角形網格，只要三角形的高度超過分割尺寸H，取邊緣的中點pi并連接以形成兩個新的三角形[12]。

Step1采用變尺度的算法，先設定某個剖分尺度H。

Step2最初用H×N(N為某個自然數)進行全局細化。

Step3在使用H×(N-1)依次到N為1的時候產生最后的剖分結果，此時形成的無約束三角網具有較好的尺度[13]。如下式所示：

(8)

Step4迭代Step 1-Step 3，直到所有三角形符合H，此時，停止剖分。

變尺度加密前后的三角網如圖9所示。

圖9 變尺度加密前后三角網

2.3 斷層數據的嵌入

(1) 斷層數據細分算法。連接所有斷層數據，形成斷層邊，細分所有斷層邊，將形成的點存入數組K[14]。

Step1從離散點集合v中取出相鄰的兩點vi和vi+1，i=(1,2，…，n-1)，首先判斷線段vivi+1是否是已構成的CDT中的三角形的一條邊，若是，則退出直接處理下一條邊，反之，進行下一步工作。

Step2計算出vi和vi+1的間距L，并在初始設置一定的步長d。如果d

(2) 插入細分后斷層數據點。細分后，將獲得的點插入網格。

Step1根據變尺度算法中的最大高度H來設置搜索半徑r。

Step2從K中隨機取一點ki，ki到三角形ti三個頂點的最短的距離賦值給di，如果di

Step3將ki與E中三角形外接圓進行判斷，不符合則標志此三角形。

Step4刪除步驟3中標記的三角形的相鄰邊，并將ki連接到每個點將得到的三角形存儲在T中。清空E，迭代Step 2、Step 3、Step 4，直到斷層數據全部插入網格。如圖10所示，首先放入新節點結點P，再決定P點該如何連接周圍點，判斷P點的相對位置，根據Delaunay三角網性質，刪除公共邊AB，連接形成新的三角形。

(a) 插入新結點P (b) 決定P如何和其他結點相連

(c) 刪除AB邊 (d) 形成新三角形圖10 插入斷層數據點

Step 3中，點定位算法具體步驟如下：

假設P在BC的右側，∠A+∠P≥π，則P在圓內，∠A+∠P<π，則P在圓外。

點P在圓內，根據余弦公式可轉化為:

-(AC×AB)×|PC|×|PB|≥(PC×PB)|AC|×|AB|

(9)

點P在圓外，則:

-(AC×AB)×|PC|×|PB|<(PC×PB)×|AC|×|AB|

(10)

如果P點在BC左側，則判斷的結論剛好相反，如圖11所示，(a)中，點ki在線段BC右側，根據式(9)計算可知，該點在三角形ABC外接圓圓內。同理可知，ki+1在外接圓圓外，ki在外接圓圓上。

(a) (b) (c)圖11 點與外接圓位置關系

(3) 斷層線段嵌入網格。將ki賦值給ki+1，存入數組K中，判斷K中線段是否是網格中線段，若不是，采用生長法嵌入約束[10]。

Step1從數組K依次取出ki、ki+1，i≤n-1，判斷kiki+1是否在網格中，如果在，繼續Step 1，否則轉至Step 2。

Step3用一個隊列Q存儲數組K中的點，在影響域內，先以kiki+1為基邊(以右邊為例)，向右擴展，搜索P，使得角Pkiki+1最大，即構成了Pkiki+1。

Step4分別以Pki、Pki+1為基礎反復Step 3繼續產生擴展，如果擴充邊(比如kiP或者ki+1P)是外邊界的邊或者已經使用了2次，那么就停止這條邊往外擴張，將三角形納入T中。

Step5迭代上述步驟直到隊列Q為空，并且完成受影響的線段的三角網的重建。

Step6搜索重建后的所有三角形，進行LOP算法優化。如果存在三角形共圓，則交換對角線。

3 算法分析

表1是常見的三角剖分各算法時間復雜度比較[15]。

綜上，本文無論是在逐點插入法、點定位算法還是約束區域的嵌入上，算法的時間復雜度均較低。

4 實 驗

文中算法在Microsoft Visual Studio 2015平臺上編寫測試通過，實現了含斷層數據的變尺度加密三角剖分。對Debug版本進行了時間性能測試，測試環境為：Window 10 64位操作系統，Intel i5-8300H CPU 3.89 GHz，16 GB 內存，GTX1060 6 GB顯卡。傳統的生長法[16]、逐點插入法[17]和含斷層數據的變尺度加密三角剖分的算法執行效率測試見表2。

表2 算法執行效率測試結果

圖13為含斷層數據的變尺度加密三角剖分效果圖。可以看出，(a)是未經過變尺度加密的三角網，不含斷層數據；(b)是經過變尺度加密后的三角剖分結果，斷層細節比較簡單；(c)是在(b)的基礎上進一步進行變尺度加密的三角網，斷層顯示較為明顯。同時從表2可以看出，本文使用的算法在三角剖分的時間效率上優于傳統算法。

(a) (b) (c)圖13 含斷層數據的變尺度加密三角剖分效果圖

5 結 語

本文通過改進傳統的逐點插入法和約束邊的嵌入，提出了一種含斷層數據的變尺度加密三角剖分算法，此算法在傳統的三角剖分算法中進行了三個方面的改進。

(1) 該算法在對傳統的凸殼生成方案上進行了改良，在其基礎上提出了“邊界收縮”，即形成高精細度且更貼合實際地層的凸包，優化了點定位算法。

(2) 該算法在傳統的加密算法上進行了改良，設定一個最長邊，進行加密后的三角網更加規整，方便貼合實際地形進行處理。

(3) 該算法優化了約束線段影響域的處理，使得含約束線段的三角剖分在時間和效率上優于傳統算法，基于本算法在三維中進行拓展，為實際建模提供了很好的方向，有較好的實用價值。