创建和编辑德劳内三角剖分
此示例说明如何使用 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
示例二:创建和绘制三维德劳内三角剖分
此示例说明如何计算和绘制三维德劳内三角剖分。
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])
要显示大型四面体网格,请使用 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')
示例六:创建地图的有约束德劳内三角剖分
加载美国本土的周长地图。构造一个受约束的代表该多边形的德劳内三角剖分。该三角剖分跨受点集的凸包约束的一个域。过滤出多边形域内的三角形并对其绘图。注意:该数据集包含重复的数据点;即两个或更多个数据点的位置相同。重复的点将被拒绝,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')
示例七:基于点云重新构造曲线
下面的示例重点介绍如何使用德劳内三角剖分从点云重新构造多边形边界。重新构造操作基于精心设计的 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
示例八:计算多边形域的近似中轴
此示例说明如何使用受约束的德劳内三角剖分创建一个多边形域的近似中轴。多边形的中轴由多边形内部最大圆的中心的轨迹定义。
构造一个域边界上采样点的有约束德劳内三角剖分。
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
示例九:按照修改的边界变换二维网格
此示例说明如何变换二维域的网格以适应对域边界的修改。
步骤 1:加载数据。要变换的网格由 trife
、xfe
和 yfe
定义,它是面-顶点格式的三角剖分。
load trimesh2d clf triplot(trife,xfe,yfe) axis equal axis([-10 310 -10 310]) axis equal xlabel('Initial Mesh','FontWeight','b')
步骤 2:构造背景三角剖分 - 由代表网格边界的点集构成的受约束的德劳内三角剖分。对于网格的每个顶点,计算用于定义与其背景三角剖分相关的位置的描述符。该描述符表示该封闭三角形以及与该三角形相关的重心坐标。
dt = delaunayTriangulation(x,y,Constraints); clf triplot(dt) axis equal axis([-10 310 -10 310]) axis equal xlabel('Background Triangulation','FontWeight','b')
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')
步骤 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')