1#version 310 es
2
3layout(local_size_x = 64) in;
4
5layout(std430, binding = 0) readonly buffer Distribution
6{
7    vec2 distribution[];
8};
9
10layout(std430, binding = 1) writeonly buffer HeightmapFFT
11{
12    uint heights[];
13};
14
15layout(binding = 2, std140) uniform UBO
16{
17    vec4 uModTime;
18};
19
20vec2 alias(vec2 i, vec2 N)
21{
22    return mix(i, i - N, greaterThan(i, 0.5 * N));
23}
24
25vec4 cmul(vec4 a, vec4 b)
26{
27    vec4 r3 = a.yxwz;
28    vec4 r1 = b.xxzz;
29    vec4 R0 = a * r1;
30    vec4 r2 = b.yyww;
31    vec4 R1 = r2 * r3;
32    return R0 + vec4(-R1.x, R1.y, -R1.z, R1.w);
33}
34
35vec2 cmul(vec2 a, vec2 b)
36{
37    vec2 r3 = a.yx;
38    vec2 r1 = b.xx;
39    vec2 R0 = a * r1;
40    vec2 r2 = b.yy;
41    vec2 R1 = r2 * r3;
42    return R0 + vec2(-R1.x, R1.y);
43}
44
45uint pack2(vec2 v)
46{
47    return packHalf2x16(v);
48}
49
50uvec2 pack4(vec4 v)
51{
52    return uvec2(packHalf2x16(v.xy), packHalf2x16(v.zw));
53}
54
55uvec2 workaround_mix(uvec2 a, uvec2 b, bvec2 sel)
56{
57    return uvec2(sel.x ? b.x : a.x, sel.y ? b.y : a.y);
58}
59
60void generate_heightmap()
61{
62    uvec2 N = gl_WorkGroupSize.xy * gl_NumWorkGroups.xy;
63    uvec2 i = gl_GlobalInvocationID.xy;
64    // Pick out the negative frequency variant.
65    uvec2 wi = workaround_mix(N - i, uvec2(0u), equal(i, uvec2(0u)));
66
67    // Pick out positive and negative travelling waves.
68    vec2 a = distribution[i.y * N.x + i.x];
69    vec2 b = distribution[wi.y * N.x + wi.x];
70
71    vec2 k = uModTime.xy * alias(vec2(i), vec2(N));
72    float k_len = length(k);
73
74    const float G = 9.81;
75
76    // If this sample runs for hours on end, the cosines of very large numbers will eventually become unstable.
77    // It is fairly easy to fix this by wrapping uTime,
78    // and quantizing w such that wrapping uTime does not change the result.
79    // See Tessendorf's paper for how to do it.
80    // The sqrt(G * k_len) factor represents how fast ocean waves at different frequencies propagate.
81    float w = sqrt(G * k_len) * uModTime.z;
82    float cw = cos(w);
83    float sw = sin(w);
84
85    // Complex multiply to rotate our frequency samples.
86    a = cmul(a, vec2(cw, sw));
87    b = cmul(b, vec2(cw, sw));
88    b = vec2(b.x, -b.y); // Complex conjugate since we picked a frequency with the opposite direction.
89    vec2 res = a + b; // Sum up forward and backwards travelling waves.
90    heights[i.y * N.x + i.x] = pack2(res);
91}
92
93void main()
94{
95    generate_heightmap();
96}
97
98