How do I extract the intersection line between a plane and a surface?

6 次查看（过去 30 天）

burges 2018-12-22

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/436917-how-do-i-extract-the-intersection-line-between-a-plane-and-a-surface

编辑： burges 2018-12-22

intersection line.png

I have 2 sets of 3D points through which I fitted a plane and a surface (see attached figure). I tried to extract the line of intersection using symbolic expressions as described in: https://uk.mathworks.com/matlabcentral/answers/301302-intersection-of-two-sfit-planes However, I am unable to visualize the intersection line as I get the error “Error using inlineeval (line 14) Error in inline expression… “. More importantly, I’m stuck as to how to get the coordinates of x,y,z points along the intersection line.

clc; clear; close all;
       
Xs = [18.5, 18, 17.5
      18.5, 18, 16.5
      19.5, 18, 17.5
      19.5, 18, 16.5
      17.5, 18, 18.5
      16.5, 18, 18.5];
 
Xn = [18.5, 18.5, 18.0
      18.5, 19.5, 18.0
      19.5, 18.5, 18.0
      19.5, 19.0, 17.5
      19.5, 19.5, 17.0
      18.0, 18.5, 18.5
      18.0, 19.5, 18.5];
  
% fit plane through Xs
[sfp,gofp,obj] = fit([Xs(:,1), Xs(:,3)],Xs(:,2), 'poly11');
plot(sfp,[Xs(:,1), Xs(:,3)],Xs(:,2))
hold on;
 
% fit surface through Xn
[sfs,gofs,obj] = fit([Xn(:,1), Xn(:,3)],Xn(:,2), 'poly22');
plot(sfs,[Xn(:,1), Xn(:,3)],Xn(:,2))
 
% obtain coefficients for yplane and ysurface
g = sfp.p10;
h = sfp.p01;
k = sfp.p00;
a = sfs.p20;
b = sfs.p02;
c = sfs.p10;
d = sfs.p01;
e = sfs.p11;
f = sfs.p00;
 
syms x z
yp = g*x + h*z + k;
ys = a*(x^2) + b*(z^2) + c*x + d*z + e*x*z + f;
 
% Solve for z in the expression yp==ys. This gives you z(x). 
z_of_x = solve(yp==ys, z);
 
% Noow substitute the value of z(x) in yp to obtain an expression for y(x)
y_of_x = simplify(subs(yp, z, z_of_x));
 
% Plot the result
figure
h1 = ezsurf(yp,[-20 20]);
set(h1,'facecolor','r','facealpha',0.5);
hold  all
h2 = ezsurf(ys,[-20 20]);
set(h2,'facecolor','b','facealpha',0.5);
hp = ezplot3(x, z_of_x, y_of_x, [-2 2]);
set(hp,'LineWidth',2,'Color','g');
ylim([-2 2]);