testcase.glsl

// 3D to 3D hash function
uvec3 hash3(uvec3 v) 
{
    v = v * 1664525u + 1013904223u;      
    v.x += v.y*v.z; v.y += v.z*v.x; v.z += v.x*v.y;
    v.x += v.y*v.z; v.y += v.z*v.x; v.z += v.x*v.y;
    return v >> 16u;
}
// build a 3D curl nosie into a 2D texture atlas

// 3D simplex weights (& corners)
// based on McEwan et al., Efficient computation of noise in GLSL, JCT 2011
mat4x3 simplexCoords(vec3 p) 
{
    // skew to tetrahedral coordinates
    vec3 tetbase = floor(p + dot(p, vec3(1./3.)));
    vec3 base = tetbase - dot(tetbase, vec3(1./6.));
    vec3 tf = p - base;

    // One of six tetrahedra: 100, 010, 001, 011, 101, 110
    // since skew is along x=y=z axis, this works in Euclidean space too
    vec3 g = step(tf.yzx, tf.xyz), h = 1. - g.zxy;
    vec3 a1 = min(g, h) - 1./6., a2 = max(g, h) - 1./3.;

    // four corners in Euclidean space
    return mat4x3(base, base + a1, base + a2, base + 0.5);
}

// simplex smoothing function
vec4 Smooth(mat4x3 f)
{
    const float scale = 1024. / 375.;       // scale factor to make noise -1..1
    vec4 d = vec4(dot(f[0], f[0]), dot(f[1], f[1]), dot(f[2], f[2]), dot(f[3], f[3]));
    vec4 s = clamp(2. * d, 0., 1.);
    
    return (1. * scale + s*(-3. * scale + s*(3. * scale - s*scale)));
}

// derivative of simplex noise smoothing function
mat3x4 dSmooth(mat4x3 f)
{
    const float scale = 1024. / 375.;       // scale factor to make noise -1..1
    vec4 d = vec4(dot(f[0], f[0]), dot(f[1], f[1]), dot(f[2], f[2]), dot(f[3], f[3]));
    vec4 s = clamp(2. * d, 0., 1.);
    
    s = -12. * scale + s*(24. * scale - s * 12. * scale);

    return mat3x4(
        s * vec4(f[0][0], f[1][0], f[2][0], f[3][0]),
        s * vec4(f[0][1], f[1][1], f[2][1], f[3][1]),
        s * vec4(f[0][2], f[1][2], f[2][2], f[3][2]));
}

// Simplex curl noise, can seamlessly wrap at multiples of 3
vec3 CurlSimplex(vec3 v)
{
    const uvec3 GMask = uvec3(0x8000, 0x4000, 0x2000);
    const vec3 GScale = 1. / vec3(0x4000, 0x2000, 0x1000);
    
    // corners of tetrahedron
    mat4x3 T = simplexCoords(v), gvec[3], fv;
    mat3x4 grad;
    
    for(int i=0; i<4; ++i) {
        fv[i] = v - T[i];
        uvec3 rand = hash3(uvec3(6. * T[i] + 0.5));
        gvec[0][i] = vec3(rand.xxx & GMask) * GScale - 1.;
        gvec[1][i] = vec3(rand.yyy & GMask) * GScale - 1.;
        gvec[2][i] = vec3(rand.zzz & GMask) * GScale - 1.;
        grad[0][i] = dot(gvec[0][i], fv[i]);
        grad[1][i] = dot(gvec[1][i], fv[i]);
        grad[2][i] = dot(gvec[2][i], fv[i]);
    }

    // blend gradients
    vec4 sv = Smooth(fv);
    mat3x4 ds = dSmooth(fv);

    // compute Jacobian, rely on compiler to get rid of unneeded elements
    mat3x3 J;
    J[0] = vec3((gvec[0] * sv) + (grad[0] * ds));
    J[1] = vec3((gvec[1] * sv) + (grad[1] * ds));
    J[2] = vec3((gvec[2] * sv) + (grad[2] * ds));

    // curl from Jacobian
    return vec3(J[1][2] - J[2][1], J[2][0]-J[0][2], J[0][1]-J[1][0]);
}