一種基于邊指針搜索及區域劃分的三角剖分算法

張俊 田慧敏

三維重建技術是根據物體的二維圖像信息恢復其三維模型,多年來一直是計算機圖形學領域的研究熱點[1?2].三角剖分是三維重建過程中最重要的基礎之一,通過將散點連接形成許多互不相交的三角形或四面體,從而使模型面化或者體化[3].經過多年探索,研究者們對于二維三角剖分的研究已經取得了很多成果,提出了多種三角網格的優化準則和構造方法.其中,Delaunay 三角剖分由于具有良好的幾何特性-可以最大化最小角[4],即能使得到的三角網格最平均,使用最廣泛[5].然而,隨著計算機圖形學技術的發展,掃描設備的精度越來越高,用于三維重建的點集對象的規模也越來越大.因此,在保證三角網格的整體質量的提前下,有必要尋找一種更高效的Delaunay 三角剖分,從而滿足三維重建的實時性.

Delaunay 三角剖分算法有很多種類,主要分為三大類:逐點插入算法[6]、分治算法[7]、三角網生長算法.其中,逐點插入法由于實現相對簡單且占用空間相對較小[8],被大量研究者所采用;分治算法在執行時間方面可以達到最優,但占用內存較大,不適用于一般的計算機平臺;三角網生長算法相對于前兩種算法效率較低,不適用于大規模點集,在實際應用中很少被采用.逐點插入法的主要流程如下:1) 構造一個包含所有插入點的超級三角形;2) 將存儲在鏈表中的點按序插入,遍歷三角形鏈表找出包含插入點的三角形(即目標三角形),利用插入點將目標三角形拆分生成新的三角形,在三角形鏈表中刪除目標三角形,完成一個點的插入;3) 利用空圓特性對新生成的三角形進行Delaunay 規則判斷;4) 循環執行2) 和3),直至完成對所有插入點的處理.經典的逐點插入法有Bowyer 算法[9]、Watson 算法[10]以及Lawson 算法[11]等.

雖然逐點插入法易于實現,但在處理大規模的散點集數據時,實時性不是很好,為了提高逐點插入法的效率,研究者們進行了大量的研究.一部分研究者認為可以對插入點進行預處理,通過改進插入點的序列來提高三角剖分的效率.Liu 等[12]提出了一種基于廣度優先搜索的確定性插入序列,使用k-d樹來構造Delaunay 三角關系,但是構造k-d 較為復雜,需耗費大量時間.在多重網格插入法[13]的基礎上,Su 等[14]提出了使用Hilbert 曲線遍歷插入點,減少了三角剖分過程中狹長三角形出現的數量,從而降低局部優化時間.Zalik 等[15]提出一種兩級均勻細分加速技術,將目標三角形問題轉化為最近點問題,在他們的方案中,從最近的Delaunay 點開始,然后從第二個最近的Delaunay 點開始,以這種遞歸的方式直到找到目標三角形.在此基礎上,Zadravec等[16]提出結合哈希表與跳躍表來尋找插入點的最近點,從而快速定位目標三角形.一些研究人員提出基于重心方向來定位目標三角形,但是,當出現分界點時,這種方法的搜索路徑可能不唯一.針對此問題,Xi 等[17]提出可以通過移動重心來避免截點,從而使目標三角形可以連續搜索.隨著計算機技術的快速發展,許多研究者們提出了并行Delaunay 三角剖分,將點集劃分成許多獨立的分區,這些分區可以同時進行Delaunay 三角剖分,Kohout 等[18]、Rong等[19]、Cuong 等[20]、Lo[21]均對此進行了研究及改進.此外,楊昊禹等[22]提出了并行動態Delaunay三角剖分算法,可以解決新增點位于原來三角網之外的情況.另外李國慶等[23]提出了基于凸多邊形的Delaunay 三角剖分算法,使用生成的凸多邊形代替傳統算法中的超級三角形.常見的生成凸包算法有Graham 掃描法[24]、分區算法[25]等.劉斌等[26]提出將主成分分析法(PCA) 與二分法結合,通過快速確定凸包邊緣點來計算平面點集的凸包.

通過對Delaunay 三角剖分逐點插入法的研究,可知在構建三角網過程中最耗時的部分為尋找插入點的目標三角形.本文在傳統算法的基礎上,將邊指針與區域劃分相結合,很大程度上提高了構建Delaunay 三角網的效率.

1 基于邊指針的搜索

在傳統的逐點插入法中,對每個插入點的目標三角形的搜索總是從三角形鏈表的表頭開始.這意味著,目標三角形若位于鏈表的表頭,則可以很快被找到;若位于鏈表的尾部,則需要遍歷整個鏈表才能被找到.隨著點的數量大幅度增加,三角形鏈表的長度也必然會大幅度增加,從而導致尋找目標三角形會越來越耗時,降低三角剖分的效率.針對此問題,本文對三角形的數據結構進行優化,提出了基于邊指針的數據結構.如圖1 所示,在△ABC中,A為起始頂點且A、B、C逆時針分布,AC為邊1,AB為邊2,BC為邊3.在三條邊上分別存在邊指針P1、P2、P3,分別指向邊指針所在的公共邊的相鄰三角形,如P1 指向與邊AC相鄰的三角形,對于P2、P3 情況類似.

圖1 三角形數據結構Fig.1 Data structure of triangle

1.1 點與三角形的位置判斷

點與三角形之間的位置關系可以根據其坐標的相對關系來判斷.對于任意三個點A(xA,yA),B(xB,yB),C(xC,yC),其中

若Flag(A,B,C)>0,則A、B、C按逆時針方向分布;若Flag(A,B,C)<0,則A、B、C按順時針方向分布;若Flag(A,B,C)=0,則A、B、C三點共線.如圖2(a) 所示,若Flag(C,A,M)、Flag(A,B,M)、Flag(B,C,M)三者的值都大于零,則點M位于△ABC的內部.若三者中一個值為零且其余二者皆為正值,則位于某條邊上,如圖2(b) 中Flag(A,B,M) 為零,Flag(C,A,M)、Flag(B,C,M)為正值,故點M位于邊AB上.若三者中至少有一個為負值,則點位于三角形的外部.如圖2(c) 所示,Flag(A,B,M)為負值,故點M位于邊AC的外側.當然,點位于三角形外側還會出現圖2(d) 與圖2(e) 的情況.在圖2(d) 中,Flag(A,B,M)、Flag(B,C,M) 皆為負值,在此情況下,由于AB為邊2,而BC為邊3,本文中優先取編號在前的邊,故認為M位于邊AB的外側.在圖2(e) 中,點M雖位于三角形外部且與邊AC共線,但是Flag(B,C,M) 為負值,故本文認為點M位于邊BC的外側.

1.2 基于邊指針搜索的路徑確定

圖2 點與三角形位置關系判斷Fig.2 Judgment of positional relationship between point and triangle

在搜索目標三角形的過程中,搜索路徑通過判斷插入點與當前被搜索三角形的位置關系來確定的.如圖3 所示,點V為插入點,假設△AIJ存儲在三角形鏈表的表頭,搜索過程如下所述.1) 檢測點V與△AIJ的位置關系,首先確定了點V位于△AIJ的外部,然后確定點V位于△AIJ的哪條邊的外側,根據式(1) 可知點V位于邊IJ外側,故下一個被搜索的三角形為△AIJ的邊IJ上的邊指針p1所指的△IEJ.2) 用同樣的方法確定接下來的被搜索三角形分別是△EFJ、△EHF,然后根據式(1)可知點V位于△EHF的內部,故△EHF為點V的目標三角形.在整個搜索過程中,共經過4 次搜索找到點V的目標三角形,搜索路徑為:△AIJ(p1)→△IEJ(p2) →△EFJ(p3) →△EHF.邊指針搜索的方式,由于在搜索目標三角形的過程中通過邊指針的方向性可以快速地定位目標三角形,使搜索路徑大大縮短,從而大大降低了三角剖分的耗時.

圖3 點V 的目標三角形的搜索路徑Fig.3 Search path for target triangle of point V

2 基于邊指針搜索及區域劃分的三角剖分

2.1 區域分塊

基于邊指針搜索目標三角形的方法,雖然可以利用插入點與三角形的位置關系來快速定位目標三角形,但是在點數量很大的情況下,目標三角形的搜索路徑可能仍然會比較長,此時依舊較為耗時.因此,在邊指針搜索的基礎上,提出了區域劃分的方案,以縮小目標三角形的搜索范圍.

設區域分塊為N × N,構造超級三角形的時候,首先遍歷所有散點找到最大的縱橫坐標值Xmax、Ymax,取兩者中較大者得到max,然后根據max構造超級三角形,其各點坐標為:

區域劃分是指將包含超級三角形的正方形劃分成N ×N個大小相同的區域,并將這些區域存儲在二維數組中.每個區域的位置坐標值X和Y對應于其在二維數組中的存儲位置,如位置坐標(1,2) 的區域塊存儲在數組area[1][2]中.對于插入點P(x,y),它的區域坐標由其本身的橫縱坐標值確定:

在式(3) 中,L指包含超級三角形的正方形的邊長.

區域可以劃分成不同的大小,對應的區域數量也會不同.對于某個確定的點集,如果區域太大,搜索范圍會變大.但是,如果區域太小,則過多的區域很難被分配到入口三角形,從而導致入口三角形的使用率下降,進一步導致搜索目標三角形的時間增加.所以對于處理不同規模的點集,最優的分塊也會不同.

對于點集V={V1,V2,···,Vn},兩個點組成一條邊,三條邊組成一個三角形,三角形通過公共邊上的邊指針相互連接組成三角形網格.基于邊指針搜索及區域劃分的Delaunay 三角剖分的數據結構如表1 所示.

表1 基于邊指針及區域劃分的算法數據結構Table 1 Data structure of algorithm based on edge-pointer and region-division

數據結構包括4 部分:點、邊、三角形和區域.每個頂點包含x和y坐標以及區域入口標志entryflag;連接頂點的邊包含頂點Vi和Vj的編號;三角形采用空間網狀結構,包含頂點Vi、Vj、Vk的編號,三角形的區域入口標志entry以及該區域的坐標X、Y,將位置相鄰的三角形連接起來的邊指針p1、p2、p3,區域坐標只有在entry=1 時才意義,如果entry=1,意味著三角形是該區域的入口三角形;區域塊包含區域是否分配入口三角形的標志valid,以及區域的入口三角形指針對于一個入口三角形,在其三個頂點中存在一個入口頂點,其值被置為1,并且入口頂點位于于該三角形所記錄的區域中.邊指針搜索及區域劃分數據結構中不同區域之間的關系如圖4 所示,不同局域會分配一個入口三角形作為該區域的搜索入口,該入口三角形再通過邊指針指向存在該區域內或其它區域內的相鄰三角形.

圖4 結構體關系圖Fig.4 Structure diagram

2.2 搜索目標三角形

實現Delaunay 三角剖分最重要也最耗時的部分就是搜索目標三角形,本文對目標三角形的搜索過程如算法1 所示.在基于邊指針搜索的基礎上,區域劃分的搜索方式為每個區域分配了入口三角形.這使得對散點的目標三角形的搜索不是從三角形鏈表的表頭開始,而是從入口三角形開始,故搜索的起始三角形更接近目標三角形,目標三角形的搜索路徑進一步被簡化.

圖5 展示了基于邊指針搜索及區域劃分算法中搜索目標三角形的過程.如圖5(a) 所示,包含超級三角形的正方形被劃分成2×2 個區域.插入點V時,首先根據點V的坐標確定其所在區域為area[1][0].由于在點V之前,點H與點K已經被插入該區域,故該區域一定存在入口三角形,假設△GCH為入口三角形,對點V的目標三角形的搜索從△GCH開始.接下來,按照第1.2 節中所述邊指針搜索的方式,確定目標三角形的搜索路徑為:△GCH(p1) →△DGH(p2) →△EDH(p3) →△EHF,目標三角形為△EHF.

如果檢測到點V位于三角形的某條邊上,則需要根據邊指針搜索另一個包含該點的一個三角形.如圖5(b) 所示,點V位于邊HF上,當找到△EHF之后,根據邊HF上的邊指針p4找到共邊三角形△FHK.故△EHF與△FHK為點V的目標三角形,并且點V 將會被插入這兩個三角形所構成的四邊形內.

如果在點V所在區域中沒有入口三角形,則從最新插入的點生成的三角形開始搜索目標三角形,搜索路徑由邊指針確定.由于該算法主要用于三維重建,而三維重建中掃描得到的連續點在位置上可能是相鄰的,因此,由最新插入點生成的三角形可能更接近當前插入點的目標三角形.

在傳統算法中,假設當前所有三角形按一下順序存儲在三角形鏈表中:△AIJ、△AJF、△IEJ、△EFJ、△IBL、△ILD、△IDE、△EDH、△EHF、△FHK、△CFK、△BGL、△DLG、△DGH、△GCH.插入散點V時,則需要依次判斷鏈表中在△EHF之前的三角形是否包含點V,可以發現一共需要經過9 次搜索才能找到目標三角形,而本文的算法只需要經過4 次搜索就能找到目標三角形,搜索深度為傳統算法的44.4%.隨著點集的規模不斷變大,傳統算法的執行時間將會呈指數增長,本文算法的優越性會越來越明顯.

圖5 搜索點V 的目標三角形Fig.5 Searching for target triangle of point V

目標三角形的定位是根據插入點與當前搜索三角形的位置關系來確定的,通過判斷插入點位于三角形的哪條邊的外側來確定下一個被搜索的三角形.每一次被搜索的三角形都比上一次的三角形更接近插入點,因此通過這種方式一定能找到目標三角形,故該算法是收斂的.

2.3 區域入口三角形及邊指針的更新

隨著點地不斷插入,區域的入口三角形也需要不斷更新以確保其不會隨著某些三角形的拆分而消失.同時,與新三角形有關的邊指針也需要被更新以確保三角形之間能夠互相找到相鄰的三角形.找到目標三角形之后,需要將點插入到目標三角形中,并使用空圓特性對新生成的三角形進行Delaunay 規則判斷及調整違背規則的三角形.在此過程中,需要更新入口三角形及邊指針,入口三角形的更新分為本區域的更新以及他區域的更新.

如圖6 所示,找到目標三角形△EHF之后,插入點V,生成了三個新的三角形(△EHV,△HFV與△FEV),并將原目標三角形△EHF從三角形鏈表中刪掉.首先,對本區域的入口三角形進行更新.如果在此之前,區域area[1][0]沒有分配到入口三角形,則選擇△EHV為該區域的入口;如果area[1][0]之前已經分配了入口三角形如第2.2 節中所提到的△GCH,則△EHV代替△GCH作為該區域的新入口,并將點V指定為入口頂點.然后,根據被刪除三角形的標志位及區域坐標信息檢測其是否為其他區域的入口,若是,則需要指定一個新的三角形來更新此區域的入口.例如,如果三角形△EHF是區域area[1][1]的入口,判斷△EHF的哪個頂點是入口頂點,根據點的入口標志位得出點F是入口頂點,因此,包含點F的△FEV成為了area[1][1]的新入口.

圖6 入口三角形及邊指針的更新Fig.6 Update of entry triangle and edge-pointer

此外,還需要更新和原三角形△EHF有關的邊指針.在圖6 中,隨著△EHF的刪除,原本指向該三角形的邊指針應該指向新的三角形,△EDH的邊EH上的邊指針p3現在應該指向△EHV,同理,△EFJ的邊EF上的邊指針p5、△FHK的邊FH上的邊指針p6現在分別指向其相鄰三角形△FEV、△HFV.

在進行Delaunay 規則判斷的過程中,區域入口及邊指針會以同樣的方式更新,此處不再贅述.

3 實驗結果及分析

為了驗證基于邊指針搜索及區域劃分的三角剖分算法的效率,本文對隨機生成的不同規模大小的點集進行三角剖分,并比較了傳統的Delaunay 三角剖分(TD)、使用CGAL 庫的三角剖分、基于邊指針搜索的三角剖分(EPD)、基于邊指針搜索和區域劃分相結合的三角剖分(EPRDD) 幾種算法之間的執行時間,此處的執行時間不包括讀取插入點數據的時間.

在進行區域劃分的時候,包含超級三角形的正方形被劃分成了15×15 個區域.本實驗的實驗環境為Intel(R) Core(TM) CPU i5-8500 @ 3.00 GHz,16 GB 內存,Windows 10,VS2013,64 位操作系統,得到的結果如表2 所示:

表2 算法的執行時間(s)Table 2 Running time of algorithms (s)

從表2 可以看出,當插入點的數量為2 萬時,基于邊指針搜索的三角剖分耗時為傳統算法的1/201,加入區域劃分之后,算法的執行時間為傳統算法1/970.當插入點的數量為10 萬時,基于邊指針的搜索算法執行時間分別為傳統算法的1/643,而基于邊指針及區域劃分算法的執行時間為傳統算法的1/4 886.此外,與CGAL 算法相比,只使用邊指針的算法優勢并不是很明顯;但是,邊指針與區域劃分相結合之后,算法的耗時明顯低于CGAL 算法.當對10 萬個點進行三角剖分的時候,CGAL 算法需要8.594 s,而基于邊指針及區域劃分的算法僅僅需要0.813 s 就能完成.

為了比較幾種算法的執行時間與散點數量的關系,繪制了圖7 所示折線圖.由于傳統算法的執行時間遠遠高于其余幾種算法,故不與其他算法作比較.從圖7 可以看出,隨著散點數量的增多,CGAL算法,基于邊指針搜索的算法與基于邊指針搜索及區域劃分算法的執行時間基本上都是呈線性增長的,CGAL 算法與基于邊指針搜索的算法的執行時間隨散點數量變化的增長速度還是較快,但是基于邊指針及區域劃分算法的執行時間的增長速度較為緩慢,斜率遠遠小于前兩種算法.可以得知,隨著插入點數量的增長,邊指針搜索與區域劃分結合的方法在時間上的優勢將會越來越大.

圖7 執行時間比較Fig.7 Comparison of run time

進一步對幾種算法的平均搜索深度進行比較,結果如表3 所示.搜索深度是指在將散點插入到Delaunay 網格中的過程中搜索三角形的次數,搜索深度越小,三角剖分的實時性越好.由于CGAL 的源代碼沒有找到,在此處未對CGAL 的搜索深度進行比較.從表中可以看出,在傳統算法中平均搜索深度與散點數目也是線性增長的,而基于邊指針搜索及區域劃分的算法的平均搜索深度很小.即使散點數量成倍增加,平均搜索深度也只是稍微增加,證明邊指針與區域劃分結合的算法實時性要高于傳統Delaunay 三角剖分算法.其原因在于,區域劃分縮小了目標三角形的搜索范圍,邊指針優化了搜索方向,最終使得搜索深度大大減小,執行時間得到降低.

表3 算法的平均搜索深度Table 3 Average search depth of algorithms

4 結論

為了滿足三維重建的實時性,針對其基礎部分三角剖分處理時間過長的問題,本文提出了一種基于邊指針搜索及區域劃分的Delaunay 三角剖分算法.在對傳統Delaunay 三角剖分算法進行研究及實現的基礎上,進一步優化了三角形存儲的數據結構及搜索過程,設計了一種通過邊指針反映三角形空間相鄰關系的三角形數據結構來取代傳統算法中的一維三角形鏈表.由于邊指針具有方向性,可以用于確定距離插入點最近的相鄰三角形,從而快速確定目標三角形的搜索方向.在此基礎上將整個超級三角形進行區域劃分,且隨著散點的插入,每個區域將會分配一個區域入口三角形,作為搜索目標三角形的起始三角形,且入口三角形在構建Delaunay 三角網的過程中會不斷更新.這使得目標三角形的搜索范圍由以前的整個三角形鏈表縮小到插入點所在區域,從而進一步縮短搜索的起始三角形與目標三角形之間的距離,大大縮小了搜索范圍,降低了搜索目標三角形的耗時,以滿足三維重建的實時性需求.實驗結果表明,該算法的執行時間隨散點數量的增長較為緩慢,且散點數量為10 萬時,所耗費的執行時間僅為0.813 s.

本文中的區域劃分為固定劃分,當散點分布密度變化大時,對于密度大的區域搜索深度可能會相應增加.為了進一步擴大算法的適用范圍,后續研究中將引入自適應算法,尋求一種能夠根據散點分布密度或者搜索深度的變化對區域劃分方式進行調整的優化算法.