map/mapTerrain.js

import { INF, MAP_H, MAP_W, SIZE, MinHeap, clamp, createMapFields, fbm, hash2, indexOf, inside, lerp, pickEntities, rand, smoothstep, valueNoise, xyOf } from "./mapUtils.js";
import {
  extractMaskBorder,
  extractRegionBorderSegments,
  generateRegionalPrefectures,
  makePrefectureMask,
  neighbors8,
} from "./mapGeneratorHelpers.js";

const ASPECT = MAP_W / MAP_H;
const SQRT2 = Math.SQRT2;

function normalizeCoord(x, y) {
  return {
    px: (x + 0.5) / MAP_W,
    py: (y + 0.5) / MAP_H,
  };
}

function distNorm(ax, ay, bx, by) {
  const dx = (ax - bx) * ASPECT;
  const dy = ay - by;
  return Math.hypot(dx, dy);
}

function rotate(dx, dy, angle) {
  const c = Math.cos(angle);
  const s = Math.sin(angle);
  return { u: dx * c + dy * s, v: -dx * s + dy * c };
}

function quantile(values, q) {
  const arr = Array.from(values).filter(Number.isFinite).sort((a, b) => a - b);
  if (!arr.length) return 0;
  const p = clamp(q) * (arr.length - 1);
  const i = Math.floor(p);
  const f = p - i;
  return lerp(arr[i], arr[Math.min(arr.length - 1, i + 1)], f);
}

function forDisk(cx, cy, radius, fn) {
  const r = Math.ceil(radius);
  for (let dy = -r; dy <= r; dy++) {
    for (let dx = -r; dx <= r; dx++) {
      const x = cx + dx;
      const y = cy + dy;
      if (!inside(x, y)) continue;
      const d = Math.hypot(dx, dy);
      if (d <= radius) fn(x, y, d);
    }
  }
}

function largestComponent(mask, allowEdgePreference = false) {
  const seen = new Uint8Array(SIZE);
  let best = [];
  let bestScore = -1;
  for (let i = 0; i < SIZE; i++) {
    if (!mask[i] || seen[i]) continue;
    const queue = [i];
    const cells = [];
    let touchesEdge = false;
    seen[i] = 1;
    for (let q = 0; q < queue.length; q++) {
      const cur = queue[q];
      cells.push(cur);
      const x = cur % MAP_W;
      const y = Math.floor(cur / MAP_W);
      if (x === 0 || y === 0 || x === MAP_W - 1 || y === MAP_H - 1) touchesEdge = true;
      for (const [nx, ny] of neighbors8(x, y)) {
        const ni = indexOf(nx, ny);
        if (!mask[ni] || seen[ni]) continue;
        seen[ni] = 1;
        queue.push(ni);
      }
    }
    const score = cells.length + (allowEdgePreference && touchesEdge ? SIZE : 0);
    if (score > bestScore) {
      bestScore = score;
      best = cells;
    }
  }
  const out = new Uint8Array(SIZE);
  for (const i of best) out[i] = 1;
  return out;
}

function distanceField(sourceMask, maxDistance = 999) {
  const dist = new Float32Array(SIZE);
  dist.fill(maxDistance);
  const heap = new MinHeap();
  for (let i = 0; i < SIZE; i++) {
    if (!sourceMask[i]) continue;
    dist[i] = 0;
    heap.push({ i, f: 0 });
  }
  while (heap.length) {
    const cur = heap.pop();
    if (!cur || cur.f > dist[cur.i] + 1e-5) continue;
    const x = cur.i % MAP_W;
    const y = Math.floor(cur.i / MAP_W);
    for (const [nx, ny] of neighbors8(x, y)) {
      const ni = indexOf(nx, ny);
      const step = (nx !== x && ny !== y) ? SQRT2 : 1;
      const nd = cur.f + step;
      if (nd >= dist[ni] || nd > maxDistance) continue;
      dist[ni] = nd;
      heap.push({ i: ni, f: nd });
    }
  }
  return dist;
}

function ridgeContribution(px, py, ridge, seed) {
  const dx = (px - ridge.x) * ASPECT;
  const dy = py - ridge.y;
  const { u, v } = rotate(dx, dy, ridge.angle);
  const half = ridge.length * 0.5;
  const along = Math.abs(u / Math.max(0.001, half));
  if (along >= 1.22) return 0;
  const taper = smoothstep(1 - clamp((along - 0.68) / 0.54));
  const wobble = (valueNoise((u + ridge.phase) * 720, (py + ridge.phase) * 720, seed + ridge.seedOffset, 14) - 0.5) * ridge.width * ridge.wobble;
  const cross = Math.abs(v + wobble);
  const core = Math.exp(-Math.pow(cross / Math.max(0.0008, ridge.width), 2.0));
  const serration = 0.82 + 0.36 * valueNoise((px + ridge.phase) * 900, (py - ridge.phase) * 900, seed + ridge.seedOffset + 71, 7.5);
  return ridge.height * core * taper * serration;
}

function ellipticalMask(px, py, system) {
  const dx = (px - system.x) * ASPECT;
  const dy = py - system.y;
  const { u, v } = rotate(dx, dy, system.angle);
  const a = Math.max(0.01, system.length * 0.5);
  const b = Math.max(0.01, system.width * 0.5);
  const r = Math.sqrt((u / a) ** 2 + (v / b) ** 2);
  return clamp(1 - smoothstep((r - 0.55) / 0.65));
}

function sampleInsideUnitDisk(seed, n) {
  const r = Math.sqrt(rand(seed, n));
  const a = rand(seed, n + 1) * Math.PI * 2;
  return { x: Math.cos(a) * r, y: Math.sin(a) * r, r };
}

export function buildTerrainTemplate(seed) {
  const mountainModeRoll = rand(seed, 12);
  const mountainMode = mountainModeRoll < 0.48 ? "range" : mountainModeRoll < 0.80 ? "mixed" : "massif";
  const mountainMassifness = mountainMode === "massif" ? 0.72 + rand(seed, 13) * 0.24 : mountainMode === "mixed" ? 0.34 + rand(seed, 14) * 0.36 : rand(seed, 15) * 0.24;
  const coastAngle = rand(seed, 21) * Math.PI * 2;
  const twoSidedCoast = rand(seed, 22) < 0.36;
  const seaRatio = 0.14 + rand(seed, 23) * 0.16;
  const mountainAngle = coastAngle + Math.PI * (0.26 + rand(seed, 24) * 0.48);
  const baseHeight = 0.52 + rand(seed, 25) * 0.46;
  const primaryLength = lerp(0.70 + rand(seed, 26) * 0.22, 0.38 + rand(seed, 27) * 0.20, mountainMassifness);
  const primaryWidth = lerp(0.15 + rand(seed, 28) * 0.13, 0.36 + rand(seed, 29) * 0.20, mountainMassifness);
  const scratchCount = Math.round(lerp(26 + rand(seed, 30) * 22, 18 + rand(seed, 31) * 18, mountainMassifness));
  // 脊梁山脈そのものを複数箇所に置く。旧版の secondary は主山脈の周囲に寄りすぎ、
  // 画面上では「単一の山塊」に見えやすかったため、独立した major system として扱う。
  const mountainSystemCount = 14 + Math.floor(rand(seed, 32) * 3); // 14〜16
  return {
    seed,
    seaRatio,
    coastAngle,
    twoSidedCoast,
    coastNoise: 0.045 + rand(seed, 33) * 0.045,
    mountainMode,
    mountainMassifness,
    mountainAngle,
    mountainBaseHeight: baseHeight,
    mountainDensity: 0.62 + rand(seed, 34) * 0.35,
    primaryMountain: {
      x: clamp(0.50 + (rand(seed, 35) - 0.5) * 0.28, 0.22, 0.78),
      y: clamp(0.50 + (rand(seed, 36) - 0.5) * 0.28, 0.22, 0.78),
      angle: mountainAngle,
      length: primaryLength,
      width: primaryWidth,
      height: baseHeight,
      scratchCount,
      massifness: mountainMassifness,
    },
    mountainSystemCount,
    secondaryCount: mountainSystemCount - 1,
    macroNoiseScale: 0.020 + rand(seed, 37) * 0.030,
    macroNoiseStrength: 0.028 + rand(seed, 38) * 0.025,
    scratchNoiseStrength: 0.018 + rand(seed, 39) * 0.022,
    roughness: 0.40 + rand(seed, 40) * 0.50,
    erosion: 0.34 + rand(seed, 41) * 0.48,
    deposition: 0.28 + rand(seed, 42) * 0.56,
    riverRichness: 0.72 + rand(seed, 43) * 0.60,
    bigRiverChance: 0.34 + rand(seed, 44) * 0.34,
    plainBias: 0.32 + rand(seed, 45) * 0.46,
  };
}

function buildMountainSystems(template, seed) {
  const systems = [];
  const targetCount = Math.max(8, template.mountainSystemCount ?? 15);
  const baseAngle = template.mountainAngle;

  // 複数の脊梁山脈システムを、画面中央ではなくマップ全域に分散配置する。
  // 5x3 / 4x4 に近い粗い格子へ jitter を入れ、さらに farthest-candidate で
  // 既存システムから離れた候補を選ぶ。これにより「中央に単一山塊」化しにくくする。
  const cols = targetCount >= 14 ? 5 : 4;
  const rows = Math.ceil(targetCount / cols);
  const cellOrder = Array.from({ length: cols * rows }, (_, i) => i)
    .map((v) => ({ v, key: rand(seed, 1000 + v * 17) }))
    .sort((a, b) => a.key - b.key)
    .map((o) => o.v);

  function gridCandidate(k, attempt) {
    const cell = cellOrder[(k + attempt * 7) % cellOrder.length];
    const cx = cell % cols;
    const cy = Math.floor(cell / cols);
    const jitterX = (rand(seed, 1100 + k * 101 + attempt * 13) - 0.5) * 0.62;
    const jitterY = (rand(seed, 1200 + k * 101 + attempt * 13) - 0.5) * 0.62;
    const x = clamp((cx + 0.5 + jitterX) / cols, 0.055, 0.945);
    const y = clamp((cy + 0.5 + jitterY) / rows, 0.055, 0.945);
    const localTurn = (rand(seed, 1300 + k * 101 + attempt) - 0.5) * Math.PI * 0.92;
    const diagonalBias = (cx / Math.max(1, cols - 1) - 0.5 + (cy / Math.max(1, rows - 1) - 0.5) * 0.35) * 0.16;
    return {
      x,
      y,
      angle: baseAngle + localTurn + diagonalBias,
    };
  }

  function randomCandidate(k, attempt) {
    return {
      x: clamp(0.055 + rand(seed, 2000 + k * 137 + attempt * 31) * 0.89, 0.055, 0.945),
      y: clamp(0.055 + rand(seed, 2100 + k * 137 + attempt * 31) * 0.89, 0.055, 0.945),
      angle: baseAngle + (rand(seed, 2200 + k * 137 + attempt) - 0.5) * Math.PI * 1.05,
    };
  }

  function candidateAt(k, attempt) {
    return attempt < 5 ? gridCandidate(k, attempt) : randomCandidate(k, attempt);
  }

  for (let k = 0; k < targetCount; k++) {
    let best = candidateAt(k, 0);
    let bestScore = -INF;
    for (let attempt = 0; attempt < 18; attempt++) {
      const c = candidateAt(k, attempt);
      let minD = 999;
      for (const s of systems) minD = Math.min(minD, distNorm(c.x, c.y, s.x, s.y));
      // 中央集中を避けるため、中心距離を少し加点する。ただし端に張り付きすぎないよう edge も見る。
      const edgeD = Math.min(c.x, c.y, 1 - c.x, 1 - c.y);
      const centerD = distNorm(c.x, c.y, 0.5, 0.5);
      const score =
        minD * 1.25 +
        centerD * 0.18 +
        Math.min(edgeD, 0.16) * 0.22 +
        rand(seed, 2300 + k * 101 + attempt) * 0.04;
      if (score > bestScore) { bestScore = score; best = c; }
    }

    const m = clamp(template.mountainMassifness + (rand(seed, 2400 + k) - 0.5) * 0.50);
    const isMassif = m > 0.58;
    const major = k < 4 || rand(seed, 2500 + k) > 0.68;
    const length = lerp(
      major ? 0.30 + rand(seed, 2600 + k) * 0.22 : 0.20 + rand(seed, 2610 + k) * 0.16,
      major ? 0.22 + rand(seed, 2620 + k) * 0.14 : 0.16 + rand(seed, 2630 + k) * 0.12,
      m
    );
    const width = lerp(
      major ? 0.055 + rand(seed, 2700 + k) * 0.060 : 0.040 + rand(seed, 2710 + k) * 0.045,
      major ? 0.120 + rand(seed, 2720 + k) * 0.090 : 0.085 + rand(seed, 2730 + k) * 0.070,
      m
    );
    const height = template.mountainBaseHeight * (
      major
        ? 0.34 + rand(seed, 2800 + k) * 0.24
        : 0.20 + rand(seed, 2810 + k) * 0.18
    );
    const scratchCount = Math.round(lerp(
      major ? 10 + rand(seed, 2900 + k) * 10 : 6 + rand(seed, 2910 + k) * 7,
      isMassif ? 8 + rand(seed, 2920 + k) * 9 : 6 + rand(seed, 2930 + k) * 7,
      m
    ));

    systems.push({
      x: best.x,
      y: best.y,
      angle: best.angle,
      length,
      width,
      height,
      scratchCount,
      massifness: m,
      role: major ? (k < 4 ? "primary" : "major") : "minor",
    });
  }

  return systems;
}

function buildScratchRidges(system, seed, systemId) {
  const ridges = [];
  const count = Math.max(6, Math.round(system.scratchCount));
  for (let i = 0; i < count; i++) {
    const p = sampleInsideUnitDisk(seed + systemId * 10000, 2000 + i * 7);
    const density = clamp(1 - p.r * 0.78);
    const localAngle = system.massifness > 0.55
      ? system.angle + (rand(seed, 2100 + i + systemId * 331) - 0.5) * Math.PI * 1.45
      : system.angle + (rand(seed, 2100 + i + systemId * 331) - 0.5) * (0.36 + system.massifness * 0.80);
    const along = p.x * system.length * 0.45;
    const cross = p.y * system.width * 0.45;
    const x = clamp(system.x + Math.cos(system.angle) * along / ASPECT + Math.cos(system.angle + Math.PI / 2) * cross / ASPECT, 0.03, 0.97);
    const y = clamp(system.y + Math.sin(system.angle) * along + Math.sin(system.angle + Math.PI / 2) * cross, 0.03, 0.97);
    const len = lerp(system.length * (0.18 + rand(seed, 2200 + i) * 0.20), system.width * (0.32 + rand(seed, 2200 + i) * 0.30), system.massifness);
    const width = lerp(0.010 + rand(seed, 2300 + i) * 0.012, 0.018 + rand(seed, 2300 + i) * 0.020, system.massifness) * (0.80 + density * 0.60);
    const height = system.height * (0.040 + density * 0.095 + rand(seed, 2400 + i) * 0.035);
    ridges.push({
      x, y,
      angle: localAngle,
      length: len,
      width,
      height,
      wobble: 1.2 + rand(seed, 2500 + i) * 2.0,
      phase: rand(seed, 2600 + i) * 10,
      seedOffset: 2700 + systemId * 997 + i * 37,
      density,
      systemId,
    });
  }
  return ridges;
}

function computeCoastLower(px, py, template, seed) {
  const axis = (px - 0.5) * Math.cos(template.coastAngle) * ASPECT + (py - 0.5) * Math.sin(template.coastAngle);
  const wave = (fbm(px * 220, py * 220, seed + 300) - 0.5) * template.coastNoise;
  const bay = (valueNoise(px * 500, py * 500, seed + 301, 22) - 0.5) * 0.055;
  const sideA = smoothstep((-axis + 0.24 + wave + bay) / 0.26);
  const sideB = template.twoSidedCoast ? smoothstep((axis + 0.20 - wave + bay * 0.7) / 0.27) : 0;
  const pressure = Math.max(sideA, sideB);
  return { pressure, signedAxis: axis };
}

function recomputeSlope(elevation, sea, slope) {
  slope.fill(0);
  for (let y = 1; y < MAP_H - 1; y++) {
    for (let x = 1; x < MAP_W - 1; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      const gx = elevation[indexOf(x + 1, y)] - elevation[indexOf(x - 1, y)];
      const gy = elevation[indexOf(x, y + 1)] - elevation[indexOf(x, y - 1)];
      slope[i] = clamp(Math.hypot(gx, gy) * 8.2);
    }
  }
}

function classifyWater(elevation, seaLevel, sea, ocean, lake) {
  sea.fill(0); ocean.fill(0); lake.fill(0);
  const water = new Uint8Array(SIZE);
  for (let i = 0; i < SIZE; i++) water[i] = elevation[i] <= seaLevel ? 1 : 0;
  const oceanMask = largestComponent(water, true);
  const seen = new Uint8Array(SIZE);
  for (let i = 0; i < SIZE; i++) {
    if (!water[i] || seen[i]) continue;
    const queue = [i];
    const cells = [];
    seen[i] = 1;
    for (let q = 0; q < queue.length; q++) {
      const cur = queue[q];
      cells.push(cur);
      const x = cur % MAP_W;
      const y = Math.floor(cur / MAP_W);
      for (const [nx, ny] of neighbors8(x, y)) {
        const ni = indexOf(nx, ny);
        if (!water[ni] || seen[ni]) continue;
        seen[ni] = 1;
        queue.push(ni);
      }
    }
    const isOcean = cells.some((ci) => oceanMask[ci]);
    if (isOcean || cells.length >= 22) {
      for (const ci of cells) {
        sea[ci] = 1;
        if (isOcean) ocean[ci] = 1;
        else lake[ci] = 1;
      }
    } else {
      for (const ci of cells) elevation[ci] = seaLevel + 0.010;
    }
  }
}

function priorityFloodFlow(elevation, sea, flowTo, filled) {
  flowTo.fill(-1);
  filled.set(elevation);
  const visited = new Uint8Array(SIZE);
  const heap = new MinHeap();
  let seedCount = 0;
  for (let i = 0; i < SIZE; i++) {
    if (sea[i]) {
      visited[i] = 1;
      heap.push({ i, f: filled[i] });
      seedCount++;
    }
  }
  if (seedCount === 0) {
    for (let i = 0; i < SIZE; i++) {
      const x = i % MAP_W;
      const y = Math.floor(i / MAP_W);
      if (x === 0 || y === 0 || x === MAP_W - 1 || y === MAP_H - 1) {
        visited[i] = 1;
        heap.push({ i, f: filled[i] });
      }
    }
  }
  while (heap.length) {
    const cur = heap.pop();
    if (!cur) continue;
    const cx = cur.i % MAP_W;
    const cy = Math.floor(cur.i / MAP_W);
    for (const [nx, ny] of neighbors8(cx, cy)) {
      const ni = indexOf(nx, ny);
      if (visited[ni]) continue;
      visited[ni] = 1;
      if (filled[ni] < filled[cur.i] + 0.00002) filled[ni] = filled[cur.i] + 0.00002;
      heap.push({ i: ni, f: filled[ni] });
    }
  }

  for (let y = 0; y < MAP_H; y++) {
    for (let x = 0; x < MAP_W; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      let best = -1;
      let bestScore = filled[i];
      for (const [nx, ny] of neighbors8(x, y)) {
        const ni = indexOf(nx, ny);
        const stepPenalty = (nx !== x && ny !== y) ? 0.000015 : 0;
        const score = filled[ni] + stepPenalty + hash2(nx, ny, 9000) * 0.000002;
        if (score < bestScore - 0.000001 || sea[ni]) {
          bestScore = score;
          best = ni;
          if (sea[ni]) break;
        }
      }
      flowTo[i] = best;
    }
  }
}

function computeFlowAccumulation(sea, flowTo, filled, flowAccum) {
  const area = new Float32Array(SIZE);
  const order = [];
  for (let i = 0; i < SIZE; i++) {
    if (sea[i]) continue;
    area[i] = 1;
    order.push(i);
  }
  order.sort((a, b) => filled[b] - filled[a]);
  for (const i of order) {
    const to = flowTo[i];
    if (to >= 0 && !sea[to]) area[to] += area[i];
  }
  let maxArea = 1;
  for (let i = 0; i < SIZE; i++) if (!sea[i]) maxArea = Math.max(maxArea, area[i]);
  for (let i = 0; i < SIZE; i++) flowAccum[i] = sea[i] ? 0 : clamp(Math.pow(area[i] / maxArea, 0.42));
  return area;
}

function traceFlowPath(start, sea, flowTo, maxSteps = 900) {
  const path = [];
  const seen = new Set();
  let i = start;
  for (let step = 0; step < maxSteps && i >= 0 && !seen.has(i); step++) {
    seen.add(i);
    const x = i % MAP_W;
    const y = Math.floor(i / MAP_W);
    path.push([x, y]);
    if (sea[i]) break;
    const next = flowTo[i];
    if (next < 0 || next === i) break;
    i = next;
  }
  return path;
}

function buildRiverNetwork(seed, template, sea, lake, elevation, slope, flowTo, flowAccum, river, erosionField) {
  river.fill(0);
  const candidates = [];
  for (let y = 2; y < MAP_H - 2; y++) {
    for (let x = 2; x < MAP_W - 2; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      const f = flowAccum[i];
      const high = clamp((elevation[i] - 0.36) * 2.2);
      const wet = valueNoise(x * 1.7, y * 1.7, seed + 12000, 24);
      const score = f * 0.80 + high * 0.25 + wet * 0.16 - slope[i] * 0.10;
      if (score > 0.24) candidates.push({ x, y, score });
    }
  }
  const desired = 44 + Math.floor(template.riverRichness * 28);
  const sources = pickEntities(candidates, { max: desired, minDistance: 6, threshold: 0.26, seed: seed + 12100, jitter: 0.035 });
  const riverPaths = [];
  for (const s of sources) {
    const path = traceFlowPath(indexOf(s.x, s.y), sea, flowTo);
    if (path.length >= 8 && path.some(([x, y], k) => k > 5 && (sea[indexOf(x, y)] || lake[indexOf(x, y)]))) riverPaths.push(path);
    else if (path.length >= 14) riverPaths.push(path);
  }

  const longPaths = riverPaths
    .map((path) => {
      let maxFlow = 0;
      let meanFlow = 0;
      for (const [x, y] of path) {
        const f = flowAccum[indexOf(x, y)];
        maxFlow = Math.max(maxFlow, f);
        meanFlow += f;
      }
      meanFlow /= Math.max(1, path.length);
      return { path, score: path.length * 0.75 + maxFlow * 90 + meanFlow * 35 };
    })
    .sort((a, b) => b.score - a.score);
  const mainCount = Math.min(longPaths.length, template.bigRiverChance > 0.56 ? 4 : 3);
  const mainSet = new Set(longPaths.slice(0, mainCount).map((p) => p.path));

  const riverThreshold = 0.26 - template.riverRichness * 0.035 - (template.bigRiverChance > 0.56 ? 0.030 : 0);
  for (let i = 0; i < SIZE; i++) {
    if (sea[i]) continue;
    const f = flowAccum[i];
    if (f > riverThreshold) river[i] = clamp((f - riverThreshold) / (0.55 - riverThreshold));
  }
  for (const path of riverPaths) {
    const main = mainSet.has(path);
    for (let k = 0; k < path.length; k++) {
      const [x, y] = path[k];
      const i = indexOf(x, y);
      if (sea[i]) continue;
      const downstream = k / Math.max(1, path.length - 1);
      const boost = main ? 0.54 + downstream * 0.42 : 0.30 + downstream * 0.22;
      river[i] = clamp(Math.max(river[i], boost + flowAccum[i] * (main ? 0.70 : 0.42)));
    }
  }

  const riverCells = [];
  for (let i = 0; i < SIZE; i++) if (!sea[i] && river[i] > 0.12) riverCells.push(i);
  for (const i of riverCells) {
    const x = i % MAP_W;
    const y = Math.floor(i / MAP_W);
    const strength = river[i];
    const width = 1.15 + clamp((strength - 0.35) * 2.2) * 1.65;
    const depth = (0.006 + strength * (0.014 + template.erosion * 0.018) + flowAccum[i] * 0.018) * (0.65 + clamp((elevation[i] - 0.28) * 1.4) * 0.55);
    forDisk(x, y, width, (nx, ny, d) => {
      const ni = indexOf(nx, ny);
      if (sea[ni]) return;
      const profile = Math.pow(Math.max(0, 1 - d / Math.max(0.1, width)), 2.15);
      const sideGuard = d === 0 ? 1 : 0.38;
      const cut = depth * profile * sideGuard;
      const floor = 0.075;
      elevation[ni] = Math.max(floor, elevation[ni] - cut);
      erosionField[ni] = clamp(erosionField[ni] + cut * 9.0);
    });
  }

  const sortedPaths = longPaths.map((p) => p.path);
  const mainRivers = sortedPaths.filter((p) => mainSet.has(p)).slice(0, mainCount);
  const tributaryRivers = sortedPaths.filter((p) => !mainSet.has(p)).slice(0, 24);
  const smallStreams = sortedPaths.slice(mainCount + 8, mainCount + 58);
  return { riverPaths: sortedPaths.slice(0, 80), mainRivers, tributaryRivers, smallStreams };
}

function deriveFields(seed, template, fields, seaLevel) {
  const {
    elevation, moisture, slope, sea, ocean, lake, river, floodplain, plain, agriculture,
    ridgeField, valleyField, visibleRavineField, surfaceTextureField, basinField,
    coastalLowland, flowAccum, erosionField, depositionField, arcSpineField,
    branchRidgeField, depositionalLowland, alluvialFanField, deltaField,
    naturalBarrierScore, portSuitability, crossingSuitability, passSuitability,
  } = fields;

  recomputeSlope(elevation, sea, slope);
  const waterDist = distanceField(sea, 80);
  const riverMask = new Uint8Array(SIZE);
  for (let i = 0; i < SIZE; i++) if (river[i] > 0.18) riverMask[i] = 1;
  const riverDist = distanceField(riverMask, 40);

  for (let y = 1; y < MAP_H - 1; y++) {
    for (let x = 1; x < MAP_W - 1; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      const e = elevation[i];
      const mean4 = (elevation[indexOf(x - 1, y)] + elevation[indexOf(x + 1, y)] + elevation[indexOf(x, y - 1)] + elevation[indexOf(x, y + 1)]) * 0.25;
      const relief = e - mean4;
      const coast = clamp((24 - waterDist[i]) / 24) * clamp((0.43 - e) * 2.4) * clamp((0.42 - slope[i]) * 2.4);
      coastalLowland[i] = coast;
      const riverNear = clamp((8 - riverDist[i]) / 8);
      const valley = clamp(flowAccum[i] * 0.82 + river[i] * 0.72 + Math.max(0, -relief) * 10.0 + riverNear * 0.38 - slope[i] * 0.22);
      valleyField[i] = clamp(Math.max(valleyField[i] * 0.30, valley));
      const ridge = clamp(arcSpineField[i] * 0.76 + branchRidgeField[i] * 0.86 + Math.max(0, relief) * 9.0 + slope[i] * 0.30 + Math.max(0, e - 0.54) * 0.88 - valleyField[i] * 0.34);
      ridgeField[i] = clamp(Math.max(ridgeField[i] * 0.30, ridge));
      basinField[i] = clamp((0.42 - slope[i]) * 1.55 + Math.max(0, -relief) * 5.0 + clamp((0.48 - e) * 1.35) - coast * 0.40 - river[i] * 0.32);
      const low = clamp((0.58 - e) * 1.55);
      const lowSlope = clamp((0.34 - slope[i]) * 2.8);
      const riverGate = clamp(riverNear * 0.72 + flowAccum[i] * 0.82 + coast * 0.70 + basinField[i] * 0.20 - ridgeField[i] * 0.40);
      plain[i] = clamp(low * lowSlope * (0.18 + template.plainBias * 0.42 + riverGate * 0.92));
      floodplain[i] = clamp(lowSlope * riverNear * (0.35 + flowAccum[i] * 0.82 + river[i] * 0.52));
      deltaField[i] = clamp(coast * river[i] * 1.4 + coast * flowAccum[i] * 0.72);
      alluvialFanField[i] = clamp(riverNear * clamp(slope[i] * 2.5) * clamp((0.58 - e) * 1.7) * clamp(ridgeField[i] * 0.8 + arcSpineField[i] * 0.4));
      depositionalLowland[i] = clamp(floodplain[i] * 0.54 + deltaField[i] * 0.66 + alluvialFanField[i] * 0.42 + coast * 0.28 + basinField[i] * 0.18);
      depositionField[i] = clamp(depositionalLowland[i] * (0.25 + template.deposition * 0.30));
      agriculture[i] = clamp(plain[i] * 0.72 + floodplain[i] * 0.42 + depositionalLowland[i] * 0.36 - slope[i] * 0.20);
      visibleRavineField[i] = clamp(visibleRavineField[i] + valleyField[i] * 0.22 + erosionField[i] * 0.35);
      surfaceTextureField[i] = clamp(surfaceTextureField[i] + slope[i] * 0.26 + visibleRavineField[i] * 0.38 + Math.max(0, relief) * 2.2);
      naturalBarrierScore[i] = clamp(ridgeField[i] * 0.82 + slope[i] * 0.44 + river[i] * 0.42 + Math.max(0, e - 0.58) * 0.34 - plain[i] * 0.24);
      portSuitability[i] = clamp(coast * (1 - slope[i]) * (0.50 + plain[i] * 0.32) - ridgeField[i] * 0.20);
      crossingSuitability[i] = clamp((1 - slope[i]) * 0.38 + plain[i] * 0.34 + floodplain[i] * 0.24 - river[i] * 0.20 - ridgeField[i] * 0.28);
      passSuitability[i] = clamp(slope[i] * 0.30 + valleyField[i] * 0.34 + clamp((0.75 - ridgeField[i]) * 0.8) + plain[i] * 0.18);
      moisture[i] = clamp(0.30 + (1 - waterDist[i] / 65) * 0.36 + valleyField[i] * 0.22 + riverNear * 0.26 - Math.max(0, e - 0.55) * 0.36 + (fbm(x * 1.6, y * 1.6, seed + 15000) - 0.5) * 0.18);
    }
  }
  for (let i = 0; i < SIZE; i++) {
    if (!sea[i]) continue;
    moisture[i] = 1;
    slope[i] = 0;
    river[i] = 0;
    ridgeField[i] = 0;
    valleyField[i] = 0;
    plain[i] = 0;
    agriculture[i] = 0;
    coastalLowland[i] = 0;
    naturalBarrierScore[i] = 0;
  }
}

function enforceLandGradient(elevation, sea, seaLevel) {
  // Keep extreme cliffs rare without flattening normal mountain relief.
  for (let pass = 0; pass < 2; pass++) {
    const next = new Float32Array(elevation);
    for (let y = 1; y < MAP_H - 1; y++) {
      for (let x = 1; x < MAP_W - 1; x++) {
        const i = indexOf(x, y);
        if (sea[i]) continue;
        let minN = elevation[i];
        let maxN = elevation[i];
        for (const [nx, ny] of neighbors8(x, y)) {
          const ni = indexOf(nx, ny);
          if (sea[ni]) continue;
          minN = Math.min(minN, elevation[ni]);
          maxN = Math.max(maxN, elevation[ni]);
        }
        const range = maxN - minN;
        if (range > 0.34) next[i] = lerp(elevation[i], (elevation[i] + minN + maxN) / 3, 0.22);
        next[i] = Math.max(next[i], seaLevel + 0.006);
      }
    }
    elevation.set(next);
  }
}

export function generateTerrainAndRivers(seed) {
  const fields = createMapFields();
  fields.visibleRavineField = new Float32Array(SIZE);
  fields.surfaceTextureField = new Float32Array(SIZE);
  const {
    elevation, moisture, slope, sea, ocean, lake, river, floodplain, plain,
    agriculture, ridgeField, valleyField, visibleRavineField, surfaceTextureField,
    basinField, coastalLowland, flowAccum, erosionField, depositionField,
    arcSpineField, branchRidgeField, depositionalLowland, alluvialFanField,
    deltaField, naturalBarrierScore, flowTo, portSuitability, crossingSuitability,
    passSuitability,
  } = fields;

  const terrainTemplate = buildTerrainTemplate(seed);
  const systems = buildMountainSystems(terrainTemplate, seed);
  const allRidges = systems.flatMap((system, id) => buildScratchRidges(system, seed, id));

  for (let y = 0; y < MAP_H; y++) {
    for (let x = 0; x < MAP_W; x++) {
      const i = indexOf(x, y);
      const { px, py } = normalizeCoord(x, y);
      const terrainLarge = (fbm(x * 0.65, y * 0.65, seed + 1) - 0.5) * 0.23;
      const terrainRegional = (valueNoise(x * 0.8, y * 0.8, seed + 2, 42) - 0.5) * 0.16;
      const { pressure: coastPressure } = computeCoastLower(px, py, terrainTemplate, seed);
      let e = 0.42 + terrainLarge + terrainRegional - coastPressure * (0.22 + terrainTemplate.deposition * 0.040);
      let mountainMaskMax = 0;
      for (let s = 0; s < systems.length; s++) {
        const system = systems[s];
        const mask = ellipticalMask(px, py, system);
        mountainMaskMax = Math.max(mountainMaskMax, mask);
        const broad = Math.pow(mask, lerp(2.0, 1.35, system.massifness)) * system.height * lerp(0.13, 0.25, system.massifness);
        e += broad;
        arcSpineField[i] = Math.max(arcSpineField[i], mask * (system.role === "minor" ? 0.42 : system.role === "primary" ? 0.86 : 0.72));
      }
      for (const ridge of allRidges) {
        const r = ridgeContribution(px, py, ridge, seed);
        if (r <= 0) continue;
        e += r;
        branchRidgeField[i] = clamp(branchRidgeField[i] + r * 5.0);
        arcSpineField[i] = clamp(Math.max(arcSpineField[i], r * 4.6));
      }
      const macro = (fbm(x * terrainTemplate.macroNoiseScale * 48, y * terrainTemplate.macroNoiseScale * 48, seed + 500) - 0.5) * 2;
      const global = (valueNoise(x * 0.23, y * 0.23, seed + 501, 38) - 0.5) * 2;
      const scratch = (fbm(x * 2.2, y * 2.2, seed + 502) - 0.5) * 2;
      e += macro * terrainTemplate.macroNoiseStrength * (0.38 + mountainMaskMax * 0.80);
      e += global * 0.020;
      e += scratch * terrainTemplate.scratchNoiseStrength * mountainMaskMax;
      elevation[i] = clamp(e, 0.025, 1.08);
      visibleRavineField[i] = clamp(Math.abs(scratch) * mountainMaskMax * 0.30 + Math.max(0, -scratch) * mountainMaskMax * 0.40);
      surfaceTextureField[i] = clamp(Math.abs(macro) * 0.12 + Math.abs(scratch) * mountainMaskMax * 0.46);
      valleyField[i] = clamp(Math.max(0, -scratch) * mountainMaskMax * 0.18);
      moisture[i] = clamp(0.45 + coastPressure * 0.28 - elevation[i] * 0.20 + (fbm(x * 1.1, y * 1.1, seed + 503) - 0.5) * 0.16);
    }
  }

  let seaLevel = quantile(elevation, terrainTemplate.seaRatio);
  seaLevel = clamp(seaLevel, 0.20, 0.47);
  classifyWater(elevation, seaLevel, sea, ocean, lake);
  recomputeSlope(elevation, sea, slope);

  const filled = new Float32Array(SIZE);
  priorityFloodFlow(elevation, sea, flowTo, filled);
  computeFlowAccumulation(sea, flowTo, filled, flowAccum);
  const { riverPaths, mainRivers, tributaryRivers, smallStreams } = buildRiverNetwork(seed, terrainTemplate, sea, lake, elevation, slope, flowTo, flowAccum, river, erosionField);
  enforceLandGradient(elevation, sea, seaLevel);
  deriveFields(seed, terrainTemplate, fields, seaLevel);

  const prefectureMask = makePrefectureMask(seed, sea, elevation, slope, river);
  const regional = generateRegionalPrefectures(seed, sea, elevation, slope, river, ridgeField, flowAccum, prefectureMask);
  const prefectureRegionId = regional.regionId;
  const regionalDebug = regional.debug;
  const regionalPrefectureBorders = extractRegionBorderSegments(prefectureRegionId, sea);
  const prefectureBorder = extractMaskBorder(prefectureMask, sea);

  let landCount = 0;
  let waterCount = 0;
  let mountainCount = 0;
  let plainCount = 0;
  let primarySpineStrength = 0;
  let spineSamples = 0;
  for (let i = 0; i < SIZE; i++) {
    if (sea[i]) { waterCount++; continue; }
    landCount++;
    if (elevation[i] > 0.60 || ridgeField[i] > 0.58) mountainCount++;
    if (plain[i] > 0.36) plainCount++;
    if (arcSpineField[i] > 0.55) { primarySpineStrength += arcSpineField[i]; spineSamples++; }
  }
  primarySpineStrength /= Math.max(1, spineSamples);
  const terrainDebug = {
    primarySpineStrength,
    riverConnectivityRate: mainRivers.length ? mainRivers.filter((path) => path.some(([x, y], k) => k > path.length * 0.45 && sea[indexOf(x, y)])).length / mainRivers.length : 0,
    smallIslandCount: 0,
    largeInlandLakeCount: lake.reduce((a, v) => a + v, 0) > 120 ? 1 : 0,
    depositionLowlandArea: depositionalLowland.reduce((a, v, i) => a + (!sea[i] && v > 0.24 ? 1 : 0), 0),
    smallStreamCount: smallStreams.length,
    erosionGullyCount: 0,
    branchRavineCount: 0,
    simpleTerrainSystem: true,
    seaRatio: waterCount / SIZE,
    landCount,
    mountainRatio: mountainCount / Math.max(1, landCount),
    plainRatio: plainCount / Math.max(1, landCount),
    mountainSystemCount: systems.length,
  };

  return {
    terrainTemplate,
    seaLevel,
    elevation,
    moisture,
    slope,
    sea,
    ocean,
    lake,
    river,
    floodplain,
    plain,
    agriculture,
    ridgeField,
    valleyField,
    visibleRavineField,
    surfaceTextureField,
    basinField,
    coastalLowland,
    flowAccum,
    erosionField,
    depositionField,
    arcSpineField,
    branchRidgeField,
    depositionalLowland,
    alluvialFanField,
    deltaField,
    naturalBarrierScore,
    portSuitability,
    crossingSuitability,
    passSuitability,
    prefectureMask,
    prefectureBorder,
    prefectureRegionId,
    regionalDebug,
    terrainDebug,
    regionalPrefectureBorders,
    riverPaths,
    mainRivers,
    tributaryRivers,
    smallStreams,
  };
}