Main Content

创建和编辑德劳内三角剖分

此示例说明如何使用 delaunayTriangulation 类创建、编辑和查询德劳内三角剖分。德劳内三角剖分是科学计算中使用最广泛的三角剖分。与三角剖分关联的属性提供了解决各种几何问题的基础。此外,还说明了如何构造受约束的德劳内三角剖分,同时提供了一个涵盖中轴计算和网格变换的应用程序。

示例一:创建和绘制二维德劳内三角剖分

此示例说明如何计算二维德劳内三角剖分,然后同时对三角剖分以及顶点和三角标签绘图。

rng default
x = rand(10,1);
y = rand(10,1);
dt = delaunayTriangulation(x,y)
dt = 
  delaunayTriangulation with properties:

              Points: [10x2 double]
    ConnectivityList: [11x3 double]
         Constraints: []

triplot(dt)

在图上显示顶点和三角标签。

hold on
vxlabels = arrayfun(@(n) {sprintf('P%d', n)}, (1:10)');
Hpl = text(x,y,vxlabels,'FontWeight','bold','HorizontalAlignment',...
   'center','BackgroundColor','none');
ic = incenter(dt);
numtri = size(dt,1);
trilabels = arrayfun(@(x) {sprintf('T%d',x)}, (1:numtri)');
Htl = text(ic(:,1),ic(:,2),trilabels,'FontWeight','bold', ...
   'HorizontalAlignment','center','Color','blue');
hold off

Figure contains an axes object. The axes object contains 22 objects of type line, text.

示例二:创建和绘制三维德劳内三角剖分

此示例说明如何计算和绘制三维德劳内三角剖分。

rng default
X = rand(10,3);
dt = delaunayTriangulation(X)
dt = 
  delaunayTriangulation with properties:

              Points: [10x3 double]
    ConnectivityList: [18x4 double]
         Constraints: []

tetramesh(dt)
view([10 20])

Figure contains an axes object. The axes object contains 18 objects of type patch.

要显示大型四面体网格,请使用 convexHull 方法计算边界三角剖分,并使用 trisurf 对其绘图。例如:

triboundary = convexHull(dt);
trisurf(triboundary, X(:,1), X(:,2), X(:,3),'FaceColor','cyan')

示例三:访问三角剖分数据结构体

可通过两种方式访问三角剖分数据结构。一种方式是通过 Triangulation 属性,另一种方式是使用索引。

基于 10 个随机点创建一个二维德劳内三角剖分。

rng default
X = rand(10,2);
dt = delaunayTriangulation(X)
dt = 
  delaunayTriangulation with properties:

              Points: [10x2 double]
    ConnectivityList: [11x3 double]
         Constraints: []

访问三角剖分数据结构体的一种方式是使用 ConnectivityList 属性。

dt.ConnectivityList
ans = 11×3

     2     8     5
     7     6     1
     3     7     8
     8     7     5
     3     8     2
     6     7     3
     7     4     5
     5     9     2
     2     9    10
     5     4     9
      ⋮

使用索引是查询三角剖分的一种简便方式。语法为 dt(i,j),其中 j 是第 i 个三角形的第 j 个顶点。需遵循标准索引规则。

使用索引查询三角剖分数据结构体。

dt(:,:)
ans = 11×3

     2     8     5
     7     6     1
     3     7     8
     8     7     5
     3     8     2
     6     7     3
     7     4     5
     5     9     2
     2     9    10
     5     4     9
      ⋮

第二个三角形为:

dt(2,:)
ans = 1×3

     7     6     1

第二个三角形的第三个顶点为:

dt(2,3)
ans = 1

前三个三角形为:

dt(1:3,:)
ans = 3×3

     2     8     5
     7     6     1
     3     7     8

示例四:编辑德劳内三角剖分以插入或删除点

下面的示例说明了如何使用基于索引的下标插入或删除点。相对于从头重新创建新的 delaunayTriangulation,编辑 delaunayTriangulation 以进行少量修改的效率更高,数据集较大时尤其如此。

基于单位正方形内的 10 个随机点构造一个德劳内三角剖分。

rng default
x = rand(10,1);
y = rand(10,1);
dt = delaunayTriangulation(x,y)
dt = 
  delaunayTriangulation with properties:

              Points: [10x2 double]
    ConnectivityList: [11x3 double]
         Constraints: []

插入 5 个其他随机点。

dt.Points(end+(1:5),:) = rand(5,2)
dt = 
  delaunayTriangulation with properties:

              Points: [15x2 double]
    ConnectivityList: [20x3 double]
         Constraints: []

替换第五个点。

dt.Points(5,:) = [0 0]
dt = 
  delaunayTriangulation with properties:

              Points: [15x2 double]
    ConnectivityList: [20x3 double]
         Constraints: []

删除第四个点。

dt.Points(4,:) = []
dt = 
  delaunayTriangulation with properties:

              Points: [14x2 double]
    ConnectivityList: [18x3 double]
         Constraints: []

示例五:创建有约束德劳内三角剖分

此示例说明如何创建一个有约束德劳内三角剖分,并说明约束的影响。

创建并绘制一个有约束德劳内三角剖分。

X = [0 0; 16 0; 16 2; 2 2; 2 3; 8 3; 8 5; 0 5];
C = [1 2; 2 3; 3 4; 4 5; 5 6; 6 7; 7 8; 8 1];
dt = delaunayTriangulation(X,C);
subplot(2,1,1)
triplot(dt)
axis([-1 17 -1 6])
xlabel('Constrained Delaunay triangulation','FontWeight','b')

用红色绘制约束边。

hold on
plot(X(C'),X(C'+size(X,1)),'-r','LineWidth',2)
hold off

现在删除约束并绘制无约束德劳内三角剖分。

dt.Constraints = [];
subplot(2,1,2)
triplot(dt)
axis([-1 17 -1 6])
xlabel('Unconstrained Delaunay triangulation','FontWeight','b')

Figure contains 2 axes objects. Axes object 1 with xlabel Constrained Delaunay triangulation contains 9 objects of type line. Axes object 2 with xlabel Unconstrained Delaunay triangulation contains an object of type line.

示例六:创建地图的有约束德劳内三角剖分

加载美国本土的周长地图。构造一个受约束的代表该多边形的德劳内三角剖分。该三角剖分跨受点集的凸包约束的一个域。过滤出多边形域内的三角形并对其绘图。注意:该数据集包含重复的数据点;即两个或更多个数据点的位置相同。重复的点将被拒绝,delaunayTriangulation 会相应地重新设置约束的格式。

clf
load usapolygon

在构成多边形边界的两个连续点之间定义一个边约束,并创建德劳内三角剖分。

nump = numel(uslon);
C = [(1:(nump-1))' (2:nump)'; nump 1];
dt = delaunayTriangulation(uslon,uslat,C);
Warning: Duplicate data points have been detected and removed.
 The Triangulation indices and constraints are defined with respect to the unique set of points in delaunayTriangulation.
Warning: Intersecting edge constraints have been split, this may have added new points into the triangulation.
io = isInterior(dt);
patch('Faces',dt(io,:),'Vertices',dt.Points,'FaceColor','r')
axis equal
axis([-130 -60 20 55])
xlabel('Constrained Delaunay Triangulation of usapolygon','FontWeight','b')

Figure contains an axes object. The axes object with xlabel Constrained Delaunay Triangulation of usapolygon contains an object of type patch.

示例七:基于点云重新构造曲线

下面的示例重点介绍如何使用德劳内三角剖分从点云重新构造多边形边界。重新构造操作基于精心设计的 Crust 算法。

参考资料:N. Amenta、M. Bern 和 D. Eppstein。The crust and the beta-skeleton: combinatorial curve reconstruction.Graphical Models and Image Processing, 60:125-135, 1998.

创建一组表示点云的点。

numpts = 192;
t = linspace( -pi, pi, numpts+1 )';
t(end) = [];
r = 0.1 + 5*sqrt( cos( 6*t ).^2 + (0.7).^2 );
x = r.*cos(t);
y = r.*sin(t);
ri = randperm(numpts);
x = x(ri);
y = y(ri);

构造该点集的德劳内三角剖分。

dt = delaunayTriangulation(x,y);
tri = dt(:,:);

将沃罗诺伊顶点的位置插入现有三角剖分中。

V = voronoiDiagram(dt);

删除无限顶点并用 unique 筛除重复点。

V(1,:) = [];
numv = size(V,1);
dt.Points(end+(1:numv),:) = unique(V,'rows');

连接成对采样点的德劳内边表示边界。

delEdges = edges(dt);
validx = delEdges(:,1) <= numpts;
validy = delEdges(:,2) <= numpts;
boundaryEdges = delEdges((validx & validy),:)';
xb = x(boundaryEdges);
yb = y(boundaryEdges);
clf
triplot(tri,x,y)
axis equal
hold on
plot(x,y,'*r')
plot(xb,yb,'-r')
xlabel('Curve reconstruction from point cloud','FontWeight','b')
hold off

Figure contains an axes object. The axes object with xlabel Curve reconstruction from point cloud contains 194 objects of type line.

示例八:计算多边形域的近似中轴

此示例说明如何使用受约束的德劳内三角剖分创建一个多边形域的近似中轴。多边形的中轴由多边形内部最大圆的中心的轨迹定义。

构造一个域边界上采样点的有约束德劳内三角剖分。

load trimesh2d
dt = delaunayTriangulation(x,y,Constraints);
inside = isInterior(dt);

构造一个三角剖分来表示域中的三角形。

tr = triangulation(dt(inside,:),dt.Points);

构造一组联接相邻三角形外心的边。附加逻辑构造一个包含这样的边的唯一组。

numt = size(tr,1);
T = (1:numt)';
neigh = neighbors(tr);
cc = circumcenter(tr);
xcc = cc(:,1);
ycc = cc(:,2);
idx1 = T < neigh(:,1);
idx2 = T < neigh(:,2);
idx3 = T < neigh(:,3);
neigh = [T(idx1) neigh(idx1,1); T(idx2) neigh(idx2,2); T(idx3) neigh(idx3,3)]';

用绿色绘制域中的三角形,用蓝色绘制域边界,用红色绘制中轴。

clf
triplot(tr,'g')
hold on
plot(xcc(neigh), ycc(neigh), '-r','LineWidth',1.5)
axis([-10 310 -10 310])
axis equal
plot(x(Constraints'),y(Constraints'),'-b','LineWidth',1.5)
xlabel('Medial Axis of Polygonal Domain','FontWeight','b')
hold off

Figure contains an axes object. The axes object with xlabel Medial Axis of Polygonal Domain contains 364 objects of type line.

示例九:按照修改的边界变换二维网格

此示例说明如何变换二维域的网格以适应对域边界的修改。

步骤 1:加载数据。要变换的网格由 trifexfeyfe 定义,它是面-顶点格式的三角剖分。

load trimesh2d
clf
triplot(trife,xfe,yfe)
axis equal
axis([-10 310 -10 310])
axis equal
xlabel('Initial Mesh','FontWeight','b')

Figure contains an axes object. The axes object with xlabel Initial Mesh contains an object of type line.

步骤 2:构造背景三角剖分 - 由代表网格边界的点集构成的受约束的德劳内三角剖分。对于网格的每个顶点,计算用于定义与其背景三角剖分相关的位置的描述符。该描述符表示该封闭三角形以及与该三角形相关的重心坐标。

dt = delaunayTriangulation(x,y,Constraints);
clf
triplot(dt)
axis equal
axis([-10 310 -10 310])
axis equal
xlabel('Background Triangulation','FontWeight','b')

Figure contains an axes object. The axes object with xlabel Background Triangulation contains an object of type line.

descriptors.tri = pointLocation(dt,xfe,yfe);
descriptors.baryCoords = cartesianToBarycentric(dt,descriptors.tri,[xfe yfe]);

步骤 3:编辑背景三角剖分以将需要的修改融入域边界中。

cc1 = [210 90];
circ1 = (143:180)';
x(circ1) = (x(circ1)-cc1(1))*0.6 + cc1(1);
y(circ1) = (y(circ1)-cc1(2))*0.6 + cc1(2);
tr = triangulation(dt(:,:),x,y);
clf
triplot(tr)
axis([-10 310 -10 310])
axis equal
xlabel('Edited Background Triangulation - Hole Size Reduced','FontWeight','b')

Figure contains an axes object. The axes object with xlabel Edited Background Triangulation - Hole Size Reduced contains an object of type line.

步骤 4:使用变形的背景三角剖分作为计算的基础,将该描述符重新转换为笛卡尔坐标。

Xnew = barycentricToCartesian(tr,descriptors.tri,descriptors.baryCoords);
tr = triangulation(trife,Xnew);
clf
triplot(tr)
axis([-10 310 -10 310])
axis equal
xlabel('Morphed Mesh','FontWeight','b')

Figure contains an axes object. The axes object with xlabel Morphed Mesh contains an object of type line.