計算幾何第五周：完美三角剖分（Delaunay Triangulation）

來自專欄計算幾何筆記6 人贊了文章

1、Point Set Triangulation

• 目標：對一系列的點進行三角剖分

當我們有一系列點集時，我們需要對這些點集數據進行結構化。對點集進行結構化的一個重要的方法就是對這些點進行三角剖分。

這裡的對角線分為兩類，如上圖所示，深綠色的外面的邊對應於原來的凸包，另外一些內部的淺綠色的線則是凸包的內對角線。

所謂的點集的三角剖分，其實就是對這個點集所對應的凸包的區域進行剖分。

• 對於點集的三角剖分並不唯一

如果一個點集的點數目是n，如果凸包的size記做h，那麼任意一個極大的對角線的集合都是由3(n-1)-h條對角線組成；這些對角線剖分出來的三角形的數目為2(n-1)-h。

每新增一個點，大致就增加三條對角線，兩個三角形。

邊翻轉（Edge Flipping）：如上面不唯一的剖分所示，從左邊的豎著的那條紅色的內對角線變成右邊橫著的內對角線的這個操作稱作邊翻轉。每進行一次邊翻轉就會得到一個不同的剖分方案。

• 剖分演算法的界

如圖所示，若我們已經有兩個做好剖分的子集，那麼根據第一章凸包演算法找上下邊界的Zigzag方法，我們就可以得到兩個子集連接後的三角剖分，時間複雜度為O(nlogn)，那麼上界為O(nlogn)。

下界：將sorting問題規約到trangulation問題。下界為O(nlogn)

將sorting的一系列輸入映射到弧上的點，然後找到輸入的最大值和最小值，再在外面引入一個點構成三角剖分的問題輸入。對於這個點與弧上的點的這些對角線就構成了三角剖分，然後可以通過這個點的DCEL結構枚舉各條邊及其對應的鄰居們，這樣就得到了這些點在圓周上的排列的次序，也就得到了sorting的輸出。這樣就將sorting問題規約到了剖分問題。

因此點集三角剖分的下界是O(nlogn)

2、Delaunay Triangulation（簡稱DT）

Delaunay三角剖分定義：平面上的點集P是一種三角剖分，使得P中沒有點嚴格處於剖分後中任意一個三角形外接圓的內部。

對偶圖：對於一個Voronoi圖，若任何兩個site之間有一條非空邊界，那麼這兩個site之間連接一條邊。

由此得到的對偶圖，就是一個三角剖分。這是由於對於一個一般的Voronoi圖（即沒有四點共圓的情況出現），它的每一個Vertex的度數都恰好是3，也就是說每一個Vertex都恰好和三條邊關聯，那麼與同一個Vertex關聯的三條邊就會圍成一個三角形。這就說明在對偶圖中每一個face都是由三條邊圍成的三角形。

Delaunay三角剖分的複雜度是O(nlogn)。因為對偶變換是線性時間，從數學上將只需要把線性個頂點，線性個邊對應過去即可，而Voronoi圖的構建是O(nlogn)的時間，因此總體的剖分演算法的複雜度也就是O(nlogn)

3、Properties

Delaunay三角剖分與Voronoi圖是對偶關係，它們的性質也是有所聯繫。

（1）、空圓性（Empty Circle）

①對於Delaunay剖分中的任何一張face的外接圓必然是空的。如下圖所示的七個圓都是空的。

②在Delaunay三角剖分中，每一條邊之所以會成為一條邊，本質上是存在一個空圓以它為弦。

這個是它的充分條件。

左圖是剖分圖，右圖是相應的Voronoi圖。

（2）、最近鄰性

任何一條連接於最近鄰之間的邊都會被Delaunay剖分所採用，因為這裡頭會存在一個以這條邊為直徑（弦）的空圓。

而按照這種最近鄰關係所生成的圖又被稱為Nearest Neighbor Graph（最近鄰圖），它是對應的Delaunay剖分的一個子圖。NNG是一個有向圖，由上圖可知最近鄰關係不是對等的。

（3）、複雜度

在二維平面中，每增加一個點三角形的數目都會大概增加2，邊數增加3。可以說在二維上的Delaunay剖分中是一個線性規模的數據結構。但在三維的情況下這兩個指標最多會達到平房的量級，更高維的空間的一般結論也會達到2^d量級。

4、Proximity Graph

• Gabriel Graph（簡稱GG）

定義：對於給定點集S，若線pq是屬於Gabriel Graph的一條邊當且僅當

|pq|2=min{|pr|2+|rq|2 | r∈S}

幾何解釋：以pq為直徑的圓是空的。如圖所示

其他等價定義：若邊pq被Gabriel Graph所採用，它的充要條件是①其他的點r與p,q兩點形成的角不是鈍角②在同樣這個點集所對應的Voronoi圖中，在p和q之間必然有一段非空的公共邊界，且pq的連線必然會和實際上存在但是看不見的那條公共邊界存在一個實在的交點。

對於第②點的解釋：在Voronoi圖中，對應的p和q是site。如果確實p和q中間的連邊會和p和q這兩個Cell的公共邊界相交於一點，那一定是中點x。這個點是Voronoi edge上的一個點，那麼以它為中心會存在一個空圓且經過p和q；反過來如果p和q之間的這條邊的確會被Gabriel Graph採用並且相應的有一個空圓，那麼也說明這個圓心是具有Voronoi edge的必要條件的一個點，即這個圓心就在某一條Voronoi edge上，且這條edge就是介於p和q所對應的Cell之間的那一段公共邊界上。

結論：Gabriel Graph中的任何一條邊也就屬於同樣的一個點集所對應的Voronoi圖中的一條邊，因此Gabriel Graph必然是Delaunay剖分的一個子圖。

構造Gabriel Graph的複雜度：在擁有Voronoi圖的基礎上可以花費O(n)時間來得到它，而構造Voronoi圖需要O(nlogn)，因此累計是O(nlogn)的時間即可構造Gabriel Graph。

• Relative Neighborhood Graph（相對鄰居圖，簡稱RNG）

定義：|pq| = min{max(|pr|, |rq|) | r∈S}

幾何解釋：若線pq屬於RNG，那麼以pq為半徑，分別以p和q為圓心作出的兩個圓的公共部分必定是空的。如圖所示的兩個圓的公共的梭形部分。

性質：RNG一定是Gabriel Graph的子圖。由上圖所示，顯然以pq為直徑的圓一定是空的。

構造RNG的複雜度：O(nlogn)

• 構造GG或RNG的下界：O(nlogn)

證明：考慮一維的情況，將這個問題規約到sorting。對於sorting輸入的一系列點，其實就相當於一維空間中GG和RNG的輸入。如圖所示，根據GG的定義，相鄰兩個點的連線必然會被GG所採用。那麼構造出來的GG的點的排列是按照緊鄰的規則來確定，而這個規則恰好就是對應的排序的結果，因此可以規約到sorting問題上。RNG圖亦然。

5、Euclidean Minimum Spanning Tree（歐氏最小生成樹，簡稱EMST）

定義：最小生成樹兩個節點之間的權重就是它們之間的歐氏距離。

如圖所示，給了一系列的點用適當的方式連接起來，使得這個連接的成本達到歐氏意義上的最小。那麼這個圖實際上是一個隱式的圖，因為所有點對之間的距離都是通過歐氏距離隱式的給出。這樣，我們的計算輸入是一幅完全圖。

由於有了歐氏距離這個條件，EMST是我們前面提到的RNG的子圖。

證明：思路是任何一條不屬於RNG的邊也必不屬於EMST。

如圖所示，圖中的pq必不被RNG所採用，假設它被EMST採用，那麼將pq剪除後會把EMST圖分為不連通的兩個部分。而r作為EMST中的一個點，它至少會和p或q中的一個點連通。不失一般性，假設r與p連通，當我們將pq剪除後再連接rq，那麼我們會得到一棵新的樹，這棵樹與之前的EMST的差別就在於pq和rq的權重不同。由於權重的定義是歐氏距離，那麼從幾何上可知rq的權重必然小於pq的權重。這說明這棵新的樹的權重比之前的最小生成樹的權重還要小，矛盾。因此不被RNG採用的邊一定不被EMST採用。

在有上述結論EMST是RNG的子圖後，那麼EMST也必然是DT的子圖。那對於EMST的構造就可以先花費O(nlogn)的時間構造DT，然後再套用Kruskal或Prim演算法即可。

總覽：

6、Euclidean Traveling Salesman Problem（簡稱ETSP）

定義：給出的所有點之間的距離採用歐氏距離，在這樣的意義下找到一條環路使得距離和最小。

已證明：ETSP依然是NP-hard問題，但有一些簡單的近似方法。我們可以在O(n)時間內從一個已經構造好的EMST中得到一個ETSP的近似解，且近似程度最多不過最優解的兩倍。

演算法：首先構造點集對應的EMST，然後利用DCEL中類似的方法構造一個初始的環路（如圖所示綠色部分）。然後從一個特定的點出發如最左側的點，沿著初始環路往前走，能往前走就走直到它試圖回到以前已經訪問過的一個點的時候，就跳過這個點並嘗試下一個點，最終會回到起點，並且總的路程（tour）會構成一個環路（如圖所示的藍色虛線）。

在沒有優化的情況下，這條tour的總權重|W| = 2 * |EMST| < 2 * |ETST|

定義：對於給定的一系列點集，找到一個三角剖分，使得它的權重之和最小。

MWT與DT並不是同一碼事，如下圖，MWT會選擇豎直的那條內對角線，而DT會選擇水平的內對角線。

MWT問題的是NP-Hard問題。

8、Construction

• 基本操作：edge flip（邊翻轉）
• 前置工作：

對於每一個三角剖分所對應的角的數目為6n-3h-6個。

一個指標：將一個三角剖分中的所有角按照從小到大排列，構成一個向量用來度量和評價這個三角剖分。如下圖所示

若採用AC作為內對角線剖分，則角度序列是 <15,30,60,75,75,105>

若採用BD作為內對角線剖分，則角度序列是<15,30,45,60,75,135>

用字典序比較這兩個角度序列，我們可以說AC＞BD

Delaunay三角剖分存在一個趨勢：儘可能使每一個三角形都勻稱均衡。如下圖所示，上圖採用的是DT，下圖則是其他方案。因為要保證空圓性，上圖剖分的三角形要比下圖更為均衡。

從數學上講，Delaunay三角剖分會儘可能使得剛才所說的角度序列極大化，也可以說是儘可能使小的角更大，又被稱為Maximizing the minimum angle.

解釋：由上圖所示，我們來考慮ad,ab,bc和cd對應的四個角，由幾何知識得：θab＜φab，θbc＜φbc，θcd＜φcd，θda＜φda。因此由DT得到的角相對而言比其他方法得到的要大。

問題：我們總共有6對角，為何只考慮其中的四對角？

個人思考：雖然四邊形內角和是360度固定，有角變大自然就有角變小，但是因為我們在進行比較的時候是從小到大按字典序比較，只要偏小的角度比其他方法得到的小角度要大，那麼後面較大的角度比其他方法得到的大角度小也無所謂了，畢竟前面就已經比較出大小了。

• 策略總結

（1）確定性演算法

①對偶：通過Voronoi圖的對偶圖可以在O(nlogn)+O(n)=O(nlogn)的時間內得到。

②分而治之

③平面掃描

（2）隨機化演算法：蠻力演算法+隨機化處理

9、RIC With Example

RIC：Randomized Incremental Construction

思想：在一個已經成形的現有框架下再增加一個點進行構造。

（1）、點定位

對於給定的一個新的點，判斷它處於哪個face中。

（2）、In-Circle test

對於新引入的這個點，判斷它的哪個三角剖分是屬於Delaunay剖分的。

當如上圖的右圖的剖分不屬於Delaunay剖分時，將該邊進行邊翻轉（edge flip）翻到左邊的那種情況。但當進行了這個操作後，要對那些粉紅色的邊進行判斷是否還是屬於Delaunay剖分，這些邊又被稱為嫌疑邊。

（3）、Frontier

對嫌疑邊做空圓檢測，若不符合，那麼繼而進行邊翻轉，再檢查其餘的嫌疑邊。這些嫌疑邊實際上是首尾相連構成一個環的，有一個中心G，是一個星形多邊形。這個星形多邊形又稱為Frontier Polygon。

（4）、Convergence

對於這個多邊形的角度序列，每做一次邊翻轉，角度序列都會有所遞增，而這個角度序列是有一個上界的，因此這個操作存在一個上界。當抵達上界時演算法結束。

10、Randomized Incremental Construction

• 基本操作

①、查找新插入的點在哪個三角形的內部

②、把新引入的點在幾何結構上縫合到其中去

③、邊翻轉

• 遞歸版本：偽代碼如下

Insert (p){ //首先做point detection,找到包含p的那個三角形abc T(a,b,c) = TriangleContaining(DT, p); //插入三條邊pa,pb,pc，構造初始的三角剖分 connect(p,a); connect(p,b); connect(p, c); //對ab, bc, ac做排查看是否符合Delaunay剖分，有必要的話需要進行交換 swapTest(p,a,b,c);}swapTest(p, a, b, c){ //針對剛剛插入的點p,然後相對於一條懷疑邊進行外接圓的空圓檢測 sTest(p, a, b); sTest(p, b, c); sTest(p, c, a);}sTest(p, a, b){ //找到對岸的點x，就是ab這條邊的twin edge所附著的你賬面的對頂的角的點 x = rightSite(a, b); //若不存在x則返回 if(!x) return; //若存在x，則判斷這個x是否包含在pab所構成的外接圓上 if(inCircle(p, a, b, x) { //如果落在裡面，那麼翻轉邊ab，連接成邊px flipEdge(a, b, p, x); //翻轉過後要繼續檢測兩條新的懷疑的邊。 sTest(p, a, x); sTest(p, x, b); }}

• 迭代版本：偽代碼如下

swapTest(p, a, b, c){ //定義一個隊列Q，最初包含三條邊 queue Q = {(a, b), (b, c), (c, a)}; while(!empty(Q)) { //取出一條邊 (a, b) = Q.dequeue(); //找到對岸的點x x = rightSite(a, b); //若不存在則演算法繼續 if(!x) continue; //若存在則判斷x是否在外接圓內 if(inCircle(p, a, b, x)) { //若在，則翻轉邊ab變成邊px flipEdge(a, b, p, x); //將邊ax和xb入隊 Q.enqueue((a, x), (x, b)); } }}

• In-Circle Test

判斷點是否在三角形外接圓內部。

若InCircle(a,b,c,p)>0，則說明這個點落在圓內部；

若InCircle(a,b,c,p)<0，則說明這個點落在圓外部；

若InCircle(a,b,c,p)=0，則說明這個點落在圓周上。

時間複雜度O(1)。

• Point Location

首先一開始的時候將所有的未插入的點進行歸類，每一類用一個桶結構（bucket）來進行存儲，具體的分類就是按照它們在已有的三角形中分別歸屬於哪個三角形。

如上圖所示，abc同屬於一個三角形內，那麼這三個點就放在一個桶里。

在進行邊翻轉的時候會導致三角形結構的變更。如下圖，一開始是水平切分成兩個三角形，那麼abc是一個桶，de是另一個桶；若發生邊翻轉，那麼這兩個桶消亡，重新生成新的兩個桶，一個桶放置bce，一個桶放置ad。

11、RIC Analysis

• 複雜度分析：主要包含邊的變化的複雜度和桶結構修正的複雜度。

後向分析法：在對於已知的結果情況下對前面的時間成本進行分析，站在一個過去完成時，隨機的在整個計算過程中的某一個時刻叫停，然後反問這個演算法在剛剛過去已經發生的這樣一步迭代過程消耗了多少時間。

需要具備的條件：①最終的結果必須是唯一確定的②在中間過程中的任何一個時刻的結果都是確定的③前面的任何一個操作都是概率均等的。

• 邊變化：對於每一次頂點的插入，邊的變化是常數O(1)的時間。

使用後向分析的條件對應：①假設所有的點都在一般性位置沒有退化，那麼最終結果是確定的②在中間的任何一個時刻的k個點的臨時局部的Delaunay剖分都是確定的③前面k個點的插入從概率上來說每個點是最後插入的概率都是相等的。雖然不知道前k個點究竟是哪個點作為最後一個點插入，但是只需要把每個點所作為最後一個元素插入所對應的時間成本平均起來就可以得到這樣一個期望值。

我們將假設是最後一個插入的點稱為點p，分析每一個點作為這個點p相應的需要做多少次邊翻轉以及相應的邊操作。

假設這個點G是作為最後一個點即p點被插入，它會引起的修改操作是：①初始插入：即引入三條新增的邊；②：邊翻轉，每一次邊翻轉會生成一條新的邊。對於這個圖而言，在第二步的時候會有兩次邊翻轉，所以總共有3+2次邊操作。

分析：因為所有的點都可能是作為最後一個點被插入，理論上而言我們應該把所有點的可能統計出來再求平均。具體的來講，假設上圖的G是作為最後一個點被插入，如果我們能把它精確所需要的時間準確的估計出來就足夠了，若這個方法能通用就可以施加到所有的點上去。因此現在我們將目光就放在一個點上。

分析這個點G，在每經過一次邊的翻轉時，都會有一條老的邊死掉並且引入一條新的邊，同時會有兩條新的邊成為需要檢查的嫌疑邊，因此對於這個Frontier而言每做過一次flip後size都是-1+2=1，這說明這個Frontier Polygon忠實的記錄了我們此前進行過多少次操作。巧合的是如果算上初始生成的三條邊這個數恰好相當於Frontier Polygon的size。注意到點G引起的edge change恰好是5次，而這個Frontier Polygon的size也是5，因此這個polygon最後的size就度量了我們此前如果G是最後一個點被插入的話它所花費的時間。而前面我們說過這個Frontier Polygon是一個星形多邊形，G就是它的核，我們就可以直接得出結論說這個點G的度數就度量了花費的時間。

結論：任何一個點作為最後一個點被插入的成本就等於它在當下這張圖中所具有的度數。那麼平均時間就是當下所有點的度數的平均值。

由歐拉公式知，從漸進的意義上來講一個平面圖的邊數大概是頂點數的3倍。而每一條邊貢獻2度，因此總共不過6n度數，分配到n個點每個點的平均度數不會超過6。也就是說在最後一步插入的過程中，從平均意義上講所需要執行的邊操作無非是O(6)即O(1)，也就是之前所說的常數時間。

• 桶變化（rebucketing）

思路：考慮某一個點，在整個演算法的執行過程中這個點從期望的角度講會參與rebucketing多少次？

結論：O(logn)

假設點pi是第i個點，即截止到它已經插入了i個點，那麼剩下n-i個點待插入。而在剩下的n-i個點鐘每一個點從期望的意義上講從概率上都是3/i的概率會參與剛剛過去的最後那次迭代中參與了rebucketing，它們所歸屬的那個三角形以及相應的bucket有變化。

對於上圖，假設在最後的那步迭代中q曾經參與一次rebucketing，它的必要條件是在當下q所屬的那個三角形必然是在剛剛最後那步迭代中最新創建的。可以再將這個條件轉化為另一個必要條件：如果這個三角形確實是剛才第i次迭代中創建的，那麼這個三角形肯定會與最新插入的頂點相關聯，也就是說新插入的這個點必然是這個三角形的三個頂點之一，即要麼a要麼b要麼c作為最後一個點引入，才會引起q的rebucketing。因此這個概率是3/i。

總體期望：

因此最後結果是O(nlogn)。