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.

Egestion

Do you have access to the Curve Fitting Toolbox? – Proteose

Unfortunately I don't, I'm stuck with vanilla MATLAB, although I have a lot of programming experience in general so I should be able to handle a fairly complex algorithm. – Egestion

If accuracy isn't your main concern, try reducing the input complexity of your data. Run kmeans or something on the initial set of 6-10k points, and then fit the plane to the exemplars. – Distributary

@Ansari: good idea. To improve performance even further, taking a random subset of points might be enough. – Biconcave

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

Roxie

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

Marji

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

Proteose

Thanks a lot, I'll look into this right away. I'm typically more from a CS background, so my surface-fitting math is a little behind. As such I'm more than happy to let somebody else's code do the job – Egestion

Hmm my problem with griddata is that in order to get it to provide me with a plane, I have to (based off of their first example) tell it to generate zq using 4 points at (-1,-1), (-1,1), (1,-1), (1,1) (using 2:-2 - the bounds of the data set in the example - for some reason only returns NaN). Unfortunately this seems to guarantee that the corners of the plane will be at (-1,-1), (-1,1), (1,-1), (1,1) and it doesn't appear to be taking any more points into consideration. If I increase the number of points, I no longer get a plane. – Egestion

@NickUdell Put some code you've tried in your answer so I can help you more. – Proteose

Answer 4

You may try the consolidator by John D'Errico. It aggregates the points within a given tolerance, this will allow to reduce the amount of data and increase the speed. You can also check John's gridfit function which is usually faster and more flexible than griddata

Antonomasia

Answer 5

If there are no outliers, just a little error in each point's location, then you can pick any 3 widely spaced points to define the plane. This is how I choose them to ensure they're far from colinear.

1. Find the 6 points with the maximum and minimum coordinate in each axis. This is a simple search through all points similar to finding an axis-aligned bounding box (AABB).

2. Identify which of those 3 pairs of points has the greatest separation in the direction of their axis (x, y or z) by subtracting their coordinates. This pair of points defines a line on the plane.

3. Find a 3rd point that's the furthest from this line. Since you don't require the absolute furthest point but just any far-away one, you don't have to do a full-blown 3D line-point-distance calculation but can project each point onto the plane normal to the long axis (x, y or z) and calculate its 2D distance from the line from 2).

Confiscatory