diff --git a/glm/ext/quaternion_exponential.inl b/glm/ext/quaternion_exponential.inl index 1365d0a9..8456c00a 100644 --- a/glm/ext/quaternion_exponential.inl +++ b/glm/ext/quaternion_exponential.inl @@ -1,3 +1,5 @@ +#include "scalar_constants.hpp" + namespace glm { template @@ -46,16 +48,30 @@ namespace glm //To deal with non-unit quaternions T magnitude = sqrt(x.x * x.x + x.y * x.y + x.z * x.z + x.w *x.w); - //Equivalent to raising a real number to a power - //Needed to prevent a division by 0 error later on - if(abs(x.w / magnitude) > static_cast(1) - epsilon() && abs(x.w / magnitude) < static_cast(1) + epsilon()) - return qua(pow(x.w, y), 0, 0, 0); + T Angle; + if(abs(x.w / magnitude) > cos_one_over_two()) + { + //Scalar component is close to 1; using it to recover angle would lose precision + //Instead, we use the non-scalar components since sin() is accurate around 0 + + //Prevent a division by 0 error later on + T VectorMagnitude = x.x * x.x + x.y * x.y + x.z * x.z; + if (glm::abs(VectorMagnitude - static_cast(0)) < glm::epsilon()) { + //Equivalent to raising a real number to a power + return qua(pow(x.w, y), 0, 0, 0); + } + + Angle = asin(sqrt(VectorMagnitude) / magnitude); + } + else + { + //Scalar component is small, shouldn't cause loss of precision + Angle = acos(x.w / magnitude); + } - T Angle = acos(x.w / magnitude); T NewAngle = Angle * y; T Div = sin(NewAngle) / sin(Angle); T Mag = pow(magnitude, y - static_cast(1)); - return qua(cos(NewAngle) * magnitude * Mag, x.x * Div * Mag, x.y * Div * Mag, x.z * Div * Mag); } diff --git a/glm/ext/quaternion_trigonometric.inl b/glm/ext/quaternion_trigonometric.inl index 0c0bf63f..06b7c4c3 100644 --- a/glm/ext/quaternion_trigonometric.inl +++ b/glm/ext/quaternion_trigonometric.inl @@ -1,8 +1,15 @@ +#include "scalar_constants.hpp" + namespace glm { template GLM_FUNC_QUALIFIER T angle(qua const& x) { + if (abs(x.w) > cos_one_over_two()) + { + return asin(sqrt(x.x * x.x + x.y * x.y + x.z * x.z)) * static_cast(2); + } + return acos(x.w) * static_cast(2); } diff --git a/glm/ext/scalar_constants.hpp b/glm/ext/scalar_constants.hpp index 4a61faa5..74e210d9 100644 --- a/glm/ext/scalar_constants.hpp +++ b/glm/ext/scalar_constants.hpp @@ -30,6 +30,10 @@ namespace glm template GLM_FUNC_DECL GLM_CONSTEXPR genType pi(); + /// Return the value of cos(1 / 2) for floating point types. + template + GLM_FUNC_DECL GLM_CONSTEXPR genType cos_one_over_two(); + /// @} } //namespace glm diff --git a/glm/ext/scalar_constants.inl b/glm/ext/scalar_constants.inl index f97c9802..b475adf8 100644 --- a/glm/ext/scalar_constants.inl +++ b/glm/ext/scalar_constants.inl @@ -15,4 +15,10 @@ namespace glm GLM_STATIC_ASSERT(std::numeric_limits::is_iec559, "'pi' only accepts floating-point inputs"); return static_cast(3.14159265358979323846264338327950288); } + + template + GLM_FUNC_QUALIFIER GLM_CONSTEXPR genType cos_one_over_two() + { + return genType(0.877582561890372716130286068203503191); + } } //namespace glm diff --git a/glm/gtc/constants.inl b/glm/gtc/constants.inl index 87f5c860..bb98c6bf 100644 --- a/glm/gtc/constants.inl +++ b/glm/gtc/constants.inl @@ -163,4 +163,5 @@ namespace glm { return genType(1.61803398874989484820458683436563811); } + } //namespace glm