MATLAB: Weighted Polynomial Surface for 3D Points

Question

Hi

I have an mx3 array of point cloud data (X,Y,Z) and a vector of weights for each point (mx1). I'm trying to create a Polynomial Surface (poly22) of the form:

f(x,y) = p00 + p10*x + p01*y + p20*x^2 + p11*x*y + p02*y^2

From this I need the value of the first coefficient, p00. At the moment I'm using the fit command as follows:

x = XYZ(:,1); % x-column

y = XYZ(:,2); % y-column

Z = XYZ(:,3); % elevation

fitsurface = fit([x,y],z,fittype('poly22'),'Weight',weight);

p00 = fitsurface(0,0);

This works well but is rather slow, considering that this is part of a function I'm working on which acts one-by-one on a point cloud of nearly 2 million points. If the points had no weight, I would use something like:

A = [ones(size(x)) x y x.^2 x.*y y.^2] \ z;

p00 = A(1);

This is much quicker than the 'fit' command; however, I have no idea how to make this work with weights involved. I have also tried the polyfitn function from the File Exchange – around 10% quicker for my dataset.

Any help finding a faster solution to this would be greatly appreciated.

Cheers, Jack

Answer 1

Weighted fitting is simpler than it might seem. Of course, since I wrote polyfitn, why would you EVER want to use anything else? :)

% I'll pick some random weights here. Your weights
% probably make more sense than mine.
W = rand(size(x));
M = [ones(size(x)) x y x.^2 x.*y y.^2];
A = (diag(w)*M)\(w.*z);
p00 = A(1);

The idea is you simply multiply every line of the least squares problem by the corresponding weight. That scales the i'th residual by w(i).

I used diag to build a matrix to scale the rows of M there. If you had a HUGE number of points, that multiply will be less efficient. If you wanted speed, it would be better to use spdiags to build the diagonal matrix with the weights as sparse diagonal. Or, I could have used bsxfun here.

This will in fact be the fastest way to compute the vector of coefficients A, I think:

A = bsxfun(@times,w,M)\(w.*z);

It will certainly be faster than polyfitn, because this does no error checking, no model building, no computation of statistics on the coefficients, or things like R^2.