Plot a plane from 3 points

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];

Answer 1

[Mathieu NOE]

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

[Alberto Acri]

as always, perfect!

[Mathieu NOE]

My pleasure !

[Alberto Acri]

@Mathieu NOE

May I ask how can I visualize the plane passing through these 3 points? It is a plane parallel to the YZ plane.

A = [12.8, 248.2111, 14];
B = [12.8, 48.2, 14];
C = [12.8, 48.2, -60];

Thank you!

[Mathieu NOE]

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


end

[Alberto Acri]

Thank you!

[Mathieu NOE]

again, my pleasure !

Answer 2

[Luca Ferro]

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.