1//  MIT License. © Ian McEwan, Stefan Gustavson, Munrocket, Johan Helsing
2
3#define_import_path noisy_bevy
4
5fn permute_3_(x: vec3<f32>) -> vec3<f32> {
6    return (((x * 34.) + 1.) * x) % vec3(289.);
7}
8
9fn step_3(edge: vec3<f32>, x: vec3<f32>) -> vec3<f32> {
10    let b = vec3(edge.x < x.x, edge.y <= x.y, edge.z <= x.z);
11    return select(vec3(0.), vec3(1.), b);
12}
13
14fn simplex_noise_2d(v: vec2<f32>) -> f32 {
15    let C = vec4(
16        0.211324865405187, // (3.0 - sqrt(3.0)) / 6.0
17        0.366025403784439, // 0.5 * (sqrt(3.0) - 1.0)
18        -0.577350269189626, // -1.0 + 2.0 * C.x
19        0.024390243902439 // 1.0 / 41.0
20    );
21
22    // first corner
23    var i = floor(v + dot(v, C.yy));
24    let x0 = v - i + dot(i, C.xx);
25
26    // other corners
27    var i1 = select(vec2(0., 1.), vec2(1., 0.), x0.x > x0.y);
28    var x12 = x0.xyxy + C.xxzz - vec4(i1, 0., 0.);
29
30    // permutations
31    i = i % vec2(289.);
32
33    let p = permute_3_(permute_3_(i.y + vec3(0., i1.y, 1.)) + i.x + vec3(0., i1.x, 1.));
34    var m = max(0.5 - vec3(dot(x0, x0), dot(x12.xy, x12.xy), dot(x12.zw, x12.zw)), vec3(0.));
35    m *= m;
36    m *= m;
37
38    // gradients: 41 points uniformly over a line, mapped onto a diamond
39    // the ring size, 17*17 = 289, is close to a multiple of 41 (41*7 = 287)
40    let x = 2. * fract(p * C.www) - 1.;
41    let h = abs(x) - 0.5;
42    let ox = floor(x + 0.5);
43    let a0 = x - ox;
44
45    // normalize gradients implicitly by scaling m
46    // approximation of: m *= inversesqrt(a0 * a0 + h * h);
47    m = m * (1.79284291400159 - 0.85373472095314 * (a0 * a0 + h * h));
48
49    // compute final noise value at P
50    let g = vec3(a0.x * x0.x + h.x * x0.y, a0.yz * x12.xz + h.yz * x12.yw);
51    return 130. * dot(m, g);
52}
53
54fn simplex_noise_2d_seeded(v: vec2<f32>, seed: f32) -> f32 {
55    let C = vec4(
56        0.211324865405187, // (3.0 - sqrt(3.0)) / 6.0
57        0.366025403784439, // 0.5 * (sqrt(3.0) - 1.0)
58        -0.577350269189626, // -1.0 + 2.0 * C.x
59        0.024390243902439 // 1.0 / 41.0
60    );
61
62    // first corner
63    var i = floor(v + dot(v, C.yy));
64    let x0 = v - i + dot(i, C.xx);
65
66    // other corners
67    var i1 = select(vec2(0., 1.), vec2(1., 0.), x0.x > x0.y);
68    var x12 = x0.xyxy + C.xxzz - vec4(i1, 0., 0.);
69
70    // permutations
71    i = i % vec2(289.);
72
73    var p = permute_3_(permute_3_(i.y + vec3(0., i1.y, 1.)) + i.x + vec3(0., i1.x, 1.));
74    p = permute_3_(p + vec3(seed));
75    var m = max(0.5 - vec3(dot(x0, x0), dot(x12.xy, x12.xy), dot(x12.zw, x12.zw)), vec3(0.));
76    m *= m;
77    m *= m;
78
79    // gradients: 41 points uniformly over a line, mapped onto a diamond
80    // the ring size, 17*17 = 289, is close to a multiple of 41 (41*7 = 287)
81    let x = 2. * fract(p * C.www) - 1.;
82    let h = abs(x) - 0.5;
83    let ox = floor(x + 0.5);
84    let a0 = x - ox;
85
86    // normalize gradients implicitly by scaling m
87    // approximation of: m *= inversesqrt(a0 * a0 + h * h);
88    m = m * (1.79284291400159 - 0.85373472095314 * (a0 * a0 + h * h));
89
90    // compute final noise value at P
91    let g = vec3(a0.x * x0.x + h.x * x0.y, a0.yz * x12.xz + h.yz * x12.yw);
92    return 130. * dot(m, g);
93}
94
95fn permute_4_(x: vec4<f32>) -> vec4<f32> {
96    return ((x * 34. + 1.) * x) % vec4<f32>(289.);
97}
98
99fn taylor_inv_sqrt_4_(r: vec4<f32>) -> vec4<f32> {
100    return 1.79284291400159 - 0.85373472095314 * r;
101}
102
103fn simplex_noise_3d(v: vec3<f32>) -> f32 {
104    let C = vec2(1. / 6., 1. / 3.);
105    let D = vec4(0., 0.5, 1., 2.);
106
107    // first corner
108    var i = floor(v + dot(v, C.yyy));
109    let x0 = v - i + dot(i, C.xxx);
110
111    // other corners
112    let g = step_3(x0.yzx, x0.xyz);
113    let l = 1. - g;
114    let i1 = min(g.xyz, l.zxy);
115    let i2 = max(g.xyz, l.zxy);
116
117    // x0 = x0 - 0. + 0. * C
118    let x1 = x0 - i1 + 1. * C.xxx;
119    let x2 = x0 - i2 + 2. * C.xxx;
120    let x3 = x0 - 1. + 3. * C.xxx;
121
122    // permutations
123    i = i % vec3(289.);
124    let p = permute_4_(permute_4_(permute_4_(
125        i.z + vec4(0., i1.z, i2.z, 1.)) +
126        i.y + vec4(0., i1.y, i2.y, 1.)) +
127        i.x + vec4(0., i1.x, i2.x, 1.)
128    );
129
130    // gradients (NxN points uniformly over a square, mapped onto an octahedron)
131    let n_ = 1. / 7.; // N=7
132    let ns = n_ * D.wyz - D.xzx;
133
134    let j = p - 49. * floor(p * ns.z * ns.z); // mod(p, N*N)
135
136    let x_ = floor(j * ns.z);
137    let y_ = floor(j - 7. * x_); // mod(j, N)
138
139    let x = x_ * ns.x + ns.yyyy;
140    let y = y_ * ns.x + ns.yyyy;
141    let h = 1. - abs(x) - abs(y);
142
143    let b0 = vec4(x.xy, y.xy);
144    let b1 = vec4(x.zw, y.zw);
145
146    let s0 = floor(b0) * 2. + 1.;
147    let s1 = floor(b1) * 2. + 1.;
148    let sh = -step(h, vec4(0.));
149
150    let a0 = b0.xzyw + s0.xzyw * sh.xxyy;
151    let a1 = b1.xzyw + s1.xzyw * sh.zzww;
152
153    var p0 = vec3(a0.xy, h.x);
154    var p1 = vec3(a0.zw, h.y);
155    var p2 = vec3(a1.xy, h.z);
156    var p3 = vec3(a1.zw, h.w);
157
158    // normalize gradients
159    let norm = taylor_inv_sqrt_4_(vec4(dot(p0, p0), dot(p1, p1), dot(p2, p2), dot(p3, p3)));
160    p0 = p0 * norm.x;
161    p1 = p1 * norm.y;
162    p2 = p2 * norm.z;
163    p3 = p3 * norm.w;
164
165    // mix final noise value
166    var m = 0.5 - vec4(dot(x0, x0), dot(x1, x1), dot(x2, x2), dot(x3, x3));
167    m = max(m, vec4(0.));
168    m *= m;
169    return 105. * dot(m * m, vec4(dot(p0, x0), dot(p1, x1), dot(p2, x2), dot(p3, x3)));
170}
171
172fn simplex_noise_3d_seeded(v: vec3<f32>, seed: vec3<f32>) -> f32 {
173    let C = vec2(1. / 6., 1. / 3.);
174    let D = vec4(0., 0.5, 1., 2.);
175
176    // first corner
177    var i = floor(v + dot(v, C.yyy));
178    let x0 = v - i + dot(i, C.xxx);
179
180    // other corners
181    let g = step_3(x0.yzx, x0.xyz);
182    let l = 1. - g;
183    let i1 = min(g.xyz, l.zxy);
184    let i2 = max(g.xyz, l.zxy);
185
186    // x0 = x0 - 0. + 0. * C
187    let x1 = x0 - i1 + 1. * C.xxx;
188    let x2 = x0 - i2 + 2. * C.xxx;
189    let x3 = x0 - 1. + 3. * C.xxx;
190
191    // permutations
192    i = i % vec3(289.);
193    let s = floor(seed + vec3(0.5));
194    let p = permute_4_(permute_4_(permute_4_(
195        i.z + vec4(0., i1.z, i2.z, 1.) + s.z) +
196        i.y + vec4(0., i1.y, i2.y, 1.) + s.y) +
197        i.x + vec4(0., i1.x, i2.x, 1.) + s.x
198    );
199
200    // gradients (NxN points uniformly over a square, mapped onto an octahedron)
201    let n_ = 1. / 7.; // N=7
202    let ns = n_ * D.wyz - D.xzx;
203
204    let j = p - 49. * floor(p * ns.z * ns.z); // mod(p, N*N)
205
206    let x_ = floor(j * ns.z);
207    let y_ = floor(j - 7. * x_); // mod(j, N)
208
209    let x = x_ * ns.x + ns.yyyy;
210    let y = y_ * ns.x + ns.yyyy;
211    let h = 1. - abs(x) - abs(y);
212
213    let b0 = vec4(x.xy, y.xy);
214    let b1 = vec4(x.zw, y.zw);
215
216    let s0 = floor(b0) * 2. + 1.;
217    let s1 = floor(b1) * 2. + 1.;
218    let sh = -step(h, vec4(0.));
219
220    let a0 = b0.xzyw + s0.xzyw * sh.xxyy;
221    let a1 = b1.xzyw + s1.xzyw * sh.zzww;
222
223    var p0 = vec3(a0.xy, h.x);
224    var p1 = vec3(a0.zw, h.y);
225    var p2 = vec3(a1.xy, h.z);
226    var p3 = vec3(a1.zw, h.w);
227
228    // normalize gradients
229    let norm = taylor_inv_sqrt_4_(vec4(dot(p0, p0), dot(p1, p1), dot(p2, p2), dot(p3, p3)));
230    p0 = p0 * norm.x;
231    p1 = p1 * norm.y;
232    p2 = p2 * norm.z;
233    p3 = p3 * norm.w;
234
235    // mix final noise value
236    var m = 0.6 - vec4(dot(x0, x0), dot(x1, x1), dot(x2, x2), dot(x3, x3));
237    m = max(m, vec4(0.));
238    m *= m;
239    return 42. * dot(m * m, vec4(dot(p0, x0), dot(p1, x1), dot(p2, x2), dot(p3, x3)));
240}
241
242// higher level concepts:
243
244/// Fractional brownian motion (fbm) based on 2d simplex noise
245fn fbm_simplex_2d(pos: vec2<f32>, octaves: i32, lacunarity: f32, gain: f32) -> f32 {
246    var sum = 0.;
247    var amplitude = 1.;
248    var frequency = 1.;
249
250    for (var i = 0; i < octaves; i+= 1) {
251        sum += simplex_noise_2d(pos * frequency) * amplitude;
252        amplitude *= gain;
253        frequency *= lacunarity;
254    }
255
256    return sum;
257}
258
259/// Fractional brownian motion (fbm) based on seeded 2d simplex noise
260fn fbm_simplex_2d_seeded(pos: vec2<f32>, octaves: i32, lacunarity: f32, gain: f32, seed: f32) -> f32 {
261    var sum = 0.;
262    var amplitude = 1.;
263    var frequency = 1.;
264
265    for (var i = 0; i < octaves; i+= 1) {
266        sum += simplex_noise_2d_seeded(pos * frequency, seed) * amplitude;
267        amplitude *= gain;
268        frequency *= lacunarity;
269    }
270
271    return sum;
272}
273
274const max_warp_iterations = 4; // Warping has diminishing returns due to the falloff param, so we don't need many iterations. Faloff makes it look more natural.
275
276struct WarpResult {
277    noise_value: f32,
278    // The history of warped coordinates, where positions[0] is the last iteration, positions[1] is second to last, etc.
279    // Can be useful for mixing colors.
280    positions: array<vec2f, max_warp_iterations>
281}
282
283/// A technique that distorts the position before feeding it to the noise
284/// inspired by https://iquilezles.org/articles/warp/
285fn fbm_simplex_2d_warp_seeded(pos_initial: vec2<f32>, octaves: i32, lacunarity: f32, gain: f32, seed: f32, warp_iterations: i32, warp_scale: vec2<f32>, falloff: f32) -> WarpResult {
286    var scale = 1.0;
287    var positions = array<vec2f, max_warp_iterations>();
288    var pos = pos_initial;
289
290    // Clamp warp_iterations to max_warp_iterations
291    let iterations = min(warp_iterations, i32(max_warp_iterations));
292
293    for (var i: i32 = 0; i < iterations; i++) {
294        pos.x += scale * warp_scale.x * fbm_simplex_2d_seeded(pos, octaves, lacunarity, gain, seed);
295        pos.y += scale * warp_scale.y * fbm_simplex_2d_seeded(pos, octaves, lacunarity, gain, seed);
296
297        // Store positions in reverse order for easier user access (last iteration at index 0)
298        let index = max_warp_iterations - 1 - i;
299        if (index < max_warp_iterations) {
300            positions[index] = pos;
301        }
302
303        scale *= falloff;
304    }
305
306    let noise_value = fbm_simplex_2d_seeded(positions[0], octaves, lacunarity, gain, seed);
307
308    return WarpResult(noise_value, positions);
309}
310
311/// Fractional brownian motion (fbm) based on 3d simplex noise
312fn fbm_simplex_3d(pos: vec3<f32>, octaves: i32, lacunarity: f32, gain: f32) -> f32 {
313    var sum = 0.;
314    var amplitude = 1.;
315    var frequency = 1.;
316
317    for (var i = 0; i < octaves; i+= 1) {
318        sum += simplex_noise_3d(pos * frequency) * amplitude;
319        amplitude *= gain;
320        frequency *= lacunarity;
321    }
322
323    return sum;
324}
325
326/// Fractional brownian motion (fbm) based on seeded 3d simplex noise
327fn fbm_simplex_3d_seeded(pos: vec3<f32>, octaves: i32, lacunarity: f32, gain: f32, seed: vec3<f32>) -> f32 {
328    var sum = 0.;
329    var amplitude = 1.;
330    var frequency = 1.;
331
332    for (var i = 0; i < octaves; i+= 1) {
333        sum += simplex_noise_3d_seeded(pos * frequency, seed) * amplitude;
334        amplitude *= gain;
335        frequency *= lacunarity;
336    }
337
338    return sum;
339}
340
341// MIT license, ported from https://github.com/bevy-interstellar/wgsl_noise
342/// Cellular noise, lower jitter makes the patern more regular
343/// The x component (F1) of the returned result represents the distance to the nearest feature point from the input position
344/// The y component (F2) represents the distance to the second nearest feature point from the input position
345fn worley_2d(pos: vec2<f32>, jitter: f32) -> vec2<f32> {
346    let k = 0.142857142857; // 1/7
347    let ko = 0.428571428571; // 3/7
348
349    // Determine the grid cell and fractional position
350    let pi = floor(pos);
351    let pf = fract(pos);
352
353    // Define offset indices for neighboring grid cells
354    let oi = vec3(-1.0, 0.0, 1.0);
355    let of_ = vec3(-0.5, 0.5, 1.5);
356
357    // Permute the grid cell indices to get unique values for each cell
358    let px = permute_3_(pi.x + oi);
359    var p = permute_3_(px.x + pi.y + oi);  // p11, p12, p13
360
361    var ox = fract(p * k) - ko;
362    var oy = (floor(p * k) % 7.0) * k - ko;
363    var dx = pf.x + 0.5 + jitter * ox;
364    var dy = pf.y - of_ + jitter * oy;
365    var d1 = dx * dx + dy * dy;  // d11, d12, d13, squared
366
367    p = permute_3_(px.y + pi.y + oi); // p21, p22, p23
368    ox = fract(p * k) - ko;
369    oy = (floor(p * k) % 7.0) * k - ko;
370    dx = pf.x - 0.5 + jitter * ox;
371    dy = pf.y - of_ + jitter * oy;
372    var d2 = dx * dx + dy * dy; // d21, d22, d23, squared
373
374    p = permute_3_(px.z + pi.y + oi); // p31, p32, p33
375    ox = fract(p * k) - ko;
376    oy = (floor(p * k) % 7.0) * k - ko;
377    dx = pf.x - 1.5 + jitter * ox;
378    dy = pf.y - of_ + jitter * oy;
379    let d3 = dx * dx + dy * dy; // d31, d32, d33, squared
380
381    // Sort out the two smallest distances (F1, F2)
382    let d1a = min(d1, d2);
383    d2 = max(d1, d2);               // Swap to keep candidates for F2
384    d2 = min(d2, d3);               // neither F1 nor F2 are now in d3
385    d1 = min(d1a, d2);              // F1 is now in d1
386    d2 = max(d1a, d2);              // Swap to keep candidates for F2
387
388    if d1.x > d1.y {                // Swap if smaller
389        let tmp = d1.x;
390        d1.x = d1.y;
391        d1.y = tmp;
392    }
393    if d1.x > d1.z {                // F1 is in d1.x
394        let tmp = d1.x;
395        d1.x = d1.z;
396        d1.z = tmp;
397    }
398
399    d1.y = min(d1.y, d2.y);         // F2 is now not in d2.yz
400    d1.z = min(d1.z, d2.z);
401    d1.y = min(d1.y, d1.z);         // nor in  d1.z
402    d1.y = min(d1.y, d2.x);         // F2 is in d1.y, we're done.
403    return sqrt(d1.xy);
404}