noise: support n-dimensional noise

Makefile.am

@@ -158,6 +158,10 @@ UTIL_FILES   =			\
     noise/basis/perlin.ipp	\
     noise/basis/runtime.cpp	\
     noise/basis/runtime.hpp	\
+    noise/basis/type/distance.cpp	\
+    noise/basis/type/distance.hpp	\
+    noise/basis/type/gradient.hpp	\
+    noise/basis/type/gradient.cpp	\
     noise/basis/value.hpp	\
     noise/basis/value.ipp	\
     noise/basis/worley.hpp	\

noise/basis/constant.cpp

@@ -19,14 +19,17 @@
 using util::noise::basis::constant;
 
 
-//-----------------------------------------------------------------------------
+///////////////////////////////////////////////////////////////////////////////
-template <typename T>
+template <size_t S, typename T>
-constant<T>::constant (seed_t _seed):
+constant<S,T>::constant (seed_t _seed):
     seed (_seed),
     value (42)
 { ; }
 
 
 ///////////////////////////////////////////////////////////////////////////////
-template struct util::noise::basis::constant<float>;
+template struct util::noise::basis::constant<2,float>;
-template struct util::noise::basis::constant<double>;
+template struct util::noise::basis::constant<2,double>;
+
+template struct util::noise::basis::constant<3,float>;
+template struct util::noise::basis::constant<3,double>;

noise/basis/constant.hpp

@@ -22,12 +22,13 @@
 #include <cstdint>
 
 namespace util { namespace noise { namespace basis {
-    template <typename T>
+    template <size_t S, typename T>
     struct constant {
         using seed_t = uint64_t;
 
         constant (seed_t);
-        T operator() (util::point<2,T>) const;
+
+        T operator() (point<S,T>) const;
 
         seed_t seed;
         T value;

noise/basis/constant.ipp

@@ -20,9 +20,9 @@
 #define __UTIL_NOISE_BASIS_CONSTANT_IPP
 
 namespace util { namespace noise { namespace basis {
-    template <typename T>
+    template <size_t S, typename T>
     T
-    constant<T>::operator() (util::point<2,T>) const
+    constant<S,T>::operator() (point<S,T>) const
     {
         return value;
     }

noise/basis/expgrad.hpp

@@ -24,10 +24,22 @@
 #include "../../point.hpp"
 #include "../../range.hpp"
 
+
+///////////////////////////////////////////////////////////////////////////////
+// modifies a standard uniforma gradient generator to follow an exponential
+// distribution, base^(-t*exp)
+//
+// 
 namespace util { namespace noise { namespace basis {
-    template <typename T, lerp_t<T> L>
+    template <
-    struct expgrad : public gradient<T,L> {
+        size_t S,   // probe point dimensionality
-        explicit expgrad <T,L> (seed_t seed, T base = (T)1.02, T exponent = T{256});
+        typename T, // probe point value_type
+        lerp_t<T> L // axis interpolation function
+    >
+    struct expgrad : public gradient<S,T,L> {
+        explicit expgrad (seed_t seed,
+                          T base = (T)1.02,
+                          T exponent = T{256});
 
         T base (void) const;
         T base (T);
@@ -36,7 +48,7 @@ namespace util { namespace noise { namespace basis {
         T exponent (T);
 
     protected:
-        vector<2,T> generate (point<2,intmax_t>) const;
+        vector<S,T> generate (pointi<S>) const;
 
         T m_base;
         T m_exponent;

noise/basis/expgrad.ipp

@@ -24,21 +24,22 @@
 
 namespace util { namespace noise { namespace basis {
     ///////////////////////////////////////////////////////////////////////////
-    template <typename T, util::noise::lerp_t<T> L>
+    template <size_t S, typename T, util::noise::lerp_t<T> L
-    expgrad<T,L>::expgrad (seed_t _seed, T _base, T _exponent):
+    >
-        gradient<T,L> (_seed),
+    expgrad<S,T,L>::expgrad (seed_t _seed, T _base, T _exponent):
+        gradient<S,T,L> (_seed),
         m_base (_base),
         m_exponent (_exponent)
     { ; }
 
 
     ///////////////////////////////////////////////////////////////////////////
-    template <typename T, util::noise::lerp_t<T> L>
+    template <size_t S, typename T, util::noise::lerp_t<T> L>
-    vector<2,T>
+    vector<S,T>
-    expgrad<T,L>::generate (point<2,intmax_t> p) const
+    expgrad<S,T,L>::generate (pointi<S> p) const
     {
         auto t = (noise::rand<float> (this->seed (), p) + 1) / 2;
         auto factor = std::pow (m_base, -t * m_exponent);
-        return factor * gradient<T,L>::generate (p);
+        return factor * gradient<S,T,L>::generate (p);
     }
 } } }

noise/basis/gradient.hpp

@@ -24,7 +24,11 @@
 
 namespace util { namespace noise { namespace basis {
     /// Perlin: interpolated value across each grid space
-    template <typename T, lerp_t<T> L>
+    template <
+        size_t S,   // probe point dimensionality
+        typename T, // probe point value_type
+        lerp_t<T> L // axis interpolation function
+    >
     struct gradient {
         gradient (seed_t);
 
@@ -32,7 +36,7 @@ namespace util { namespace noise { namespace basis {
         seed_t seed (seed_t);
 
     protected:
-        vector<2,T> generate (point<2,intmax_t>) const;
+        vector<S,T> generate (pointi<S>) const;
 
         seed_t m_seed;
     };

noise/basis/gradient.ipp

@@ -23,36 +23,36 @@
 
 namespace util { namespace noise { namespace basis {
     ///////////////////////////////////////////////////////////////////////////
-    template <typename T, util::noise::lerp_t<T> L>
+    template <size_t S, typename T, util::noise::lerp_t<T> L>
-    gradient<T,L>::gradient (seed_t _seed):
+    gradient<S,T,L>::gradient (seed_t _seed):
         m_seed (_seed)
     { ; }
 
 
     //-------------------------------------------------------------------------
-    template <typename T, lerp_t<T> L>
+    template <size_t S, typename T, lerp_t<T> L>
     seed_t
-    gradient<T,L>::seed (void) const
+    gradient<S,T,L>::seed (void) const
     {
         return m_seed;
     }
 
 
     //-------------------------------------------------------------------------
-    template <typename T, lerp_t<T> L>
+    template <size_t S, typename T, lerp_t<T> L>
     seed_t
-    gradient<T,L>::seed (seed_t _seed)
+    gradient<S,T,L>::seed (seed_t _seed)
     {
         return m_seed = _seed;
     }
 
 
     //-------------------------------------------------------------------------
-    template <typename T, util::noise::lerp_t<T> L>
+    template <size_t S, typename T, util::noise::lerp_t<T> L>
-    vector<2,T>
+    vector<S,T>
-    gradient<T,L>::generate (point<2,intmax_t> p) const
+    gradient<S,T,L>::generate (pointi<S> p) const
     {
-        return noise::rand<2,T> (m_seed, p);
+        return noise::rand<vector,T> (m_seed, p);
     }
 } } }
 

noise/basis/patch.hpp

@@ -17,16 +17,19 @@
 #ifndef __UTIL_NOISE_BASIS_PATCH_HPP
 #define __UTIL_NOISE_BASIS_PATCH_HPP
 
+#include "./type/distance.hpp"
+
 #include "../fwd.hpp"
 #include "../../point.hpp"
 
 namespace util { namespace noise { namespace basis {
-    template <typename T>
+    template <size_t S, typename T>
-    struct patch {
+    struct patch : public type::distance<S,2> {
         patch (seed_t, T width = 0);
         
         range<T> bounds (void) const;
-        T operator() (point2<T>) const;
+
+        T operator() (point<S,T>) const;
         
         seed_t seed (void) const;
         seed_t seed (seed_t);
@@ -35,9 +38,6 @@ namespace util { namespace noise { namespace basis {
         T width (T);
 
     private:
-        point2<T> centroid (util::point2i) const;
-        T generate (util::point2i) const;
-
         static constexpr T THRESHOLD = 1 - T(0.999);
 
         T m_width;

noise/basis/patch.ipp

@@ -26,8 +26,8 @@
 
 namespace util { namespace noise { namespace basis {
     ///////////////////////////////////////////////////////////////////////////
-    template <typename T>
+    template <size_t S, typename T>
-    patch<T>::patch (seed_t _seed, T _width):
+    patch<S,T>::patch (seed_t _seed, T _width):
         m_width (_width),
         m_power (exactly_zero (_width)
                  ? std::numeric_limits<T>::infinity ()
@@ -37,43 +37,40 @@ namespace util { namespace noise { namespace basis {
 
 
     ///////////////////////////////////////////////////////////////////////////
-    template <typename T>
+    template <size_t S, typename T>
     range<T>
-    patch<T>::bounds (void) const
+    patch<S,T>::bounds (void) const
     {
         return { T{0}, T{1} };
     }
 
 
     ///////////////////////////////////////////////////////////////////////////
-    template <typename T>
+    template <size_t S, typename T>
     T
-    patch<T>::operator () (point2<T> p) const
+    patch<S,T>::operator () (point<S,T> p) const
     {
+        static const size_t COUNT = type::distance<S,2>::OFFSET_SIZE;
+
         // extract integer and fractional parts. be careful to always round down
         auto p_int = floor (p).template cast<intmax_t> ();
         auto p_rem = (p - p_int).template as<point> ();
 
-        static const util::vector2i OFFSETS[] = {
+        // find the distances to each neighbour's centroid.
-                                    {  0, -2 },
+        util::point<S,T> centres[COUNT];
-                        { -1, -1 }, {  0, -1 }, {  1, -1 },
+        std::transform (std::begin (this->OFFSETS),
-            { -2,  0 }, { -1,  0 }, {  0,  0 }, {  1,  0 }, {  2,  0 },
+                        std::end   (this->OFFSETS),
-                        { -1,  1 }, {  0,  1 }, {  1,  1 },
+                        std::begin (centres),
-                                    {  0,  2 },
+                        [this,p_int] (auto i) { return (noise::rand<point,T> (m_seed, p_int + i) + 1) / 2  + i; });
-        };
-
-        static const size_t COUNT = elems (OFFSETS);
-
-        // find the distances to each neighbour's centroid
-        util::point2<T> centres[COUNT];
-        for (size_t i = 0; i < COUNT; ++i)
-            centres[i] = centroid (p_int + OFFSETS[i]) + OFFSETS[i];
 
         T distances[COUNT];
-        for (size_t i = 0; i < COUNT; ++i)
+        std::transform (std::begin (centres),
-            distances[i] = std::sqrt (util::distance2 (p_rem, centres[i]));
+                        std::end   (centres),
+                        std::begin (distances),
+                        [p_rem] (auto i) { return util::distance (p_rem, i); });
 
-        // sort the distances, using indices so we can use 'offsets' to generate values
+        // sort the distances using indices so we can reuse indices into
+        // 'OFFSETS' to generate the random patch values
         unsigned indices[COUNT];
         std::iota (std::begin (indices), std::end (indices), 0);
 
@@ -83,12 +80,13 @@ namespace util { namespace noise { namespace basis {
             return distances[a] < distances[b];
         });
 
-        // calculate normalisation constants for the 9 nearest points. the
+        // setup normalisation for the distances to the nearest points. the
         // neighbourhood size is implicitly specified by the 1.5 unit maximum
         // distance.
-        constexpr auto MAX_DISTANCE = 2.1213203435596424f; // std::hypot (1.5f, 1.5f);
+        const size_t hi_off = pow(3,S);
+
         const auto lo = distances[indices[0     ]];
-        const auto hi = std::min (distances[indices[COUNT-1]], MAX_DISTANCE);
+        const auto hi = distances[indices[hi_off]];
 
         T out  = 0.f;
         T sumw = 0.f;
@@ -96,9 +94,13 @@ namespace util { namespace noise { namespace basis {
         // sum the weight values of each neighbour. weight by a function of
         // the distance. we use an power function which allows a known width
         // to blend.
-        for (size_t i = 0; i < COUNT && distances[indices[i]] <= MAX_DISTANCE; ++i)
+        for (size_t i = 0; i < hi_off; ++i)
         {
-            auto v = generate (p_int + OFFSETS[indices[i]]);
+            auto v = util::noise::rand<T> (
+                m_seed,
+                p_int + this->OFFSETS[indices[i]]
+            );
+
             auto d = (distances[indices[i]] - lo) / (hi - lo);
             auto w = std::pow (1 - d, m_power);
 
@@ -111,36 +113,36 @@ namespace util { namespace noise { namespace basis {
 
 
     ///////////////////////////////////////////////////////////////////////////
-    template <typename T>
+    template <size_t S, typename T>
     seed_t
-    patch<T>::seed (void) const
+    patch<S,T>::seed (void) const
     {
         return m_seed;
     }
 
 
     //-------------------------------------------------------------------------
-    template <typename T>
+    template <size_t S, typename T>
     seed_t
-    patch<T>::seed (util::noise::seed_t _seed)
+    patch<S,T>::seed (util::noise::seed_t _seed)
     {
         return m_seed = _seed;
     }
 
 
     ///////////////////////////////////////////////////////////////////////////
-    template <typename T>
+    template <size_t S, typename T>
     T
-    patch<T>::width (void) const
+    patch<S,T>::width (void) const
     {
         return m_width;
     }
 
 
     //-------------------------------------------------------------------------
-    template <typename T>
+    template <size_t S, typename T>
     T
-    patch<T>::width (T _width)
+    patch<S,T>::width (T _width)
     {
         m_width = _width;
         m_power = exactly_zero (_width)
@@ -149,36 +151,4 @@ namespace util { namespace noise { namespace basis {
 
         return m_width;
     }
-
-
-    ///////////////////////////////////////////////////////////////////////////
-    template <typename T>
-    util::point2<T>
-    patch<T>::centroid (util::point2i p) const
-    {
-        using util::hash::murmur2::mix;
-
-        auto u = mix (m_seed, mix (uint64_t (p.x), uint64_t (p.y)));
-        auto v = mix (u, m_seed);
-
-        auto r = util::point<2,T> {
-            (u & 0xffff) / T{0xffff},
-            (v & 0xffff) / T{0xffff}
-        };
-
-        CHECK_LIMIT (r, T{0}, T{1});
-        return r;
-    }
-
-
-    //-------------------------------------------------------------------------
-    template <typename T>
-    T
-    patch<T>::generate (util::point2i p) const
-    {
-        using util::hash::murmur2::mix;
-
-        auto u = mix (m_seed, mix (uint64_t (p.x), uint64_t (p.y)));
-        return (u & 0xffff) / T{0xffff};
-    }
 } } }

noise/basis/perlin.hpp

@@ -19,6 +19,7 @@
 #define __UTIL_NOISE_BASIS_PERLIN_HPP
 
 #include "./gradient.hpp"
+#include "./type/gradient.hpp"
 
 #include "../fwd.hpp"
 #include "../../point.hpp"
@@ -27,18 +28,21 @@
 namespace util { namespace noise { namespace basis {
     /// Perlin: interpolated value across each grid space
     template <
+        size_t S,       // point point dimensionality
         typename T,     // arithmetic and result value_type, must be floating point
         lerp_t<T> L,    // gradient interpolation function
         template <      // gradient provider class, must provide generate(point_t)
+            size_t,
             typename,
             lerp_t<T>
         > class G = gradient
     >
-    struct perlin : public G<T,L> {
+    struct perlin : public G<S,T,L>, public type::gradient<S> {
         perlin (seed_t);
 
         range<T> bounds (void) const;
-        T operator() (point<2,T>) const;
+
+        T operator() (point<S,T>) const;
     };
 } } }
 

noise/basis/perlin.ipp

@@ -20,19 +20,21 @@
 #define __UTIL_NOISE_BASIS_PERLIN_IPP
 
 #include "../rand.hpp"
+#include "../../types.hpp"
+
 
 namespace util { namespace noise { namespace basis {
     ///////////////////////////////////////////////////////////////////////////
-    template <typename T, util::noise::lerp_t<T> L, template <typename,lerp_t<T>> class G>
+    template <size_t S, typename T, util::noise::lerp_t<T> L, template <size_t,typename,lerp_t<T>> class G>
-    perlin<T,L,G>::perlin (seed_t _seed):
+    perlin<S,T,L,G>::perlin (seed_t _seed):
-        G<T,L> (_seed)
+        G<S,T,L> (_seed)
     { ; }
 
 
     //-------------------------------------------------------------------------
-    template <typename T, util::noise::lerp_t<T> L, template <typename,lerp_t<T>> class G>
+    template <size_t S, typename T, util::noise::lerp_t<T> L, template <size_t,typename,lerp_t<T>> class G>
     util::range<T>
-    perlin<T,L,G>::bounds (void) const
+    perlin<S,T,L,G>::bounds (void) const
     {
         return {
             -std::sqrt (T{2}) / 2,
@@ -42,37 +44,38 @@ namespace util { namespace noise { namespace basis {
 
 
     //-------------------------------------------------------------------------
-    template <typename T, util::noise::lerp_t<T> L, template <typename,lerp_t<T>> class G>
+    template <size_t S, typename T, util::noise::lerp_t<T> L, template <size_t,typename,lerp_t<T>> class G>
     T
-    perlin<T,L,G>::operator() (util::point<2,T> p) const
+    perlin<S,T,L,G>::operator() (util::point<S,T> p) const
     {
         // extract integer and fractional parts. be careful to always round down
         auto p_int = floor (p).template cast<intmax_t> ();
         auto p_rem = p - p_int;
 
         // generate the corner positions
-        auto p0 = p_int + util::vector<2,intmax_t> { 0, 0 };
+        pointi<S> p_[pow(2,S)];
-        auto p1 = p_int + util::vector<2,intmax_t> { 1, 0 };
+        std::transform (std::begin (this->CORNERS), std::end (this->CORNERS),
-        auto p2 = p_int + util::vector<2,intmax_t> { 0, 1 };
+                        std::begin (p_),
-        auto p3 = p_int + util::vector<2,intmax_t> { 1, 1 };
+                        [p_int] (auto i) { return i + p_int; });
 
         // generate the corner gradients
-        auto g0 = G<T,L>::generate (p0);
+        vector<S,T> g_[pow(2,S)];
-        auto g1 = G<T,L>::generate (p1);
+        std::transform (std::begin (p_), std::end (p_),
-        auto g2 = G<T,L>::generate (p2);
+                        std::begin (g_),
-        auto g3 = G<T,L>::generate (p3);
+                        [this] (auto i) { return this->generate (i); });
 
         // compute the dot products
-        T v0 = dot (g0, p - p0);
+        T v_[pow(2,S)];
-        T v1 = dot (g1, p - p1);
+        for (size_t i = 0; i < elems (v_); ++i)
-        T v2 = dot (g2, p - p2);
+            v_[i] = dot (g_[i], p - p_[i]);
-        T v3 = dot (g3, p - p3);
 
         // interpolate the results
-        auto L0 = L (v0, v1, p_rem.x);
+        T l_[pow(2,S)];
-        auto L1 = L (v2, v3, p_rem.x);
+        std::copy (std::begin (v_), std::end (v_), std::begin (l_));
-        auto L_ = L (L0, L1, p_rem.y);
 
-        return L_;
+        for (size_t i = S; i; --i)
+            for (size_t j = 0; j < std::pow(2,i); j += 2)
+                l_[j / 2] = L (l_[j], l_[j+1], p_rem[S-i]);
+        return l_[0];
     }
 } } }

noise/basis/runtime.cpp

@@ -17,4 +17,4 @@
 #include "runtime.hpp"
 
 
-template struct util::noise::basis::runtime<float>;
+template struct util::noise::basis::runtime<2,float>;

noise/basis/runtime.hpp

@@ -25,7 +25,7 @@
 #include <memory>
 
 namespace util { namespace noise { namespace basis {
-    template <typename T>
+    template <size_t S, typename T>
     struct runtime {
     public:
         runtime (seed_t) {}
@@ -34,7 +34,7 @@ namespace util { namespace noise { namespace basis {
         runtime& operator= (const runtime&) = delete;
 
         // basis functions
-        T operator () (util::point<2,T> p) const { return (*m_child) (p); }
+        T operator () (util::point<S,T> p) const { return (*m_child) (p); }
         range<T> bounds (void) const { return m_child->bounds (); }
 
         seed_t seed (void) const { return m_child->seed (); }
@@ -43,7 +43,7 @@ namespace util { namespace noise { namespace basis {
     private:
         struct base {
             virtual ~base () = default;
-            virtual T operator() (util::point<2,T>) const = 0;
+            virtual T operator() (util::point<S,T>) const = 0;
             virtual range<T> bounds (void) const = 0;
             virtual seed_t seed (void) const = 0;
             virtual seed_t seed (seed_t) = 0;
@@ -53,7 +53,7 @@ namespace util { namespace noise { namespace basis {
         struct child : public base {
             template <typename ...Args>
             child (seed_t _seed, Args&& ...args): data (_seed, std::forward<Args> (args)...) { }
-            virtual T operator() (util::point<2,T> p) const override { return data (p); }
+            virtual T operator() (util::point<S,T> p) const override { return data (p); }
             virtual range<T> bounds (void) const override { return data.bounds (); }
             virtual seed_t seed (void) const override { return data.seed (); }
             virtual seed_t seed (seed_t _seed) override { return data.seed (_seed); }

noise/basis/type/distance.cpp

@@ -0,0 +1,43 @@
+#include "./distance.hpp"
+#include "../../../extent.hpp"
+
+
+///////////////////////////////////////////////////////////////////////////////
+template <size_t S, size_t R>
+static const std::array<util::vectori<S>, util::pow(R*2+1,S)>
+generate (void)
+{
+    static const util::extent_range<
+        S,typename util::vectori<S>::value_type
+    > area (util::extent<S,typename util::vectori<S>::value_type> {R*2+1});
+
+    std::array<
+        util::vectori<S>,
+        util::pow(R*2+1,S)
+    > out;
+
+    std::transform (area.begin (), area.end (),
+                    out.begin (),
+                    [] (auto i) { return i.template as<util::vector> () - R; });
+
+    return out;
+}
+
+
+///////////////////////////////////////////////////////////////////////////////
+template <size_t S, size_t R>
+const std::array<
+    util::vectori<S>,
+    util::pow(R*2+1,S)
+>
+util::noise::basis::type::distance<S,R>::OFFSETS = generate<S,R> ();
+
+
+//-----------------------------------------------------------------------------
+template struct util::noise::basis::type::distance<1,1>;
+template struct util::noise::basis::type::distance<2,1>;
+template struct util::noise::basis::type::distance<3,1>;
+
+template struct util::noise::basis::type::distance<1,2>;
+template struct util::noise::basis::type::distance<2,2>;
+template struct util::noise::basis::type::distance<3,2>;

noise/basis/type/distance.hpp

@@ -0,0 +1,22 @@
+
+#ifndef __UTIL_NOISE_BASIS_TYPE_HPP
+#define __UTIL_NOISE_BASIS_TYPE_HPP
+
+#include "../../../vector.hpp"
+#include "../../../maths.hpp"
+
+#include <array>
+
+namespace util { namespace noise { namespace basis { namespace type {
+    template <size_t S, size_t R>
+    struct distance {
+    protected:
+        static constexpr size_t OFFSET_SIZE = util::pow(R*2+1,S);
+
+        static const std::array<
+            vectori<S>, util::pow(R*2+1,S)
+        > OFFSETS;
+    };
+} } } }
+
+#endif

noise/basis/type/gradient.cpp

@@ -0,0 +1,44 @@
+/*
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ *
+ * Copyright 2015 Danny Robson <danny@nerdcruft.net>
+ */
+
+#include "./gradient.hpp"
+
+
+///////////////////////////////////////////////////////////////////////////////
+template <size_t S>
+static std::array<util::vectori<S>,util::pow(2,S)>
+generate (void)
+{
+    std::array<util::vectori<S>,util::pow(2,S)> out;
+
+    for (size_t i = 0; i < util::pow(2,S); ++i)
+        for (size_t s = 0; s < S; ++s)
+            out[i][s] = (i >> s) & 1;
+
+    return out;
+}
+
+
+///////////////////////////////////////////////////////////////////////////////
+template <size_t S>
+const std::array<util::vectori<S>,util::pow(2,S)>
+util::noise::basis::type::gradient<S>::CORNERS = generate<S> ();
+
+//-----------------------------------------------------------------------------
+template struct util::noise::basis::type::gradient<1>;
+template struct util::noise::basis::type::gradient<2>;
+template struct util::noise::basis::type::gradient<3>;
+template struct util::noise::basis::type::gradient<4>;

noise/basis/type/gradient.hpp

@@ -0,0 +1,36 @@
+/*
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ *
+ * Copyright 2015 Danny Robson <danny@nerdcruft.net>
+ */
+
+#ifndef __UTIL_NOISE_BASIS_TYPE_GRADIENT_HPP
+#define __UTIL_NOISE_BASIS_TYPE_GRADIENT_HPP
+
+#include "../../../point.hpp"
+#include "../../../maths.hpp"
+
+#include <array>
+
+
+///////////////////////////////////////////////////////////////////////////////
+namespace util { namespace noise { namespace basis { namespace type {
+    template <size_t S>
+    struct gradient {
+
+    protected:
+        static const std::array<vectori<S>,util::pow(2,S)> CORNERS;
+    };
+} } } }
+
+#endif

noise/basis/value.hpp

@@ -18,18 +18,21 @@
 #ifndef __UTIL_NOISE_BASIS_VALUE_HPP
 #define __UTIL_NOISE_BASIS_VALUE_HPP
 
+#include "./type/gradient.hpp"
+
 #include "../fwd.hpp"
 #include "../../range.hpp"
 #include "../../point.hpp"
 
 namespace util { namespace noise { namespace basis {
     /// Single value per grid space
-    template <typename T, lerp_t<T>>
+    template <size_t S, typename T, lerp_t<T>>
-    struct value {
+    struct value : public type::gradient<S> {
         value (seed_t);
 
         range<T> bounds (void) const;
-        T operator() (util::point<2,T>) const;
+
+        T operator() (util::point<S,T>) const;
 
         seed_t seed (void) const;
         seed_t seed (seed_t);

noise/basis/value.ipp

@@ -20,68 +20,72 @@
 #define __UTIL_NOISE_BASIS_VALIE_IPP
 
 #include "../rand.hpp"
+#include "../../types.hpp"
+
 
 namespace util { namespace noise { namespace basis {
     ///////////////////////////////////////////////////////////////////////////
-    template <typename T, util::noise::lerp_t<T> L>
+    template <size_t S, typename T, util::noise::lerp_t<T> L>
-    value<T,L>::value (seed_t _seed):
+    value<S,T,L>::value (seed_t _seed):
         m_seed (_seed)
     { ; }
 
 
     //-------------------------------------------------------------------------
-    template <typename T, util::noise::lerp_t<T> L>
+    template <size_t S, typename T, util::noise::lerp_t<T> L>
     util::range<T>
-    value<T,L>::bounds (void) const
+    value<S,T,L>::bounds (void) const
     {
         return { -1, 1 };
     }
 
 
     //-------------------------------------------------------------------------
-    template <typename T, lerp_t<T> L>
+    template <size_t S, typename T, lerp_t<T> L>
     seed_t
-    value<T,L>::seed (void) const
+    value<S,T,L>::seed (void) const
     {
         return m_seed;
     }
 
 
     //-------------------------------------------------------------------------
-    template <typename T, lerp_t<T> L>
+    template <size_t S, typename T, lerp_t<T> L>
     seed_t
-    value<T,L>::seed (seed_t _seed)
+    value<S,T,L>::seed (seed_t _seed)
     {
         return m_seed = _seed;
     }
 
 
     //-------------------------------------------------------------------------
-    template <typename T, util::noise::lerp_t<T> L>
+    template <size_t S, typename T, util::noise::lerp_t<T> L>
     T
-    value<T,L>::operator() (util::point<2,T> p) const
+    value<S,T,L>::operator() (util::point<S,T> p) const
     {
         // extract integer and fractional parts. be careful to always round down
         auto p_int = floor (p).template cast<intmax_t> ();
         auto p_rem = p - p_int;
 
         // generate the corner points
-        auto p0 = p_int + util::vector<2,intmax_t> { 0, 0 };
+        std::array<pointi<S>,pow(2,S)> p_;
-        auto p1 = p_int + util::vector<2,intmax_t> { 1, 0 };
+        std::transform (std::begin (this->CORNERS), std::end (this->CORNERS),
-        auto p2 = p_int + util::vector<2,intmax_t> { 0, 1 };
+                        std::begin (p_),
-        auto p3 = p_int + util::vector<2,intmax_t> { 1, 1 };
+                        [p_int] (auto i) { return p_int + i; });
 
-        // Generate the four corner values
+        // Generate the corner values
-        T g0 = noise::rand<float> (m_seed, p0);
+        std::array<T,pow(2,S)> g_;
-        T g1 = noise::rand<float> (m_seed, p1);
+        std::transform (std::begin (p_), std::end (p_),
-        T g2 = noise::rand<float> (m_seed, p2);
+                        std::begin (g_),
-        T g3 = noise::rand<float> (m_seed, p3);
+                        [this] (auto i) { return noise::rand<float> (m_seed, i); });
 
         // Interpolate on one dimension, then the other.
-        auto l0 = L (g0, g1, p_rem.x);
+        T l_[pow(2,S)];
-        auto l1 = L (g2, g3, p_rem.x);
+        std::copy (std::begin (g_), std::end (g_), std::begin (l_));
-        auto l_ = L (l0, l1, p_rem.y);
 
-        return l_;
+        for (size_t i = S; i; --i)
+            for (size_t j = 0; j < std::pow(2,i); j += 2)
+                l_[j / 2] = L (l_[j], l_[j+1], p_rem[S-i]);
+        return l_[0];
     }
 } } }

noise/basis/worley.hpp

@@ -18,23 +18,27 @@
 #ifndef __UTIL_NOISE_BASIS_WORLEY_HPP
 #define __UTIL_NOISE_BASIS_WORLEY_HPP
 
+#include "./type/distance.hpp"
+
 #include "../fwd.hpp"
 #include "../../point.hpp"
 #include "../../range.hpp"
 
 namespace util { namespace noise { namespace basis {
-    template <typename T, size_t F = 0>
+    template <size_t S, typename T, size_t F = 0>
-    struct worley {
+    struct worley : public type::distance<S,1> {
         worley (seed_t);
 
         range<T> bounds (void) const;
-        T operator() (util::point<2,T>) const;
+
+        T operator() (util::point<S,T>) const;
 
         seed_t seed (void) const;
         seed_t seed (seed_t);
 
     private:
-        point<2,T> generate (point<2,intmax_t>) const;
+        point<S,T> generate (point<S,intmax_t>) const;
+
         seed_t m_seed;
     };
 } } }

noise/basis/worley.ipp

@@ -25,67 +25,58 @@
 
 namespace util { namespace noise { namespace basis {
     ///////////////////////////////////////////////////////////////////////////
-    template <typename T, size_t F>
+    template <size_t S, typename T, size_t F>
-    worley<T,F>::worley (seed_t _seed):
+    worley<S,T,F>::worley (seed_t _seed):
         m_seed (_seed)
     { ; }
 
 
     //-------------------------------------------------------------------------
-    template <typename T, size_t F>
+    template <size_t S, typename T, size_t F>
     util::range<T>
-    worley<T,F>::bounds (void) const
+    worley<S,T,F>::bounds (void) const
     {
         return { 0.0, 1.5 };
     }
 
 
     //-------------------------------------------------------------------------
-    template <typename T, size_t F>
+    template <size_t S, typename T, size_t F>
     seed_t
-    worley<T,F>::seed (void) const
+    worley<S,T,F>::seed (void) const
     {
         return m_seed;
     }
 
 
     //-------------------------------------------------------------------------
-    template <typename T, size_t F>
+    template <size_t S, typename T, size_t F>
     seed_t
-    worley<T,F>::seed (seed_t _seed)
+    worley<S,T,F>::seed (seed_t _seed)
     {
         return m_seed = _seed;
     }
 
 
     //-------------------------------------------------------------------------
-    template <typename T, size_t F>
+    template <size_t S, typename T, size_t F>
     T
-    worley<T,F>::operator() (util::point<2,T> p) const
+    worley<S,T,F>::operator() (point<S,T> p) const
     {
         // extract integer and fractional parts. be careful to always round down
         auto p_int = floor (p).template cast<intmax_t> ();
         auto p_rem = (p - p_int).template as<point> ();
 
         // setup an array of distances
-        constexpr size_t RADIUS = 1;
+        constexpr size_t COUNT = type::distance<S,1>::OFFSET_SIZE;
-        constexpr size_t COUNT = pow2 (RADIUS * 2 + 1);
+        T distances[COUNT];
-        T distances[COUNT] = { std::numeric_limits<T>::quiet_NaN () };
-        T *cursor = distances;
 
-        // record the distances to each candidate point
+        std::transform (std::begin (this->OFFSETS), std::end (this->OFFSETS),
-        for (signed y_off = -signed(RADIUS); y_off <= signed(RADIUS) ; ++y_off) {
+                        distances,
-            for (signed x_off = -signed(RADIUS); x_off <= signed(RADIUS); ++x_off) {
+                        [p_int,p_rem,this] (auto i) {
-                auto off = vector<2,intmax_t> {x_off, y_off};
+            auto q = this->generate (p_int + i);
-                auto pos = generate (p_int + off);
+            return distance2 (q + i, p_rem);
-
+        });
-                CHECK_LIMIT (pos.x, T{0}, T{1});
-                CHECK_LIMIT (pos.y, T{0}, T{1});
-
-                *cursor = distance2 (pos + off, p_rem);
-                cursor++;
-            }
-        }
 
         // find the f'th lowest value
         static_assert (F < COUNT, "worley order must be less than search radius");
@@ -96,21 +87,10 @@ namespace util { namespace noise { namespace basis {
 
 
     //////////////////////////////////////////////////////////////////////////
-    template <typename T, size_t F>
+    template <size_t S, typename T, size_t F>
-    point<2,T>
+    point<S,T>
-    worley<T,F>::generate (point<2,intmax_t> p) const
+    worley<S,T,F>::generate (point<S,intmax_t> p) const
     {
-        using util::hash::murmur2::mix;
+        return (noise::rand<util::point,T> (m_seed, p) + 1) / 2;
-
-        auto u = mix (m_seed, mix (uint64_t (p.x), uint64_t (p.y)));
-        auto v = mix (u, m_seed);
-
-        auto r = util::point<2,T> {
-            (u & 0xffff) / T{0xffff},
-            (v & 0xffff) / T{0xffff}
-        };
-
-        CHECK_LIMIT (r, T{0}, T{1});
-        return r;
     }
 } } }

noise/fractal/base.hpp

@@ -32,7 +32,11 @@ namespace util { namespace noise { namespace fractal {
     /// lacunarity: per octave frequency scaling factor
     /// amplitude: maximum absolute value of the noise
     /// gain: per octave amplitude scaling factor. typically 1/f.
-    template <typename T, typename B>
+    template <
+        size_t S, // probe point dimensionality
+        typename T, // probe point value_type
+        typename B  // noise basis function
+    >
     struct base {
         using seed_t = uint64_t;
 

noise/fractal/base.ipp

@@ -23,8 +23,8 @@
 
 namespace util { namespace noise { namespace fractal {
     ///////////////////////////////////////////////////////////////////////////
-    template <typename T, typename B>
+    template <size_t S, typename T, typename B>
-    base<T,B>::base (seed_t _seed,
+    base<S,T,B>::base (seed_t _seed,
                        unsigned _octaves,
                        T _H,
                        T _frequency,
@@ -51,27 +51,27 @@ namespace util { namespace noise { namespace fractal {
 
 
     ///////////////////////////////////////////////////////////////////////////
-    template <typename T, typename B>
+    template <size_t S, typename T, typename B>
     unsigned
-    base<T,B>::octaves (unsigned _octaves)
+    base<S,T,B>::octaves (unsigned _octaves)
     {
         return m_octaves = _octaves;
     }
 
 
     //-------------------------------------------------------------------------
-    template <typename T, typename B>
+    template <size_t S, typename T, typename B>
     constexpr unsigned
-    base<T,B>::octaves (void) const
+    base<S,T,B>::octaves (void) const
     {
         return m_octaves;
     }
 
 
     //-------------------------------------------------------------------------
-    template <typename T, typename B>
+    template <size_t S, typename T, typename B>
     T
-    base<T,B>::H (T _h)
+    base<S,T,B>::H (T _h)
     {
         m_H = _h;
         m_invAH = std::pow (m_amplitude, -m_H);
@@ -81,63 +81,63 @@ namespace util { namespace noise { namespace fractal {
 
 
     //-------------------------------------------------------------------------
-    template <typename T, typename B>
+    template <size_t S, typename T, typename B>
     constexpr T
-    base<T,B>::H (void) const
+    base<S,T,B>::H (void) const
     {
         return m_H;
     }
 
 
     //-------------------------------------------------------------------------
-    template <typename T, typename B>
+    template <size_t S, typename T, typename B>
     T
-    base<T,B>::frequency (T _frequency)
+    base<S,T,B>::frequency (T _frequency)
     {
         return m_frequency = _frequency;
     }
 
 
     //-------------------------------------------------------------------------
-    template <typename T, typename B>
+    template <size_t S, typename T, typename B>
     constexpr T
-    base<T,B>::frequency (void) const
+    base<S,T,B>::frequency (void) const
     {
         return m_frequency;
     }
 
 
     //-------------------------------------------------------------------------
-    template <typename T, typename B>
+    template <size_t S, typename T, typename B>
     T
-    base<T,B>::lacunarity (T _lacunarity)
+    base<S,T,B>::lacunarity (T _lacunarity)
     {
         return m_lacunarity = _lacunarity;
     }
 
 
     //-------------------------------------------------------------------------
-    template <typename T, typename B>
+    template <size_t S, typename T, typename B>
     constexpr T
-    base<T,B>::lacunarity (void) const
+    base<S,T,B>::lacunarity (void) const
     {
         return m_lacunarity;
     }
 
 
     //-------------------------------------------------------------------------
-    template <typename T, typename B>
+    template <size_t S, typename T, typename B>
     constexpr T
-    base<T,B>::amplitude (void) const
+    base<S,T,B>::amplitude (void) const
     {
         return m_amplitude;
     }
 
 
     //-------------------------------------------------------------------------
-    template <typename T, typename B>
+    template <size_t S, typename T, typename B>
     T
-    base<T,B>::amplitude (T _amplitude)
+    base<S,T,B>::amplitude (T _amplitude)
     {
         m_amplitude = _amplitude;
         m_invAH = std::pow (m_amplitude, -m_H);
@@ -146,18 +146,18 @@ namespace util { namespace noise { namespace fractal {
 
 
     //-------------------------------------------------------------------------
-    template <typename T, typename B>
+    template <size_t S, typename T, typename B>
     constexpr T
-    base<T,B>::gain (void) const
+    base<S,T,B>::gain (void) const
     {
         return m_gain;
     }
 
 
     //-------------------------------------------------------------------------
-    template <typename T, typename B>
+    template <size_t S, typename T, typename B>
     T
-    base<T,B>::gain (T _gain)
+    base<S,T,B>::gain (T _gain)
     {
         m_gain = _gain;
         m_invGH = std::pow (_gain, m_H);
@@ -166,36 +166,36 @@ namespace util { namespace noise { namespace fractal {
 
 
     //-------------------------------------------------------------------------
-    template <typename T, typename B>
+    template <size_t S, typename T, typename B>
-    typename base<T,B>::seed_t
+    typename base<S,T,B>::seed_t
-    base<T,B>::seed (seed_t _seed)
+    base<S,T,B>::seed (seed_t _seed)
     {
         return m_basis.seed (_seed);
     }
 
 
     //-------------------------------------------------------------------------
-    template <typename T, typename B>
+    template <size_t S, typename T, typename B>
-    typename base<T,B>::seed_t
+    typename base<S,T,B>::seed_t
-    base<T,B>::seed (void) const
+    base<S,T,B>::seed (void) const
     {
         return m_basis.seed ();
     }
 
 
     //-------------------------------------------------------------------------
-    template <typename T, typename B>
+    template <size_t S, typename T, typename B>
     const B&
-    base<T,B>::basis (void) const
+    base<S,T,B>::basis (void) const
     {
         return m_basis;
     }
 
 
     //-------------------------------------------------------------------------
-    template <typename T, typename B>
+    template <size_t S, typename T, typename B>
     B&
-    base<T,B>::basis (void)
+    base<S,T,B>::basis (void)
     {
         return m_basis;
     }

noise/fractal/fbm.hpp

@@ -33,9 +33,13 @@ namespace util { namespace noise { namespace fractal {
     /// lacunarity: per octave frequency scaling factor
     /// amplitude: maximum absolute value of the noise
     /// gain: per octave amplitude scaling factor. typically 1/f.
-    template <typename T, typename B>
+    template <
-    struct fbm : public base<T,B> {
+        size_t S,   // probe point dimensionality
-        using seed_t = typename base<T,B>::seed_t;
+        typename T, // probe point value_type
+        typename B  // generating basis function
+    >
+    struct fbm : public base<S,T,B> {
+        using seed_t = typename base<S,T,B>::seed_t;
 
         static constexpr unsigned DEFAULT_OCTAVES = 8;
         static constexpr T DEFAULT_H          = 1;
@@ -53,7 +57,7 @@ namespace util { namespace noise { namespace fractal {
              T gain);
         fbm (seed_t);
 
-        T operator() (util::point<2,T>) const;
+        T operator() (point<S,T>) const;
     };
 } } }
 

noise/fractal/fbm.ipp

@@ -23,15 +23,15 @@
 
 namespace util { namespace noise { namespace fractal {
     ///////////////////////////////////////////////////////////////////////////
-    template <typename T, typename B>
+    template <size_t S, typename T, typename B>
-    fbm<T,B>::fbm (seed_t _seed,
+    fbm<S,T,B>::fbm (seed_t _seed,
                      unsigned _octaves,
                      T        _H,
                      T        _frequency,
                      T        _lacunarity,
                      T        _amplitude,
                      T        _gain):
-        base<T,B> (_seed,
+        base<S,T,B> (_seed,
                      _octaves,
                      _H,
                      _frequency,
@@ -42,9 +42,9 @@ namespace util { namespace noise { namespace fractal {
 
 
     //-------------------------------------------------------------------------
-    template <typename T, typename B>
+    template <size_t S, typename T, typename B>
-    fbm<T,B>::fbm (seed_t _seed):
+    fbm<S,T,B>::fbm (seed_t _seed):
-        fbm<T,B> (_seed,
+        fbm<S,T,B> (_seed,
                     DEFAULT_OCTAVES,
                     DEFAULT_H,
                     DEFAULT_FREQUENCY,
@@ -55,9 +55,9 @@ namespace util { namespace noise { namespace fractal {
 
 
     ///////////////////////////////////////////////////////////////////////////
-    template <typename T, typename B>
+    template <size_t S, typename T, typename B>
     T
-    fbm<T,B>::operator() (util::point<2,T> p) const
+    fbm<S,T,B>::operator() (point<S,T> p) const
     {
         T total = 0;
         T scale = this->m_invAH;
@@ -67,7 +67,7 @@ namespace util { namespace noise { namespace fractal {
         for (size_t i = 0; i < this->m_octaves; ++i) {
             total += this->m_basis (p) * scale;
 
-            p += T{1};
+            p += PI<float>;
             p *= this->m_lacunarity;
             scale *= this->m_invGH;
         }

noise/fractal/hetero.hpp

@@ -25,9 +25,13 @@
 namespace util { namespace noise { namespace fractal {
     ///////////////////////////////////////////////////////////////////////
     /// Heterogeneous procedural terrain fucntion: stats by altitude method
-    template <typename T, typename B>
+    template <
-    struct hetero : public base<T,B> {
+        size_t S,   // probe point dimensionality
-        using seed_t = typename base<T,B>::seed_t;
+        typename T, // probe point value_type
+        typename B  // generating basis function
+    >
+    struct hetero : public base<S,T,B> {
+        using seed_t = typename base<S,T,B>::seed_t;
 
         static constexpr unsigned DEFAULT_OCTAVES = 6;
         static constexpr T DEFAULT_H = T(0.75);
@@ -59,7 +63,7 @@ namespace util { namespace noise { namespace fractal {
         constexpr T offset (void) const;
         T offset (T);
 
-        T operator() (util::point<2,T>) const;
+        T operator() (point<S,T>) const;
 
     private:
         T m_offset;

noise/fractal/hetero.ipp

@@ -21,15 +21,15 @@
 
 namespace util { namespace noise { namespace fractal {
     ///////////////////////////////////////////////////////////////////////////
-    template <typename T, typename B>
+    template <size_t S, typename T, typename B>
-    hetero<T,B>::hetero(seed_t _seed,
+    hetero<S,T,B>::hetero(seed_t _seed,
                           unsigned _octaves,
                           T _H,
                           T _frequency,
                           T _lacunarity,
                           T _amplitude,
                           T _gain):
-        hetero<T,B> (_seed,
+        hetero<S,T,B> (_seed,
                        _octaves,
                        _H,
                        _frequency,
@@ -41,8 +41,8 @@ namespace util { namespace noise { namespace fractal {
 
 
     //-------------------------------------------------------------------------
-    template <typename T, typename B>
+    template <size_t S, typename T, typename B>
-    hetero<T,B>::hetero(seed_t _seed,
+    hetero<S,T,B>::hetero(seed_t _seed,
                           unsigned _octaves,
                           T _H,
                           T _frequency,
@@ -50,7 +50,7 @@ namespace util { namespace noise { namespace fractal {
                           T _amplitude,
                           T _gain,
                           T _offset):
-        base<T,B> (_seed,
+        base<S,T,B> (_seed,
                      _octaves,
                      _H,
                      _frequency,
@@ -62,9 +62,9 @@ namespace util { namespace noise { namespace fractal {
 
 
     //-------------------------------------------------------------------------
-    template <typename T, typename B>
+    template <size_t S, typename T, typename B>
-    hetero<T,B>::hetero (seed_t _seed):
+    hetero<S,T,B>::hetero (seed_t _seed):
-        hetero<T,B> (_seed,
+        hetero<S,T,B> (_seed,
                        DEFAULT_OCTAVES,
                        DEFAULT_H,
                        DEFAULT_FREQUENCY,
@@ -76,27 +76,27 @@ namespace util { namespace noise { namespace fractal {
 
 
     ///////////////////////////////////////////////////////////////////////////
-    template <typename T, typename B>
+    template <size_t S, typename T, typename B>
     constexpr T
-    hetero<T,B>::offset (void) const
+    hetero<S,T,B>::offset (void) const
     {
         return m_offset;
     }
 
 
     //-------------------------------------------------------------------------
-    template <typename T, typename B>
+    template <size_t S, typename T, typename B>
     T
-    hetero<T,B>::offset (T _offset)
+    hetero<S,T,B>::offset (T _offset)
     {
         return m_offset = _offset;
     }
 
 
     ///////////////////////////////////////////////////////////////////////////
-    template <typename T, typename B>
+    template <size_t S, typename T, typename B>
     T
-    hetero<T,B>::operator() (util::point<2,T> p) const
+    hetero<S,T,B>::operator() (point<S,T> p) const
     {
         T scale = this->m_invAH;
 
@@ -115,7 +115,7 @@ namespace util { namespace noise { namespace fractal {
 
             result += increment;
 
-            p += T{1};
+            p += PI<T>;
             p *= this->m_lacunarity;
         }
 

noise/fractal/hmf.hpp

@@ -24,9 +24,13 @@
 namespace util { namespace noise { namespace fractal {
     ///////////////////////////////////////////////////////////////////////
     /// Musgrave's "Hybrid MultiFractal"
-    template <typename T, typename B>
+    template <
-    struct hmf : public base<T,B> {
+        size_t S,   // probe point dimensionality
-        using seed_t = typename base<T,B>::seed_t;
+        typename T, // probe point value_type
+        typename B  // generating basis function
+    >
+    struct hmf : public base<S,T,B> {
+        using seed_t = typename base<S,T,B>::seed_t;
 
         // H should be fairly low due to the decreasing weight parameter in eval
         static constexpr unsigned DEFAULT_OCTAVES = 6;
@@ -48,7 +52,7 @@ namespace util { namespace noise { namespace fractal {
 
         hmf (seed_t);
 
-        T operator() (point<2,T>) const;
+        T operator() (point<S,T>) const;
 
     private:
         T m_offset;

noise/fractal/hmf.ipp

@@ -22,8 +22,8 @@
 
 namespace util { namespace noise { namespace fractal {
     ///////////////////////////////////////////////////////////////////////////
-    template <typename T, typename B>
+    template <size_t S, typename T, typename B>
-    hmf<T,B>::hmf (seed_t _seed,
+    hmf<S,T,B>::hmf (seed_t _seed,
                      unsigned _octaves,
                      T _H,
                      T _frequency,
@@ -31,7 +31,7 @@ namespace util { namespace noise { namespace fractal {
                      T _amplitude,
                      T _gain,
                      T _offset):
-        base<T,B> (_seed,
+        base<S,T,B> (_seed,
                      _octaves,
                      _H,
                      _frequency,
@@ -43,9 +43,9 @@ namespace util { namespace noise { namespace fractal {
 
 
     //-------------------------------------------------------------------------
-    template <typename T, typename B>
+    template <size_t S, typename T, typename B>
-    hmf<T,B>::hmf (seed_t _seed):
+    hmf<S,T,B>::hmf (seed_t _seed):
-        hmf<T,B> (_seed,
+        hmf<S,T,B> (_seed,
                     DEFAULT_OCTAVES,
                     DEFAULT_H,
                     DEFAULT_FREQUENCY,
@@ -57,9 +57,9 @@ namespace util { namespace noise { namespace fractal {
 
 
     ///////////////////////////////////////////////////////////////////////////
-    template <typename T, typename B>
+    template <size_t S, typename T, typename B>
     T
-    hmf<T,B>::operator() (util::point<2,T> p) const
+    hmf<S,T,B>::operator() (point<S,T> p) const
     {
         T scale = this->m_invAH;
 
@@ -77,7 +77,7 @@ namespace util { namespace noise { namespace fractal {
             weight  = min (weight, T{1});
 
             scale *= this->m_invGH;
-            p += T{1};
+            p += PI<T>;
             p *= this->m_lacunarity;
         }
 

noise/fractal/rmf.hpp

@@ -32,9 +32,13 @@ namespace util { namespace noise { namespace fractal {
     /// lacunarity: incremental octave frequency scaling factor
     /// amplitude: value scaling factor for the base octave
     /// gain: incremental octave value scaling factor
-    template <typename T, typename B>
+    template <
-    struct rmf : public base<T,B> {
+        size_t S,   // probe point dimensionality
-        using seed_t = typename base<T,B>::seed_t;
+        typename T, // probe point value_type
+        typename B  // generating basis function
+    >
+    struct rmf : public base<S,T,B> {
+        using seed_t = typename base<S,T,B>::seed_t;
 
         static constexpr unsigned DEFAULT_OCTAVES = 5;
         static constexpr T DEFAULT_H = 1;
@@ -55,7 +59,7 @@ namespace util { namespace noise { namespace fractal {
 
         rmf (seed_t);
 
-        T operator() (util::point<2,T>) const;
+        T operator() (point<S,T>) const;
 
         T offset (void) const;
         T offset (T);

noise/fractal/rmf.ipp

@@ -22,8 +22,8 @@
 
 namespace util { namespace noise { namespace fractal {
     ///////////////////////////////////////////////////////////////////////////
-    template <typename T, typename B>
+    template <size_t S, typename T, typename B>
-    rmf<T,B>::rmf (seed_t _seed,
+    rmf<S,T,B>::rmf (seed_t _seed,
                      unsigned _octaves,
                      T        _H,
                      T        _offset,
@@ -31,7 +31,7 @@ namespace util { namespace noise { namespace fractal {
                      T        _lacunarity,
                      T        _amplitude,
                      T        _gain):
-        base<T,B> (_seed,
+        base<S,T,B> (_seed,
                      _octaves,
                      _H,
                      _frequency,
@@ -43,9 +43,9 @@ namespace util { namespace noise { namespace fractal {
 
 
     //-------------------------------------------------------------------------
-    template <typename T, typename B>
+    template <size_t S, typename T, typename B>
-    rmf<T,B>::rmf (seed_t _seed):
+    rmf<S,T,B>::rmf (seed_t _seed):
-        rmf<T,B> (_seed,
+        rmf<S,T,B> (_seed,
                     DEFAULT_OCTAVES,
                     DEFAULT_H,
                     DEFAULT_FREQUENCY,
@@ -59,9 +59,9 @@ namespace util { namespace noise { namespace fractal {
     ///////////////////////////////////////////////////////////////////////////
     // we use the name 'amplitude' instead of musgrave's 'gain'.
     // assumes basis distribution [-1,1] and offset ~= 1
-    template <typename T, typename B>
+    template <size_t S, typename T, typename B>
     T
-    rmf<T,B>::operator() (util::point<2,T> p) const
+    rmf<S,T,B>::operator() (point<S,T> p) const
     {
         T scale = this->m_invAH;
 
@@ -92,7 +92,7 @@ namespace util { namespace noise { namespace fractal {
 
             scale *= this->m_invGH;
 
-            p += T{1};
+            p += PI<T>;
             p *= this->m_lacunarity;
         }
 
@@ -101,18 +101,18 @@ namespace util { namespace noise { namespace fractal {
 
 
     ///////////////////////////////////////////////////////////////////////////
-    template <typename T, typename B>
+    template <size_t S, typename T, typename B>
     T
-    rmf<T,B>::offset (void) const
+    rmf<S,T,B>::offset (void) const
     {
         return m_offset;
     }
 
 
     //-------------------------------------------------------------------------
-    template <typename T, typename B>
+    template <size_t S, typename T, typename B>
     T
-    rmf<T,B>::offset (T _offset)
+    rmf<S,T,B>::offset (T _offset)
     {
         return m_offset = _offset;
     }

noise/fractal/runtime.cpp

@@ -17,4 +17,7 @@
 #include "runtime.hpp"
 #include "../basis/runtime.hpp"
 
-template struct util::noise::fractal::runtime<float,util::noise::basis::runtime<float>>;
+template struct util::noise::fractal::runtime<
+    2,float,
+    util::noise::basis::runtime<2,float>
+>;

noise/fractal/runtime.hpp

@@ -25,7 +25,11 @@
 #include <memory>
 
 namespace util { namespace noise { namespace fractal {
-    template <typename T, typename B>
+    template <
+        size_t S,   // probe point dimensionality
+        typename T, // probe point value_type
+        typename B  // generating basis function
+    >
     struct runtime {
     public:
         using seed_t = uint64_t;
@@ -37,7 +41,7 @@ namespace util { namespace noise { namespace fractal {
         runtime& operator= (const runtime&) = delete;
 
         // basis functions
-        T operator () (util::point<2,T> p) const { return (*m_child) (p); }
+        T operator () (point<S,T> p) const { return (*m_child) (p); }
 
         unsigned octaves (void) const { return m_child->octaves (); }
         unsigned octaves (unsigned _octaves) { return m_child->octaves (_octaves); }
@@ -67,7 +71,7 @@ namespace util { namespace noise { namespace fractal {
         struct base {
             virtual ~base () = default;
 
-            virtual T operator() (util::point<2,T>) = 0;
+            virtual T operator() (util::point<S,T>) = 0;
 
             virtual unsigned octaves (void) const = 0;
             virtual unsigned octaves (unsigned) = 0;
@@ -116,7 +120,7 @@ namespace util { namespace noise { namespace fractal {
                      _gain)
             { ; }
  
-            T operator() (util::point<2,T> p) override { return data (p); }
+            T operator() (util::point<S,T> p) override { return data (p); }
 
             unsigned octaves (void) const override { return data.octaves (); }
             unsigned octaves (unsigned _octaves) override { return data.octaves (_octaves); }

noise/rand.hpp

@@ -19,13 +19,24 @@
 
 namespace util { namespace noise {
     /// generate a uniform random floating point in the range [-1, 1] from a seed and vector
-    template <typename T, typename V>
+    template <
-    T
+        typename U,
-    rand (uint64_t seed, V value);
+        size_t S,
+        typename T,
+        template <size_t,typename> class Q
+    >
+    U
+    rand (uint64_t seed, Q<S,T> value);
 
-    template <size_t N, typename T, typename V>
+    template <
-    vector<N,T>
+        template <size_t,typename> class R,
-    rand (uint64_t seed, V value);
+        typename U,
+        size_t S,
+        typename T,
+        template <size_t,typename> class Q
+    >
+    R<S,U>
+    rand (uint64_t seed, Q<S,T> value);
 } }
 
 #include "rand.ipp"

noise/rand.ipp

@@ -24,37 +24,53 @@
 
 
 namespace util { namespace noise {
-    template <typename T, typename V>
+    //-------------------------------------------------------------------------
-    T
+    template <
-    rand (uint64_t seed, V value)
+        typename U,
+        size_t S,
+        typename T,
+        template <size_t,typename> class Q
+    >
+    U
+    rand (uint64_t seed, Q<S,T> query)
     {
+        static_assert (std::is_integral<T>::value,
+                       "mixing only works on integral types");
+
         uint64_t accum = seed;
-        for (auto i: value)
+        for (auto i: query)
             accum = hash::murmur2::mix (accum, i);
 
-        T out = (accum & 0xFFFF) / T{0xFFFF};
+        U result = (accum & 0xFFFF) / U{0xFFFF};
 
-        out *= 2;
+        result *= 2;
-        out -= 1;
+        result -= 1;
 
-        return out;
+        return result;
     }
 
 
-    template <size_t N, typename T, typename V>
+    //-------------------------------------------------------------------------
-    vector<N,T>
+    template <
-    rand (uint64_t seed, V value)
+        template <size_t,typename> class R,
+        typename U,
+        size_t S,
+        typename T,
+        template <size_t,typename> class Q
+    >
+    R<S,U>
+    rand (uint64_t seed, Q<S,T> query)
     {
         uint64_t accum = seed;
-        for (auto i: value)
+        for (auto i: query)
             accum = hash::murmur2::mix (accum, i);
 
-        vector<N,T> out;
+        R<S,U> result;
-        for (auto &i: out) {
+        for (auto &i: result) {
-            i = (accum & 0xFFFF) / T{0xFFFF};
+            i = (accum & 0xFFFF) / U{0xFFFF};
             accum = hash::murmur2::mix (accum, seed);
         }
 
-        return out * 2 - 1;
+        return result * 2 - 1;
     }
 } }

noise/turbulence.hpp

@@ -29,20 +29,36 @@ namespace util { namespace noise {
     /// assumes the pertubation function is roughly symetrical around 0.
     /// nothing will explode if it isn't, but you'll see strong directional
     /// artefacts with higher scaling factors.
-    template <typename T, typename D, typename P>
+    template <
+        size_t   S, // dimension
+        typename T, // value_type
+        typename D, // data fractal
+        typename P  // pertubation fractal
+    >
     struct turbulence {
+        static constexpr auto dimension = S;
+        using value_type = T;
+
         using seed_t = uint64_t;
-        turbulence (seed_t, vector<2,T> scale);
+
+        turbulence (seed_t, vector<S,T> scale);
 
         seed_t seed (seed_t);
         seed_t seed (void) const;
 
-        constexpr T operator() (point<2,T>) const;
+        constexpr T operator() (point<S,T>) const;
 
         D data;
-        P perturb[2];
 
-        vector<2,T> scale;
+        // XXX: use a union to defer initialization of pertubation fractals in
+        // the constructor. i know this is horrible, but there's no time to
+        // write the proper generator constructor to pass out the seeds.
+        union {
+            char _[0];
+            P perturb[S];
+        };
+
+        vector<S,T> scale;
     };
 } }
 

noise/turbulence.ipp

@@ -24,50 +24,48 @@
 
 namespace util { namespace noise {
     ///////////////////////////////////////////////////////////////////////////
-    template <typename T, typename D, typename P>
+    template <size_t S, typename T, typename D, typename P>
-    turbulence<T,D,P>::turbulence (seed_t _seed,
+    turbulence<S,T,D,P>::turbulence (seed_t _seed,
-                                   vector<2,T> _scale):
+                                     vector<S,T> _scale):
         data (_seed),
-        perturb {
-            hash::wang (_seed),
-            hash::wang (hash::wang (_seed))
-        },
         scale (_scale)
-    { ; }
+    {
+        for (auto &p: perturb)
+            new (&p) P (_seed = hash::wang (_seed));
+    }
 
     ////////////////////////////////////////////////////////////////////////////
-    template <typename T, typename D, typename P>
+    template <size_t S, typename T, typename D, typename P>
-    typename turbulence<T,D,P>::seed_t
+    typename turbulence<S,T,D,P>::seed_t
-    turbulence<T,D,P>::seed (void) const
+    turbulence<S,T,D,P>::seed (void) const
     {
         return data.seed ();
     }
 
 
     //-------------------------------------------------------------------------
-    template <typename T, typename D, typename P>
+    template <size_t S, typename T, typename D, typename P>
-    typename turbulence<T,D,P>::seed_t
+    typename turbulence<S,T,D,P>::seed_t
-    turbulence<T,D,P>::seed (seed_t _seed)
+    turbulence<S,T,D,P>::seed (seed_t _seed)
     {
         auto ret = _seed;
+        data.seed (_seed);
 
-        data.seed (_seed);          _seed = hash::wang (_seed);
+        for (size_t i = 0; i < S; ++i)
-        perturb[0].seed (_seed);    _seed = hash::wang (_seed);
+            perturb[i].seed (_seed = hash::wang (_seed));
-        perturb[1].seed (_seed);    _seed = hash::wang (_seed);
 
         return ret;
     }
 
 
     ///////////////////////////////////////////////////////////////////////////
-    template <typename T, typename D, typename P>
+    template <size_t S, typename T, typename D, typename P>
     constexpr T
-    turbulence<T,D,P>::operator() (point<2,T> p) const
+    turbulence<S,T,D,P>::operator() (point<S,T> p) const
     {
-        vector<2,T> n = {
+        vector<S,T> n;
-            perturb[0] (p),
+        for (size_t i = 0; i < S; ++i)
-            perturb[1] (p)
+            n[i] = perturb[i] (p);
-        };
 
         // scale by the data frequency so that we match scale
         return data (p + n * scale / data.frequency ());

tools/noise.cpp

@@ -25,11 +25,13 @@
 #include "region.hpp"
 
 
+constexpr size_t S = 3;
+
 ///////////////////////////////////////////////////////////////////////////////
-template struct util::noise::fractal::fbm<float, util::noise::basis::perlin<float,util::lerp::cubic>>;
+template struct util::noise::fractal::fbm<S,float, util::noise::basis::perlin<S,float,util::lerp::cubic>>;
-template struct util::noise::fractal::hmf<float, util::noise::basis::value<float,util::lerp::cubic>>;
+template struct util::noise::fractal::hmf<S,float, util::noise::basis::value<S,float,util::lerp::cubic>>;
-template struct util::noise::fractal::rmf<float, util::noise::basis::constant<float>>;
+template struct util::noise::fractal::rmf<S,float, util::noise::basis::constant<S,float>>;
-template struct util::noise::fractal::hetero<float, util::noise::basis::worley<float,2>>;
+template struct util::noise::fractal::hetero<S,float, util::noise::basis::worley<S,float,S>>;
 
 
 ///////////////////////////////////////////////////////////////////////////////
@@ -221,30 +223,30 @@ main (int argc, char **argv)
 #endif
 
     util::noise::turbulence<
-        float,
+        S,float,
         util::noise::fractal::runtime<
-            float,
+            S,float,
-            util::noise::basis::runtime<float>
+            util::noise::basis::runtime<S,float>
         >,
         util::noise::fractal::fbm<
-            float,
+            S,float,
             util::noise::basis::perlin<
-                float,
+                S,float,
                 util::lerp::cubic
             >
         >
-    > t (seed, { turbulence, turbulence });
+    > t (seed, util::vectorf<S> (turbulence));
 
     auto &f = t.data;
 
     switch (fractal) {
         using namespace util::noise;
 
-        case FBM:    f.reset<fractal::fbm   <float,basis::runtime<float>>> (seed); break;
+        case FBM:    f.reset<fractal::fbm<S,float,basis::runtime<S,float>>> (seed); break;
-        case HMF:    f.reset<fractal::hmf   <float,basis::runtime<float>>> (seed); break;
+        case HMF:    f.reset<fractal::hmf<S,float,basis::runtime<S,float>>> (seed); break;
-        case RMF:    f.reset<fractal::rmf   <float,basis::runtime<float>>> (seed); break;
+        case RMF:    f.reset<fractal::rmf<S,float,basis::runtime<S,float>>> (seed); break;
         case HETERO: {
-            auto &child = f.reset<fractal::hetero<float,basis::runtime<float>>> (seed);
+            auto &child = f.reset<fractal::hetero<S,float,basis::runtime<S,float>>> (seed);
             if (!std::isnan (offset))
                 child.offset (offset);
             break;
@@ -260,11 +262,11 @@ main (int argc, char **argv)
 
         case PERLIN: {
             switch (lerp) {
-                case LINEAR:    b.reset<basis::perlin<float,util::lerp::linear>> (seed); break;
+                case LINEAR:    b.reset<basis::perlin<S,float,util::lerp::linear>> (seed); break;
-                case CUBIC:     b.reset<basis::perlin<float,util::lerp::cubic>> (seed);    break;
+                case CUBIC:     b.reset<basis::perlin<S,float,util::lerp::cubic>> (seed);    break;
-                case QUINTIC:   b.reset<basis::perlin<float,util::lerp::quintic>> (seed);  break;
+                case QUINTIC:   b.reset<basis::perlin<S,float,util::lerp::quintic>> (seed);  break;
-                case COSINE:    b.reset<basis::perlin<float,util::lerp::cosine>> (seed);   break;
+                case COSINE:    b.reset<basis::perlin<S,float,util::lerp::cosine>> (seed);   break;
-                case TRUNC:     b.reset<basis::perlin<float,util::lerp::trunc>> (seed);    break;
+                case TRUNC:     b.reset<basis::perlin<S,float,util::lerp::trunc>> (seed);    break;
 
             default:
                 unreachable ();
@@ -274,11 +276,11 @@ main (int argc, char **argv)
 
         case EXPDIST: {
             switch (lerp) {
-                case LINEAR:    b.reset<basis::perlin<float,util::lerp::linear,basis::expgrad>> (seed); break;
+                case LINEAR:    b.reset<basis::perlin<S,float,util::lerp::linear,basis::expgrad>> (seed); break;
-                case CUBIC:     b.reset<basis::perlin<float,util::lerp::cubic,basis::expgrad>> (seed);    break;
+                case CUBIC:     b.reset<basis::perlin<S,float,util::lerp::cubic,basis::expgrad>> (seed);    break;
-                case QUINTIC:   b.reset<basis::perlin<float,util::lerp::quintic,basis::expgrad>> (seed);  break;
+                case QUINTIC:   b.reset<basis::perlin<S,float,util::lerp::quintic,basis::expgrad>> (seed);  break;
-                case COSINE:    b.reset<basis::perlin<float,util::lerp::cosine,basis::expgrad>> (seed);   break;
+                case COSINE:    b.reset<basis::perlin<S,float,util::lerp::cosine,basis::expgrad>> (seed);   break;
-                case TRUNC:     b.reset<basis::perlin<float,util::lerp::trunc,basis::expgrad>> (seed);    break;
+                case TRUNC:     b.reset<basis::perlin<S,float,util::lerp::trunc,basis::expgrad>> (seed);    break;
 
             default:
                 unreachable ();
@@ -288,11 +290,11 @@ main (int argc, char **argv)
 
         case VALUE: {
             switch (lerp) {
-                case LINEAR:    b.reset<basis::value<float,util::lerp::linear>> (seed); break;
+                case LINEAR:    b.reset<basis::value<S,float,util::lerp::linear>> (seed); break;
-                case CUBIC:     b.reset<basis::value<float,util::lerp::cubic>> (seed);    break;
+                case CUBIC:     b.reset<basis::value<S,float,util::lerp::cubic>> (seed);    break;
-                case QUINTIC:   b.reset<basis::value<float,util::lerp::quintic>> (seed);  break;
+                case QUINTIC:   b.reset<basis::value<S,float,util::lerp::quintic>> (seed);  break;
-                case COSINE:    b.reset<basis::value<float,util::lerp::cosine>> (seed);   break;
+                case COSINE:    b.reset<basis::value<S,float,util::lerp::cosine>> (seed);   break;
-                case TRUNC:     b.reset<basis::value<float,util::lerp::trunc>> (seed);    break;
+                case TRUNC:     b.reset<basis::value<S,float,util::lerp::trunc>> (seed);    break;
 
             default:
                 unreachable ();
@@ -301,12 +303,12 @@ main (int argc, char **argv)
         }
 
         case WORLEY: {
-            b.reset<util::noise::basis::worley<float>> (seed);
+            b.reset<util::noise::basis::worley<S,float>> (seed);
             break;
         }
 
         case PATCH: {
-            b.reset<util::noise::basis::patch<float>> (seed, width);
+            b.reset<util::noise::basis::patch<S,float>> (seed, width);
             break;
         }
 
@@ -322,9 +324,8 @@ main (int argc, char **argv)
     if (!std::isnan (amplitude)) f.amplitude (amplitude);
     if (!std::isnan (gain)) f.gain (gain);
 
-
+    for (auto &p: t.perturb)
-    t.perturb[0].frequency ( scale / res.w);
+        p.frequency (scale / res.w);
-    t.perturb[1].frequency ( scale / res.w);
 
     util::image::buffer<float> img (res);
 
@@ -336,8 +337,10 @@ main (int argc, char **argv)
 
     {
         for (size_t y = 0; y < res.h; ++y)
-            for (size_t x = 0; x < res.w; ++x)
+            for (size_t x = 0; x < res.w; ++x) {
-                img[{x, y}] = t (util::point2f {float (x), float (y)} + OFFSET);
+                util::point2f p (x, y);
+                img[{x, y}] = t ((p + OFFSET).redim<S> ());
+            }
     }
 
     // working on the assumption that all octave images are based on summation,
@@ -350,8 +353,10 @@ main (int argc, char **argv)
         auto prev = img.clone ();
 
         for (size_t y = 0; y < res.h; ++y)
-            for (size_t x = 0; x < res.w; ++x)
+            for (size_t x = 0; x < res.w; ++x) {
-                prev[{x,y}] = t (util::point2f {float (x), float (y)} + OFFSET);
+                util::point2f p (x, y);
+                prev[{x,y}] = t ((p + OFFSET).redim<S> ());
+            }
 
         CHECK_EQ (img.stride (), prev.stride ());
         for (size_t i = 0; i < img.size (); ++i)