From 33469f5bb2f4eb393fdac4e5bd7ed72b6b6868af Mon Sep 17 00:00:00 2001
From: Alexander Gavrilov <angavrilov@gmail.com>
Date: Wed, 2 Oct 2013 18:55:41 +0400
Subject: [PATCH] Use a better hash function in the Perlin noise generator.

---
 library/include/modules/PerlinNoise.inc | 104 ++++++++++--------------
 library/include/modules/Random.h        |  16 ++--
 2 files changed, 52 insertions(+), 68 deletions(-)

diff --git a/library/include/modules/PerlinNoise.inc b/library/include/modules/PerlinNoise.inc
index 573968a95..6b7fec17b 100644
--- a/library/include/modules/PerlinNoise.inc
+++ b/library/include/modules/PerlinNoise.inc
@@ -13,7 +13,11 @@ namespace Random {
 template<class T>
 inline T s_curve(T t)
 {
-    return t * t * (3 - 2*t);
+    // Classical function
+    //return t * t * (3 - 2*t);
+
+    // 2002 version from http://mrl.nyu.edu/~perlin/paper445.pdf
+    return t * t * t * (t * (t * 6 - 15) + 10);
 }
 
 template<class T>
@@ -22,50 +26,46 @@ inline T lerp(T s, T a, T b)
     return a + s * (b-a);
 }
 
+// Dot product of VSIZE vectors pointed by pa, pb
+
+template<class T, unsigned i>
+struct DotProduct {
+    static inline T eval(T *pa, T *pb);
+};
+template<class T>
+struct DotProduct<T,0> {
+    static inline T eval(T *pa, T *pb) { return pa[0]*pb[0]; }
+};
+template<class T, unsigned i>
+inline T DotProduct<T,i>::eval(T *pa, T *pb) {
+    return DotProduct<T,i-1>::eval(pa, pb) + pa[i]*pb[i];
+}
+
 // Templates used to force unrolling and inlining of the loops
 
 template<class T, unsigned VSIZE, unsigned BITS, class IDXT>
 template<unsigned mask>
-struct PerlinNoise<T,VSIZE,BITS,IDXT>::Impl<mask,0> {
+struct PerlinNoise<T,VSIZE,BITS,IDXT>::Impl<mask,-1> {
     typedef typename PerlinNoise<T,VSIZE,BITS,IDXT>::Temp Temp;
-    static inline T dot(T *pa, T *pb);
-    static inline void setup(const T *pv, Temp *pt);
+    static inline void setup(PerlinNoise<T,VSIZE,BITS,IDXT> *, const T *, Temp *) {}
     static inline T eval(PerlinNoise<T,VSIZE,BITS,IDXT> *self, Temp *pt, unsigned idx, T *pq);
 };
 
-// Dot product of VSIZE vectors pointed by pa, pb
-
-template<class T, unsigned VSIZE, unsigned BITS, class IDXT>
-template<unsigned mask>
-inline T PerlinNoise<T,VSIZE,BITS,IDXT>::Impl<mask,0>::dot(T *pa, T *pb)
-{
-    return pa[0] * pb[0];
-}
-
-template<class T, unsigned VSIZE, unsigned BITS, class IDXT>
-template<unsigned mask, unsigned i>
-inline T PerlinNoise<T,VSIZE,BITS,IDXT>::Impl<mask,i>::dot(T *pa, T *pb)
-{
-    return Impl<mask,i-1>::dot(pa, pb) + pa[i] * pb[i];
-}
-
 // Initialization of the temporaries from input coordinates
 
 template<class T, unsigned VSIZE, unsigned BITS, class IDXT>
-template<unsigned mask>
-inline void PerlinNoise<T,VSIZE,BITS,IDXT>::Impl<mask,0>::setup(const T *pv, Temp *pt)
-{
-    T t = std::floor(*pv);
-    pt->s = s_curve(pt->r0 = *pv - t);
-    pt->b0 = unsigned(int32_t(t)) & mask;
-}
+template<unsigned mask, int i>
+inline void PerlinNoise<T,VSIZE,BITS,IDXT>::Impl<mask,i>::setup(
+    PerlinNoise<T,VSIZE,BITS,IDXT> *self, const T *pv, Temp *pt
+) {
+    Impl<mask,i-1>::setup(self, pv, pt);
 
-template<class T, unsigned VSIZE, unsigned BITS, class IDXT>
-template<unsigned mask, unsigned i>
-inline void PerlinNoise<T,VSIZE,BITS,IDXT>::Impl<mask,i>::setup(const T *pv, Temp *pt)
-{
-    Impl<mask,i-1>::setup(pv,pt);
-    Impl<mask,0>::setup(pv+i,pt+i);
+    T t = std::floor(pv[i]);
+    pt[i].s = s_curve(pt[i].r0 = pv[i] - t);
+
+    unsigned b = unsigned(int32_t(t));
+    pt[i].b0 = self->idxmap[i][b & mask];
+    pt[i].b1 = self->idxmap[i][(b+1) & mask];
 }
 
 // Main recursion. Uses tables from self and pt.
@@ -73,32 +73,22 @@ inline void PerlinNoise<T,VSIZE,BITS,IDXT>::Impl<mask,i>::setup(const T *pv, Tem
 
 template<class T, unsigned VSIZE, unsigned BITS, class IDXT>
 template<unsigned mask>
-inline T PerlinNoise<T,VSIZE,BITS,IDXT>::Impl<mask, 0>::eval(
+inline T PerlinNoise<T,VSIZE,BITS,IDXT>::Impl<mask, -1>::eval(
     PerlinNoise<T,VSIZE,BITS,IDXT> *self, Temp *pt, unsigned idx, T *pq
 ) {
-    pq[0] = pt[0].r0;
-    idx += pt[0].b0;
-    T u = Impl<mask,VSIZE-1>::dot(pq, self->gradients[idx]);
-
-    pq[0] -= 1;
-    idx += 1;
-    T v = Impl<mask,VSIZE-1>::dot(pq, self->gradients[idx]);
-
-    return lerp(pt[0].s, u, v);
+    return DotProduct<T,VSIZE-1>::eval(pq, self->gradients[idx]);
 }
 
 template<class T, unsigned VSIZE, unsigned BITS, class IDXT>
-template<unsigned mask, unsigned i>
+template<unsigned mask, int i>
 inline T PerlinNoise<T,VSIZE,BITS,IDXT>::Impl<mask,i>::eval(
     PerlinNoise<T,VSIZE,BITS,IDXT> *self, Temp *pt, unsigned idx, T *pq
 ) {
     pq[i] = pt[i].r0;
-    idx += pt[i].b0;
-    T u = Impl<mask,i-1>::eval(self, pt, self->idxmap[idx], pq);
+    T u = Impl<mask,i-1>::eval(self, pt, idx ^ pt[i].b0, pq);
 
     pq[i] -= 1;
-    idx += 1;
-    T v = Impl<mask,i-1>::eval(self, pt, self->idxmap[idx], pq);
+    T v = Impl<mask,i-1>::eval(self, pt, idx ^ pt[i].b1, pq);
 
     return lerp(pt[i].s, u, v);
 }
@@ -114,19 +104,13 @@ void PerlinNoise<T,VSIZE,BITS,IDXT>::init(MersenneRNG &rng)
     for (unsigned i = 0; i < TSIZE; i++)
         rng.unitvector(gradients[i], VSIZE);
 
-    // Random permutation table
-    for (unsigned i = 0; i < TSIZE; i++)
-        idxmap[i] = i;
-
-    rng.permute(idxmap, TSIZE);
-
-    // Extended part of the table to avoid doing bitwise ops
-    for (unsigned i = TSIZE; i < TSIZE_EXT; i++)
+    // Random permutation tables
+    for (unsigned j = 0; j < VSIZE; j++)
     {
-        for (unsigned j = 0; j < VSIZE; j++)
-            gradients[i][j] = gradients[i-TSIZE][j];
+        for (unsigned i = 0; i < TSIZE; i++)
+            idxmap[j][i] = i;
 
-        idxmap[i] = idxmap[i-TSIZE];
+        rng.permute(idxmap[j], TSIZE);
     }
 }
 
@@ -138,7 +122,7 @@ T PerlinNoise<T,VSIZE,BITS,IDXT>::eval(const T coords[VSIZE])
     // Temporary used to build the current offset vector
     T q[VSIZE];
 
-    Impl<TSIZE-1,VSIZE-1>::setup(coords, tmp);
+    Impl<TSIZE-1,VSIZE-1>::setup(this, coords, tmp);
 
     return Impl<TSIZE-1,VSIZE-1>::eval(this, tmp, 0, q);
 }
diff --git a/library/include/modules/Random.h b/library/include/modules/Random.h
index 67a9f904b..25d694ad7 100644
--- a/library/include/modules/Random.h
+++ b/library/include/modules/Random.h
@@ -109,6 +109,9 @@ namespace Random
     /*
      * Classical Perlin noise function in template form.
      * http://mrl.nyu.edu/~perlin/doc/oscar.html#noise
+     *
+     * Using an improved hash function from:
+     * http://www.cs.utah.edu/~aek/research/noise.pdf
      */
 
     template<class T, unsigned VSIZE, unsigned BITS = 8, class IDXT = uint8_t>
@@ -116,21 +119,18 @@ namespace Random
     {
         // Size of randomness tables
         static const unsigned TSIZE = 1<<BITS;
-        // Extended size with repeated data to avoid bitwise masking
-        static const unsigned TSIZE_EXT = 2*(TSIZE+1);
 
-        T gradients[TSIZE_EXT][VSIZE];
-        IDXT idxmap[TSIZE_EXT];
+        T gradients[TSIZE][VSIZE];
+        IDXT idxmap[VSIZE][TSIZE];
 
         // Templates used to unwind and inline recursion and loops
         struct Temp {
             T r0, s;
-            unsigned b0;
+            unsigned b0, b1;
         };
-        template<unsigned mask, unsigned i>
+        template<unsigned mask, int i>
         struct Impl {
-            static inline T dot(T *pa, T *pb);
-            static inline void setup(const T *pv, Temp *pt);
+            static inline void setup(PerlinNoise<T,VSIZE,BITS,IDXT> *self, const T *pv, Temp *pt);
             static inline T eval(PerlinNoise<T,VSIZE,BITS,IDXT> *self, Temp *pt, unsigned idx, T *pq);
         };