Fast plane fitting to many points
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
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