// Copyright (C) 2022 Rémi Bèges - Jérôme "Lynix" Leclercq (lynix680@gmail.com) // This file is part of the "Nazara Engine - Math module" // For conditions of distribution and use, see copyright notice in Config.hpp #include #include #include #include #include #include #include #include #include #include namespace Nz { /*! * \ingroup math * \class Nz::Quaternion * \brief Math class that represents an element of the quaternions * * \remark The quaternion is meant to be 'unit' to represent rotations in a three dimensional space */ /*! * \brief Constructs a Quaternion object from its components * * \param W W component * \param X X component * \param Y Y component * \param Z Z component */ template Quaternion::Quaternion(T W, T X, T Y, T Z) { Set(W, X, Y, Z); } /*! * \brief Constructs a Quaternion object from an angle * * \param angle Angle representing a 2D rotation */ template template Quaternion::Quaternion(const Angle& angle) { Set(angle); } /*! * \brief Constructs a Quaternion object from a EulerAngles * * \param angles Easier representation of rotation of space * * \see EulerAngles */ template Quaternion::Quaternion(const EulerAngles& angles) { Set(angles); } /*! * \brief Constructs a Quaternion object from an angle and a direction * * \param angle Angle to rotate along the axis * \param axis Vector3 which represents a direction, no need to be normalized */ template Quaternion::Quaternion(RadianAngle angle, const Vector3& axis) { Set(angle, axis); } /*! * \brief Constructs a Quaternion object from an array of four elements * * \param quat[4] quat[0] is W component, quat[1] is X component, quat[2] is Y component and quat[3] is Z component */ template Quaternion::Quaternion(const T quat[4]) { Set(quat); } /* template Quaternion::Quaternion(const Matrix3& mat) { Set(mat); } */ /*! * \brief Constructs a Quaternion object from another type of Quaternion * * \param quat Quaternion of type U to convert to type T */ template template Quaternion::Quaternion(const Quaternion& quat) { Set(quat); } /*! * \brief Computes the w component of the quaternion to make it unit * \return A reference to this quaternion */ template Quaternion& Quaternion::ComputeW() { T t = T(1.0) - SquaredMagnitude(); if (t < T(0.0)) w = T(0.0); else w = -std::sqrt(t); return *this; } /*! * \brief Returns the rotational conjugate of this quaternion * \return A reference to this quaternion * * The conjugate of a quaternion represents the same rotation in the opposite direction about the rotational axis * * \see GetConjugate */ template Quaternion& Quaternion::Conjugate() { x = -x; y = -y; z = -z; return *this; } /*! * \brief Calculates the dot (scalar) product with two quaternions * \return The value of the dot product * * \param quat The other quaternion to calculate the dot product with */ template T Quaternion::DotProduct(const Quaternion& quat) const { return w * quat.w + x * quat.x + y * quat.y + z * quat.z; } /*! * \brief Gets the rotational conjugate of this quaternion * \return A new quaternion which is the conjugate of this quaternion * * The conjugate of a quaternion represents the same rotation in the opposite direction about the rotational axis * * \see Conjugate */ template Quaternion Quaternion::GetConjugate() const { Quaternion quat(*this); quat.Conjugate(); return quat; } /*! * \brief Gets the inverse of this quaternion * \return A new quaternion which is the inverse of this quaternion * * \remark If this quaternion is (0, 0, 0, 0), then it returns (0, 0, 0, 0) * * \see Inverse */ template Quaternion Quaternion::GetInverse() const { Quaternion quat(*this); quat.Inverse(); return quat; } /*! * \brief Gets the normalization of this quaternion * \return A new quaternion which is the normalization of this quaternion * * \param length Optional argument to obtain the length's ratio of the quaternion and the unit-length * * \remark If this quaternion is (0, 0, 0, 0), then it returns (0, 0, 0, 0) and length is 0 * * \see Normalize */ template Quaternion Quaternion::GetNormal(T* length) const { Quaternion quat(*this); quat.Normalize(length); return quat; } /*! * \brief Inverts this quaternion * \return A reference to this quaternion which is now inverted * * \remark If this quaternion is (0, 0, 0, 0), then it returns (0, 0, 0, 0) * * \see GetInverse */ template Quaternion& Quaternion::Inverse() { T norm = SquaredMagnitude(); if (norm > T(0.0)) { T invNorm = T(1.0) / std::sqrt(norm); w *= invNorm; x *= -invNorm; y *= -invNorm; z *= -invNorm; } return *this; } /*! * \brief Makes the quaternion (1, 0, 0, 0) * \return A reference to this vector with components (1, 0, 0, 0) * * \see Unit */ template Quaternion& Quaternion::MakeIdentity() { return Set(T(1.0), T(0.0), T(0.0), T(0.0)); } /*! * \brief Makes this quaternion to the rotation required to rotate direction Vector3 from to direction Vector3 to * \return A reference to this vector which is the rotation needed * * \param from Initial vector * \param to Target vector * * \remark Vectors are not required to be normalized * * \see RotationBetween */ template Quaternion& Quaternion::MakeRotationBetween(const Vector3& from, const Vector3& to) { // Based on: http://lolengine.net/blog/2013/09/18/beautiful-maths-quaternion-from-vectors T norm = std::sqrt(from.GetSquaredLength() * to.GetSquaredLength()); Vector3 crossProduct = from.CrossProduct(to); Set(norm + from.DotProduct(to), crossProduct.x, crossProduct.y, crossProduct.z); return Normalize(); } /*! * \brief Makes the quaternion (0, 0, 0, 0) * \return A reference to this vector with components (0, 0, 0, 0) * * \see Zero */ template Quaternion& Quaternion::MakeZero() { return Set(T(0.0), T(0.0), T(0.0), T(0.0)); } /*! * \brief Calculates the magnitude (length) of the quaternion * \return The magnitude * * \see SquaredMagnitude */ template T Quaternion::Magnitude() const { return std::sqrt(SquaredMagnitude()); } /*! * \brief Normalizes the current quaternion * \return A reference to this quaternion which is now normalized * * \param length Optional argument to obtain the length's ratio of the quaternion and the unit-length * * \remark If the quaternion is (0, 0, 0, 0), then it returns (0, 0, 0, 0) and length is 0 * * \see GetNormal */ template Quaternion& Quaternion::Normalize(T* length) { T norm = std::sqrt(SquaredMagnitude()); if (norm > T(0.0)) { T invNorm = T(1.0) / norm; w *= invNorm; x *= invNorm; y *= invNorm; z *= invNorm; } if (length) *length = norm; return *this; } /*! * \brief Sets the components of the quaternion * \return A reference to this quaternion * * \param W W component * \param X X component * \param Y Y component * \param Z Z component */ template Quaternion& Quaternion::Set(T W, T X, T Y, T Z) { w = W; x = X; y = Y; z = Z; return *this; } /*! * \brief Sets this quaternion from a 2D rotation specified by an Angle * \return A reference to this quaternion * * \param angle 2D angle * * \see Angle */ template template Quaternion& Quaternion::Set(const Angle& angle) { return Set(angle.ToQuaternion()); } /*! * \brief Sets this quaternion from rotation specified by Euler angle * \return A reference to this quaternion * * \param angles Easier representation of rotation of space * * \see EulerAngles */ template Quaternion& Quaternion::Set(const EulerAngles& angles) { return Set(angles.ToQuaternion()); } /*! * \brief Sets this quaternion from rotation specified by axis and angle * \return A reference to this quaternion * * \param angle Angle to rotate along the axis * \param axis Vector3 which represents a direction, no need to be normalized */ template Quaternion& Quaternion::Set(RadianAngle angle, const Vector3& axis) { angle /= T(2.0); Vector3 normalizedAxis = axis.GetNormal(); auto sincos = angle.GetSinCos(); w = sincos.second; x = normalizedAxis.x * sincos.first; y = normalizedAxis.y * sincos.first; z = normalizedAxis.z * sincos.first; return Normalize(); } /*! * \brief Sets the components of the quaternion from an array of four elements * \return A reference to this quaternion * * \param quat[4] quat[0] is W component, quat[1] is X component, quat[2] is Y component and quat[3] is Z component */ template Quaternion& Quaternion::Set(const T quat[4]) { w = quat[0]; x = quat[1]; y = quat[2]; z = quat[3]; return *this; } /*! * \brief Sets the components of the quaternion from another type of Quaternion * \return A reference to this quaternion * * \param quat Quaternion of type U to convert its components */ template template Quaternion& Quaternion::Set(const Quaternion& quat) { w = T(quat.w); x = T(quat.x); y = T(quat.y); z = T(quat.z); return *this; } /*! * \brief Calculates the squared magnitude (length) of the quaternion * \return The squared magnitude * * \see Magnitude */ template T Quaternion::SquaredMagnitude() const { return w * w + x * x + y * y + z * z; } /*! * \brief Returns the "roll angle" of this quaternion * \return Roll rotation * * \remark This function only has sense when quaternion only represents a "roll rotation" */ template RadianAngle Quaternion::To2DAngle() const { T siny_cosp = T(2.0) * (w * z + x * y); T cosy_cosp = T(1.0) - T(2.0) * (y * y + z * z); return std::atan2(siny_cosp, cosy_cosp); } /*! * \brief Converts this quaternion to Euler angles representation * \return EulerAngles which is the representation of this rotation * * \remark Rotation are "left-handed" */ template EulerAngles Quaternion::ToEulerAngles() const { T test = x * y + z * w; if (test > T(0.499)) // singularity at north pole return EulerAngles(DegreeAngle(T(0.0)), RadianAngle(T(2.0) * std::atan2(x, w)), DegreeAngle(T(90.0))); if (test < T(-0.499)) // singularity at south pole return EulerAngles(DegreeAngle(T(0.0)), RadianAngle(T(-2.0) * std::atan2(x, w)), DegreeAngle(T(-90.0))); return EulerAngles(RadianAngle(std::atan2(T(2.0) * x * w - T(2.0) * y * z, T(1.0) - T(2.0) * x * x - T(2.0) * z * z)), RadianAngle(std::atan2(T(2.0) * y * w - T(2.0) * x * z, T(1.0) - T(2.0) * y * y - T(2.0) * z * z)), RadianAngle(std::asin(T(2.0) * test))); } /*! * \brief Gives a string representation * \return A string representation of the object: "Quaternion(w | x, y, z)" */ template std::string Quaternion::ToString() const { std::ostringstream ss; ss << *this; return ss.str(); } /*! * \brief Adds the components of the quaternion with other quaternion * \return A quaternion where components are the sum of this quaternion and the other one * * \param quat The other quaternion to add components with */ template Quaternion Quaternion::operator+(const Quaternion& quat) const { Quaternion result; result.w = w + quat.w; result.x = x + quat.x; result.y = y + quat.y; result.z = z + quat.z; return result; } /*! * \brief Multiplies of the quaternion with other quaternion * \return A quaternion which is the product of those two according to operator* in quaternions * * \param quat The other quaternion to multiply with */ template Quaternion Quaternion::operator*(const Quaternion& quat) const { Quaternion result; result.w = w * quat.w - x * quat.x - y * quat.y - z * quat.z; result.x = w * quat.x + x * quat.w + y * quat.z - z * quat.y; result.y = w * quat.y + y * quat.w + z * quat.x - x * quat.z; result.z = w * quat.z + z * quat.w + x * quat.y - y * quat.x; return result; } /*! * \brief Apply the quaternion to the Vector3 * \return A Vector3f which is the vector rotated by this quaternion * * \param vec The vector to multiply with */ template Vector3 Quaternion::operator*(const Vector3& vec) const { Vector3 quatVec(x, y, z); Vector3 uv = quatVec.CrossProduct(vec); Vector3 uuv = quatVec.CrossProduct(uv); uv *= T(2.0) * w; uuv *= T(2.0); return vec + uv + uuv; } /*! * \brief Multiplies the components of the quaternion with a scalar * \return A quaternion where components are the product of this quaternion and the scalar * * \param scale The scalar to multiply components with */ template Quaternion Quaternion::operator*(T scale) const { return Quaternion(w * scale, x * scale, y * scale, z * scale); } /*! * \brief Divides the quaternion with other quaternion * \return A quaternion which is the quotient of those two according to operator* in quaternions * * \param quat The other quaternion to divide with */ template Quaternion Quaternion::operator/(const Quaternion& quat) const { return quat.GetConjugate() * (*this); } /*! * \brief Adds the components of the quaternion with other quaternion * \return A reference to this quaternion where components are the sum of this quaternion and the other one * * \param quat The other quaternion to add components with */ template Quaternion& Quaternion::operator+=(const Quaternion& quat) { return operator=(operator+(quat)); } /*! * \brief Multiplies of the quaternion with other quaternion * \return A reference to this quaternion which is the product of those two according to operator* in quaternions * * \param quat The other quaternion to multiply with */ template Quaternion& Quaternion::operator*=(const Quaternion& quat) { return operator=(operator*(quat)); } /*! * \brief Multiplies the components of the quaternion with a scalar * \return A reference to this quaternion where components are the product of this quaternion and the scalar * * \param scale The scalar to multiply components with */ template Quaternion& Quaternion::operator*=(T scale) { return operator=(operator*(scale)); } /*! * \brief Divides the quaternion with other quaternion * \return A reference to this quaternion which is the quotient of those two according to operator* in quaternions * * \param quat The other quaternion to divide with */ template Quaternion& Quaternion::operator/=(const Quaternion& quat) { return operator=(operator/(quat)); } /*! * \brief Compares the quaternion to other one * \return true if the quaternions are the same * * \param quat Other quaternion to compare with */ template bool Quaternion::operator==(const Quaternion& quat) const { return NumberEquals(w, quat.w) && NumberEquals(x, quat.x) && NumberEquals(y, quat.y) && NumberEquals(z, quat.z); } /*! * \brief Compares the quaternion to other one * \return false if the quaternions are the same * * \param quat Other quaternion to compare with */ template bool Quaternion::operator!=(const Quaternion& quat) const { return !operator==(quat); } /*! * \brief Shorthand for the quaternion (1, 0, 0, 0) * \return A quaternion with components (1, 0, 0, 0) * * \see MakeIdentity */ template Quaternion Quaternion::Identity() { Quaternion quaternion; quaternion.MakeIdentity(); return quaternion; } /*! * \brief Interpolates the quaternion to other one with a factor of interpolation * \return A new quaternion which is the interpolation of two quaternions * * \param from Initial quaternion * \param to Target quaternion * \param interpolation Factor of interpolation * * \remark interpolation is meant to be between 0 and 1, other values are potentially undefined behavior * \remark With NAZARA_DEBUG, a NazaraError is thrown and Zero() is returned * * \see Lerp, Slerp */ template Quaternion Quaternion::Lerp(const Quaternion& from, const Quaternion& to, T interpolation) { #ifdef NAZARA_DEBUG if (interpolation < T(0.0) || interpolation > T(1.0)) { NazaraError("Interpolation must be in range [0..1] (Got " + NumberToString(interpolation) + ')'); return Zero(); } #endif Quaternion interpolated; interpolated.w = Nz::Lerp(from.w, to.w, interpolation); interpolated.x = Nz::Lerp(from.x, to.x, interpolation); interpolated.y = Nz::Lerp(from.y, to.y, interpolation); interpolated.z = Nz::Lerp(from.z, to.z, interpolation); return interpolated; } template Quaternion Quaternion::LookAt(const Vector3& forward, const Vector3& up) { // From https://gamedev.stackexchange.com/questions/53129/quaternion-look-at-with-up-vector Vector3 forward_w(1, 0, 0); Vector3 axis = Vector3::CrossProduct(forward, forward_w); RadianAngle angle = std::acos(Vector3::DotProduct(forward, forward_w)); Vector3 third = Vector3::CrossProduct(axis, forward_w); if (Vector3::DotProduct(third, forward) < 0) angle = -angle; Quaternion q1 = Quaternion(angle, axis); Vector3 up_l = q1 * up; Vector3 right = Vector3::Normalize(Vector3::CrossProduct(forward, up)); Vector3 up_w = Vector3::Normalize(Vector3::CrossProduct(right, forward)); Vector3 axis2 = Vector3::CrossProduct(up_l, up_w); RadianAngle angle2 = std::acos(Vector3::DotProduct(forward, forward_w)); Quaternion q2 = Quaternion(angle2, axis2); return q2 * q1; } /*! * \brief Gives the normalized quaternion * \return A normalized quaternion from the quat * * \param quat Quaternion to normalize * \param length Optional argument to obtain the length's ratio of the vector and the unit-length * * \see GetNormal */ template Quaternion Quaternion::Normalize(const Quaternion& quat, T* length) { return quat.GetNormal(length); } /*! * \brief Gets the rotation required to rotate direction Vector3 from to direction Vector3 to * \return A quaternion which is the rotation needed between those two Vector3 * * \param from Initial vector * \param to Target vector * * \see MakeRotationBetween */ template Quaternion Quaternion::RotationBetween(const Vector3& from, const Vector3& to) { Quaternion quaternion; quaternion.MakeRotationBetween(from, to); return quaternion; } template Quaternion Quaternion::Mirror(Quaternion quat, const Vector3& axis) { float x = std::copysign(T(1.0), axis.x); float y = std::copysign(T(1.0), axis.y); float z = std::copysign(T(1.0), axis.z); quat.x = y * z * quat.x; quat.y = x * z * quat.y; quat.z = x * y * quat.z; return quat; } /*! * \brief Interpolates spherically the quaternion to other one with a factor of interpolation * \return A new quaternion which is the interpolation of two quaternions * * \param from Initial quaternion * \param to Target quaternion * \param interpolation Factor of interpolation * * \remark interpolation is meant to be between 0 and 1, other values are potentially undefined behavior * \remark With NAZARA_DEBUG, a NazaraError is thrown and Zero() is returned * * \see Lerp */ template Quaternion Quaternion::Slerp(const Quaternion& from, const Quaternion& to, T interpolation) { #ifdef NAZARA_DEBUG if (interpolation < T(0.0) || interpolation > T(1.0)) { NazaraError("Interpolation must be in range [0..1] (Got " + std::to_string(interpolation) + ')'); return Zero(); } #endif Quaternion q; T cosOmega = from.DotProduct(to); if (cosOmega < T(0.0)) { // We invert everything q.Set(-to.w, -to.x, -to.y, -to.z); cosOmega = -cosOmega; } else q.Set(to); T k0, k1; if (cosOmega > T(0.9999)) { // Linear interpolation to avoid division by zero k0 = T(1.0) - interpolation; k1 = interpolation; } else { T sinOmega = std::sqrt(T(1.0) - cosOmega*cosOmega); T omega = std::atan2(sinOmega, cosOmega); // To avoid two divisions sinOmega = T(1.0)/sinOmega; k0 = std::sin((T(1.0) - interpolation) * omega) * sinOmega; k1 = std::sin(interpolation*omega) * sinOmega; } Quaternion result(k0 * from.w, k0 * from.x, k0 * from.y, k0 * from.z); return result += q * k1; } /*! * \brief Shorthand for the quaternion (0, 0, 0, 0) * \return A quaternion with components (0, 0, 0, 0) * * \see MakeZero */ template Quaternion Quaternion::Zero() { Quaternion quaternion; quaternion.MakeZero(); return quaternion; } /*! * \brief Serializes a Quaternion * \return true if successfully serialized * * \param context Serialization context * \param quat Input Quaternion */ template bool Serialize(SerializationContext& context, const Quaternion& quat, TypeTag>) { if (!Serialize(context, quat.x)) return false; if (!Serialize(context, quat.y)) return false; if (!Serialize(context, quat.z)) return false; if (!Serialize(context, quat.w)) return false; return true; } /*! * \brief Unserializes a Quaternion * \return true if successfully unserialized * * \param context Serialization context * \param quat Output Quaternion */ template bool Unserialize(SerializationContext& context, Quaternion* quat, TypeTag>) { if (!Unserialize(context, &quat->x)) return false; if (!Unserialize(context, &quat->y)) return false; if (!Unserialize(context, &quat->z)) return false; if (!Unserialize(context, &quat->w)) return false; return true; } } /*! * \brief Output operator * \return The stream * * \param out The stream * \param quat The quaternion to output */ template std::ostream& operator<<(std::ostream& out, const Nz::Quaternion& quat) { return out << "Quaternion(" << quat.w << " | " << quat.x << ", " << quat.y << ", " << quat.z << ')'; } #include