How to find intersection between a straight line and a curve?

I want to find the coordinates for the two intersection points between the horisontal line and the graph. I think the problem is that there is no matching value between the line and the elements in the vector, so when I try to use intersect() I get an empty set. Any idea how this can be done? Thx .

On the y-axis I have a vector with intensity levels, and on the x-axis position. The data is from a file.

az = load('profile_AZ_7cm_Full.mat');
el = load('profile_EL_7cm_Full.mat');

x = az.Position(1,1:end); 
intensity_az = 10*log10(sum(az.RF.^2));
normx_intensity = intensity_az-max(intensity_az);
y = el.Position(2,1:end);
intensity_el = 10*log10(sum(el.RF.^2));
normy_intensity = intensity_el-max(intensity_el);

l = -6;
isx = interp1(normx_intensity,l);

figure;
subplot(2,1,1)
plot(x*1e3,normx_intensity)
hold on
yline(l)
title('Azimuth scan of beam profile')
xlabel('x-axis (mm)')
ylabel('Intensity (dB)')

Answer 1

[Star Strider]

Another approach —

x = linspace(-15, 15, 150);
y = sinc(0.5*x);
l = 0.6;

[~,idx] = max(y);
xq(1) = interp1(y(1:idx),x(1:idx), l);
xq(2) = interp1(y(idx+1:end), x(idx+1:end), l);

figure
plot(x, y)
hold on
plot(x, 0.6*ones(size(x)), ':k')
plot(xq, [1 1]*l, 'xr', 'MarkerSize',10)
hold off
text(xq, [1 1]*l, compose('\\leftarrow x = %+.3f',xq), 'Horiz','left', 'Vert','middle', 'Rotation',-90)

.

[Davide Masiello]

Nice! And definitely faster than what I suggested.

But I think this works well for the specific case given in the question, right?

I was thinking of a more general approach when I gave my answer.

[Star Strider]

Thank you!

My solution is intended to work specifically for this problem, so yes.

To get multiple intersections in different problems, I always begin with:

zxi = find(diff(sign(y - l)));

This creates ‘zero-crossings’ where ‘y’ is approximately equal to ‘l’ (in the context of the current problem), and returns their indices. This then guides the rest of the code:

for k = 1:numel(zxi)

idxrng = max([1, zxi(k)-2]) : min([numel(y), zxi(k)+2]); % Index Range For Interpolation

xi(k) = interp1(y(idxrng), x(idxrng), l); % High-Precision Estimate Of Intersection Coordinate

end

Testing it:

x = linspace(-15, 15);
y = sinc(0.5*x);

l = 0.05;

zxi = find(diff(sign(y - l)));
for k = 1:numel(zxi)
   idxrng = max([1, zxi(k)-2]) : min([numel(y), zxi(k)+2]);                % Index Range For Interpolation
   xi(k) = interp1(y(idxrng), x(idxrng), l);                               % High-Precision Estimate Of Intersection Coordinate
end   

figure
plot(x, y)
hold on
plot(x, ones(size(x))*l, ':k')
plot(xi, ones(size(xi))*l, 'xr', 'MarkerSize',7.5)
hold off

This gets around the problem with all interpolation methods (even hand-coded original versions) in that they need to be given a specific region over which to interpolate. After that, the interpolation proceeds straight away! (The loop is obviously necessary.) The index range ‘idxrng’ calculation traps situations where a zero-crossing is at the beginning or end of the sequence so it will not reference points that do not exist. The accuracy of the interpolation to calculate the intersections is in part a function of the resolution of the data, so initially interpolating to a finer resolution (more data points) is occasionally required.

.

[Davide Masiello]

That's great, thanks very much for the detailed explanation!

[Star Strider]

My pleasure!

Feel free to use that approach in your answers as well. (We generally learn from each other here, so that is not considered ‘cheating’, and is instead ‘learning’.) For example the ‘zxi’ code is an improvement on an approach I previously used that someone else developed, so I just adopted it to use with my code.

[Priyam Mishra]

Sir, I have same kind of problem in which in which blue color line is the cone tip data which is plotted along the depth. A linear fitted line with the same data is also shown. I want to find the values of depth (y-axis) where they are intersecting to each other. I want to find out all the points of intersection. Please help me Sir.

Thanks for your kind response.

Priyam Mishra

[Star Strider]

Change the variables —

y = linspace(-15, 15);
x = sinc(0.5*y);

l = 0.05;

zxi = find(diff(sign(x - l)));
for k = 1:numel(zxi)
   idxrng = max([1, zxi(k)-2]) : min([numel(y), zxi(k)+2]);                % Index Range For Interpolation
   yi(k) = interp1(x(idxrng), y(idxrng), l);                               % High-Precision Estimate Of Intersection Coordinate
end   

figure
plot(x, y)
hold on
plot(ones(size(x))*l, y, ':k')
plot(ones(size(yi))*l,  yi, 'xr', 'MarkerSize',7.5)
hold off

.

Answer 2

[Davide Masiello]

I am a bit ovewhelmed by your code, so I will give you a qualitative example instead of modifying it.

x = [x1,x2,x3,___,xn];      % Your array of x data
y = [y1,y2,y3,___,yn];      % Your array of y data
c = c0;                     % The y value of your straight horizontal line

f = @(z) interp1(x,y-c,z);  % Creates anonymous function representing the difference between curve and straight line

xi = fzero(f,x0);           % Finds zero of the function (i.e., the intersection point)

The value of xi depends on the the initial condition x0, so you might want to loop the fzero over several x0 values to find multiple intersection points. It could be like this:

x0 = linspace(x1,xn,5);
for i = 1:5
    xi(i) = fzero(f,x0(i));           % Finds zero of the function (i.e., the intersection point)
end

Then you can do

xi = unique(xi);

to remove repeated values.

Hope it makes sense.