Fit a plane on 3D-point data, orthagonal to an existing plane?

21 views (last 30 days)
Dear MATLAB users,
I am stuck on a problem I am working with, and I really hope you guys/gals can help me.
Suppose I have two different sets of 3D point data:
Data1 = rand(100,3) .* [100 100 10];
Data2 = rand(100,3) .* [100 10 100];
figure; hold on; axis equal;
plot3(Data1(:,1),Data1(:,2),Data1(:,3),'.')
plot3(Data2(:,1),Data2(:,2),Data2(:,3),'.')
Now suppose I want to fit a plane through my first dataset:
SurfFit = fit([Data1(:,1),Data1(:,2)], Data1(:,3),'poly11');
This goes according to plan. I can see that MATLAB fits a plane (Poly11: Z = p00 + p10*x + p01*y) through my data. I can see that for the three coefficients it gets certain values:
Linear model Poly11:
val(x,y) = p00 + p10*x + p01*y
Coefficients (with 95% confidence bounds):
p00 = 6.429 (4.915, 7.943)
p10 = -0.01357 (-0.03247, 0.005329)
p01 = -0.005913 (-0.02611, 0.01429)
Obviously, due to randomness, your values might be different..
Now, I want to fit a surface (poly11) through my second dataset, but it should be orthogonal to my first plane. This means that the following statement must be met:
p00(1)*p00(2) + p10(1)*p10(2) + p01(1)*p01(2) = 0
Where (1) refers to the first surface, and (2) to the second surface.
At this point I don't know how to proceed, since I am not that familiar with surface fitting. I know that you can make your own custom equations to fit a surface on, but I do not know how to combine the Poly11 surface with the above equation to make it fit my data. Is it even possible at all? Any help would be appreciated!
Kind regards,
Mark
  2 Comments
William Rose
William Rose on 3 May 2021
As you probably know, the best-fit plane you have found for the first dataset is the plane that minimizes the errors in the z-direction. That s good, if you know the x- and y-values precisely, and only z is subject to measurement error or modeling error. If you want to minimize the perpendicular distance between each point and the plane, you need a different approach.
Let us assume that you like plane 1, which fits data set 1. You define plane 1 by the equation
You are right that perpendicularity (is that a word?) of the second plane means that the normals to the two plane must be orthogonal, i.e. the inner product of the normals is zero. But the equation you gave for acomplishing this is not correct, because it is not true that the normal to plane 1 is (). It's more complicated.
Your second plane will have two degrees of freedom: where it intersects some axis, such as the z-axis, and the direction of the normal, which could be expressed ans an angle from 0 to 2*pi, ltying within or parallel to plane 1.
I don't have time to get farther now, but you mght want to think about what I've posted so far.
Matt J
Matt J on 3 May 2021
If you want to minimize the perpendicular distance between each point and the plane, you need a different approach.
Yes, and you can use planarFit() in this FEX submission package to do so,

Sign in to comment.

Accepted Answer

Matt J
Matt J on 3 May 2021
Edited: Matt J on 3 May 2021
I think if you want to do this properly, the Curve Fitting Toolbox won't be enough. You really need to jointly fit both Data sets using fmincon, using nonlinear constraints.
D1=Data1-mean(Data1);
D2=Data2-mean(Data2);
C=blkdiag(D1,D2);
fun=@(x) norm(C*x(:)).^2;
nonlcon=@(x) deal([], [dot(x(1:3),x(4:6)) ; norm(x)^2-1] );
[~,~,V1]=svd(D1,0); [~,~,V2]=svd(D2,0);
x=fmincon(fun, [V1(:,end);V2(:,end)], [],[],[],[],[],[],nonlcon);
normal1=x(1:3)/norm(x(1:3)); %plane normals
normal2=x(4:6)/norm(x(4:6));
  4 Comments
Mark Kamps
Mark Kamps on 4 May 2021
You, sir, are a genius!
Thank you very much for your participation! This works as intented!
William Rose
William Rose on 4 May 2021
@Mark Kamps, I agree that @Matt J has provided the ideal solution. It uses singular value decomposition to find the best fit plane subject to the constraint that the plane be normal to the first plane. There is a lot of subtlety in his solution.
You may wonder how the plane found by his method relates to the equation you used to charaterize a plane, which was
First: the normal to the plane you have described is
This normal is not a unit vector, but that is OK: the normal does not have to be a unit vector. You can prove that this is in fact the normal, by noting that another way to decribe a plane is that it is the set of points P=[x,y,z] that satisfy
dot(P-P0,n)=0, which is equivalent to saying
where is a point (any point) in the plane, and n=[A,B,C] is the normal. You can re-arrange to show that and .
The solution of @Matt J makes plane 2 go through the centroid of the second set of points. He finds the plane that minimzes the perpendicular errors, subject to the constraint that it be perpendicular to plane 1. The error for each point is the dot product of the unit normal with the vector from the point to the centroid of points. (The centroid is in the plane, by definition.*) Therefore @Matt J must insure that the normal has unit length, and he does, by passing an additional constraint to fmincon(), to make it so.
When the fit for plane 2 is done, you have P0=mean(Data2)=[x0,y0,z0], and you have normal2=[A,B,C], which has unit length. Therefore
and
which has the form of the equation for the plane in your original post.
*It is possible that you could get a slightly better fit for plane 2 by not forcing it to go through the centroid of Data2. In the case where the plane is not constrained to be normal to another plane, the best-fit plane goes exactly thorugh the centroid.

Sign in to comment.

More Answers (0)

Products


Release

R2019a

Community Treasure Hunt

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

Start Hunting!