85 views (last 30 days)

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

% 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.

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)

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)

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)

(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

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

p123*abc + d

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)

Or if you prefer...

sum(abc.^2)

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)

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.

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)

d = -points(2,:)*abc

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)

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.

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

Matt J
on 21 Jan 2021

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

Start Hunting!
## 0 Comments

Sign in to comment.