MATLAB实现三角剖分（Delaunay）算法

最新推荐文章于 2023-07-17 09:21:44 发布

知行流浪

最新推荐文章于 2023-07-17 09:21:44 发布

阅读量2w

点赞数 22

分类专栏： MATLAB MATLAB菜鸟之路文章标签： matlab 三角剖分 Delaunay

本文链接：https://blog.csdn.net/zengxiantao1994/article/details/77407750

版权

MATLAB菜鸟之路同时被 2 个专栏收录

15 篇文章 67 订阅

订阅专栏

MATLAB

13 篇文章 29 订阅

订阅专栏

三角剖分定义

【定义】三角剖分：假设V是二维实数域上的有限点集，边e是由点集中的点作为端点构成的封闭线段，E为e的集合。那么该点集V的一个三角剖分T = (V,E)是一个平面图G，该平面图满足条件：

1、除了端点，平面图中的边不包含点集中的任何点。

2、没有相交边。// 边和边没有交叉点。

3、平面图中所有的面都是三角面，且所有三角面的合集是散点集V的凸包。

// 凸包的概念：用不严谨的话来讲，给定二维平面上的点集，凸包就是将最外层的点连接起来构成的凸多边型，它能包含点集中所有的点。

Delaunay三角剖分的定义

在实际中运用的最多的三角剖分是Delaunay三角剖分，它是一种特殊的三角剖分。先从Delaunay边说起：

【定义】Delaunay边：假设E中的一条边e (两个端点为a，b)，e若满足下列条件，则称之为Delaunay边：存在一个圆经过a，b两点，圆内（注意是圆内，圆上最多三点共圆）不含点集V中任何其他的点，这一特性又称空圆特性。

【定义】Delaunay三角剖分：如果点集V的一个三角剖分T只包含Delaunay边，那么该三角剖分称为Delaunay三角剖分。

Delaunay三角剖分的准则

要满足Delaunay三角剖分的定义，必须符合两个重要的准则：

1、空圆特性：Delaunay三角网是唯一的（任意四点不能共圆），在Delaunay三角形网中任一三角形的外接圆范围内不会有其它点存在。如下图所示：

2、最大化最小角特性：在散点集可能形成的三角剖分中，Delaunay三角剖分所形成的三角形的最小角最大。从这个意义上讲，Delaunay三角网是"最接近于规则化的"三角网。具体的说是在两个相邻的三角形构成凸四边形的对角线，在相互交换后，两个内角的最小角不再增大。如下图所示：

Delaunay三角剖分的特性

以下是Delaunay剖分所具备的优异特性：

1、最接近：以最接近的三点形成三角形，且各线段（三角行的边）皆不相交。

2、唯一性：不论从区域何处开始构建，最终都将得到一致的结果。

3、最优性：任意两个相邻三角形构成的凸四边形的对角线如果可以互换，那么两个三角形六个内角中最小角度不会变化。

4、最规则：如果将三角网中的每个三角形的最小角进行升序排列，则Delaunay三角网的排列得到的数值最大。

5、区域性：新增、删除、移动某一个顶点只会影响邻近的三角形。

6、具有凸边形的外壳：三角网最外层的边界形成一个凸多边形的外壳。

局部最优化处理

理论上为了构造Delaunay三角网，Lawson提出的局部优化过程LOP(Local Optimization Procedure)，一般三角网经过LOP处理，即可确保为Delaunay三角网，其基本做法如下所示：

1、将两个具有共同边的三角形合成一个多边形。

2、以最大空圆准则作检查，看其第四个顶点是否在三角形的外接圆内。

3、如果在，修正对角线即将对角线对调，即完成局部优化过程的处理。

LOP处理过程如下图所示：

Delaunay剖分的算法

Delaunay剖分是一种三角剖分的标准，实现它有多种算法。最常用的是逐点插入法。

Lawson算法的基本步骤是：

1、构造一个超级三角形，包含所有散点，放入三角形链表。

2、将点集中的散点依次插入，在三角形链表中找出其外接圆包含插入点的三角形（称为该点的影响三角形），删除影响三角形的公共边，将插入点同影响三角形的全部顶点连接起来，从而完成一个点在Delaunay三角形链表中的插入。

3、根据优化准则对局部新形成的三角形进行优化。将形成的三角形放入Delaunay三角形链表。

4、循环执行上述第2步，直到所有散点插入完毕。

这一算法的关键的第2步图示如下：

伪代码表述

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

MATLAB实现

clc;
clear;
close all;

rand('state', 0);
node = 8;
x = rand(1,node);
y = rand(1,node);
% delaunay是MATLAB中三角剖分的函数，返回的TRI是三角形的矩阵
% TRI的每一行表示三角形的三个点
TRI = delaunay(x,y);

% 绘图
figure; 
xmin = min(x(:)); xmax = max(x(:));
ymin = min(y(:)); ymax = max(y(:));
xl = xmax - xmin; yl = ymax - ymin;
axis([xmin-xl*0.1, xmax+xl*0.1,...
    ymin-yl*0.1, ymax+yl*0.1]);
hold on;

n = size(TRI, 1);
for i = 1 : n
    t1 = TRI(i, :);
    for j = 1 : length(t1)-1
        xt = [x(t1(j)) x(t1(j+1))];
        yt = [y(t1(j)) y(t1(j+1))];
        plot(xt, yt, 'k-', 'LineWidth', 2);
        pause(0.1);
    end
    xt = [x(t1(end)) x(t1(1))];
    yt = [y(t1(end)) y(t1(1))];
    plot(xt, yt, 'k-', 'LineWidth', 2);
    pause(0.1);
end

W = zeros(node);
for i = 1 : n
    for j = 1 : length(TRI(i, :))-1
        W(TRI(i, j), TRI(i, j+1)) = 1;
        W(TRI(i, j+1), TRI(i, j)) = 1;
    end
    W(TRI(i, end), TRI(i, 1)) = 1;
    W(TRI(i, 1), TRI(i, end)) = 1;
end
for i = 1 : node
    for j = 1 : node
        if ~W(i, j)
            W(i, j) = 10000;
        end
    end
end