Modernize sampler in testrender (AcademySoftwareFoundation#1534)

The sampler implementation in testrender was very close to the 2020 article from Brent Burley. This patch
catches up to best practices, including a faster owen scrambler from Nathan Vegdahl and recognizing that
the index scrambling can be achived by the same owen scrambler method (which saves a bunch of logic and runs
faster). Finally, avoid the redundant bit reversals and use compiler intrinsics when available to benefit
from hardware instructions (for ARM or GPUs). Finally the integer hash constants were updated to the latest
from the hash prospector project.

While more optimizations are certainly possible, this is a reasonable balance between quality and speed.

Update reference images for new sampler

Signed-off-by: Chris Kulla <[email protected]>

src/testrender/sampling.h

@@ -168,102 +168,80 @@ struct MIS {
     }
 };
 
-// Simple stratified progressive sampling using owen scrambled sobol points.
-// Code is written for clarity and simplicity over maximum speed.
+// "Practical Hash-based Owen Scrambling" - Brent Burley - JCGT 2020
+//    https://jcgt.org/published/0009/04/01/
 struct Sampler {
     Sampler(int px, int py, int si)
-        : seed(((px & 2047) << 22) | ((py & 2047) << 11)), si(si)
+        : seed(((px & 2047) << 22) | ((py & 2047) << 11)), index(reversebits(si))
     {
         assert(si < (1 << 24));
     }
 
     Vec3 get()
     {
         static const uint32_t zmatrix[24] = {
-            // 2^24 precision
-            0x800000u, 0xc00000u, 0x600000u, 0x900000u, 0xe80000u, 0x5c0000u,
-            0x8e0000u, 0xc50000u, 0x688000u, 0x9cc000u, 0xee6000u, 0x559000u,
-            0x806800u, 0xc09c00u, 0x60ee00u, 0x905500u, 0xe88080u, 0x5cc0c0u,
-            0x8e6060u, 0xc59090u, 0x6868e8u, 0x9c9c5cu, 0xeeee8eu, 0x5555c5u,
+            // 2^24 precision (reversed)
+            0x000001u, 0x000003u, 0x000006u, 0x000009u, 0x000017u, 0x00003au,
+            0x000071u, 0x0000a3u, 0x000116u, 0x000339u, 0x000677u, 0x0009aau,
+            0x001601u, 0x003903u, 0x007706u, 0x00aa09u, 0x010117u, 0x03033au,
+            0x060671u, 0x0909a3u, 0x171616u, 0x3a3939u, 0x717777u, 0xa3aaaau
         };
         seed += 4;  // advance depth for next call
-        uint32_t index = progressive_permute(si, hash(seed - 4));
-        uint32_t px = 0, py = 0, pz = 0, dx = 0x800000u, dy = 0x800000u;
-        for (int c = 0; index; c++, index >>= 1) {
-            if (index & 1) {
-                px ^= dx;
-                py ^= dy;
-                pz ^= zmatrix[c];
-            }
-            dx >>= 1;
-            dy ^= dy >> 1;
-        }  // scramble and scale by 2^-24
-        return { owen_scramble(px, seed - 3) * 5.96046448e-08f,
-                 owen_scramble(py, seed - 2) * 5.96046448e-08f,
-                 owen_scramble(pz, seed - 1) * 5.96046448e-08f };
+        uint32_t scrambled_index = owen_scramble(index, hash(seed - 4)) & 0xFFFFFF;
+        uint32_t result_x = scrambled_index; // already reversed
+        uint32_t result_y = 0;
+        uint32_t result_z = 0;
+        uint32_t ymatrix = 1;
+        for (int c = 0; c < 24; c++) {
+            uint32_t bit = (scrambled_index >> c) & 1;
+            result_y ^= bit * ymatrix;
+            result_z ^= bit * zmatrix[c];
+            ymatrix ^= ymatrix << 1; // generate procedurally instead of storing this
+        }
+        // scramble results and scale by 2^-24 to guarantee equally spaced values in [0,1)
+        return { (owen_scramble(result_x, hash(seed - 3)) >> 8) * 5.96046448e-8f,
+                 (owen_scramble(result_y, hash(seed - 2)) >> 8) * 5.96046448e-8f,
+                 (owen_scramble(result_z, hash(seed - 1)) >> 8) * 5.96046448e-8f };
     }
 
 private:
-    uint32_t seed, si;
+    uint32_t seed, index;
 
     static uint32_t hash(uint32_t s)
     {
-        // https://nullprogram.com/blog/2018/07/31/
+        // https://github.com/skeeto/hash-prospector
         s ^= s >> 16;
-        s *= 0x7feb352du;
+        s *= 0x21f0aaadu;
+        s ^= s >> 15;
+        s *= 0xd35a2d97u;
         s ^= s >> 15;
-        s *= 0x846ca68bu;
-        s ^= s >> 16;
         return s;
     }
-    static uint32_t progressive_permute(uint32_t si, uint32_t p)
-    {
-        // shuffle order of points in power of 2 blocks
-        if (si < 4)
-            return cmj_permute(si, 4, p);
-        uint32_t l = si;
-        l          = l | (l >> 1);
-        l          = l | (l >> 2);
-        l          = l | (l >> 4);
-        l          = l | (l >> 8);
-        l          = l | (l >> 16);
-        l          = l - (l >> 1);
-        return cmj_permute(si - l, l, p) + l;
-    }
-    static inline uint32_t cmj_permute(uint32_t i, uint32_t l, uint32_t p)
+
+    static uint32_t reversebits(uint32_t x)
     {
-        // in-place random permutation (power of 2), see:
-        // "Correlated Multi-Jittered Sampling" by "Andrew Kensler"
-        const uint32_t w = l - 1;
-        assert((l & w) == 0);
-        i ^= p;
-        i *= 0xe170893d;
-        i ^= p >> 16;
-        i ^= (i & w) >> 4;
-        i ^= p >> 8;
-        i *= 0x0929eb3f;
-        i ^= p >> 23;
-        i ^= (i & w) >> 1;
-        i *= 1 | p >> 27;
-        i *= 0x6935fa69;
-        i ^= (i & w) >> 11;
-        i *= 0x74dcb303;
-        i ^= (i & w) >> 2;
-        i *= 0x9e501cc3;
-        i ^= (i & w) >> 2;
-        i *= 0xc860a3df;
-        i &= w;
-        i ^= i >> 5;
-        return (i + p) & w;
+#if defined(__clang__)
+        return __builtin_bitreverse32(x);
+#else
+        x = (x << 16) | (x >> 16);
+        x = ((x & 0x00ff00ff) << 8) | ((x & 0xff00ff00) >> 8);
+        x = ((x & 0x0f0f0f0f) << 4) | ((x & 0xf0f0f0f0) >> 4);
+        x = ((x & 0x33333333) << 2) | ((x & 0xcccccccc) >> 2);
+        x = ((x & 0x55555555) << 1) | ((x & 0xaaaaaaaa) >> 1);
+        return x;
+#endif
     }
+
     static uint32_t owen_scramble(uint32_t p, uint32_t s)
     {
-        for (uint32_t m = 1u << 23; m; m >>= 1) {
-            s = hash(s);  // randomize state
-            p ^= s & m;   // flip output (depending on state)
-            s ^= p & m;   // flip state  (depending on output)
-        }
-        return p;
+        // https://psychopath.io/post/2021_01_30_building_a_better_lk_hash
+        // assumes reversed input
+        p ^= p * 0x3d20adea;
+        p += s;
+        p *= (s>> 16) | 1;
+        p ^= p * 0x05526c56;
+        p ^= p * 0x53a22864;
+        return reversebits(p);
     }
 };