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

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

[William Rose]

@Mark Kamps,

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]

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,

https://www.mathworks.com/matlabcentral/fileexchange/87584-object-oriented-tools-for-fitting-conics-and-quadrics

Answer 1

[Matt J]

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

[Mark Kamps]

@Matt J

Thank you for your reply! I am trying my best to understand what happens and how to interpret the results.

I understand the plane normals now give a proper fit of the data. However, this still means I have to construct the planes on the normal vectors? Since that still gives one free parameters that determines the perpendicular distance to each point? Am I correct?

Furthermore, by performing the above method, both planes are fit simultainously, given an overal best fit for both planes. However. I am specifically looking for a fit on a second plane, based on a prior fit of the first plane.

So; I want to fit a plane to a single data set, given te best perpendicular fit. (@William Rose and @Matt J, you both were obviously right that my best-fit approach was not correct). After I have fitted this plane, I want to fit a second plane on the second data set, with the condition that it is orthogonal to the first plane. Even if this means that the overal fit is sub-optimal compared to your proposed method.

Again, thank you for the help, I really appreciate it!

[Matt J]

Since that still gives one free parameters that determines the perpendicular distance to each point? Am I correct?

Yes, but that's the easy part. Once you have the normals, you can solve for the last parameter just by plugging in a known point from each plane (probably mean(Data1) and mean(Data2)) into the plane equations.

I want to fit a second plane on the second data set, with the condition that it is orthogonal to the first plane. Even if this means that the overal fit is sub-optimal compared to your proposed method.

If you're sure that's what you want, only small modifications are required:

D1=Data1-mean(Data1);
D2=Data2-mean(Data2);

[~,~,V1]=svd(D1,0); 
    z=V1(:,end);  %unconstrained normal to first plane
    
[~,~,V2]=svd(D2,0); 
   x0=V2(:,end);  %unconstrained normal to 2nd plane (initial guess)

fun=@(x) norm(D2*x(:)).^2;
nonlcon=@(x) deal([], [dot(x,z) ; norm(x)^2-1] );

x=fmincon(fun, x0,  [],[],[],[],[],[],nonlcon);

normal1=z/norm(z);
normal2=x/norm(x); %plane normals

[Mark Kamps]

You, sir, are a genius!

Thank you very much for your participation! This works as intented!

[William Rose]

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