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

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!

[Rena Berman]

(Answers Dev) Restored edit

Answer 1

[Jan]

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.

[J. D.]

Wonderful! Thanks for your help!

[Jan]

@J. D.: You are welcome. I've adjusted your code, so it was useful to find it here. I think this clarifies the former discussion.

Answer 2

[Cameron B]

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