Can I determine the points (x,y) with a single line looping back and cross its self?

I have a plot(x,y) where the data loop back and crosses at a certain point (x,y). Can Matlab tell me where it intersect?

[Walter Roberson]

Is the graph in the form of an image, or is it a list of coordinates with linear connection between coordinates assumed?

[A C]

It is a line on a plot(x,y) that loops around and crosses itself.

[Walter Roberson]

If Roger Stafford were still active he could probably come up with a nice approach...

At one point, Roger posted code to determine whether any two line segments intersected. For a smallish number of points you could probably afford to compare every segment to every other segment. However that would clearly be O(n^2), which would not work so well with larger lines.

For larger lines you would want to preprocess to see if the segments were in the same general area before bothering to compare them.

You could do something like divide the x range up into a number of bins. For each segment, mark all of the bins from its minimum x to maximum x as occupied by that segment number. This would generate a set of segment numbers for each bin, but is something that can be done in less than n^2 time (the time for any one segment is proportionate to the range of x it spans, expressed as number of bins. Once you have a list of bin numbers for each segment, accumarray can be used to create the sets for each bin.)

Once you have a set of segment numbers for each bin, you can proceed to compare the segments within each bin to each other.

The worst case scenario for this is still n^2 but the best case is more like linear.

There just might be a way with lower cost using a variation of https://en.m.wikipedia.org/wiki/Bresenham%27s_line_algorithm

Answer 1

[David Goodmanson]

Hi A C

Good one.

I'm assuming that the x and y arrays are such that x(1:end),y(1:end) traces out a 'continuous' (so to speak) path. The code below outputs a small matrix of indices, such as

ind =
    23    24   114   115
   133   134   224   225

with one row for each crossing occurrence. The first row indicates that the line between (x,y) points 23 and 24 intersects the line between (x,y) points 114 and 115, and so forth. With that information there are many ways to find the exact crossing point. One of them is (after constructing points a,b,c,d, for the mth instance of line crossing),

abcd = [x(ind(m,:)); y(ind(m,:))];
a = abcd(:,1);  b = abcd(:,2);  c= abcd(:,3);  d = abcd(:,4);  % four column vectors of x,y coordinates

M = [b-a,d-c];
crossingpoint = (1/2)*( a+c + M*[1 0;0 -1]*inv(M)*(c-a) )

As Walter mentioned, without special effort the search will be O(n^2) and it's hard to get more O(n^2) than this code is. For 10000 points it takes only about 1.8 sec. on my PC, but uses about 5 GB of memory. And as bad as O(n^2) is going up, it's really great going back down. So it seems to be good for small to medium size sets of points.

The plot shows the crossing segment endpoints in color but does not show the exact crossing point calculation.

% make up some data
t = linspace(4,21,1000);
x = 1.4*cos(t) + .28*t;     % cycloid example
y = sin(t);

N = length(t)
% segments down, points across
xsegs = repmat(diff(x)',1,N);
ysegs = repmat(diff(y)',1,N);
xd = x-x(1:end-1)';                 % point displacements from the first point in each segment
yd = y-y(1:end-1)';
s = sign(xsegs.*yd - ysegs.*xd);    % denotes which side of the segment each point is on
sd = abs(diff(s,[],2));             % find line crossings, sign difference = +-2
[n1 n2] = find(triu(sd)==2 & triu(sd')==2);
ind = [n1 n1+1 n2  n2+1];

%plot
figure(1)
plot(x,y); grid on
hold on
n1n2 = [n1;n2]                       % mashup of all the indices
plot(x(n1n2(:)),y(n1n2(:)),'o')
plot(x(n1n2(:)+1),y(n1n2(:)+1),'o')
hold off

The code uses the principle that

for an extended line through points a and b, if points c and d are on opposite sides, and

for an extended line through points c and d, if points a and b are on opposite sides, then segment ab crosses segment cd.