// Procedural ground-cover blades as a MESH shader (Vulkan, VK_EXT_mesh_shader): the same blades
// grass.vert draws - the same places, density, thinning, culling, size, sway and lighting normal -
// but a blade the tests reject emits nothing, where the instanced path still runs all eight of its
// vertices to a degenerate position.
//
// One dispatch covers a chunk of tiles (grass.ludic): work group y is the tile's place in u_tiles,
// work group x a batch of BLADES blade indices within that tile. Everything past the tile's own
// count is skipped.
layout(local_size_x = 1, local_size_y = 1, local_size_z = 1) in;
const int BLADES = 16;                 // blades an invocation may emit
const int ROWS = 5;                    // grass_blade_mesh(5): two vertices a row, four quads
layout(triangles, max_vertices = 128, max_primitives = 96) out;

uniform mat4 u_view;
uniform mat4 u_proj;
uniform mat4 u_vp;
uniform vec3 u_cam_pos;
uniform float u_time;
uniform sampler2D u_ts_height;
// the baked terrain normal: x and z in RG16F, y rebuilt (a terrain normal always points up)
uniform sampler2D u_ter_normal;
uniform vec2 u_ts_origin;
uniform float u_ts_half;
uniform sampler2D u_ortho;
uniform float u_ortho_on;
// terpage.glsl, pasted (the reference copy and its rules are there)
// ---- the page table ----------------------------------------------------------------------
#ifndef TP_HELPERS
#define TP_HELPERS
uniform sampler2D u_tp_page;   // R32F, a texel a tile: slot + 1, 0 = not resident
uniform float u_tp_on;         // 1 while paging; 0 reads the whole maps exactly as before
uniform vec4 u_tp_dims;        // tiles a side, height texels a tile, photo texels a tile, full height res
// the layer holding a full-map uv, or -1: paging off, off the map, or not resident
float tpSlot(vec2 uv) {
  if (u_tp_on < 0.5 || uv.x < 0.0 || uv.y < 0.0 || uv.x >= 1.0 || uv.y >= 1.0) return -1.0;
  ivec2 t = min(ivec2(floor(uv * u_tp_dims.x)), ivec2(int(u_tp_dims.x) - 1));
  return texelFetch(u_tp_page, t, 0).r - 1.0;
}
// a full-map uv inside its tile's layer of k texels a side, past the one-texel border
vec2 tpUV(vec2 uv, float k) { return (1.0 + fract(uv * u_tp_dims.x) * k) / (k + 2.0); }
#endif
// ---- height (vertex and mesh stages) -----------------------------------------------------
#ifndef TP_HEIGHT
#define TP_HEIGHT
#ifndef TP_HMAP
#define TP_HMAP u_ts_height
#endif
uniform sampler2DArray u_tp_h;
// bilinear height, fine where a tile is resident. Outside a fragment stage texture() IS level 0
// (glslang emits the same explicit-lod sample), so this is the old read bit for bit.
float terHeight(vec2 uv) {
  float s = tpSlot(uv);
  if (s >= 0.0) return textureLod(u_tp_h, vec3(tpUV(uv, u_tp_dims.y), s), 0.0).r;
  return textureLod(TP_HMAP, uv, 0.0).r;
}
// The B-spline through four bilinear taps (heightSmooth). Paging on, the texel is the FULL
// map's - the coarse map's own size would halve the surface's detail - and each tap finds its
// own tile, so a patch straddling two tiles, or a tile and the coarse map, stays one surface.
float terHeightSmooth(vec2 uv) {
  vec2 res = u_tp_on > 0.5 ? vec2(u_tp_dims.w) : vec2(textureSize(TP_HMAP, 0));
  vec2 t = uv * res - 0.5;
  vec2 f = fract(t);
  vec2 i = floor(t);
  vec2 w0 = (1.0 - f) * (1.0 - f) * (1.0 - f) / 6.0;
  vec2 w1 = (4.0 - 6.0 * f * f + 3.0 * f * f * f) / 6.0;
  vec2 w3 = f * f * f / 6.0;
  vec2 w2 = 1.0 - w0 - w1 - w3;
  vec2 s0 = w0 + w1, s1 = w2 + w3;
  vec2 o0 = (i - 1.0 + w1 / s0 + 0.5) / res;
  vec2 o1 = (i + 1.0 + w3 / s1 + 0.5) / res;
  return (terHeight(vec2(o0.x, o0.y)) * s0.x + terHeight(vec2(o1.x, o0.y)) * s1.x) * s0.y
       + (terHeight(vec2(o0.x, o1.y)) * s0.x + terHeight(vec2(o1.x, o1.y)) * s1.x) * s1.y;
}
#endif
// ---- the baked normal (x, z; y is rebuilt by the caller) ----------------------------------
#ifndef TP_NORMAL
#define TP_NORMAL
uniform sampler2DArray u_tp_nrm;
vec2 terNormalXZ(vec2 uv) {
#ifdef TP_FRAGMENT
  if (u_tp_on < 0.5) return texture(u_ter_normal, uv).rg;
  // taken before the per-pixel branch: a neighbour on another tile must not decide the level
  vec2 dx = dFdx(uv), dy = dFdy(uv);
#endif
  float s = tpSlot(uv);
  if (s >= 0.0) return textureLod(u_tp_nrm, vec3(tpUV(uv, u_tp_dims.y), s), 0.0).rg;
#ifdef TP_FRAGMENT
  return textureGrad(u_ter_normal, uv, dx, dy).rg;
#else
  return textureLod(u_ter_normal, uv, 0.0).rg;
#endif
}
#endif
// ---- the photograph's measures (declares no sampler: a stage reading u_ortho blurred needs no tile)
#ifndef TP_ORTHO_LOD
#define TP_ORTHO_LOD
// the FULL photograph's texels a side, for filters that weigh by the texel (orthoSmooth)
vec2 tpOrthoRes() { return u_tp_on > 0.5 ? vec2(u_tp_dims.z * u_tp_dims.x) : vec2(textureSize(u_ortho, 0)); }
// A blurred read (textureLod 1..2.5) names a level of the FULL photograph; on the coarse map
// the same detail is that many levels lower, so the ground's far colour does not soften.
float tpOrthoLod(float lod) {
  if (u_tp_on < 0.5) return lod;
  return max(lod - log2(u_tp_dims.z * u_tp_dims.x / float(textureSize(u_ortho, 0).x)), 0.0);
}
#endif
vec3 terNormal(vec2 uv) { vec2 xz = terNormalXZ(uv);  return vec3(xz.x, sqrt(max(1.0 - dot(xz, xz), 0.0)), xz.y); }
uniform float u_lake_level;
uniform float u_sea_level;
uniform vec4 u_lake;
uniform float u_snow_line;
uniform float u_wind;
uniform int u_tile_cells;    // 16 m cells per tile side
uniform vec4 u_tiles[256];   // per tile of the dispatch: corner x, corner z, indices per cell
uniform float u_s0;
uniform float u_d0;
uniform float u_radius;
uniform float u_px;          // one pixel in radians (grass.vert)
uniform int u_dbg;
uniform vec4 u_gk_ortho;     // the grass kind's (grass_bind_kind): w its cell (m); the rest is grass.vert's
out vec3 v_wpos[];
out vec3 v_nrm[];
out vec2 v_uv[];
out float v_seed[];
out vec2 v_rot[];
out float v_hull[];
out float v_quake[];   // grass does not quake; written so model.frag can read it

// The kind's cell, as the CPU's tiles count it. Nothing else of the kind reaches this path yet: its
// shape, sizes and density test are an older copy of grass.vert's, kept as they draw today.
#define CELL u_gk_ortho.w
float hash1(vec2 p) { return fract(sin(dot(p, vec2(127.1, 311.7))) * 43758.5453123); }
uint pcg(uint v) { uint s = v * 747796405u + 2891336453u; uint w = ((s >> ((s >> 28u) + 4u)) ^ s) * 277803737u; return (w >> 22u) ^ w; }
float bladeHash(ivec2 cell, int j, int k) {
  uint h = pcg(uint(cell.x + 32768) * 73856093u ^ uint(cell.y + 32768) * 19349663u ^ uint(j) * 83492791u ^ uint(k) * 2654435761u);
  return float(h) * (1.0 / 4294967295.0);
}

void main() {
  int ti = int(gl_WorkGroupID.y);
  vec2 tile = u_tiles[ti].xy;
  int per_cell = int(u_tiles[ti].z + 0.5);
  int total = min(per_cell * u_tile_cells * u_tile_cells, 65535);   // the instanced path's cap on a tile
  int first = int(gl_WorkGroupID.x) * BLADES;
  int nv = 0;
  int np = 0;
  for (int b = 0; b < BLADES; b++) {
    int i = first + b;
    if (per_cell <= 0 || i >= total) break;
    int c = i / per_cell;
    int j = i - c * per_cell;
    vec2 cell = tile + vec2(float(c % u_tile_cells), float(c / u_tile_cells)) * CELL;
    vec2 cid = floor(cell / CELL + 0.5);
    ivec2 ci = ivec2(cid);
    float fj = float(j);
    vec2 hv = vec2(bladeHash(ci, j, 0), bladeHash(ci, j, 1));
    vec2 xz = cell + hv * CELL;
    float dist = length(xz - u_cam_pos.xz);
    if (dist >= u_radius) continue;
    float spacing = u_s0 * (1.0 + dist / u_d0);
    float count = CELL * CELL / (spacing * spacing) * (1.0 - smoothstep(u_radius * 0.7, u_radius, dist));
    if (fj >= count) continue;
    float life = 1.0 - smoothstep(0.8, 1.0, fj / max(count, 1.0));
    vec2 huv = (xz - u_ts_origin) / (2.0 * u_ts_half) + 0.5;
    if (huv.x < 0.0 || huv.x > 1.0 || huv.y < 0.0 || huv.y > 1.0) continue;
    vec4 ht = vec4(terHeight(huv));   // only .r is read
    vec4 croot = u_vp * vec4(xz.x, ht.r, xz.y, 1.0);
    if (croot.w < -1.0 || abs(croot.x) > croot.w * 1.25 + 1.5 || abs(croot.y) > croot.w * 1.4 + 1.5) continue;
    vec3 gn = terNormal(huv);
    float h3 = bladeHash(ci, j, 2), h4 = bladeHash(ci, j, 3);
    float wl = u_sea_level;
    if (u_lake.z > 0.0) { vec2 q = (xz - u_lake.xy) / u_lake.zw;  if (dot(q, q) < 1.0) wl = max(wl, u_lake_level); }
    float ok = (1.0 - smoothstep(0.30, 0.55, 1.0 - gn.y)) * smoothstep(0.0, 0.6, ht.r - wl - 0.15) * smoothstep(u_snow_line - 80.0, u_snow_line - 200.0, ht.r);
    if (u_ortho_on > 0.5) {
      vec3 oc = textureLod(u_ortho, huv, tpOrthoLod(1.5)).rgb;
      ok *= 0.25 + 0.75 * smoothstep(0.0, 0.02, oc.g - max(oc.r, oc.b));
    }
    if (h4 > ok) continue;
    bool far = dist > 300.0;
    float h = far ? ht.r : terHeightSmooth(huv);
    if (far) h += 0.03;
    if ((u_dbg & 1) != 0) h += 0.3;
    float seed = hv.x * 0.7 + hv.y * 0.3;
    float ang = hv.y * 6.2831853;
    float s = sin(ang), c_ = cos(ang);
    float grow = spacing / u_s0;
    float tall = mix(0.18, 0.42, h3) * mix(0.8, 1.2, hash1(cid * 0.1)) * (1.0 + 0.35 * smoothstep(1.0, 12.0, grow)) * life;
    float bw = 0.028 * mix(1.0, 0.45 * grow, smoothstep(1.0, 4.0, grow));
    if (far) { bw = max(bw, spacing * 0.35); tall = min(tall, spacing * 0.3); }
    bw = max(bw, dist * u_px * 1.1);
    float gust = sin(xz.x * 0.09 + u_time * 1.1) * 0.5 + sin(xz.y * 0.13 - u_time * 0.8 + xz.x * 0.05) * 0.5;
    float ph = u_time * 1.7 + seed * 6.2831 + xz.x * 0.05 + xz.y * 0.07;
    float sway = (sin(ph) * 0.6 + sin(ph * 2.3 + 1.0) * 0.4 + gust) * u_wind;
    // the ground's frame, shared by the blade's vertices
    vec3 up = vec3(0.0, 1.0, 0.0);
    vec3 k = cross(up, gn);
    float sk = length(k), ck = gn.y;
    bool tilt = sk > 1e-4;
    if (tilt) k /= sk;
    vec3 n = vec3(0.0, 0.3, 1.0);
    n = normalize(vec3(c_ * n.x + s * n.z, n.y, -s * n.x + c_ * n.z));
    if (tilt) n = normalize(n * ck + cross(k, n) * sk + k * dot(k, n) * (1.0 - ck));
    n = normalize(mix(n, gn, smoothstep(2.0, 12.0, dist)));
    float hull = (dist > 2.0 || far || (u_dbg & 2) != 0) ? -1.0 : 1.0;
    // the blade: grass_blade_mesh(5)'s vertices, placed as grass.vert places them
    for (int r = 0; r < ROWS; r++) {
      float t = float(r) / float(ROWS - 1);
      // A COPY of grass_blade_mesh's profile in grass.ludic - narrow at the sheath,
      // widest a fifth of the way up, a fine point, and arched. Change both together.
      float taper = max((0.50 + 0.33 * min(t / 0.22, 1.0)) * (1.0 - t * t * t), 0.05);
      float bend = t * t * 0.52;
      for (int sd = 0; sd < 2; sd++) {
        vec3 a_pos = vec3((float(sd) - 0.5) * taper, t, bend);
        vec2 a_uv = vec2(float(sd), t);
        vec3 p = vec3(a_pos.x * bw, a_pos.y * tall, a_pos.z * tall * (0.6 + 0.8 * h4));
        float hgt = max(p.y, 0.0);
        p.x += sway * hgt * hgt * 0.35;
        p.z += sway * hgt * hgt * 0.15 * cos(ph * 0.7);
        p = vec3(c_ * p.x + s * p.z, p.y, -s * p.x + c_ * p.z);
        if (tilt) p = p * ck + cross(k, p) * sk + k * dot(k, p) * (1.0 - ck);
        vec3 w = vec3(xz.x, h - 0.02, xz.y) + p;
        vec4 clip = u_proj * u_view * vec4(w, 1.0);
        clip.z = (clip.z + clip.w) * 0.5;               // OpenGL's depth range to Vulkan's, as the vertex wrapper does
        int o = nv + r * 2 + sd;
        gl_MeshVerticesEXT[o].gl_Position = clip;
        v_wpos[o] = w;
        v_nrm[o] = n;
        v_uv[o] = far ? vec2(a_uv.x, 0.45 + 0.2 * a_uv.y) : a_uv;
        v_seed[o] = seed;
        v_rot[o] = vec2(s, c_);
        v_hull[o] = hull;
        v_quake[o] = 0.0;
      }
    }
    for (int q = 0; q < ROWS - 1; q++) {
      uint bb = uint(nv + q * 2);
      gl_PrimitiveTriangleIndicesEXT[np] = uvec3(bb, bb + 1u, bb + 2u);
      gl_PrimitiveTriangleIndicesEXT[np + 1] = uvec3(bb + 1u, bb + 3u, bb + 2u);
      np += 2;
    }
    nv += ROWS * 2;
  }
  SetMeshOutputsEXT(uint(nv), uint(np));
}