Weighted Polynomial Surface for 3D Points

Question

Jack Williams 2015-4-4

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/196103-weighted-polynomial-surface-for-3d-points

评论： Chad Greene 2019-11-15

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

John D'Errico 2015-4-4

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/196103-weighted-polynomial-surface-for-3d-points#answer_173961

编辑：John D'Errico 2015-4-4

在 MATLAB Online 中打开

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.

4 个评论
显示 2更早的评论隐藏 2更早的评论

Jack Williams 2015-4-6

编辑：Jack Williams 2015-4-6

John - many thanks indeed for your explanation. This makes sense to me now.

All the best,

Jack

Chad Greene 2019-11-15

在 MATLAB Online 中打开

Note: Since John posted his answer, Matlab now has implicit expansion, meaning

A = (diag(w)*M)\(w.*z);

can be written simply as

A = (w.*M)\(w.*z);

and memory is no longer an issue by this method, even with a long vector of weights.

One point I'm not clear about, though, is whether perhaps it should actually be the square root of weights. I would love an explanation of whether we should or shouldn't take the square root of formal weights when calculating least squares polynomial fits.

请先登录，再进行评论。