Delaunay三角剖分算法改進與對比分析

袁小翠　吳祿慎　陳華偉

(南昌大學機電工程學院　江西 南昌 330031)

?

Delaunay三角剖分算法改進與對比分析

袁小翠吳祿慎陳華偉

(南昌大學機電工程學院江西 南昌 330031)

針對Delaunay算法的計算速度問題，從數據結構和算法兩個方面加以改進。對Delaunay三角剖分的代數拓撲分析，設計一種順序存貯的Hash數據結構，實現臨時單純形對象的快速和順序存取、查詢、插入和刪除等操作；以單純形邊對象的活性分析為核心，以Hash數據結構進行操作，消去生長法的遞歸過程；此外，提出基于微切平面的生長法，將基于空間四面體的空球搜索降維至局部二維的空圓搜索。對汽車擋泥板和兔子模型進行三角剖分實驗，實驗結果表明，消去遞歸的生長法和基于微切平面的生長法和傳統(tǒng)的生長法三角剖分效果相同，但是計算速度比傳統(tǒng)方法效率更高。

Delaunay三角剖分生長法半空間隱式曲面

0　引　言

逆向工程通過三維光學掃描儀獲取物體形貌數據，該數據通常是海量的三維點，在數據處理中，需要經過三角化、光滑曲面重建才能為正向CAD軟件所采用。點云三角化網格重建是數據由離散點域向面域求解的第一步，是后續(xù)光滑曲面重構的基礎。其重構方法主要有基于Delaunay三角剖分的方法[1,2]、基于區(qū)域擴張的方法[3,4]和基于實體布爾差運算的方法[5,6]，其中生長法是最常用的Delaunay三角剖分法[7]。目前，已有諸多學者開始采用隱式曲面法[9-11]構造三角剖分。盡管如此，作為三角剖分的基礎，Delaunay法仍然集成于主流算法和商品化軟件之中，尤其是在國內，Delaunay三角網格構造仍然是一個熱點研究問題。作為一個三角網格剖分的基礎和共性問題，本文對Delaunay三角剖分中的生長法進行分析，針對Delaunay算法的計算速度問題，設計一種順序存貯的Hash數據結構，實現臨時單純形對象的快速存取、插入和刪除等操作。此外，本文還提出基于微切平面的生長法，將基于空間四面體的空球搜索降維至局部二維的空圓搜索，提高三角剖分的時間效率。

1　Delaunay三角剖分

代數拓撲中，n-單純形(n-simplex)拓撲等價于n-球面，每個n-單純形是n維帶邊界流形。n-單純形是三角形和四面體的一種泛化，一個n維單純形是指包含n+1個節(jié)點的凸多面體。其精確定義是：n-simplex是某個n維以上的歐氏空間Rd(d≥n)中的n+1個仿射無關的點集的凸包。例如，0-simplex就是點，1-simplex就是線段，2-simplex就是三角形，3-simplex就是四面體。將構成單純形的各點稱為單純形的頂點(vertex)。進一步地，將單純形的k+1個頂點的凸包組合所構成的單純形稱為k維面(k-face)，如點是0-face，邊是1-face，三角面是2-face。

拓撲空間可以通過單純形組合構造，人們希望把一個拓撲對象剖分成許多個小的單純形，并要求任何兩個相鄰單純形的相交部分仍是一個單純形——這種剖分稱為單純剖分。在曲面情形，就是三角剖分。

R3空間中面的剖分是2-simplex，即三角形的組合，則相鄰單純形的相交部分為1-face，即為邊。三角化過程就是以邊為基礎，逐步向外擴展生成三角形，新三角形的邊將作為繼續(xù)生長的基礎。三角形生長應遵循一定的規(guī)則，Delaunay三角剖分規(guī)定了兩項重要規(guī)則：

1) 空圓特性：任一Delaunay三角形的外接圓內不會有其他點存在，如圖1(a)所示；

2) 最大化最小角特性：在點集P可能形成的三角剖分中，Delaunay剖分形成的三角形的最小角最大，如圖1(b)所示。

圖1　Delaunay三角剖分中的兩個重要特性

2　生長法

圖2　生長法三角剖分流程圖

生長法[12]是Delaunay三角剖分最直接的方法，算法設計直接由Delaunay三角剖分定義而來，算法簡單明了。根據拓撲定義，Rn空間可由n維(及n維以下)單純形組合構造，相鄰兩個單純形的相交部分為n-1維面。生長法是以面為基礎，在生長規(guī)則的控制下，為面加入頂點，構造新的單純形。生長法三角剖分流程如圖2所示。

傳統(tǒng)的生長法[13，14]采用遞歸過程實現，Triangulation算法為遞歸生長法的C語言偽代碼：

FunctionTriangulation(P,f0)/*P為輸入點集，f0為輸入單純形面*/{

v = FindNextVertex(P,f0);

s = BuildSimplex(f0,v);

print(s);

//打印輸出

for each f∈s

Triangulation(P,f);

}

Triangulation算法以一輸入面f0為基礎，使用FindNextVertex函數查找一個頂點v，使用f0和v構造單純形s。接下來遍歷s中的每個面f，重新調用Triangulation構造下一個單純形。print語句打印輸出新構建的單純形，如果需要存貯輸出，則只需設置一個列表，順序存貯即可。

3　改進生長法

Triangulation算法是一個遞歸算法，容易理解和實現，但是大數據量計算的算法時域和空域特性差，因此需進行消去遞歸。據此，本文從數據結構上改進遞歸生長法，設計了一個(n-1)維面列表FList，用于動態(tài)存貯各面，從而消除遞歸。

3.1數據結構設計

數據結構設計上，隊列能夠實現順序存取和快速查找，但是刪除操作效率低。為了提高算法性能，尤其是實現隨機位的快速搜索和刪除操作，采用哈希表進行數據操作。常規(guī)哈希表中數據的存貯位是隨機的，因此還必須設計順序存貯的哈希表結構。本文使用C++標準庫的list容器作為順序存貯器，結合Boost庫的非排序類容器unordered中的unordered_set哈希模板，構造順序存貯的哈希表。

unordered_set模板主要有三個參數：T1、T2和T3分別代表存貯對象、Hash編碼函數(hash_value)和等值判斷函數(equal)。其中hash_value函數為一個存貯對象生成哈希鍵值(key)，當兩個對象有相同的key值而產生訪問沖突時，就會調用equal函數對兩個對象作進一步比較。為了提高程序的可讀性和可維護性，通常在對象類中定義hash_value和equal函數，然后以函數指針方式傳入Hash表中進行調用。但unordered_set不能直接以函數指針為模板參數，因此將函數指針封裝為結構體，以結構體參數形式傳入。

數據結構定義：

template

struct StruHashFunc

{

size_t operator()(Tp1) const

{

return pHashFunc(p1);

//hash_value函數以一個對象為參數

}

};

template

struct StruEqualFunc

{

bool operator()(Tp1,Tp2) const

{

return pEqualFunc(p1,p2);

//equal函數以兩個對象為參數

}

};

typedef StruHashFuncTHashFunc;

typedef StruEqualFuncTEqualFunc;

接下來，就可以自定義順序哈希表。為了便于用戶使用，以list為外部容器，用于輸入輸出。在數據結構內部使用hash容器存貯對象并生成對應key值，用于搜索定位，兩者內部的讀寫操作要保證一致。如下為自定義接口IListedHash的過程：

1) 以對象類為模板參數定義該接口類：

ctemplateclass IListedHash;

2) 以TValue為list參數，定義list接口類：

typedef std::listIList;

3) 定義哈希表接口類：

typedefboost::unordered_set;

IHash;

4) 定義迭代器，其中hash迭代器為內部迭代器，list迭代器為外部迭代器：

typedef typenameI Hash::iteratoriterator_in;

typedef typename IList::iteratoriterator;

5) 定義內部數據存貯變量：

IListlistdata和IHash hashdata;

6) 定義與list類似的各項操作，主要有insert、find、erase等操作和begin、end等迭代器操作方法，在實現上保證listdata和hashdata數據的同步讀寫。

采用以上設計的哈希表進行數據操作，對Triangulation算法進行消去遞歸，改進后為Triangulation2算法。

Triangulation2算法的偽代碼為：

FunctionTriangulation2(P,FList)

/*P為輸入點集，FList為輸入面列表*/

{

whileFList≠

{

f0 = ExtractHash(FList);

v = FindNextVertex(P,f0);

s = BuildSimplex(f0,v);print(s);

//打印輸出

for each f∈s

UpdateHash(FList,f);

}

}

UpdateHash(FList,f)

{

if f∈FList

DeleteHash(FList,f);

else

InsertHash(FList,f);

}

上述代碼表示中有三個哈希表的操作函數：ExtactHash、DeleteHash和InsertHash，它們分別對面列表執(zhí)行讀取、刪除和插入操作，自定義UpdateHash函數對面進行更新操作。為了避免重復生長，讀取f0的ExtractHash函數中包含了f0的DeleteHash操作。算法以面列表FList為輸入，循環(huán)判斷是否為空，每次循環(huán)從FList中讀取一個面f0，構建單純形s，遍歷s中的每個面f，如果已存在于列表，則從FList中刪除，否則，插入f0，算法直至FList為空時結束。

3.2三角剖分的詳細過程

三角化算法的核心語句是FindNextVertex函數中滿足Delaunay規(guī)則的點的搜索。為了敘述方便，本文給出如下兩個定義：

定義1由多個三角形或者兩條邊界邊完全封閉的點為死點；反之，未被完全封閉的點稱為活點。

定義2相鄰三角形的共有邊或邊界邊為死邊；其他邊稱為活邊。

設點集為P，當前遍歷面為邊e(v0,v1)，其中v0，v1∈P是e的兩端點，以邊e為基礎，采用以下生長規(guī)則搜索第三點v2：

?只對鄰域點進行搜索

此規(guī)則用于確定總體搜索范圍。設點v0和v1的鄰域分別為Neib0和Neib1，則搜索范圍確定為Neib01=Neib0∪Neib1。

?只對活點進行搜索

圖3　點和邊的活性

在圖3中，粗線邊v5v7和邊v7v8為已確定的邊界邊。以v0為端點的邊，以及邊v5v4、v4v7、v4v8均為死邊，兩條邊界邊v5v7、v7v8亦為死邊，其他邊均為活邊。顯然，死邊是當前搜索狀態(tài)下的內部邊，不能作為基礎邊再行生長，需要從面列表FList中刪除；相對應地，活邊作為外部邊，應作為生長邊進入FList。

為了準確表述，用點或邊的使用頻次f表示點或邊的活性。在構網過程中，一條邊只能為邊界邊，或者為相鄰兩個三角形的共有邊，因此邊的活性取值范圍為{0,1,2}。點的使用頻次與其構邊的次數有關，設與點v關聯(lián)的邊集為E∈FList，則f的計數規(guī)則是：如果E=?，則置f=-1；否則f=E.size。由于E是FList的子集，因此f值和E的增減都在FList中體現：v的每次構邊都將通過InsertHash加入FList，相應有f++；每次提取邊構建三角形時都將使用ExtractHash從FList中刪除邊，相應有f--。

?生長點必須位于近鄰三角形的右側

此規(guī)則用于限定生長方向。為了定義三角形的生長方向，一般統(tǒng)一定義為逆時針(或順時針)方向為三角面正向。這也延伸定義了構面邊為有向邊，其方向與面的方向一致。如圖4所示，已構三角形面v0v1v2為逆時針方向，三條邊v0v1、v1v2、v2v0亦為逆時針方向。此時，若對邊v1v2搜索新的生長點，為了避免出現交叉三角形，必須限定生長點位于邊v1v2的“右側”。二維情形下，邊v1v2將平面劃分為左右兩個半平面；三維空間中，假想經過邊v1v2創(chuàng)建一個與三角面v0v1v2垂直的平面，該平面將整個空間劃分為兩個半空間，分別記為HS-和HS+。邊v1v2是兩個半空間的一條交界邊，為了滿足三角面方向統(tǒng)一的要求，在搜索空間HS+中，對邊v1v2進行擴展時，應取其反向邊v2v1。本規(guī)則限定了生長點必須位于HS+空間，該規(guī)則下，圖3的備選點集{v3,v4,v5}中只有{v4,v5}為合法生長點集。

圖4　半空間搜索

?限定三角形的最小角度Amin

此規(guī)則避免出現狹長三角形，用于滿足Delaunay最大化最小角特性。推薦Amin=5°～30°。

?三角形外接圓半徑與當前邊的比值大于限定值Rmin

此規(guī)則用于外接圓半徑，即縮小待建三角形的影響范圍，旨在滿足Delaunay空圓特性。

3.3基于微切平面的生長法

在三維點云空間中，對鄰域點的搜索實際是基于四面體的空球搜索，搜索速度和精度以及編程難度比二維點云大。對此，考慮將三維空間降維至二維，將三維點投影至微切平面，在二維平面內構建三角形。之后通過空間點與投影點的映射關系將平面三角網映射回空間，從而實現三角形的空間構網。算法可簡單描述為以下四個步驟：

1) 由當前遍歷面為邊e(v0,v1)的鄰域Neib01=Neib0∪Neib1構建微切平面T，構造UV二維坐標系；

4　實例分析

采用C++語言編程在Microsoft visual studio 2008上實現以上算法，并且調用OpenGL庫函數顯示點云，實驗所用的計算機配置為Intel Core 2.30 GHz CPU，1.19 GB內存。對汽車擋泥板、兔子、馬鞍模型進行三角剖分，其點云數量分別是：6408、34 834、56 780個。為便于敘述，將Triangulation算法稱為方法1，Triangulation2為方法2，基于微切平面的生長法為方法3。方法1、方法2三角剖分思路相同，不同的是數據結構操作以及遞歸過程。三種方法對三種模型的三角化結果和耗時比較如表1和圖5-圖7所示。圖5-圖7中(a)是點云模型，(b)是方法1、方法2三角剖分結果，(c)是方法3結果。從表1和圖5-圖7中得出，方法1-方法3的三角化結果相同，但從耗時比較上來看，方法2三角剖分速度比方法1快，尤其是當點云數據較大時，方法2耗時幾乎是方法1的一半。此外，本文提出的基于微切平面的生長法即方法3耗時更小，與方法2相同數據結構實現方法3時耗時比傳統(tǒng)的生長法(方法2)耗時更小，對大數據點云模型，本文的生長法三角剖分效率更高。

表1　3個模型的實驗數據

圖5　擋泥板模型三角化結果

圖6　兔子模型三角化結果

圖7　馬鞍模型三角化結果

5　結　語

本文對Delaunay三角剖分算法中的生長法進行了詳細的分析，設計了一種順序存貯的Hash數據結構，實現臨時單純形對象的快速存取、插入和刪除等操作，采用Hash數據結構消去遞歸的生長法執(zhí)行時間更短。此外，本文提出了基于微切平面的生長法，將基于空間四面體的空球搜索降維至局部二維的空圓搜索，同樣采用Hash數據結構對基于微切平面法的生長法進行操作。微切平面生長法三角剖分效果與傳統(tǒng)方法效果相同，但微切平面生長法速度更快。對數據量為30 KB的兔子模型，三種方法的計算時間分別是8.5、4.3、3.2 s。

[1] Digne J,Cohen-Steiner D,Alliez P,et al.Feature-preserving surface reconstruction and simplification from defect-laden point sets[J].Journal of Mathematical Imaging and Vision,2014,48(2):369-382.

[2] Lhuillier M.2-Manifold Tests for 3D Delaunay Triangulation-Based Surface Reconstruction[J].Journal of Mathematical Imaging and Vision,2014,51(1):98-105.

[3] Gong W,Liu Y J,Tang K,et al.Approximate Delaunay mesh reconstruction and quality estimation from point samples[J].Journal of Computational and Applied Mathematics,2015,274(15):23-34.

[4] Gálvez A,Iglesias A.Particle swarm optimization for non-uniform rational B-spline surface reconstruction from clouds of 3D data points[J].Information Sciences,2012,192(1):174-192.

[5] Boissonnat J D.Geometric structures for three-dimensional shape representation[J].ACM Transactions on Graphics (TOG),1984,3(4):266-286.

[6] Edelsbrunner H,Mücke E P.Three-dimensional alpha shapes[J].ACM Transactions on Graphics (TOG),1994,13(1):43-72.

[7] Coakley E S,Rokhlin V.A fast divide-and-conquer algorithm for computing the spectra of real symmetric tridiagonal matrices[J].Applied and Computational Harmonic Analysis,2013,34(3):379-414.

[8] Xu Y,Liu L,Gotsman C,et al.Capacity-constrained Delaunay triangulation for point distributions[J].Computers & Graphics,2011,35(3):510-516.

[9] Carr J C,Beatson R K,Cherrie J B,et al.Reconstruction and representation of 3D objects with radial basis functions[C]//Proceedings of the 28th annual conference on Computer graphics and interactive techniques,2001:67-76.

[10] 趙建東,康寶生,康健超,等.改進的基于徑向基函數的曲面重建算法[J].西北大學學報:自然科學版,2012,42(5):744-748.

[11] 方林聰,汪國昭.基于徑向基函數的曲面重建算法[J].浙江大學學報:工學版,2010,44(4):728-731.

[12] 余杰,呂品,鄭昌文.Delaunay 三角網構建方法比較研究[J].中國圖象圖形學報,2010,15(8):1158-1167.

[13] 袁舒，楊烜.基于 Delaunay 三角網生長法的并行圖像插值方法[J].計算機應用與軟件,2012,29(3):62-68.

[14] 孟憲海,成文迪，徐博，等.基于 Voronoi 取小鄰近點集的 Delaunay 三角化方法[J].圖學學報,2013,34(6):36-41.

IMPROVEMENT AND COMPARATIVE ANALYSIS OF DELAUNAY TRIANGULATION ALGORITHM

Yuan XiaocuiWu LushenChen Huawei

(School of Mechanical and Electrical Engineering,Nanchang University,Nanchang 330031,Jiangxi,China)

Aiming at the problem of computation speed of Delaunay algorithm,we improved it from two aspects of data structure and algorithm.After analysing the algebraic topology of Delaunay triangulation,we designed a sequentially stored Hash data structure so that the temporary simplex object can be accessed,queried,inserted and deleted rapidly and in order.We took the activity analysis of edge object of simplex as the core,and employed Hash data structure to operate the recursion process of incremental algorithm removal.Besides we proposed the micro-tangent plane-based incremental algorithm,which decreases the dimension of the space tetrahedron-based blank sphere searching to local two-dimension blank circle searching.The triangulation experiments were carried out on car fender and bunny model,experimental results demonstrated that the triangulation results of recursion-removed and micro-tangent plane-based Incremental algorithms were the same as that of traditional Incremental algorithm,but the computation speed was higher than the traditional method.

Delaunay triangulationIncremental algorithmHalf spaceImplicit surface

2015-04-22。國家自然科學基金項目(51365037，51065021)。袁小翠，博士生，主研領域：逆向工程與數字圖像處理。吳祿慎，教授。陳華偉，講師。

TP391

A

10.3969/j.issn.1000-386x.2016.09.039