MATLAB小技巧（2）Delaunay三角剖分算法

前言

MATLAB进行图像处理相关的学习是非常友好的，可以从零开始，对基础的图像处理都已经有了封装好的许多可直接调用的函数，这个系列文章的话主要就是介绍一些大家在MATLAB中常用一些概念函数进行例程演示！

Delaunay三角剖分的概念：给定平面一组离散点，构造一个三角剖分，使得每个三角形的外接圆内不包含第四个点。在三维情形，我们构造四面体剖分，使得每个四面体的外接球不包含第五个点。Delaunay三角剖分使得最小内角最大化，从而提高数值模拟的稳定性，因此成为网格生成算法的根基。

Delaunay三角剖分算法于1970年代被发明出来，由此诞生了“计算几何”这一学科。从此，计算几何领域的主旋律一直是Delaunay三角剖分。经过数十年的发展，Delaunay三角剖分算法日趋成熟，目前自动四面体网格生成问题已经被解决，并且在工业界得到广泛应用。

计算几何领域的翘楚- Herbert Edelsbrunner 基于Delaunay三角剖分技巧开创了知名的Geomagic公司，主要应用于点云重建。Geomagic于2013年被3D Systems收购，成为3D打印工业不可或缺的软件工具。

一. MATLAB仿真一

clc;
clear;
close all;

rand('state',0);
x = rand(1,10);
y = rand(1,10);
TRI = delaunay(x,y);
figure(1),triplot(TRI,x,y)
axis([0 1 0 1]);
hold on;
plot(x,y,'or');
hold off

二. 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

三. MATLAB仿真三

clc; clear; close all;
rand('state',0);
x = rand(1,15);
y = rand(1,15);%随机数
TRI = delaunay(x,y);%三角剖分
triplot(TRI,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;
for i = 1 : size(TRI, 1)%遍历三角形个数
    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

四. 小结

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

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

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

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

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

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

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

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

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

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

三角剖分不仅可以用于作图，在多目标跟踪领域也有相应的应用，例如远距离的多个弱小目标，由于距离较远，在足够的采样精度下，经过三角剖分定位可以大致确定不同目标的相对位置，从而实现多目标的标定跟踪。每天学一个MATLAB小知识，大家一起来学习进步阿！