LYGIA Shader Library

psrdnoise (lygia/generative/psrdnoise)

Tiling simplex flow noise in 2-D and 3-D from https://github.com/stegu/psrdnoise "vec2/3 x" is the point to evaluate, "vec2/3 period" is the desired periods "float alpha" is the rotation (in radians) for the swirling gradients. The "float" return value is the noise value n, the "out gradient" argument returns the first order derivatives, and "out dg" returns the second order derivatives as (dn2/dx2, dn2/dy2, dn2/dxy) For the original code please visit: https://github.com/stegu/psrdnoise

Dependencies:

Use:

<float> psrdnoise(<vec2|vec3> x, <vec2|vec3> period, <float> alpha, out <vec2|vec3> gradient [, out <vec3> dg])

Check it on Github



#ifndef FNC_PSRFNOISE
#define FNC_PSRFNOISE

float psrdnoise(vec2 x, vec2 period, float alpha, out vec2 gradient) {

    // Transform to simplex space (axis-aligned hexagonal grid)
    vec2 uv = vec2(x.x + x.y*0.5, x.y);

    // Determine which simplex we're in, with i0 being the "base"
    vec2 i0 = floor(uv);
    vec2 f0 = fract(uv);
    // o1 is the offset in simplex space to the second corner
    float cmp = step(f0.y, f0.x);
    vec2 o1 = vec2(cmp, 1.0-cmp);

    // Enumerate the remaining simplex corners
    vec2 i1 = i0 + o1;
    vec2 i2 = i0 + vec2(1.0, 1.0);

    // Transform corners back to texture space
    vec2 v0 = vec2(i0.x - i0.y * 0.5, i0.y);
    vec2 v1 = vec2(v0.x + o1.x - o1.y * 0.5, v0.y + o1.y);
    vec2 v2 = vec2(v0.x + 0.5, v0.y + 1.0);

    // Compute vectors from v to each of the simplex corners
    vec2 x0 = x - v0;
    vec2 x1 = x - v1;
    vec2 x2 = x - v2;

    vec3 iu = vec3(0.0);
    vec3 iv = vec3(0.0);
    vec3 xw = vec3(0.0);
    vec3 yw = vec3(0.0);

    // Wrap to periods, if desired
    if(any(greaterThan(period, vec2(0.0)))) {
        xw = vec3(v0.x, v1.x, v2.x);
        yw = vec3(v0.y, v1.y, v2.y);
        if(period.x > 0.0)
            xw = mod(vec3(v0.x, v1.x, v2.x), period.x);
        if(period.y > 0.0)
            yw = mod(vec3(v0.y, v1.y, v2.y), period.y);
        // Transform back to simplex space and fix rounding errors
        iu = floor(xw + 0.5*yw + 0.5);
        iv = floor(yw + 0.5);
    } else { // Shortcut if neither x nor y periods are specified
        iu = vec3(i0.x, i1.x, i2.x);
        iv = vec3(i0.y, i1.y, i2.y);
    }

    // Compute one pseudo-random hash value for each corner
    vec3 hash = mod(iu, 289.0);
    hash = mod((hash*51.0 + 2.0)*hash + iv, 289.0);
    hash = mod((hash*34.0 + 10.0)*hash, 289.0);

    // Pick a pseudo-random angle and add the desired rotation
    vec3 psi = hash * 0.07482 + alpha;
    vec3 gx = cos(psi);
    vec3 gy = sin(psi);

    // Reorganize for dot products below
    vec2 g0 = vec2(gx.x,gy.x);
    vec2 g1 = vec2(gx.y,gy.y);
    vec2 g2 = vec2(gx.z,gy.z);

    // Radial decay with distance from each simplex corner
    vec3 w = 0.8 - vec3(dot(x0, x0), dot(x1, x1), dot(x2, x2));
    w = max(w, 0.0);
    vec3 w2 = w * w;
    vec3 w4 = w2 * w2;

    // The value of the linear ramp from each of the corners
    vec3 gdotx = vec3(dot(g0, x0), dot(g1, x1), dot(g2, x2));

    // Multiply by the radial decay and sum up the noise value
    float n = dot(w4, gdotx);

    // Compute the first order partial derivatives
    vec3 w3 = w2 * w;
    vec3 dw = -8.0 * w3 * gdotx;
    vec2 dn0 = w4.x * g0 + dw.x * x0;
    vec2 dn1 = w4.y * g1 + dw.y * x1;
    vec2 dn2 = w4.z * g2 + dw.z * x2;
    gradient = 10.9 * (dn0 + dn1 + dn2);

    // Scale the return value to fit nicely into the range [-1,1]
    return 10.9 * n;
}

float psrdnoise(vec2 x, vec2 period, float alpha, out vec2 gradient, out vec3 dg) {

    // Transform to simplex space (axis-aligned hexagonal grid)
    vec2 uv = vec2(x.x + x.y*0.5, x.y);

    // Determine which simplex we're in, with i0 being the "base"
    vec2 i0 = floor(uv);
    vec2 f0 = fract(uv);
    // o1 is the offset in simplex space to the second corner
    float cmp = step(f0.y, f0.x);
    vec2 o1 = vec2(cmp, 1.0-cmp);

    // Enumerate the remaining simplex corners
    vec2 i1 = i0 + o1;
    vec2 i2 = i0 + vec2(1.0, 1.0);

    // Transform corners back to texture space
    vec2 v0 = vec2(i0.x - i0.y * 0.5, i0.y);
    vec2 v1 = vec2(v0.x + o1.x - o1.y * 0.5, v0.y + o1.y);
    vec2 v2 = vec2(v0.x + 0.5, v0.y + 1.0);

    // Compute vectors from v to each of the simplex corners
    vec2 x0 = x - v0;
    vec2 x1 = x - v1;
    vec2 x2 = x - v2;

    vec3 iu, iv;
    vec3 xw, yw;

    // Wrap to periods, if desired
    if(any(greaterThan(period, vec2(0.0)))) {
        xw = vec3(v0.x, v1.x, v2.x);
        yw = vec3(v0.y, v1.y, v2.y);
        if(period.x > 0.0)
            xw = mod(vec3(v0.x, v1.x, v2.x), period.x);
        if(period.y > 0.0)
            yw = mod(vec3(v0.y, v1.y, v2.y), period.y);
        // Transform back to simplex space and fix rounding errors
        iu = floor(xw + 0.5*yw + 0.5);
        iv = floor(yw + 0.5);
    } else { // Shortcut if neither x nor y periods are specified
        iu = vec3(i0.x, i1.x, i2.x);
        iv = vec3(i0.y, i1.y, i2.y);
    }

    // Compute one pseudo-random hash value for each corner
    vec3 hash = mod(iu, 289.0);
    hash = mod((hash*51.0 + 2.0)*hash + iv, 289.0);
    hash = mod((hash*34.0 + 10.0)*hash, 289.0);

    // Pick a pseudo-random angle and add the desired rotation
    vec3 psi = hash * 0.07482 + alpha;
    vec3 gx = cos(psi);
    vec3 gy = sin(psi);

    // Reorganize for dot products below
    vec2 g0 = vec2(gx.x,gy.x);
    vec2 g1 = vec2(gx.y,gy.y);
    vec2 g2 = vec2(gx.z,gy.z);

    // Radial decay with distance from each simplex corner
    vec3 w = 0.8 - vec3(dot(x0, x0), dot(x1, x1), dot(x2, x2));
    w = max(w, 0.0);
    vec3 w2 = w * w;
    vec3 w4 = w2 * w2;

    // The value of the linear ramp from each of the corners
    vec3 gdotx = vec3(dot(g0, x0), dot(g1, x1), dot(g2, x2));

    // Multiply by the radial decay and sum up the noise value
    float n = dot(w4, gdotx);

    // Compute the first order partial derivatives
    vec3 w3 = w2 * w;
    vec3 dw = -8.0 * w3 * gdotx;
    vec2 dn0 = w4.x * g0 + dw.x * x0;
    vec2 dn1 = w4.y * g1 + dw.y * x1;
    vec2 dn2 = w4.z * g2 + dw.z * x2;
    gradient = 10.9 * (dn0 + dn1 + dn2);

    // Compute the second order partial derivatives
    vec3 dg0, dg1, dg2;
    vec3 dw2 = 48.0 * w2 * gdotx;
    // d2n/dx2 and d2n/dy2
    dg0.xy = dw2.x * x0 * x0 - 8.0 * w3.x * (2.0 * g0 * x0 + gdotx.x);
    dg1.xy = dw2.y * x1 * x1 - 8.0 * w3.y * (2.0 * g1 * x1 + gdotx.y);
    dg2.xy = dw2.z * x2 * x2 - 8.0 * w3.z * (2.0 * g2 * x2 + gdotx.z);
    // d2n/dxy
    dg0.z = dw2.x * x0.x * x0.y - 8.0 * w3.x * dot(g0, x0.yx);
    dg1.z = dw2.y * x1.x * x1.y - 8.0 * w3.y * dot(g1, x1.yx);
    dg2.z = dw2.z * x2.x * x2.y - 8.0 * w3.z * dot(g2, x2.yx);
    dg = 10.9 * (dg0 + dg1 + dg2);

    // Scale the return value to fit nicely into the range [-1,1]
    return 10.9 * n;
}

float psrdnoise(vec2 x, vec2 period, float alpha) {
    vec2 g = vec2(0.0);
    return psrdnoise(x, period, alpha, g);
}

float psrdnoise(vec2 x, vec2 period) {
    return psrdnoise(x, period, 0.0);
}

float psrdnoise(vec2 x) {
    return psrdnoise(x, vec2(0.0));
}

float psrdnoise(vec3 x, vec3 period, float alpha, out vec3 gradient) {

#ifndef PSRDNOISE_PERLIN_GRID
    // Transformation matrices for the axis-aligned simplex grid
    const mat3 M = mat3(0.0, 1.0, 1.0,
                        1.0, 0.0, 1.0,
                        1.0, 1.0, 0.0);

    const mat3 Mi = mat3(-0.5, 0.5, 0.5,
                            0.5,-0.5, 0.5,
                            0.5, 0.5,-0.5);
#endif

    vec3 uvw = vec3(0.0);

    // Transform to simplex space (tetrahedral grid)
#ifndef PSRDNOISE_PERLIN_GRID
    // Use matrix multiplication, let the compiler optimise
    uvw = M * x;
 #else
    // Optimised transformation to uvw (slightly faster than
    // the equivalent matrix multiplication on most platforms)
    uvw = x + dot(x, vec3(0.33333333));
 #endif

    // Determine which simplex we're in, i0 is the "base corner"
    vec3 i0 = floor(uvw);
    vec3 f0 = fract(uvw); // coords within "skewed cube"

    // To determine which simplex corners are closest, rank order the
    // magnitudes of u,v,w, resolving ties in priority order u,v,w,
    // and traverse the four corners from largest to smallest magnitude.
    // o1, o2 are offsets in simplex space to the 2nd and 3rd corners.
    vec3 g_ = step(f0.xyx, f0.yzz); // Makes comparison "less-than"
    vec3 l_ = 1.0 - g_;             // complement is "greater-or-equal"
    vec3 g = vec3(l_.z, g_.xy);
    vec3 l = vec3(l_.xy, g_.z);
    vec3 o1 = min( g, l );
    vec3 o2 = max( g, l );

    // Enumerate the remaining simplex corners
    vec3 i1 = i0 + o1;
    vec3 i2 = i0 + o2;
    vec3 i3 = i0 + vec3(1.0);

    vec3 v0 = vec3(0.0);
    vec3 v1 = vec3(0.0);
    vec3 v2 = vec3(0.0);
    vec3 v3 = vec3(0.0);

    // Transform the corners back to texture space
#ifndef PSRDNOISE_PERLIN_GRID
    v0 = Mi * i0;
    v1 = Mi * i1;
    v2 = Mi * i2;
    v3 = Mi * i3;
#else
    // Optimised transformation (mostly slightly faster than a matrix)
    v0 = i0 - dot(i0, vec3(1.0/6.0));
    v1 = i1 - dot(i1, vec3(1.0/6.0));
    v2 = i2 - dot(i2, vec3(1.0/6.0));
    v3 = i3 - dot(i3, vec3(1.0/6.0));
#endif

    // Compute vectors to each of the simplex corners
    vec3 x0 = x - v0;
    vec3 x1 = x - v1;
    vec3 x2 = x - v2;
    vec3 x3 = x - v3;

    if(any(greaterThan(period, vec3(0.0)))) {
        // Wrap to periods and transform back to simplex space
        vec4 vx = vec4(v0.x, v1.x, v2.x, v3.x);
        vec4 vy = vec4(v0.y, v1.y, v2.y, v3.y);
        vec4 vz = vec4(v0.z, v1.z, v2.z, v3.z);
        // Wrap to periods where specified
        if(period.x > 0.0) vx = mod(vx, period.x);
        if(period.y > 0.0) vy = mod(vy, period.y);
        if(period.z > 0.0) vz = mod(vz, period.z);
        // Transform back
#ifndef PSRDNOISE_PERLIN_GRID
        i0 = M * vec3(vx.x, vy.x, vz.x);
        i1 = M * vec3(vx.y, vy.y, vz.y);
        i2 = M * vec3(vx.z, vy.z, vz.z);
        i3 = M * vec3(vx.w, vy.w, vz.w);
#else
        v0 = vec3(vx.x, vy.x, vz.x);
        v1 = vec3(vx.y, vy.y, vz.y);
        v2 = vec3(vx.z, vy.z, vz.z);
        v3 = vec3(vx.w, vy.w, vz.w);
        // Transform wrapped coordinates back to uvw
        i0 = v0 + dot(v0, vec3(1.0/3.0));
        i1 = v1 + dot(v1, vec3(1.0/3.0));
        i2 = v2 + dot(v2, vec3(1.0/3.0));
        i3 = v3 + dot(v3, vec3(1.0/3.0));
#endif
        // Fix rounding errors
        i0 = floor(i0 + 0.5);
        i1 = floor(i1 + 0.5);
        i2 = floor(i2 + 0.5);
        i3 = floor(i3 + 0.5);
    }

    // Compute one pseudo-random hash value for each corner
    vec4 hash = permute( permute( permute( 
                vec4(i0.z, i1.z, i2.z, i3.z ))
                + vec4(i0.y, i1.y, i2.y, i3.y ))
                + vec4(i0.x, i1.x, i2.x, i3.x ));

    // Compute generating gradients from a Fibonacci spiral on the unit sphere
    vec4 theta = hash * 3.883222077;  // 2*pi/golden ratio
    vec4 sz    = hash * -0.006920415 + 0.996539792; // 1-(hash+0.5)*2/289
    vec4 psi   = hash * 0.108705628 ; // 10*pi/289, chosen to avoid correlation

    vec4 Ct = cos(theta);
    vec4 St = sin(theta);
    vec4 sz_prime = sqrt( 1.0 - sz*sz ); // s is a point on a unit fib-sphere

    vec4 gx = vec4(0.0);
    vec4 gy = vec4(0.0);
    vec4 gz = vec4(0.0);

    // Rotate gradients by angle alpha around a pseudo-random ortogonal axis
#ifdef PSRDNOISE_FAST_ROTATION
    // Fast algorithm, but without dynamic shortcut for alpha = 0
    vec4 qx = St;         // q' = norm ( cross(s, n) )  on the equator
    vec4 qy = -Ct; 
    vec4 qz = vec4(0.0);

    vec4 px =  sz * qy;   // p' = cross(q, s)
    vec4 py = -sz * qx;
    vec4 pz = sz_prime;

    psi += alpha;         // psi and alpha in the same plane
    vec4 Sa = sin(psi);
    vec4 Ca = cos(psi);

    gx = Ca * px + Sa * qx;
    gy = Ca * py + Sa * qy;
    gz = Ca * pz + Sa * qz;
#else
    // Slightly slower algorithm, but with g = s for alpha = 0, and a
    // useful conditional speedup for alpha = 0 across all fragments
    if(alpha != 0.0) {
        vec4 Sp = sin(psi);          // q' from psi on equator
        vec4 Cp = cos(psi);

        vec4 px = Ct * sz_prime;     // px = sx
        vec4 py = St * sz_prime;     // py = sy
        vec4 pz = sz;

        vec4 Ctp = St*Sp - Ct*Cp;    // q = (rotate( cross(s,n), dot(s,n))(q')
        vec4 qx = mix( Ctp*St, Sp, sz);
        vec4 qy = mix(-Ctp*Ct, Cp, sz);
        vec4 qz = -(py*Cp + px*Sp);

        vec4 Sa = vec4(sin(alpha));       // psi and alpha in different planes
        vec4 Ca = vec4(cos(alpha));

        gx = Ca * px + Sa * qx;
        gy = Ca * py + Sa * qy;
        gz = Ca * pz + Sa * qz;
    }
    else {
        gx = Ct * sz_prime;  // alpha = 0, use s directly as gradient
        gy = St * sz_prime;
        gz = sz;  
    }
#endif

    // Reorganize for dot products below
    vec3 g0 = vec3(gx.x, gy.x, gz.x);
    vec3 g1 = vec3(gx.y, gy.y, gz.y);
    vec3 g2 = vec3(gx.z, gy.z, gz.z);
    vec3 g3 = vec3(gx.w, gy.w, gz.w);

    // Radial decay with distance from each simplex corner
    vec4 w = 0.5 - vec4(dot(x0,x0), dot(x1,x1), dot(x2,x2), dot(x3,x3));
    w = max(w, 0.0);
    vec4 w2 = w * w;
    vec4 w3 = w2 * w;

    // The value of the linear ramp from each of the corners
    vec4 gdotx = vec4(dot(g0,x0), dot(g1,x1), dot(g2,x2), dot(g3,x3));

    // Multiply by the radial decay and sum up the noise value
    float n = dot(w3, gdotx);

    // Compute the first order partial derivatives
    vec4 dw = -6.0 * w2 * gdotx;
    vec3 dn0 = w3.x * g0 + dw.x * x0;
    vec3 dn1 = w3.y * g1 + dw.y * x1;
    vec3 dn2 = w3.z * g2 + dw.z * x2;
    vec3 dn3 = w3.w * g3 + dw.w * x3;
    gradient = 39.5 * (dn0 + dn1 + dn2 + dn3);

    // Scale the return value to fit nicely into the range [-1,1]
    return 39.5 * n; 
}

float psrdnoise(vec3 x, vec3 period, float alpha, out vec3 gradient, out vec3 dg, out vec3 dg2) {

#ifndef PSRDNOISE_PERLIN_GRID
    // Transformation matrices for the axis-aligned simplex grid
    const mat3 M = mat3(0.0, 1.0, 1.0,
                        1.0, 0.0, 1.0,
                        1.0, 1.0, 0.0);

    const mat3 Mi = mat3(-0.5, 0.5, 0.5,
                            0.5,-0.5, 0.5,
                            0.5, 0.5,-0.5);
#endif

    vec3 uvw = vec3(0.0);

    // Transform to simplex space (tetrahedral grid)
#ifndef PSRDNOISE_PERLIN_GRID
    // Use matrix multiplication, let the compiler optimise
    uvw = M * x;
#else
    // Optimised transformation to uvw (slightly faster than
    // the equivalent matrix multiplication on most platforms)
    uvw = x + dot(x, vec3(0.3333333));
#endif

    // Determine which simplex we're in, i0 is the "base corner"
    vec3 i0 = floor(uvw);
    vec3 f0 = fract(uvw); // coords within "skewed cube"

    // To determine which simplex corners are closest, rank order the
    // magnitudes of u,v,w, resolving ties in priority order u,v,w,
    // and traverse the four corners from largest to smallest magnitude.
    // o1, o2 are offsets in simplex space to the 2nd and 3rd corners.
    vec3 g_ = step(f0.xyx, f0.yzz); // Makes comparison "less-than"
    vec3 l_ = 1.0 - g_;             // complement is "greater-or-equal"
    vec3 g = vec3(l_.z, g_.xy);
    vec3 l = vec3(l_.xy, g_.z);
    vec3 o1 = min( g, l );
    vec3 o2 = max( g, l );

    // Enumerate the remaining simplex corners
    vec3 i1 = i0 + o1;
    vec3 i2 = i0 + o2;
    vec3 i3 = i0 + vec3(1.0);

    vec3 v0, v1, v2, v3;

    // Transform the corners back to texture space
#ifndef PSRDNOISE_PERLIN_GRID
    v0 = Mi * i0;
    v1 = Mi * i1;
    v2 = Mi * i2;
    v3 = Mi * i3;
#else
    // Optimised transformation (mostly slightly faster than a matrix)
    v0 = i0 - dot(i0, vec3(1.0/6.0));
    v1 = i1 - dot(i1, vec3(1.0/6.0));
    v2 = i2 - dot(i2, vec3(1.0/6.0));
    v3 = i3 - dot(i3, vec3(1.0/6.0));
#endif

    // Compute vectors to each of the simplex corners
    vec3 x0 = x - v0;
    vec3 x1 = x - v1;
    vec3 x2 = x - v2;
    vec3 x3 = x - v3;

    if(any(greaterThan(period, vec3(0.0)))) {
        // Wrap to periods and transform back to simplex space
        vec4 vx = vec4(v0.x, v1.x, v2.x, v3.x);
        vec4 vy = vec4(v0.y, v1.y, v2.y, v3.y);
        vec4 vz = vec4(v0.z, v1.z, v2.z, v3.z);
        // Wrap to periods where specified
        if(period.x > 0.0) vx = mod(vx, period.x);
        if(period.y > 0.0) vy = mod(vy, period.y);
        if(period.z > 0.0) vz = mod(vz, period.z);
        // Transform back
#ifndef PSRDNOISE_PERLIN_GRID
        i0 = M * vec3(vx.x, vy.x, vz.x);
        i1 = M * vec3(vx.y, vy.y, vz.y);
        i2 = M * vec3(vx.z, vy.z, vz.z);
        i3 = M * vec3(vx.w, vy.w, vz.w);
#else
        v0 = vec3(vx.x, vy.x, vz.x);
        v1 = vec3(vx.y, vy.y, vz.y);
        v2 = vec3(vx.z, vy.z, vz.z);
        v3 = vec3(vx.w, vy.w, vz.w);
        // Transform wrapped coordinates back to uvw
        i0 = v0 + dot(v0, vec3(0.3333333));
        i1 = v1 + dot(v1, vec3(0.3333333));
        i2 = v2 + dot(v2, vec3(0.3333333));
        i3 = v3 + dot(v3, vec3(0.3333333));
#endif
        // Fix rounding errors
        i0 = floor(i0 + 0.5);
        i1 = floor(i1 + 0.5);
        i2 = floor(i2 + 0.5);
        i3 = floor(i3 + 0.5);
    }

    // Compute one pseudo-random hash value for each corner
    vec4 hash = permute( permute( permute( 
                vec4(i0.z, i1.z, i2.z, i3.z ))
                + vec4(i0.y, i1.y, i2.y, i3.y ))
                + vec4(i0.x, i1.x, i2.x, i3.x ));

    // Compute generating gradients from a Fibonacci spiral on the unit sphere
    vec4 theta = hash * 3.883222077;  // 2*pi/golden ratio
    vec4 sz    = hash * -0.006920415 + 0.996539792; // 1-(hash+0.5)*2/289
    vec4 psi   = hash * 0.108705628 ; // 10*pi/289, chosen to avoid correlation

    vec4 Ct = cos(theta);
    vec4 St = sin(theta);
    vec4 sz_prime = sqrt( 1.0 - sz*sz ); // s is a point on a unit fib-sphere

    vec4 gx, gy, gz;

    // Rotate gradients by angle alpha around a pseudo-random ortogonal axis
#ifdef PSRDNOISE_FAST_ROTATION
    // Fast algorithm, but without dynamic shortcut for alpha = 0
    vec4 qx = St;         // q' = norm ( cross(s, n) )  on the equator
    vec4 qy = -Ct; 
    vec4 qz = vec4(0.0);

    vec4 px =  sz * qy;   // p' = cross(q, s)
    vec4 py = -sz * qx;
    vec4 pz = sz_prime;

    psi += alpha;         // psi and alpha in the same plane
    vec4 Sa = sin(psi);
    vec4 Ca = cos(psi);

    gx = Ca * px + Sa * qx;
    gy = Ca * py + Sa * qy;
    gz = Ca * pz + Sa * qz;
    #else
    // Slightly slower algorithm, but with g = s for alpha = 0, and a
    // strong conditional speedup for alpha = 0 across all fragments
    if(alpha != 0.0) {
        vec4 Sp = sin(psi);          // q' from psi on equator
        vec4 Cp = cos(psi);

        vec4 px = Ct * sz_prime;     // px = sx
        vec4 py = St * sz_prime;     // py = sy
        vec4 pz = sz;

        vec4 Ctp = St*Sp - Ct*Cp;    // q = (rotate( cross(s,n), dot(s,n))(q')
        vec4 qx = mix( Ctp*St, Sp, sz);
        vec4 qy = mix(-Ctp*Ct, Cp, sz);
        vec4 qz = -(py*Cp + px*Sp);

        vec4 Sa = vec4(sin(alpha));       // psi and alpha in different planes
        vec4 Ca = vec4(cos(alpha));

        gx = Ca * px + Sa * qx;
        gy = Ca * py + Sa * qy;
        gz = Ca * pz + Sa * qz;
    }
    else {
        gx = Ct * sz_prime;  // alpha = 0, use s directly as gradient
        gy = St * sz_prime;
        gz = sz;  
    }
#endif

    // Reorganize for dot products below
    vec3 g0 = vec3(gx.x, gy.x, gz.x);
    vec3 g1 = vec3(gx.y, gy.y, gz.y);
    vec3 g2 = vec3(gx.z, gy.z, gz.z);
    vec3 g3 = vec3(gx.w, gy.w, gz.w);

    // Radial decay with distance from each simplex corner
    vec4 w = 0.5 - vec4(dot(x0,x0), dot(x1,x1), dot(x2,x2), dot(x3,x3));
    w = max(w, 0.0);
    vec4 w2 = w * w;
    vec4 w3 = w2 * w;

    // The value of the linear ramp from each of the corners
    vec4 gdotx = vec4(dot(g0,x0), dot(g1,x1), dot(g2,x2), dot(g3,x3));

    // Multiply by the radial decay and sum up the noise value
    float n = dot(w3, gdotx);

    // Compute the first order partial derivatives
    vec4 dw = -6.0 * w2 * gdotx;
    vec3 dn0 = w3.x * g0 + dw.x * x0;
    vec3 dn1 = w3.y * g1 + dw.y * x1;
    vec3 dn2 = w3.z * g2 + dw.z * x2;
    vec3 dn3 = w3.w * g3 + dw.w * x3;
    gradient = 39.5 * (dn0 + dn1 + dn2 + dn3);

    // Compute the second order partial derivatives
    vec4 dw2 = 24.0 * w * gdotx;
    vec3 dga0 = dw2.x * x0 * x0 - 6.0 * w2.x * (gdotx.x + 2.0 * g0 * x0);
    vec3 dga1 = dw2.y * x1 * x1 - 6.0 * w2.y * (gdotx.y + 2.0 * g1 * x1);
    vec3 dga2 = dw2.z * x2 * x2 - 6.0 * w2.z * (gdotx.z + 2.0 * g2 * x2);
    vec3 dga3 = dw2.w * x3 * x3 - 6.0 * w2.w * (gdotx.w + 2.0 * g3 * x3);
    dg = 35.0 * (dga0 + dga1 + dga2 + dga3); // (d2n/dx2, d2n/dy2, d2n/dz2)
    vec3 dgb0 = dw2.x * x0 * x0.yzx - 6.0 * w2.x * (g0 * x0.yzx + g0.yzx * x0);
    vec3 dgb1 = dw2.y * x1 * x1.yzx - 6.0 * w2.y * (g1 * x1.yzx + g1.yzx * x1);
    vec3 dgb2 = dw2.z * x2 * x2.yzx - 6.0 * w2.z * (g2 * x2.yzx + g2.yzx * x2);
    vec3 dgb3 = dw2.w * x3 * x3.yzx - 6.0 * w2.w * (g3 * x3.yzx + g3.yzx * x3);
    dg2 = 39.5 * (dgb0 + dgb1 + dgb2 + dgb3); // (d2n/dxy, d2n/dyz, d2n/dxz)

    // Scale the return value to fit nicely into the range [-1,1]
    return 39.5 * n;
}

float psrdnoise(vec3 x, vec3 period, float alpha) {
    vec3 g = vec3(0.0);
    return psrdnoise(x, period, alpha, g);
}

float psrdnoise(vec3 x, vec3 period) {
    return psrdnoise(x, period, 0.0);
}

float psrdnoise(vec3 x) {
    return psrdnoise(x, vec3(0.0));
}
#endif

Dependencies:

Use:

<float> psrdnoise(<float2|float3> x, <float2|float3> period, <float> alpha, out <float2|float3> gradient [, out <float3> dg])

Check it on Github



#ifndef FNC_PSRFNOISE
#define FNC_PSRFNOISE

float psrdnoise(float2 x, float2 period, float alpha, out float2 gradient) {

    // Transform to simplex space (axis-aligned hexagonal grid)
    float2 uv = float2(x.x + x.y*0.5, x.y);

    // Determine which simplex we're in, with i0 being the "base"
    float2 i0 = floor(uv);
    float2 f0 = frac(uv);
    // o1 is the offset in simplex space to the second corner
    float cmp = step(f0.y, f0.x);
    float2 o1 = float2(cmp, 1.0-cmp);

    // Enumerate the remaining simplex corners
    float2 i1 = i0 + o1;
    float2 i2 = i0 + float2(1.0, 1.0);

    // Transform corners back to texture space
    float2 v0 = float2(i0.x - i0.y * 0.5, i0.y);
    float2 v1 = float2(v0.x + o1.x - o1.y * 0.5, v0.y + o1.y);
    float2 v2 = float2(v0.x + 0.5, v0.y + 1.0);

    // Compute vectors from v to each of the simplex corners
    float2 x0 = x - v0;
    float2 x1 = x - v1;
    float2 x2 = x - v2;

    float3 iu = float3(0.0, 0.0, 0.0);
    float3 iv = float3(0.0, 0.0, 0.0);
    float3 xw = float3(0.0, 0.0, 0.0);
    float3 yw = float3(0.0, 0.0, 0.0);

    // Wrap to periods, if desired
    if( any(greaterThan(period, float2(0.0, 0.0))) ) {
        xw = float3(v0.x, v1.x, v2.x);
        yw = float3(v0.y, v1.y, v2.y);
        if(period.x > 0.0)
            xw = mod(float3(v0.x, v1.x, v2.x), period.x);
        if(period.y > 0.0)
            yw = mod(float3(v0.y, v1.y, v2.y), period.y);
        // Transform back to simplex space and fix rounding errors
        iu = floor(xw + 0.5*yw + 0.5);
        iv = floor(yw + 0.5);
    } else { // Shortcut if neither x nor y periods are specified
        iu = float3(i0.x, i1.x, i2.x);
        iv = float3(i0.y, i1.y, i2.y);
    }

    // Compute one pseudo-random hash value for each corner
    float3 hash = mod(iu, 289.0);
    hash = mod((hash*51.0 + 2.0)*hash + iv, 289.0);
    hash = mod((hash*34.0 + 10.0)*hash, 289.0);

    // Pick a pseudo-random angle and add the desired rotation
    float3 psi = hash * 0.07482 + alpha;
    float3 gx = cos(psi);
    float3 gy = sin(psi);

    // Reorganize for dot products below
    float2 g0 = float2(gx.x,gy.x);
    float2 g1 = float2(gx.y,gy.y);
    float2 g2 = float2(gx.z,gy.z);

    // Radial decay with distance from each simplex corner
    float3 w = 0.8 - float3(dot(x0, x0), dot(x1, x1), dot(x2, x2));
    w = max(w, 0.0);
    float3 w2 = w * w;
    float3 w4 = w2 * w2;

    // The value of the linear ramp from each of the corners
    float3 gdotx = float3(dot(g0, x0), dot(g1, x1), dot(g2, x2));

    // Multiply by the radial decay and sum up the noise value
    float n = dot(w4, gdotx);

    // Compute the first order partial derivatives
    float3 w3 = w2 * w;
    float3 dw = -8.0 * w3 * gdotx;
    float2 dn0 = w4.x * g0 + dw.x * x0;
    float2 dn1 = w4.y * g1 + dw.y * x1;
    float2 dn2 = w4.z * g2 + dw.z * x2;
    gradient = 10.9 * (dn0 + dn1 + dn2);

    // Scale the return value to fit nicely into the range [-1,1]
    return 10.9 * n;
}

float psrdnoise(float2 x, float2 period, float alpha, out float2 gradient, out float3 dg) {

    // Transform to simplex space (axis-aligned hexagonal grid)
    float2 uv = float2(x.x + x.y*0.5, x.y);

    // Determine which simplex we're in, with i0 being the "base"
    float2 i0 = floor(uv);
    float2 f0 = frac(uv);
    // o1 is the offset in simplex space to the second corner
    float cmp = step(f0.y, f0.x);
    float2 o1 = float2(cmp, 1.0-cmp);

    // Enumerate the remaining simplex corners
    float2 i1 = i0 + o1;
    float2 i2 = i0 + float2(1.0, 1.0);

    // Transform corners back to texture space
    float2 v0 = float2(i0.x - i0.y * 0.5, i0.y);
    float2 v1 = float2(v0.x + o1.x - o1.y * 0.5, v0.y + o1.y);
    float2 v2 = float2(v0.x + 0.5, v0.y + 1.0);

    // Compute vectors from v to each of the simplex corners
    float2 x0 = x - v0;
    float2 x1 = x - v1;
    float2 x2 = x - v2;

    float3 iu, iv;
    float3 xw, yw;

    // Wrap to periods, if desired
    if(any(greaterThan(period, float2(0.0, 0.0)))) {
        xw = float3(v0.x, v1.x, v2.x);
        yw = float3(v0.y, v1.y, v2.y);
        if(period.x > 0.0)
            xw = mod(float3(v0.x, v1.x, v2.x), period.x);
        if(period.y > 0.0)
            yw = mod(float3(v0.y, v1.y, v2.y), period.y);
        // Transform back to simplex space and fix rounding errors
        iu = floor(xw + 0.5*yw + 0.5);
        iv = floor(yw + 0.5);
    } else { // Shortcut if neither x nor y periods are specified
        iu = float3(i0.x, i1.x, i2.x);
        iv = float3(i0.y, i1.y, i2.y);
    }

    // Compute one pseudo-random hash value for each corner
    float3 hash = mod(iu, 289.0);
    hash = mod((hash*51.0 + 2.0)*hash + iv, 289.0);
    hash = mod((hash*34.0 + 10.0)*hash, 289.0);

    // Pick a pseudo-random angle and add the desired rotation
    float3 psi = hash * 0.07482 + alpha;
    float3 gx = cos(psi);
    float3 gy = sin(psi);

    // Reorganize for dot products below
    float2 g0 = float2(gx.x,gy.x);
    float2 g1 = float2(gx.y,gy.y);
    float2 g2 = float2(gx.z,gy.z);

    // Radial decay with distance from each simplex corner
    float3 w = 0.8 - float3(dot(x0, x0), dot(x1, x1), dot(x2, x2));
    w = max(w, 0.0);
    float3 w2 = w * w;
    float3 w4 = w2 * w2;

    // The value of the linear ramp from each of the corners
    float3 gdotx = float3(dot(g0, x0), dot(g1, x1), dot(g2, x2));

    // Multiply by the radial decay and sum up the noise value
    float n = dot(w4, gdotx);

    // Compute the first order partial derivatives
    float3 w3 = w2 * w;
    float3 dw = -8.0 * w3 * gdotx;
    float2 dn0 = w4.x * g0 + dw.x * x0;
    float2 dn1 = w4.y * g1 + dw.y * x1;
    float2 dn2 = w4.z * g2 + dw.z * x2;
    gradient = 10.9 * (dn0 + dn1 + dn2);

    // Compute the second order partial derivatives
    float3 dg0, dg1, dg2;
    float3 dw2 = 48.0 * w2 * gdotx;
    // d2n/dx2 and d2n/dy2
    dg0.xy = dw2.x * x0 * x0 - 8.0 * w3.x * (2.0 * g0 * x0 + gdotx.x);
    dg1.xy = dw2.y * x1 * x1 - 8.0 * w3.y * (2.0 * g1 * x1 + gdotx.y);
    dg2.xy = dw2.z * x2 * x2 - 8.0 * w3.z * (2.0 * g2 * x2 + gdotx.z);
    // d2n/dxy
    dg0.z = dw2.x * x0.x * x0.y - 8.0 * w3.x * dot(g0, x0.yx);
    dg1.z = dw2.y * x1.x * x1.y - 8.0 * w3.y * dot(g1, x1.yx);
    dg2.z = dw2.z * x2.x * x2.y - 8.0 * w3.z * dot(g2, x2.yx);
    dg = 10.9 * (dg0 + dg1 + dg2);

    // Scale the return value to fit nicely into the range [-1,1]
    return 10.9 * n;
}

float psrdnoise(float2 x, float2 period, float alpha) {
    float2 g = float2(0.0, 0.0);
    return psrdnoise(x, period, alpha, g);
}

float psrdnoise(float2 x, float2 period) {
    return psrdnoise(x, period, 0.0);
}

float psrdnoise(float2 x) {
    return psrdnoise(x, float2(0.0, 0.0));
}

float psrdnoise(float3 x, float3 period, float alpha, out float3 gradient) {

#ifndef PSRDNOISE_PERLIN_GRID
    // Transformation matrices for the axis-aligned simplex grid
    const float3x3 M = float3x3(0.0, 1.0, 1.0,
                                1.0, 0.0, 1.0,
                                1.0, 1.0, 0.0);

    const float3x3 Mi = float3x3(-0.5, 0.5, 0.5,
                                0.5,-0.5, 0.5,
                                0.5, 0.5,-0.5);
#endif

    float3 uvw = float3(0.0, 0.0, 0.0);

    // Transform to simplex space (tetrahedral grid)
#ifndef PSRDNOISE_PERLIN_GRID
    // Use matrix multiplication, let the compiler optimise
    uvw = mul(M, x);
 #else
    // Optimised transformation to uvw (slightly faster than
    // the equivalent matrix multiplication on most platforms)
    uvw = x + dot(x, float3(0.33333333, 0.33333333, 0.33333333));
 #endif

    // Determine which simplex we're in, i0 is the "base corner"
    float3 i0 = floor(uvw);
    float3 f0 = frac(uvw); // coords within "skewed cube"

    // To determine which simplex corners are closest, rank order the
    // magnitudes of u,v,w, resolving ties in priority order u,v,w,
    // and traverse the four corners from largest to smallest magnitude.
    // o1, o2 are offsets in simplex space to the 2nd and 3rd corners.
    float3 g_ = step(f0.xyx, f0.yzz); // Makes comparison "less-than"
    float3 l_ = 1.0 - g_;             // complement is "greater-or-equal"
    float3 g = float3(l_.z, g_.xy);
    float3 l = float3(l_.xy, g_.z);
    float3 o1 = min( g, l );
    float3 o2 = max( g, l );

    // Enumerate the remaining simplex corners
    float3 i1 = i0 + o1;
    float3 i2 = i0 + o2;
    float3 i3 = i0 + float3(1.0, 1.0, 1.0);

    float3 v0 = float3(0.0, 0.0, 0.0);
    float3 v1 = float3(0.0, 0.0, 0.0);
    float3 v2 = float3(0.0, 0.0, 0.0);
    float3 v3 = float3(0.0, 0.0, 0.0);

    // Transform the corners back to texture space
#ifndef PSRDNOISE_PERLIN_GRID
    v0 = mul(Mi, i0);
    v1 = mul(Mi, i1);
    v2 = mul(Mi, i2);
    v3 = mul(Mi, i3);
#else
    // Optimised transformation (mostly slightly faster than a matrix)
    v0 = i0 - dot(i0, float3(1.0/6.0));
    v1 = i1 - dot(i1, float3(1.0/6.0));
    v2 = i2 - dot(i2, float3(1.0/6.0));
    v3 = i3 - dot(i3, float3(1.0/6.0));
#endif

    // Compute vectors to each of the simplex corners
    float3 x0 = x - v0;
    float3 x1 = x - v1;
    float3 x2 = x - v2;
    float3 x3 = x - v3;

    if(any(greaterThan(period, float3(0.0, 0.0, 0.0)))) {
        // Wrap to periods and transform back to simplex space
        float4 vx = float4(v0.x, v1.x, v2.x, v3.x);
        float4 vy = float4(v0.y, v1.y, v2.y, v3.y);
        float4 vz = float4(v0.z, v1.z, v2.z, v3.z);
        // Wrap to periods where specified
        if(period.x > 0.0) vx = mod(vx, period.x);
        if(period.y > 0.0) vy = mod(vy, period.y);
        if(period.z > 0.0) vz = mod(vz, period.z);
        // Transform back
#ifndef PSRDNOISE_PERLIN_GRID
        i0 = mul(M, float3(vx.x, vy.x, vz.x));
        i1 = mul(M, float3(vx.y, vy.y, vz.y));
        i2 = mul(M, float3(vx.z, vy.z, vz.z));
        i3 = mul(M, float3(vx.w, vy.w, vz.w));
#else
        v0 = float3(vx.x, vy.x, vz.x);
        v1 = float3(vx.y, vy.y, vz.y);
        v2 = float3(vx.z, vy.z, vz.z);
        v3 = float3(vx.w, vy.w, vz.w);
        // Transform wrapped coordinates back to uvw
        i0 = v0 + dot(v0, float3(1.0/3.0));
        i1 = v1 + dot(v1, float3(1.0/3.0));
        i2 = v2 + dot(v2, float3(1.0/3.0));
        i3 = v3 + dot(v3, float3(1.0/3.0));
#endif
        // Fix rounding errors
        i0 = floor(i0 + 0.5);
        i1 = floor(i1 + 0.5);
        i2 = floor(i2 + 0.5);
        i3 = floor(i3 + 0.5);
    }

    // Compute one pseudo-random hash value for each corner
    float4 hash = permute( permute( permute( 
                float4(i0.z, i1.z, i2.z, i3.z ))
                + float4(i0.y, i1.y, i2.y, i3.y ))
                + float4(i0.x, i1.x, i2.x, i3.x ));

    // Compute generating gradients from a Fibonacci spiral on the unit sphere
    float4 theta = hash * 3.883222077;  // 2*pi/golden ratio
    float4 sz    = hash * -0.006920415 + 0.996539792; // 1-(hash+0.5)*2/289
    float4 psi   = hash * 0.108705628 ; // 10*pi/289, chosen to avoid correlation

    float4 Ct = cos(theta);
    float4 St = sin(theta);
    float4 sz_prime = sqrt( 1.0 - sz*sz ); // s is a point on a unit fib-sphere

    float4 gx = float4(0.0, 0.0, 0.0, 0.0);
    float4 gy = float4(0.0, 0.0, 0.0, 0.0);
    float4 gz = float4(0.0, 0.0, 0.0, 0.0);

    // Rotate gradients by angle alpha around a pseudo-random ortogonal axis
#ifdef PSRDNOISE_FAST_ROTATION
    // Fast algorithm, but without dynamic shortcut for alpha = 0
    float4 qx = St;         // q' = norm ( cross(s, n) )  on the equator
    float4 qy = -Ct; 
    float4 qz = float4(0.0, 0.0, 0.0, 0.0);

    float4 px =  sz * qy;   // p' = cross(q, s)
    float4 py = -sz * qx;
    float4 pz = sz_prime;

    psi += alpha;         // psi and alpha in the same plane
    float4 Sa = sin(psi);
    float4 Ca = cos(psi);

    gx = Ca * px + Sa * qx;
    gy = Ca * py + Sa * qy;
    gz = Ca * pz + Sa * qz;
#else
    // Slightly slower algorithm, but with g = s for alpha = 0, and a
    // useful conditional speedup for alpha = 0 across all fragments
    if(alpha != 0.0) {
        float4 Sp = sin(psi);          // q' from psi on equator
        float4 Cp = cos(psi);

        float4 px = Ct * sz_prime;     // px = sx
        float4 py = St * sz_prime;     // py = sy
        float4 pz = sz;

        float4 Ctp = St*Sp - Ct*Cp;    // q = (rotate( cross(s,n), dot(s,n))(q')
        float4 qx = lerp( Ctp*St, Sp, sz);
        float4 qy = lerp(-Ctp*Ct, Cp, sz);
        float4 qz = -(py*Cp + px*Sp);

        float sin_a = sin(alpha);
        float cos_a = cos(alpha);
        float4 Sa = float4(sin_a, sin_a, sin_a, sin_a);       // psi and alpha in different planes
        float4 Ca = float4(cos_a, cos_a, cos_a, cos_a);

        gx = Ca * px + Sa * qx;
        gy = Ca * py + Sa * qy;
        gz = Ca * pz + Sa * qz;
    }
    else {
        gx = Ct * sz_prime;  // alpha = 0, use s directly as gradient
        gy = St * sz_prime;
        gz = sz;  
    }
#endif

    // Reorganize for dot products below
    float3 g0 = float3(gx.x, gy.x, gz.x);
    float3 g1 = float3(gx.y, gy.y, gz.y);
    float3 g2 = float3(gx.z, gy.z, gz.z);
    float3 g3 = float3(gx.w, gy.w, gz.w);

    // Radial decay with distance from each simplex corner
    float4 w = 0.5 - float4(dot(x0,x0), dot(x1,x1), dot(x2,x2), dot(x3,x3));
    w = max(w, 0.0);
    float4 w2 = w * w;
    float4 w3 = w2 * w;

    // The value of the linear ramp from each of the corners
    float4 gdotx = float4(dot(g0,x0), dot(g1,x1), dot(g2,x2), dot(g3,x3));

    // Multiply by the radial decay and sum up the noise value
    float n = dot(w3, gdotx);

    // Compute the first order partial derivatives
    float4 dw = -6.0 * w2 * gdotx;
    float3 dn0 = w3.x * g0 + dw.x * x0;
    float3 dn1 = w3.y * g1 + dw.y * x1;
    float3 dn2 = w3.z * g2 + dw.z * x2;
    float3 dn3 = w3.w * g3 + dw.w * x3;
    gradient = 39.5 * (dn0 + dn1 + dn2 + dn3);

    // Scale the return value to fit nicely into the range [-1,1]
    return 39.5 * n;

}

float psrdnoise(float3 x, float3 period, float alpha, out float3 gradient, out float3 dg, out float3 dg2) {

#ifndef PSRDNOISE_PERLIN_GRID
    // Transformation matrices for the axis-aligned simplex grid
    const float3x3 M = float3x3(0.0, 1.0, 1.0,
                        1.0, 0.0, 1.0,
                        1.0, 1.0, 0.0);

    const float3x3 Mi = float3x3(-0.5, 0.5, 0.5,
                            0.5,-0.5, 0.5,
                            0.5, 0.5,-0.5);
#endif

    float3 uvw = float3(0.0, 0.0, 0.0);

    // Transform to simplex space (tetrahedral grid)
#ifndef PSRDNOISE_PERLIN_GRID
    // Use matrix multiplication, let the compiler optimise
    uvw = mul(M, x);
#else
    // Optimised transformation to uvw (slightly faster than
    // the equivalent matrix multiplication on most platforms)
    uvw = x + dot(x, float3(0.3333333, 0.3333333, 0.3333333));
#endif

    // Determine which simplex we're in, i0 is the "base corner"
    float3 i0 = floor(uvw);
    float3 f0 = frac(uvw); // coords within "skewed cube"

    // To determine which simplex corners are closest, rank order the
    // magnitudes of u,v,w, resolving ties in priority order u,v,w,
    // and traverse the four corners from largest to smallest magnitude.
    // o1, o2 are offsets in simplex space to the 2nd and 3rd corners.
    float3 g_ = step(f0.xyx, f0.yzz); // Makes comparison "less-than"
    float3 l_ = 1.0 - g_;             // complement is "greater-or-equal"
    float3 g = float3(l_.z, g_.xy);
    float3 l = float3(l_.xy, g_.z);
    float3 o1 = min( g, l );
    float3 o2 = max( g, l );

    // Enumerate the remaining simplex corners
    float3 i1 = i0 + o1;
    float3 i2 = i0 + o2;
    float3 i3 = i0 + float3(1.0, 1.0, 1.0);

    float3 v0, v1, v2, v3;

    // Transform the corners back to texture space
#ifndef PSRDNOISE_PERLIN_GRID
    v0 = mul(Mi, i0);
    v1 = mul(Mi, i1);
    v2 = mul(Mi, i2);
    v3 = mul(Mi, i3);
#else
    // Optimised transformation (mostly slightly faster than a matrix)
    v0 = i0 - dot(i0, float3(0.166666667, 0.166666667, 0.166666667));
    v1 = i1 - dot(i1, float3(0.166666667, 0.166666667, 0.166666667));
    v2 = i2 - dot(i2, float3(0.166666667, 0.166666667, 0.166666667));
    v3 = i3 - dot(i3, float3(0.166666667, 0.166666667, 0.166666667));
#endif

    // Compute vectors to each of the simplex corners
    float3 x0 = x - v0;
    float3 x1 = x - v1;
    float3 x2 = x - v2;
    float3 x3 = x - v3;

    if(any(greaterThan(period, float3(0.0, 0.0, 0.0)))) {
        // Wrap to periods and transform back to simplex space
        float4 vx = float4(v0.x, v1.x, v2.x, v3.x);
        float4 vy = float4(v0.y, v1.y, v2.y, v3.y);
        float4 vz = float4(v0.z, v1.z, v2.z, v3.z);
        // Wrap to periods where specified
        if(period.x > 0.0) vx = mod(vx, period.x);
        if(period.y > 0.0) vy = mod(vy, period.y);
        if(period.z > 0.0) vz = mod(vz, period.z);
        // Transform back
#ifndef PSRDNOISE_PERLIN_GRID
        i0 = mul(M, float3(vx.x, vy.x, vz.x));
        i1 = mul(M, float3(vx.y, vy.y, vz.y));
        i2 = mul(M, float3(vx.z, vy.z, vz.z));
        i3 = mul(M, float3(vx.w, vy.w, vz.w));
#else
        v0 = float3(vx.x, vy.x, vz.x);
        v1 = float3(vx.y, vy.y, vz.y);
        v2 = float3(vx.z, vy.z, vz.z);
        v3 = float3(vx.w, vy.w, vz.w);
        // Transform wrapped coordinates back to uvw
        i0 = v0 + dot(v0, float3(0.3333333, 0.3333333, 0.3333333));
        i1 = v1 + dot(v1, float3(0.3333333, 0.3333333, 0.3333333));
        i2 = v2 + dot(v2, float3(0.3333333, 0.3333333, 0.3333333));
        i3 = v3 + dot(v3, float3(0.3333333, 0.3333333, 0.3333333));
#endif
        // Fix rounding errors
        i0 = floor(i0 + 0.5);
        i1 = floor(i1 + 0.5);
        i2 = floor(i2 + 0.5);
        i3 = floor(i3 + 0.5);
    }

    // Compute one pseudo-random hash value for each corner
    float4 hash = permute( permute( permute( 
                float4(i0.z, i1.z, i2.z, i3.z ))
                + float4(i0.y, i1.y, i2.y, i3.y ))
                + float4(i0.x, i1.x, i2.x, i3.x ));

    // Compute generating gradients from a Fibonacci spiral on the unit sphere
    float4 theta = hash * 3.883222077;  // 2*pi/golden ratio
    float4 sz    = hash * -0.006920415 + 0.996539792; // 1-(hash+0.5)*2/289
    float4 psi   = hash * 0.108705628 ; // 10*pi/289, chosen to avoid correlation

    float4 Ct = cos(theta);
    float4 St = sin(theta);
    float4 sz_prime = sqrt( 1.0 - sz*sz ); // s is a point on a unit fib-sphere

    float4 gx, gy, gz;

    // Rotate gradients by angle alpha around a pseudo-random ortogonal axis
#ifdef PSRDNOISE_FAST_ROTATION
    // Fast algorithm, but without dynamic shortcut for alpha = 0
    float4 qx = St;         // q' = norm ( cross(s, n) )  on the equator
    float4 qy = -Ct; 
    float4 qz = float4(0.0, 0.0, 0.0, 0.0);

    float4 px =  sz * qy;   // p' = cross(q, s)
    float4 py = -sz * qx;
    float4 pz = sz_prime;

    psi += alpha;         // psi and alpha in the same plane
    float4 Sa = sin(psi);
    float4 Ca = cos(psi);

    gx = Ca * px + Sa * qx;
    gy = Ca * py + Sa * qy;
    gz = Ca * pz + Sa * qz;
    #else
    // Slightly slower algorithm, but with g = s for alpha = 0, and a
    // strong conditional speedup for alpha = 0 across all fragments
    if(alpha != 0.0) {
        float4 Sp = sin(psi);          // q' from psi on equator
        float4 Cp = cos(psi);

        float4 px = Ct * sz_prime;     // px = sx
        float4 py = St * sz_prime;     // py = sy
        float4 pz = sz;

        float4 Ctp = St*Sp - Ct*Cp;    // q = (rotate( cross(s,n), dot(s,n))(q')
        float4 qx = lerp( Ctp*St, Sp, sz);
        float4 qy = lerp(-Ctp*Ct, Cp, sz);
        float4 qz = -(py*Cp + px*Sp);

        float sin_a = sin(alpha);
        float cos_a = cos(alpha);
        float4 Sa = float4(sin_a, sin_a, sin_a, sin_a);       // psi and alpha in different planes
        float4 Ca = float4(cos_a, cos_a, cos_a, cos_a);

        gx = Ca * px + Sa * qx;
        gy = Ca * py + Sa * qy;
        gz = Ca * pz + Sa * qz;
    }
    else {
        gx = Ct * sz_prime;  // alpha = 0, use s directly as gradient
        gy = St * sz_prime;
        gz = sz;  
    }
#endif

    // Reorganize for dot products below
    float3 g0 = float3(gx.x, gy.x, gz.x);
    float3 g1 = float3(gx.y, gy.y, gz.y);
    float3 g2 = float3(gx.z, gy.z, gz.z);
    float3 g3 = float3(gx.w, gy.w, gz.w);

    // Radial decay with distance from each simplex corner
    float4 w = 0.5 - float4(dot(x0,x0), dot(x1,x1), dot(x2,x2), dot(x3,x3));
    w = max(w, 0.0);
    float4 w2 = w * w;
    float4 w3 = w2 * w;

    // The value of the linear ramp from each of the corners
    float4 gdotx = float4(dot(g0,x0), dot(g1,x1), dot(g2,x2), dot(g3,x3));

    // Multiply by the radial decay and sum up the noise value
    float n = dot(w3, gdotx);

    // Compute the first order partial derivatives
    float4 dw = -6.0 * w2 * gdotx;
    float3 dn0 = w3.x * g0 + dw.x * x0;
    float3 dn1 = w3.y * g1 + dw.y * x1;
    float3 dn2 = w3.z * g2 + dw.z * x2;
    float3 dn3 = w3.w * g3 + dw.w * x3;
    gradient = 39.5 * (dn0 + dn1 + dn2 + dn3);

    // Compute the second order partial derivatives
    float4 dw2 = 24.0 * w * gdotx;
    float3 dga0 = dw2.x * x0 * x0 - 6.0 * w2.x * (gdotx.x + 2.0 * g0 * x0);
    float3 dga1 = dw2.y * x1 * x1 - 6.0 * w2.y * (gdotx.y + 2.0 * g1 * x1);
    float3 dga2 = dw2.z * x2 * x2 - 6.0 * w2.z * (gdotx.z + 2.0 * g2 * x2);
    float3 dga3 = dw2.w * x3 * x3 - 6.0 * w2.w * (gdotx.w + 2.0 * g3 * x3);
    dg = 35.0 * (dga0 + dga1 + dga2 + dga3); // (d2n/dx2, d2n/dy2, d2n/dz2)
    float3 dgb0 = dw2.x * x0 * x0.yzx - 6.0 * w2.x * (g0 * x0.yzx + g0.yzx * x0);
    float3 dgb1 = dw2.y * x1 * x1.yzx - 6.0 * w2.y * (g1 * x1.yzx + g1.yzx * x1);
    float3 dgb2 = dw2.z * x2 * x2.yzx - 6.0 * w2.z * (g2 * x2.yzx + g2.yzx * x2);
    float3 dgb3 = dw2.w * x3 * x3.yzx - 6.0 * w2.w * (g3 * x3.yzx + g3.yzx * x3);
    dg2 = 39.5 * (dgb0 + dgb1 + dgb2 + dgb3); // (d2n/dxy, d2n/dyz, d2n/dxz)

    // Scale the return value to fit nicely into the range [-1,1]
    return 39.5 * n;
}

float psrdnoise(float3 x, float3 period, float alpha) {
    float3 g = float3(0.0, 0.0, 0.0);
    return psrdnoise(x, period, alpha, g);
}

float psrdnoise(float3 x, float3 period) {
    return psrdnoise(x, period, 0.0);
}

float psrdnoise(float3 x) {
    return psrdnoise(x, float3(0.0, 0.0, 0.0));
}
#endif

License

Copyright 2021-2023 by Stefan Gustavson and Ian McEwan. Published under the terms of the MIT license: https://opensource.org/license/mit/

LYGIA is dual-licensed under the Prosperity License and the Patron License for sponsors and contributors.

Sponsors and contributors are automatically added to the Patron License and they can ignore the any non-commercial rule of the Prosperity Licensed software (please take a look to the exception).

It's also possible to get a permanent comercial license hook to a single and specific version of LYGIA.

Get the latest news and releases

Sign up for the news letter bellow, joing the LYGIA's channel on Discord or follow the Github repository