math: optimise Perlin noise by parsing the hypercube using a Gray code pattern.

src/lol/math/noise/perlin.h

@@ -34,51 +34,59 @@ public:
     /* Evaluate noise at a given point */
     inline float eval(vec_t<float, N> position) const
     {
-        int const cells = 1 << N;
-
-        /* Retrieve the containing hypercube origin */
+        /* Compute the containing hypercube origin */
         vec_t<int, N> origin;
         for (int i = 0; i < N; ++i)
             origin[i] = (int)position[i] - (position[i] < 0);
 
         vec_t<float, N> delta = position - (vec_t<float, N>)origin;
 
+        /* Apply a smooth step to delta and store it in “t”. */
         vec_t<float, N> t = delta;
-        /* DEBUG: original Perlin noise polynomial */
-        //t = (vec_t<float, N>(3.f) - 2.f * t) * t * t;
-        /* Improved polynomial (with null second derivative) */
         t = ((6.f * t - vec_t<float, N>(15.f))
                   * t + vec_t<float, N>(10.f)) * (t * t * t);
+        /* DEBUG: original Perlin noise polynomial */
+        //t = (vec_t<float, N>(3.f) - 2.f * t) * t * t;
 
-        /* Compute all gradient contributions */
-        array<float> values;
-        values.Resize(cells);
-        for (int i = 0; i < cells; ++i)
+        /* Premultiply and predivide (1-t)/t and t/(1-t) into “u” and “v”. */
+        vec_t<float, N> u, v;
+        float multiplier = 1.f;
+        for (int bit = 0; bit < N; ++bit)
         {
-            /* FIXME: find another way to travel the hypercube that
-             * causes fewer bit flips, and compute the return value
-             * directly without allocating an array value. */
-            vec_t<float, N> v = delta;
-            vec_t<int, N> corner = origin;
-            for (int bit = 0; bit < N; ++bit)
-                if ((1 << bit) & i)
-                {
-                    ++corner[bit];
-                    v[bit] = v[bit] - 1.f;
-                }
-
-            values[i] = dot(v, this->get_gradient(corner));
+            /* Avoid divisions by zero near the hypercube boundaries. */
+            float f = clamp(t[bit], 0.001f, 0.999f);
+
+            multiplier *= (1.f - f);
+            u[bit] = (1.f - f) / f;
+            v[bit] = f / (1.f - f);
         }
 
-        /* Interpolate all contributions together */
-        for (int bit = N; bit--; )
+        float ret = 0.f;
+
+        /* Compute all gradient contributions, for each of the 2^N corners
+         * of the hypercube. */
+        for (int i = 0; i < (1 << N); ++i)
         {
-            for (int i = 0; i < (1 << bit); ++i)
-                values[i] = (1.f - t[bit]) * values[i]
-                                  + t[bit] * values[i + (1 << bit)];
+            /* Accumulate Perlin noise */
+            ret += multiplier * dot(delta, this->get_gradient(origin));
+
+            /* Don’t use the binary pattern for “i” but use its Gray code
+             * “j” instead, so we know we only have one component to alter
+             * in “origin” and in “delta”. We know which bit was flipped by
+             * looking at “k”, the Gray code for the next value of “i”. */
+            int j = i ^ (i >> 1);
+            int k = (i + 1) ^ ((i + 1) >> 1);
+
+            int bit = 0;
+            while ((j ^ k) > (1 << bit))
+                ++bit;
+
+            origin[bit] += j > k ? -1 : 1;
+            delta[bit] += j > k ? 1.f : -1.f;
+            multiplier *= (j > k ? u : v)[bit];
         }
 
-        return values[0];
+        return sqrt(2.f) * ret;
     }
 };