BoundingOrientedBox::CreateFromPoints solver suffers from numeric problems

[walbourn]

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.

[walbourn]

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;
}

[anishsingh935]

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 vertices[8] =
{
Vector3{-59.372719f, 26.978844f, 39.425350f},
Vector3{-29.475891f, 116.858276f, 7.364998f},
Vector3{1.924920f, -16.858200f, -26.308300f},
Vector3{13.762627f, 26.978811f, 107.625198f},
Vector3{31.821747f, 73.021233f, -58.368652f,},
Vector3{43.659454f, 116.858246f, 75.564850f},
Vector3{75.060265f, -16.858232f, 41.891552f},
Vector3{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 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 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 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 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;
}