How to find the intersection points of a plot that contains multiple straight lines?

2 views (last 30 days)
Hi,
I have a simple simulation in which I spawn small lines with random position and direction. As time propagates, each line grows linearly from both its end points (denoted in the code as U and D) at a rate "lambda".
A plot figure is then presented upon each time-step, plotting the different lines simultaneously.
I would like to check, upon each time-step, if an intersection happens, and, if does happen, what's the coordinate of intersection and between which two lines did this intersection happen?
I tried InterX and polyxpoly; I don't think they can be applied here. This is my code:
clear all
close all
N=3; %number of spawns
initpos_x=rand(1,N)-0.5;
D_endpos_x=initpos_x;
U_endpos_x=initpos_x;
initpos_y=rand(1,N)-0.5;
D_endpos_y=initpos_y;
U_endpos_y=initpos_y;
direction=pi*rand(1,N)-pi/2;
figure
scatter(initpos_x,initpos_y)
xlim([-0.5 0.5])
ylim([-0.5 0.5])
drawnow
hold on
n=1000000; %number of time steps
lambda=0.001; %rate of growth.
for t=2:n
D_endpos_x=D_endpos_x-lambda*cos(direction);
U_endpos_x=U_endpos_x+lambda*cos(direction);
D_endpos_y=D_endpos_y-lambda*sin(direction);
U_endpos_y=U_endpos_y+lambda*sin(direction);
h=plot([D_endpos_x ; U_endpos_x],[D_endpos_y ; U_endpos_y]);
set(h, {'color'}, num2cell(0.8*jet(N),2));
drawnow;
end
Thank you!

Accepted Answer

Jan
Jan on 9 May 2022
Edited: Jan on 9 May 2022
N = 5; %number of spawns
X1 = rand(1,N) - 0.5;
X2 = X1;
Y1 = rand(1,N) - 0.5;
Y2 = Y1;
D = pi * rand(1,N) - pi/2;
sD = sin(D);
cD = cos(D);
figure;
axes('XLim', [-1, 1], 'YLim', [-1, 1], 'NextPlot', 'add');
hLine = line([X1; X2], [Y1; Y2]);
hDot = [];
set(hLine, {'color'}, num2cell(0.8*jet(N), 2));
n = 100; % number of time steps
g = 0.01; % rate of growth.
for t = 1:n
X1 = X1 - g * cD;
X2 = X2 + g * cD;
Y1 = Y1 - g * sD;
Y2 = Y2 + g * sD;
for k = 1:numel(hLine)
set(hLine(k), 'XData', [X1(k); X2(k)], 'YData', [Y1(k); Y2(k)]);
end
% Check for intersection now:
% [xi, yi] = polyxpoly(X1, Y1, X2, Y2); % If you have the mapping toolbox
% Or try your own line intersection code:
[xi, yi] = mySegmentCrossing(X1, Y1, X2, Y2);
delete(hDot);
hDot = plot(xi, yi, 'ro');
pause(0.02);
end
function [xi, yi] = mySegmentCrossing(X1, Y1, X2, Y2)
% (C) 2022, Jan, Heidelberg, CC BY-SA 3.0
n = numel(X1);
Result = cell(2, n);
for k = 1:n - 1
x1 = X1(k);
x2 = X2(k);
x3 = X1(k+1:n); % check crossing with following segments only
x4 = X2(k+1:n);
y1 = Y1(k);
y2 = Y2(k);
y3 = Y1(k+1:n);
y4 = Y2(k+1:n);
x13 = x1 - x3;
y13 = y1 - y3;
x34 = x3 - x4;
y34 = y3 - y4;
u = (x13 .* (y1-y2) - y13 .* (x1-x2)) ./ ...
((x1-x2) .* y34 - (y1-y2) .* x34);
t = (x13 .* y34 - y13 .* x34) ./ ...
((x1-x2) .* y34 - (y1-y2) .* x34);
c = (u >= 0 & u <= 1.0) & (t >= 0 & t <= 1.0); % TRUE for crossings
% Take mean of both cross points to reduce noise:
Result{1, k} = ((x3(c) - u(c) .* x34(c)) + (x1 + t(c) .* (x2-x1))) / 2;
Result{2, k} = ((y3(c) - u(c) .* y34(c)) + (y1 + t(c) .* (y2-y1))) / 2;
end
xi = cat(2, Result{1, :});
yi = cat(2, Result{2, :});
end
A further improvement is to avoid redrawing all found crossings:
delete(hDot); % Simple and expensive
hDot = plot(xi, yi, 'ro');
Use scatter once and update the XData and YData only as for hLine.

More Answers (1)

Cameron B
Cameron B on 6 Jan 2020
This currently only works for the three lines. There's a way to do it for multiple "spawns" but I didn't have time to look at that.
clear all
close all
N=3; %number of spawns
initpos_x=rand(1,N)-0.5;
D_endpos_x=initpos_x;
U_endpos_x=initpos_x;
initpos_y=rand(1,N)-0.5;
D_endpos_y=initpos_y;
U_endpos_y=initpos_y;
direction=pi*rand(1,N)-pi/2;
figure
plot(initpos_x,initpos_y,'o');
%xlim([-0.5 0.5])
%ylim([-0.5 0.5])
drawnow
hold on
n=1000000; %number of time steps
lambda=0.008; %rate of growth.
D_pt2_end_x = D_endpos_x-lambda*cos(direction);
U_pt2_end_x = U_endpos_x+lambda*cos(direction);
D_pt2_end_y = D_endpos_y-lambda*sin(direction);
U_pt2_end_y = U_endpos_y+lambda*sin(direction);
poly1 = polyfit([D_pt2_end_x(1), U_pt2_end_x(1)],[D_pt2_end_y(1), U_pt2_end_y(1)],1);
poly2 = polyfit([D_pt2_end_x(2), U_pt2_end_x(2)],[D_pt2_end_y(2), U_pt2_end_y(2)],1);
poly3 = polyfit([D_pt2_end_x(3), U_pt2_end_x(3)],[D_pt2_end_y(3), U_pt2_end_y(3)],1);
x12 = (poly2(2)-poly1(2))/(poly1(1)-poly2(1));
x13 = (poly3(2)-poly1(2))/(poly1(1)-poly3(1));
x23 = (poly3(2)-poly2(2))/(poly2(1)-poly3(1));
for t=2:n
D_endpos_x=D_endpos_x-lambda*cos(direction);
U_endpos_x=U_endpos_x+lambda*cos(direction);
D_endpos_y=D_endpos_y-lambda*sin(direction);
U_endpos_y=U_endpos_y+lambda*sin(direction);
h=plot([D_endpos_x ; U_endpos_x],[D_endpos_y ; U_endpos_y]);
legend('start','1','2','3')
if D_endpos_x(1) <= x12 && U_endpos_x(1) >= x12 && D_endpos_x(2) <= x12 && U_endpos_x(2) >= x12
plot(x12,poly1(1)*x12+poly1(2),'ko');
end
if D_endpos_x(1) <= x13 && U_endpos_x(1) >= x13 && D_endpos_x(3) <= x13 && U_endpos_x(3) >= x13
plot(x13,poly1(1)*x13+poly1(2),'ro');
end
if D_endpos_x(2) <= x23 && U_endpos_x(2) >= x23 && D_endpos_x(3) <= x23 && U_endpos_x(3) >= x23
plot(x23,poly2(1)*x23+poly2(2),'bo');
end
set(h, {'color'}, num2cell(0.8*jet(N),2));
drawnow;
end

Categories

Find more on Function Creation in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!