lygia/generative/psrdnoise)Tiling simplex flow noise in 2-D and 3-D from https://github.com/stegu/psrdnoise\">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\">https://github.com/stegu/psrdnoise
Dependencies:
Use:
<float> psrdnoise(<vec2|vec3> x, <vec2|vec3> period, <float> alpha, out <vec2|vec3> gradient [, out <vec3> dg])
#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:
lygia/math/mod289.glsllygia/math/permute.glsllygia/math/taylorInvSqrt.glsllygia/math/grad4.glsllygia/math/mod.glsllygia/math/greaterThan.glslUse:
<float> psrdnoise(<float2|float3> x, <float2|float3> period, <float> alpha, out <float2|float3> gradient [, out <float3> dg])
#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
static const float3x3 M = float3x3(0.0, 1.0, 1.0,
1.0, 0.0, 1.0,
1.0, 1.0, 0.0);
static 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
static const float3x3 M = float3x3(0.0, 1.0, 1.0,
1.0, 0.0, 1.0,
1.0, 1.0, 0.0);
static 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
Dependencies:
lygia/math/greaterThan.glsllygia/math/mod.glsllygia/math/mod289.glsllygia/math/permute.glsllygia/math/taylorInvSqrt.glsllygia/math/grad4.glslUse:
<float> psrdnoise(<float2|float3> x, <float2|float3> period, <float> alpha, out <float2|float3> gradient [, out <float3> dg])
#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 = fract(uv);
// o1 is the offset 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);
float3 iv = float3(0.0);
float3 xw = float3(0.0);
float3 yw = float3(0.0);
// Wrap to periods, if desired
if(any(greaterThan(period, float2(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 = fract(uv);
// o1 is the offset 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)))) {
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);
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));
}
float psrdnoise(float3 x, float3 period, float alpha, out float3 gradient) {
#ifndef PSRDNOISE_PERLIN_GRID
// Transformation matrices for the axis-aligned simplex grid
const matrix<float, 3, 3> M = matrix<float, 3, 3>(0.0, 1.0, 1.0,
1.0, 0.0, 1.0,
1.0, 1.0, 0.0);
const matrix<float, 3, 3> Mi = matrix<float, 3, 3>(-0.5, 0.5, 0.5,
0.5,-0.5, 0.5,
0.5, 0.5,-0.5);
#endif
float3 uvw = float3(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, float3(0.33333333));
#endif
// Determine which simplex we're in, i0 is the "base corner"
float3 i0 = floor(uvw);
float3 f0 = fract(uvw); // coords with"skewed cube"
// To determine which simplex corners are closest, rank order the
// magnitudes of u,v,w, resolving ties priority order u,v,w,
// and traverse the four corners from largest to smallest magnitude.
// o1, o2 are offsets 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);
float3 v0 = float3(0.0);
float3 v1 = float3(0.0);
float3 v2 = float3(0.0);
float3 v3 = float3(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, 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)))) {
// 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 = M * float3(vx.x, vy.x, vz.x);
i1 = M * float3(vx.y, vy.y, vz.y);
i2 = M * float3(vx.z, vy.z, vz.z);
i3 = 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);
float4 gy = float4(0.0);
float4 gz = float4(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);
float4 px = sz * qy; // p' = cross(q, s)
float4 py = -sz * qx;
float4 pz = sz_prime;
psi += alpha; // psi and alpha 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 = mix( Ctp*St, Sp, sz);
float4 qy = mix(-Ctp*Ct, Cp, sz);
float4 qz = -(py*Cp + px*Sp);
float4 Sa = float4(sin(alpha)); // psi and alpha different planes
float4 Ca = float4(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
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 matrix<float, 3, 3> M = matrix<float, 3, 3>(0.0, 1.0, 1.0,
1.0, 0.0, 1.0,
1.0, 1.0, 0.0);
const matrix<float, 3, 3> Mi = matrix<float, 3, 3>(-0.5, 0.5, 0.5,
0.5,-0.5, 0.5,
0.5, 0.5,-0.5);
#endif
float3 uvw = float3(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, float3(0.3333333));
#endif
// Determine which simplex we're in, i0 is the "base corner"
float3 i0 = floor(uvw);
float3 f0 = fract(uvw); // coords with"skewed cube"
// To determine which simplex corners are closest, rank order the
// magnitudes of u,v,w, resolving ties priority order u,v,w,
// and traverse the four corners from largest to smallest magnitude.
// o1, o2 are offsets 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);
float3 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, 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)))) {
// 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 = M * float3(vx.x, vy.x, vz.x);
i1 = M * float3(vx.y, vy.y, vz.y);
i2 = M * float3(vx.z, vy.z, vz.z);
i3 = 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));
i1 = v1 + dot(v1, float3(0.3333333));
i2 = v2 + dot(v2, float3(0.3333333));
i3 = v3 + dot(v3, float3(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);
float4 px = sz * qy; // p' = cross(q, s)
float4 py = -sz * qx;
float4 pz = sz_prime;
psi += alpha; // psi and alpha 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 = mix( Ctp*St, Sp, sz);
float4 qy = mix(-Ctp*Ct, Cp, sz);
float4 qz = -(py*Cp + px*Sp);
float4 Sa = float4(sin(alpha)); // psi and alpha different planes
float4 Ca = float4(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
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);
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));
}
#endif
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.
Sign up for the news letter bellow, joing the LYGIA's channel on Discord or follow the Github repository