1

I need to implement the algorithm described below. Everything is fine until the eigenvalues computation. I'm completely new to them and I found a lot of very complicated paper on the net. Is it possible that this 'best-fitting plane' case requires such a complex approach? Where can I find straightforward implementation of a suitable eigenvalue solver for this?

I am also interested in different approaches to the 'best-fitting plane' problem.

Thanks.

xm = mean(x); ym = mean(y); zm = mean(z);

u = x - xm; v = y - ym; w = z - zm;

M = [sum(u.^2),sum(u.*v),sum(u.*w); ...
     sum(v.*u),sum(v.^2),sum(v.*w); ...
     sum(w.*u),sum(w.*v),sum(w.^2)];

[V,D] = eig(M);

Now choose the smallest eigenvalue in D. The corresponding eigenvector in V gives the coefficients, v1, v2, v3, in the best-fitting plane:

v1*(X-xm) + v2*(Y-ym) + v3*(Z-zm) = 0

By the way, if that smallest eigenvalue should turn out to be zero, that is your indication that the four points are coplanar and your plane is an exact fit.

abenci
  • 185
  • Yes, this is the correct approach for the best-fitting plane. To get useful responses where to find a suitable eigenvalue solver, you should probably tell us more about your computing environment and requirements. – joriki Nov 14 '11 at 09:43
  • 1
    Bleh, I don't recommend this method, though. Better to use SVD instead of the eigendecomposition for these matters... – J. M. ain't a mathematician Nov 14 '11 at 09:52
  • @J.M.: I don't understand -- since $M$ is symmetric, isn't that the same thing? – joriki Nov 14 '11 at 09:58
  • @joriki: theoretically, yes. Computationally, a bad idea, since the cross-product matrix's condition number is the square of the condition number of the design matrix. See this previous answer of mine, where I show a simple example due to Läuchli. – J. M. ain't a mathematician Nov 14 '11 at 10:56
  • @J.M.: Ah, I see, sorry -- I thought you were suggesting to apply SVD instead of eigendecomposition to the covariance matrix :-) – joriki Nov 14 '11 at 11:15
  • I have a SVD matrix solver inside our product, can SVD help to solve this problem? We use it for solving certain Ax = b system of equations. – abenci Nov 14 '11 at 13:16
  • @joriky: Everything readable as C/C++ would be fine. – abenci Nov 14 '11 at 15:02
  • @devdept: I don't get pinged if you misspell my name :-) Yes, as J. M. rightly pointed out, this can be done using SVD. You can follow J. M.'s link; if you don't understand the answer, this Wikipedia section might help. – joriki Nov 14 '11 at 19:07
  • @joriki: Sorry ;-). I checked my SVD implementation and accepts two matrices or one matrix and one vector as input. The problem above only prepares the M matrix, what can I do? I also found an internal method called void svdcmp(double[,] a, out double[] w, out double[,] v), maybe this could help? – abenci Nov 16 '11 at 16:00

0 Answers0