Update flake nodes to standard cellnoise functions (#2842)

This changelist updates the implementations of `flake2d` and `flake3d` to use standard cellnoise functions across all shading languages, including a new implementation in OSL.  The following specific changes are included:

- Replace the custom random number generator with standard cellnoise functions in all languages, using scalar cellnoise for density and priority rejection and vector cellnoise for flake rotation.
- Optimize `mx_cell_noise_vec3` in GLSL and `mx_cell_noise_float3` in MDL to share a single bjmix across all three output channels, matching the existing optimization in the native OSL version.
- Add OSL implementations of the `flake2d` and `flake3d` nodes, removing corresponding exclusions in the render test suite.
- Fix the GLSL rotation matrix to use correct sine and cosine terms in the Arvo fast rotation construction, aligning with the MDL implementation.
- Tighten the bounding sphere early-rejection test to the exact circumsphere of the flake cube.

Fixes issue #2819.

Author: jstone-lucasfilm

libraries/stdlib/genglsl/lib/mx_flake.glsl

@@ -1,25 +1,4 @@
-uint mx_flake_hash(uint seed, uint i)
-{
-    return (i ^ seed) * 1075385539u; 
-}
-
-uint mx_flake_init_seed(ivec3 i)
-{
-    return mx_flake_hash(mx_flake_hash(mx_flake_hash(0, i.x), i.y), i.z);
-}
-
-uint mx_flake_xorshift32(uint seed)
-{
-    seed ^= seed << 13;
-    seed ^= seed >> 17;
-    seed ^= seed << 5;
-    return seed;
-}
-
-float mx_uint_to_01(uint x)
-{
-    return float(x) / float(0xffffffffu);  // scale to [0, 1)
-}
+#include "mx_noise.glsl"
 
 // "Fast Random Rotation Matrices" by James Arvo, Graphics Gems3 P.117
 vec3 mx_rotate_flake(vec3 p, vec3 i)
@@ -35,12 +14,12 @@ vec3 mx_rotate_flake(vec3 p, vec3 i)
 
     float s_theta = sin(theta);
     float c_theta = cos(theta);
-    float sx = vx * s_theta - vy * s_theta;
-    float sy = vx * c_theta + vy * c_theta;
+    float sx = vx * c_theta - vy * s_theta;
+    float sy = vx * s_theta + vy * c_theta;
 
     mat3 m = mat3(
-        vx * sx - s_theta, vx * sy - s_theta, vx * vz,
-        vy * sx + c_theta, vy * sy - c_theta, vy * vz,
+        vx * sx - c_theta, vx * sy - s_theta, vx * vz,
+        vy * sx + s_theta, vy * sy - c_theta, vy * vz,
         vz * sx          , vz * sy          , 1.0 - z
     );
 
@@ -77,51 +56,41 @@ void mx_flake(
 
     vec3 P = position / vec3(size);
     vec3 base_P = floor(P);
-    ivec3 base_P_int = ivec3(base_P);
 
     // flake priority in [0..1], 0: no flake, flakes with higher priority shadow flakes "below" them
     float flake_priority = 0.0;
-    uint flake_seed = 0;
+    vec3 flake_cell = vec3(0.0);
 
-    // Examine the 3×3×3 lattice neighborhood around the sample cell. Flakes are seeded at cell
-    // centers but can overlap adjacent cells by up to flake_diameter, so neighbors may contribute
-    // at the sample position. For each neighbor we deterministically generate a seed, reject it
-    // by the density probability, compute a per-flake priority, and test the rotated, centered
-    // flake position for overlap. The highest-priority overlapping flake is selected.
+    // Examine the 3x3x3 neighborhood of cells around the sample position, selecting the
+    // highest-priority overlapping flake.
     for (int i = -1; i < 2; ++i)
     {
         for (int j = -1; j < 2; ++j)
         {
             for (int k = -1; k < 2; ++k)
             {
-                uint seed = mx_flake_init_seed(base_P_int + ivec3(i, j, k));
+                vec3 cell_pos = base_P + vec3(i, j, k);
 
-                seed = mx_flake_xorshift32(seed);
-                if (mx_uint_to_01(seed) > probability)
+                vec3 PP = P - cell_pos - vec3(0.5);
+                if (dot(PP, PP) >= flake_diameter * flake_diameter * 3.0)
                     continue;
 
-                seed = mx_flake_xorshift32(seed);
-                float priority = mx_uint_to_01(seed);
-                if (priority < flake_priority)
+                if (mx_cell_noise_float(cell_pos) > probability)
                     continue;
 
-                vec3 flake_P = base_P + vec3(i, j, k) + vec3(0.5);
-                vec3 PP = P - flake_P;
-                if (dot(PP, PP) >= flake_diameter * flake_diameter * 4.0)
+                float priority = mx_cell_noise_float(vec4(cell_pos, 3.0));
+                if (priority < flake_priority)
                     continue;
 
-                vec3 rot;
-                seed = mx_flake_xorshift32(seed); rot.x = mx_uint_to_01(seed);
-                seed = mx_flake_xorshift32(seed); rot.y = mx_uint_to_01(seed);
-                seed = mx_flake_xorshift32(seed); rot.z = mx_uint_to_01(seed);
+                vec3 rot = mx_cell_noise_vec3(cell_pos);
                 PP = mx_rotate_flake(PP, rot);
 
                 if (abs(PP.x) <= flake_diameter &&
                     abs(PP.y) <= flake_diameter &&
                     abs(PP.z) <= flake_diameter)
                 {
                     flake_priority = priority;
-                    flake_seed = seed;
+                    flake_cell = cell_pos;
                 }
             }
         }
@@ -138,12 +107,12 @@ void mx_flake(
     }
 
     // create a flake normal by importance sampling a microfacet distribution with given roughness
-    uint seed = flake_seed;
-    float xi0 = mx_uint_to_01(seed); seed = mx_flake_xorshift32(seed);
-    float xi1 = mx_uint_to_01(seed); seed = mx_flake_xorshift32(seed);
+    vec3 flake_noise = mx_cell_noise_vec3(vec4(flake_cell, 2.0));
+    float xi0 = flake_noise.x;
+    float xi1 = flake_noise.y;
 
-    id = int(seed);  // not ideal but MaterialX does not support unsigned integer type
-    rand = mx_uint_to_01(seed);
+    rand = flake_noise.z;
+    id = int(rand * 2147483647.0);
     presence = flake_priority;
 
     float phi = M_PI * 2.0 * xi0;

libraries/stdlib/genglsl/lib/mx_noise.glsl

@@ -378,10 +378,16 @@ vec3 mx_cell_noise_vec3(vec3 p)
     int ix = mx_floor(p.x);
     int iy = mx_floor(p.y);
     int iz = mx_floor(p.z);
+    uint a, b, c;
+    a = b = c = uint(0xdeadbeef) + (4u << 2u) + 13u;
+    a += uint(ix);
+    b += uint(iy);
+    c += uint(iz);
+    mx_bjmix(a, b, c);
     return vec3(
-            mx_bits_to_01(mx_hash_int(ix, iy, iz, 0)),
-            mx_bits_to_01(mx_hash_int(ix, iy, iz, 1)),
-            mx_bits_to_01(mx_hash_int(ix, iy, iz, 2))
+            mx_bits_to_01(mx_bjfinal(a,      b, c)),
+            mx_bits_to_01(mx_bjfinal(a + 1u, b, c)),
+            mx_bits_to_01(mx_bjfinal(a + 2u, b, c))
     );
 }
 
@@ -391,10 +397,17 @@ vec3 mx_cell_noise_vec3(vec4 p)
     int iy = mx_floor(p.y);
     int iz = mx_floor(p.z);
     int iw = mx_floor(p.w);
+    uint a, b, c;
+    a = b = c = uint(0xdeadbeef) + (5u << 2u) + 13u;
+    a += uint(ix);
+    b += uint(iy);
+    c += uint(iz);
+    mx_bjmix(a, b, c);
+    a += uint(iw);
     return vec3(
-            mx_bits_to_01(mx_hash_int(ix, iy, iz, iw, 0)),
-            mx_bits_to_01(mx_hash_int(ix, iy, iz, iw, 1)),
-            mx_bits_to_01(mx_hash_int(ix, iy, iz, iw, 2))
+            mx_bits_to_01(mx_bjfinal(a, b,      c)),
+            mx_bits_to_01(mx_bjfinal(a, b + 1u, c)),
+            mx_bits_to_01(mx_bjfinal(a, b + 2u, c))
     );
 }
 

libraries/stdlib/genosl/lib/mx_flake.osl

@@ -0,0 +1,129 @@
+// "Fast Random Rotation Matrices" by James Arvo, Graphics Gems3 P.117
+vector mx_rotate_flake(vector p, vector i)
+{
+    float theta = M_PI * 2.0 * i[0];
+    float phi = M_PI * 2.0 * i[1];
+    float z = i[2] * 2.0;
+
+    float r = sqrt(z);
+    float vx = sin(phi) * r;
+    float vy = cos(phi) * r;
+    float vz = sqrt(2.0 - z);
+
+    float s_theta = sin(theta);
+    float c_theta = cos(theta);
+    float sx = vx * c_theta - vy * s_theta;
+    float sy = vx * s_theta + vy * c_theta;
+
+    // Matrix-vector multiply (matching GLSL mat3 column-major layout)
+    vector result;
+    result[0] = (vx * sx - c_theta) * p[0] + (vy * sx + s_theta) * p[1] + (vz * sx) * p[2];
+    result[1] = (vx * sy - s_theta) * p[0] + (vy * sy - c_theta) * p[1] + (vz * sy) * p[2];
+    result[2] = (vx * vz) * p[0] + (vy * vz) * p[1] + (1.0 - z) * p[2];
+
+    return result;
+}
+
+// compute a flake probability for a given flake coverage density x
+float mx_flake_density_to_probability(float x)
+{
+    // constants for numerical fitted curve to observed flake noise density behavior
+    float a = -26.19771808;
+    float b = 26.39663835;
+    float c = 85.53857017;
+    float d = -102.35069432;
+    float e = -101.42634862;
+    float f = 118.45082288;
+    float xx = x * x;
+
+    return (a * xx + b * x) / (c * xx * x + d * xx + e * x + f);
+}
+
+void mx_flake(
+    float size,
+    float roughness,
+    float coverage,
+    vector position,
+    vector N,
+    vector tangent,
+    vector bitangent,
+    output int id,
+    output float rand,
+    output float presence,
+    output vector flakenormal
+)
+{
+    float probability = mx_flake_density_to_probability(clamp(coverage, 0.0, 1.0));
+    float flake_diameter = 1.5 / sqrt(3.0);
+
+    vector P = position / vector(size);
+    point base_P = point(floor(P));
+
+    // flake priority in [0..1], 0: no flake, flakes with higher priority shadow flakes "below" them
+    float flake_priority = 0.0;
+    point flake_cell = point(0);
+
+    // Examine the 3x3x3 neighborhood of cells around the sample position, selecting the
+    // highest-priority overlapping flake.
+    for (int i = -1; i < 2; ++i)
+    {
+        for (int j = -1; j < 2; ++j)
+        {
+            for (int k = -1; k < 2; ++k)
+            {
+                point cell_pos = base_P + vector(i, j, k);
+
+                vector PP = P - vector(cell_pos) - vector(0.5);
+                if (dot(PP, PP) >= flake_diameter * flake_diameter * 3.0)
+                    continue;
+
+                if ((float) cellnoise(cell_pos) > probability)
+                    continue;
+
+                float priority = (float) cellnoise(cell_pos, 3.0);
+                if (priority < flake_priority)
+                    continue;
+
+                vector rot = (vector) cellnoise(cell_pos);
+                PP = mx_rotate_flake(PP, rot);
+
+                if (abs(PP[0]) <= flake_diameter &&
+                    abs(PP[1]) <= flake_diameter &&
+                    abs(PP[2]) <= flake_diameter)
+                {
+                    flake_priority = priority;
+                    flake_cell = cell_pos;
+                }
+            }
+        }
+    }
+
+    if (flake_priority <= 0.0)
+    {
+        // no flake
+        id = 0;
+        rand = 0.0;
+        presence = 0.0;
+        flakenormal = N;
+        return;
+    }
+
+    // create a flake normal by importance sampling a microfacet distribution with given roughness
+    vector flake_noise = cellnoise(flake_cell, 2.0);
+    float xi0 = flake_noise[0];
+    float xi1 = flake_noise[1];
+
+    rand = flake_noise[2];
+    id = (int)(rand * 2147483647.0);
+    presence = flake_priority;
+
+    float phi = M_PI * 2.0 * xi0;
+    float tan_theta = roughness * roughness * sqrt(xi1) / sqrt(1.0 - xi1);  // GGX
+    float sin_theta = tan_theta / sqrt(1.0 + tan_theta * tan_theta);
+    float cos_theta = sqrt(1.0 - sin_theta * sin_theta);
+
+    flakenormal = tangent * cos(phi) * sin_theta +
+                  bitangent * sin(phi) * sin_theta +
+                  N * cos_theta;
+    flakenormal = normalize(flakenormal);
+}

libraries/stdlib/genosl/mx_flake2d.osl

@@ -0,0 +1,20 @@
+#include "lib/mx_flake.osl"
+
+void mx_flake2d(
+    float size,
+    float roughness,
+    float coverage,
+    vector2 texcoord,
+    vector N,
+    vector tangent,
+    vector bitangent,
+    output int id,
+    output float rand,
+    output float presence,
+    output vector flakenormal
+)
+{
+    // reuse the 3d flake implementation. this could be optimized by using a 2d flake implementation
+    vector position = vector(texcoord.x, texcoord.y, 0.0);
+    mx_flake(size, roughness, coverage, position, N, tangent, bitangent, id, rand, presence, flakenormal);
+}

libraries/stdlib/genosl/mx_flake3d.osl

@@ -0,0 +1,18 @@
+#include "lib/mx_flake.osl"
+
+void mx_flake3d(
+    float size,
+    float roughness,
+    float coverage,
+    vector position,
+    vector N,
+    vector tangent,
+    vector bitangent,
+    output int id,
+    output float rand,
+    output float presence,
+    output vector flakenormal
+)
+{
+    mx_flake(size, roughness, coverage, position, N, tangent, bitangent, id, rand, presence, flakenormal);
+}

libraries/stdlib/genosl/stdlib_genosl_impl.mtlx

@@ -176,6 +176,12 @@
   <implementation name="IM_worleynoise3d_vector2_genosl" nodedef="ND_worleynoise3d_vector2" file="mx_worleynoise3d_vector2.osl" function="mx_worleynoise3d_vector2" target="genosl" />
   <implementation name="IM_worleynoise3d_vector3_genosl" nodedef="ND_worleynoise3d_vector3" file="mx_worleynoise3d_vector3.osl" function="mx_worleynoise3d_vector3" target="genosl" />
 
+  <!-- <flake2d> -->
+  <implementation name="IM_flake2d_genosl" nodedef="ND_flake2d" file="mx_flake2d.osl" function="mx_flake2d" target="genosl" />
+
+  <!-- <flake3d> -->
+  <implementation name="IM_flake3d_genosl" nodedef="ND_flake3d" file="mx_flake3d.osl" function="mx_flake3d" target="genosl" />
+
   <!-- ======================================================================== -->
   <!-- Geometric nodes                                                          -->
   <!-- ======================================================================== -->