Combinations of triangles that form a convex hull

Question

I have a set of triangles. I am looking for a way to find all combinations of these triangles that constitute a convex hull when joined together. The convex hull should be empty, ie. no points inside the convex hull only on the edge. And only triangles that share a side can be joined together, ie. no gaps in the union.

Example: The following points gives 12 triangles (Delaunay Triangulation).

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));

I am looking for the largest convex hulls, so the convex hull should include as many triangles as possible. But if I have all the possible combinations I can easily filter out those with less triangles. In the example from above I should end up with these six convex hulls:

I am guessing that I should use that every triangle has at most three neighbouring triangles (one for each side). And that I should check whether the sum of the angle is less or equal 180 degrees at the points of intersect. This would insure that the union is convex -- see figure below. (The angle might also be exactly 360 degrees if several triangles form a full circle).

---

Triangle angles:

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

I reformat the information such that rows are points and columns are triangles:

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

I loop through all the edges and calculate the angles on both sides of the edge (edges angles). This way I am able to find pairs of triangles that form a convex set (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;

However I am unsure how to proceed if I want all combinations, or the largest convex hull. And I am unsure how to account for the special case where several triangles from a full circle (ie. 360 degrees).

Answer 1

Hopefully some of your concerns can be answered using this lines of code

Convexity of polygon

Assume you have a set of triangles, then you should remove all points inside that polygon, i.e. do only look at points at the border of the polygon. You can check whether or not a point is inside the polygon using the inpolygon function

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

The condition, that all angles need to be less that 180° is in this case i assume sufficient. You could also just create a convex hull, and check whether or not the set is identical. Use convhull

K = convhull(X,Y)

Is subset of larger one?

One question you addressed is, how you do check if one (convex) polygon is a subset of another (convex) polygon. Let's assume the following; you have a target polygon and you know all triangles in that polygon. Next you also know all other (convex) polygons containing at least one of those polygons. Then, you can check if one polygon is a subset of another by using again 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

The find all condition

Here I would use the - say - poor mans approach. Simply brute-force loop through all combinations and check for convexity. A list of all triangles is given to you by the DelaunayTriangulation. I.e.

DT(:,:)

gives you all triangles and its points.

The delaunayTriangulation class

You are already using the delaunayTriangulation class, why not use methods of that class? You have basically all you need

1. Check if triangulation is subset of other triangulation
2. Get boundary points of triangulation.
3. Get the convex hull of a triangulation.

Do you need more?

Answer 2

This answer can be solved with a simple recursive algorithm:

1. Start with any triangle
2. Find all triangles that share an edge with this triangle
3. For each of these triangles, check if this set is convex.
   • Stop if not convex
   • Add triangle connected to last added triangle to set that still needs to be tested

So this algorithm is recursive, as it greedily tries to add more triangles to the set until it is not convex any more. The code below gives this result (I omitted all trivial (1 triangle) answers.

The check if there are no interior points in the convex hull is a bit naive: construct convex hull and see if all points are on it.

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