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

export function generateTerrainAndRivers(seed) {
  let prefectureMask;
  let prefectureBorder;

  const {
    elevation,
    moisture,
    slope,
    sea,
    river,
    floodplain,
    plain,
    agriculture,
    ridgeField,
    valleyField,
    basinField,
    coastalLowland,
    flowAccum,
    erosionField,
    depositionField,
    flowTo,
    portSuitability,
    crossingSuitability,
    passSuitability,
  } = createMapFields();

  const coastAngle = rand(seed, 11) * Math.PI * 2;
  const coastX = Math.cos(coastAngle);
  const coastY = Math.sin(coastAngle);
  const coastThreshold = 0.22 + rand(seed, 12) * 0.22;
  const coastStrength = 0.15 + rand(seed, 13) * 0.23;

  const seaLevel = 0.285;

  const mountainBlobs = Array.from({ length: 2 + Math.floor(rand(seed, 98) * 3) }, (_, i) => ({
    x: rand(seed, 100 + i) * MAP_W,
    y: rand(seed, 200 + i) * MAP_H,
    r: 10 + rand(seed, 300 + i) * 24,
    h: 0.08 + rand(seed, 400 + i) * 0.16,
  }));

  const ridgeBands = Array.from({ length: 5 + Math.floor(rand(seed, 97) * 4) }, (_, i) => ({
    x: rand(seed, 1500 + i) * MAP_W,
    y: rand(seed, 1600 + i) * MAP_H,
    angle: rand(seed, 1700 + i) * Math.PI * 2,
    width: 3 + rand(seed, 1800 + i) * 7,
    length: 42 + rand(seed, 1900 + i) * 92,
    h: 0.11 + rand(seed, 2000 + i) * 0.22,
  }));

  for (let y = 0; y < MAP_H; y++) {
    for (let x = 0; x < MAP_W; x++) {
      const nx = x / (MAP_W - 1) - 0.5;
      const ny = y / (MAP_H - 1) - 0.5;
      const i = indexOf(x, y);

      const warpX = (fbm(x * 0.62 + 180, y * 0.62 - 90, seed + 3101) - 0.5) * 13;
      const warpY = (fbm(x * 0.62 - 70, y * 0.62 + 210, seed + 3201) - 0.5) * 13;
      const wx = x + warpX;
      const wy = y + warpY;

      let mountains = 0;
      for (const blob of mountainBlobs) {
        const d = Math.hypot(wx - blob.x, wy - blob.y) / blob.r;
        mountains += Math.exp(-d * d * 2.35) * blob.h;
      }

      let ridges = 0;
      for (const ridge of ridgeBands) {
        const dx = wx - ridge.x;
        const dy = wy - ridge.y;
        const along = dx * Math.cos(ridge.angle) + dy * Math.sin(ridge.angle);
        const perp = -dx * Math.sin(ridge.angle) + dy * Math.cos(ridge.angle);
        const lengthFade = smoothstep(1 - Math.abs(along) / ridge.length);
        const serration = 0.72 + valueNoise(wx + along * 0.15, wy + perp * 0.15, seed + 2220, 8) * 0.56;
        ridges += Math.exp(-(perp * perp) / (ridge.width * ridge.width)) * lengthFade * ridge.h * serration;
      }

      const directionalCoast = nx * coastX + ny * coastY;
      const coastWave = (fbm(wx * 0.72, wy * 0.72, seed + 2222) - 0.5) * 0.12 + (valueNoise(wx, wy, seed + 2233, 18) - 0.5) * 0.08;
      const coastLower = smoothstep((directionalCoast + coastWave - coastThreshold) / 0.26);
      // Four terrain-noise bands from continental structure to fine surface roughness.
      const terrainLarge = fbm(wx * 0.36 + 40, wy * 0.36 - 60, seed + 710);
      const terrainRegional = fbm(wx * 0.95 + 80, wy * 0.95 - 20, seed + 777);
      const terrainLocal = fbm(wx * 2.05 + 17, wy * 2.05 - 31, seed + 1777);
      const terrainFine = valueNoise(wx * 2.9 + 11, wy * 2.9 - 19, seed + 2444, 4.5);
      const fineDissection = Math.abs(terrainLocal - 0.5) * 0.08 + Math.abs(terrainFine - 0.5) * 0.035;
      const basin = 0.1 * Math.sin((nx * 3.1 + ny * 1.7 + rand(seed, 15)) * Math.PI) - 0.045 * Math.cos((nx * 5.2 - ny * 3.6 + rand(seed, 16)) * Math.PI);
      const rawElevation =
        0.30 * terrainLarge +
        0.235 * terrainRegional +
        0.105 * terrainLocal +
        0.055 * terrainFine +
        mountains * 0.54 +
        ridges * 1.22 +
        basin +
        fineDissection -
        coastLower * (coastStrength + 0.19) +
        0.055;

      elevation[i] = clamp(0.5 + (rawElevation - 0.5) * 1.26);
      ridgeField[i] = clamp(ridges * 4.8 + Math.max(0, mountains - 0.10) * 0.95 + fineDissection * 2.0);
      basinField[i] = clamp(Math.max(0, -basin) * 3.0 + (1 - coastLower) * Math.max(0, 0.42 - elevation[i]) * 0.7);
      moisture[i] = clamp(0.44 * fbm(wx + 400, wy - 200, seed + 333) + 0.18 * valueNoise(wx, wy, seed + 343, 11) + 0.22 * (1 - Math.abs(ny * 1.7)) + 0.28 * coastLower - Math.max(0, elevation[i] - 0.62) * 0.22);
    }
  }

  for (let y = 0; y < MAP_H; y++) {
    for (let x = 0; x < MAP_W; x++) {
      const i = indexOf(x, y);
      const nx = x / (MAP_W - 1) - 0.5;
      const ny = y / (MAP_H - 1) - 0.5;
      const directionalCoast = nx * coastX + ny * coastY;
      const coastNoise = (fbm(x * 0.95, y * 0.95, seed + 2222) - 0.5) * 0.14 + (valueNoise(x, y, seed + 2233, 13) - 0.5) * 0.08;
      const oceanSide = directionalCoast + coastNoise > coastThreshold + 0.055;
      if (elevation[i] < seaLevel || oceanSide) sea[i] = 1;
      if (sea[i]) elevation[i] = Math.min(elevation[i], seaLevel - 0.018 + hash2(x, y, seed + 2311) * 0.012);
    }
  }

  // Align coastal elevation with the sea mask. This prevents artificial one-cell cliffs
  // when the directional coastline cuts through a high terrain cell.
  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 nearestSea = INF;
      for (let dy = -7; dy <= 7; dy++) {
        for (let dx = -7; dx <= 7; dx++) {
          const nx = x + dx;
          const ny = y + dy;
          if (!inside(nx, ny) || !sea[indexOf(nx, ny)]) continue;
          nearestSea = Math.min(nearestSea, Math.hypot(dx, dy));
        }
      }
      if (nearestSea <= 7) {
        const coastalCap = seaLevel + 0.018 + nearestSea * 0.028 + Math.max(0, fbm(x * 1.4, y * 1.4, seed + 2350) - 0.5) * 0.022;
        elevation[i] = Math.min(elevation[i], coastalCap);
        coastalLowland[i] = clamp(1 - nearestSea / 7);
      }
    }
  }

  for (let y = 1; y < MAP_H - 1; y++) {
    for (let x = 1; x < MAP_W - 1; x++) {
      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[indexOf(x, y)] = clamp(Math.sqrt(gx * gx + gy * gy) * 10.5);
    }
  }

  const landOrder = [];
  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 low = i;
      let best = elevation[i] + 0.012 * hash2(x, y, seed + 2468);
      let localMean = 0;
      let localMax = elevation[i];
      let localMin = elevation[i];
      let nCount = 0;
      for (const [nx, ny] of neighbors8(x, y)) {
        const ni = indexOf(nx, ny);
        const ev = elevation[ni];
        localMean += ev;
        localMax = Math.max(localMax, ev);
        localMin = Math.min(localMin, ev);
        nCount++;
        const directed = ev + 0.008 * hash2(nx, ny, seed + 2469);
        if (directed < best || sea[ni]) {
          best = directed;
          low = ni;
        }
      }
      if (low !== i) flowTo[i] = low;
      localMean /= Math.max(1, nCount);
      const hollow = Math.max(0, localMean - elevation[i]);
      const relief = localMax - localMin;
      valleyField[i] = clamp(hollow * 8.4 + Math.max(0, 0.42 - elevation[i]) * 0.32 + moisture[i] * 0.08 - ridgeField[i] * 0.18);
      basinField[i] = clamp(basinField[i] + hollow * 2.4 + (relief < 0.055 && elevation[i] < 0.55 ? 0.18 : 0));
      flowAccum[i] = 0.7 + moisture[i] * 0.7 + valleyField[i] * 0.55;
      landOrder.push(i);
    }
  }
  landOrder.sort((a, b) => elevation[b] - elevation[a]);
  for (const i of landOrder) {
    const to = flowTo[i];
    if (to >= 0 && to !== i) flowAccum[to] += flowAccum[i] * 0.82;
  }
  let maxFlowAccum = 0;
  for (let i = 0; i < SIZE; i++) if (!sea[i]) maxFlowAccum = Math.max(maxFlowAccum, flowAccum[i]);
  if (maxFlowAccum > 0) {
    for (let i = 0; i < SIZE; i++) flowAccum[i] = clamp(flowAccum[i] / maxFlowAccum);
  }
  for (let i = 0; i < SIZE; i++) {
    if (!sea[i]) valleyField[i] = clamp(valleyField[i] * 0.68 + Math.pow(flowAccum[i], 0.55) * 0.48);
  }

  // First-order fluvial shaping: cut valley floors on steep/high-flow cells and
  // deposit gently in coastal lowlands and basin floors. This gives visible
  // river valleys without destroying the macro terrain structure.
  const shapedElevation = 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;
      const flow = Math.pow(flowAccum[i], 0.46);
      const incisionNoise = 0.82 + hash2(x, y, seed + 8120) * 0.36;
      const steepValley = clamp(flow * (0.058 + slope[i] * 0.21 + ridgeField[i] * 0.046) * incisionNoise);
      const lateralCut = clamp(Math.pow(flowAccum[i], 0.66) * valleyField[i] * 0.078);
      const lowSettling = clamp(flow * (coastalLowland[i] * 0.036 + basinField[i] * 0.020 + (elevation[i] < 0.40 ? 0.012 : 0)) * (1 - slope[i] * 0.82));
      erosionField[i] = steepValley + lateralCut;
      depositionField[i] = lowSettling;
      shapedElevation[i] = clamp(elevation[i] - steepValley - lateralCut + lowSettling * 0.72, seaLevel + 0.006, 1);
    }
  }
  elevation.set(shapedElevation);

  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.sqrt(gx * gx + gy * gy) * 10.5);
      valleyField[i] = clamp(valleyField[i] + erosionField[i] * 2.1 + depositionField[i] * 0.8 - ridgeField[i] * 0.06);
      basinField[i] = clamp(basinField[i] + depositionField[i] * 1.6);
    }
  }

  const sourceCandidates = [];
  for (let y = 4; y < MAP_H - 4; y++) {
    for (let x = 4; x < MAP_W - 4; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      const score = elevation[i] * 0.38 + moisture[i] * 0.24 + ridgeField[i] * 0.08 + flowAccum[i] * 0.56 + valleyField[i] * 0.28 + hash2(x, y, seed + 9000) * 0.06;
      if (elevation[i] > 0.40 && elevation[i] < 0.82 && moisture[i] > 0.28 && flowAccum[i] > 0.020 && ridgeField[i] < 0.88) sourceCandidates.push({ x, y, score });
    }
  }

  const sources = pickEntities(sourceCandidates, {
    max: 20 + Math.floor(rand(seed, 910) * 28),
    minDistance: 8,
    threshold: 0.53 + rand(seed, 911) * 0.11,
    seed,
  });

  function nearestWaterGoal(from) {
    let bestSea = null;
    let bestScore = INF;
    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;
        const d = Math.hypot(x - from.x, y - from.y);
        const score = d - coastalLowland[indexOf(Math.max(0, Math.min(MAP_W - 1, from.x)), Math.max(0, Math.min(MAP_H - 1, from.y)))] * 2;
        if (score < bestScore) {
          bestScore = score;
          bestSea = { x, y };
        }
      }
    }
    return bestSea;
  }

  function riverRouteCost(x, y, cx, cy) {
    const i = indexOf(x, y);
    const ci = indexOf(cx, cy);
    if (sea[i]) return 0.18;
    const uphill = Math.max(0, elevation[i] - elevation[ci]);
    const downhill = Math.max(0, elevation[ci] - elevation[i]);
    if (!sea[i] && uphill > 0.035 && flowAccum[i] < flowAccum[ci] + 0.015) return INF;
    return Math.max(
      0.18,
      1 +
        uphill * 86 +
        slope[i] * 0.38 +
        elevation[i] * 0.42 -
        downhill * 2.1 -
        valleyField[i] * 0.92 -
        flowAccum[i] * 0.72 -
        moisture[i] * 0.18 -
        coastalLowland[i] * 0.22
    );
  }

  function forceRiverToWater(path) {
    if (!path.length) return path;
    const [ex, ey] = path[path.length - 1];
    if (sea[indexOf(ex, ey)]) return path;
    const goal = nearestWaterGoal({ x: ex, y: ey });
    if (!goal) return path;
    const startElevation = elevation[indexOf(ex, ey)];
    const tail = aStar({ x: ex, y: ey }, goal, (x, y, cx, cy) => {
      const i = indexOf(x, y);
      const ci = indexOf(cx, cy);
      if (!sea[i] && elevation[i] > Math.max(startElevation + 0.045, elevation[ci] + 0.030)) return INF;
      return riverRouteCost(x, y, cx, cy);
    });
    if (tail.length <= 2) return path;
    return path.concat(tail.slice(1));
  }

  function confluenceAnglePenalty(nx, ny, dx, dy, lengthSoFar) {
    if (lengthSoFar < 7 || river[indexOf(nx, ny)] < 0.24) return 0;
    let best = 0.16;
    const inLen = Math.hypot(dx, dy) || 1;
    for (const [rx, ry] of neighbors8(nx, ny)) {
      if (river[indexOf(rx, ry)] < 0.22) continue;
      const rdx = rx - nx;
      const rdy = ry - ny;
      const cos = clamp((dx * rdx + dy * rdy) / Math.max(0.001, inLen * Math.hypot(rdx, rdy)), -1, 1);
      const angle = Math.acos(cos);
      const shallow = angle < 0.45 ? 0.28 : 0;
      best = Math.min(best, Math.abs(angle - Math.PI * 0.62) * 0.045 + shallow);
    }
    return best;
  }

  function traceRiverPath(startX, startY, bonusSeed = 0) {
    let x = startX;
    let y = startY;
    let lastDx = 0;
    let lastDy = 0;
    const path = [];
    const seen = new Set();
    let accum = 0;

    for (let step = 0; step < 600; step++) {
      const i = indexOf(x, y);
      if (seen.has(i)) break;
      seen.add(i);
      path.push([x, y]);
      river[i] += 0.44 + path.length / 160 + flowAccum[i] * 0.55;
      accum += river[i] + flowAccum[i];
      if (sea[i]) break;

      let best = null;
      let bestValue = INF;
      const currentElevation = elevation[i];
      const preferred = flowTo[i];

      for (const [nx, ny] of neighbors8(x, y)) {
        const ni = indexOf(nx, ny);
        const dx = nx - x;
        const dy = ny - y;
        const drop = currentElevation - elevation[ni];
        const uphill = Math.max(0, -drop);
        if (!sea[ni] && uphill > 0.032 && flowAccum[ni] < flowAccum[i] + 0.018) continue;
        let surrounding = 0;
        let surroundingCount = 0;
        for (const [vx, vy] of neighbors8(nx, ny)) {
          surrounding += elevation[indexOf(vx, vy)];
          surroundingCount++;
        }
        const valley = Math.max(0, surrounding / Math.max(1, surroundingCount) - elevation[ni]);
        const sameDirection = lastDx || lastDy ? (dx * lastDx + dy * lastDy) / Math.max(0.001, Math.hypot(dx, dy) * Math.hypot(lastDx, lastDy)) : 0;
        const straightPenalty = Math.max(0, sameDirection) * 0.075;
        const turnPenalty = sameDirection < -0.35 ? 0.24 : 0;
        const sideSwing = Math.abs(dx * lastDy - dy * lastDx);
        const meanderPhase = Math.sin((path.length + bonusSeed * 0.013) * 0.73) * 0.5 + 0.5;
        const meander = sideSwing * (0.032 + meanderPhase * 0.026);
        const flowBonus = ni === preferred ? 0.62 : 0;
        const junctionPenalty = confluenceAnglePenalty(nx, ny, dx, dy, path.length);
        const noise = (hash2(nx, ny, seed + bonusSeed + step * 11) - 0.5) * 0.04;
        const value =
          elevation[ni] * 1.45 +
          uphill * 88 -
          Math.max(0, drop) * 2.05 -
          valley * 1.05 -
          valleyField[ni] * 1.72 -
          flowAccum[ni] * 0.94 -
          moisture[ni] * 0.14 -
          coastalLowland[ni] * 0.28 -
          (river[ni] > 0 ? 0.22 : 0) -
          flowBonus +
          slope[ni] * 0.04 +
          straightPenalty +
          turnPenalty +
          junctionPenalty * 1.35 -
          meander +
          noise -
          (sea[ni] ? 0.6 : 0);

        if (value < bestValue) {
          bestValue = value;
          best = [nx, ny, dx, dy];
        }
      }
      if (!best) break;
      x = best[0];
      y = best[1];
      lastDx = best[2];
      lastDy = best[3];
    }

    const forced = forceRiverToWater(path);
    if (forced.length > path.length) {
      for (const [rx, ry] of forced.slice(path.length)) {
        const ri = indexOf(rx, ry);
        river[ri] += 0.32 + flowAccum[ri] * 0.4;
        accum += river[ri] + flowAccum[ri];
      }
    }
    return { path: forced, accum };
  }

  function traceSmallStreamPath(startX, startY, bonusSeed = 0) {
    let x = startX;
    let y = startY;
    let lastDx = 0;
    let lastDy = 0;
    const path = [];
    const seen = new Set();
    for (let step = 0; step < 160; step++) {
      const i = indexOf(x, y);
      if (seen.has(i)) break;
      seen.add(i);
      path.push([x, y]);
      river[i] += 0.12 + flowAccum[i] * 0.18;
      if ((river[i] > 0.48 && path.length > 5) || sea[i]) break;
      let best = null;
      let bestValue = INF;
      for (const [nx, ny] of neighbors8(x, y)) {
        const ni = indexOf(nx, ny);
        const dx = nx - x;
        const dy = ny - y;
        const drop = elevation[i] - elevation[ni];
        const sameDirection = lastDx || lastDy ? (dx * lastDx + dy * lastDy) / Math.max(0.001, Math.hypot(dx, dy) * Math.hypot(lastDx, lastDy)) : 0;
        const value = elevation[ni] * 1.2 + Math.max(0, -drop) * 26 - Math.max(0, drop) * 1.4 - valleyField[ni] * 1.15 - flowAccum[ni] * 0.55 - moisture[ni] * 0.12 + Math.max(0, sameDirection) * 0.04 - Math.abs(dx * lastDy - dy * lastDx) * 0.018 + (hash2(nx, ny, seed + bonusSeed + step * 13) - 0.5) * 0.05;
        if (value < bestValue) { bestValue = value; best = [nx, ny, dx, dy]; }
      }
      if (!best) break;
      x = best[0];
      y = best[1];
      lastDx = best[2];
      lastDy = best[3];
    }
    return path;
  }

  const riverPaths = [];
  const riverScores = [];
  for (const source of sources) {
    const { path, accum } = traceRiverPath(source.x, source.y, 0);
    if (path.length > 6) {
      riverPaths.push(path);
      riverScores.push(path.length + accum * 0.18);
    }
  }

  const preliminaryMainRiverCells = new Set(riverPaths.slice().sort((a, b) => b.length - a.length).slice(0, 5).flatMap((path) => path.map(([x, y]) => `${x},${y}`)));
  const tributarySources = pickEntities(sourceCandidates
    .filter((p) => !preliminaryMainRiverCells.has(`${p.x},${p.y}`))
    .map((p) => ({ ...p, score: p.score + flowAccum[indexOf(p.x, p.y)] * 0.75 + valleyField[indexOf(p.x, p.y)] * 0.24 })), {
      max: 14 + Math.floor(rand(seed, 915) * 20),
      minDistance: 6,
      threshold: 0.45,
      seed: seed + 916,
      jitter: 0.02,
    });
  for (const source of tributarySources) {
    const { path, accum } = traceRiverPath(source.x, source.y, 4000 + source.x * 7 + source.y * 11);
    if (path.length > 8) {
      riverPaths.push(path);
      riverScores.push(path.length * 0.7 + accum * 0.12);
    }
  }

  const streamPaths = [];
  const streamSources = pickEntities(sourceCandidates
    .map((p) => ({ ...p, score: valleyField[indexOf(p.x, p.y)] * 0.46 + flowAccum[indexOf(p.x, p.y)] * 0.36 + moisture[indexOf(p.x, p.y)] * 0.18 + hash2(p.x, p.y, seed + 918) * 0.05 }))
    .filter((p) => p.score > 0.18), {
      max: 22 + Math.floor(rand(seed, 919) * 20),
      minDistance: 4,
      threshold: 0.18,
      seed: seed + 919,
      jitter: 0.015,
    });
  for (const source of streamSources) {
    const path = traceSmallStreamPath(source.x, source.y, 7000 + source.x * 5 + source.y * 17);
    if (path.length > 4) streamPaths.push(path);
  }

  if (riverPaths.length === 0 && sourceCandidates.length > 0) {
    const fallback = sourceCandidates.slice().sort((a, b) => b.score - a.score)[0];
    let bestSea = null;
    let bestSeaDist = INF;
    for (let y = 0; y < MAP_H; y++) {
      for (let x = 0; x < MAP_W; x++) {
        if (!sea[indexOf(x, y)]) continue;
        const d = Math.hypot(x - fallback.x, y - fallback.y);
        if (d < bestSeaDist) {
          bestSeaDist = d;
          bestSea = { x, y };
        }
      }
    }
    if (bestSea) {
      const fallbackPath = aStar(fallback, bestSea, (x, y, cx, cy) => {
        const i = indexOf(x, y);
        const ci = indexOf(cx, cy);
        if (sea[i]) return 0.25;
        const uphill = Math.max(0, elevation[i] - elevation[ci]) * 24;
        const downhill = Math.max(0, elevation[ci] - elevation[i]) * 1.8;
        return Math.max(0.24, 1 + uphill + slope[i] * 0.7 + elevation[i] * 0.8 - downhill - Math.min(0.55, river[i] * 0.1));
      });
      if (fallbackPath.length > 6) {
        let accum = 0;
        for (const [x, y] of fallbackPath) {
          const i = indexOf(x, y);
          river[i] += 0.42;
          accum += river[i];
        }
        riverPaths.push(fallbackPath);
        riverScores.push(fallbackPath.length + accum * 0.18);
      }
    }
  }

  function sanitizeDownhillRiverPath(path, tolerance = 0.040) {
    if (!path || path.length < 2) return path || [];
    const out = [path[0]];
    for (let k = 1; k < path.length; k++) {
      const [px, py] = out[out.length - 1];
      const [x, y] = path[k];
      const pi = indexOf(px, py);
      const i = indexOf(x, y);
      if (!sea[i] && elevation[i] > elevation[pi] + tolerance) break;
      out.push(path[k]);
      if (sea[i]) break;
    }
    return out.length >= 2 ? out : [];
  }
  function trimMountainHeadwaters(path) {
    if (!path || path.length < 4) return path || [];
    let start = 0;
    while (start < path.length - 3) {
      const [x, y] = path[start];
      const i = indexOf(x, y);
      if (sea[i]) break;
      if (elevation[i] <= 0.72 && (valleyField[i] >= 0.18 || flowAccum[i] >= 0.05)) break;
      start++;
    }
    return path.slice(start);
  }
  for (let r = 0; r < riverPaths.length; r++) riverPaths[r] = sanitizeDownhillRiverPath(trimMountainHeadwaters(riverPaths[r]), 0.032);
  for (let r = riverPaths.length - 1; r >= 0; r--) if (riverPaths[r].length < 2) riverPaths.splice(r, 1);
  for (let r = 0; r < streamPaths.length; r++) streamPaths[r] = sanitizeDownhillRiverPath(trimMountainHeadwaters(streamPaths[r]), 0.026);
  for (let r = streamPaths.length - 1; r >= 0; r--) if (streamPaths[r].length < 2) streamPaths.splice(r, 1);
  river.fill(0);
  for (const path of riverPaths) {
    for (let k = 0; k < path.length; k++) {
      const [x, y] = path[k];
      const i = indexOf(x, y);
      river[i] += 0.42 + k / 170 + flowAccum[i] * 0.55;
    }
  }
  for (const path of streamPaths) {
    for (let k = 0; k < path.length; k++) {
      const [x, y] = path[k];
      const i = indexOf(x, y);
      river[i] += 0.11 + flowAccum[i] * 0.18;
    }
  }

  const expandedRiver = new Float32Array(river);
  for (let y = 1; y < MAP_H - 1; y++) {
    for (let x = 1; x < MAP_W - 1; x++) {
      const i = indexOf(x, y);
      if (river[i] <= 0) continue;
      for (const [nx, ny] of neighbors8(x, y)) {
        expandedRiver[indexOf(nx, ny)] = Math.max(expandedRiver[indexOf(nx, ny)], river[i] * 0.35);
      }
    }
  }
  river.set(expandedRiver);

  // Second fluvial pass uses the actual traced river network. Main channels cut
  // visible V-shaped valleys; lower reaches accumulate alluvial deposits.
  const fluvialElevation = 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] || river[i] <= 0.02) continue;
      const r = clamp(river[i] / 3.4);
      const channelCut = clamp(Math.pow(r, 0.55) * (0.060 + slope[i] * 0.145 + ridgeField[i] * 0.038));
      const valleyWiden = clamp(Math.pow(r, 0.72) * (0.020 + Math.max(0, elevation[i] - seaLevel) * 0.058 + valleyField[i] * 0.040));
      const alluvium = clamp(Math.pow(r, 0.72) * (coastalLowland[i] * 0.030 + basinField[i] * 0.020 + (slope[i] < 0.10 ? 0.010 : 0)));
      erosionField[i] = clamp(erosionField[i] + channelCut + valleyWiden);
      depositionField[i] = clamp(depositionField[i] + alluvium);
      fluvialElevation[i] = clamp(elevation[i] - channelCut - valleyWiden + alluvium, seaLevel + 0.005, 1);
      valleyField[i] = clamp(valleyField[i] + r * 0.62 + channelCut * 6.4);
      basinField[i] = clamp(basinField[i] + alluvium * 3.2);
    }
  }
  // Lateral valley carving around the traced river network deepens valleys and
  // makes ridge/valley contrast legible at the map scale.
  for (const path of riverPaths) {
    for (const [rx, ry] of path) {
      const ri = indexOf(rx, ry);
      const r = clamp(river[ri] / 3.0);
      const radius = r > 0.48 ? 2 : 1;
      for (let dy = -radius; dy <= radius; dy++) {
        for (let dx = -radius; dx <= radius; dx++) {
          const nx = rx + dx;
          const ny = ry + dy;
          if (!inside(nx, ny)) continue;
          const ni = indexOf(nx, ny);
          if (sea[ni]) continue;
          const d = Math.hypot(dx, dy);
          if (d > radius || d === 0) continue;
          const weight = (radius + 0.35 - d) / (radius + 0.35);
          const carve = Math.max(0, weight) * (0.008 + r * 0.026) * Math.max(0.45, slope[ni] + 0.22);
          fluvialElevation[ni] = clamp(fluvialElevation[ni] - carve, seaLevel + 0.005, 1);
          erosionField[ni] = clamp(erosionField[ni] + carve * 3.0);
          valleyField[ni] = clamp(valleyField[ni] + carve * 12.0);
        }
      }
    }
  }

  // Restore rugged summit relief after strong river incision. This prevents highlands
  // from becoming unnaturally flat or visually concave while keeping valleys cut.
  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 high = clamp((fluvialElevation[i] - 0.62) / 0.26);
      const summit = high * clamp(ridgeField[i] * 1.4 - flowAccum[i] * 0.8);
      const rugged = (valueNoise(x * 2.1 + 19, y * 2.1 - 23, seed + 9661, 3.2) - 0.5) * 0.035;
      const uplift = summit * (0.018 + Math.max(0, rugged));
      if (uplift > 0) {
        fluvialElevation[i] = clamp(fluvialElevation[i] + uplift, seaLevel + 0.005, 1);
        erosionField[i] = Math.max(0, erosionField[i] - uplift * 0.6);
      }
    }
  }

  elevation.set(fluvialElevation);

  // Broad alluvial/coastal/basin plains. The plain score alone is not enough;
  // the elevation surface must also be locally calm, otherwise every lowland
  // still reads as rugged terrain. Smooth only low, wet depositional cells and
  // leave ridges/headwaters untouched.
  for (let pass = 0; pass < 4; pass++) {
    const nextElevation = new Float32Array(elevation);
    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 lowland = clamp(
          coastalLowland[i] * 0.72 +
          basinField[i] * 0.54 +
          valleyField[i] * 0.34 +
          Math.pow(flowAccum[i], 0.58) * 0.24 -
          ridgeField[i] * 0.62 -
          Math.max(0, elevation[i] - 0.54) * 1.65 -
          slope[i] * 0.74
        );
        if (lowland <= 0.12) continue;
        let sum = 0;
        let weight = 0;
        for (let dy = -2; dy <= 2; dy++) {
          for (let dx = -2; dx <= 2; dx++) {
            const nx = x + dx;
            const ny = y + dy;
            const ni = indexOf(nx, ny);
            if (sea[ni]) continue;
            const d = Math.hypot(dx, dy);
            if (d > 2.25) continue;
            const compatible = clamp(1 - Math.abs(elevation[ni] - elevation[i]) / 0.11);
            const w = compatible / (1 + d);
            sum += elevation[ni] * w;
            weight += w;
          }
        }
        if (weight <= 0) continue;
        const localMean = sum / weight;
        const terrace = Math.round(localMean * 42) / 42;
        const target = lerp(localMean, terrace, 0.28);
        nextElevation[i] = clamp(lerp(elevation[i], target, lowland * 0.42), seaLevel + 0.006, 1);
        if (lowland > 0.55) {
          depositionField[i] = clamp(depositionField[i] + lowland * 0.018);
          erosionField[i] = Math.max(0, erosionField[i] - lowland * 0.012);
        }
      }
    }
    elevation.set(nextElevation);
  }

  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.sqrt(gx * gx + gy * gy) * 11.2);
    }
  }

  // Re-trim visible river paths after fluvial reshaping changes local elevation.
  for (let r = 0; r < riverPaths.length; r++) riverPaths[r] = sanitizeDownhillRiverPath(trimMountainHeadwaters(riverPaths[r]), 0.028);
  for (let r = riverPaths.length - 1; r >= 0; r--) if (riverPaths[r].length < 2) riverPaths.splice(r, 1);
  for (let r = 0; r < streamPaths.length; r++) streamPaths[r] = sanitizeDownhillRiverPath(trimMountainHeadwaters(streamPaths[r]), 0.022);
  for (let r = streamPaths.length - 1; r >= 0; r--) if (streamPaths[r].length < 2) streamPaths.splice(r, 1);

  const mainRivers = riverPaths
    .map((p, i) => ({ path: p, score: riverScores[i] }))
    .sort((a, b) => b.score - a.score)
    .slice(0, Math.min(6, riverPaths.length))
    .map((x) => x.path);

  if (mainRivers.length === 0 && riverPaths.length > 0) mainRivers.push(riverPaths[0]);
  if (mainRivers.length === 0) {
    let start = null;
    let startScore = -INF;
    for (let y = 4; y < MAP_H - 4; y++) {
      for (let x = 4; x < MAP_W - 4; x++) {
        const i = indexOf(x, y);
        if (sea[i]) continue;
        const score = elevation[i] * 0.55 + moisture[i] * 0.35 - slope[i] * 0.15;
        if (score > startScore) {
          startScore = score;
          start = { x, y };
        }
      }
    }
    if (start) {
      let goal = null;
      let goalDist = INF;
      for (let y = 0; y < MAP_H; y++) {
        for (let x = 0; x < MAP_W; x++) {
          if (!sea[indexOf(x, y)]) continue;
          const d = Math.hypot(x - start.x, y - start.y);
          if (d < goalDist) {
            goalDist = d;
            goal = { x, y };
          }
        }
      }
      if (goal) {
        const fallbackPath = aStar(start, goal, (x, y, cx, cy) => {
          const i = indexOf(x, y);
          const ci = indexOf(cx, cy);
          if (sea[i]) return 0.2;
          const uphillBias = Math.max(0, elevation[i] - elevation[ci]) * 22;
          const downhillBias = Math.max(0, elevation[ci] - elevation[i]) * 1.7;
          return Math.max(0.25, 1 + uphillBias + slope[i] * 0.65 + elevation[i] * 0.8 - downhillBias);
        });
        if (fallbackPath.length > 4) {
          riverPaths.push(fallbackPath);
          mainRivers.push(fallbackPath);
          for (const [x, y] of fallbackPath) river[indexOf(x, y)] += 0.4;
        }
      }
    }
  }

  const mainRiverCells = new Set(mainRivers.flatMap((path) => path.map(([x, y]) => `${x},${y}`)));
  const tributaryRivers = riverPaths.filter((path) => path.some(([x, y]) => !mainRiverCells.has(`${x},${y}`)) && !mainRivers.includes(path));
  const smallStreams = streamPaths.filter((path) => path.length >= 5);

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

  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 low = 1 - clamp((elevation[i] - 0.28) / 0.4);
      const flat = 1 - slope[i];
      const valleyPlain = valleyField[i] * 0.44 + basinField[i] * 0.36 + coastalLowland[i] * 0.55;
      plain[i] = clamp(low * 0.44 + flat * 0.58 + valleyPlain - ridgeField[i] * 0.28 - (elevation[i] > 0.62 ? 0.48 : 0));

      let nearRiver = 0;
      for (let dy = -4; dy <= 4; dy++) {
        for (let dx = -4; dx <= 4; dx++) {
          const nx = x + dx;
          const ny = y + dy;
          if (!inside(nx, ny)) continue;
          nearRiver = Math.max(nearRiver, river[indexOf(nx, ny)] / (1 + Math.hypot(dx, dy)));
        }
      }

      const fan = clamp(valleyField[i] * (1 - coastalLowland[i]) * (elevation[i] > 0.34 && elevation[i] < 0.58 ? 0.9 : 0.35) * (1 - slope[i] * 0.55));
      floodplain[i] = clamp(nearRiver * plain[i] * 0.92 + coastalLowland[i] * nearRiver * 0.22);
      agriculture[i] = clamp(plain[i] * 0.58 + fan * 0.26 + basinField[i] * 0.2 + moisture[i] * 0.14 + clamp(nearRiver) * 0.32 - slope[i] * 0.34 - ridgeField[i] * 0.18 - floodplain[i] * 0.06);
    }
  }

  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;
      let seaNear = 0;
      let riverNear = 0;
      let sheltered = 0;

      for (let dy = -5; dy <= 5; dy++) {
        for (let dx = -5; dx <= 5; dx++) {
          const nx = x + dx;
          const ny = y + dy;
          if (!inside(nx, ny)) continue;
          const d = Math.hypot(dx, dy);
          if (sea[indexOf(nx, ny)]) seaNear += 1 / (1 + d);
          riverNear = Math.max(riverNear, river[indexOf(nx, ny)] / (1 + d));
        }
      }

      for (let dy = -2; dy <= 2; dy++) {
        for (let dx = -2; dx <= 2; dx++) {
          const nx = x + dx;
          const ny = y + dy;
          if (inside(nx, ny) && !sea[indexOf(nx, ny)]) sheltered += 1;
        }
      }

      const isDelta = riverNear > 0.22 && coastalLowland[i] > 0.18;
      const bayShelter = sheltered * 0.012 + seaNear * 0.055 + coastalLowland[i] * 0.16;
      portSuitability[i] = clamp(bayShelter + riverNear * 0.24 + (isDelta ? 0.22 : 0) + plain[i] * 0.08 - slope[i] * 0.48 - ridgeField[i] * 0.16);
    }
  }

  for (let y = 3; y < MAP_H - 3; y++) {
    for (let x = 3; x < MAP_W - 3; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      const r = river[i];
      if (r < 0.2 || r > 1.85) continue;
      let bankPlain = 0;
      for (const [nx, ny] of neighbors8(x, y)) bankPlain += plain[indexOf(nx, ny)];
      crossingSuitability[i] = clamp(r * 0.34 + (bankPlain / 8) * 0.54 + valleyField[i] * 0.18 - slope[i] * 0.55 - floodplain[i] * 0.06);
    }
  }

  for (let y = 4; y < MAP_H - 4; y++) {
    for (let x = 4; x < MAP_W - 4; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      const e = elevation[i];
      if (e < 0.43 || e > 0.82) continue;
      const ewHigh = (elevation[indexOf(x - 3, y)] + elevation[indexOf(x + 3, y)]) / 2;
      const nsHigh = (elevation[indexOf(x, y - 3)] + elevation[indexOf(x, y + 3)]) / 2;
      const diagLow = Math.min(
        elevation[indexOf(x - 3, y - 3)],
        elevation[indexOf(x + 3, y + 3)],
        elevation[indexOf(x - 3, y + 3)],
        elevation[indexOf(x + 3, y - 3)]
      );
      passSuitability[i] = clamp((Math.max(ewHigh, nsHigh) - e) * 2.2 + (e - diagLow) * 0.55 + valleyField[i] * 0.28 - ridgeField[i] * 0.18 - slope[i] * 0.2);
    }
  }


  return {
    elevation,
    moisture,
    slope,
    sea,
    river,
    floodplain,
    plain,
    agriculture,
    ridgeField,
    valleyField,
    basinField,
    coastalLowland,
    flowAccum,
    erosionField,
    depositionField,
    portSuitability,
    crossingSuitability,
    passSuitability,
    prefectureMask,
    prefectureBorder,
    prefectureRegionId,
    regionalDebug,
    regionalPrefectureBorders,
    riverPaths,
    mainRivers,
    tributaryRivers,
    smallStreams,
  };
}