二维Delaunay（德洛内）三角网剖分的matlab实现

大家好，欢迎来到IT知识分享网。

二维Delaunay（德洛内）三角网的matlab实现

1.Delaunay三角网的概述

德洛内(Delaunay)三角网的定义: 它是一系列相连的但不重叠的三角形的集合, 而且这些三角形的外接圆不包含这个面域的其他任何点。它具有两个特征：

(1) 每个德洛内(Delaunay) 三角形的外接圆不包含面内的其他任何点, 称之为德洛内(Delaunay) 三角网的空外接圆性质, 这个特征已经作为创建德洛内(Delaunay) 三角网的一项判别标准;

(2) 它的另一个性质最大最小角性质: 每两个相邻的三角形构成的凸四边形的对角线，在相互交换后，六个内角的最小角不再增大。

Delaunay 三角网的优点是结构良好, 数据结构简单, 数据冗余度小, 存储效率高, 与不规则的地面特征和谐一致,可以表示线性特征和迭加任意形状的区域边界, 易于更新,可适应各种分布密度的数据等; 它的局限性是, 算法实现比较复杂和困难, 但现在已经有了较多成熟的实现算法。

其中Voronoi图（泰森多边形）的画法：https://blog.csdn.net/weixin_42943114/article/details/82319332
彩色Voronoi图的画法：https://blog.csdn.net/weixin_42943114/article/details/82461228

2.Delaunay三角网的算法

算法的伪代码如下：

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=1: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的每一个点。再之后和上述伪代码中一样，严格按照算法要求走就行了。

4.原算法的一些问题和改进思路

原blog的算法其实存在一个比较大的bug，就是边缘连线并不能保证是一个凸型，也就是说边缘的Delaunay三角剖分不完整。

可以看到，至少左边可以补齐2个三角形作为外边三角形，使得图形变成凸边形，而且每个点都最大限度的得到了利用。

本文的算法思路是获取边缘三角形信息，生成边缘点列表，按逆时针排序（或顺时针）。然后隔一个点连一条线，判断这条线是否和原图形相交，如果不想交，则合法，生成新三角形。这么做的缺点就是不能保证边缘得到的图形一定是Delaunay三角形，但是根据几何判定应该说大部分情况都符合Delaunay三角形，只有像上图这种凹陷点涉及到4个或以上的情况有概率不是。

这个bug的补救其准备在下一篇文章voronoi图（泰森多边形）的文章里应用。
更新：由于在下一篇介绍对阅读造成了很大的困扰，我决定在这里介绍一下：

4.1原算法的问题

%画出最大包含的三角形 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]; bigtri=10*bigtri; 

这里最后加上bigtri=10*bigtri;来扩大三角形，可以减小这个发生的概率。

4.2凸包检测

由于Delaunay三角形构成的图案一定是凸多边形，这里为了避免会出现凹多边形的结果出现，将外缘凹陷的边缘进行补齐，来完善Delaunay三角形。这部分代码如下：

%凸包监测 %思路是先找出边缘点（三角形只有1个或2个的），顺便整出一个三角形相互关系图，以后用。 %然后顺时针，依次隔一个点连接出一条线段，如果这个和之前的线段相交，则不算；如果不交，则记录出三角形 %更新完了以后，再监测一遍，直到没有新的为止。 t_w=0; while t_w==0 [~,border_point,~]=makebordertri(trimat); border_point=[border_point;border_point(1,:)]; temp_edgemat=[]; temp_trimat=[]; for j=1:size(border_point,1)-1 tempboderedge=[border_point(j,1),border_point(j+1,2)]; tempboderdot=border_point(j,2); %寻找带tempboderdot的所有边 tempdotex=edgemat(logical(sum(edgemat==tempboderdot,2)),:); %删除相邻边 tempdotex(ismember(tempdotex,[tempboderdot,tempboderedge(1)],'rows'),:)=[]; tempdotex(ismember(tempdotex,[tempboderedge(1),tempboderdot],'rows'),:)=[]; tempdotex(ismember(tempdotex,[tempboderdot,tempboderedge(2)],'rows'),:)=[]; tempdotex(ismember(tempdotex,[tempboderedge(2),tempboderdot],'rows'),:)=[]; %检测tempdotex是否为空，如果是证明不用相连 t_N=size(tempdotex,1); t_t=0; if t_N>0 %依次检测是否相交，只要有一个相交就不算；如果都不想交，则相连 for k=1:t_N if tempdotex(k,1)==tempboderdot t_xdotno4=tempdotex(k,2); else t_xdotno4=tempdotex(k,1); end tt_xdotno4=xdot(t_xdotno4,:)-xdot(tempboderdot,:); xdotno4=xdot(tempboderdot,:)+tt_xdotno4/sqrt(sum(tt_xdotno4.^2))*(sqrt((xmax-xmin)^2+(ymax-ymin)^2)); panduan=crossornot(xdot(tempboderedge(1),:),xdot(tempboderedge(2),:),xdot(tempboderdot,:),xdotno4); if panduan==1 t_t=t_t+1; break end end %t_t大于0说明有相交的线,略过 if t_t==0 temp_edgemat=[temp_edgemat;tempboderedge]; temp_trimat=[temp_trimat;[tempboderedge,tempboderdot]]; break end end end trimat=[trimat;temp_trimat]; edgemat=[edgemat;temp_edgemat]; %删除重复的三角形 trimat=sort(trimat,2); trimat=unique(trimat,'stable','rows'); if j==size(border_point,1)-1 t_w=1; end end 

4.3改进效果