Math: fix angle() returning NaN for close arguments.

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
 

src/Magnum/Math/Complex.h

@@ -67,14 +67,18 @@ template<class T> inline T dot(const Complex<T>& a, const Complex<T>& 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<T>&, const Quaternion<T>&),
     @ref angle(const Vector<size, FloatingPoint>&, const Vector<size, FloatingPoint>&)
 */
 template<class T> inline Rad<T> angle(const Complex<T>& normalizedA, const Complex<T>& normalizedB) {
     CORRADE_ASSERT(normalizedA.isNormalized() && normalizedB.isNormalized(),
         "Math::angle(): complex numbers" << normalizedA << "and" << normalizedB << "are not normalized", {});
-    return Rad<T>(std::acos(dot(normalizedA, normalizedB)));
+    return Rad<T>(std::acos(clamp(dot(normalizedA, normalizedB), T(-1), T(1))));
 }
 
 /**

src/Magnum/Math/Quaternion.h

@@ -63,14 +63,18 @@ template<class T> inline T dot(const Quaternion<T>& a, const Quaternion<T>& 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<T>&, const Complex<T>&),
     @ref angle(const Vector<size, FloatingPoint>&, const Vector<size, FloatingPoint>&)
  */
 template<class T> inline Rad<T> angle(const Quaternion<T>& normalizedA, const Quaternion<T>& normalizedB) {
     CORRADE_ASSERT(normalizedA.isNormalized() && normalizedB.isNormalized(),
         "Math::angle(): quaternions" << normalizedA << "and" << normalizedB << "are not normalized", {});
-    return Rad<T>{std::acos(dot(normalizedA, normalizedB))};
+    return Rad<T>{std::acos(clamp(dot(normalizedA, normalizedB), T(-1), T(1)))};
 }
 
 /** @relatesalso Quaternion

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

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

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

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<T>&, const Complex<T>&),
     @ref angle(const Quaternion<T>&, const Quaternion<T>&)
@@ -125,7 +129,7 @@ typename std::enable_if<std::is_floating_point<FloatingPoint>::value, Rad<Floati
 angle(const Vector<size, FloatingPoint>& normalizedA, const Vector<size, FloatingPoint>& normalizedB) {
     CORRADE_ASSERT(normalizedA.isNormalized() && normalizedB.isNormalized(),
         "Math::angle(): vectors" << normalizedA << "and" << normalizedB << "are not normalized", {});
-    return Rad<FloatingPoint>(std::acos(dot(normalizedA, normalizedB)));
+    return Rad<FloatingPoint>(std::acos(clamp(dot(normalizedA, normalizedB), FloatingPoint(-1), FloatingPoint(1))));
 }
 
 /**