From 498f4c59e1186a46cd8b30547d9383a68498a8a5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Vladim=C3=ADr=20Vondru=C5=A1?= Date: Tue, 7 Apr 2020 19:52:05 +0200 Subject: [PATCH] Math: fix angle() returning NaN for close arguments. --- doc/changelog.dox | 2 ++ src/Magnum/Math/Complex.h | 6 +++++- src/Magnum/Math/Quaternion.h | 6 +++++- src/Magnum/Math/Test/ComplexTest.cpp | 13 +++++++++++++ src/Magnum/Math/Test/QuaternionTest.cpp | 13 +++++++++++++ src/Magnum/Math/Test/VectorTest.cpp | 13 +++++++++++++ src/Magnum/Math/Vector.h | 6 +++++- 7 files changed, 56 insertions(+), 3 deletions(-) diff --git a/doc/changelog.dox b/doc/changelog.dox index 60743d8bb..99b95fc82 100644 --- a/doc/changelog.dox +++ b/doc/changelog.dox @@ -461,6 +461,8 @@ See also: - @ref Shaders::MeshVisualizer3D accidentally didn't enable @glsl noperspective @ce interpolation on desktop, resulting in minor wireframe rendering artifacts +- @ref Math::angle() got fixed to not produce NaN results for vectors, + complex numbers or quaternions very close to each other @subsection changelog-latest-deprecated Deprecated APIs diff --git a/src/Magnum/Math/Complex.h b/src/Magnum/Math/Complex.h index 7b2ae415c..b8b7b63a0 100644 --- a/src/Magnum/Math/Complex.h +++ b/src/Magnum/Math/Complex.h @@ -67,6 +67,10 @@ template inline T dot(const Complex& a, const Complex& b) { Expects that both complex numbers are normalized. @f[ \theta = \arccos \left( \frac{Re(c_0 \cdot c_1))}{|c_0| |c_1|} \right) = \arccos (a_0 a_1 + b_0 b_1) @f] + +To avoid numerical issues when two complex numbers are very close to each +other, the dot product is clamped to the @f$ [-1, +1] @f$ range before being +passed to @f$ \arccos @f$. @see @ref Complex::isNormalized(), @ref angle(const Quaternion&, const Quaternion&), @ref angle(const Vector&, const Vector&) @@ -74,7 +78,7 @@ Expects that both complex numbers are normalized. @f[ template inline Rad angle(const Complex& normalizedA, const Complex& normalizedB) { CORRADE_ASSERT(normalizedA.isNormalized() && normalizedB.isNormalized(), "Math::angle(): complex numbers" << normalizedA << "and" << normalizedB << "are not normalized", {}); - return Rad(std::acos(dot(normalizedA, normalizedB))); + return Rad(std::acos(clamp(dot(normalizedA, normalizedB), T(-1), T(1)))); } /** diff --git a/src/Magnum/Math/Quaternion.h b/src/Magnum/Math/Quaternion.h index 81ae3cc8a..be9fe94ac 100644 --- a/src/Magnum/Math/Quaternion.h +++ b/src/Magnum/Math/Quaternion.h @@ -63,6 +63,10 @@ template inline T dot(const Quaternion& a, const Quaternion& b) { Expects that both quaternions are normalized. @f[ \theta = \arccos \left( \frac{p \cdot q}{|p| |q|} \right) = \arccos(p \cdot q) @f] + +To avoid numerical issues when two complex numbers are very close to each +other, the dot product is clamped to the @f$ [-1, +1] @f$ range before being +passed to @f$ \arccos @f$. @see @ref Quaternion::isNormalized(), @ref angle(const Complex&, const Complex&), @ref angle(const Vector&, const Vector&) @@ -70,7 +74,7 @@ Expects that both quaternions are normalized. @f[ template inline Rad angle(const Quaternion& normalizedA, const Quaternion& normalizedB) { CORRADE_ASSERT(normalizedA.isNormalized() && normalizedB.isNormalized(), "Math::angle(): quaternions" << normalizedA << "and" << normalizedB << "are not normalized", {}); - return Rad{std::acos(dot(normalizedA, normalizedB))}; + return Rad{std::acos(clamp(dot(normalizedA, normalizedB), T(-1), T(1)))}; } /** @relatesalso Quaternion diff --git a/src/Magnum/Math/Test/ComplexTest.cpp b/src/Magnum/Math/Test/ComplexTest.cpp index 71780b035..3dc6e55b4 100644 --- a/src/Magnum/Math/Test/ComplexTest.cpp +++ b/src/Magnum/Math/Test/ComplexTest.cpp @@ -89,6 +89,7 @@ struct ComplexTest: Corrade::TestSuite::Tester { void invertedNormalizedNotNormalized(); void angle(); + void angleNormalizedButOver1(); void angleNotNormalized(); void rotation(); void matrix(); @@ -142,6 +143,7 @@ ComplexTest::ComplexTest() { &ComplexTest::invertedNormalizedNotNormalized, &ComplexTest::angle, + &ComplexTest::angleNormalizedButOver1, &ComplexTest::angleNotNormalized, &ComplexTest::rotation, &ComplexTest::matrix, @@ -442,6 +444,17 @@ void ComplexTest::angle() { CORRADE_COMPARE(Math::angle(a, -a), 180.0_degf); } +void ComplexTest::angleNormalizedButOver1() { + /* This complex *is* normalized, but its length is larger than 1, which + would cause acos() to return a NaN. Ensure it's clamped to correct range + before passing it there. */ + Complex a{1.0f + Math::TypeTraits::epsilon()/2, 0.0f}; + CORRADE_VERIFY(a.isNormalized()); + + CORRADE_COMPARE(Math::angle(a, a), 0.0_radf); + CORRADE_COMPARE(Math::angle(a, -a), 180.0_degf); +} + void ComplexTest::angleNotNormalized() { std::ostringstream out; Error redirectError{&out}; diff --git a/src/Magnum/Math/Test/QuaternionTest.cpp b/src/Magnum/Math/Test/QuaternionTest.cpp index 07579df75..70289280a 100644 --- a/src/Magnum/Math/Test/QuaternionTest.cpp +++ b/src/Magnum/Math/Test/QuaternionTest.cpp @@ -94,6 +94,7 @@ struct QuaternionTest: Corrade::TestSuite::Tester { void rotation(); void rotationNotNormalized(); void angle(); + void angleNormalizedButOver1(); void angleNotNormalized(); void matrix(); void matrixNotOrthogonal(); @@ -175,6 +176,7 @@ QuaternionTest::QuaternionTest() { &QuaternionTest::rotation, &QuaternionTest::rotationNotNormalized, &QuaternionTest::angle, + &QuaternionTest::angleNormalizedButOver1, &QuaternionTest::angleNotNormalized, &QuaternionTest::matrix, &QuaternionTest::matrixNotOrthogonal, @@ -511,6 +513,17 @@ void QuaternionTest::angle() { Corrade::TestSuite::Compare::around(0.0005_radf)); } +void QuaternionTest::angleNormalizedButOver1() { + /* This quaternion *is* normalized, but its length is larger than 1, which + would cause acos() to return a NaN. Ensure it's clamped to correct range + before passing it there. */ + Quaternion a{{1.0f + Math::TypeTraits::epsilon()/2, 0.0f, 0.0f}, 0.0f}; + CORRADE_VERIFY(a.isNormalized()); + + CORRADE_COMPARE(Math::angle(a, a), 0.0_radf); + CORRADE_COMPARE(Math::angle(a, -a), 180.0_degf); +} + void QuaternionTest::angleNotNormalized() { std::ostringstream out; Error redirectError{&out}; diff --git a/src/Magnum/Math/Test/VectorTest.cpp b/src/Magnum/Math/Test/VectorTest.cpp index 55ea310a4..7c1ffcce9 100644 --- a/src/Magnum/Math/Test/VectorTest.cpp +++ b/src/Magnum/Math/Test/VectorTest.cpp @@ -108,6 +108,7 @@ struct VectorTest: Corrade::TestSuite::Tester { void flipped(); void angle(); + void angleNormalizedButOver1(); void angleNotNormalized(); void subclassTypes(); @@ -181,6 +182,7 @@ VectorTest::VectorTest() { &VectorTest::flipped, &VectorTest::angle, + &VectorTest::angleNormalizedButOver1, &VectorTest::angleNotNormalized, &VectorTest::subclassTypes, @@ -583,6 +585,17 @@ void VectorTest::angle() { Corrade::TestSuite::Compare::around(0.0005_radf)); } +void VectorTest::angleNormalizedButOver1() { + /* This vector *is* normalized, but its length is larger than 1, which + would cause acos() to return a NaN. Ensure it's clamped to correct range + before passing it there. */ + Vector3 a{1.0f + Math::TypeTraits::epsilon()/2, 0.0f, 0.0f}; + CORRADE_VERIFY(a.isNormalized()); + + CORRADE_COMPARE(Math::angle(a, a), 0.0_radf); + CORRADE_COMPARE(Math::angle(a, -a), 180.0_degf); +} + void VectorTest::angleNotNormalized() { std::ostringstream out; Error redirectError{&out}; diff --git a/src/Magnum/Math/Vector.h b/src/Magnum/Math/Vector.h index 2ed8aa18d..83060ca9b 100644 --- a/src/Magnum/Math/Vector.h +++ b/src/Magnum/Math/Vector.h @@ -112,6 +112,10 @@ Expects that both vectors are normalized. Enabled only for floating-point types. @f[ \theta = \arccos \left( \frac{\boldsymbol a \cdot \boldsymbol b}{|\boldsymbol a| |\boldsymbol b|} \right) = \arccos (\boldsymbol a \cdot \boldsymbol b) @f] + +To avoid numerical issues when two vectors are very close to each other, the +dot product is clamped to the @f$ [-1, +1] @f$ range before being passed to +@f$ \arccos @f$. @see @ref Vector::isNormalized(), @ref angle(const Complex&, const Complex&), @ref angle(const Quaternion&, const Quaternion&) @@ -125,7 +129,7 @@ typename std::enable_if::value, Rad& normalizedA, const Vector& normalizedB) { CORRADE_ASSERT(normalizedA.isNormalized() && normalizedB.isNormalized(), "Math::angle(): vectors" << normalizedA << "and" << normalizedB << "are not normalized", {}); - return Rad(std::acos(dot(normalizedA, normalizedB))); + return Rad(std::acos(clamp(dot(normalizedA, normalizedB), FloatingPoint(-1), FloatingPoint(1)))); } /**