shaders/comp/generate_height.comp

   1 #version 310 es
   2
   3 layout(local_size_x = 64) in;
   4
   5 layout(std430, binding = 0) readonly buffer Distribution
   6 {
   7     vec2 distribution[];
   8 };
   9
  10 layout(std430, binding = 1) writeonly buffer HeightmapFFT
  11 {
  12     uint heights[];
  13 };
  14
  15 layout(binding = 2, std140) uniform UBO
  16 {
  17     vec4 uModTime;
  18 };
  19
  20 vec2 alias(vec2 i, vec2 N)
  21 {
  22     return mix(i, i - N, greaterThan(i, 0.5 * N));
  23 }
  24
  25 vec4 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
  35 vec2 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
  45 uint pack2(vec2 v)
  46 {
  47     return packHalf2x16(v);
  48 }
  49
  50 uvec2 pack4(vec4 v)
  51 {
  52     return uvec2(packHalf2x16(v.xy), packHalf2x16(v.zw));
  53 }
  54
  55 uvec2 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
  60 void 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
  93 void main()
  94 {
  95     generate_heightmap();
  96 }
  97