647 lines
No EOL
21 KiB
GLSL
647 lines
No EOL
21 KiB
GLSL
#include "../math/mod289.glsl"
|
|
#include "../math/permute.glsl"
|
|
#include "../math/taylorInvSqrt.glsl"
|
|
#include "../math/grad4.glsl"
|
|
|
|
/*
|
|
contributors: [Stefan Gustavson, Ian McEwan]
|
|
description: |
|
|
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
|
|
|
|
use: <float> psrdnoise(<vec2|vec3> x, <vec2|vec3> period, <float> alpha, out <vec2|vec3> gradient [, out <vec3> dg])
|
|
options:
|
|
- PSRDNOISE_PERLIN_GRID
|
|
- PSRDNOISE_FAST_ROTATION
|
|
license: |
|
|
Copyright 2021-2023 by Stefan Gustavson and Ian McEwan.
|
|
Published under the terms of the MIT license:
|
|
https://opensource.org/license/mit/
|
|
*/
|
|
|
|
#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 |