How to find intersections in contour plot?

I've been trying to find the points where this contour intersects the line y = -3, but it's not working. I tried just counting the intersection points and got 0, when it should clearly be 2. How do I fix this?

This is my entire code btw, in case you want to try running it. The part where I'm trying to find the intersections starts at line 51, where I've divided the code sections.

clear all
close all
clc
%Constant
rho = 4420;          %kg/m^3
Cp = 550;            %J/kg?K
T0 = 303.15;         %K
A = 0.5;             %[Absorbtivity] 
k = 7.2;             %W/m/K
alpha = 2.96*10^-6;  %m^2/s
D = alpha;
P = 500;             %W
u = 0.5;               %m/s
Tm = 1933;        %K
d_laser = 0.0001;    %m
r_laser = d_laser/2; %m
a = r_laser;
p = D/(u*a);
%Define
n = 100;
% x = linspace(-0.00015,0.00025,n);
% z = linspace(-0.000005,0,n);
x = linspace(-0.00025,0.0012,n);
z = linspace(-0.0001,0,n);
%Normalized
x_nor = x/a;
z_nor = z/(D*a/u).^0.5;
[x_mesh,z_mesh] = meshgrid(x_nor,z_nor);
y_mesh = 0;
fun = @(t) exp((-z_mesh.^2./(4*t))-((y_mesh.^2+(x_mesh-t).^2)./(4*p.*t+1)))./((4.*p.*t+1).*sqrt(t));
g = integral(fun,0,Inf,'ArrayValued',true);
% Ts = (A*P)/(pi*rho*Cp*sqrt(D*u*(a^3)))
Ts = Tm/0.16713
g2 = Tm/Ts

figure
hold on
[C,h] = contourf(x_mesh,z_mesh,g, 'ShowText', 'On')
[C1,h1] = contour(x_mesh,z_mesh,g, 'ShowText', 'On',LevelList=0.16713,Color='red')
[C2,h2] = contour(x_mesh,z_mesh,g, 'ShowText', 'On',LevelList = g2,Color='white')
hold off
clabel(C,h,'Color','white')
clabel(C1,h1,'Color','white')
clabel(C2,h2,'Color','white')
% xticks(-5:5:25)
% yticks(-10:2:0)
% axis([-5 25 -10 0])
title('x\primez\prime plane')
xlabel('x\prime')
ylabel('z\prime')
colorbar

%Extract data from contour plot
contourdata = contourc(x_nor,z_nor,g,[g2 g2])

% Initialize arrays to store the boundary coordinates
        x_prime = [contourdata(1,2:size(contourdata,2))]
        z_prime = [contourdata(2,2:size(contourdata,2))]
        outlineplot = plot(x_prime,z_prime)
        x_outline = outlineplot.XData
        z_outline = outlineplot.YData
        outline = [x,z]
        z_line = yline(-3)

% Find the indices where the curve crosses the line y = 3
crossing_indices = find(outlineplot == -3)
% intersections = intersect(outline, x_line)
% Count the number of intersections
num_intersections = numel(crossing_indices)

% Display the result
disp(['The curve intersects the line z = -3 at ', num2str(num_intersections), ' point(s).']);

[Dyuman Joshi]

You are assuming that the plot data will have a point with y-coordinate equals to -3, which is not the case.

Also, you are comparing a graphics object, a plot (i.e. outlineplot variable), with a numerical value via this line of code, which does not make sense.

crossing_indices = find(outlineplot == -3)

And as a plot is not equal to a number, the output of the comparison is 0 thus find() returns an empty array.

If you want to compare, use the YData from the plot to compare.

There are some FEX files you can look into - Fast and Robust Curve Intersections and Curve Intersections

If you have the Mapping toolbox, polyxpoly would be helpful.

[Sarinya]

Thank you so much! I'm still a beginner so I guess I'm not well-versed on the different types of data and how you're supposed to work with them yet. This was very helpful.

Answer 1

[Chunru]

rho = 4420;          %kg/m^3
Cp = 550;            %J/kg?K
T0 = 303.15;         %K
A = 0.5;             %[Absorbtivity] 
k = 7.2;             %W/m/K
alpha = 2.96*10^-6;  %m^2/s
D = alpha;
P = 500;             %W
u = 0.5;               %m/s
Tm = 1933;        %K
d_laser = 0.0001;    %m
r_laser = d_laser/2; %m
a = r_laser;
p = D/(u*a);
%Define
n = 100;
% x = linspace(-0.00015,0.00025,n);
% z = linspace(-0.000005,0,n);
x = linspace(-0.00025,0.0012,n);
z = linspace(-0.0001,0,n);
%Normalized
x_nor = x/a;
z_nor = z/(D*a/u).^0.5;
[x_mesh,z_mesh] = meshgrid(x_nor,z_nor);
y_mesh = 0;
fun = @(t) exp((-z_mesh.^2./(4*t))-((y_mesh.^2+(x_mesh-t).^2)./(4*p.*t+1)))./((4.*p.*t+1).*sqrt(t));
g = integral(fun,0,Inf,'ArrayValued',true);
% Ts = (A*P)/(pi*rho*Cp*sqrt(D*u*(a^3)))
Ts = Tm/0.16713

Ts = 1.1566e+04

g2 = Tm/Ts

g2 = 0.1671

figure
hold on
[C,h] = contourf(x_mesh,z_mesh,g, 'ShowText', 'On')

C = 2×525

0.2000 -1.3514 -1.3179 -1.2815 -1.2421 -1.1992 -1.1919 -1.1694 -1.1407 -1.1099 -1.0767 -1.0410 -1.0024 -0.9609 -0.9160 -0.8990 -0.8787 -0.8452 -0.8095 -0.7713 -0.7305 -0.6867 -0.6398 -0.6061 -0.5941 -0.5556 -0.5148 -0.4715 -0.4255 -0.3767 177.0000 0 -0.0587 -0.1174 -0.1761 -0.2348 -0.2446 -0.2936 -0.3523 -0.4110 -0.4697 -0.5284 -0.5871 -0.6458 -0.7045 -0.7262 -0.7632 -0.8220 -0.8807 -0.9394 -0.9981 -1.0568 -1.1155 -1.1559 -1.1742 -1.2329 -1.2916 -1.3504 -1.4091 -1.4678

h = 
  Contour with properties:

    EdgeColor: [0 0 0]
    LineStyle: '-'
    LineWidth: 0.5000
    FaceColor: 'flat'
    LevelList: [1.3251e-13 0.2000 0.4000 0.6000 0.8000 1 1.2000 1.4000 1.6000 1.8000]
        XData: [100×100 double]
        YData: [100×100 double]
        ZData: [100×100 double]

  Show all properties

[C1,h1] = contour(x_mesh,z_mesh,g, 'ShowText', 'On',LevelList=0.16713,Color='red')

C1 = 2×203

0.1671 -1.4112 -1.3830 -1.3524 -1.3192 -1.2831 -1.2440 -1.2014 -1.1919 -1.1712 -1.1431 -1.1128 -1.0803 -1.0452 -1.0073 -0.9664 -0.9223 -0.8990 -0.8835 -0.8510 -0.8163 -0.7791 -0.7394 -0.6968 -0.6511 -0.6061 -0.6032 -0.5661 -0.5268 -0.4851 202.0000 0 -0.0587 -0.1174 -0.1761 -0.2348 -0.2936 -0.3523 -0.3652 -0.4110 -0.4697 -0.5284 -0.5871 -0.6458 -0.7045 -0.7632 -0.8220 -0.8518 -0.8807 -0.9394 -0.9981 -1.0568 -1.1155 -1.1742 -1.2329 -1.2871 -1.2916 -1.3504 -1.4091 -1.4678

h1 = 
  Contour with properties:

    EdgeColor: [1 0 0]
    LineStyle: '-'
    LineWidth: 0.5000
    FaceColor: 'none'
    LevelList: 0.1671
        XData: [100×100 double]
        YData: [100×100 double]
        ZData: [100×100 double]

  Show all properties

[C2,h2] = contour(x_mesh,z_mesh,g, 'ShowText', 'On',LevelList = g2,Color='white')

C2 = 2×203

0.1671 -1.4112 -1.3830 -1.3524 -1.3192 -1.2831 -1.2440 -1.2014 -1.1919 -1.1712 -1.1431 -1.1128 -1.0803 -1.0452 -1.0073 -0.9664 -0.9223 -0.8990 -0.8835 -0.8510 -0.8163 -0.7791 -0.7394 -0.6968 -0.6511 -0.6061 -0.6032 -0.5661 -0.5268 -0.4851 202.0000 0 -0.0587 -0.1174 -0.1761 -0.2348 -0.2936 -0.3523 -0.3652 -0.4110 -0.4697 -0.5284 -0.5871 -0.6458 -0.7045 -0.7632 -0.8220 -0.8518 -0.8807 -0.9394 -0.9981 -1.0568 -1.1155 -1.1742 -1.2329 -1.2871 -1.2916 -1.3504 -1.4091 -1.4678

h2 = 
  Contour with properties:

    EdgeColor: [1 1 1]
    LineStyle: '-'
    LineWidth: 0.5000
    FaceColor: 'none'
    LevelList: 0.1671
        XData: [100×100 double]
        YData: [100×100 double]
        ZData: [100×100 double]

  Show all properties

hold off
clabel(C,h,'Color','white')
clabel(C1,h1,'Color','white')
clabel(C2,h2,'Color','white')
% xticks(-5:5:25)
% yticks(-10:2:0)
% axis([-5 25 -10 0])
title('x\primez\prime plane')
xlabel('x\prime')
ylabel('z\prime')
colorbar

%Extract data from contour plot
contourdata = contourc(x_nor,z_nor,g,[g2 g2])

contourdata = 2×203

0.1671 -1.4112 -1.3830 -1.3524 -1.3192 -1.2831 -1.2440 -1.2014 -1.1919 -1.1712 -1.1431 -1.1128 -1.0803 -1.0452 -1.0073 -0.9664 -0.9223 -0.8990 -0.8835 -0.8510 -0.8163 -0.7791 -0.7394 -0.6968 -0.6511 -0.6061 -0.6032 -0.5661 -0.5268 -0.4851 202.0000 0 -0.0587 -0.1174 -0.1761 -0.2348 -0.2936 -0.3523 -0.3652 -0.4110 -0.4697 -0.5284 -0.5871 -0.6458 -0.7045 -0.7632 -0.8220 -0.8518 -0.8807 -0.9394 -0.9981 -1.0568 -1.1155 -1.1742 -1.2329 -1.2871 -1.2916 -1.3504 -1.4091 -1.4678

% Initialize arrays to store the boundary coordinates
        x_prime = [contourdata(1,2:size(contourdata,2))]

x_prime = 1×202

-1.4112 -1.3830 -1.3524 -1.3192 -1.2831 -1.2440 -1.2014 -1.1919 -1.1712 -1.1431 -1.1128 -1.0803 -1.0452 -1.0073 -0.9664 -0.9223 -0.8990 -0.8835 -0.8510 -0.8163 -0.7791 -0.7394 -0.6968 -0.6511 -0.6061 -0.6032 -0.5661 -0.5268 -0.4851 -0.4408

        z_prime = [contourdata(2,2:size(contourdata,2))]

z_prime = 1×202

0 -0.0587 -0.1174 -0.1761 -0.2348 -0.2936 -0.3523 -0.3652 -0.4110 -0.4697 -0.5284 -0.5871 -0.6458 -0.7045 -0.7632 -0.8220 -0.8518 -0.8807 -0.9394 -0.9981 -1.0568 -1.1155 -1.1742 -1.2329 -1.2871 -1.2916 -1.3504 -1.4091 -1.4678 -1.5265

        outlineplot = plot(x_prime,z_prime)

outlineplot = 
  Line with properties:

              Color: [0 0.4470 0.7410]
          LineStyle: '-'
          LineWidth: 0.5000
             Marker: 'none'
         MarkerSize: 6
    MarkerFaceColor: 'none'
              XData: [-1.4112 -1.3830 -1.3524 -1.3192 -1.2831 -1.2440 -1.2014 -1.1919 -1.1712 -1.1431 -1.1128 -1.0803 -1.0452 -1.0073 -0.9664 -0.9223 -0.8990 -0.8835 -0.8510 -0.8163 -0.7791 … ]
              YData: [0 -0.0587 -0.1174 -0.1761 -0.2348 -0.2936 -0.3523 -0.3652 -0.4110 -0.4697 -0.5284 -0.5871 -0.6458 -0.7045 -0.7632 -0.8220 -0.8518 -0.8807 -0.9394 -0.9981 -1.0568 -1.1155 … ]

  Show all properties

        x_outline = outlineplot.XData

x_outline = 1×202

-1.4112 -1.3830 -1.3524 -1.3192 -1.2831 -1.2440 -1.2014 -1.1919 -1.1712 -1.1431 -1.1128 -1.0803 -1.0452 -1.0073 -0.9664 -0.9223 -0.8990 -0.8835 -0.8510 -0.8163 -0.7791 -0.7394 -0.6968 -0.6511 -0.6061 -0.6032 -0.5661 -0.5268 -0.4851 -0.4408

        z_outline = outlineplot.YData

z_outline = 1×202

0 -0.0587 -0.1174 -0.1761 -0.2348 -0.2936 -0.3523 -0.3652 -0.4110 -0.4697 -0.5284 -0.5871 -0.6458 -0.7045 -0.7632 -0.8220 -0.8518 -0.8807 -0.9394 -0.9981 -1.0568 -1.1155 -1.1742 -1.2329 -1.2871 -1.2916 -1.3504 -1.4091 -1.4678 -1.5265

        outline = [x,z]

outline = 1×200

1.0e+00 * -0.0003 -0.0002 -0.0002 -0.0002 -0.0002 -0.0002 -0.0002 -0.0001 -0.0001 -0.0001 -0.0001 -0.0001 -0.0001 -0.0001 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0000 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0002 0.0002

        z_line = yline(-3)

z_line = 
  ConstantLine with properties:

    InterceptAxis: 'y'
            Value: -3
            Color: [0.1500 0.1500 0.1500]
        LineStyle: '-'
        LineWidth: 0.5000
            Label: ''
      DisplayName: ''

  Show all properties

% Find the indices where the curve crosses the line y = 3
crossing_indices = find(outlineplot == -3)

crossing_indices =

     []

% intersections = intersect(outline, x_line)
% Count the number of intersections
num_intersections = numel(crossing_indices)

num_intersections = 0

% Display the result
disp(['The curve intersects the line z = -3 at ', num2str(num_intersections), ' point(s).']);

The curve intersects the line z = -3 at 0 point(s).

% Intersection of polygon and line
figure
plot(x_prime, z_prime);
hold on
yline(-3)
p = polyshape(x_prime, z_prime);
l = [-5 -3; 15 -3];
[in out] = intersect(p, l);
in(:, 1)   % the intersect point

ans = 2×1

1.0354 11.6586

plot(in(:, 1), [-3 -3], 'ro')

[Sarinya]

Thank you so much!

Answer 2

[Sandeep Mishra]

Hello Sarinya,

I understand that you want to find the intersection points of you contour plot for y==-3, But you are using the find operator with plot object, so you are getting an empty array as result.

Also in your code the outlinePlot.YData conatins double data type & it doesn't contain exact -3.0000 value so you will not get your desried result.

You can try the below code to get your desired result.

% Indices of value > -3.01
ind1 = find(outlineplot.YData>-3.01)

ind1 = 1×122

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

% Indices of value < -2.99
ind2 = find(outlineplot.YData<-2.99)

ind2 = 1×82

60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89

% Indices of value > -3.01 and < -2.99
ind = intersect(ind1, ind2)

ind = 1×2

60 141

% X coordinates
xCoord = outlineplot.XData(ind)

xCoord = 1×2

1.0276 11.6709

You can refer to the below documentation to learn more about "find" in MATLAB.

https://in.mathworks.com/help/matlab/ref/find.html