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

三角剖分定义

        【定义】三角剖分：假设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);endW = 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;endfor i = 1 : node    for j = 1 : node        if ~W(i, j)            W(i, j) = 10000;        end    endend