quaternion: update to use coord framework

Makefile.am

@@ -435,6 +435,7 @@ TEST_BIN = 			\
     test/point			\
     test/polynomial		\
     test/pool			\
+    test/quaternion		\
     test/rand/buckets		\
     test/random			\
     test/range 			\

coord/names.hpp

@@ -20,11 +20,23 @@
 namespace util { namespace coord {
     ///////////////////////////////////////////////////////////////////////
     // tags for accessor names
+    //
+    // colours
     struct rgba { };
     struct hsv  { };
+
+    // physical positions
     struct xyzw { };
+
+    // texture coords
     struct stpq { };
+
+    // physical dimensions
     struct whd  { };
+
+    // quaternions
+    struct wxyz { };
+    struct abcd { };
 } }
 
 #endif

coord/ops.hpp

@@ -30,6 +30,7 @@ namespace util {
     template <size_t,typename> struct extent;
     template <size_t,typename> struct vector;
     template <size_t,typename> struct colour;
+    template <size_t,typename> struct quaternion;
 
     ///////////////////////////////////////////////////////////////////////
     // operation traits
@@ -56,6 +57,7 @@ namespace util {
     template <> struct is_coord<extent>     : std::true_type { };
     template <> struct is_coord<vector>     : std::true_type { };
     template <> struct is_coord<colour>     : std::true_type { };
+    template <> struct is_coord<quaternion> : std::true_type { };
 
     template <template <size_t,typename> class K>
     constexpr bool

coord/store.hpp

@@ -168,6 +168,17 @@ namespace util { namespace coord {
             struct { T w,h,d; };
         };
     };
+
+
+    ///////////////////////////////////////////////////////////////////////////
+    template <typename T>
+    struct store<4,T,wxyz,abcd> {
+        union {
+            T data[4];
+            struct { T w,x,y,z; };
+            struct { T a,b,c,d; };
+        };
+    };
 } }
 
 #endif

quaternion.cpp

@@ -15,75 +15,52 @@
  */
 
 
-#include "quaternion.hpp"
+#include "./quaternion.hpp"
 
-#include "debug.hpp"
+#include "./debug.hpp"
-#include "maths.hpp"
+#include "./maths.hpp"
 
+#include <iostream>
 
 //-----------------------------------------------------------------------------
 using util::quaternion;
-using util::matrix4;
 
 
 //-----------------------------------------------------------------------------
-template<> const quaternion< float> quaternion< float>::IDENTITY = { 1, 0, 0, 0 };
+template<> const quaternion<4,  float> quaternion<4,  float>::IDENTITY = { 1, 0, 0, 0 };
-template<> const quaternion<double> quaternion<double>::IDENTITY = { 1, 0, 0, 0 };
+template<> const quaternion<4, double> quaternion<4, double>::IDENTITY = { 1, 0, 0, 0 };
 
 
 ///////////////////////////////////////////////////////////////////////////////
-template <typename T>
+template <size_t S, typename T>
-quaternion<T>::quaternion (T _a, T _b, T _c, T _d):
+bool
-    a (_a),
+quaternion<S,T>::is_normalised (void) const
-    b (_b),
+{
-    c (_c),
+    return almost_equal (T{1}, magnitude2 ());
-    d (_d)
+}
-{ ; }
 
 
 //-----------------------------------------------------------------------------
-template <typename T>
+template <size_t S, typename T>
-quaternion<T>::quaternion (T _a):
+quaternion<S,T>
-    a (_a),
+quaternion<S,T>::rotation (const T radians, const vector<3,T> axis)
-    b (T{}),
+{
-    c (T{}),
-    d (T{})
-{ ; }
-
-
-//-----------------------------------------------------------------------------
-template <typename T>
-quaternion<T>::quaternion ():
-    quaternion (T{}, T{}, T{}, T{})
-{ ; }
-
-
-//-----------------------------------------------------------------------------
-template <typename T>
-quaternion<T>::quaternion (vector3<T> v):
-    quaternion (0, v.x, v.y, v.z)
-{ ; }
-
-
-///////////////////////////////////////////////////////////////////////////////
-template <typename T>
-quaternion<T>
-quaternion<T>::rotation (const T radians, const vector<3,T> axis) {
     CHECK (axis.is_normalised ());
 
+    auto w   = std::cos (radians / 2);
+    auto xyz = std::sin (radians / 2) * axis;
+
     return {
-        std::cos (radians / 2),
+        w, xyz.x, xyz.y, xyz.z
-        std::sin (radians / 2) * axis.x,
-        std::sin (radians / 2) * axis.y,
-        std::sin (radians / 2) * axis.z
     };
 }
 
 
 //-----------------------------------------------------------------------------
-template <typename T>
+template <size_t S, typename T>
-quaternion<T>
+quaternion<S,T>
-quaternion<T>::rotation (vector<3,T> src, vector<3,T> dst) {
+quaternion<S,T>::rotation (const vector<3,T> src, const vector<3,T> dst)
+{
     auto v = util::cross (src, dst);
 
     return {
@@ -96,196 +73,148 @@ quaternion<T>::rotation (vector<3,T> src, vector<3,T> dst) {
 
 
 ///////////////////////////////////////////////////////////////////////////////
-template <typename T>
+template <size_t S, typename T>
 T
-quaternion<T>::norm2 (void) const
+quaternion<S,T>::magnitude2 (void) const
 {
-    return a * a + b * b + c * c + d * d;
+    return this->a * this->a +
+           this->b * this->b +
+           this->c * this->c +
+           this->d * this->d;
 }
 
 
 //-----------------------------------------------------------------------------
-template <typename T>
+template <size_t S, typename T>
 T
-quaternion<T>::norm (void) const
+quaternion<S,T>::magnitude (void) const
 {
-    return std::sqrt (norm2 ());
+    return std::sqrt (magnitude2 ());
 }
 
 
 //-----------------------------------------------------------------------------
-template <typename T>
+template <size_t S, typename T>
-quaternion<T>
+quaternion<S,T>
-quaternion<T>::normalised (void) const
+quaternion<S,T>::normalised (void) const
 {
-    return *this / norm ();
+    return *this / magnitude ();
 }
 
 
 ///////////////////////////////////////////////////////////////////////////////
-template <typename T>
+template <size_t S, typename T>
-quaternion<T>
+quaternion<S,T>
-quaternion<T>::operator- (void) const
+util::operator* (const quaternion<S,T> a, const quaternion<S,T> b)
-{
-    return { -a, -b, -c, -d };
-}
-
-
-//-----------------------------------------------------------------------------
-template <typename T>
-quaternion<T>
-quaternion<T>::conjugate (void) const
-{
-    return { a, -b, -c, -d };
-}
-
-
-///////////////////////////////////////////////////////////////////////////////
-template <typename T>
-quaternion<T>
-quaternion<T>::operator+ (const quaternion<T> q) const
 {
     return {
-        a + q.a,
+        a.w * b.w - a.x * b.x - a.y * b.y - a.z * b.z,
-        b + q.b,
+        a.w * b.x + a.x * b.w + a.y * b.z - a.z * b.y,
-        c + q.c,
+        a.w * b.y - a.x * b.z + a.y * b.w + a.z * b.x,
-        d + q.d
+        a.w * b.z + a.x * b.y - a.y * b.x + a.z * b.w,
+    };
+}
+
+template quaternion<4,float> util::operator* (quaternion<4,float>, quaternion<4,float>);
+
+
+//-----------------------------------------------------------------------------
+template <size_t S, typename T>
+quaternion<S,T>
+util::operator/ (const quaternion<S,T> a, const quaternion<S,T> b)
+{
+    CHECK (a.is_normalised ());
+    CHECK (b.is_normalised ());
+
+    return quaternion<S,T> {
+          a.w * b.w + a.x * b.x + a.y * b.y + a.z * b.z,
+        - a.w * b.x + a.x * b.w + a.y * b.z - a.z * b.y,
+        - a.w * b.y - a.x * b.z + a.y * b.w + a.z * b.x,
+        - a.w * b.z + a.x * b.y - a.y * b.x + a.z * b.w,
+    };
+}
+
+template quaternion<4,float> util::operator/ (quaternion<4,float>, quaternion<4,float>);
+
+
+///////////////////////////////////////////////////////////////////////////////
+template <size_t S, typename T>
+quaternion<S,T>
+util::operator* (const quaternion<S,T> q, const T t)
+{
+    return {
+        q.w * t,
+        q.x * t,
+        q.y * t,
+        q.z * t
     };
 }
 
 
 //-----------------------------------------------------------------------------
-template <typename T>
+template <size_t S, typename T>
-quaternion<T>
+quaternion<S,T>
-quaternion<T>::operator- (const quaternion<T> q) const
+util::operator/ (const quaternion<S,T> q, const T t)
 {
     return {
-        a - q.a,
+        q.w / t,
-        b - q.b,
+        q.x / t,
-        c - q.c,
+        q.y / t,
-        d - q.d
+        q.z / t
     };
 }
 
 
-//-----------------------------------------------------------------------------
-template <typename T>
-quaternion<T>
-quaternion<T>::operator* (const quaternion<T> q) const {
-    return {
-        a * q.a - b * q.b - c * q.c - d * q.d,
-        a * q.b + b * q.a + c * q.d - d * q.c,
-        a * q.c - b * q.d + c * q.a + d * q.b,
-        a * q.d + b * q.c - c * q.b + d * q.a,
-    };
-}
-
-
-//-----------------------------------------------------------------------------
-template <typename T>
-quaternion<T>
-quaternion<T>::operator/ (const quaternion<T> q) const
-{
-    auto n = q.norm2 ();
-
-    return {
-        (  a * q.a + b * q.b + c * q.c + d * q.d) / n,
-        (- a * q.b + b * q.a + c * q.d - d * q.c) / n,
-        (- a * q.c - b * q.d + c * q.a + d * q.b) / n,
-        (- a * q.d + b * q.c - c * q.b + d * q.a) / n
-    };
-}
-
 ///////////////////////////////////////////////////////////////////////////////
-template <typename T>
+template <size_t S, typename T>
-quaternion<T>
+util::matrix4<T>
-quaternion<T>::operator+ (const T t) const
+quaternion<S, T>::as_matrix (void) const
 {
-    return { a + t, b, c, d };
+    CHECK (is_normalised ());
-}
 
-
+    const T wx = this->w * this->x, wy = this->w * this->y, wz = this->w * this->z;
-//-----------------------------------------------------------------------------
+    const T xx = this->x * this->x, xy = this->x * this->y, xz = this->x * this->z;
-template <typename T>
+    const T yy = this->y * this->y, yz = this->y * this->z, zz = this->z * this->z;
-quaternion<T>
-quaternion<T>::operator- (const T t) const
-{
-    return { a - t, b, c, d };
-}
-
-
-//-----------------------------------------------------------------------------
-template <typename T>
-quaternion<T>
-quaternion<T>::operator* (const T t) const
-{
-    return { a * t, b, c, d };
-}
-
-
-//-----------------------------------------------------------------------------
-template <typename T>
-quaternion<T>
-quaternion<T>::operator/ (const T t) const
-{
-    return { a / t, b, c, d };
-}
-
-
-///////////////////////////////////////////////////////////////////////////////
-template <typename T>
-bool
-quaternion<T>::operator== (const quaternion<T> rhs) const
-{
-    return almost_equal (a, rhs.a) &&
-           almost_equal (b, rhs.b) &&
-           almost_equal (c, rhs.c) &&
-           almost_equal (d, rhs.d);
-}
-
-
-//-----------------------------------------------------------------------------
-template <typename T>
-bool
-quaternion<T>::operator!= (const quaternion<T> rhs) const
-{
-    return !(*this == rhs);
-}
-
-
-///////////////////////////////////////////////////////////////////////////////
-template <typename T>
-matrix4<T>
-quaternion<T>::rotation_matrix (void) const
-{
-    CHECK_EQ (T{1}, norm ());
-
-    const T wx = w * x, wy = w * y, wz = w * z;
-    const T xx = x * x, xy = x * y, xz = x * z;
-    const T yy = y * y, yz = y * z, zz = z * z;
 
     return { {
         { 1 - 2 * (yy - zz),     2 * (xy - wz),     2 * (xz + wy), 0 },
         {     2 * (xy + wz), 1 - 2 * (xx - zz),     2 * (yz - wx), 0 },
-        {     2 * (xz - wy),     2 * (yz + wx), 1 - 2 * (xx - yy), 0 },
+        {     2 * (xz - wy),     2 * (yz + wx), 1 - 2 * (xx + yy), 0 },
         {                 0,                 0,                 0, 1 }
     } };
 }
 
 
 ///////////////////////////////////////////////////////////////////////////////
-template <typename T>
+template <size_t S, typename T>
 std::ostream&
-util::operator<< (std::ostream &os, quaternion<T> q)
+util::operator<< (std::ostream &os, const quaternion<S,T> q)
 {
-    os << q.w << ' ' << q.x << "i " << q.y << "j " << q.z << 'k';
+    return os << "[" << q.w << ", " << q.x << ", " << q.y << ", " << q.z << "]";
-    return os;
 }
 
 
 ///////////////////////////////////////////////////////////////////////////////
-template struct util::quaternion<float>;
+template struct util::quaternion<4,float>;
-template struct util::quaternion<double>;
+template struct util::quaternion<4,double>;
 
-template std::ostream& util::operator<< (std::ostream&, quaternion<float>);
+
-template std::ostream& util::operator<< (std::ostream&, quaternion<double>);
+template std::ostream& util::operator<< (std::ostream&, quaternion<4,float>);
+template std::ostream& util::operator<< (std::ostream&, quaternion<4,double>);
+
+
+///////////////////////////////////////////////////////////////////////////////
+namespace util { namespace debug {
+    template <size_t S, typename T>
+    struct validator<quaternion<S,T>> {
+        static bool is_valid (const quaternion<S,T> &q)
+        {
+            return q.is_normalised ();
+        }
+    };
+} }
+
+
+//-----------------------------------------------------------------------------
+template bool util::debug::is_valid(const quaternion<4,float>&);
+template bool util::debug::is_valid(const quaternion<4,double>&);

quaternion.hpp

@@ -17,6 +17,8 @@
 #ifndef __UTIL_QUATERNION_HPP
 #define __UTIL_QUATERNION_HPP
 
+#include "coord.hpp"
+
 #include "vector.hpp"
 #include "matrix.hpp"
 
@@ -24,53 +26,46 @@
 
 
 namespace util {
-    template <typename T>
+    // quaternions must be 4 elements, but we include a size parameter so it
-    struct quaternion {
+    // fits with the generic coord infrastructure more easily.
-        ///////////////////////////////////////////////////////////////////////
+    //
-        union {
+    // specifically: 
-            struct { T w, x, y, z; };
+    // large regions of base code require a template template parameter with
-            struct { T a, b, c, d; };
+    // size and type arguments, which is annoying to work around for this one
-            T data[4];
+    // case.
-        };
+    //
+    // we protect against invalid instantiations through static_assert
+    template <size_t S, typename T>
+    struct quaternion : public coord::base<4,T,quaternion,coord::wxyz,coord::abcd> {
+        static_assert (S == 4, "quaternions must be 4 elements");
 
-        static const quaternion IDENTITY;
+        using coord::base<S,T,::util::quaternion,::util::coord::wxyz,::util::coord::abcd>::base;
 
         static quaternion rotation (T radians, vector<3,T> axis);
         static quaternion rotation (vector<3,T> src, vector<3,T> dst);
 
-        quaternion (T a, T b, T c, T d);
+        T magnitude  (void) const;
-        quaternion (T a);
+        T magnitude2 (void) const;
-        quaternion ();
-        quaternion (vector<3,T>);
 
-        T norm  (void) const;
+        bool is_normalised (void) const;
-        T norm2 (void) const;
         quaternion normalised (void) const;
 
-        quaternion operator- (void) const;
+        matrix4<T> as_matrix (void) const;
-        quaternion conjugate (void) const;
 
-        quaternion operator+ (const quaternion) const;
+        static const quaternion IDENTITY;
-        quaternion operator- (const quaternion) const;
-        quaternion operator* (const quaternion) const;
-        quaternion operator/ (const quaternion) const;
-
-        quaternion operator+ (const T) const;
-        quaternion operator- (const T) const;
-        quaternion operator* (const T) const;
-        quaternion operator/ (const T) const;
-
-        bool operator== (const quaternion) const;
-        bool operator!= (const quaternion) const;
-
-        matrix4<T> rotation_matrix (void) const;
     };
 
-    typedef quaternion<float> quaternionf;
+    template <size_t S, typename T> quaternion<S,T> operator* (const quaternion<S,T>, const quaternion<S,T>);
+    template <size_t S, typename T> quaternion<S,T> operator/ (const quaternion<S,T>, const quaternion<S,T>);
 
-    template <typename T>
+    template <size_t S, typename T> quaternion<S,T> operator* (const quaternion<S,T>, const T);
+    template <size_t S, typename T> quaternion<S,T> operator/ (const quaternion<S,T>, const T);
+
+    typedef quaternion<4,float> quaternionf;
+
+    template <size_t S, typename T>
     std::ostream&
-    operator<< (std::ostream&, quaternion<T>);
+    operator<< (std::ostream&, quaternion<S,T>);
 }
 
 

test/quaternion.cpp

@@ -0,0 +1,90 @@
+#include <quaternion.hpp>
+
+#include "tap.hpp"
+#include "types.hpp"
+
+using util::quaternion;
+using util::quaternionf;
+
+
+///////////////////////////////////////////////////////////////////////////////
+int
+main (void)
+{
+    util::TAP::logger tap;
+
+    tap.expect_eq (
+        quaternionf::IDENTITY.magnitude (), 1.f,
+        "identity magnitude is unit"
+    );
+
+    tap.expect_eq (
+        quaternionf::IDENTITY * quaternionf::IDENTITY,
+        quaternionf::IDENTITY,
+        "identity multiplication with identity"
+    );
+
+    {
+        auto val = quaternionf (2, 3, 4, 7).normalised ();
+
+        tap.expect_eq (
+            val * quaternionf::IDENTITY,
+            val,
+            "identity multiplication with quaternion constant"
+        );
+    }
+
+    {
+        util::vector4f a_v { 2, -11,  5, -17};
+        util::vector4f b_v { 3,  13, -7, -19};
+
+        auto a = a_v.normalised ().as<quaternion> ();
+        auto b = b_v.normalised ().as<quaternion> ();
+        auto c = quaternionf {
+            -0.27358657116960006f,
+            -0.43498209092420004f,
+            -0.8443769181970001f,
+            -0.15155515799559996f,
+        };
+
+        tap.expect_eq  (a * b, c, "multiplication");
+        tap.expect_neq (b * a, c, "multiplication is not commutative");
+    }
+
+    tap.expect_eq (
+        quaternionf::IDENTITY.as_matrix (),
+        util::matrix4f::IDENTITY,
+        "identity quaternion to matrix"
+    );
+
+    {
+        static const struct {
+            float mag;
+            util::vector3f axis;
+        } ROTATIONS[] = {
+            { 0.f, { 1.f, 0.f, 0.f } },
+            { 1.f, { 0.f, 1.f, 0.f } },
+            { 0.f, { 0.f, 0.f, 1.f } },
+        };
+
+        for (size_t i = 0; i < elems (ROTATIONS); ++i) {
+            const auto &r = ROTATIONS[i];
+
+            tap.expect_eq (quaternionf::rotation    (r.mag, r.axis).as_matrix (),
+                           util::matrix4f::rotation (r.mag, r.axis),
+                           "single basis rotation %zu", i);
+        }
+
+        auto q_total = quaternionf::IDENTITY;
+        auto m_total = util::matrix4f::IDENTITY;
+
+        for (auto r: ROTATIONS) {
+            q_total = q_total.rotation (r.mag, r.axis) * q_total;
+            m_total = m_total.rotation (r.mag, r.axis) * m_total;
+        }
+
+        tap.expect_eq (q_total.as_matrix (), m_total, "chained single axis rotations");
+    }
+
+    return tap.status ();
+}