magnum/src/Math/Quaternion.h

#ifndef Magnum_Math_Quaternion_h
#define Magnum_Math_Quaternion_h
/*
    Copyright © 2010, 2011, 2012 Vladimír Vondruš <mosra@centrum.cz>

    This file is part of Magnum.

    Magnum is free software: you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License version 3
    only, as published by the Free Software Foundation.

    Magnum is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
    GNU Lesser General Public License version 3 for more details.
*/

/** @file
 * @brief Class Magnum::Math::Quaternion
 */

#include <cmath>
#include <Utility/Assert.h>
#include <Utility/Debug.h>

#include "Math/MathTypeTraits.h"
#include "Math/Matrix.h"
#include "Math/Vector3.h"

namespace Magnum { namespace Math {

/**
@brief %Quaternion
@tparam T   Underlying data type

Represents 3D rotation.
@see Magnum::Quaternion, DualQuaternion, Matrix4
*/
template<class T> class Quaternion {
    public:
        typedef T Type;                         /**< @brief Underlying data type */

        /**
         * @brief Dot product
         *
         * @f[
         *      p \cdot q = \boldsymbol p_V \cdot \boldsymbol q_V + p_S q_S
         * @f]
         * @see dot() const
         */
        inline static T dot(const Quaternion<T>& a, const Quaternion<T>& b) {
            /** @todo Use four-component SIMD implementation when available */
            return Vector3<T>::dot(a.vector(), b.vector()) + a.scalar()*b.scalar();
        }

        /**
         * @brief Angle between normalized quaternions
         *
         * Expects that both quaternions are normalized. @f[
         *      \theta = acos \left( \frac{p \cdot q}{|p| |q|} \right) = acos(p \cdot q)
         * @f]
         * @see Complex::angle(), Vector::angle()
         */
        inline static Rad<T> angle(const Quaternion<T>& normalizedA, const Quaternion<T>& normalizedB) {
            CORRADE_ASSERT(MathTypeTraits<T>::equals(normalizedA.dot(), T(1)) && MathTypeTraits<T>::equals(normalizedB.dot(), T(1)),
                           "Math::Quaternion::angle(): quaternions must be normalized", Rad<T>(std::numeric_limits<T>::quiet_NaN()));
            return Rad<T>(angleInternal(normalizedA, normalizedB));
        }

        /**
         * @brief Linear interpolation of two quaternions
         * @param normalizedA   First quaternion
         * @param normalizedB   Second quaternion
         * @param t             Interpolation phase (from range @f$ [0; 1] @f$)
         *
         * Expects that both quaternions are normalized. @f[
         *      q_{LERP} = \frac{(1 - t) q_A + t q_B}{|(1 - t) q_A + t q_B|}
         * @f]
         * @see slerp(), Math::lerp()
         */
        inline static Quaternion<T> lerp(const Quaternion<T>& normalizedA, const Quaternion<T>& normalizedB, T t) {
            CORRADE_ASSERT(MathTypeTraits<T>::equals(normalizedA.dot(), T(1)) && MathTypeTraits<T>::equals(normalizedB.dot(), T(1)),
                           "Math::Quaternion::lerp(): quaternions must be normalized",
                           Quaternion<T>({}, std::numeric_limits<T>::quiet_NaN()));
            return ((T(1) - t)*normalizedA + t*normalizedB).normalized();
        }

        /**
         * @brief Spherical linear interpolation of two quaternions
         * @param normalizedA   First quaternion
         * @param normalizedB   Second quaternion
         * @param t             Interpolation phase (from range @f$ [0; 1] @f$)
         *
         * Expects that both quaternions are normalized. @f[
         *      q_{SLERP} = \frac{sin((1 - t) \theta) q_A + sin(t \theta) q_B}{sin \theta}
         *      ~~~~~~~~~~
         *      \theta = acos \left( \frac{q_A \cdot q_B}{|q_A| \cdot |q_B|} \right) = acos(q_A \cdot q_B)
         * @f]
         * @see lerp()
         */
        inline static Quaternion<T> slerp(const Quaternion<T>& normalizedA, const Quaternion<T>& normalizedB, T t) {
            CORRADE_ASSERT(MathTypeTraits<T>::equals(normalizedA.dot(), T(1)) && MathTypeTraits<T>::equals(normalizedB.dot(), T(1)),
                           "Math::Quaternion::slerp(): quaternions must be normalized",
                           Quaternion<T>({}, std::numeric_limits<T>::quiet_NaN()));
            T a = angleInternal(normalizedA, normalizedB);
            return (std::sin((T(1) - t)*a)*normalizedA + std::sin(t*a)*normalizedB)/std::sin(a);
        }

        /**
         * @brief Rotation quaternion
         * @param angle             Rotation angle (counterclockwise)
         * @param normalizedAxis    Normalized rotation axis
         *
         * Expects that the rotation axis is normalized. @f[
         *      q = [\boldsymbol a \cdot sin \frac \theta 2, cos \frac \theta 2]
         * @f]
         * @see rotationAngle(), rotationAxis(), DualQuaternion::rotation(),
         *      Matrix4::rotation(), Complex::rotation(), Vector3::xAxis(),
         *      Vector3::yAxis(), Vector3::zAxis()
         */
        inline static Quaternion<T> rotation(Rad<T> angle, const Vector3<T>& normalizedAxis) {
            CORRADE_ASSERT(MathTypeTraits<T>::equals(normalizedAxis.dot(), T(1)),
                           "Math::Quaternion::rotation(): axis must be normalized", {});

            return {normalizedAxis*std::sin(T(angle)/2), std::cos(T(angle)/2)};
        }

        /**
         * @brief Default constructor
         *
         * Creates unit quaternion. @f[
         *      q = [\boldsymbol 0, 1]
         * @f]
         */
        inline constexpr /*implicit*/ Quaternion(): _scalar(T(1)) {}

        /**
         * @brief Construct quaternion from vector and scalar
         *
         * @f[
         *      q = [\boldsymbol v, s]
         * @f]
         */
        inline constexpr /*implicit*/ Quaternion(const Vector3<T>& vector, T scalar): _vector(vector), _scalar(scalar) {}

        /**
         * @brief Construct quaternion from vector
         *
         * To be used in transformations later. @f[
         *      q = [\boldsymbol v, 0]
         * @f]
         * @see transformVector(), transformVectorNormalized()
         */
        inline constexpr explicit Quaternion(const Vector3<T>& vector): _vector(vector), _scalar(T(0)) {}

        /** @brief Equality comparison */
        inline bool operator==(const Quaternion<T>& other) const {
            return _vector == other._vector && MathTypeTraits<T>::equals(_scalar, other._scalar);
        }

        /** @brief Non-equality comparison */
        inline bool operator!=(const Quaternion<T>& other) const {
            return !operator==(other);
        }

        /** @brief %Vector part */
        inline constexpr Vector3<T> vector() const { return _vector; }

        /** @brief %Scalar part */
        inline constexpr T scalar() const { return _scalar; }

        /**
         * @brief Rotation angle of unit quaternion
         *
         * Expects that the quaternion is normalized. @f[
         *      \theta = 2 \cdot acos q_S
         * @f]
         * @see rotationAxis(), rotation(), DualQuaternion::rotationAngle()
         */
        inline Rad<T> rotationAngle() const {
            CORRADE_ASSERT(MathTypeTraits<T>::equals(dot(), T(1)),
                           "Math::Quaternion::rotationAngle(): quaternion must be normalized",
                           Rad<T>(std::numeric_limits<T>::quiet_NaN()));
            return Rad<T>(T(2)*std::acos(_scalar));
        }

        /**
         * @brief Rotation axis of unit quaternion
         *
         * Expects that the quaternion is normalized. Returns either unit-length
         * vector for valid rotation quaternion or NaN vector for
         * default-constructed quaternion. @f[
         *      \boldsymbol a = \frac{\boldsymbol q_V}{\sqrt{1 - q_S^2}}
         * @f]
         * @see rotationAngle(), rotation()
         */
        inline Vector3<T> rotationAxis() const {
            CORRADE_ASSERT(MathTypeTraits<T>::equals(dot(), T(1)),
                           "Math::Quaternion::rotationAxis(): quaternion must be normalized",
                           {});
            return _vector/std::sqrt(1-pow2(_scalar));
        }

        /**
         * @brief Convert quaternion to rotation matrix
         *
         * @see DualQuaternion::matrix(), Matrix4::from(const Matrix<3, T>&, const Vector3<T>&)
         */
        Matrix<3, T> matrix() const {
            return {
                Vector<3, T>(T(1) - 2*pow2(_vector.y()) - 2*pow2(_vector.z()),
                    2*_vector.x()*_vector.y() + 2*_vector.z()*_scalar,
                        2*_vector.x()*_vector.z() - 2*_vector.y()*_scalar),
                Vector<3, T>(2*_vector.x()*_vector.y() - 2*_vector.z()*_scalar,
                    T(1) - 2*pow2(_vector.x()) - 2*pow2(_vector.z()),
                        2*_vector.y()*_vector.z() + 2*_vector.x()*_scalar),
                Vector<3, T>(2*_vector.x()*_vector.z() + 2*_vector.y()*_scalar,
                    2*_vector.y()*_vector.z() - 2*_vector.x()*_scalar,
                        T(1) - 2*pow2(_vector.x()) - 2*pow2(_vector.y()))
            };
        }

        /**
         * @brief Add and assign quaternion
         *
         * The computation is done in-place. @f[
         *      p + q = [\boldsymbol p_V + \boldsymbol q_V, p_S + q_S]
         * @f]
         */
        inline Quaternion<T>& operator+=(const Quaternion<T>& other) {
            _vector += other._vector;
            _scalar += other._scalar;
            return *this;
        }

        /**
         * @brief Add quaternion
         *
         * @see operator+=()
         */
        inline Quaternion<T> operator+(const Quaternion<T>& other) const {
            return Quaternion<T>(*this) += other;
        }

        /**
         * @brief Negated quaternion
         *
         * @f[
         *      -q = [-\boldsymbol q_V, -q_S]
         * @f]
         */
        inline Quaternion<T> operator-() const {
            return {-_vector, -_scalar};
        }

        /**
         * @brief Subtract and assign quaternion
         *
         * The computation is done in-place. @f[
         *      p - q = [\boldsymbol p_V - \boldsymbol q_V, p_S - q_S]
         * @f]
         */
        inline Quaternion<T>& operator-=(const Quaternion<T>& other) {
            _vector -= other._vector;
            _scalar -= other._scalar;
            return *this;
        }

        /**
         * @brief Subtract quaternion
         *
         * @see operator-=()
         */
        inline Quaternion<T> operator-(const Quaternion<T>& other) const {
            return Quaternion<T>(*this) -= other;
        }

        /**
         * @brief Multiply with scalar and assign
         *
         * The computation is done in-place. @f[
         *      q \cdot a = [\boldsymbol q_V \cdot a, q_S \cdot a]
         * @f]
         */
        inline Quaternion<T>& operator*=(T scalar) {
            _vector *= scalar;
            _scalar *= scalar;
            return *this;
        }

        /**
         * @brief Multiply with scalar
         *
         * @see operator*=(T)
         */
        inline Quaternion<T> operator*(T scalar) const {
            return Quaternion<T>(*this) *= scalar;
        }

        /**
         * @brief Divide with scalar and assign
         *
         * The computation is done in-place. @f[
         *      \frac q a = [\frac {\boldsymbol q_V} a, \frac {q_S} a]
         * @f]
         */
        inline Quaternion<T>& operator/=(T scalar) {
            _vector /= scalar;
            _scalar /= scalar;
            return *this;
        }

        /**
         * @brief Divide with scalar
         *
         * @see operator/=(T)
         */
        inline Quaternion<T> operator/(T scalar) const {
            return Quaternion<T>(*this) /= scalar;
        }

        /**
         * @brief Multiply with quaternion
         *
         * @f[
         *      p q = [p_S \boldsymbol q_V + q_S \boldsymbol p_V + \boldsymbol p_V \times \boldsymbol q_V,
         *             p_S q_S - \boldsymbol p_V \cdot \boldsymbol q_V]
         * @f]
         */
        inline Quaternion<T> operator*(const Quaternion<T>& other) const {
            return {_scalar*other._vector + other._scalar*_vector + Vector3<T>::cross(_vector, other._vector),
                    _scalar*other._scalar - Vector3<T>::dot(_vector, other._vector)};
        }

        /**
         * @brief Dot product of the quaternion
         *
         * Should be used instead of length() for comparing quaternion length
         * with other values, because it doesn't compute the square root. @f[
         *      q \cdot q = \boldsymbol q_V \cdot \boldsymbol q_V + q_S^2
         * @f]
         * @see dot(const Quaternion&, const Quaternion&)
         */
        inline T dot() const {
            return dot(*this, *this);
        }

        /**
         * @brief %Quaternion length
         *
         * See also dot() const which is faster for comparing length with other
         * values. @f[
         *      |q| = \sqrt{q \cdot q}
         * @f]
         */
        inline T length() const {
            return std::sqrt(dot());
        }

        /** @brief Normalized quaternion (of unit length) */
        inline Quaternion<T> normalized() const {
            return (*this)/length();
        }

        /**
         * @brief Conjugated quaternion
         *
         * @f[
         *      q^* = [-\boldsymbol q_V, q_S]
         * @f]
         */
        inline Quaternion<T> conjugated() const {
            return {-_vector, _scalar};
        }

        /**
         * @brief Inverted quaternion
         *
         * See invertedNormalized() which is faster for normalized
         * quaternions. @f[
         *      q^{-1} = \frac{q^*}{|q|^2} = \frac{q^*}{q \cdot q}
         * @f]
         */
        inline Quaternion<T> inverted() const {
            return conjugated()/dot();
        }

        /**
         * @brief Inverted normalized quaternion
         *
         * Equivalent to conjugated(). Expects that the quaternion is
         * normalized. @f[
         *      q^{-1} = \frac{q^*}{|q|^2} = q^*
         * @f]
         * @see inverted()
         */
        inline Quaternion<T> invertedNormalized() const {
            CORRADE_ASSERT(MathTypeTraits<T>::equals(dot(), T(1)),
                           "Math::Quaternion::invertedNormalized(): quaternion must be normalized",
                           Quaternion<T>({}, std::numeric_limits<T>::quiet_NaN()));
            return conjugated();
        }

        /**
         * @brief Rotate vector with quaternion
         *
         * See transformVectorNormalized(), which is faster for normalized
         * quaternions. @f[
         *      v' = qvq^{-1} = q [\boldsymbol v, 0] q^{-1}
         * @f]
         * @see Quaternion(const Vector3&), vector(), Matrix4::transformVector(),
         *      DualQuaternion::transformPoint(), Complex::transformVector()
         */
        inline Vector3<T> transformVector(const Vector3<T>& vector) const {
            return ((*this)*Quaternion<T>(vector)*inverted()).vector();
        }

        /**
         * @brief Rotate vector with normalized quaternion
         *
         * Faster alternative to transformVector(), expects that the quaternion
         * is normalized. @f[
         *      v' = qvq^{-1} = qvq^* = q [\boldsymbol v, 0] q^*
         * @f]
         * @see Quaternion(const Vector3&), vector(), Matrix4::transformVector(),
         *      DualQuaternion::transformPointNormalized(), Complex::transformVector()
         */
        inline Vector3<T> transformVectorNormalized(const Vector3<T>& vector) const {
            CORRADE_ASSERT(MathTypeTraits<T>::equals(dot(), T(1)),
                           "Math::Quaternion::transformVectorNormalized(): quaternion must be normalized",
                           Vector3<T>(std::numeric_limits<T>::quiet_NaN()));
            return ((*this)*Quaternion<T>(vector)*conjugated()).vector();
        }

    private:
        /* Used to avoid including Functions.h */
        inline constexpr static T pow2(T value) {
            return value*value;
        }

        /* Used in angle() and slerp() (no assertions) */
        inline static T angleInternal(const Quaternion<T>& normalizedA, const Quaternion<T>& normalizedB) {
            return std::acos(dot(normalizedA, normalizedB));
        }

        Vector3<T> _vector;
        T _scalar;
};

/** @relates Quaternion
@brief Multiply scalar with quaternion

Same as Quaternion::operator*(T) const.
*/
template<class T> inline Quaternion<T> operator*(T scalar, const Quaternion<T>& quaternion) {
    return quaternion*scalar;
}

/** @relates Quaternion
@brief Divide quaternion with number and invert

@f[
    \frac a q = [\frac a {\boldsymbol q_V}, \frac a {q_S}]
@f]
@see Quaternion::operator/()
*/
template<class T> inline Quaternion<T> operator/(T scalar, const Quaternion<T>& quaternion) {
    return {scalar/quaternion.vector(), scalar/quaternion.scalar()};
}

/** @debugoperator{Magnum::Math::Quaternion} */
template<class T> Corrade::Utility::Debug operator<<(Corrade::Utility::Debug debug, const Quaternion<T>& value) {
    debug << "Quaternion({";
    debug.setFlag(Corrade::Utility::Debug::SpaceAfterEachValue, false);
    debug << value.vector().x() << ", " << value.vector().y() << ", " << value.vector().z() << "}, " << value.scalar() << ")";
    debug.setFlag(Corrade::Utility::Debug::SpaceAfterEachValue, true);
    return debug;
}

/* Explicit instantiation for commonly used types */
#ifndef DOXYGEN_GENERATING_OUTPUT
extern template Corrade::Utility::Debug MAGNUM_EXPORT operator<<(Corrade::Utility::Debug, const Quaternion<float>&);
#ifndef MAGNUM_TARGET_GLES
extern template Corrade::Utility::Debug MAGNUM_EXPORT operator<<(Corrade::Utility::Debug, const Quaternion<double>&);
#endif
#endif

}}

#endif