# Looking for where line intersects on itself

4 views (last 30 days)

Show older comments

### Accepted Answer

Bruno Luong
on 15 Jul 2022

Edited: Bruno Luong
on 17 Jul 2022

If you interpolate your data by line segments it cross it selft at least 135 times

data = readtable('data.csv');

xy=table2array(data)';

X=polyselfx(xy);

fprintf('Number of intersections = %d\n', size(X,2))

figure

plot(xy(1,:),xy(2,:))

hold on

plot(X(1,:),X(2,:),'rx')

function X = polyselfx(P)

% Empty buffer

X = zeros(2,0);

filled = 0;

sizec = 0;

for n=2:size(P,2)-1

cn = seg2poly(P(:,n:n+1), P(:,1:n-1));

m = size(cn,2);

filled = filled+m;

% Buffer too small

if sizec < filled

sizec = max(filled, 2*sizec);

X(2,sizec) = 0;

end

% Store the result

X(:,filled+(-m+1:0)) = cn;

end

% remove the tails

X(:,filled+1:end) = [];

end % polyselfx

%%

function [cross_X, cross_t1, cross_t2, cross_i2] = seg2poly(s1, P)

a = s1(:,1);

M = P-a;

b = s1(:,2)-a;

nb2 = b(1)^2+b(2)^2;

% Check if the points are on the left/right side

x = [b(2) -b(1)]*M; % 1 x n

sx = sign(x);

% x -coordinates has opposite signs

crossx = sx(1:end-1).*sx(2:end) <= 0;

ix = find(crossx);

% cross point to the y-axis (along the segment)

x1 = x(ix);

x2 = x(ix+1);

d = b.'/nb2;

y1 = d*M(:,ix);

y2 = d*M(:,ix+1);

dx = x2-x1;

t1 = (y1.*x2-y2.*x1)./dx;

% Check if the cross point is inside the segment

ind = t1>0 & t1<=1;

if any(ind)

cross_t1 = t1(ind);

cross_t2= -x1(ind)./dx(ind);

cross_X = a + b*cross_t1;

cross_i2 = ix(ind);

else

cross_X = zeros(2,0);

cross_t1 = zeros(1,0);

cross_t2 = zeros(1,0);

cross_i2 = [];

end

end % seg2poly

##### 5 Comments

Bruno Luong
on 18 Jul 2022

Edited: Bruno Luong
on 18 Jul 2022

Here we go I modify the code and you get in idx the (first) indexes of the line segments and the fractional parametriec coordinates t of the intersections

data = readtable('data.csv');

xy=table2array(data)';

[X,idx,t]=polyselfx(xy);

fprintf('Number of intersections = %d\n', size(X,2))

figure

plot(xy(1,:),xy(2,:))

hold on

plot(X(1,:),X(2,:),'rx')

function [X, loc, t] = polyselfx(P)

% Crossing points polygons P

%

% INPUTS:

% P: two-row arrays, each column is a vertices

% OUTPUTS:

% X is two-row array, each column is an crossing point

% loc: two-row array, which edges the crossing point belong?

% first row corresponds to P1, second row to P2

% edge#1 is P(:,[1 2]), edge#2 is P(:,[2 3]), ... etc

% t: floating parametric of the crossing point

% Empty buffer

X = zeros(2,0);

loc = zeros(2,0);

t = zeros(2,0);

filled = 0;

sizec = 0;

for n=2:size(P,2)-1

[cn, t1, t2, i2] = seg2poly(P(:,n:n+1), P(:,1:n-1));

m = size(cn,2);

filled = filled+m;

% Buffer too small

if sizec < filled

sizec = max(filled, 2*sizec);

X(2,sizec) = 0;

loc(2,sizec) = 0;

t(2,sizec) = 0;

end

% Store the result

ifill = filled+(-m+1:0);

X(:,ifill) = cn;

loc(1,ifill) = n;

loc(2,ifill) = i2;

t(1,ifill) = t1;

t(2,ifill) = t2;

end

% remove the tails

X(:,filled+1:end) = [];

loc(:,filled+1:end) = [];

t(:,filled+1:end) = [];

end % polyselfx

%%

function [cross_X, cross_t1, cross_t2, cross_i2] = seg2poly(s1, P)

a = s1(:,1);

M = P-a;

b = s1(:,2)-a;

nb2 = b(1)^2+b(2)^2;

% Check if the points are on the left/right side

x = [b(2) -b(1)]*M; % 1 x n

sx = sign(x);

% x -coordinates has opposite signs

crossx = sx(1:end-1).*sx(2:end) <= 0;

ix = find(crossx);

% cross point to the y-axis (along the segment)

x1 = x(ix);

x2 = x(ix+1);

d = b.'/nb2;

y1 = d*M(:,ix);

y2 = d*M(:,ix+1);

dx = x2-x1;

t1 = (y1.*x2-y2.*x1)./dx;

% Check if the cross point is inside the segment

ind = t1>0 & t1<=1;

if any(ind)

cross_t1 = t1(ind);

cross_t2= -x1(ind)./dx(ind);

cross_X = bsxfun(@plus, a, b*cross_t1);

cross_i2 = ix(ind);

else

cross_X = zeros(2,0);

cross_t1 = zeros(1,0);

cross_t2 = zeros(1,0);

cross_i2 = [];

end

end % seg2poly

### More Answers (0)

### See Also

### Community Treasure Hunt

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

Start Hunting!