二維Delaunay（德洛內）三角網的matlab實現

1.Delaunay三角網的概述

德洛內(Delaunay)三角網的定義: 它是一系列相連的但不重疊的三角形的集合, 而且這些三角形的外接圓不包含這個面域的其他任何點。它具有兩個特徵：

(1) 每個德洛內(Delaunay) 三角形的外接圓不包含面內的其他任何點, 稱之為德洛內(Delaunay) 三角網的空外接圓性質, 這個特徵已經作為建立德洛內(Delaunay) 三角網的一項判別標準;

(2) 它的另一個性質最大最小角性質: 每兩個相鄰的三角形構成的凸四邊形的對角線，在相互交換後，六個內角的最小角不再增大。

Delaunay 三角網的優點是結構良好, 資料結構簡單, 資料冗餘度小, 儲存效率高, 與不規則的地面特徵和諧一致,可以表示線性特徵和迭加任意形狀的區域邊界, 易於更新,可適應各種分佈密度的資料等; 它的侷限性是, 演算法實現比較複雜和困難, 但現在已經有了較多成熟的實現演算法。

以上直接複製百度百科的。
下面是一個自己做的Delaunay三角網的劃分的演示例項

2.Delaunay三角網的演算法

本文的演算法採用的是改進型的駐點插入演算法，也就是Bowyer-Watson演算法。由於前面具體的演算法概述和參考的文章大致相同，所以就不再多說了（其實就是懶得複製貼上了）。
這個演算法採用的比較廣泛，效果比較好，而且在原網格中插入新點會更方便，但是缺點是不太好理解。

主要思路是：
1。先構建一個超級三角形，把所有點包圍起來
2。之後依次一個一個插入點。(a)在網格中插入新點 (b)相關三角形畫外接圓，記錄下所有外接圓包含新點的三角形 (c)刪除影響三角形的公共邊 (d)重新規劃新邊
3。記錄三角形至網格中
4。重複步2，直至所有點插入
5。刪除超級三角形以及和其相關的三角形。

演算法的虛擬碼如下：

input: 頂點列表(vertices)　　　　　　　　　　　　　　　　　　　　  　//vertices為外部生成的隨機或亂序頂點列表
output:已確定的三角形列表(triangles)
　　　　初始化頂點列表
　　　　建立索引列表(indices = new Array(vertices.length))　　　　//indices陣列中的值為0,1,2,3,......,vertices.length-1
　　　　基於vertices中的頂點x座標對indices進行sort　　  　　　　　  //sort後的indices值順序為頂點座標x從小到大排序（也可對y座標，本例中針對x座標）
 

　　　　確定超級三角形
　　　　將超級三角形儲存至未確定三角形列表（temp triangles）
　　　　將超級三角形push到triangles列表
　　　　遍歷基於indices順序的vertices中每一個點　　　　　　　　　  　//基於indices後，則頂點則是由x從小到大出現
　　　　　　初始化邊快取陣列（edge buffer）
　　　　　　遍歷temp triangles中的每一個三角形
　　　　　　　　計算該三角形的圓心和半徑
　　　　　　　　如果該點在外接圓的右側
　　　　　　　　　　則該三角形為Delaunay三角形，儲存到triangles
　　　　　　　　　　並在temp裡去除掉
　　　　　　　　　　跳過
　　　　　　　　如果該點在外接圓外（即也不是外接圓右側）
　　　　　　　　　　則該三角形為不確定        　　　　　　　　　     //後面會在問題中討論
　　　　　　　　　　跳過
　　　　　　　　如果該點在外接圓內
　　　　　　　　　　則該三角形不為Delaunay三角形
　　　　　　　　　　將三邊儲存至edge buffer
　　　　　　　　　　在temp中去除掉該三角形
　　　　　　對edge buffer進行去重
　　　　　　將edge buffer中的邊與當前的點進行組合成若干三角形並儲存至temp triangles中
　　　　將triangles與temp triangles進行合併
　　　　除去與超級三角形有關的三角形
end

這篇文章（三角剖分演算法(delaunay) http://www.cnblogs.com/zhiyishou/p/4430017.html）有一個用三個點的例項來演示這個程式執行的結果，如果可以的話，也推薦用三個點或4個點來檢測一下自己程式執行的效果。

3。Delaunay三角剖分的matlab實現

下面是matlab實現的原始碼，之後會解釋一下實現的思路

clear
N=1000;
%點隨機
xdot=rand(N,2);
%點按圓形隨機
%r=rand(N,1).^0.3;
%theta=rand(N,1)*2*pi;
%xdot=[r.*cos(theta),r.*sin(theta)];
%點按圓形規則
% r=0:0.001:1;
% r=r.^0.8;
% theta=0:0.1:1000*0.1;
% theta=theta.^1.2;
% xdot=[(r.*cos(theta))',(r.*sin(theta))'];


%1Delaulay三角形的構建

%整理點，遵循從左到右，從上到下的順序
xdot=sortrows(xdot,[1 2]);

%畫出最大包含的三角形
xmin=min(xdot(:,1));xmax=max(xdot(:,1));
ymin=min(xdot(:,2));ymax=max(xdot(:,2));
bigtri=[(xmin+xmax)/2-(xmax-xmin)*1.5,ymin-(xmax-xmin)*0.5;...
    (xmin+xmax)/2,ymax+(ymax-ymin)+(xmax-xmin)*0.5;...
   (xmin+xmax)/2+(xmax-xmin)*1.5,ymin-(xmax-xmin)*0.5];

xdot=[bigtri;xdot];%點集
edgemat=[1 2 xdot(1,:) xdot(2,:);...
    2 3 xdot(2,:) xdot(3,:);1 3 xdot(1,:) xdot(3,:)];%邊集，每個點包含2個點，4個座標值
trimat=[1 2 3];%三角集,每個三角包含3個點
temp_trimat=[1 2 3];
for j=4:N+3
    pointtemp=xdot(j,:);%迴圈每一個點
    deltemp=[];%初始化刪除temp_trimat的點
    temp_edgemat=[];%初始化臨時邊
    for k=1:size(temp_trimat,1)%迴圈每一個temp_trimat的三角形
        panduan=whereispoint(xdot(temp_trimat(k,1),:),...
            xdot(temp_trimat(k,2),:),xdot(temp_trimat(k,3),:),pointtemp);%判斷點在圓內0、圓外1、圓右側2
        switch panduan
            case 0
                %點在圓內
                %則該三角形不為Delaunay三角形
                temp_edge=maketempedge(temp_trimat(k,1),temp_trimat(k,2),temp_trimat(k,3),j,xdot);%把三條邊暫時存放於臨時邊矩陣
                temp_edgemat=[temp_edgemat;temp_edge];
                deltemp=[deltemp,k];
                ;
            case 1
                %點在圓外，pass
                ;
            case 2
                %點在圓右
                %則該三角形為Delaunay三角形，儲存到triangles
                trimat=[trimat;temp_trimat(k,:)];%新增到正式三角形中
                deltemp=[deltemp,k];
                %並在temp裡去除掉
                %別忘了把正式的邊也新增進去
                edgemat=[edgemat;makeedge(temp_trimat(k,1),temp_trimat(k,2),temp_trimat(k,3),xdot)];%遵循12,13,23的順序
                edgemat=unique(edgemat,'stable','rows');          
        end


    %三角迴圈結束    
    end

    %除去上述步驟中的臨時三角形
    temp_trimat(deltemp,:)=[];
    temp_trimat(~all(temp_trimat,2),:)=[];
    %對temp_edgemat去重複
    temp_edgemat=unique(temp_edgemat,'stable','rows');
    %將edge buffer中的邊與當前的點進行組合成若干三角形並儲存至temp triangles中
    temp_trimat=[temp_trimat;maketemptri(temp_edgemat,xdot,j)];
    k=k;


%點迴圈結束
end

%合併temptri
trimat=[trimat;temp_trimat];
edgemat=[edgemat;temp_edgemat];
%刪除大三角形
deltemp=[];
for j=1:size(trimat,1)
    if ismember(1,trimat(j,:))||ismember(2,trimat(j,:))||ismember(3,trimat(j,:))
        deltemp=[deltemp,j];
    end
end
trimat(deltemp,:)=[];

%凸包監測
%思路是先找出邊緣點（三角形只有1個或2個的），順便整出一個三角形相互關係圖，以後用。
%然後順時針，依次隔一個點連接出一條線段，如果這個和之前的線段相交，則不算；如果不交，則記錄出三角形
%更新完了以後，再監測一遍，直到沒有新的為止。


figure(2)
hold on
% plot(xdot(:,1),xdot(:,2),'ko')
for j=4:size(trimat,1)
    plot([xdot(trimat(j,1),1),xdot(trimat(j,2),1)],[xdot(trimat(j,1),2),xdot(trimat(j,2),2)],'k-')
    plot([xdot(trimat(j,1),1),xdot(trimat(j,3),1)],[xdot(trimat(j,1),2),xdot(trimat(j,3),2)],'k-')
    plot([xdot(trimat(j,3),1),xdot(trimat(j,2),1)],[xdot(trimat(j,3),2),xdot(trimat(j,2),2)],'k-')
end
hold off




%判斷點在三角形外接圓的哪個部分
function panduan=whereispoint(xy1,xy2,xy3,xy0)
%判斷點在三角形外接圓的哪個部分
[a,b,r2]=maketricenter(xy1,xy2,xy3);
x0=xy0(1);y0=xy0(2);
if a+sqrt(r2)<x0
    %x0在圓的右側
    panduan=2;
elseif (x0-a)^2+(y0-b)^2<r2
    %x0在圓內
    panduan=0;
else
    %在圓外
    panduan=1;
end
end

%做出三角形三點與內部1點之間的線段
function temp_edge=maketempedge(dot1,dot2,dot3,dot0,xdot)
%做出連線點與三角形之間的線
%每行包含2個點，4個座標值，共3行
%xy1和xy0組成線段
temp_edge=zeros(3,6);
if xdot(dot1,1)<xdot(dot0,1)
    temp_edge(1,:)=[dot1,dot0,xdot(dot1,:),xdot(dot0,:)];
elseif xdot(dot1,1)==xdot(dot0,1)
    if xdot(dot1,2)<xdot(dot0,2)
        temp_edge(1,:)=[dot1,dot0,xdot(dot1,:),xdot(dot0,:)];
    else
        temp_edge(1,:)=[dot0,dot1,xdot(dot0,:),xdot(dot1,:)];
    end
else
    temp_edge(1,:)=[dot0,dot1,xdot(dot0,:),xdot(dot1,:)];
end
%xy2和xy0組成線段
if xdot(dot2,1)<xdot(dot0,1)
    temp_edge(2,:)=[dot2,dot0,xdot(dot2,:),xdot(dot0,:)];
elseif xdot(dot2,1)==xdot(dot0,1)
    if xdot(dot2,2)<xdot(dot0,2)
        temp_edge(2,:)=[dot2,dot0,xdot(dot2,:),xdot(dot0,:)];
    else
        temp_edge(2,:)=[dot0,dot2,xdot(dot0,:),xdot(dot2,:)];
    end
else
    temp_edge(2,:)=[dot0,dot2,xdot(dot0,:),xdot(dot2,:)];
end
%xy3和xy0組成線段
if xdot(dot3,1)<xdot(dot0,1)
    temp_edge(3,:)=[dot3,dot0,xdot(dot3,:),xdot(dot0,:)];
elseif xdot(dot3,1)==xdot(dot0,1)
    if xdot(dot3,2)<xdot(dot0,2)
        temp_edge(3,:)=[dot3,dot0,xdot(dot3,:),xdot(dot0,:)];
    else
        temp_edge(3,:)=[dot0,dot3,xdot(dot0,:),xdot(dot3,:)];
    end
else
    temp_edge(3,:)=[dot0,dot3,xdot(dot0,:),xdot(dot3,:)];
end

end

%做出一些列固定點發散的線段外點組成的三角形
function temp_trimat=maketemptri(temp_edgemat,xdot,dot0)
%將edge buffer中的邊與當前的點進行組合成若干三角形
%temp_edgemat是新邊，x是中心點
%思路是計算各個邊對應角度，然後排序相連

A=temp_edgemat(:,1:2);
pointline=A(A~=dot0);
N=length(pointline);
pointaxe=xdot(pointline,:);
img_pointaxe=pointaxe(:,1)+1i*pointaxe(:,2);
d_img_pointaxe=img_pointaxe-xdot(dot0,1)-1i*xdot(dot0,2);
angle_d_img_pointaxe=angle(d_img_pointaxe);
[~,index]=sort(angle_d_img_pointaxe);
index=[index;index(1)];%排序，然後依次串起來
temp_trimat=zeros(N,3);
for j=1:N
    temp_trimat(j,:)=[pointline(index(j)),pointline(index(j+1)),dot0];
end


end

%將三個點構成3條邊
function edgemat=makeedge(dot1,dot2,dot3,xdot)
%將dot1 2 3這三個點構成三條邊
%每行包含2個點，4個座標值，共3行
edgemat=zeros(3,6);
%點12
if xdot(dot1,1)<xdot(dot2,1)
    edgemat(1,:)=[dot1,dot2,xdot(dot1,:),xdot(dot2,:)];
elseif xdot(dot1,1)==xdot(dot2,1)
    if xdot(dot1,2)<xdot(dot2,2)
        edgemat(1,:)=[dot1,dot2,xdot(dot1,:),xdot(dot2,:)];
    else
        edgemat(1,:)=[dot2,dot1,xdot(dot2,:),xdot(dot1,:)];
    end
else
    edgemat(1,:)=[dot2,dot1,xdot(dot2,:),xdot(dot1,:)];
end
%點13
if xdot(dot1,1)<xdot(dot3,1)
    edgemat(2,:)=[dot1,dot3,xdot(dot1,:),xdot(dot3,:)];
elseif xdot(dot1,1)==xdot(dot3,1)
    if xdot(dot1,2)<xdot(dot3,2)
        edgemat(2,:)=[dot1,dot3,xdot(dot1,:),xdot(dot3,:)];
    else
        edgemat(2,:)=[dot3,dot1,xdot(dot3,:),xdot(dot1,:)];
    end
else
    edgemat(2,:)=[dot3,dot1,xdot(dot3,:),xdot(dot1,:)];
end
%點23
if xdot(dot3,1)<xdot(dot2,1)
    edgemat(3,:)=[dot3,dot2,xdot(dot3,:),xdot(dot2,:)];
elseif xdot(dot3,1)==xdot(dot2,1)
    if xdot(dot3,2)<xdot(dot2,2)
        edgemat(3,:)=[dot3,dot2,xdot(dot3,:),xdot(dot2,:)];
    else
        edgemat(3,:)=[dot2,dot3,xdot(dot2,:),xdot(dot3,:)];
    end
else
    edgemat(3,:)=[dot2,dot3,xdot(dot2,:),xdot(dot3,:)];
end
% edgemat
end

%求三角形外接圓圓心
function [a,b,r2]=maketricenter(xy1,xy2,xy3)
x1=xy1(1);y1=xy1(2);
x2=xy2(1);y2=xy2(2);
x3=xy3(1);y3=xy3(2);
a=((y2-y1)*(y3*y3-y1*y1+x3*x3-x1*x1)-(y3-y1)*(y2*y2-y1*y1+x2*x2-x1*x1))/(2.0*((x3-x1)*(y2-y1)-(x2-x1)*(y3-y1)));
b=((x2-x1)*(x3*x3-x1*x1+y3*y3-y1*y1)-(x3-x1)*(x2*x2-x1*x1+y2*y2-y1*y1))/(2.0*((y3-y1)*(x2-x1)-(y2-y1)*(x3-x1)));
r2=(x1-a)*(x1-a)+(y1-b)*(y1-b);
end

首先是隨機建立一個xdot點矩陣，用來生成網格點。

然後對點進行先x後y的整理。這個要求是改進後演算法本身的要求，因為演算法要求點的序號要從左到右依次增大，這樣有助於節約相關點判斷的數量。

然後畫出最大包含的三角形，這裡要注意單純的增大座標並不意味著能把所有點包住，可能會有點落在斜邊外。建議先求出包含所有點矩形，然後在新增一個斜邊斜率約束的思考考慮。

之後生成包含大三角形的點集xdot，序號1 2 3對應大三角形，每行包含兩個座標值。邊集edgemat的每行包含2個點序號和4個座標值。三角集trimat的每行只包含三個點的序號值。

之後初始化temp_trimat=[1,2,3]，開始遍歷xdot的每一個點。再之後和上述虛擬碼中一樣，嚴格按照演算法要求走就行了。

接下來介紹一下程式中用到的內建函式
panduan=whereispoint(xy1,xy2,xy3,xy0)
是用來判斷點在三角形外接圓的哪個部分，具體實現是求出三角形點的圓心和半徑，然後利用簡單的幾何做判斷。

temp_edge=maketempedge(dot1,dot2,dot3,dot0,xdot)
是用來做出三角形dot1,dot2,dot3和內部一個點dot0之間的連線。這個函式看起來比較亂，但實際上只是想用判斷去整理一下邊的序號，最後結果看起來更整齊一些。

temp_trimat=maketemptri(temp_edgemat,xdot,dot0)
這個函式是做出一些列固定點發散的線段外點組成的三角形，對應之前演算法裡的小圖(d)。
形象的可以看下圖，即輸入為中間這些黑線，輸出為藍線包圍的這4個三角形。具體思路是求出各個線段的角度，然後按照角度排序依次相連。我試過用判斷是否交叉來做，但是出現了下圖這個反例所以又重新修改了程式碼。

edgemat=makeedge(dot1,dot2,dot3,xdot)
這個函式是將dot1,dot2,dot3三個邊生成三條邊，雖然沒什麼難度，但為了整齊，我還是寫了一大堆。我覺得這個在之後的優化中可以刪掉精簡一些。

實現的最終效果圖如下，這是1000點的圓形隨機點圖：

4。原演算法的一些問題和改進思路

原blog的演算法其實存在一個比較大的bug，就是邊緣連線並不能保證是一個凸型，也就是說邊緣的Delaunay三角剖分不完整。

這個現象在邊緣點比較少的時候比較突出，比如下圖就是隻設定了10個點的隨機圖

可以看到，至少左邊可以補齊2個三角形作為外邊三角形，使得圖形變成凸邊形，而且每個點都最大限度的得到了利用。

本文的演算法思路是獲取邊緣三角形資訊，生成邊緣點列表，按逆時針排序（或順時針）。然後隔一個點連一條線，判斷這條線是否和原圖形相交，如果不想交，則合法，生成新三角形。這麼做的缺點就是不能保證邊緣得到的圖形一定是Delaunay三角形，但是根據幾何判定應該說大部分情況都符合Delaunay三角形，只有像上圖這種凹陷點涉及到4個或以上的情況有概率不是。

這個bug的補救其準備在下一篇文章voronoi圖（泰森多邊形）的文章裡應用。