MATLAB求解抛物面和多条线段的交点

Question

MATLAB求解 抛物面 和 多条线段 的交点，三维图，含图像。
大佬们求代码和讲解。

Answer 1

同学你好，先看求交点的原理：

然后求这个参数里t的符号解算代码：

syms alpha beta gamma 
syms a b c x1 y1 z1 t
eq = z1 + gamma*t == ...
    a*(x1+alpha*t)^2 + ...
    +b*(y1+beta*t)^2 + c;
t = solve(eq,t)

给出求交点的函数

function [x0,y0,z0] = crossPoint(a,b,c,x1,y1,z1,x2,y2,z2)
% 输入参数abc为抛物面的三个参数
% x1,y1,z1和x2,y2,z2分别是点P1和点P2的坐标【P1和P2构成线段】
% 输出参数x0,y0,z0分别是求得的抛物面和线段的交点坐标
% 若无交点：输出x0，y0，z0为空
% 有一个交点则输出x0，y0，z0各有一个值
% 若两个交点，那么x0，y0，z0各有两个值
alpha =  x2 - x1;%计算alpha、beta、gamma
beta  =  y2 - y1;
gamma =  z2 - z1;
p = gamma^2-4*a*alpha^2*c-4*b*beta^2*c...
    +4*a*alpha^2*z1 + 4*b*beta^2*z1...
    -4*a*alpha*gamma*x1 - 4*b*beta*gamma*y1...
    -4*a*b*beta^2*x1^2 - 4*a*alpha^2*b*y1^2+8*a*b*alpha*beta*x1*y1;% 计算sigma1右边根号里的内容
x0 = []; y0=[]; z0 = []; % 首先取空值
t1 = 2; t2=2; % 这里随便给了t的初始值
if(p==0) % 根号里的内容为0，说明t有一个解，一个交点
    sig = sqrt(p);
    t1 =  -(sig-gamma+2*a*alpha*x1+2*b*beta*y1)/(2*(a*alpha^2+b*beta^2));
elseif(p>0) %根号里的内容大于0，说明t有两个解，两个交点
    sig = sqrt(p);
    t1 =   -(sig-gamma+2*a*alpha*x1+2*b*beta*y1)/(2*(a*alpha^2+b*beta^2));
    t2 =  -(-sig-gamma+2*a*alpha*x1+2*b*beta*y1)/(2*(a*alpha^2+b*beta^2));
end
if(t1<=1 && t1>=0) % 对t1进行区间判断，如果落在0-1之间，就可以
    x0 = [x0; x1 + t1*alpha];
    y0 = [y0; y1 + t1*beta];
    z0 = [z0; z1 + t1*gamma];
end
if(t2<=1 && t2>=0) % 对t2进行区间判断，如果落在0-1之间，就可以
    x0 = [x0; x1 + t2*alpha];
    y0 = [y0; y1 + t2*beta];
    z0 = [z0; z1 + t2*gamma];
end
end

为了调用这个函数，我们假设了抛物面的参数和随机的线段参数，设置了10条随机线段跟这个抛物面进行交点求解，主函数如下：

% 多个相交情况
% 首先要给出 a b c 和多条线段的 x1,y1,z1 以及 x2,y2,z2
N = 10; % 10条线段,我这个例子取随机线段x和y在-4到4之间取值，z在0到5之间取值
x1 = rand(N,1)*8-4;
y1 = rand(N,1)*8-4;
z1 = rand(N,1)*5;
x2 = rand(N,1)*8-4;
y2 = rand(N,1)*8-4;
z2 = rand(N,1)*5;
a = 1/8; b = 1/8; c=1; % 给出抛物面参数
f = @(x,y) a*x.^2+b*y.^2+c; % 这个是抛物面的方程
[X, Y] = meshgrid(linspace(-4,4,101));%准备画抛物面的数据
Z = f(X,Y); %准备画抛物面的数据

figure(1);clf
mesh(X,Y,Z); % 画抛物面
for i = 1:1:N
hold on
plot3([x1(i),x2(i)],[y1(i),y2(i)],[z1(i),z2(i)],'r-o','markerfacecolor','r','markersize',5); %画线段 红色
[x0,y0,z0] = crossPoint(a,b,c,x1(i),y1(i),z1(i),x2(i),y2(i),z2(i));% 调用函数求交点
plot3(x0,y0,z0,'bo','markerfacecolor','b','markersize',5)% 画交点 蓝色
axis equal
end

你的N可以设置为10，这个时候画出来图像如下（因为是随机的，所以不唯一，仅供参考）：

你的N可以设置为100，这个时候画出来图像如下（不唯一，仅供参考）：

当然，自己也可以设置一些坐标数据来画交点。希望能够帮助你！