Plot a plane from 3 points

Question

Alberto Acri 2023-3-15

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1929300-plot-a-plane-from-3-points

评论： Mathieu NOE 2023-3-16

Hi, I want to create a plane that goes through points A, B, C. I am using the "plot_line" function found on "file exchange" below:

%% ====== FUNCTION ======
function [normal, d, X, Y, Z] = plot_line_deviazione(p1, p2, p3, x_min, x_max, y_min, y_max)
% This function plots a line from three points. 
% I/P arguments: 
%   p1, p2, p3 eg, p1 = [x y z]
% 
%
% O/P is: 
% normal: it contains a,b,c coeff , normal = [a b c]
% d : coeff
normal = cross(p1 - p2, p1 - p3);
d = p1(1)*normal(1) + p1(2)*normal(2) + p1(3)*normal(3);
d = -d;
x = x_min:x_max; y = y_min:y_max;
[X,Y] = meshgrid(x,y);
Z = (-d - (normal(1)*X) - (normal(2)*Y))/normal(3);
mesh(X,Y,Z)

Can anyone tell me why it is not shown in the plotting? I have this code:

A = [91.0049, 222.2614, -12.2419];
B = [12.8, 48.2, -5.6631];
C = [12.8, 48.2, -53.7279];
figure
plot3(A(:,1),A(:,2),A(:,3),'r.','Markersize',25)
hold on
plot3(B(:,1),B(:,2),B(:,3),'r.','Markersize',25)
plot3(C(:,1),C(:,2),C(:,3),'r.','Markersize',25)
[normal0, d0, X0, Y0, Z0] = plot_line(A, B, C, -100, 100, -100, 100);  % new plane
hold off
grid off
xlabel('x')
ylabel('y')
zlabel('z')
view([15,50,30])
axis equal

NOTE: using these three points shows the plane

A = [-62.6412, 32.6849, -195.9077];
B = [-36.5, -33.1755, -2];
C = [-36.5, 34.2726, -2];

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Mathieu NOE 2023-3-16

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1929300-plot-a-plane-from-3-points#answer_1194145

在 MATLAB Online 中打开

hello Alberto

welcome back

I was first a bit puzzled because you speak of "line" whereas we are talking here "plane" that goes through 3 points. just a minor comment

also I don't see wich Fex function you refer to (missing link) but again this is not important

more important :

1 in your code you don't call the function plot_line_deviazione that you posted but this one : plot_line.

[normal0, d0, X0, Y0, Z0] = plot_line(A, B, C, -100, 100, -100, 100);  % new plane

I supposed that it's actually plot_line_deviazione that you wanted to use in your main code

now let's see what happen when you test your code with the following data (representing a triangle in the vertical plane)

A = [91.0049, 222.2614, -12.2419];

B = [12.8, 48.2, -5.6631];

C = [12.8, 48.2, -53.7279];

the cross product is a vector in the horizontal plane , so it's z coordinate is zero

normal = 1.0e+03 *( 8.3662 -3.7589 0)

therefore in this line you are making a division by zero resulting in Z being Inf :

Z = (-d - (normal(1)*X) - (normal(2)*Y))/normal(3);

so you have to check if your cross product has one zero value and make sure you are not making a division by zero when you generate the mesh.

for this specific case you have to create first the x and z plan mesh and then compute the Y result (see below new code of your function)

all the best

% 3 points on vertical plane
A = [91.0049, 222.2614, -12.2419];
B = [12.8, 48.2, -5.6631];
C = [12.8, 48.2, -53.7279];
% % 3 points on inclined plane
% A = [-62.6412, 32.6849, -195.9077];
% B = [-36.5, -33.1755, -2];
% C = [-36.5, 34.2726, -2];
figure(1)
plot3(A(:,1),A(:,2),A(:,3),'r.','Markersize',25)
hold on
plot3(B(:,1),B(:,2),B(:,3),'r.','Markersize',25)
plot3(C(:,1),C(:,2),C(:,3),'r.','Markersize',25)
[normal0, d0, X, Y, Z] = plot_line_deviazione(A, B, C, -100, 100, -100, 100);  % new plane
mesh(X,Y,Z)
hold off
grid off
xlabel('x')
ylabel('y')
zlabel('z')
view([15,50,30])
axis equal
%% ====== FUNCTION ======
function [normal, d, X, Y, Z] = plot_line_deviazione(p1, p2, p3, x_min, x_max, y_min, y_max)
% This function plots a line from three points. 
% I/P arguments: 
%   p1, p2, p3 eg, p1 = [x y z]
% 
%
% O/P is: 
% normal: it contains a,b,c coeff , normal = [a b c]
% d : coeff
normal = cross(p1 - p2, p1 - p3);
d = p1(1)*normal(1) + p1(2)*normal(2) + p1(3)*normal(3);
d = -d;
if abs(normal(3))>eps % cross product is non zero in Z direction ( your original code)
    x = x_min:x_max;
    y = y_min:y_max;
    [X,Y] = meshgrid(x,y);
    Z = (-d - (normal(1)*X) - (normal(2)*Y))/(normal(3));
else % cross product is zero in Z direction => modified code
    x = x_min:x_max;
    z = min([p1(3) p2(3) p3(3)]):max([p1(3) p2(3) p3(3)]);
    [X,Z] = meshgrid(x,z);
    Y = (-d - (normal(1)*X) - (normal(3)*Z))/(normal(2));
end
end

6 个评论
显示 4更早的评论隐藏 4更早的评论

Mathieu NOE 2023-3-16

在 MATLAB Online 中打开

ok

this code should be covering all cases now (I knew the previous release was not !)

I made quite a few changes...

% plane parallel to the YZ plane. % tested OK

A = [12.8, 248.2111, 14];

B = [12.8, 48.2, 14];

C = [12.8, 48.2, -60];

% % 3 points on vertical plane % tested OK

% A = [91.0049, 222.2614, -12.2419];

% B = [12.8, 48.2, -5.6631];

% C = [12.8, 48.2, -53.7279];

% % 3 points on inclined plane % tested OK

% A = [-62.6412, 32.6849, -195.9077];

% B = [-36.5, -33.1755, -2];

% C = [-36.5, 34.2726, -2];

x_min = floor(min([A(1) B(1) C(1)]))-10;

x_max = ceil(max([A(1) B(1) C(1)]))+10;

y_min = floor(min([A(2) B(2) C(2)]))-10;

y_max = ceil(max([A(2) B(2) C(2)]))+10;

z_min = floor(min([A(3) B(3) C(3)]))-10;

z_max = ceil(max([A(3) B(3) C(3)]))+10;

[normal0, d0, X, Y, Z] = plot_line_deviazione(A, B, C, x_min, x_max, y_min, y_max, z_min, z_max); % new plane

figure(1)

mesh(X,Y,Z);

hold on

plot3(A(:,1),A(:,2),A(:,3),'r.','Markersize',25)

plot3(B(:,1),B(:,2),B(:,3),'r.','Markersize',25)

plot3(C(:,1),C(:,2),C(:,3),'r.','Markersize',25)

axis equal

view([15,50,30])

xlim([x_min x_max]);

ylim([y_min y_max]);

zlim([z_min z_max]);

hold off

grid off

xlabel('x')

ylabel('y')

zlabel('z')

%% ====== FUNCTION ======

function [normal, d, X, Y, Z] = plot_line_deviazione(p1, p2, p3, x_min, x_max, y_min, y_max, z_min, z_max)

% This function plots a plane going through three points.

% I/P arguments:

% p1, p2, p3 eg, p1 = [x y z]

%

% O/P is:

% normal: it contains a,b,c coeff , normal = [a b c]

% d : coeff

normal = cross(p1 - p2, p1 - p3);

n1 = normal(1);

n2 = normal(2);

n3 = normal(3);

d = p1(1)*n1 + p1(2)*n2 + p1(3)*n3;

% d = -d; % sign has been changed in the lines below

x = x_min:x_max;

y = y_min:y_max;

z = z_min:z_max;

% case a)

if abs(n1) > eps

[Y,Z] = meshgrid(y,z);

X = (d - (n2*Y) - (n3*Z))/n1;

end

% case b)

if abs(n2) > eps

[X,Z] = meshgrid(x,z);

Y = (d - (n1*X) - (n3*Z))/n2;

end

% case c)

if abs(n3) > eps

[X,Y] = meshgrid(x,y);

Z = (d - (n1*X) - (n2*Y))/n3;

end

Alberto Acri 2023-3-16

Thank you!

Mathieu NOE 2023-3-16

again, my pleasure !

请先登录，再进行评论。

Answer 2

Luca Ferro 2023-3-16

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1929300-plot-a-plane-from-3-points#answer_1194130

编辑：Luca Ferro 2023-3-16

basically after debugging i got this results:

mesh(X,Y,Z) fails because Z is a 201x201 whose elements are ALL Inf.

This is caused by its declaration Z = (-d - (normal(1)*X) - (normal(2)*Y))/normal(3) in which the last term normale(3) is 0.

The last term is 0 since in cross(p1 - p2, p1 - p3), the two differences inside are 78.2049 174.0614 41.4860 and 78.2049 174.0614 -6.5788. As you can see, using cross product definition of ||a|| ||b|| sin(theta), the theta is 0 due to their alignment, resulting in the 0.

The latter 3 points work because the normal vector has a non zero third element, 0.1763 to be exact, so using it to divide is not problematic.