MATLAB: Efficient method to find all intersections of triangulation edges

edgeintersectionMATLABmeshtriangulation

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

Best Answer

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.
Related Question