an even more stable procedure

This commit is contained in:
Gael Guennebaud 2011-02-03 11:25:34 +01:00
parent 5beb2f4f0d
commit 1eae6d0fb9

View File

@ -153,11 +153,18 @@ void eigen33(const Matrix& mat, Matrix& evecs, Vector& evals)
// tmp.diagonal().array() -= evals(2); // tmp.diagonal().array() -= evals(2);
// evecs.col(2) = tmp.row(0).cross(tmp.row(1)).normalized(); // evecs.col(2) = tmp.row(0).cross(tmp.row(1)).normalized();
// a slighlty more stable version: // a more stable version:
if((evals(2)-evals(0))<=Eigen::NumTraits<Scalar>::epsilon())
{
evecs.setIdentity();
}
else
{
Matrix tmp; Matrix tmp;
tmp = scaledMat; tmp = scaledMat;
tmp.diagonal ().array () -= evals (2); tmp.diagonal ().array () -= evals (2);
evecs.col (2) = tmp.row (0).cross (tmp.row (1)).normalized (); evecs.col (2) = tmp.row (0).cross (tmp.row (1)).normalized ();
tmp = scaledMat; tmp = scaledMat;
tmp.diagonal ().array () -= evals (1); tmp.diagonal ().array () -= evals (1);
evecs.col(1) = tmp.row (0).cross(tmp.row (1)); evecs.col(1) = tmp.row (0).cross(tmp.row (1));
@ -166,7 +173,11 @@ void eigen33(const Matrix& mat, Matrix& evecs, Vector& evals)
evecs.col(1) = evecs.col(2).unitOrthogonal(); evecs.col(1) = evecs.col(2).unitOrthogonal();
else else
evecs.col(1) /= n1; evecs.col(1) /= n1;
// make sure that evecs[1] is orthogonal to evecs[2]
evecs.col(1) = evecs.col(2).cross(evecs.col(1).cross(evecs.col(2))).normalized();
evecs.col(0) = evecs.col(2).cross(evecs.col(1)); evecs.col(0) = evecs.col(2).cross(evecs.col(1));
}
// Rescale back to the original size. // Rescale back to the original size.
evals *= scale; evals *= scale;