5 views (last 30 days)

Show older comments

Walter Roberson
on 8 Mar 2020

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

David Goodmanson
on 10 Mar 2020

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.

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

Start Hunting!