coord/ops: use free functions for normalisations

coord/ops.hpp

@@ -52,6 +52,17 @@ namespace util {
         using result_t = typename result<A,B>::type;
 
 
+        //---------------------------------------------------------------------
+        template <template <size_t,typename> class K>
+        struct has_norm : public std::false_type { };
+
+        template <> struct has_norm<vector>     : public std::true_type { };
+        template <> struct has_norm<quaternion> : public std::true_type { };
+
+        template <template <size_t,typename> class K>
+        constexpr auto has_norm_v = has_norm<K>::value;
+
+
         //---------------------------------------------------------------------
         template <template <size_t,typename> class K>
         struct has_scalar_op : public std::false_type { };
@@ -407,7 +418,79 @@ namespace util {
     }
 
 
+    ///////////////////////////////////////////////////////////////////////////
+    template <
+        size_t S,
+        typename T,
+        template <size_t,typename> class K,
+        typename = std::enable_if_t<coord::has_norm_v<K>,void>
+    >
+    constexpr
+    T
+    norm2 (const K<S,T> &k)
+    {
+        T sum = T{0};
+        for (auto &t: k)
+            sum += t * t;
+        return sum;
+    }
+
+
     //-------------------------------------------------------------------------
+    template <
+        size_t S,
+        typename T,
+        template <size_t,typename> class K,
+        typename = std::enable_if_t<
+            coord::has_norm_v<K>,
+            void
+        >
+    >
+    constexpr
+    T
+    norm (const K<S,T> &k)
+    {
+        return std::sqrt (norm2 (k));
+    }
+
+
+    //-------------------------------------------------------------------------
+    template <
+        size_t S,
+        typename T,
+        template <size_t,typename> class K,
+        typename = std::enable_if_t<
+            coord::has_norm_v<K>,
+            void
+        >
+    >
+    constexpr
+    K<S,T>
+    normalised (const K<S,T> &k)
+    {
+        return k / norm (k);
+    }
+
+
+    //-------------------------------------------------------------------------
+    template <
+        size_t S,
+        typename T,
+        template <size_t,typename> class K,
+        typename = std::enable_if_t<
+            coord::has_norm_v<K>,
+            void
+        >
+    >
+    constexpr
+    bool
+    is_normalised (const K<S,T> &k)
+    {
+        return almost_equal (norm2 (k), T{1});
+    }
+
+
+    ///////////////////////////////////////////////////////////////////////////
     template <
         size_t S,
         typename T,

geom/cylinder.cpp

@@ -28,10 +28,10 @@ cylinder<S,T>::includes (util::point<S,T> p_) const
     auto p_0 = p_ - p0;
 
     auto l = dot (p10, p_0);
-    if (l < 0 || l > p10.magnitude2 ())
+    if (l < 0 || l > norm2 (p10))
         return false;
 
-    auto d = dot (p10, p10) - l * l * p10.magnitude2 ();
+    auto d = dot (p10, p10) - l * l * norm2 (p10);
     if (d > radius * radius)
         return false;
 

geom/plane.cpp

@@ -29,7 +29,7 @@ plane<S,T>::plane (point<S,T>  _p,
     p (_p),
     n (_n)
 {
-    CHECK_EQ (n.magnitude2 (), T{1});
+    CHECK (is_normalised (n));
 }
 
 

geom/ray.cpp

@@ -31,7 +31,7 @@ ray<S,T>::ray (util::point<S,T> _origin,
     origin (_origin),
     direction (_direction)
 {
-    CHECK (direction.is_normalised ());
+    CHECK (is_normalised (direction));
 }
 
 
@@ -43,7 +43,7 @@ ray<S,T>::make (util::point<S,T> origin,
 {
     return {
         origin,
-        (target - origin).normalised ()
+        normalised (target - origin)
     };
 }
 

matrix.cpp

@@ -333,8 +333,8 @@ matrix<S,T>::look_at (util::point<3,T> eye,
                       util::point<3,T> centre,
                       util::vector<3,T> up)
 {
-    const auto f = (centre - eye).normalise ();
+    const auto f = normalised (centre - eye);
-    const auto s = cross (f, up).normalise ();
+    const auto s = normalised (cross (f, up));
     const auto u = cross (s, f);
 
     return { {
@@ -396,7 +396,7 @@ template <size_t S, typename T>
 matrix4<T>
 matrix<S,T>::rotation (T angle, util::vector<3,T> about)
 {
-    about.normalise ();
+    CHECK (is_normalised (about));
 
     T c = std::cos (angle);
     T s = std::sin (angle);

quaternion.cpp

@@ -19,6 +19,7 @@
 
 #include "./debug.hpp"
 #include "./maths.hpp"
+#include "./vector.hpp"
 
 #include <iostream>
 
@@ -33,19 +34,10 @@ template<> const quaternion<4, double> quaternion<4, double>::IDENTITY = { 1, 0,
 
 ///////////////////////////////////////////////////////////////////////////////
 template <size_t S, typename T>
-bool
-quaternion<S,T>::is_normalised (void) const
-{
-    return almost_equal (T{1}, magnitude2 ());
-}
-
-
-//-----------------------------------------------------------------------------
-template <size_t S, typename T>
 quaternion<S,T>
 quaternion<S,T>::rotation (const T radians, const vector<3,T> axis)
 {
-    CHECK (axis.is_normalised ());
+    CHECK (is_normalised (axis));
 
     auto w   = std::cos (radians / 2);
     auto xyz = std::sin (radians / 2) * axis;
@@ -72,36 +64,6 @@ quaternion<S,T>::rotation (const vector<3,T> src, const vector<3,T> dst)
 }
 
 
-///////////////////////////////////////////////////////////////////////////////
-template <size_t S, typename T>
-T
-quaternion<S,T>::magnitude2 (void) const
-{
-    return this->a * this->a +
-           this->b * this->b +
-           this->c * this->c +
-           this->d * this->d;
-}
-
-
-//-----------------------------------------------------------------------------
-template <size_t S, typename T>
-T
-quaternion<S,T>::magnitude (void) const
-{
-    return std::sqrt (magnitude2 ());
-}
-
-
-//-----------------------------------------------------------------------------
-template <size_t S, typename T>
-quaternion<S,T>
-quaternion<S,T>::normalised (void) const
-{
-    return *this / magnitude ();
-}
-
-
 ///////////////////////////////////////////////////////////////////////////////
 template <size_t S, typename T>
 quaternion<S,T>
@@ -123,8 +85,8 @@ 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 (is_normalised (a));
-    CHECK (b.is_normalised ());
+    CHECK (is_normalised (b));
 
     return quaternion<S,T> {
           a.w * b.w + a.x * b.x + a.y * b.y + a.z * b.z,
@@ -142,7 +104,7 @@ template <size_t S, typename T>
 util::matrix4<T>
 quaternion<S, T>::as_matrix (void) const
 {
-    CHECK (is_normalised ());
+    CHECK (is_normalised (*this));
 
     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;
@@ -181,7 +143,7 @@ namespace util { namespace debug {
     struct validator<quaternion<S,T>> {
         static bool is_valid (const quaternion<S,T> &q)
         {
-            return q.is_normalised ();
+            return is_normalised (q);
         }
     };
 } }

quaternion.hpp

@@ -44,12 +44,6 @@ namespace util {
         static quaternion rotation (T radians, vector<3,T> axis);
         static quaternion rotation (vector<3,T> src, vector<3,T> dst);
 
-        T magnitude  (void) const;
-        T magnitude2 (void) const;
-
-        bool is_normalised (void) const;
-        quaternion normalised (void) const;
-
         matrix4<T> as_matrix (void) const;
 
         static const quaternion IDENTITY;

test/coord.cpp

@@ -1,5 +1,6 @@
 #include "tap.hpp"
 #include "point.hpp"
+#include "vector.hpp"
 
 #include "coord/iostream.hpp"
 
@@ -30,5 +31,8 @@ main (void)
         "unary point boolean negation has type bool"
     );
 
+    auto vec = util::vector4f (0.5f);
+    t.expect_eq (vec, util::normalised (vec), "normalisation of normalised vector");
+
     return t.status ();
 }

test/quaternion.cpp

@@ -14,7 +14,7 @@ main (void)
     util::TAP::logger tap;
 
     tap.expect_eq (
-        quaternionf::IDENTITY.magnitude (), 1.f,
+        norm (quaternionf::IDENTITY), 1.f,
         "identity magnitude is unit"
     );
 
@@ -25,7 +25,7 @@ main (void)
     );
 
     {
-        auto val = quaternionf (2, 3, 4, 7).normalised ();
+        auto val = normalised (quaternionf (2, 3, 4, 7));
 
         tap.expect_eq (
             val * quaternionf::IDENTITY,
@@ -35,11 +35,17 @@ main (void)
     }
 
     {
+        // construct two quaternions for a trivial multiplication check.
+        //
+        // we use vector as an intermediate form so that normalisation can be
+        // done before the quaternion code is invoked (which may consider the
+        // values to be invalid and throw an exception given it's geared more
+        // towards rotations than general maths).
         util::vector4f a_v { 2, -11,  5, -17};
         util::vector4f b_v { 3,  13, -7, -19};
 
-        auto a = a_v.normalised ().as<quaternion> ();
+        auto a = normalised (a_v).as<quaternion> ();
-        auto b = b_v.normalised ().as<quaternion> ();
+        auto b = normalised (b_v).as<quaternion> ();
         auto c = quaternionf {
             -0.27358657116960006f,
             -0.43498209092420004f,

test/vector.cpp

@@ -47,7 +47,7 @@ test_polar (util::TAP::logger &tap)
         auto in_cart = t.cartesian;
         auto to_cart = util::polar_to_cartesian (t.polar);
 
-        tap.expect_lt ((in_cart - to_cart).magnitude (), 0.00001f, "%s", t.desc);
+        tap.expect_lt (norm (in_cart - to_cart), 0.00001f, "%s", t.desc);
 
         // Compare polar representations. Make sure to normalise them first.
         auto in_polar = t.polar;
@@ -91,11 +91,11 @@ test_euler (util::TAP::logger &tap)
     for (auto i: TESTS) {
         auto trip = util::from_euler (util::to_euler (i.dir));
         auto diff = i.dir - trip;
-        auto norm = diff.magnitude ();
+        auto mag  = norm (diff);
 
         // trig functions reduce precision above almost_equal levels, so we
         // hard code a fairly low bound here instead.
-        tap.expect_lt (norm, 1e-7, "euler round-trip error, %s", i.name);
+        tap.expect_lt (mag, 1e-7, "euler round-trip error, %s", i.name);
     }
 }
 
@@ -108,8 +108,8 @@ main ()
     test_polar (tap);
     test_euler (tap);
 
-    tap.expect (!util::vector3f::ZERO.is_normalised (), "zero isn't normalised");
+    tap.expect (!is_normalised (util::vector3f::ZERO), "zero isn't normalised");
-    tap.expect (!util::vector3f::UNIT.is_normalised (), "unit is normalised");
+    tap.expect (!is_normalised (util::vector3f::UNIT), "unit is normalised");
 
     return tap.status ();
 }

vector.cpp

@@ -33,28 +33,6 @@ using util::vector3d;
 ///////////////////////////////////////////////////////////////////////////////
 template <size_t S, typename T>
 T
-util::vector<S,T>::magnitude (void) const
-{
-    // TODO: this should not truncate for integral types
-    return static_cast<T> (std::sqrt (magnitude2 ()));
-}
-
-
-//-----------------------------------------------------------------------------
-template <size_t S, typename T>
-T
-util::vector<S,T>::magnitude2 (void) const
-{
-    T total { 0 };
-    for (size_t i = 0; i < S; ++i)
-        total += pow2 (this->data[i]);
-    return total;
-}
-
-
-//-----------------------------------------------------------------------------
-template <size_t S, typename T>
-T
 util::vector<S,T>::difference (vector<S,T> rhs) const
 {
     // TODO: change the signature to ensure it does not truncate
@@ -74,44 +52,6 @@ util::vector<S,T>::difference2 (vector<S,T> rhs) const
 }
 
 
-///////////////////////////////////////////////////////////////////////////////
-template <size_t S, typename T>
-bool
-vector<S,T>::is_normalised (void) const
-{
-    return almost_equal (magnitude2 (), T{1});
-}
-
-
-//-----------------------------------------------------------------------------
-template <size_t S, typename T>
-util::vector<S,T>&
-util::vector<S,T>::normalise (void)
-{
-    T mag = magnitude ();
-
-    for (size_t i = 0; i < S; ++i)
-        this->data[i] /= mag;
-
-    return *this;
-}
-
-
-//-----------------------------------------------------------------------------
-template <size_t S, typename T>
-util::vector<S,T>
-util::vector<S,T>::normalised (void) const
-{
-    T mag = magnitude ();
-    util::vector<S,T> out;
-
-    for (size_t i = 0; i < S; ++i)
-        out.data[i] = this->data[i] / mag;
-
-    return out;
-}
-
-
 ///////////////////////////////////////////////////////////////////////////////
 template <typename T>
 vector<2,T>
@@ -157,10 +97,10 @@ template <typename T>
 vector<2,T>
 util::to_euler (vector<3,T> vec)
 {
-    vec.normalise ();
+    CHECK (is_normalised (vec));
 
     return {
-        std::acos  (vec.y / vec.magnitude ()),
+         std::acos  (vec.y),
         -std::atan2 (vec.z, vec.x),
     };
 }
@@ -204,7 +144,7 @@ template <typename T>
 vector<3,T>
 util::cartesian_to_spherical (vector<3,T> c)
 {
-    T mag = c.magnitude ();
+    T mag = norm (c);
 
     return vector<3,T> {
         mag,

vector.hpp

@@ -32,18 +32,9 @@ namespace util {
         // vector size
         bool is_zero (void) const;
 
-        T magnitude  (void) const;
-        T magnitude2 (void) const;
-
         T difference  (vector<S,T>) const;
         T difference2 (vector<S,T>) const;
 
-        // normalisation
-        bool is_normalised (void) const;
-
-        vector<S,T>& normalise  (void);
-        vector<S,T>  normalised [[gnu::warn_unused_result]] (void) const;
-
         // representations
         template <size_t D> vector<D,T> homog (void) const;