math: shitloads of tweaks, optimisations, fixes and comments to the

simplex noise code.

demos/test/simplex.cpp

@@ -17,44 +17,108 @@
 
 using namespace lol;
 
+/* Constants to tweak */
+float const zoom = 0.03f;
+int const octaves = 1;
+
 int main(int argc, char **argv)
 {
     UNUSED(argc, argv);
 
-    float const zoom = 0.03f;
-    int const octaves = 10;
-
     /* Create an image */
     ivec2 const size(1280, 720);
     Image img(size);
     array2d<vec4> &data = img.Lock2D<PixelFormat::RGBA_F32>();
 
-    /* Fill image with simplex noise */
+    /* Declare plenty of allocators */
     simplex_interpolator<2> s2;
     simplex_interpolator<3> s3;
+    simplex_interpolator<4> s4;
+    simplex_interpolator<5> s5;
+    simplex_interpolator<6> s6;
     simplex_interpolator<7> s7;
+    simplex_interpolator<8> s8;
+    simplex_interpolator<9> s9;
+    simplex_interpolator<10> s10;
+    simplex_interpolator<11> s11;
+    simplex_interpolator<12> s12;
 
+    /* Fill image with simplex noise */
     for (int j = 0; j < size.y; ++j)
     for (int i = 0; i < size.x; ++i)
     {
+        int cell = j / (size.y / 2) * 3 + i / (size.x / 3);
+
         float x = (float)i, y = (float)j;
-        float sum = 0.f;
-        int maxoctave = (j < size.y / 2) ? 1 : (1 << octaves);
+        float sum = 0.f, coeff = 0.f;
+        int const max_k = 1 << octaves;
 
-        for (int k = 1; k <= maxoctave; k *= 2)
+        for (int k = 1; k < max_k; k *= 2)
         {
-            if (i < size.x / 3)
-                sum += 0.5f / k * s2.Interp(zoom * k * vec2(x, y));
-            else if (i < size.x * 2 / 3)
-                sum += 0.5f / k * s3.Interp(zoom * k * vec3(x, y, 0.f));
-            else
-                sum += 0.5f / k * s7.Interp(zoom * k * vec7(x, 0.f, 0.f, 0.f, y, 0.f, 0.f));
+            float t = 0.f;
+
+            switch (cell)
+            {
+            case 0:
+                t = s2.Interp(zoom * k * vec2(x, y));
+                break;
+            case 1:
+                t = s3.Interp(zoom * k * vec3(x, 1.f, y));
+                break;
+            case 2:
+                t = s4.Interp(zoom * k * vec4(x, 1.f, y, 2.f));
+                break;
+            case 3:
+                t = s6.Interp(zoom * k * vec6(x, 1.f, 2.f, y, 3.f, 4.f));
+                break;
+            case 4:
+                t = s8.Interp(zoom * k * vec8(x, 1.f, 2.f, 3.f, y, 4.f,
+                                              5.f, 6.f));
+                break;
+            case 5:
+                t = s12.Interp(zoom * k * vec12(x, 1.f, 2.f, 3.f, 4.f, 5.f,
+                                                y, 6.f, 7.f, 8.f, 9.f, 10.f));
+                break;
+            default:
+                break;
+            }
+
+            sum += 1.f / k * t;
+            coeff += 1.f / k;
         }
 
-        float c = saturate(0.5f + sum);
+        float c = saturate(0.5f + 0.5f * sum / coeff);
         data[i][j] = vec4(c, c, c, 1.f);
+        //data[i][j] = Color::HSVToRGB(vec4(c, 1.0f, 0.5f, 1.f));
     }
 
+#if 0
+    /* 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)
+    {
+        auto putpixel = [&](ivec2 p, vec4 c = vec4(1.f, 0.f, 1.f, 1.f))
+        {
+            if (p.x >= 0 && p.x < size.x - 1 && p.y >= 0 && p.y < size.y - 1)
+                data[p.x][p.y] = c;
+        };
+
+        ivec2 coord = ivec2(i / zoom * dx + j / zoom * dy);
+
+        vec2 g = s2.GetGradient((i + 0.1f) * dx + (j + 0.1f) * dy);
+        for (int t = 0; t < 40; ++t)
+            putpixel(coord + (ivec2)(g * (t / 2.f)), vec4(0.f, 1.f, 0.f, 1.f));
+
+        putpixel(coord);
+        putpixel(coord + ivec2(1, 0));
+        putpixel(coord + ivec2(0, 1));
+        putpixel(coord + ivec2(1, 1));
+    }
+#endif
+
     /* Save image */
     img.Unlock2D(data);
     img.Save("simplex.png");

src/lol/math/simplex_interpolator.h

@@ -17,6 +17,27 @@
 namespace lol
 {
 
+/*
+ * Simplex noise in dimension N
+ * ----------------------------
+ *
+ *  The N-dimensional regular hypercube can be split into N! simplices that
+ * all have the main diagonal as a shared edge.
+ *  - number of simplices: N!
+ *  - number of vertices per simplex: N+1
+ *  - number of edges: N(N+1)/2
+ *  - minimum edge length: 1           (hypercube edges, e.g. [1,0,0,…,0])
+ *  - maximum edge length: sqrt(N)     (hypercube diagonal, i.e. [1,1,1,…,1])
+ *
+ *  We skew the simplicial grid along the main diagonal by a factor f =
+ * sqrt(N+1), which means the diagonal of the initial parallelotope has
+ * length sqrt(N/(N+1)). The edges of that parallelotope have length
+ * sqrt(N/(N+1)), too. A formula for the maximum edge length was found
+ * empirically:
+ *  - minimum edge length: sqrt(N/(N+1))   (parallelotope edges and diagonal)
+ *  - maximum edge length: sqrt(floor((N+1)²/4)/(N+1))
+ */
+
 template<int N>
 class simplex_interpolator
 {
@@ -24,165 +45,218 @@ public:
     simplex_interpolator(int seed = 0)
       : m_seed(seed)
     {
-        /* 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 ↴
-         */
-
-        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;
-
+#if 0
+        // Print some debug information
+        printf("Simplex Noise of Dimension %d\n", N);
+
+        long long int n = 1; for (int i = 1; i <= N; ++i) n *= i;
+        printf(" - each hypercube cell has %lld simplices "
+               "with %d vertices and %d edges\n", n, N + 1, N * (N + 1) / 2);
+
+        vec_t<float, N> diagonal(1.f);
+        printf(" - length of hypercube edges is 1, diagonal is %f\n",
+               length(diagonal));
+        printf(" - length of parallelotope edges and diagonal is %f\n",
+               length(unskew(diagonal)));
+
+        vec_t<float, N> vertices[N + 1];
+        vertices[0] = vec_t<float, N>(0.f);
         for (int i = 0; i < N; ++i)
         {
-            for (int j = 0; j < N; ++j)
+            vertices[i + 1] = vertices[i];
+            vertices[i + 1][i] += 1.f;
+        }
+        float minlength = 10.0f;
+        float maxlength = 0.0f;
+        for (int i = 0; i < N + 1; ++i)
+            for (int j = i + 1; j < N + 1; ++j)
             {
-                m_base[i][j] = (i == j ? a : b);
-                m_base_inverse[i][j] = (i == j ? c : d);
+                float l = length(unskew(vertices[i] - vertices[j]));
+                minlength = min(minlength, l);
+                maxlength = max(maxlength, l);
             }
-        }
+        printf(" - edge lengths between %f and %f\n", minlength, maxlength);
+        printf(" - predicted: %f and %f\n", sqrt((float)N/(N+1)),
+                                            sqrt((N+1)*(N+1)/4/(float)(N+1)));
+        printf("\n");
+#endif
     }
 
-    // Single interpolation
+    /* Single interpolation */
     inline float Interp(vec_t<float, N> position) const
     {
-        // Computing position in simplex referential
-        vec_t<float, N> simplex_position = m_base_inverse * position;
-
-        // Retrieving the closest floor point and decimals
+        // Retrieve the containing hypercube origin and associated decimals
         vec_t<int, N> origin;
         vec_t<float, N> pos;
+        get_origin(skew(position), origin, pos);
 
-        this->ExtractFloorDecimal(simplex_position, origin, pos);
-
-        vec_t<int, N> index_order = this->GetIndexOrder(pos);
-
-        // Resetting decimal point in regular orthonormal coordinates
-        pos = m_base * pos;
+        return get_noise(origin, pos);
+    }
 
-        return this->LastInterp(origin, pos, index_order);
+    /* Only for debug purposes: return the gradient vector */
+    inline vec_t<float, N> GetGradient(vec_t<float, N> position) const
+    {
+        vec_t<int, N> origin;
+        vec_t<float, N> pos;
+        get_origin(skew(position), origin, pos);
+        return get_gradient(origin);
     }
 
 protected:
-    inline float LastInterp(vec_t<int, N> origin,
-                            vec_t<float, N> const & pos,
-                            vec_t<int, N> const & index_order) const
+    inline float get_noise(vec_t<int, N> origin,
+                           vec_t<float, N> const & pos) const
     {
-        /* “corner” will traverse the simplex along its edges in
-         * orthonormal coordinates. */
-        vec_t<float, N> corner(0);
-        float result = 0;
+        /* For a given position [0…1]^N inside a regular N-hypercube, find
+         * the N-simplex which contains that position, and return a path
+         * along the hypercube edges from (0,0,…,0) to (1,1,…,1) which
+         * uniquely describes that simplex. */
+        vec_t<int, N> traversal_order;
+        for (int i = 0; i < N; ++i)
+            traversal_order[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[traversal_order[i]] < pos[traversal_order[j]])
+                    std::swap(traversal_order[i], traversal_order[j]);
+
+
+        /* Get the position in world coordinates, too */
+        vec_t<float, N> world_pos = unskew(pos);
+
+        /* “corner” will traverse the simplex along its edges in world
+         * coordinates. */
+        vec_t<float, N> world_corner(0.f);
+        float result = 0.f, coeff = 0.f;
 
         for (int i = 0; i < N + 1; ++i)
         {
-            vec_t<float, N> const delta = pos - corner;
-            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
-             * holds for higher dimensions. */
-            float dist = 1.0f - 4.f / 3.f * sqlength(delta);
-            if (dist > 0)
+#if 1
+            // FIXME: In “Noise Hardware” (2-17) Perlin uses 0.6 but Gustavson
+            // makes an exception for dimension 2 and uses 0.5 instead.
+            // I think m should actually increase with the dimension count.
+            float const m = N == 2 ? 0.5f : 0.6f;
+            float d = m - sqlength(world_pos - world_corner);
+#else
+            // DEBUG: this is the linear contribution of each vertex
+            // in the skewed simplex. Unfortunately it creates artifacts.
+            float d = ((i == 0) ? 1.f : pos[traversal_order[i - 1]])
+                    - ((i == N) ? 0.f : pos[traversal_order[i]]);
+#endif
+
+            if (d > 0)
             {
-                dist *= dist;
-                result += dist * dist * dot(gradient, delta);
+                vec_t<float, N> const &gradient = get_gradient(origin);
+
+                // Perlin uses 8d⁴ whereas Gustavson uses d⁴ and a final
+                // multiplication factor at the end.
+                //d = 8.f * d * d * d * d;
+                d = d * d * d * d;
+
+                //d = (3.f - 2.f * d) * d * d;
+                //d = ((6 * d - 15) * d + 10) * d * d * d;
+
+                result += d * dot(gradient, world_pos - world_corner);
+                coeff += d;
             }
 
             if (i < N)
             {
-                corner += m_base[index_order[i]];
-                origin[index_order[i]] += 1;
+                vec_t<float, N> v(0.f);
+                v[traversal_order[i]] = 1.f;
+                world_corner += unskew(v);
+                origin[traversal_order[i]] += 1;
             }
         }
 
-        // FIXME: another paper uses the value 70 for dimension 2, 32 for
-        // dimension 3, and 27 for dimension 4; find where this comes from.
-        return result * 70.f / 16.f;
-    }
-
-    /* For a given position [0…1]^n inside a square/cube/hypercube etc.,
-     * find the simplex which contains that position, and return the path
-     * from (0,0,…,0) to (1,1,…,1) that describes that simplex. */
-    inline vec_t<int, N> GetIndexOrder(vec_t<float, N> const & pos) const
-    {
-        vec_t<int, N> result;
-        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]])
-                    std::swap(result[i], result[j]);
-
-        return result;
+        // FIXME: Gustavson uses the value 70 for dimension 2, 32 for
+        // dimension 3, and 27 for dimension 4; find where this comes from
+        // and maybe find a reasonable formula.
+        //return 70.f * result / coeff;
+        float const k = N == 2 ? 70 : N == 3 ? 32 : N == 4 ? 27 : 20;
+        return k * result;
     }
 
-    inline vec_t<float, N> GetGradient(vec_t<int, N> origin) const
+    inline vec_t<float, N> get_gradient(vec_t<int, N> origin) const
     {
         /* 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,
+            111, 14, 180, 186, 221, 114, 219, 79, 66, 46, 152, 81, 246, 200,
+            141, 172, 85, 244, 112, 92, 34, 106, 218, 205, 236, 7, 121, 115,
+            109, 131, 10, 96, 188, 148, 17, 107, 94, 182, 235, 163, 143, 63,
+            248, 202, 52, 154, 37, 241, 53, 129, 25, 159, 242, 38, 171, 213,
+            6, 203, 255, 193, 42, 209, 28, 176, 210, 60, 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, 11, 75, 185, 73, 93, 31, 30, 44,
+            122, 173, 139, 91, 136, 162, 194, 41, 56, 101, 68, 69, 211, 151,
+            97, 55, 83, 33, 50, 119, 156, 149, 208, 157, 253, 247, 161, 133,
+            230, 166, 225, 204, 224, 13, 110, 123, 142, 64, 65, 155, 215, 9,
+            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, 120,
+            21, 245, 118, 243, 82, 183, 238, 150, 158, 61, 4, 177, 146, 153,
+            117, 249, 254, 233, 90, 222, 207, 48, 15, 18, 20, 16, 47, 0, 51,
+            165, 138, 127, 169, 72, 1, 201, 145, 191, 192, 239, 49, 19, 160,
+            226, 228, 84, 181, 251, 36, 87, 22, 43, 70, 45, 105, 5, 189, 95,
+            40, 196, 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,
         };
 
-        /* 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);
-        };
+        /* Generate 2^(N+2) random vectors, but at least 2^5 (32) and not
+         * more than 2^20 (~ 1 million). */
+        int const gradient_count = 1 << min(max(N + 2, 5), 20);
 
-        static vec_t<float, N> const gradients[16] =
+        static auto build_gradients = [&]()
         {
-            v(), v(), v(), v(), v(), v(), v(), v(),
-            v(), v(), v(), v(), v(), v(), v(), v(),
+            array<vec_t<float, N>> ret;
+            for (int k = 0; k < gradient_count; ++k)
+            {
+                vec_t<float, N> v;
+                for (int i = 0; i < N; ++i)
+                    v[i] = rand(-1.f, 1.f);
+                ret << normalize(v);
+            }
+            return ret;
         };
 
+        static array<vec_t<float, N>> const gradients = build_gradients();
+
         int idx = m_seed;
         for (int i = 0; i < N; ++i)
-            idx = shuffle[(idx + origin[i]) & 255];
+            idx ^= shuffle[(idx + origin[i]) & 255];
+
+        idx &= (gradient_count - 1);
+#if 0
+        // DEBUG: only output a few gradients
+        if (idx > 2)
+            return vec_t<float, N>(0);
+#endif
+        return gradients[idx];
+    }
 
-        return gradients[(idx ^ (idx >> 4)) & 15];
+    static inline vec_t<float, N> skew(vec_t<float, N> const &v)
+    {
+        /* Quoting Perlin in “Hardware Noise” (2-18):
+         *   The “skew factor” f should be set to f = sqrt(N+1), so that
+         *   the point (1,1,...1) is transformed to the point (f,f,...f). */
+        float const sum = dot(v, vec_t<float, N>(1));
+        float const f = sqrt(1.f + N);
+        return v + vec_t<float, N>(sum * (f - 1) / N);
+    }
+
+    static inline vec_t<float, N> unskew(vec_t<float, N> const &v)
+    {
+        float const sum = dot(v, vec_t<float, N>(1));
+        float const f = sqrt(1.f + N);
+        return v + vec_t<float, N>(sum * (1 / f - 1) / N);
     }
 
     /* For a given world position, extract grid coordinates (origin) and
      * the corresponding delta position (pos). */
-    inline void ExtractFloorDecimal(vec_t<float, N> const & simplex_position,
-                                    vec_t<int, N> & origin,
-                                    vec_t<float, N> & pos) const
+    inline void get_origin(vec_t<float, N> const & simplex_position,
+                           vec_t<int, N> & origin, vec_t<float, N> & pos) const
     {
         // Finding floor point index
         for (int i = 0; i < N; ++i)
@@ -192,7 +266,7 @@ protected:
         pos = simplex_position - (vec_t<float, N>)origin;
     }
 
-    mat_t<float, N, N> m_base, m_base_inverse;
+    /* A user-provided random seed. Defaults to zero. */
     int m_seed;
 };