Best fit of a surface

90 views (last 30 days)
Michele
Michele on 31 Jan 2014
Commented: Matt J on 1 Feb 2014
Dear all, i have a matrix A(188x215) but is has many pics as you can see in the image
i would like to visualize a plane that best fits this matrix.
So far i made it with some for cycles but i do think that there can be a better way to solve that.
Any help would be appreciated

Accepted Answer

arich82
arich82 on 31 Jan 2014
If you're just looking for a best fit plane (for some definition of 'best'), you can use the backslash operator, just like finding the the least-squares fit for a line.
For a line, you have data x and y, but unknown coeffcients m abnd b:
y1 = m*x1 + b;
y2 = m*x2 + b;
y3 = m*x3 + b;
...
This represents an overdetermined system X*M = Y, where we can solve M = [b; m]
X = [ ...
1, x1; ...
1, x2; ...
1, x3; ...
];
Y = [y1; y2; y3;];
M = X\Y;
For a plane, we have the general equation
0 = a*(x - x0) + b*(y - y0) + c*(z - z0)
z = (1/c)*(a*x0 + b*y0 + c*z0) - ((a/c)*x + *b/c)*y);
Let C = (1/c)*(a*x0 + b*y0 + c*z0), A = -a/c, B = -b/c, then for M = [C, A, B], you again have the overdetermined equation Z = X*M, which can be solved via backslash:
X = [ ...
1, x1, y1; ...
1, x2, y2; ...
1, x3, y3; ...
];
Z = [z1; z2; z3;];
M = X\Z;
In general, if you have vectors of data x(:), y(:), z(:) of length n, then you can create the data variables as X = [ones(n, 1), x(:), y(:)]. Below is a full example:
m = 188;
n = 215;
[y, x] = meshgrid(1:n, 1:m);
z = sind(10*x + 5*y) + log(x.*y); % filler data; this would be your A matrix
X = [ones(m*n, 1), x(:), y(:)];
M = X\z(:);
figure;
surf(x, y, z, 'EdgeColor', 'none'); % original wavy data
hold all;
surf(x, y, reshape(X*M, m, n), 'EdgeColor', 'none', 'FaceColor', 'b'); % blue fitted-plane
Alternatively, if you just want to get rid of the peaks, you might consider filtering your data (e.g. filter2).
Hope this helps.
--Andy
  4 Comments
arich82
arich82 on 1 Feb 2014
Matt J--
Do you have a reference thread for the QR stability comment? Or a pathological case where mldivide would fail in fitting a simple plane? I'm genuinely interested.
Once upon a time, the Matlab help documentation described mldivide as using QR for non-square matrices, though perhaps there was some option for a more stable decomp that wasn't used by default. (It now appears there's just a brief comment under "Systems of Linear Equations" about mldivide and qr; it would be nice MW didn't gut their help docs...) The mldivide doc no longer offers any real info, and the new example in qr shows a "Q-less" solution to a system of linear equations, but doesn't offer any insights on what might be different with the backslash operator.
Matt J
Matt J on 1 Feb 2014
Once upon a time, the Matlab help documentation described mldivide as using QR for non-square matrices,
It still does,
so I guess it doesn't matter if all you want is the fit. If you wanted an error analysis of the fit, like polyfit provides for example, it helps to do the QR decomp explicitly and reuse it for the error analysis...

Sign in to comment.

More Answers (1)

Matt J
Matt J on 31 Jan 2014
  1 Comment
Michele
Michele on 31 Jan 2014
Thank you for the link! I had a look and there are many interersting infos!

Sign in to comment.

Categories

Find more on Polynomials in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!