Intersection points of generated lines and polygon

I have a script (below) that generates two plots: a set of points that represent a cross_section of a 3D surface (blue) and a set of fixed, set angle lines (red)

vv = Cross_sections;
cla
hold on
for k = 1:1 %length(Cross_sections)   
    x = [vv{k}(:,1); vv{k}(1,1)];
    y = [vv{k}(:,2); vv{k}(1,2)];
    z = [vv{k}(:,3); vv{k}(1,3)];
    plot3(x,y,z,'.-b');
   
    N = length(x);
        % use 3 points to calculate normal vector of crossection plane
    i1 = 1;                 % first index            
    i2 = round(N/3);        % second index
    i3 = round(2*N/3);      % third index
    v1 = [x(i1)-x(i2) y(i1)-y(i2) z(i1)-z(i2)];
    v2 = [x(i1)-x(i3) y(i1)-y(i3) z(i1)-z(i3)];
    v = cross(v1,v2);       % normal vector of a plane
    a = v(1); b = v(2); c = v(3);
    
        % use normal vector to find approximate centroid
	x0 = (max(x)+min(x))/2;     % approximate centroid X
    y0 = (max(y)+min(y))/2;     % approximate centroid Y
        
    % use a(px-x0) + b(py-y0) + c(pz-z0) == 0 
        % equation to find centroid Z
        
        % choose first point for (px, py, pz)
    z0 = (a*(x(1)-x0) + b*(y(1)-y0))/c + z(1);
    
   
    % circle
    R = 10;
    ac  = sqrt(a*a+c*c);
    abc = sqrt(a*a+b*b+c*c);
    t1 = linspace(0,2*pi,62);
    ct = cos(t1);
    st = sin(t1);
    x1 = x0 + R/ac*(c*ct-a*b*st/abc);
    y1 = y0 + R*ac/abc*st;
    z1 = z0 - R/ac*(a*ct+b*c*st/abc); 

        % draw lines
    for i = 1:length(x1)
        plot3([x0 x1(i)],[y0 y1(i)],[z0 z1(i)],'.-r')
    end
end
hold off

Question: Could someone please give me some guidance so as to generate or incorporate a script into MY script in order to generate the intersection points of the blue lines with the red lines. Perhaps the intersection points to be stored and also plotted on the same figure as an asterix or so.

Thank you!!

[Arvin Rasoulzadeh]

I noticed that your "pencil of lines" and your "blue curve" are planar. Therefore, you can find the carrier plane of them and with a suitable rotation (and translation) make it coincide with xy-plane. From then you can use polyxpoly command in the following way:

for i = 1:length(x1)

[p1{i},p2{i}] = polyxpoly([x0, x1(i)],[y0,y1(i)],x,y);

end

pp1 = cell2mat(p1);

pp2 = cell2mat(p2);

hold on

axis equal

plot(pp1',pp2',zeros(numel(pp1),1),'g-');

hold off

Which the result will be:

Answer 1

[KSSV]

The best thing for you today is this function: https://in.mathworks.com/matlabcentral/fileexchange/22441-curve-intersections?focused=5165138&tab=function

[Devrim Tugberk]

I tried this option and my P always ends up being empty suggesting that there is no interseciton but clearly there is. Maybe I am using the arguments that are suppose to go in the function wrong? Some guidance maybe? Thanks

[KSSV]

This will definitely work.......need the data to help your further.

[Devrim Tugberk]

Here it is, thank you i appreciate it

[KSSV]

load('cross_secitions.mat') ;

vv = Cross_sections;
cla
hold on
for k = 1:1 %length(Cross_sections)
    x = [vv{k}(:,1); vv{k}(1,1)];
    y = [vv{k}(:,2); vv{k}(1,2)];
    z = [vv{k}(:,3); vv{k}(1,3)];
    plot3(x,y,z,'-b');
    L1 = [x' ; y' ;z'] ; 
    N = length(x);
    % use 3 points to calculate normal vector of crossection plane
    i1 = 1;                 % first index
    i2 = round(N/3);        % second index
    i3 = round(2*N/3);      % third index
    v1 = [x(i1)-x(i2) y(i1)-y(i2) z(i1)-z(i2)];
    v2 = [x(i1)-x(i3) y(i1)-y(i3) z(i1)-z(i3)];
    v = cross(v1,v2);       % normal vector of a plane
    a = v(1); b = v(2); c = v(3);
    
    % use normal vector to find approximate centroid
    x0 = (max(x)+min(x))/2;     % approximate centroid X
    y0 = (max(y)+min(y))/2;     % approximate centroid Y
    
    % use a(px-x0) + b(py-y0) + c(pz-z0) == 0
    % equation to find centroid Z
    
    % choose first point for (px, py, pz)
    z0 = (a*(x(1)-x0) + b*(y(1)-y0))/c + z(1);
    
    
    % circle
    R = 10;
    ac  = sqrt(a*a+c*c);
    abc = sqrt(a*a+b*b+c*c);
    t1 = linspace(0,2*pi,62);
    ct = cos(t1);
    st = sin(t1);
    x1 = x0 + R/ac*(c*ct-a*b*st/abc);
    y1 = y0 + R*ac/abc*st;
    z1 = z0 - R/ac*(a*ct+b*c*st/abc);
    % draw lines
    P = zeros([],2) ;
    for i = 1:length(x1)
        plot3([x0 x1(i)],[y0 y1(i)],[z0 z1(i)],'.-r')
        L2= [x0 x1(i) ; y0 y1(i) ; z0 z1(i)] ; 
        % GEt Intersection
        P(i,:) = InterX(L1,L2)' ; 
    end
end
% GEt Z for intersection points 
F = scatteredInterpolant(x,y,z) ; 
P(:,3) = F(P(:,1),P(:,2)) ; 
plot3(P(:,1),P(:,2),P(:,3),'*k')
hold off

[ARUN]

Hello,

I am facing similar kinda problem, trying to use InterX but not able get the intersection point as expected.

Below is my case. I need to find the intersection points of red curve on each of the horizontal greenlines.

# Plotting code(datapoints atttached)

P1 = [Xin,Yin];

P2 = [Xout,Yout];

Xa = [P1(:,1),P2(:,1)];

Xb = [P1(:,2),P2(:,2)];

Xa1 = Xa';

Xb1 = Xb';

X1 = X';

Y1 = Y';

plot(Xa1,Xb1,'g');

hold on

plot(X1,Y1)

Would be of great help if you can suggest me on fixing this problem. Thank you !