0

我有一组三角形。我正在寻找一种方法来找到这些三角形的所有组合，这些三角形在连接在一起时构成凸包。凸包应该是空的，即。仅在边缘上的凸包内没有点。并且只有共享一条边的三角形可以连接在一起，即。工会没有缺口。

示例：以下几点给出了 12 个三角形（Delaunay 三角剖分）。

xy = [3.3735 0.7889; -0.1072 -3.4814; -3.9732 4.1955; -5 5; 5 5; 5 -5; -5 -5];
DT = delaunayTriangulation(xy);
triplot(DT);
%The coordinates for each triangle -- each row is a triangle.
TRIX = reshape(DT.Points(DT.ConnectivityList, 1), size(DT.ConnectivityList));
TRIY = reshape(DT.Points(DT.ConnectivityList, 2), size(DT.ConnectivityList));

12个三角形

我正在寻找最大的凸包，所以凸包应该包含尽可能多的三角形。但如果我拥有所有可能的组合，我可以轻松过滤掉三角形较少的组合。在上面的示例中，我最终应该得到这六个凸包：

6个凸包

我猜我应该使用每个三角形最多有三个相邻的三角形（每边一个）。而且我应该检查相交点的角度总和是否小于或等于 180 度。这将确保联合是凸的——见下图。（如果几个三角形形成一个完整的圆，角度也可能正好是 360 度）。

5个凸包

三角角：

% Vectors connecting points
diffx = diff([TRIX TRIX(:,1)], [], 2); diffy = diff([TRIY TRIY(:,1)], [], 2); diffxy = [diffx(:) diffy(:)];
% Norm
normxy = reshape( arrayfun(@(row) norm(diffxy(row,:)), 1:size(diffxy,1)), size(DT.ConnectivityList));
nominator = repmat(sum(normxy.^2, 2), 1, 3) - 2*normxy.^2;
denominator = 2 * repmat(prod(normxy, 2), 1, 3)./normxy;
% Angle
tri_angle = acosd(nominator./denominator);
tri_angle = circshift(tri_angle, [0 -1]); % Shift so the angles match the correct point.

我重新格式化信息，使行是点，列是三角形：

n_tri = size(TRIX,1); % Number of triangles
% Adjacency matrix connecting points (rows) with triangles (columns).
adj_points = zeros(size(xy,1), n_tri);
adj_angle = NaN(size(adj_points));
for point =1:size(xy,1)
    idx = find(DT.ConnectivityList == point);
    [a_tri, ~] = ind2sub(size(DT.ConnectivityList), idx);
    adj_points(point,a_tri) = 1;
    adj_angle(point,a_tri) = tri_angle(idx);
end

我遍历所有边缘并计算边缘两侧的角度（edges angles）。通过这种方式，我能够找到形成凸集 ( adj_convex) 的三角形对：

DT_edges = edges(DT); % All edges in the Delaunay triangulation
% Adjacency matrix connecting edges (rows) with triangles (columns).
adj_edge = logical(adj_points(DT_edges(:,1),:) .* adj_points(DT_edges(:,2),:));

edgesangles = NaN(size(DT_edges));
adj = zeros(n_tri); % Adjacency matrix indicating which triangles are neighbours.
adj_convex = zeros(n_tri);
for edge=1:size(DT_edges,1)
    % The angles on either side of the edge.
    tri = adj_edge(edge,:);
    t = adj_angle(DT_edges(edge,:), tri );
    edgesangles(edge,:) = sum(t, 2);
    tri_idx = find(tri);
    adj(tri_idx,tri_idx) = 1;
    adj_convex(tri_idx,tri_idx) = prod(edgesangles(edge,:) <= 180);
end
convexedges = (edgesangles <= 180);
% Set diagonals to zero.
adj(logical(eye(n_tri))) = 0;
adj_convex(logical(eye(n_tri))) = 0;

但是，如果我想要所有组合或最大的凸包，我不确定如何进行。而且我不确定如何解释一个完整的圆中有几个三角形（即360度）的特殊情况。

4

2 回答 2

2

这个答案可以用一个简单的递归算法来解决：

  1. 从任何三角形开始
  2. 查找与该三角形共享边的所有三角形
  3. 对于这些三角形中的每一个，检查这个集合是否是凸的。
    • 如果不凸则停止
    • 添加连接到最后添加的三角形的三角形以设置仍需要测试的设置

所以这个算法是递归的，因为它贪婪地尝试向集合中添加更多的三角形，直到它不再是凸的。下面的代码给出了这个结果（我省略了所有琐碎的（1 个三角形）答案。

检查凸包中是否没有内部点有点幼稚：构造凸包并查看是否所有点都在上面。

在此处输入图像描述

function [convexTriangleSets,DT] = triangles()
    % inputs
    xy = [3.3735 0.7889; -0.1072 -3.4814; -3.9732 4.1955; -5 5; 5 5; 5 -5; -5 -5];
    DT = delaunayTriangulation(xy);

    function convexTriangleSets = testAddTriangle(iTriangle,attachedTriangleIndices,includedTriangleIndices)
        % add triangle
        includedTriangleIndices(end+1) = iTriangle;

        % naive test if allpoint of set of triangles are on convex hull
        nodes = unique(DT(includedTriangleIndices,:));
        coords = DT.Points(nodes,:);
        ch = convexHull(delaunayTriangulation(coords)); 
        allNodesOnConvexHull = length(nodes) == length(ch)-1;

        if ~allNodesOnConvexHull
            convexTriangleSets = {};
            return 
        end

        % find triangles connected to iTriangle
        currentTriangle = DT.ConnectivityList(iTriangle,:)';

        attachedCell = DT.edgeAttachments(currentTriangle([1 2 3]),currentTriangle([2 3 1]));
        attachedRow = unique([attachedTriangleIndices,attachedCell{:}]);
        attachedTriangleIndices = attachedRow(~ismember(attachedRow,includedTriangleIndices));

        % recursively try to expand connected triangles
        convexTriangleSets = {sort(includedTriangleIndices)};
        for ii = 1:length(attachedTriangleIndices)
            convexTriangleSets = [convexTriangleSets,...
                testAddTriangle(attachedTriangleIndices(ii),...
                                attachedTriangleIndices,...
                                includedTriangleIndices)]; %#ok<AGROW>
        end
    end

    includedTriangleIndices = [];
    attachedTriangleIndices = [];
    convexTriangleSets = {};
    for iTriangle = 1:DT.size
        convexTriangleSets = [convexTriangleSets,...
            testAddTriangle(iTriangle,attachedTriangleIndices,includedTriangleIndices)]; %#ok<AGROW>
    end

    % filter single triangles
    convexTriangleSets(cellfun(@length,convexTriangleSets) == 1) = [];

    % filter unique sets; convert to string because matlab cannot unique a cell array
    [~,c] = unique(cellfun(@(x) sprintf('%d,',x),convexTriangleSets,'UniformOutput',false));

    convexTriangleSets = convexTriangleSets(c);

    % plot result
    n = ceil(sqrt(length(convexTriangleSets)));
    for kk = 1:length(convexTriangleSets)
        subplot(n,n,kk)
        triplot(DT,'k')
        hold on
        patch('faces',DT(convexTriangleSets{kk},:), 'vertices', DT.Points, 'FaceColor','r');
    end
end
于 2016-02-12T14:48:32.147 回答
1

希望您的一些问题可以使用这行代码得到解答

多边形的凸性

假设您有一组三角形，那么您应该删除该多边形内的所有点，即只查看多边形边界处的点。您可以使用inpolygon函数检查点是否在多边形内

[in,on] = inpolygon(xq,yq,xv,yv);

在这种情况下，我假设所有角度都必须小于 180° 的条件是足够的。您也可以只创建一个凸包，并检查该集合是否相同。使用convhull

K = convhull(X,Y)

是更大的子集吗？

您解决的一个问题是，您如何检查一个（凸）多边形是否是另一个（凸）多边形的子集。让我们假设以下内容；您有一个目标多边形，并且您知道该多边形中的所有三角形。接下来，您还知道包含至少一个这些多边形的所有其他（凸）多边形。然后，您可以再次使用检查一个多边形是否是另一个多边形的子集inpolygon

x %// target polygon x
y
xt %// test polygon x
yt
in = inpolygon(x, v, xt, yt);
if sum(in)==length(x)
   %// target polygon is subset of testpolygon
end

查找所有条件

在这里，我将使用-say- poor mans 方法。只需暴力循环所有组合并检查凸度。DelaunayTriangulation 为您提供了所有三角形的列表。IE

DT(:,:)

给你所有的三角形和它的点。

delaunayTriangulation 类

您已经在使用 delaunayTriangulation 类，为什么不使用该类的方法呢？你基本上有你需要的一切

  1. 检查三角剖分是否是其他三角剖分的子集
  2. 获取三角剖分的边界点。
  3. 获取三角剖分的凸包。

你需要更多吗？

于 2016-02-12T07:34:37.833 回答