筆記：Delaunay三角剖分(Delaunay Triangulation)相關知識

最近接觸到計算Delaunay三角剖分的問題，也算是計算幾何的一個經典問題了。按照別人的演算法，也自己實現了個（原始碼下載），發現點集大的時候，程式計算起來特慢。後來分析發現，別人程式號稱的都是O(nlogn)的，我的卻成了O(n*n)的，演算法都是一樣，後來才發現是資料結構的問題，看來程式＝演算法＋資料結構，有道理。閑著，就整理了些相關知識，組織如下：

1.Delaunay三角剖分&Voronoi圖定義
2.計算Delaunay三角剖分的演算法及分析
3.例子程式&代碼

大話
點集的三角剖分（Triangulation），對數值分析（比如有限元分析）以及圖形學來說，都是極為重要的一項預先處理技術。
尤其是Delaunay三角剖分，由於其獨特性，關於點集的很多種幾何圖都和Delaunay三角剖分相關，如Voronoi圖，EMST樹，Gabriel圖等。Delaunay三角剖分有幾個很好的特性：
1.最大化最小角，“最接近於規則化的“的三角網。
2.唯一性（任意四點不能共圓）。

概念及定義
二維實數域(二維平面)上的三角剖分

定義1：假設V是二維實數域上的有限點集，邊e是由點集中的點作為端點構成的封閉線段, E為e的集合。
那麼該點集V的一個三角剖分T=(V,E)是一個平面圖G，該平面圖滿足條件：
1.除了端點，平面圖中的邊不包含點集中的任何點。
2.沒有相交邊。
3.平面圖中所有的面都是三角面，且所有三角面的合集就是點集V的凸包。
那什麼是Delaunay三角剖分呢?不過是一種特殊的三角剖分罷了。從Delaunay邊說起。

Delaunay邊
定義2：假設E中的一條邊e（兩個端點為a,b），e若滿足下列條件，則稱之為Delaunay邊：
存在一個圓經過a,b兩點，圓內不含點集V中任何的點，這一特性又稱空圓特性。

Delaunay三角剖分
定義3：如果點集V的一個三角剖分T只包含Delaunay邊，那麼該三角剖分稱為Delaunay三角剖分。

我們看幾個圖：

由上面的圖引出Delaunay三角剖分的另一種定義：
定義4：假設T為V的任一三角剖分，則T是V的一個Delaunay三角剖分，當前僅當T中的每個三角形的外接圓的內部不包含V中任何的點。

Voronoi圖

定義5：V的Voronoi圖是由多邊形地區的集合（有些地區可能不是閉合的），該地區僅含點集中的一個點v，地區中的任何位置到v的距離都比該位置到點集中其它所有點的距離短。

由Voronoi圖和Delaunay三角剖分的關係，可以引出另一個Delaunay三角剖分的定義：
定義6：將Voronoi圖相鄰地區（共邊的地區）中的點串連起來構成的圖，稱為Delaunay三角剖分。
如：

概念部分到此，下面看看怎麼求Delaunay三角剖分。

計算Delaunay三角剖分

問題1：計算二維Delaunay三角剖分
問題輸入：二維實數域上的點集V
問題輸出：Delaunay三角剖分DT = (V, E).

演算法

由不同的定義對應有不同的演算法。用得較多的是基於定義3或4的演算法。
目前常用的演算法又分為好幾種，被不同的傢伙發現。什麼掃描線法（Sweepline），隨機增量法（Incremental），分治法（Divide and Conquer）啊等等。

隨機增量法（Incremental）

其中，隨機增量法較為簡單，遵循增量法的一貫思路，即按照隨機的順序依次插入點集中的點，在整個過程中都要維護並更新一個與當前點集對應的Delaunay三角剖分。考慮插入vi點的情況，由於前面插入所有的點v1,v2,...,vi-1構成的DT(v1,v2,...,vi-1)已經是Delaunay三角剖分，只需要考慮插入vi點後引起的變化，並作調整使得DT(v1,v2,...,vi-1) U vi成為新的Delaunay三角剖分DT(v1,v2,...,vi)。
(插入調整過程：首先確定vi落在哪個三角形中（或邊上），然後將vi與三角形三個頂點串連起來構成三個三角形（或與共邊的兩個三角形的對頂點串連起來構成四個三角形），由於新產生的邊以及原來的邊可能不是或不再是Delaunay邊，故進行邊翻轉來調整使之都成為Delaunay邊，從而得出DT(v1,v2,...,vi)。)

其Pseudocode如下：

Algorithm IncrementalDelaunay(V)
Input: 由n個點組成的二維點集V
Output: Delaunay三角剖分DT
1.add a appropriate triangle boudingbox to contain V ( such as: we can use triangle abc, a=(0, 3M), b=(-3M,-3M), c=(3M, 0), M=Max({|x1|,|x2|,|x3|,...} U {|y1|,|y2|,|y3|,...}))
2.initialize DT(a,b,c) as triangle abc
3.for i <- 1 to n
    do (Insert(DT(a,b,c,v1,v2,...,vi-1), vi))  
4.remove the boundingbox and relative triangle which cotains any vertex of triangle abc from DT(a,b,c,v1,v2,...,vn) and return DT(v1,v2,...,vn).

Algorithm Insert(DT(a,b,c,v1,v2,...,vi-1), vi)
1.find the triangle vavbvc which contains vi // FindTriangle()
2.if (vi located at the interior of vavbvc) 
3.    then add triangle vavbvi, vbvcvi and vcvavi into DT // UpdateDT()
    FlipTest(DT, va, vb, vi)
    FlipTest(DT, vb, vc, vi)
    FlipTest(DT, vc, va, vi)
4.else if (vi located at one edge (E.g. edge vavb) of vavbvc)
5.    then add triangle vavivc, vivbvc, vavdvi and vivdvb into DT (here, d is the third vertex of triangle which contains edge vavb) // UpdateDT()
    FlipTest(DT, va, vd, vi)
    FlipTest(DT, vc, va, vi)
    FlipTest(DT, vd, vb, vi)
    FlipTest(DT, vb, vc, vi)
6.return DT(a,b,c,v1,v2,...,vi)
   
Algorithm FlipTest(DT(a,b,c,v1,v2,...,vi), va, vb, vi)
1.find the third vertex (vd) of triangle which contains edge vavb // FindThirdVertex()
2.if(vi is in circumcircle of abd)  // InCircle()
3.    then remove edge vavb, add new edge vivd into DT // UpdateDT()
    FlipTest(DT, va, vd, vi)
    FlipTest(DT, vd, vb, vi)

Algorithm InCircle(va, vb, vd, vc)
if |a^b^, c^b^, d^b^| > 0, vd is NOT in the circumcircle of abc //a^座標為(ax,ay,ax*ax+ay*ay),a^b^為向量
if |a^b^, c^b^, d^b^| == 0, vd is on the circumcircle of abc
if |a^b^, c^b^, d^b^| < 0, vd is in the circumcircle of abc

Note: Details of InCircle():
如：

拋物線P (z = x*x + y*y), 現有x-y平面上a,b,c三點，現在要測試點d是否在a,b,c的外接圓內，將a,b,c,d四點lift up,與P相交於a^,b^,c^,d^，則有：
case 1：如果d正好位於外接圓上，則d^必在由a^,b^,c^三點構成的平面H上。此刻向量a^b^, c^b^, d^b^構成的矩陣的行列式det==0
case 2：如果d位於外接圓內，即d'，則d'^必在由a^,b^,c^三點構成的平面H下面。即向量a^b^, c^b^, d^b^構成的矩陣的行列式det<0
case 3：如果d位於外接圓外，即d''，則d''^必在由a^,b^,c^三點構成的平面H上面。即向量a^b^, c^b^, d^b^構成的矩陣的行列式det>0

複雜度分析

問題的規模為點集中的點的總個數n（沒有重合的點），迴圈內的基本的操作有：
1.尋找插入點所在的三角形(FindTriangle())
2.測試點是否在外接圓內(InCircle())
3.更新三角網(UpdateDT())
4.尋找共測試邊的第三頂點(FindThirdVertex())

考慮最壞情況：

時間複雜度T ＝ T(addboudingbox()) + Sum(T(insert(i), i=1,2,...,n)) + T(removeboundingbox)
因為addboudingbox()和removeboundingbox不隨n變化，是個常量。T(addboudingbox()) = O(1), T(removeboundingbox()) = O(1).

T = Sum(T(insert(i), i=1,2,...,n)) + O(1) + O(1).
 
考慮插入第i個的點的時間:
 T(insert(i)) = T(FindTriangle(i)) + T(UpdateDT(i)) + K*T(FlipTest(i)) (K為常數)

故T = Sum(T(FindTriangle(i)), i=1,2,..,n) + Sum(T(UpdateTD(i)), i=1,2,..,n) + K*Sum(T(FlipTest(i)), i=1,2,..,n)

挨個考慮：
FindTriangle(i)是要找出包含第i個點的三角形，由歐拉公式知道，平面圖的面數F是O(n), n為頂點數，故此時總的三角形數是O(i)的。所以問題相當於在O(i)個記錄中尋找目標記錄，如果不藉助特殊的資料結構，按照一般順序尋找，需要O(i)的時間。
T(FindTriangle(i)) = O(i),故：Sum(T(FindTriangle(i)), i=1,2,..,n) = O(n*n)

UpdateTD(i)是更新三角網資料結構，插入和刪除三角形到當前的三角網是個常量操作，因為已經知道插入和刪除的位置。
T(UpdateDT(i)) = O(1),故：Sum(T(UpdateTD(i)), i=1,2,..,n) = O(n)

FlipTest(i)比較複雜些，是個遞迴過程。細分為：T(FlipTest(i)) = T(FindThirdVertex(i)) + T(InCircle(i)) + T(UpdateDT(i)) + 2*T(FlipTest(j))。(這裡，用j來區分不同的深度)
因為InCircle(i)是測試點是否在外接圓內，是個常量操作，不隨點數變化而變化。故T(InCircle(i)) = O(1), 又T(UpdateDT(i) = O(1)(見上)
FindThirdVertex(i)是尋找目標點，在O(i)個記錄中尋找目標記錄，如果不藉助特殊的資料結構，按照一般順序尋找需要O(i)的時間。T(FindThirdVertex(i)) = O(i).
剩下的是遞迴調用FlipTest的過程，不過還好，因為FlipTest的遞迴深度只有常數級(設為K)。正比於點在三角網中的度數(degree)。故FlipTest(i)最多調用的次數為2*2*,...,2 = K',還是常數。
故T(FlipTest(i)) = O(i) + O(1) + O(1) + K'*(O(i) + O(1) + O(1)) = O(i) + O(1) + O(1) .
Sum(T(FlipTest(i)), i=1,2,..,n) = O(n*n) + O(n) + O(n)

綜上，最壞情況下演算法總時間複雜度 T = O(n*n) + O(n) + K*(O(n*n) + O(n) + O(n)) + O(1) + O(1) = O(n*n)
其中，關鍵的操作是FindTriangle()和FindThirdVertex()，都是n*n次的。

在網上很多資料說隨機增量法是O(nlogn)的，但是分析下來，卻是O(n*n)的。後來看到別人的實現，原來是用的別的資料結構來儲存三角網，減少了FindTriangle()和FindThirdVertex()的複雜度。使得某次尋找三角形和共邊三角形的第三頂點能在O(logn)，而非O(n)的時間內實現。這樣，總的尋找的時間為 O(log1)+O(log2)+,...+O(logn) = O(nlogn)。程式=演算法+資料結構，看來一點沒錯。
比如說用DAG，Quad-edge等，都能達到O(nlogn)的複雜度。

分治法（Divide and Conquer）

據說是現在實際表現最好的。

掃描線法（Sweepline）

有空再看看其演算法思路。

例子程式&代碼

網上由很多代碼關於計算Delaunay的庫，如Triangle, Qhull,看的煩，自己寫了個，很簡陋，基本上就是隨機增量法虛擬碼的直譯，還沒時間最佳化。

有空再整理下約束性Delaunay三角剖分相關知識。

to be continued...