DirectXMath
DirectXMath copied to clipboard
microsoft
BoundingOrientedBox::CreateFromPoints solver suffers from numeric problems
Report from Dave Eberly about problems with the oriented box SolveCubic helper function.
The SolveCubic is called to compute the eigenvalues of a real-valued symmetric 3x3 matrix. Such matrices ALWAYS have 3 real eigenvalues, so SolveCubic is returning a theoretically incorrect result. The problem has to do with the precision of 'float' numbers. The cubic-solver approach can be made more robust by preconditioning the matrix.
Root finding using the closed-form algebraic equations for low-degree polynomials is known to be numerically ill conditioned. It is better to use iterative methods. For 3x3 symmetric matrices, a specialized iterative method can be very fast, extremely accurate, and robust.
int main()
{
// Cube with vertices at (0,0,0), (100,0,0), (0,100,0), (100,100,0),
// (0,0,100), (100,0,100), (0,100,100), and (100,100,100) rotated by
// an arbitrary rotation matrix.
Vector3<float> vertices[8] =
{
Vector3<float>{-59.372719f, 26.978844f, 39.425350f},
Vector3<float>{-29.475891f, 116.858276f, 7.364998f},
Vector3<float>{1.924920f, -16.858200f, -26.308300f},
Vector3<float>{13.762627f, 26.978811f, 107.625198f},
Vector3<float>{31.821747f, 73.021233f, -58.368652f,},
Vector3<float>{43.659454f, 116.858246f, 75.564850f},
Vector3<float>{75.060265f, -16.858232f, 41.891552f},
Vector3<float>{104.957092f, 73.021202f, 9.831200f}
};
DirectX::BoundingOrientedBox box;
DirectX::BoundingOrientedBox::CreateFromPoints(box, 8,
(const DirectX::XMFLOAT3*)&vertices[0], 3 * sizeof(float));
// box.Center = (22.7921867, 50.0000229, 24.6282730)
// box.Extents = (82.1649017, 66.8582535, 82.9969254)
// box.Orientation = (-0.000000000, 0.000000000, -0.000000000, 1.00000000)
//
// Internally, the parameters passed to SolveCubic (DirectXCollision.inl) are
// e = -59999.9922
// f = 1.19999974e+009
// g = -7.99999787e+012
// Local parameters are
// p = 128.000000
// q = 1048576.00
// h = 2.74877972e+011
// Because 'h > 0.0', the values t, u, and v are all set to zero and the
// function returns 'false'. The comment is "only one real root". The
// eigenvectors are then set to the default (1,0,0), (0,1,0), and (0,0,1).
// These are not correct should someone actually want the eigenvectors of
// the matrix (rather than some bounding box of the points).
//
// The SolveCubic is called to compute the eigenvalues of a real-valued symmetric
// 3x3 matrix. Such matrices ALWAYS have 3 real eigenvalues, so SolveCubic
// is returning a theoretically incorrect result. The problem has to do with
// the precision of 'float' numbers. The cubic-solver approach can be made
// more robust by preconditioning the matrix. For example, see
// www.geometrictools.com/.../EigenSymmetric3x3.pdf
//
// Root finding using the closed-form algebraic equations for low-degree
// polynomials is known to be numerically ill conditioned. It is better
// to use iterative methods. For 3x3 symmetric matrices, a specialized
// iterative method can be very fast, extremely accurate, and robust.
Vector3<float> mean{ 0.0f, 0.0f, 0.0f };
for (int i = 0; i < 8; ++i)
{
mean += vertices[i];
}
mean /= 8.0f;
float covar00 = 0.0f, covar01 = 0.0f, covar02 = 0.0f;
float covar11 = 0.0f, covar12 = 0.0f, covar22 = 0.0f;
for (int i = 0; i < 8; ++i)
{
Vector3<float> diff = vertices[i] - mean;
covar00 += diff[0] * diff[0];
covar01 += diff[0] * diff[1];
covar02 += diff[0] * diff[2];
covar11 += diff[1] * diff[1];
covar12 += diff[1] * diff[2];
covar22 += diff[2] * diff[2];
}
float e = -(covar00 + covar11 + covar22);
float f =
covar00 * covar11 +
covar11 * covar22 +
covar22 * covar00 -
covar01 * covar01 -
covar02 * covar02 -
covar12 * covar12;
float g =
covar01 * covar01 * covar22 +
covar02 * covar02 * covar11 +
covar12 * covar12 * covar00 -
covar01 * covar12 * covar02 * 2.0f -
covar00 * covar11 * covar22;
float p = f - e * e / 3.0f;
float q = g - e * f / 3.0f + e * e * e * 2.0f / 27.0f;
float h = q * q / 4.0f + p * p * p / 27.0f;
// To show how bad the floating-point errors are, compute
// e, f, g, p, q, and h using exact rational arithmetic.
typedef BSRational<UIntegerAP32> Rational;
Rational rcovar00(covar00);
Rational rcovar01(covar01);
Rational rcovar02(covar02);
Rational rcovar11(covar11);
Rational rcovar12(covar12);
Rational rcovar22(covar22);
Rational re = -(rcovar00 + rcovar11 + rcovar22);
Rational rf =
rcovar00 * rcovar11 +
rcovar11 * rcovar22 +
rcovar22 * rcovar00 -
rcovar01 * rcovar01 -
rcovar02 * rcovar02 -
rcovar12 * rcovar12;
Rational rg =
rcovar01 * rcovar01 * rcovar22 +
rcovar02 * rcovar02 * rcovar11 +
rcovar12 * rcovar12 * rcovar00 -
rcovar01 * rcovar12 * rcovar02 * Rational(2) -
rcovar00 * rcovar11 * rcovar22;
Rational rp = rf - re * re / Rational(3);
Rational rq = rg - re * rf / Rational(3) + re * re * re * Rational(2) / Rational(27);
Rational rh = rq * rq / Rational(4) + rp * rp * rp / Rational(27);
double drp, drq, drh;
drp = (double)rp; // -6.4158812165260315e-006
drq = (double)rq; // 1.9736035028472543e-009
drh = (double)rh; // -8.8077160212578135e-018
// Now let's compute the roots using the formulas in SolveCubic but with
// more accurate values for p, q, and h.
p = (float)drp;
q = (float)drq;
h = (float)drh;
float d = sqrtf(q * q / 4.0f - h);
float rc;
if (d < 0)
rc = -powf(-d, 1.0f / 3.0f);
else
rc = powf( d, 1.0f / 3.0f );
float theta = acos( -q / ( 2.0f * d ) );
float costh3 = cos( theta / 3.0f );
float sinth3 = sqrtf( 3.0f ) * sin( theta / 3.0f );
float t, u, v;
t = 2.0f * rc * costh3 - e / 3.0f; // 20000.00
u = -rc * ( costh3 + sinth3 ) - e / 3.0f; // 19999.9961
v = -rc * ( costh3 - sinth3 ) - e / 3.0f; // 19999.9980
// You can find the iterative solver at
// www.geometrictools.com/.../GteSymmetricEigensolver3x3.h
// based on the documentation
// www.geometrictools.com/.../RobustEigenSymmetric3x3.pdf
// None of this is intellectual property. The algorithms are based on the
// book "Matrix Computations" by Gene Golub and Charles van Loan, and their
// work is based on research papers from the 1960s. You are welcome to grab
// my code and use convert it to your DirectX math style. If Microsoft
// lawyers are uncomfortable with this, I'll be glad to provide a legal
// document that says you can use the code freely.
SymmetricEigensolver3x3<float> es;
std::array<float, 3> eval;
std::array<std::array<float, 3>, 3> evec;
es(covar00, covar01, covar02, covar11, covar12, covar22, false, +1,
eval, evec);
// eval = { 19999.9941, 19999.9961, 19999.9980 }
// evec[0] = { 0.0722742379, 0.997384787, 0.000000000 }
// evec[1] = { -0.913869500, 0.0662224069, -0.400571018 }
// evec[2] = { -0.399523437, 0.0289509650, 0.916265726 }
// Theoretically, the eigenvalues are all 20000 and the eigenvectors are
// the axes of the rotation matrix used to rotate the cube (see my initial
// statements). However, my solver is given the floating-point computations
// for the covariance matrix entries, so some rounding errors have already
// occurred when I execute my solver.
return 0;
}
int main()
{
// Cube with vertices at (0,0,0), (100,0,0), (0,100,0), (100,100,0),
// (0,0,100), (100,0,100), (0,100,100), and (100,100,100) rotated by
// an arbitrary rotation matrix.
Vector3
