Fast plane fitting to many points

Question

I'm looking to fit a plane to a set of ~ 6-10k 3D points. I'm looking to do this as fast as possible, and accuracy is not the highest concern (frankly the plane can be off by +-10 degrees in any of the cardinal axes).

My current approach is to use best of best fit, but it's incredibly slow (I'm hoping to extract planes at a rate of about 10-50k times each time I run the algorithm, and at this rate it would finish in weeks, as opposed to hours) as it works on all possible combinations of 6000 points, so ~35,000,000,000 iterations, and frankly it has a much higher accuracy than I need.

Does anybody know of any weaker plane-fitting techniques that might speed my algorithm up considerably?

EDIT:

I've managed to get the number of iterations down to ~42k by creating planes at each possible 3D angle (stepping through at 5 degrees each time) and testing the existing points against these to find the best plane, instead of fitting planes to the points I have.

I'm sure there's something to be gained here by divide and conquering too, although I worry I could jump straight past the best plane.

Answer 1

Use the standard plane equation Ax + By + Cz + D = 0, and write the equation as a matrix multiplication. P is your unknown 4x1 [A;B;C;D]

g = [x y z 1];  % represent a point as an augmented row vector
g*P = 0;        % this point is on the plane

Now expand this to all your actual points, an Nx4 matrix G. The result is no longer exactly 0, it's the error you're trying to minimize.

G*P = E;   % E is a Nx1 vector

So what you want is the closest vector to the null-space of G, which can be found from the SVD. Let's test:

% Generate some test data
A = 2;
B = 3;
C = 2.5;
D = -1;

G = 10*rand(100, 2);  % x and y test points
% compute z from plane, add noise (zero-mean!)
G(:,3) = -(A*G(:,1) + B*G(:,2) + D) / C + 0.1*randn(100,1);

G(:,4) = ones(100,1);   % augment your matrix

[u s v] = svd(G, 0);
P = v(:,4);             % Last column is your plane equation

OK, remember that P can vary by a scalar. So just to show that we match:

scalar = 2*P./P(1);
P./scalar

ans = 2.0000 3.0038 2.5037 -0.9997

Answer 2

In computer vision a standard way is to use RANSAC or MSAC, in your case;

1. Take 3 random points from the population
2. Calculate the plane defined by the 3 points
3. Sum the errors (distance to plane) for all of the points to that plane.
4. Keep the 3 points that show the smallest sum of errors (and fall within a threshold).
5. Repeat N iterations (see RANSAC theory to choose N, may I suggest 50?)

http://en.wikipedia.org/wiki/RANSAC

Answer 3

It looks like griddata might be what you want. The link has an example in it.

If this doesn't work, maybe check out gridfit on the MATLAB File Exchange. It's made to match a more general case than griddata.

You probably don't want to be rolling your own surface fitting, as there's several well-documented tools out there.

Take the example from griddata:

x = % some values 
y = % some values
z = % function values to fit to

ti = % this range should probably be greater than or equal to your x,y test values
[xq,yq] = meshgrid(ti,ti);
zq = griddata(x,y,z,xq,yq,'linear'); % NOTE: linear will fit to a plane!
Plot the gridded data along with the scattered data.

mesh(xq,yq,zq), hold
plot3(x,y,z,'o'), hold off