史曉楠 王夢凡 王子童

(西安科技大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院 陜西 西安 710054)

0 引 言

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

本文主要研究的是包含斷層的數(shù)據(jù)如何較高效率且高質(zhì)量的實(shí)現(xiàn)Delaunay三角剖分的問題?；趦刹椒?，本文首先使用逐點(diǎn)插值方法和可變尺度加密算法生成無約束DT(Delaunay Triangulation)網(wǎng)格，然后先用細(xì)分算法細(xì)分?jǐn)鄬舆?，避免了狹長三角形的出現(xiàn)，再結(jié)合生長法和分治法思想構(gòu)建帶斷層約束的CDT(Constrained Delaunay Triangulation)網(wǎng)格。對比其他常見三角剖分方法，本文方法可以較好地解決斷層線段的難以合理的插入的問題，有效提高了構(gòu)網(wǎng)質(zhì)量和效率。

1 算法整體流程

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

圖1 算法流程圖

2 關(guān)鍵算法設(shè)計(jì)

2.1 無約束三角網(wǎng)的形成

(1)

(2)

(3)

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

(4)

將所有離散點(diǎn)按照上述規(guī)則投影到網(wǎng)格，如圖2所示。

圖2 虛擬網(wǎng)格

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

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

圖3 凸殼初步建立

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

查詢出線段右側(cè)距離最大的點(diǎn)，可以結(jié)合網(wǎng)格和三角形面積判斷。已知三點(diǎn)坐標(biāo)，由面積公式得面積最大對應(yīng)的點(diǎn)(x3,y3)：

(5)

Step3當(dāng)鏈表中頂點(diǎn)個數(shù)等于每次沒有插入點(diǎn)的次數(shù)，停止操作，初始凸包成功建立。

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

Step1設(shè)置最大凸殼邊界約束值。

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

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

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

圖4 “高精細(xì)度”凸殼

(4) 初始三角剖分。任取凸殼內(nèi)一點(diǎn)與凸殼點(diǎn)進(jìn)行連接，構(gòu)成三角網(wǎng)的雛形。隨機(jī)選取凸殼內(nèi)任意一點(diǎn)，依次連接該點(diǎn)和圖中所有的點(diǎn)進(jìn)行圖5所示的初始三角剖分。

圖5 初始三角剖分

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

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

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

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

(6)

反之，若出現(xiàn)如下所示的不等關(guān)系，表示點(diǎn)在外部。

abs(S1)

(7)

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

圖6 面積法判斷點(diǎn)的位置

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

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

Step4LOP算法優(yōu)化[11]。

Step5重復(fù)進(jìn)行Step 1，當(dāng)所有點(diǎn)都完成插入時，三角網(wǎng)構(gòu)建完成。

圖8是逐點(diǎn)插入法的流程框架。

圖8 逐點(diǎn)插入法框架圖

2.2 變尺度加密算法

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

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

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

Step2最初用H×N(N為某個自然數(shù))進(jìn)行全局細(xì)化。

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

(8)

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

變尺度加密前后的三角網(wǎng)如圖9所示。

圖9 變尺度加密前后三角網(wǎng)

2.3 斷層數(shù)據(jù)的嵌入

(1) 斷層數(shù)據(jù)細(xì)分算法。連接所有斷層數(shù)據(jù)，形成斷層邊，細(xì)分所有斷層邊，將形成的點(diǎn)存入數(shù)組K[14]。

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

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

(2) 插入細(xì)分后斷層數(shù)據(jù)點(diǎn)。細(xì)分后，將獲得的點(diǎn)插入網(wǎng)格。

Step1根據(jù)變尺度算法中的最大高度H來設(shè)置搜索半徑r。

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

Step3將ki與E中三角形外接圓進(jìn)行判斷，不符合則標(biāo)志此三角形。

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

(a) 插入新結(jié)點(diǎn)P (b) 決定P如何和其他結(jié)點(diǎn)相連

(c) 刪除AB邊 (d) 形成新三角形圖10 插入斷層數(shù)據(jù)點(diǎn)

Step 3中，點(diǎn)定位算法具體步驟如下：

假設(shè)P在BC的右側(cè)，∠A+∠P≥π，則P在圓內(nèi)，∠A+∠P<π，則P在圓外。

點(diǎn)P在圓內(nèi)，根據(jù)余弦公式可轉(zhuǎn)化為:

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

(9)

點(diǎn)P在圓外，則:

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

(10)

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

(a) (b) (c)圖11 點(diǎn)與外接圓位置關(guān)系

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

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

Step3用一個隊(duì)列Q存儲數(shù)組K中的點(diǎn)，在影響域內(nèi)，先以kiki+1為基邊(以右邊為例)，向右擴(kuò)展，搜索P，使得角Pkiki+1最大，即構(gòu)成了Pkiki+1。

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

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

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

3 算法分析

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

綜上，本文無論是在逐點(diǎn)插入法、點(diǎn)定位算法還是約束區(qū)域的嵌入上，算法的時間復(fù)雜度均較低。

4 實(shí) 驗(yàn)

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

表2 算法執(zhí)行效率測試結(jié)果

圖13為含斷層數(shù)據(jù)的變尺度加密三角剖分效果圖?？梢钥闯觯?a)是未經(jīng)過變尺度加密的三角網(wǎng)，不含斷層數(shù)據(jù)；(b)是經(jīng)過變尺度加密后的三角剖分結(jié)果，斷層細(xì)節(jié)比較簡單；(c)是在(b)的基礎(chǔ)上進(jìn)一步進(jìn)行變尺度加密的三角網(wǎng)，斷層顯示較為明顯。同時從表2可以看出，本文使用的算法在三角剖分的時間效率上優(yōu)于傳統(tǒng)算法。

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

5 結(jié) 語

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

(1) 該算法在對傳統(tǒng)的凸殼生成方案上進(jìn)行了改良，在其基礎(chǔ)上提出了“邊界收縮”，即形成高精細(xì)度且更貼合實(shí)際地層的凸包，優(yōu)化了點(diǎn)定位算法。

(2) 該算法在傳統(tǒng)的加密算法上進(jìn)行了改良，設(shè)定一個最長邊，進(jìn)行加密后的三角網(wǎng)更加規(guī)整，方便貼合實(shí)際地形進(jìn)行處理。

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