Efficient method to find all intersections of triangulation edges?

Hi,

I have a triangulation object which is a list of vertices and connection list to create edges which form triangular mesh. I am looking for an efficient solution to find all missing intersections between all edges. In the plot below green dots are original vertices and red crosses are intersections I would like to find:

Example code is below:

clc
clearvars

% 3x3 1 internal triangle
[xx, yy] = meshgrid(1:3, 1:3);
V = [[xx(:), yy(:)]; 1.5 1.75; 2 2.5; 2.5 1.75];
V(4, :) = [];
F = [1 2 4; 2 3 5; 2 4 5; 4 6 7; 4 5 7; 5 7 8; 9 10 11];

% Create triangulation object
tr = triangulation(F, V);

% Plot
fig = figure(298);
clf(fig);
set(fig, 'Color', [1 1 1], 'WindowState', 'Normal')
ax = axes(fig, 'OuterPosition', [0 0 1 1], 'Clipping', 'off', 'XLim', [xx(1)-0.1 xx(end)+0.1], 'YLim', [yy(1)-0.1 yy(end)+0.1], 'Box', 'on', ...
    'XGrid', 'on', 'XMinorGrid', 'on', 'YGrid', 'on', 'YMinorGrid', 'on', 'DataAspectRatio', [1 1 1]);
p = patch(ax, 'faces', F, 'vertices' ,V, 'FaceColor', [0 0.5 1], 'Marker', '.', 'MarkerEdgeColor', 'g', 'MarkerSize', 18);

% Find and add all intersections
warning('This needs to be replaced with a function/algorithm')
missingIntersections = [1.75 1.75; 2.25 1.75; 5/3 2; 7/3 2]; % <-- this is a 'manual' solution

line(ax, missingIntersections(:,1), missingIntersections(:,2), 'LineStyle', 'none', 'Marker', 'x', 'Color', 'r', 'LineWidth', 2, 'MarkerSize', 16)

Are there any functions available that could find all intersections efficiently? I could loop through all edges and check if it intersects with any other edge but that would be quite slow and I imagine there are more efficient methods for that. The actual application of this is much bigger and thus looping through is way too slow.

Many Thanks,

Titas

Answer 1

[Matt J]

The closed form parametric solution for the intersection of two edges is easily derived using the Symbolic Toolbox:

syms xa xb xc xd ya yb yc yd s t
equ1=[xa;ya]+s*[xb-xa;yb-ya]; %edge from [xa,ya] to [xb,yb]
equ2=[xc;yc]+t*[xd-xc;yd-yc]; %edge from [xc,yc] to [xd,yd]
 
sol=solve(equ1==equ2,[s,t]);

s=vectorize(sol.s)

s = '(xa.*yc - xc.*ya - xa.*yd + xd.*ya + xc.*yd - xd.*yc)./(xa.*yc - xc.*ya - xa.*yd - xb.*yc + xc.*yb + xd.*ya + xb.*yd - xd.*yb)'

t=vectorize(sol.t)

t = '-(xa.*yb - xb.*ya - xa.*yc + xc.*ya + xb.*yc - xc.*yb)./(xa.*yc - xc.*ya - xa.*yd - xb.*yc + xc.*yb + xd.*ya + xb.*yd - xd.*yb)'

As you can see, these expressions are quite vectorizable. If you gather vectors xa,ya,xb,yb, etc... of vertex coordinates, you can compute s and t for all edge intersections in a single statement. Of course, a solution (s,t) only corresponds to a physical intersection if , so you will have to do some post-processing to discard non-physical solutions. But that step, too, is eminently vectorizable.