Plot a plane from 3 points

108 vues (au cours des 30 derniers jours)
Alberto Acri
Alberto Acri le 15 Mar 2023
Commenté : Mathieu NOE le 16 Mar 2023
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];

Réponse acceptée

Mathieu NOE
Mathieu NOE le 16 Mar 2023
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 commentaires
Alberto Acri
Alberto Acri le 16 Mar 2023
Thank you!
Mathieu NOE
Mathieu NOE le 16 Mar 2023
again, my pleasure !

Connectez-vous pour commenter.

Plus de réponses (1)

Luca Ferro
Luca Ferro le 16 Mar 2023
Modifié(e) : Luca Ferro le 16 Mar 2023
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.

Produits


Version

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by