1#define_import_path bevy_gaussian_splatting::noise
2
3
4//  MIT License. © Ian McEwan, Stefan Gustavson, Munrocket, Johan Helsing
5//
6fn mod289_2d(x: vec2<f32>) -> vec2<f32> {
7    return x - floor(x * (1.0 / 289.0)) * 289.0;
8}
9
10fn mod289_3d(x: vec3<f32>) -> vec3<f32> {
11    return x - floor(x * (1.0 / 289.0)) * 289.0;
12}
13
14fn permute_3d(x: vec3<f32>) -> vec3<f32> {
15    return mod289_3d(((x * 34.0) + 1.0) * x);
16}
17
18//  MIT License. © Ian McEwan, Stefan Gustavson, Munrocket
19fn simplex_2d(v: vec2<f32>) -> f32 {
20    let C = vec4(
21        0.211324865405187, // (3.0-sqrt(3.0))/6.0
22        0.366025403784439, // 0.5*(sqrt(3.0)-1.0)
23        -0.577350269189626, // -1.0 + 2.0 * C.x
24        0.024390243902439 // 1.0 / 41.0
25    );
26
27    // First corner
28    var i = floor(v + dot(v, C.yy));
29    let x0 = v - i + dot(i, C.xx);
30
31    // Other corners
32    var i1 = select(vec2(0., 1.), vec2(1., 0.), x0.x > x0.y);
33
34    // x0 = x0 - 0.0 + 0.0 * C.xx ;
35    // x1 = x0 - i1 + 1.0 * C.xx ;
36    // x2 = x0 - 1.0 + 2.0 * C.xx ;
37    var x12 = x0.xyxy + C.xxzz;
38    x12.x = x12.x - i1.x;
39    x12.y = x12.y - i1.y;
40
41    // Permutations
42    i = mod289_2d(i); // Avoid truncation effects in permutation
43
44    var p = permute_3d(permute_3d(i.y + vec3(0., i1.y, 1.)) + i.x + vec3(0., i1.x, 1.));
45    var m = max(0.5 - vec3(dot(x0, x0), dot(x12.xy, x12.xy), dot(x12.zw, x12.zw)), vec3(0.));
46    m *= m;
47    m *= m;
48
49    // Gradients: 41 points uniformly over a line, mapped onto a diamond.
50    // The ring size 17*17 = 289 is close to a multiple of 41 (41*7 = 287)
51    let x = 2. * fract(p * C.www) - 1.;
52    let h = abs(x) - 0.5;
53    let ox = floor(x + 0.5);
54    let a0 = x - ox;
55
56    // Normalize gradients implicitly by scaling m
57    // Approximation of: m *= inversesqrt( a0*a0 + h*h );
58    m *= 1.79284291400159 - 0.85373472095314 * (a0 * a0 + h * h);
59
60    // Compute final noise value at P
61    let g = vec3(a0.x * x0.x + h.x * x0.y, a0.yz * x12.xz + h.yz * x12.yw);
62    return 130. * dot(m, g);
63}
64
65
66// MIT License. © Stefan Gustavson, Munrocket
67fn permute4_(x: vec4<f32>) -> vec4<f32> { return ((x * 34.0 + 1.0) * x) % vec4<f32>(289.0); }
68fn taylorInvSqrt4_(r: vec4<f32>) -> vec4<f32> { return 1.79284291400159 - 0.85373472095314 * r; }
69fn fade3_(t: vec3<f32>) -> vec3<f32> { return t * t * t * (t * (t * 6.0 - 15.0) + 10.0); }
70
71fn perlin_3d(P: vec3<f32>) -> f32 {
72    var Pi0 : vec3<f32> = floor(P); // Integer part for indexing
73    var Pi1 : vec3<f32> = Pi0 + vec3<f32>(1.0); // Integer part + 1
74    Pi0 = Pi0 % vec3<f32>(289.0);
75    Pi1 = Pi1 % vec3<f32>(289.0);
76    let Pf0 = fract(P); // Fractional part for interpolation
77    let Pf1 = Pf0 - vec3<f32>(1.0); // Fractional part - 1.
78    let ix = vec4<f32>(Pi0.x, Pi1.x, Pi0.x, Pi1.x);
79    let iy = vec4<f32>(Pi0.yy, Pi1.yy);
80    let iz0 = Pi0.zzzz;
81    let iz1 = Pi1.zzzz;
82
83    let ixy = permute4_(permute4_(ix) + iy);
84    let ixy0 = permute4_(ixy + iz0);
85    let ixy1 = permute4_(ixy + iz1);
86
87    var gx0: vec4<f32> = ixy0 / 7.0;
88    var gy0: vec4<f32> = fract(floor(gx0) / 7.0) - 0.5;
89    gx0 = fract(gx0);
90    var gz0: vec4<f32> = vec4<f32>(0.5) - abs(gx0) - abs(gy0);
91    var sz0: vec4<f32> = step(gz0, vec4<f32>(0.0));
92    gx0 = gx0 + sz0 * (step(vec4<f32>(0.0), gx0) - 0.5);
93    gy0 = gy0 + sz0 * (step(vec4<f32>(0.0), gy0) - 0.5);
94
95    var gx1: vec4<f32> = ixy1 / 7.0;
96    var gy1: vec4<f32> = fract(floor(gx1) / 7.0) - 0.5;
97    gx1 = fract(gx1);
98    var gz1: vec4<f32> = vec4<f32>(0.5) - abs(gx1) - abs(gy1);
99    var sz1: vec4<f32> = step(gz1, vec4<f32>(0.0));
100    gx1 = gx1 - sz1 * (step(vec4<f32>(0.0), gx1) - 0.5);
101    gy1 = gy1 - sz1 * (step(vec4<f32>(0.0), gy1) - 0.5);
102
103    var g000: vec3<f32> = vec3<f32>(gx0.x, gy0.x, gz0.x);
104    var g100: vec3<f32> = vec3<f32>(gx0.y, gy0.y, gz0.y);
105    var g010: vec3<f32> = vec3<f32>(gx0.z, gy0.z, gz0.z);
106    var g110: vec3<f32> = vec3<f32>(gx0.w, gy0.w, gz0.w);
107    var g001: vec3<f32> = vec3<f32>(gx1.x, gy1.x, gz1.x);
108    var g101: vec3<f32> = vec3<f32>(gx1.y, gy1.y, gz1.y);
109    var g011: vec3<f32> = vec3<f32>(gx1.z, gy1.z, gz1.z);
110    var g111: vec3<f32> = vec3<f32>(gx1.w, gy1.w, gz1.w);
111
112    let norm0 = taylorInvSqrt4_(
113        vec4<f32>(dot(g000, g000), dot(g010, g010), dot(g100, g100), dot(g110, g110)));
114    g000 = g000 * norm0.x;
115    g010 = g010 * norm0.y;
116    g100 = g100 * norm0.z;
117    g110 = g110 * norm0.w;
118    let norm1 = taylorInvSqrt4_(
119        vec4<f32>(dot(g001, g001), dot(g011, g011), dot(g101, g101), dot(g111, g111)));
120    g001 = g001 * norm1.x;
121    g011 = g011 * norm1.y;
122    g101 = g101 * norm1.z;
123    g111 = g111 * norm1.w;
124
125    let n000 = dot(g000, Pf0);
126    let n100 = dot(g100, vec3<f32>(Pf1.x, Pf0.yz));
127    let n010 = dot(g010, vec3<f32>(Pf0.x, Pf1.y, Pf0.z));
128    let n110 = dot(g110, vec3<f32>(Pf1.xy, Pf0.z));
129    let n001 = dot(g001, vec3<f32>(Pf0.xy, Pf1.z));
130    let n101 = dot(g101, vec3<f32>(Pf1.x, Pf0.y, Pf1.z));
131    let n011 = dot(g011, vec3<f32>(Pf0.x, Pf1.yz));
132    let n111 = dot(g111, Pf1);
133
134    var fade_xyz: vec3<f32> = fade3_(Pf0);
135    let temp = vec4<f32>(f32(fade_xyz.z)); // simplify after chrome bug fix
136    let n_z = mix(vec4<f32>(n000, n100, n010, n110), vec4<f32>(n001, n101, n011, n111), temp);
137    let n_yz = mix(n_z.xy, n_z.zw, vec2<f32>(f32(fade_xyz.y))); // simplify after chrome bug fix
138    let n_xyz = mix(n_yz.x, n_yz.y, fade_xyz.x);
139    return 2.2 * n_xyz;
140}
141
142
143//  <https://www.shadertoy.com/view/Xd23Dh>
144//  by Inigo Quilez
145fn hash_23_(p: vec2<f32>) -> vec3<f32> {
146    let q = vec3<f32>(dot(p, vec2<f32>(127.1, 311.7)),
147        dot(p, vec2<f32>(269.5, 183.3)),
148        dot(p, vec2<f32>(419.2, 371.9)));
149    return fract(sin(q) * 43758.5453);
150}
151
152fn voro_2d(x: vec2<f32>, u: f32, v: f32) -> f32 {
153    let p = floor(x);
154    let f = fract(x);
155    let k = 1.0 + 63.0 * pow(1. - v, 4.0);
156    var va: f32 = 0.0;
157    var wt: f32 = 0.0;
158    for(var j: i32 = -2; j <= 2; j = j + 1) {
159      for(var i: i32 = -2; i <= 2; i = i + 1) {
160        let g = vec2<f32>(f32(i), f32(j));
161        let o = hash_23_(p + g) * vec3<f32>(u, u, 1.0);
162        let r = g - f + o.xy;
163        let d = dot(r, r);
164        let ww = pow(1. - smoothstep(0.0, 1.414, sqrt(d)), k);
165        va = va + o.z * ww;
166        wt = wt + ww;
167      }
168    }
169    return va / wt;
170}
171
172
173
174fn nrand(n: vec2<f32>) -> f32 {
175    return fract(sin(dot(n, vec2<f32>(12.9898, 4.1414))) * 43758.5453);
176}
177
178fn noise_2d(n: vec2<f32>) -> f32 {
179    let d = vec2<f32>(0.0, 1.0);
180    let b = floor(n);
181    let f = smoothstep(vec2<f32>(0.0), vec2<f32>(1.0), fract(n));
182    return mix(mix(nrand(b), nrand(b + d.yx), f.x), mix(nrand(b + d.xy), nrand(b + d.yy), f.x), f.y);
183}
184
185
186// https://www.shadertoy.com/view/MlVSzw
187const ALPHA: f32 = 0.14;
188const INV_ALPHA: f32 = 7.14285714286;
189const K: f32 = 0.08912676813;
190
191fn inv_error_function(x: f32) -> f32 {
192    let y: f32 = log(1.0 - x*x);
193    let z: f32 = K + 0.5 * y;
194    return sqrt(sqrt(z*z - y * INV_ALPHA) - z) * sign(x);
195}
196
197// expects n to be in ~[0, 1]
198fn gaussian_rand(n: vec2<f32>) -> f32 {
199    let x: f32 = nrand(n * 13.7);
200
201    return inv_error_function(x * 2.0 - 1.0) * 0.3;
202}