Skip to content

Commit

Permalink
Merge pull request #2473 from pleroy/Diagonalize
Browse files Browse the repository at this point in the history
Deal with a degenerate case when diagonalizing a form
pleroy authored Feb 18, 2020
2 parents 585223b + a398933 commit b5dd83a
Showing 3 changed files with 38 additions and 3 deletions.
3 changes: 2 additions & 1 deletion geometry/symmetric_bilinear_form.hpp
Original file line number Diff line number Diff line change
@@ -66,7 +66,8 @@ class SymmetricBilinearForm {

// The eigensystem for a form is described by (1) the form in its eigenbasis,
// which gives the eigenvalues; and (2) a rotation from the current basis to
// the eigenbasis, which gives the eigenvectors.
// the eigenbasis, which gives the eigenvectors. The eigenvalues are in
// increasing order.
template<typename Eigenframe>
struct Eigensystem {
SymmetricBilinearForm<Scalar, Eigenframe, Multivector> form;
16 changes: 14 additions & 2 deletions geometry/symmetric_bilinear_form_body.hpp
Original file line number Diff line number Diff line change
@@ -129,6 +129,7 @@ template<typename Eigenframe>
typename SymmetricBilinearForm<Scalar, Frame, Multivector>::
template Eigensystem<Eigenframe>
SymmetricBilinearForm<Scalar, Frame, Multivector>::Diagonalize() const {
Scalar const zero;
R3x3Matrix<Scalar> const& A = matrix_;
auto const I = R3x3Matrix<double>::Identity();

@@ -138,20 +139,31 @@ typename SymmetricBilinearForm<Scalar, Frame, Multivector>::
Scalar const q = A.Trace() / 3;
R3x3Matrix<Scalar> const A_minus_qI = A - q * I;
Scalar const p = Sqrt((A_minus_qI * A_minus_qI).Trace() / 6);

if (p == zero) {
// A is very close to q * I.
SymmetricBilinearForm<Scalar, Eigenframe, Multivector> form(
R3x3Matrix<Scalar>({q, zero, zero},
{zero, q, zero},
{zero, zero, q}));
return {form, Rotation<Eigenframe, Frame>::Identity()};
}

R3x3Matrix<double> const B = A_minus_qI / p;
double const det_B = B.Determinant();
Angle const θ = ArcCos(det_B * 0.5);
Angle const θ = ArcCos(std::clamp(det_B, -2.0, 2.0) * 0.5);
double const β₀ = 2 * Cos(θ / 3);
double const β₁ = 2 * Cos((θ + 2 * π * Radian) / 3);
double const β₂ = 2 * Cos((θ + 4 * π * Radian) / 3);
std::array<Scalar, 3> αs = {p * β₀ + q, p * β₁ + q, p * β₂ + q};
// We expect αs[1] <= αs[2] <= αs[0] here, but sorting ensures that we are
// correct irrespective of numerical errors.
std::sort(αs.begin(), αs.end());
Scalar const& α₀ = αs[0];
Scalar const& α₁ = αs[1];
Scalar const& α₂ = αs[2];

// The form in its diagonal basis.
Scalar const zero;
SymmetricBilinearForm<Scalar, Eigenframe, Multivector> form(
R3x3Matrix<Scalar>({α₀, zero, zero},
{zero, α₁, zero},
22 changes: 22 additions & 0 deletions geometry/symmetric_bilinear_form_test.cpp
Original file line number Diff line number Diff line change
@@ -359,6 +359,28 @@ TEST_F(SymmetricBilinearFormTest, Diagonalize) {
VanishesBefore(1, 2),
AlmostEquals(1, 0)));
}

// A degenerate case.
{
auto const f = MakeSymmetricBilinearForm<World>(R3x3Matrix<double>(
{2, 0, 0},
{0, 2, 0},
{0, 0, 2}));
auto const f_eigensystem = f.Diagonalize<Eigenworld>();

Vector<double, Eigenworld> const e₀({1, 0, 0});
Vector<double, Eigenworld> const e₁({0, 1, 0});
Vector<double, Eigenworld> const e₂({0, 0, 1});
EXPECT_THAT(f_eigensystem.form(e₀, e₀), AlmostEquals(2 * Metre, 0));
EXPECT_THAT(f_eigensystem.form(e₀, e₁), AlmostEquals(0 * Metre, 0));
EXPECT_THAT(f_eigensystem.form(e₀, e₂), AlmostEquals(0 * Metre, 0));
EXPECT_THAT(f_eigensystem.form(e₁, e₀), AlmostEquals(0 * Metre, 0));
EXPECT_THAT(f_eigensystem.form(e₁, e₁), AlmostEquals(2 * Metre, 0));
EXPECT_THAT(f_eigensystem.form(e₁, e₂), AlmostEquals(0 * Metre, 0));
EXPECT_THAT(f_eigensystem.form(e₂, e₀), AlmostEquals(0 * Metre, 0));
EXPECT_THAT(f_eigensystem.form(e₂, e₁), AlmostEquals(0 * Metre, 0));
EXPECT_THAT(f_eigensystem.form(e₂, e₂), AlmostEquals(2 * Metre, 0));
}
}

} // namespace internal_symmetric_bilinear_form

0 comments on commit b5dd83a

Please sign in to comment.