当前位置: 代码迷 >> 综合 >> js-提取离散点边界(二维+三维)MATLAB Delaunay
  详细解决方案

js-提取离散点边界(二维+三维)MATLAB Delaunay

热度:76   发布时间:2023-12-23 06:13:50.0

 先看三维效果图:

在这里插入图片描述

 

二维算例,参考
https://blog.csdn.net/qq_38991255/article/details/81223437
借鉴该实例,修改后的三维边界提取代码见文后

该算法的总体思路如下:
1、利用 delaunay 函数,对所有数据点进行 Delaunay 三角剖分处理,delaunay 函数的返回值是一个 N * 3 的矩阵,其中 N 为剖分出的三角形个数,3 为每个三角形的三个端点的序号。
2、根据 triangles 矩阵,提取出所有 delaunay 三角剖分时所连接的边,依次扫描 triangles 矩阵的每一行,将 delaunay 三角剖分时所连接的边添加到一个新的矩阵中,最后构成一个 M * 2 的矩阵,其中 M 是一共所连接的边的条数。
3、显然,最小凸多边形上的边应该仅在以上矩阵中出现一次,因此,将以上矩阵中那些出现次数超过一次的边全部去掉,最后保留的便是最小凸多边形的边。
4、根据最小凸多边形的边,很容易得到构成最小凸多边形的结点的顺序,从而解决问题。
输入参数 points 是一个 2 * P 矩阵, P 为数据点的个数,第一行是这些数据点对应的 x 坐标,第二行是对应的 y 坐标;输出参数 polygon 是一个 2 * Q 矩阵, Q 为凸多边形的顶点个数(首尾相连),第一行是这些顶点对应的 x 坐标,第二行是对应的 y 坐标。代码实现如下:

 function polygon= minimal_convex_polygon(points)% 进行 delaunay 三角剖分,将所有连接了的边保存在矩阵 lines 中triangles = sort(delaunay(points(1, :), points(2, :)), 2);lines = zeros(size(triangles, 1) * 3, 2);for i = 1:size(triangles, 1)      % 对三个点进行组合,转换成三条线存放在lines中 lines(3 * i - 2,:) = [triangles(i, 1), triangles(i, 2)];lines(3 * i - 1,:) = [triangles(i, 1), triangles(i, 3)];lines(3 * i,:) = [triangles(i, 2), triangles(i, 3)];end% 去掉 lines 中出现次数超过一次的边[~, IA] = unique(lines, 'rows');          %IA表示,取出lines中只出现一次行的索引,波浪线~不是忽略该参数lines = setdiff(lines(IA, :), lines(setdiff(1:size(lines, 1), IA), :), 'rows');% % 只保留出现一次的边。C=setdiff(A,B)函数返回在向量A中却不在向量B中的元素,并且C中不包含重复元素,并且从小到大排序% % 跟踪 lines 中的数据点,将凸多边形的顶点编号保存在 seqs 中seqs = zeros(size(lines, 1) + 1,1);seqs(1:2) = lines(1, :);lines(1, :) = [];for i = 3:size(seqs)pos = find(lines == seqs(i - 1));row = rem(pos - 1, size(lines, 1)) + 1;col = ceil(pos / size(lines, 1));seqs(i) = lines(row, 3 - col);lines(row, :) = [];end% 根据 seqs , 得到凸多边形顶点坐标polygon = points(:, seqs);
end

 

定义了实现函数,下面进行调用:

P = gallery('uniformdata',[30 2],0);
Pp=P';% 调用polygon
plot(Pp(1,:),Pp(2,:), '*r', 'LineWidth', 4);        % Pp第一行为x坐标,第二行为y坐标
polygon = minimal_convex_polygon(Pp);
hold on;
plot(polygon(1, :), polygon(2, :), 'r','LineWidth', 5);%显示三角剖分
DT = delaunay(P);
triplot(DT,P(:,1),P(:,2))
axis equal

二维图呈现:

在这里插入图片描述

 

基于以上思路,修改的三维离散点边界提取算法如下:


P = gallery('uniformdata',30,3,5);
points=P;
triangles = sort(delaunay(points(:,1), points(:,2),points(:,3)), 2);   %生成四面体四个点的索引
faces = zeros(size(triangles, 1) * 4, 3);        %三维为面 for i = 1:size(triangles, 1)        %对四个点进行组合,转换成四个面存放在faces中  faces(4 * i-3,:)=[triangles(i, 1), triangles(i, 2), triangles(i, 3)];faces(4 * i - 2,:) = [triangles(i, 1), triangles(i, 2), triangles(i, 4)];faces(4 * i - 1,:) = [triangles(i, 1), triangles(i, 3), triangles(i, 4)];faces(4 * i,:) = [triangles(i, 2), triangles(i, 3), triangles(i, 4)];end% 去掉 faces中出现次数超过一次的面
[~, IA] = unique(faces, 'rows');      %IA表示,取出faces中只出现一次行的索引(注意,其中包含重复行索引,但只会取一次)
faces = setdiff(faces(IA, :), faces(setdiff(1:size(faces, 1), IA), :),'rows');  %只保留出现一次的面。C=setdiff(A,B)函数返回在向量A中却不在向量B中的元素,并且C中不包含重复元素,并且从小到大排序
trisurf(faces,P(:,1),P(:,2),P(:,3),'FaceAlpha',0.5);hold on
plot3(P(:,1),P(:,2),P(:,3),'.','MarkerSize',10);   
axis equal

结果如下:
在这里插入图片描述

 

  相关解决方案