MATLAB Answers

Solve plane equation with 3 points and additional condition

85 views (last 30 days)
link
link on 21 Jan 2021
Commented: link on 21 Jan 2021
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
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.

  0 Comments

Sign in to comment.

Accepted Answer

John D'Errico
John D'Errico on 21 Jan 2021
Edited: John D'Errico on 21 Jan 2021
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.

  4 Comments

Show 1 older comment
link
link on 21 Jan 2021
The null command does not seems to work for me (R2020b)
mp = mean(points,1)
abc = null(points - mp)
John D'Errico
John D'Errico on 21 Jan 2021
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.

Sign in to comment.

More Answers (1)

Matt J
Matt J on 21 Jan 2021
Edited: Matt J on 21 Jan 2021
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

  6 Comments

Show 3 older comments
link
link on 21 Jan 2021
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
Matt J on 21 Jan 2021
There is always rounding error in a floating point computer... It is something to be accepted, not fixed.
link
link on 21 Jan 2021
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

Sign in to comment.

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by