How can i find the intersection of two lines between points? One of the lines isnt well behaved

I need to find the intersection values of x for these two lines between output points

[Adam Danz]

Check out this function on the exchange

https://www.mathworks.com/matlabcentral/fileexchange/11837-fast-and-robust-curve-intersections

[Guillaume]

What's a badly behaved line? Sings loudly outside your house at 2AM? Is slightly curvy at the ends?

[John D'Errico]

Guillaume - I wish I could upvote comments. :)

Answer 1

[Star Strider]

Try this:

x = linspace(0, 1, 40);                                     % Create Data
y1 = 1.0592*ones(size(x));                                  % Create Data
y1(18) = 1.0251;                                            % Create Data
y2 = 0.093*x + 1;

xi = linspace(min(x), max(x), numel(x)*100);                % Interpolation Vector
y1i = interp1(x, y1, xi);                                   % Interpolate
y2i = interp1(x, y2, xi);                                   % Interpolate
yd = y1i - y2i;                                             % Difference
zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0);        % Returns Approximate Zero-Crossing Indices Of Argument Vector
ydi = zci(yd);                                              % Zero-Crossing Indices

for k = 1:2
    ix = (ydi(k)-1):(ydi(k)+1);
    xv = xi(ix);
    yv = yd(ix);
    B = [xv(:), ones(size(xv(:)))] \  yv(:);
    xval(k) = -B(2)/B(1);                                   % Intersection X-Coordinates <—
end

y2r = [x(:) ones(size(y2(:)))] \ y2(:);                     % Regression For ‘y2’
yval = [xval(:) ones(2,1)] * y2r;                           % Intersection Y-Coordinates <—

figure
plot(xi, y1i)
hold on
plot(xi, y2i)
plot(xval, yval, 'pg')
hold off
grid

Here ‘xval’ are the x-values of the intersection. This computes the associated ‘yval’ vector and plots the points as well.

[Alex Earle-Sodhi]

Thank you, this is really helpful!

However this graph comes from a procedure which i would like to iterate and can have more than 2 intersections, I see that changing for k = 1:2 to for k = 1:numel(ydi) i can incorporate all of the intersections.

However it also includes point 4000 and then xv cannot find a point 4000+1 is there any way I can avoid this and allow any number of intersections to occur?

[Star Strider]

As always, my pleasure!

One of the problems with simple anonymous function implementation of ‘zci’ is that if the end of a vector has a different sign than the beginning of the vector, it will pick that up as a zero-crossing, because of the ‘wrap-around’ effect of circshift operation. Logic to detect that and correct for it, in this context is:

if yd(1)*yd(end) < 1                                        % Check For ‘Wrap-Around’ Effect
    ydi = ydi(1:end-1)
end

and the slightly revised version of my code then becomes:

x = linspace(0, 1, 40);                                     % Create Data
y1 = 1.0592*ones(size(x));                                  % Create Data
y1(18) = 1.0251;                                            % Create Data
y2 = 0.093*x + 1;

xi = linspace(min(x), max(x), numel(x)*100);                % Interpolation Vector
y1i = interp1(x, y1, xi);                                   % Interpolate
y2i = interp1(x, y2, xi);                                   % Interpolate
yd = y1i - y2i;                                             % Difference
zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0);        % Returns Approximate Zero-Crossing Indices Of Argument Vector
ydi = zci(yd);                                              % Zero-Crossing Indices
if yd(1)*yd(end) < 1                                        % Check For ‘Wrap-Around’ Effect
    ydi = ydi(1:end-1)
end

for k = 1:numel(ydi)
    ix = (ydi(k)-1):(ydi(k)+1);
    xv = xi(ix);
    yv = yd(ix);
    B = [xv(:), ones(size(xv(:)))] \  yv(:);
    xval(k) = -B(2)/B(1);                                   % Intersection X-Coordinates <—
end

y2r = [x(:) ones(size(y2(:)))] \ y2(:);                     % Regression For ‘y2’
yval = [xval(:) ones(size(xval(:)))] * y2r;                 % Intersection Y-Coordinates <—

figure
plot(xi, y1i)
hold on
plot(xi, y2i)
plot(xval, yval, 'pg')
hold off
grid

This detects, interpolates, and plots all the valid intersections.

[Alex Earle-Sodhi]

Simply wonderful!

I thought it might have been something to do with that function only via process of elimination aha

I shall apply this to my code,

Thanks alot!

[Star Strider]

Thank you!

As always, my pleasure!