最優三角剖分演算法

⑴ 請問什麼是雙向插入法

逐點插入法進行構網，分析該演算法中制約運行速度的因素，在三角網拓撲關系、三角形的查找、LOP演算法（Local Optimization Procere）等方面進行了優化處理，使之有了較高的執行效率。
1 數據結構
在採用逐點插入進行Delaunay三角剖分的過程中，存在大量的查詢操作，利用數據結構提供三角形之間的拓撲信息，能夠有效地提高演算法效率。邊結構應包含點與三角形的信息，三角形結構應包含點與邊的信息，再考慮到構網過程中動態的數據更新，演算法中採用了雙向鏈表結構，以方便於剖分過程中新舊邊（三角形）的添加、刪除工作。
2 演算法原理
逐點插入法是Lawson在1977年提出的，該演算法思路簡單，易於編程實現。基本原理為：對於已有的三角網路，向其中插入一點，該點與包含它的三角形三個頂點相連，形成三個新的三角形，然後逐個對它們進行空外接圓檢測，通過交換對角線的方法來保證所形成的三角網為Delaunay三角網。
3 演算法基本步驟

a: 點在三角形內，b:點在三角形邊上。
在按a進行剖分時，添加了三條新邊及三個新三角形NT，刪除了舊三角形T。同樣，在b的情況中，記錄點所在的邊E，根據拓撲關系，找出該邊的左右相鄰三角形T1，T2，添加四條新邊和四個新三角形NT，刪除T1，T2以及邊E。
構網的關鍵是對新剖分出的三角形用LOP演算法進行優化，其基本過程為：新三角形與周圍的三角形構成共用同一條對角線的四邊形，逐個對四邊形中的兩個三角形進行空外接圓檢測，如果滿足空外接圓准則，則略過。如果不滿足，則用另一對角線替換現有對角線，在交換對角線後，又會產生相應的新四邊形，繼續進行空外接圓檢測。這種方法可連續進行,直到全部滿足空外接准則為止。這個過程是隨著數據點P的逐次插入，與三角剖分同時進行的，可以通過遞歸調用優化程序來實現。

⑵ 點集的Delaunay三角剖分方法

3.2.1.1 基本理論

B.Delaunay於1934年提出了Delaunay三角網格的概念，它是Voronoi圖(簡稱V圖)的幾何對偶圖，具有嚴格的數學定義和完備的理論基礎。

圖3.1 Voronoi圖(虛線)及對應的Delaunay三角剖分(實線)

3.2.1.1.1 Voronoi圖

假設V={v₁，v₂，…，v_N}，N≥3是歐幾里得平面上的一個點集，並且這些點不共線，四點不共圓。用d(v_i，v_j)表示點v_i與v_j間的歐幾里得距離。

設x為平面上的點，則:

區域V_(i)={x∈E²d(x，v_i)≤d(x，v_j)，j=1，2，…，N，j≠i}稱為Voronoi多邊形，也稱為該點的鄰域。點集中所有點的Voronoi多邊形組成Voronoi圖，如圖3.1所示。

平面上的Voronoi圖可以看做是點集V中的每個點作為生長核，以相同的速率向外擴張，直到彼此相遇為止而在平面上形成的圖形。除最外層的點形成開放的區域外，其餘每個點都形成一個凸多邊形。

3.2.1.1.2 Delaunay三角剖分

Delaunay三角形網格為V圖的幾何對偶圖。在二維平面中，點集中若無四點共圓，則該點集V圖中每個頂點恰好是3個邊的公共頂點，並且是3個Voronoi多邊形的公共頂點;上述3個Voronoi多邊形所對應的點集中的點連成的三角形稱為與該Voronoi頂點對應的Delaunay三角形，如圖3.1所示。如果一個二維點集中有四點共圓的情況，此時，這些點對應的Voronoi多邊形共用一個Voronoi頂點，這個公共的Voronoi頂點對應多於3個Voronoi多邊形，也就是對應於點集中多於3個的點;這些點可以連成多於一個的三角形。此時，可以任意將上述幾個點形成的凸包劃分為若干三角形，這些三角形也稱為和這個Voronoi頂點對應的Delaunay三角形。

所有與Voronoi頂點對應的Delaunay三角形就構成了Delaunay三角剖分。當無退化情況(四點共圓)出現時，點集的Delaunay三角剖分是唯一的。

3.2.1.1.3 Delaunay三角剖分的特性

Delaunay三角剖分具有兩個重要特性:

(1)最小角最大化特性:即要求三角形的最小內角盡量最大，具體地說是指在兩個相鄰的三角形構成凸四邊形的對角線，在相互交換後，6個內角的最小角不再增大，並且使三角形盡量接近等邊。

(2)空外接圓特性:即三角形的外接圓中不包含其他三角形的頂點(任意四點不能共圓)，該特性保證了最鄰近的點構成三角形，使三角形的邊長之和盡量最小。

3.2.1.2 常用演算法

Delaunay三角剖分方法是目前最流行的通用的全自動網格生成方法之一。比較有效的Delaunay三角剖分演算法有分治演算法、逐點插入法和三角網生長法等(Tsai，1993)，其中逐點插入法由於其演算法的簡潔性且易於實現，因而獲得廣泛的應用。其主要思路是先構建一個包含點集或區域的初始網格，再依次向初始網格中插入點，最後形成Delaunay三角剖分。

採用逐點插入法建立Delaunay三角網的演算法思想最初是由Lawson於1977年提出的(Lawson，1977)，Bowyer和Watson等先後對該演算法進行了發展和完善(Bowyer，1981;Watson，1981)。目前涌現出的大量逐點插入法中，主要為以Lawson演算法代表的對角線交換演算法和以Bowyer-Watson演算法代表的空外接圓法。

3.2.1.2.1 Lawson演算法

Lawson演算法的主要思想是將要插入的數據點逐一插入到一個已存在的Delaunay三角網內，然後再用局部優化演算法(Local Optimization Procere，LOP)優化使其滿足Delau-nay三角網的要求，其主要步驟如下:

圖3.7 Bowyer-Watson演算法剖分實例

⑶ 　靈活約束的四面體化

研究幾何模型所固有的物理特性的基礎步驟是網格化。基於簡單形體（即二維的三角形和三維的四面體）的結構化網格對自然物體的建模非常有意義，因為它們對填充由簡單邊界（即二維情形中的邊和三維情形中的三角形）所定義的區域的物性非常靈活。

自從提出了三角剖分一個多邊形的O（n）最優演算法（Chazelle,1991）後，從計算幾何學觀點來看，二維情況現在已可以很好地處理了。二維情況的其他重要演算法有Chew（1987）的約束Delaunay三角剖分方法和Ruppert（1993）的限制縱橫比三角剖分方法，前者給出一種來考慮給定連接關系的點集的三角剖分方法，後者可提供一個具有「好」的元素的網格。

所以，給定一個可在其端點彼此相交的邊界集合定義的二維模型，便可能利用三角形來填充所定義的全部區域。圖2.13顯示一個二維三角剖分模型的例子。

圖2.13二維模型的三角形化（Joel Conraud,1995）

●輸入角點；——輸入邊界

如果所有給定的邊界都在剖分結果中出現，這樣一個模型的三角剖分被稱作約束三角剖分；而如果一些邊界被分割後出現在剖分結果中，則被稱為匹配三角剖分（即在這些邊界上增加稱為「Steiner點」的一些點，用來分割它們）。為了考慮特定的連接，Chew修改了對Delaunay三角剖分的「空球標准」的定義，在文獻（Nackman and Srinivason,1991）中提出了一種在輸入的邊界上增加一些點，用來實現一種保持Delaunay特點的匹配三角剖分方法。

三維情況。三維情況要難得多，所有在二維中的一系列命題不再有效：

●有n條邊的多邊形的三角剖分有（n-2）個三角形；但一個多面體的四面體剖分具有不同數目的剖分結果。

●所有的多邊形都能被三角剖分；但有些多面體，如果其邊界拓撲不改變就不能被四面體填充；換句話說，不能對它們進行約束三維三角剖分，而只能進行匹配三維三角剖分的計算。

●二維點集的Delaunay三角剖分對所有可能的三角剖分的最小角最大化；但三維點集的Delaunay三角剖分做不到。

文獻（Ruppert and Seidel,1992）對非凸多面體的四面體化所涉及的困難有一個更詳細的描述。因此，如果有人想對一個涉及節點、邊和三角形的幾何結構四面體化，一般情況下需要增加Steiner點，它們的數目可能是很大的，導致一個復雜得多的四面體化。本節給出一個類似於匹配三維三角剖分的方法，稱為靈活約束的四面體化方法，這將大大減少增加點的數目。

問題。首先給出一些定義，然後描述我們感興趣的問題。

定義1：我們定義一個基於簡單邊界的模型M為三維空間中三角剖分定義的曲面的集合。曲面之一是閉合的，並且確定了我們研究的區域，而其他曲面有的確定了子區域，有的確定了不連續體（即面具有不閉合的邊界）。它們可以沿著三角形邊界彼此相交。

定義2：T是M四面體化的點集，當且僅當M的每個頂點是T的四面體頂點。

定義3：T是M的一個約束四面體化，當且僅當每個三角形（邊和頂點）是T的四面體表面（邊和頂點）。正如上面所寫的，不總是有一個與模型相聯系的四面體化。

定義4：T是M的一個匹配四面體化，當且僅當M的每個三角形：

●可以是T的四面體的表面；

●或者與T中四面體的表面有一個相關部分，這種情況意味著已增加了Steiner點。

從一個四面體化點集開始，怎樣修改使它成為M的一個可接受的四面體網格呢？

初步觀察。這個工作從以下觀察開始：

我們想對圖2.14所示的錐體作四面體化，在這個棱錐的表面，有三角形ABD、ABC。如果我們使用Delaunay三維點集三角剖分演算法（例如文獻Watson,1981中的方法）產生一個初始解，則由於d（C,D）＜d（A,B）（在一個Delaunay三角剖分中，對於存在4個共面點的情況，生成的兩個三角形遵循對角線最短原則），在點A和B之間連接的情況不會出現在棱錐的表面拓撲中，因此也就不會出現在四面體剖分網格中。

圖2.14錐體的Delaunay四面體化（Joel Conraud,1995）

如果我們忽略在共面區域這兩種三角剖分的不同之處，則生成的兩個四面體的邊界定義了由六個輸入三角形描述的同一個多面體。所以，即使拓撲結構不同，幾何形態是相同的。在一個靈活約束的四面體化方法中，上面的四面體網可被視作輸入的可接受的近似。

有時，一個僅基於模型節點剖分而成的四面體網格中的三角形，與一個基於邊界剖分而成的三角形網格不匹配。在研究這種情況時，則可發現兩種另外的方案：

（1）雖然，邊界三角形與四面體網格邊界不相交，但幾何形狀不相似。如果圖2.15中的點D沿ACB面的垂線向遠離E的方向移動，這個錐體被變形為一個凹多面體。按以前的情形生成兩個四面體，但這時在輸入三角形和四面體形成的多面體表面之間出現一個間隙。在圖2.15中的截面上可看出問題。

圖2.15（a）改變D的位置（Joel Conraud,1995）

圖2.15（b）（E,D,C）被兩個四面體共用（Joel Conraud,1995）

（2）邊界三角形與一個四面體邊界相交。在圖2.16的例子中，三角形（C,B,E）與四面體的邊（A,D）相交。要麼在交點處插入一個點，要麼進行網格的局部更新：可用四面體（A,B,C,E）和（D,E,B,C）代替四面體（A,D,C,E）、（A,D,E,B）、（A,D,B,C），產生一個具有8個輸入三角形的網格作為四面體的表面，但該表面不再是Delaunay的了。

圖2.16（a）由8個三角形確定的多面體（Joel Conraud,1995）

圖2.16（b）（A,B,C,D,E,F）的Delaunay四面體化（Joel Conraud,1995）

靈活約束四面體化的目標是僅當上述的第二種情況出現時，去組合測試模型的四面體化的點集並局部修改它。對付第一種情況（不論移不移D點），則用前面的方法，通過局部網格修改或增加點來處理（George等，1991）。

下面，我們給出靈活約束四面體化的如下定義：

定義5：T是一個靈活約束四面體化，當且僅當M的每個三角形：

●可以是T的一個四面體表面；

●或者與T中四面體的表面有一個重合部分；

●或者屬於一組n個相連的三角形G_Tgr1，通過它我們可以聯系一組n個四面體表面G_Tetra，這組四面體與相互連接的三角形組有相同的表面。

也就是說，靈活約束的四面體化只要求約束面與四面體表面的一部分幾何上重合，拓撲可以不同。

下面描述從一個模型的四面體化點集來構建一個靈活約束四面體化的演算法。

演算法。從一個簡單邊界的模型M和一個用M中的點通過無約束四面體剖分形成的四面體集T_{Point Set}（M）開始，演算法分為四步：

（1）計算無約束四面體剖分T_{Point Set}（M）中與M中的三角形不匹配的三角形集合S_UnMatched

for每個M中的t

if non-matched（t,M,T_{Point Set}）/＊用判斷函數判斷＊/

then S_UnMatched←S_UnMatched∪｛t}

endif

endfor

這里non-matched（）函數定義如下：

（t是基於三個頂點：v₁,v₂和v₃的三角形）

TetraFaces←φ

If edge（v₁,v₂）不存在於T_{Point Set}中Retur false

else TetraFaces←TetraFaces U tetrafaces（edge（v₁,v₂））

If edge（v₂,v₃）不存在於T_{Point Set}中Return false

else TetraFaces←TetraFaces U tetrafaces（edge（v₂,v₃）

If edge（v₃,v₁）不存在於T_{Point Set}中Return false

else TetraFaces←TetraFaces U tetrafaces（edge（v₃,v₁））

Return#TetraFaces=1

（2）將S_UnMatched分為相連的三角形子集

S_Connex←φ

for S_UnMatched中的每個三角形t

Connex←｛t和M中的鄰近t的三角形以及S_UnMatched的成員｝

S_UnMatched←S_UnMatchedConnex/＊把找到的三角形從S_UnMatched中去掉＊/

S_Connex←S_Connex∪｛Connex}/＊把找到的三角形合並到S_Connex中去＊/

endfor

（3）將每個S_Connex的子集分為由邊界既存在於M又存在於T_{Point Set}約束的三角形的子集

S_{BoundedConnex}←φ

for S_Connex中的三角形Connex的每一個集合

for Connex的每個三角形t

BoundedConnex←｛t和Connex中鄰近t的三角形以及既在M中又在T_{Point Set}沒有共享邊界的三角形｝

Connex←ConnexBoundedConnex/＊把找到的三角形從Connex中去掉＊/

S_{BoundedConnex}←S_{BoundedConnex}∪｛BoundedConnex｝

/＊把找到的三角形合並到S_{BoundedConnex}中去＊/

endfor

endfor

（4）對於每個子集BoundedConnex，其數據項的個數與基於同一邊界的四面體表面集合的數據項個數進行對比

for S_{BoundedConnex}的每一個三角形集合BoundedConnex do

① 計算集合edges（BoundedConnex,M）和edges（BoundedConnex,T_{Point Set}）

edges（BoundedConnex,M）←φ

edges（BoundedConnex,T_{Point Set}）←φ

BoundedConnexVertices←Vertices（BoundedConnex）

for BoundedConnexVertices的每一個頂點vdo

for基於v的T_{Point Set}的每一條邊e do

opposite=opposite（e,v）

if opposite是BoundedConnexVertices的成員

then edges（BoundedConnex,T_{Point Set}）←

edges（BoundedConnex,T_{Point Set}）Ue

ifM中在e和opposite之間存在邊

then edges（BoundedConnex,M）←

edges（BoundedConnex,M）Ue

endif

endif

endfor

endfor

② 計算共享edges（BoundedConnex,M）的四面體表面集合TetraFaces

TetraFaces←φ

for edges（BoundedConnex,M）的每一個邊e do

for基於e的T_{Point Set}的每一個四面體表面tfdo

for i from 1 to 3 do

neighber_edge←edge（tf,i）

if neighber_edge≠e

and member of edges（BoundedConnex,M）

then TetraFaces←TetraFaceUtf

endif

endfor

endfor

endfor

③ 計算在BoundedConnex的邊界上的邊集合BorderEdges

BorderEdges←φ

for BoundedConnex的每一個三角形t do

for i from 1 to 3 do

edge←edge（t,i）

trgl←neighber（t,i）

if trgl不是BoundedConnex的成員（或不存在）

then BorderEdges←BorderEdgesUedge

endif

endfor

endfor

④ BorderEdges定義一條閉合曲線，TetraFaces的某些元素在其外部，刪除它們

for TetraFaces的每一個四面體表面tf do

for i from 1 to 3 do

e←edge（tf,i）

if e不是BorderEdges的成員

then

for基於e的每一個四面體表面tf′ do

iftf′≠tf並且tf′為TetraFaces成員

then TetraFaces

←TetraFaces

{tf′以及TetraFaces中tf′的鄰接面

直到遇到BorderEdges的邊界）

endif

endfor

endif

endfor

endfor

⑤ 對比＃TetraFaces和＃BoundedConnex

if#TetraFaces=BoundedConnex

then OK

else局部網格修改（或增加Stiener點或同時進行）

endif

endfor

下一節我們將給出上述實現的結果。

實驗結果。上面所述演算法的一個C+＋版本已被結合到GOCAD軟體中，GOCAD軟體是一個自然物體的三維建模工具，特別適用於石油工業中地下情況的建模。

從模型點集的一個Delaunay三維三角剖分（在浮點運算中使用Watson的演算法。由Schroeder等，1990和George等，1992中所提到了同類型檢查方法使得數據結構的有效性在增加點後保持不變）開始，上一節的組合演算法使我們可以鑒別在T_{Point Set}和M之間的間隙。此時，還未應用上一節後面所講的局部網格修改。在後處理過程中，通過在區域內增加一些點來使最終的網格可用。

基本的例子。如果M被簡化為定義一個立方體的閉合三角剖分曲面，立方體的每個面分成兩個三角形，我們可得到以下結果（圖2.17）。

圖2.17正方體的四面體化（Joel Conraud,1995）

小窗口中顯示的四面體化的網格，元素在其中心周圍收縮，以便讓用戶看到結構內部。這個立方體的三個可見面的三角剖分，在約束面上和四面體邊界上是不同的。但兩個對象定義的幾何形體是一樣的。

實例。彩色圖版1.1是一個描繪具有超過8000個節點和17000個三角形的地質體的簡單三角剖分曲面。如果這個面的網格質量好的話，是不需要對四面體網格進行修改的。

彩圖版1.2中，M定義一個閉合的曲面（下面沒顯示），它和另外三個曲面定義四個地質層位。

彩色圖版2.1表示四面體網格的一個切片，我們可以看到三個曲面定義的界面的「傷疤」。

在具有256M的HP9000/800 K200上實現第二、第三個例子需要10min和12min的CPU時間，Delaunay三維三角剖分和M的三角形與四面體表面相比較花費的時間不及總時間的1/4。

討論。靈活約束的四面體化方法對地下情況建模的目的來說，提供了一個可接受的基於邊界模型的近似：在某些情況下，以用盡可能少的元素得到四面體化網格比考慮全部的三角形更重要。如果有人想在表面三角形和四面體T_{Point Set}的邊界表面之間得到1對1的關聯，則有必要從T_{Point Set}中提取約束和「靈活約束」的元素，從而得到一個新的基於邊界的模型M′，它與M幾何形狀相同，但拓撲不同。第二類型靈活約束改造的T_{Point Set}是一個M′約束的四面體化。

下一步工作將包括比較靈活約束四面體化方法與在T_{Point Set}中與M不匹配的每個三角形增加Steiner點的方法的計算效率比較。

⑷ 求一個用C++寫的Delaunay三角剖分間接實現Voronoi圖的代碼。最好有演算法說明謝謝！！ 急用！！

#include<iostream>
#include<cmath>
using namespace std;
#define N 30
typedef struct //定義點的結構體
{
int x,y;
}Point;
class point
{
private:
Point *v;
public:
int distance(Point i,Point j); //計算兩點的距離
int w(Point i,Point j,Point k); //計算三條邊的長度之和
void minWeightTriangulation(int n,int t[][N],int s[][N]); //用動態規劃計算最優值
void print(int s[][N],int i,int j); //輸出
};
int point::distance(Point i,Point j)
{
int s=(i.x-j.x)*(i.x-j.x)+(i.y-i.y)*(i.y-i.y);
return sqrt(s);
}
int point::w(Point i,Point j,Point k)
{
return distance(i,j)+distance(j,k)+distance(i,k);
}
void point::minWeightTriangulation(int n,int t[][N],int s[][N]) //用動態規劃計算最優值
{
int i=0;
int r=0;
int k=0;
for(i=1;i<=n;i++) t[i][i]=0;
for(r=2;r<=n;r++)
for(i=1;i<=n-r+1;i++)
{
int j=i+r-1;
t[i][j]=t[i+1][j]+w(v[i-1],v[i],v[j]);
s[i][j]=i;
for(k=i+1;k<j;k++)
{
int u=t[i][k]+t[k+1][j]+w(v[i-1],v[k],v[j]);
if(u<t[i][j])
{
t[i][j]=u;
s[i][j]=k;
}
}
}
}
void point::print(int s[][N],int i,int j)
{
if(i==j)
return;
print(s,i,s[i][j]);
print(s,s[i][j]+1,j);
cout<<"三角行：v"<<i-1<<"v"<<s[i][j]<<"v"<<j<<endl;
}
int main()
{
int n,i;
Point v[N]={0,0};
point triangle;
int t[N][N],s[N][N];
cout<<"輸入多邊形的頂點數：";
cin>>n;
for(i=0;i<n;i++)
{
cout<<"輸入第"<<i+1<<"點的坐標：";
cin>>v[i].x>>v[i].y;
}
triangle.minWeightTriangulation(n,t,s);
triangle.print(s,1,n);
return 0;
}