simplex: replace the N-dimensional gradient array with a simple seed.

demos/test/simplex.cpp

@@ -21,51 +21,37 @@ int main(int argc, char **argv)
 {
     UNUSED(argc, argv);
 
-    int const period = 10;
-    float const zoom = 0.02f;
-
-    /* Create some gradients */
-    array2d<vec2> gradients(vec_t<ptrdiff_t, 2>({period, period}));
-    for (int i = 0 ; i < period ; ++i)
-        for (int j = 0 ; j < period ; ++j)
-            gradients[i][j] = normalize(vec2(rand(-1.f, 1.f), rand(-1.f, 1.f)));
-    simplex_interpolator<2> s;
-    s.SetGradients(gradients);
+    float const zoom = 0.03f;
+    int const octaves = 10;
 
     /* Create an image */
-    ivec2 const size(800, 600);
+    ivec2 const size(1280, 720);
     Image img(size);
     array2d<vec4> &data = img.Lock2D<PixelFormat::RGBA_F32>();
 
     /* Fill image with simplex noise */
+    simplex_interpolator<2> s2;
+    simplex_interpolator<3> s3;
+    simplex_interpolator<4> s4;
+
     for (int j = 0; j < size.y; ++j)
     for (int i = 0; i < size.x; ++i)
     {
-        float p = 0.5f * s.Interp(zoom * vec2(i, j));
-#if 1
-        for (int k = 2; k < 128; k *= 2)
-            p += 0.5f / k * s.Interp(zoom * k * vec2(i, j));
-            //p += 0.5f / k * lol::abs(s.Interp(zoom * k * vec2(i, j)));
-#endif
-        data[i][j] = vec4(saturate(0.5f + p), 0.f, 0.f, 1.f);
-    }
+        float sum = 0.f;
+        int maxoctave = (j < size.y / 2) ? 1 : (1 << octaves);
 
-    /* Mark simplex vertices */
-    vec2 diagonal = normalize(vec2(1.f));
-    vec2 dx = mat2::rotate(60.f) * diagonal;
-    vec2 dy = mat2::rotate(-60.f) * diagonal;
-    for (int j = -100; j < 100; ++j)
-    for (int i = -100; i < 100; ++i)
-    {
-        ivec2 coord = ivec2(i / zoom * dx + j / zoom * dy);
-        if (coord.x >= 0 && coord.x < size.x - 1
-             && coord.y >= 0 && coord.y < size.y - 1)
+        for (int k = 1; k <= maxoctave; k *= 2)
         {
-            data[coord.x][coord.y] = vec4(1.f, 0.f, 1.f, 1.f);
-            data[coord.x + 1][coord.y] = vec4(1.f, 0.f, 1.f, 1.f);
-            data[coord.x][coord.y + 1] = vec4(1.f, 0.f, 1.f, 1.f);
-            data[coord.x + 1][coord.y + 1] = vec4(1.f, 0.f, 1.f, 1.f);
+            if (i < size.x / 3)
+                sum += 0.5f / k * s2.Interp(zoom * k * vec2(i, j));
+            else if (i < size.x * 2 / 3)
+                sum += 0.5f / k * s3.Interp(zoom * k * vec3(i, j, 0.0f));
+            else
+                sum += 0.5f / k * s4.Interp(zoom * k * vec4(i, j, 0.0f, 0.0f));
         }
+
+        float c = saturate(0.5f + sum);
+        data[i][j] = vec4(c, c, c, 1.f);
     }
 
     /* Save image */

src/lol/math/simplex_interpolator.h

@@ -12,6 +12,8 @@
 
 #pragma once
 
+#include <functional>
+
 namespace lol
 {
 
@@ -19,21 +21,37 @@ template<int N>
 class simplex_interpolator
 {
 public:
-
-    simplex_interpolator() :
-        m_gradients()
+    simplex_interpolator(int seed = 0)
+      : m_seed(seed)
     {
-        this->InitBase();
-    }
+        /* Matrix coordinate transformation to skew simplex referential is done
+         * by inverting the base matrix M which is written as follows:
+         *
+         *  M = | a b b b … |        M^(-1) = | c d d d … |
+         *      | b a b b … |                 | d c d d … |
+         *      | b b a b … |                 | d d c d … |
+         *      | …         |                 | …         |
+         *
+         *  where a, b, c, d are computed below ↴
+         */
 
-    inline void SetGradients(arraynd<N, vec_t<float, N> > const & gradients)
-    {
-        this->m_gradients = gradients;
-    }
+        float b = (1.f - lol::sqrt(N + 1.f)) / lol::sqrt((float)N * N * N);
+        float a = b + lol::sqrt((N + 1.f) / N);
 
-    inline arraynd<N, vec_t<float, N> > const & GetGradients() const
-    {
-        return this->m_gradients;
+        // determinant of matrix M
+        float det = a * (a + (N - 2) * b) - (b * b) * (N - 1);
+
+        float c = (a + (N - 2) * b) / det;
+        float d = -b / det;
+
+        for (int i = 0; i < N; ++i)
+        {
+            for (int j = 0; j < N; ++j)
+            {
+                m_base[i][j] = (i == j ? a : b);
+                m_base_inverse[i][j] = (i == j ? c : d);
+            }
+        }
     }
 
     // Single interpolation
@@ -57,20 +75,19 @@ public:
     }
 
 protected:
-
     inline float LastInterp(vec_t<int, N> origin,
                             vec_t<float, N> const & pos,
                             vec_t<int, N> const & index_order) const
     {
-        // “corner” will traverse the simplex along its edges in
-        // orthonormal coordinates
+        /* “corner” will traverse the simplex along its edges in
+         * orthonormal coordinates. */
         vec_t<float, N> corner(0);
         float result = 0;
 
         for (int i = 0; i < N + 1; ++i)
         {
             vec_t<float, N> const delta = pos - corner;
-            vec_t<float, N> const &gradient = this->m_gradients[origin];
+            vec_t<float, N> const &gradient = GetGradient(origin);
 
             /* We use 4/3 because the maximum radius of influence for
              * a given corner is sqrt(3/4). FIXME: check whether this
@@ -85,8 +102,7 @@ protected:
             if (i < N)
             {
                 corner += m_base[index_order[i]];
-                origin[index_order[i]] = this->Mod(origin[index_order[i]] + 1,
-                                                   index_order[i]);
+                origin[index_order[i]] += 1;
             }
         }
 
@@ -104,6 +120,7 @@ protected:
         for (int i = 0; i < N; ++i)
             result[i] = i;
 
+        /* Naïve bubble sort — enough for now */
         for (int i = 0; i < N; ++i)
             for (int j = i + 1; j < N; ++j)
                 if (pos[result[i]] < pos[result[j]])
@@ -112,11 +129,53 @@ protected:
         return result;
     }
 
-    inline int Mod(int value, int index) const
+    inline vec_t<float, N> GetGradient(vec_t<int, N> origin) const
     {
-        int const dim = this->m_gradients.GetSize()[index];
-        int const ret = value % dim;
-        return ret >= 0 ? ret : ret + dim;
+        /* Quick shuffle table:
+         * strings /dev/urandom | grep . -nm256 | sort -k2 -t: | sed 's|:.*|,|'
+         * Then just replace “256” with “0”. */
+        static int const shuffle[256] =
+        {
+            111, 14, 180, 186, 221, 114, 17, 79, 66, 46, 11, 81, 246, 200, 141,
+            172, 85, 244, 112, 92, 34, 106, 218, 205, 236, 7, 121, 115, 109,
+            131, 10, 96, 188, 148, 219, 107, 94, 182, 235, 163, 143, 213, 248,
+            202, 52, 154, 37, 241, 53, 129, 25, 60, 242, 38, 171, 63, 203, 255,
+            193, 6, 42, 209, 28, 176, 210, 159, 54, 144, 3, 71, 89, 116, 12,
+            237, 67, 216, 252, 178, 174, 164, 98, 234, 32, 26, 175, 24, 130,
+            128, 113, 99, 212, 62, 152, 75, 185, 73, 93, 31, 30, 151, 122, 173,
+            139, 91, 136, 162, 194, 208, 56, 101, 68, 69, 211, 44, 97, 55, 83,
+            33, 50, 119, 156, 149, 41, 157, 253, 247, 161, 47, 230, 166, 225,
+            204, 224, 13, 110, 123, 142, 64, 65, 155, 215, 120, 197, 140, 58,
+            77, 214, 126, 195, 179, 220, 232, 125, 147, 8, 39, 187, 27, 217,
+            100, 134, 199, 88, 206, 231, 250, 74, 2, 135, 9, 245, 118, 21, 243,
+            82, 183, 238, 87, 158, 61, 4, 177, 146, 153, 117, 249, 254, 233,
+            90, 222, 207, 48, 15, 18, 20, 239, 133, 0, 165, 138, 127, 169, 72,
+            1, 201, 145, 191, 192, 16, 49, 19, 95, 226, 228, 84, 181, 251, 36,
+            150, 22, 43, 70, 45, 105, 5, 189, 160, 196, 40, 59, 57, 190, 80,
+            104, 167, 78, 124, 103, 240, 184, 170, 137, 29, 23, 223, 108, 102,
+            86, 198, 227, 35, 229, 76, 168, 132, 51,
+        };
+
+        /* 16 random vectors; this should be enough for small dimensions */
+        auto v = [&]()
+        {
+            vec_t<float, N> ret;
+            for (int i = 0; i < N; ++i)
+                ret[i] = rand(-1.f, 1.f);
+            return normalize(ret);
+        };
+
+        static vec_t<float, N> const gradients[16] =
+        {
+            v(), v(), v(), v(), v(), v(), v(), v(),
+            v(), v(), v(), v(), v(), v(), v(), v(),
+        };
+
+        int idx = m_seed;
+        for (int i = 0; i < N; ++i)
+            idx = shuffle[(idx + origin[i]) & 255];
+
+        return gradients[(idx ^ (idx >> 4)) & 15];
     }
 
     /* For a given world position, extract grid coordinates (origin) and
@@ -131,46 +190,11 @@ protected:
 
         // Extracting decimal part from simplex sample
         pos = simplex_position - (vec_t<float, N>)origin;
-
-        // Never exceed the size of the gradients and loop on it
-        for (int i = 0; i < N; ++i)
-            origin[i] = this->Mod(origin[i], i);
-    }
-
-    inline void InitBase()
-    {
-        /*  Matrix coordinate transformation to skew simplex referential is done
-            by inversing the base matrix M which is written as follows:
-
-            M = | a b b b … |        M^(-1) = | c d d d … |
-                | b a b b … |                 | d c d d … |
-                | b b a b … |                 | d d c d … |
-                | …         |                 | …         |
-
-            where a, b, c, d are computed below ↴
-        */
-
-        float b = (1.f - lol::sqrt(N + 1.f)) / lol::sqrt((float)N * N * N);
-        float a = b + lol::sqrt((N + 1.f) / N);
-
-        // determinant of matrix M
-        float det = a * (a + (N - 2) * b) - (b * b) * (N - 1);
-
-        float c = (a + (N - 2) * b) / det;
-        float d = -b / det;
-
-        for (int i = 0; i < N; ++i)
-        {
-            for (int j = 0; j < N; ++j)
-            {
-                m_base[i][j] = (i == j ? a : b);
-                m_base_inverse[i][j] = (i == j ? c : d);
-            }
-        }
     }
 
     mat_t<float, N, N> m_base, m_base_inverse;
-    arraynd<N, vec_t<float, N> > m_gradients;
+    int m_seed;
 };
 
 }
+