Best fitting a plane (ax + by + cz = 0) where c could be 0

I'm having trouble with some plane fitting in Matlab 2016. I have a set of points I know to fit to a plane, which varies in angle. When using the polynomial function (z = ax + by + c), I get good fitting for most planes, but when the plane is vertical, it has trouble fitting. I have tried the custom equation z = (ax+by+c)/d, but Matlab has struggled to fit the plane, as d would be infinite for a fully vertical plane. Any ideas?

Answer 1

[John D'Errico]

This is an errors in variables problem anyway. So a standard regression will yield biased estimates. The trick is to use an orthogonal regression. You can do that using a SVD.

But first, do you really want to force the plane through the origin? The model that you pose does exactly that. It has no constant term, therefore, the point (x,y,z) = (0,0,0) must lie in the plane. The data that you show does not seem to fit that either, but it is difficult to see since I cannot rotate your plot.

So, ASSUMING that you really wish to fit the model...

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

where a,b,c,d are the unknowns, and the fit is an orthogonal regression, then use the svd as follows:

1. Store your data in the form of an nx3 array. I'll call it xyz.

2. Compute the mean of your data as a 1x3 vector, muxyz. You will need this later.

3. Subtract off the mean vector from each data point. bsxfun does this well.

4. Call svd on the mean subtracted data. The singular vector (in V) corresponding to the smallest singular value contains the coefficients a,b,c.

5. Recover the constant d.

So, assuming you have data in the array xyz (I'll make some up.)

xyz = randn(100,3) * [1 2 3;1 1 1;-1 0 1] + randn(100,3)/5 + 5;

muxyz = mean(xyz,1);
xyzsub = bsxfun(@minus,xyz,muxyz);
[U,S,V] = svd(xyzsub,0);
abc = V(:,3);
d = -dot(muxyz,abc);

abc
abc =
      -0.425079062281321
       0.807539049132476
      -0.408886873029998

d
d =
       0.112055404899681

Note that the vector abc has the property that c (or a or b) MAY be zero, if your data justifies that result.

The magnitude of the singular values will be a measure of how well the data fits a plane. So look at diag(S). If the third element of that vector is significantly smaller than the others, then a plane is a good fit. If the second one is also seriously smaller than the first, then a line was probably a better choice than trying to fit a plane.

diag(S)
ans =
        41.6967917874865
        15.4981087558357
        1.84059008683425

So the third element was small, as we wanted to see. That means the deviations from the regression plane were small.

If your goal really is to force the plane through the origin, then you do not subtract off the mean, in which case d would become zero.

Edit: I might as well add some explanation as to why this works.

Think of the rows of xyz as a set of points in space. You can represent any row as a linear combination of the columns of V. Those diagonal elements of S tell you something about the variability the point set has in the direction of the corresponding basis vector (column of V). The resulting plane chosen by svd will be the one with the SMALLEST possible variability out of plane.

Ok, so if the third singular value is relatively small, then it means the points lie essentially in a plane. That vector V(:,3) will be the normal vector to the plane of interest, but that is how a plane is defined, by the normal vector to the plane. So SVD works very nicely here.

By the way, you will sometimes find people telling you to use principle components to do the orthogonal regression. That is effectively what we did here.

[Alistair Rice]

Many thanks John. I know that the fitted plane should include the point (25,25,25). Would your above explanation cover this?

[Star Strider]

@John — Excellent, clear discussion. Bookmarked it.

+1

[John D'Errico]

(thanks)

If the point [25,25,25] must lie in the plane, then use it as the point to subtract off. Just don't compute the mean.

muxyz = [25 25 25];

Then just do the rest, subtracting the supplied point from your data instead of a computed mean.

To be honest, I would do the computation twice, once by computing the mean as I showed, then with the point you have supplied. Now, compare the singular values. If the smallest singular value is significantly smaller when you computed the mean, then it suggests that [25,25,25] is not actually a point on the plane. So you can at least verify your assertion that this point does lie in the plane. Think of this as a lack of fit test on the assertion that [25,25,25] lies in the plane.