Solve plane equation with 3 points and additional condition

I am quite new to matlab and only now some of the basic features, so every new problem cause trouble for me.

Given:

p1 = [x1 y1 z1]
p2 = [x2 y2 z2]
p3 = [x3 y3 z3]

Wanted: I need to find a plane

p = [a b c d]'

which satisfies the followind conditions

1. 
2. 

I found this approach, which builds a linear system

% should d = 1 or d = 0?
[p1 1; p2 1; p3 1] * [a b c d]' = 0

But solving this leads to the zero vector and also it does not incorporate the second condition.

Answer 1

[John D'Errico]

This is actually pretty simple. I think you may be missing how simple it is.

In three dimensions, 3 points determine a plane. True. At least as long as they are not collinear, this is so.

But, as you have written them, the coefficients of the plane you will generate are not set in stone. You can scale them arbitrarily.

That is, if it is true that

a*x + b*y + c*z + d == 0

then it is also true that

a/2*x + b/2*y + c/2*z + d/2 == 0

In fact, we can choose any non-zero constant k, and scale all of the coefficients of that plane, multiplying them all by k, as you wrote it. So first, I'll show you how to compute the basic plane equation coefficients.

Given three points, I'll store them as rows of the array p123. Since I have seen no actual numbers, I'll create them randomly.

format long g
p123 = randn(3,3)

p123 = 3×3

0.76501453055664 1.00436795844351 2.18328686621106 0.630610395275719 -0.480810118977791 -0.284734549172793 -0.287533856243955 0.808317469346715 0.716379900768767

So each row of p123 is a point in R^3. LEARN TO USE ARRAYS IN MATLAB!

How can we define a plane? The solution I like to use is to use the normal vector to the plane, and one point that lies in the plane.

We can arbitrarily choose any of those three points as the point in the plane, or I suppose, i could pick the average of all three points. I'll pich the mean.

P = mean(p123,1)

P = 1×3

0.369363689862801 0.44395843627081 0.871644072602345

Now we can compute the normal vector simple using a cross product, or I can use the function null. It creates a vector of length 3 containing the desired coefficients [a,b;c], so I'll just call it abc.

abc = null(p123 - P)

abc = 3×1

-0.511062968255164 -0.723897858039545 0.463450680875517

(That line of code assumes you have release R2016b or later.)

A nice feature of the use of null here is it AUTOMATICALLY normalizes the vector of coefficients [a,b,c] such that a^2+b^2+c^2==1.

d = -p123(1,:)*abc

d = 
         0.106184625853324

Now we can test to see if these coefficients satisfy your goals.

p123*abc + d

ans = 3×1

1.11022302462516e-16 2.77555756156289e-16 1.11022302462516e-16

That should yield zero, or at least as close as we can come in terms of floating point arithmetic on double precision numbers. I'm satisfied with that result.

norm(abc)

ans = 
                         1

Or if you prefer...

sum(abc.^2)

ans = 
                         1

as I am with the latter.

Finally, I did say we could compute the normal vector using a cross product.

N = cross(p123(2,:) - p123(1,:),p123(3,:) - p123(1,:)).';
N = N./norm(N)

N = 3×1

0.511062968255164 0.723897858039545 -0.463450680875517

So the same vector resulted. There could have been an arbitrary sign swap, but that would not have been relevant. In fact, we did see a sign swap of the coefficients as I computed it using the cross product.

[link]

Wow, thank you so much for this detailed description! I really appreciate the time you put into this to explain the solution.

[link]

The null command does not seems to work for me (R2020b)

mp = mean(points,1)
abc = null(points - mp)

[John D'Errico]

Sorry. I guess I was being lazy. As it turns out, you have a subtle numerical problem, where the null function decided that your matrix was not actually singular, after the mean vector subtraction. This should work just as well, since I can arbitrarily use the first point as the point in the plane.

points = [284.36 -47.258 -752.47; ...
    282.64 -47.382 -752.19; ...
    283.3 -46.523 -751.25];
P = points(1,:);
abc = null(points(2:3,:) - P)

abc = 3×1

0.1548 -0.7810 0.6050

d = -points(2,:)*abc

d = 374.3359

Because I pass a 2x3 array to null, it will not fail for the reason it failed in your case. Does it satisfy the resquirements?

norm(abc)

ans = 1

points*abc + d
ans =
            0
            0
   5.6843e-14

So all three points lie in the plane. That last 5.6e-14 is again floating point trash that we can ignore.

[link]

Thanks! You saved me a lot of time

Answer 2

[Matt J]

You can use my planefit() utility

[abc,d]=planefit([p1;p2;p3].');
d=-d;

function varargout=planefit(xyz)
%Fit 3D plane to given (x,y,z) data
%
%    [a,b,c,d] =  planefit(xyz)   
%     [abc,d]  =  planefit(xyz)
%      hcoeff  =  planefit(xyz)
%           
%IN:        
%           
%    xyz:   3xN input array with coordinates organized in columns
%           [x1,x2,x3...;y1,y2,y3,...;z1,z2,z3,...]
%                 
%OUT:       
%           
%    [a,b,c,d] : Coefficients of fitted plane equation a*x+b*y+c*z=d
%    [abc,d]   : Coefficients in 2-argument form where abc=[a,b,c]
%     hcoeff   : Homogeneous coefficient vector hcoeff=[a,b,c,-d]


 if size(xyz,1)~=3
    error 'Input xyz matrix must be 3xN' 
 end

 xyz=xyz.';
 
 mu=mean(xyz,1);
 
 [~,~,V]=svd(xyz-mu,0);
 
 normal=V(:,end).';
 d=normal*mu';
 
 switch nargout
     
     case {0,1}
 
      varargout={[normal,-d]};
      
     case 2
         
      varargout={ normal, d};
      
     case 4
       
      varargout=[num2cell(normal),{d}]; 
         
     otherwise
         
      error 'Must be 1,2, or 4 output args'
          
 end
 
end

[link]

Thank you very much. I will try this

[link]

So I tested your algorithm and it seems like the second condition ist not always satisfied. Do know how I can include this in the computation of the coefficients?

[Matt J]

So I tested your algorithm and it seems like the second condition ist not always satisfied

For example...?

[link]

Ok, maybe I am wrong. E.g. m is the input argument, p is the result of your algorithm

I check with ~eq(c2, 1) if it is not zero, and if so, I print the line "NOT ONE". But obiously it is 1, maybe with some rounding error. Do you know how to fix that?

Or can I assume, that is always true

[Matt J]

There is always rounding error in a floating point computer... It is something to be accepted, not fixed.

[link]

Yes sure. Sorry, I was imprecise. I wanted to know, if I could then apply a different test for equality. But now know that this is not even necessary. So just ignore my question. Thank you anyway for your help