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,
  makePrefectureMask,
  neighbors8,
} from "./mapGeneratorHelpers.js";
import { buildNaturalCompartments } from "./adminRegions.js";
import { createRectContext, createRectTerrainFields, rectIndexOf, rectInside, rectNeighbors8, rectQuantile } from "./rectContext.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 softCapElevation(e, start = 0.91, cap = 1.08) {
  if (e <= start) return e;
  const over = e - start;
  const span = Math.max(0.001, cap - start);
  // Hard clipping made high mountains become flat mesas. This keeps peaks high,
  // but compresses only the excess so local relief survives near the top.
  return start + span * (1 - Math.exp(-over / span));
}

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, (v + ridge.phase * 0.37) * 980, 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((u + ridge.phase) * 900, (v - ridge.phase * 0.41) * 1200, 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 along = clamp((u / a + 1) * 0.5);
  const widthWave = 1
    + (system.widthVariance ?? 0.28) * Math.sin((along + (system.phase ?? 0)) * Math.PI * 2.0)
    + (system.widthVariance ?? 0.28) * 0.50 * Math.sin((along * 2.7 + (system.phase ?? 0) * 1.7) * Math.PI * 2.0);
  const endTaper = lerp(0.60, 1.0, Math.sin(along * Math.PI));
  const b = Math.max(0.012, system.width * 0.5 * clamp(widthWave, 0.62, 1.60) * endTaper);
  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 };
}

const TERRAIN_TYPES = [
  {
    id: "tohoku_spine",
    label: "東北型・長大脊梁",
    coastStyle: "parallel_spine",
    mountainMode: "range",
    massifnessRange: [0.06, 0.26],
    seaRatioRange: [0.13, 0.23],
    twoSidedChance: 0.96,
    mountainOffsetRange: [0.47, 0.53],
    baseHeightRange: [0.56, 0.82],
    primaryLengthRange: [0.76, 0.96],
    primaryWidthRange: [0.18, 0.30],
    systemCountRange: [14, 18],
    beltCountRange: [4, 5],
    angleSpread: 0.18,
    crossSpread: 0.58,
    lengthScale: 1.22,
    widthScale: 1.16,
    heightScale: 0.86,
    coastStrength: 0.90,
    plainBiasRange: [0.16, 0.34],
    riverRichnessRange: [0.70, 1.18],
    bigRiverChanceRange: [0.22, 0.46],
  },
  {
    id: "chubu_mountain",
    label: "中部型・交差高山地",
    coastStyle: "outer_coast",
    mountainMode: "massif",
    massifnessRange: [0.42, 0.74],
    seaRatioRange: [0.000, 0.030],
    twoSidedChance: 0.10,
    mountainOffsetRange: [0.16, 0.36],
    baseHeightRange: [0.76, 1.12],
    primaryLengthRange: [0.62, 0.92],
    primaryWidthRange: [0.30, 0.58],
    systemCountRange: [16, 20],
    beltCountRange: [3, 5],
    angleSpread: 0.92,
    crossSpread: 0.82,
    lengthScale: 1.34,
    widthScale: 1.30,
    heightScale: 1.32,
    coastStrength: 0.30,
    plainBiasRange: [0.08, 0.24],
    riverRichnessRange: [0.74, 1.14],
    bigRiverChanceRange: [0.30, 0.60],
  },
  {
    id: "oceanic_archipelago",
    label: "Ocean",
    autoSelectable: false,
    coastStyle: "oceanic_archipelago",
    mountainMode: "mixed",
    massifnessRange: [0.04, 0.28],
    seaRatioRange: [0.955, 0.985],
    twoSidedChance: 1.0,
    mountainOffsetRange: [0.25, 0.55],
    baseHeightRange: [0.20, 0.46],
    primaryLengthRange: [0.20, 0.52],
    primaryWidthRange: [0.08, 0.24],
    systemCountRange: [7, 14],
    beltCountRange: [2, 4],
    angleSpread: 0.70,
    crossSpread: 0.90,
    lengthScale: 0.78,
    widthScale: 0.82,
    heightScale: 0.58,
    coastStrength: 2.10,
    plainBiasRange: [0.12, 0.34],
    riverRichnessRange: [0.18, 0.52],
    bigRiverChanceRange: [0.02, 0.12],
  },
  {
    id: "setouchi_inland_sea",
    label: "瀬戸内型・内海多島",
    coastStyle: "inland_sea",
    mountainMode: "mixed",
    massifnessRange: [0.34, 0.62],
    seaRatioRange: [0.40, 0.56],
    twoSidedChance: 0.96,
    mountainOffsetRange: [0.22, 0.34],
    baseHeightRange: [0.46, 0.78],
    primaryLengthRange: [0.52, 0.78],
    primaryWidthRange: [0.20, 0.38],
    systemCountRange: [20, 28],
    beltCountRange: [4, 5],
    angleSpread: 0.34,
    crossSpread: 0.86,
    lengthScale: 0.92,
    widthScale: 1.02,
    heightScale: 0.78,
    coastStrength: 1.52,
    plainBiasRange: [0.26, 0.50],
    riverRichnessRange: [0.58, 0.96],
    bigRiverChanceRange: [0.18, 0.42],
  },
  {
    id: "kanto_alluvial",
    label: "関東・濃尾型・大河川平野",
    coastStyle: "open_bay",
    mountainMode: "range",
    massifnessRange: [0.10, 0.30],
    seaRatioRange: [0.08, 0.17],
    twoSidedChance: 0.10,
    mountainOffsetRange: [0.28, 0.46],
    baseHeightRange: [0.48, 0.82],
    primaryLengthRange: [0.38, 0.62],
    primaryWidthRange: [0.16, 0.30],
    systemCountRange: [7, 11],
    beltCountRange: [2, 3],
    angleSpread: 0.34,
    crossSpread: 0.62,
    lengthScale: 0.92,
    widthScale: 1.10,
    heightScale: 0.84,
    coastStrength: 0.68,
    plainBiasRange: [0.70, 0.96],
    riverRichnessRange: [1.18, 1.58],
    bigRiverChanceRange: [0.80, 0.98],
  },
  {
    id: "mixed_archipelago",
    label: "混合型・列島変化",
    coastStyle: "mixed_archipelago",
    mountainMode: "mixed",
    massifnessRange: [0.16, 0.72],
    seaRatioRange: [0.13, 0.29],
    twoSidedChance: 0.42,
    mountainOffsetRange: [0.18, 0.40],
    baseHeightRange: [0.68, 1.18],
    primaryLengthRange: [0.46, 0.82],
    primaryWidthRange: [0.20, 0.48],
    systemCountRange: [14, 17],
    beltCountRange: [3, 4],
    angleSpread: 0.50,
    crossSpread: 0.74,
    lengthScale: 1.00,
    widthScale: 1.18,
    heightScale: 1.25,
    coastStrength: 0.96,
    plainBiasRange: [0.22, 0.52],
    riverRichnessRange: [0.72, 1.26],
    bigRiverChanceRange: [0.34, 0.68],
  },
];

function pickTerrainType(seed, requestedType = "auto") {
  if (requestedType && requestedType !== "auto") {
    const normalizedType = requestedType === "touhoku_spine" ? "tohoku_spine" : requestedType;
    const selected = TERRAIN_TYPES.find((type) => type.id === normalizedType);
    if (selected) return selected;
  }
  // Auto excludes explicitly manual/special-purpose templates such as Ocean.
  // The selection remains uniform over ordinary land-bearing templates.
  const autoTypes = TERRAIN_TYPES.filter((type) => type.autoSelectable !== false);
  const index = Math.floor(rand(seed, 10001) * autoTypes.length) % autoTypes.length;
  return autoTypes[index];
}

function rangeValue(seed, salt, [lo, hi]) {
  return lo + rand(seed, salt) * (hi - lo);
}

function rangeInt(seed, salt, [lo, hi]) {
  return Math.round(lo + rand(seed, salt) * (hi - lo));
}

export function buildTerrainTemplate(seed, options = {}) {
  const terrainType = pickTerrainType(seed, options.terrainType || options.generationType || "auto");
  const mountainMode = terrainType.mountainMode === "mixed"
    ? (rand(seed, 12) < 0.42 ? "range" : rand(seed, 13) < 0.72 ? "mixed" : "massif")
    : terrainType.mountainMode;
  let mountainMassifness = rangeValue(seed, 14, terrainType.massifnessRange);
  if (mountainMode === "range") mountainMassifness *= 0.70;
  if (mountainMode === "massif") mountainMassifness = clamp(mountainMassifness + 0.12);

  let coastAngle = rand(seed, 21) * Math.PI * 2;
  let twoSidedCoast = rand(seed, 22) < terrainType.twoSidedChance;
  const seaRatio = rangeValue(seed, 23, terrainType.seaRatioRange);
  let mountainAngle = coastAngle + Math.PI * (rangeValue(seed, 24, terrainType.mountainOffsetRange));
  if (terrainType.id === "tohoku_spine") {
    // 東北型の脊梁は南北/東西だけでなく斜め軸も許容する。
    // 海岸勾配は脊梁軸に概ね直交させるが、山脈自体の向きは独立に選ぶ。
    const axisChoices = [0, Math.PI / 2, Math.PI / 4, -Math.PI / 4, Math.PI * 0.35, Math.PI * 0.65];
    mountainAngle = axisChoices[Math.floor(rand(seed, 2101) * axisChoices.length) % axisChoices.length] + (rand(seed, 2102) - 0.5) * 0.24;
    coastAngle = mountainAngle - Math.PI / 2 + (rand(seed, 2104) - 0.5) * 0.12;
    twoSidedCoast = true;
  }
  const baseHeight = rangeValue(seed, 25, terrainType.baseHeightRange);
  const primaryLength = rangeValue(seed, 26, terrainType.primaryLengthRange);
  const primaryWidth = rangeValue(seed, 28, terrainType.primaryWidthRange);
  const scratchCount = Math.round(lerp(24 + rand(seed, 30) * 20, 16 + rand(seed, 31) * 18, mountainMassifness));
  const mountainSystemCount = rangeInt(seed, 32, terrainType.systemCountRange);
  const mountainBeltCount = rangeInt(seed, 46, terrainType.beltCountRange);

  return {
    seed,
    terrainType: terrainType.id,
    terrainTypeLabel: terrainType.label,
    coastStyle: terrainType.coastStyle,
    seaRatio,
    coastAngle,
    twoSidedCoast,
    coastNoise: 0.040 + rand(seed, 33) * 0.056,
    coastStrength: terrainType.coastStrength,
    mountainMode,
    mountainMassifness,
    mountainAngle,
    mountainAngleSpread: terrainType.angleSpread,
    mountainCrossSpread: terrainType.crossSpread,
    mountainLengthScale: terrainType.lengthScale,
    mountainWidthScale: terrainType.widthScale,
    mountainHeightScale: terrainType.heightScale,
    mountainBeltCount,
    mountainBaseHeight: baseHeight,
    mountainDensity: 0.58 + rand(seed, 34) * 0.39,
    primaryMountain: {
      x: clamp(0.50 + (rand(seed, 35) - 0.5) * 0.36, 0.18, 0.82),
      y: clamp(0.50 + (rand(seed, 36) - 0.5) * 0.36, 0.18, 0.82),
      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: rangeValue(seed, 43, terrainType.riverRichnessRange),
    bigRiverChance: rangeValue(seed, 44, terrainType.bigRiverChanceRange),
    plainBias: rangeValue(seed, 45, terrainType.plainBiasRange),
  };
}

function buildMountainSystems(template, seed) {
  const systems = [];
  const targetCount = Math.max(8, template.mountainSystemCount ?? 15);
  const baseAngle = template.mountainAngle;
  const angleSpread = template.mountainAngleSpread ?? 0.42;
  const crossSpread = template.mountainCrossSpread ?? 0.70;
  const lengthScale = template.mountainLengthScale ?? 1;
  const widthScale = template.mountainWidthScale ?? 1;
  const heightScale = template.mountainHeightScale ?? 1;
  const isChubu = template.terrainType === "chubu_mountain";
  const isTohoku = template.terrainType === "tohoku_spine";

  // 山脈システムは完全ランダムではなく、複数の広い造山帯に沿って配置する。
  // これにより「方向性はそこそこ揃う」が、「中央一点に集まらない」分布になる。
  const beltCount = Math.max(1, template.mountainBeltCount ?? (targetCount >= 15 ? 4 : 3));
  const belts = [];
  for (let b = 0; b < beltCount; b++) {
    const t = beltCount === 1 ? 0 : (b / (beltCount - 1) - 0.5);
    const angle = baseAngle + (rand(seed, 1000 + b) - 0.5) * angleSpread;
    const axisX = Math.cos(angle);
    const axisY = Math.sin(angle);
    const crossX = Math.cos(angle + Math.PI / 2);
    const crossY = Math.sin(angle + Math.PI / 2);
    const crossOffset = t * crossSpread + (rand(seed, 1010 + b) - 0.5) * (0.10 + crossSpread * 0.08);
    const alongShift = (rand(seed, 1020 + b) - 0.5) * 0.22;
    belts.push({
      angle,
      x: clamp(0.50 + axisX * alongShift / ASPECT + crossX * crossOffset / ASPECT, 0.08, 0.92),
      y: clamp(0.50 + axisY * alongShift + crossY * crossOffset, 0.08, 0.92),
      lengthBias: 0.82 + rand(seed, 1030 + b) * 0.32,
      heightBias: 0.82 + rand(seed, 1040 + b) * 0.42,
    });
  }

  function beltCandidate(k, attempt) {
    const beltIndex = (k + Math.floor(k / beltCount)) % beltCount;
    const belt = belts[beltIndex];
    const perBelt = Math.ceil(targetCount / beltCount);
    const ordinal = Math.floor(k / beltCount);
    const baseT = perBelt <= 1 ? 0 : ordinal / (perBelt - 1) - 0.5;
    const alongJitter = (rand(seed, 1100 + k * 79 + attempt * 11) - 0.5) * (attempt < 4 ? 0.15 : 0.28);
    const crossJitter = (rand(seed, 1200 + k * 79 + attempt * 11) - 0.5) * (attempt < 4 ? crossSpread * 0.22 : crossSpread * 0.40);
    const along = (baseT + alongJitter) * 1.03 * belt.lengthBias;
    const cross = crossJitter;
    const axisX = Math.cos(belt.angle);
    const axisY = Math.sin(belt.angle);
    const crossX = Math.cos(belt.angle + Math.PI / 2);
    const crossY = Math.sin(belt.angle + Math.PI / 2);
    return {
      x: clamp(belt.x + axisX * along / ASPECT + crossX * cross / ASPECT, 0.045, 0.955),
      y: clamp(belt.y + axisY * along + crossY * cross, 0.045, 0.955),
      angle: belt.angle + (rand(seed, 1300 + k * 79 + attempt) - 0.5) * angleSpread * 0.82,
      beltIndex,
    };
  }

  function edgeAwareScore(c) {
    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);
    return Math.min(edgeD, 0.16) * 0.16 + centerD * 0.08;
  }

  if (isTohoku) {
    const centralAngle = baseAngle + (rand(seed, 3330) - 0.5) * 0.06;
    const centralAlong = (rand(seed, 3331) - 0.5) * 0.10;
    const centralCross = (rand(seed, 3332) - 0.5) * 0.045;
    systems.push({
      x: clamp(0.50 + Math.cos(centralAngle) * centralAlong / ASPECT + Math.cos(centralAngle + Math.PI / 2) * centralCross / ASPECT, 0.12, 0.88),
      y: clamp(0.50 + Math.sin(centralAngle) * centralAlong + Math.sin(centralAngle + Math.PI / 2) * centralCross, 0.12, 0.88),
      angle: centralAngle,
      length: (0.78 + rand(seed, 3333) * 0.18) * lengthScale,
      width: (0.12 + rand(seed, 3334) * 0.07) * widthScale,
      height: template.mountainBaseHeight * heightScale * (0.50 + rand(seed, 3335) * 0.15),
      scratchCount: Math.round(20 + rand(seed, 3336) * 10),
      massifness: clamp((template.mountainMassifness ?? 0.16) * 0.55),
      role: "central-primary",
      beltIndex: 0,
      widthVariance: 0.30 + rand(seed, 3337) * 0.46,
      phase: rand(seed, 3338),
    });

    if (rand(seed, 3339) < 0.72) {
      const side = rand(seed, 3340) < 0.5 ? -1 : 1;
      systems.push({
        x: clamp(0.50 + Math.cos(centralAngle) * (centralAlong + side * 0.18) / ASPECT + Math.cos(centralAngle + Math.PI / 2) * (centralCross + side * 0.028) / ASPECT, 0.10, 0.90),
        y: clamp(0.50 + Math.sin(centralAngle) * (centralAlong + side * 0.18) + Math.sin(centralAngle + Math.PI / 2) * (centralCross + side * 0.028), 0.10, 0.90),
        angle: centralAngle + (rand(seed, 3341) - 0.5) * 0.08,
        length: (0.48 + rand(seed, 3342) * 0.20) * lengthScale,
        width: (0.075 + rand(seed, 3343) * 0.055) * widthScale,
        height: template.mountainBaseHeight * heightScale * (0.33 + rand(seed, 3344) * 0.13),
        scratchCount: Math.round(12 + rand(seed, 3345) * 8),
        massifness: clamp((template.mountainMassifness ?? 0.16) * 0.65),
        role: "central-secondary",
        beltIndex: 0,
        widthVariance: 0.24 + rand(seed, 3346) * 0.36,
        phase: rand(seed, 3347),
      });
    }
  }

  for (let k = systems.length; k < targetCount; k++) {
    let best = beltCandidate(k, 0);
    let bestScore = -INF;
    for (let attempt = 0; attempt < 12; attempt++) {
      const c = beltCandidate(k, attempt);
      let minD = 999;
      for (const s of systems) minD = Math.min(minD, distNorm(c.x, c.y, s.x, s.y));
      const score = minD * 1.05 + edgeAwareScore(c) + rand(seed, 2300 + k * 101 + attempt) * 0.035;
      if (score > bestScore) { bestScore = score; best = c; }
    }

    const belt = belts[best.beltIndex];
    const m = clamp(template.mountainMassifness + (rand(seed, 2400 + k) - 0.5) * 0.38);
    const isMassif = m > 0.58;
    const major = k < beltCount || rand(seed, 2500 + k) > 0.72;
    const lengthBaseRange = isChubu
      ? (major ? [0.48, 0.78] : [0.34, 0.58])
      : isTohoku
        ? (major ? [0.44, 0.72] : [0.28, 0.48])
        : (major ? [0.32, 0.52] : [0.21, 0.36]);
    const lengthMassifRange = isChubu
      ? (major ? [0.38, 0.58] : [0.28, 0.44])
      : (major ? [0.23, 0.35] : [0.17, 0.27]);
    const length = lerp(
      lengthBaseRange[0] + rand(seed, 2600 + k) * (lengthBaseRange[1] - lengthBaseRange[0]),
      lengthMassifRange[0] + rand(seed, 2620 + k) * (lengthMassifRange[1] - lengthMassifRange[0]),
      m
    ) * belt.lengthBias * lengthScale;
    const widthRangeA = major ? [0.082, 0.170] : [0.060, 0.120];
    const widthRangeB = major ? [0.150, 0.260] : [0.110, 0.200];
    const width = lerp(
      widthRangeA[0] + rand(seed, 2700 + k) * (widthRangeA[1] - widthRangeA[0]),
      widthRangeB[0] + rand(seed, 2720 + k) * (widthRangeB[1] - widthRangeB[0]),
      m
    ) * widthScale * (0.82 + rand(seed, 2740 + k) * 0.46);
    const heightBase = template.mountainBaseHeight * belt.heightBias * heightScale * (
      major
        ? 0.44 + rand(seed, 2800 + k) * 0.28
        : 0.26 + rand(seed, 2810 + k) * 0.20
    );
    const height = isChubu ? heightBase * 0.82 : heightBase;
    const scratchCount = Math.round(lerp(
      major ? 9 + rand(seed, 2900 + k) * 9 : 6 + rand(seed, 2910 + k) * 6,
      isMassif ? 8 + rand(seed, 2920 + k) * 8 : 6 + rand(seed, 2930 + k) * 6,
      m
    ));

    systems.push({
      x: best.x,
      y: best.y,
      angle: best.angle,
      length,
      width,
      height,
      scratchCount,
      massifness: m,
      role: major ? (k < beltCount ? "primary" : "major") : "minor",
      beltIndex: best.beltIndex,
      widthVariance: 0.18 + rand(seed, 3100 + k) * 0.46,
      phase: rand(seed, 3200 + k),
    });
  }

  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.014 + rand(seed, 2300 + i) * 0.018, 0.024 + rand(seed, 2300 + i) * 0.026, system.massifness) * (0.85 + density * 0.70);
    const height = system.height * (0.060 + density * 0.128 + rand(seed, 2400 + i) * 0.050);
    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 angle = template.coastAngle;
  const axis = (px - 0.5) * Math.cos(angle) * ASPECT + (py - 0.5) * Math.sin(angle);
  const cross = -(px - 0.5) * Math.sin(angle) * ASPECT + (py - 0.5) * Math.cos(angle);
  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 islandNoise = (fbm(px * 420 + 17, py * 420 - 11, seed + 302) - 0.5) * 0.032;
  let pressure = 0;

  if (template.coastStyle === "oceanic_archipelago") {
    // 海洋型: 外洋を強く取り、島列・湾・水道が点在する低標高圧を作る。
    const radial = distNorm(px, py, 0.5, 0.5);
    const outerOcean = smoothstep((radial - 0.30 + wave * 1.15 + bay * 0.85) / 0.18) * 0.96;
    const diagonalChannel = smoothstep((0.13 - Math.abs(cross + wave * 0.82 + islandNoise * 0.70)) / 0.14) * 0.70;
    const openSide = smoothstep((-axis + 0.12 + wave + bay) / 0.23) * 0.82;
    const islandGaps = clamp((valueNoise(px * 780 + 31, py * 780 - 19, seed + 303, 20) - 0.42) * 1.25) * 0.22;
    pressure = clamp(Math.max(outerOcean, diagonalChannel, openSide) + islandGaps);
  } else if (template.coastStyle === "inland_sea") {
    // 瀬戸内型は旧来の大きな内海＋両岸海岸線に戻す。
    // 海面比率はテンプレート側で高めに保ち、微細な島ノイズではなく
    // 連続した水道形状で海を増やす。
    const sideA = smoothstep((-axis + 0.28 + wave + bay) / 0.28);
    const sideB = smoothstep((axis + 0.26 - wave + bay * 0.7) / 0.29);
    const channel = smoothstep((0.082 - Math.abs(cross + wave * 0.68 + islandNoise * 0.55)) / 0.112) * 0.98;
    pressure = Math.max(sideA, sideB, channel);
  } else if (template.coastStyle === "parallel_spine") {
    // 東北型: 左右端または上下端に海を置く。海岸線は脊梁山脈とおおよそ平行。
    // 内海的な中央水路は作らない。
    const edgeA = smoothstep((-axis - 0.26 + wave * 0.42 + bay * 0.35) / 0.20);
    const edgeB = template.twoSidedCoast ? smoothstep((axis - 0.26 - wave * 0.42 + bay * 0.25) / 0.22) * 0.86 : 0;
    pressure = Math.max(edgeA, edgeB);
  } else if (template.coastStyle === "outer_coast") {
    // 中部型: 外縁海を中心にし、内陸へ海が入り込みすぎないようにする。
    const radial = distNorm(px, py, 0.5, 0.5);
    const outer = smoothstep((radial - 0.44 + wave * 0.8 + bay * 0.5) / 0.24);
    const side = smoothstep((-axis + 0.30 + wave) / 0.30) * 0.45;
    pressure = Math.max(outer, side);
  } else if (template.coastStyle === "open_bay") {
    // 関東・濃尾型: 一方向に開いた湾と、その背後の沖積平野を作りやすくする。
    const openSide = smoothstep((-axis + 0.29 + wave + bay) / 0.25);
    const bayMouth = smoothstep((0.22 - Math.abs(cross + wave * 0.8)) / 0.25) * smoothstep((-axis + 0.16 + bay) / 0.22) * 0.68;
    pressure = Math.max(openSide, bayMouth);
  } else {
    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 outerBite = smoothstep((distNorm(px, py, 0.5, 0.5) - 0.54 + islandNoise) / 0.22) * 0.25;
    pressure = Math.max(sideA, sideB, outerBite);
  }

  return { pressure: clamp(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 buildWatershedId(sea, flowTo, flowAccum) {
  const outletKey = new Int32Array(SIZE);
  outletKey.fill(-1);
  const ids = new Int32Array(SIZE);
  ids.fill(-1);
  const outletToId = new Map();
  const quant = 12;

  function keyForOutlet(i) {
    if (i < 0) return -1;
    const x = i % MAP_W;
    const y = Math.floor(i / MAP_W);
    const qx = Math.floor(x / quant);
    const qy = Math.floor(y / quant);
    return qy * 1000 + qx;
  }

  function resolve(start) {
    if (start < 0 || sea[start]) return -1;
    if (outletKey[start] >= 0) return outletKey[start];
    const chain = [];
    const seen = new Set();
    let i = start;
    let key = -1;
    for (let guard = 0; guard < SIZE && i >= 0; guard++) {
      if (sea[i]) { key = keyForOutlet(i); break; }
      if (outletKey[i] >= 0) { key = outletKey[i]; break; }
      if (seen.has(i)) { key = keyForOutlet(i); break; }
      seen.add(i);
      chain.push(i);
      const to = flowTo[i];
      if (to < 0 || to === i) { key = keyForOutlet(i); break; }
      // Major channels should be a watershed's spine, not a sequence of tiny
      // drainage labels. Continue to the coast/outlet even after hitting them.
      i = to;
    }
    if (key < 0 && chain.length) key = keyForOutlet(chain[chain.length - 1]);
    for (const ci of chain) outletKey[ci] = key;
    return key;
  }

  for (let i = 0; i < SIZE; i++) {
    if (sea[i]) continue;
    const key = resolve(i);
    if (key < 0) continue;
    if (!outletToId.has(key)) outletToId.set(key, outletToId.size);
    ids[i] = outletToId.get(key);
  }

  // Very small coastal outlet labels create noisy slivers. Merge them into the
  // strongest neighbouring watershed so natural compartments remain basin-scale.
  const counts = new Int32Array(outletToId.size || 1);
  for (let i = 0; i < SIZE; i++) if (ids[i] >= 0) counts[ids[i]]++;
  const minArea = 18;
  for (let pass = 0; pass < 2; pass++) {
    for (let y = 0; y < MAP_H; y++) {
      for (let x = 0; x < MAP_W; x++) {
        const i = indexOf(x, y);
        const id = ids[i];
        if (id < 0 || counts[id] >= minArea) continue;
        const choices = new Map();
        for (const [nx, ny] of [[x + 1, y], [x - 1, y], [x, y + 1], [x, y - 1]]) {
          if (!inside(nx, ny)) continue;
          const nid = ids[indexOf(nx, ny)];
          if (nid >= 0 && nid !== id) choices.set(nid, (choices.get(nid) || 0) + 1 + (flowAccum[indexOf(nx, ny)] || 0));
        }
        let best = -1, bestScore = -1;
        for (const [nid, score] of choices) if (score > bestScore) { bestScore = score; best = nid; }
        if (best >= 0) { counts[id]--; counts[best]++; ids[i] = best; }
      }
    }
  }
  return ids;
}

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 lineCellsBetween(ax, ay, bx, by) {
  const cells = [];
  const steps = Math.max(Math.abs(bx - ax), Math.abs(by - ay), 1);
  for (let t = 0; t <= steps; t++) {
    const x = Math.round(ax + (bx - ax) * t / steps);
    const y = Math.round(ay + (by - ay) * t / steps);
    if (inside(x, y) && (!cells.length || cells[cells.length - 1][0] !== x || cells[cells.length - 1][1] !== y)) cells.push([x, y]);
  }
  return cells;
}

function meanderRiverPath(path, seed, salt, sea, lake, elevation, slope) {
  if (!path || path.length < 18) return path;
  // Add a visible but controlled meander before valley incision.  The meander
  // diameter is about 3-10 cells.  Mountain reaches are no longer damped;
  // the same low-frequency bend model is applied throughout the course.
  const lengthFactor = clamp(path.length / 190);
  const controlStep = Math.max(4, Math.round(8 - lengthFactor * 3));
  const phase = hash2(seed + salt, 17) * Math.PI * 2;
  const waveCells = 14 + Math.floor(hash2(seed, salt + 31) * 16); // broad wavelength
  const secondaryCells = 28 + Math.floor(hash2(seed + 3, salt + 53) * 24);
  const maxDiameter = 3 + hash2(seed + 5, salt + 71) * 7;
  const baseAmp = maxDiameter * 0.5;
  const controls = [];
  for (let k = 0; k < path.length; k += controlStep) controls.push(k);
  if (controls[controls.length - 1] !== path.length - 1) controls.push(path.length - 1);

  const displacedControls = [];
  for (const k of controls) {
    const [x, y] = path[k];
    if (k === 0 || k === path.length - 1) { displacedControls.push([x, y]); continue; }
    const [px, py] = path[Math.max(0, k - controlStep * 2)];
    const [nx0, ny0] = path[Math.min(path.length - 1, k + controlStep * 2)];
    const tx = nx0 - px;
    const ty = ny0 - py;
    const len = Math.hypot(tx, ty) || 1;
    const normalX = -ty / len;
    const normalY = tx / len;
    const primary = Math.sin(k / waveCells * Math.PI * 2 + phase);
    const secondary = Math.sin(k / secondaryCells * Math.PI * 2 + phase * 0.43) * 0.35;
    const amp = baseAmp * 0.82 * (primary + secondary);
    let mx = Math.round(x + normalX * amp);
    let my = Math.round(y + normalY * amp);
    if (!inside(mx, my)) { displacedControls.push([x, y]); continue; }
    const mi = indexOf(mx, my);
    if (sea[mi] && !lake[mi] && k < path.length - controlStep) { displacedControls.push([x, y]); continue; }
    displacedControls.push([mx, my]);
  }

  const out = [];
  for (let c = 0; c < displacedControls.length - 1; c++) {
    const [ax, ay] = displacedControls[c];
    const [bx, by] = displacedControls[c + 1];
    const line = lineCellsBetween(ax, ay, bx, by);
    for (const cell of line) {
      if (out.length && out[out.length - 1][0] === cell[0] && out[out.length - 1][1] === cell[1]) continue;
      const ci = indexOf(cell[0], cell[1]);
      if (sea[ci] && !lake[ci] && c < displacedControls.length - 3) continue;
      out.push(cell);
    }
  }
  return out.length >= Math.max(8, path.length * 0.42) ? out : path;
}

function scoreRiverPathForDedup(path, flowAccum) {
  let maxFlow = 0;
  let meanFlow = 0;
  for (const [x, y] of path) {
    const f = flowAccum[indexOf(x, y)] || 0;
    maxFlow = Math.max(maxFlow, f);
    meanFlow += f;
  }
  meanFlow /= Math.max(1, path.length);
  return path.length * 0.75 + maxFlow * 90 + meanFlow * 35;
}

function dedupeRiverPaths(paths, flowAccum, sea, lake) {
  // Multiple traces often share, or run one cell beside, the same downstream
  // trunk.  Trim later traces at the first near-confluence so visually only one
  // river occupies a channel, while true tributaries remain visible upstream.
  const sorted = paths
    .map((path) => ({ path, score: scoreRiverPathForDedup(path, flowAccum) }))
    .sort((a, b) => b.score - a.score);
  const occupied = new Uint8Array(SIZE);
  const accepted = [];
  const nearOccupied = (x, y) => {
    for (let dy = -1; dy <= 1; dy++) {
      for (let dx = -1; dx <= 1; dx++) {
        const nx = x + dx;
        const ny = y + dy;
        if (!inside(nx, ny)) continue;
        if (occupied[indexOf(nx, ny)]) return true;
      }
    }
    return false;
  };
  const markNear = (x, y) => {
    for (let dy = -1; dy <= 1; dy++) {
      for (let dx = -1; dx <= 1; dx++) {
        const nx = x + dx;
        const ny = y + dy;
        if (inside(nx, ny)) occupied[indexOf(nx, ny)] = 1;
      }
    }
  };
  for (const item of sorted) {
    const path = item.path;
    if (!path || path.length < 8) continue;
    let joinAt = -1;
    let nearRun = 0;
    for (let k = 0; k < path.length; k++) {
      const [x, y] = path[k];
      if (nearOccupied(x, y) && k > 6) {
        nearRun++;
        if (nearRun >= 2) { joinAt = Math.max(6, k - 1); break; }
      } else {
        nearRun = 0;
      }
    }
    const trimmed = joinAt >= 0 ? path.slice(0, Math.min(path.length, joinAt + 1)) : path;
    let uniqueCells = 0;
    for (const [x, y] of trimmed) if (!nearOccupied(x, y)) uniqueCells++;
    if (trimmed.length < 8 || uniqueCells < Math.max(5, Math.min(18, trimmed.length * 0.38))) continue;
    accepted.push(trimmed);
    for (const [x, y] of trimmed) {
      const i = indexOf(x, y);
      if (!sea[i] || lake[i]) markNear(x, y);
    }
    if (accepted.length >= 62) break;
  }
  return accepted;
}

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 = 26 + Math.floor(template.riverRichness * 15);
  const sources = pickEntities(candidates, { max: desired, minDistance: 8, threshold: 0.27, seed: seed + 12100, jitter: 0.035 });
  const riverPaths = [];
  for (const s of sources) {
    const traced = traceFlowPath(indexOf(s.x, s.y), sea, flowTo);
    const path = meanderRiverPath(traced, seed, s.x * 4096 + s.y, sea, lake, elevation, slope);
    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 visibleRiverPaths = dedupeRiverPaths(riverPaths, flowAccum, sea, lake);
  const longPaths = visibleRiverPaths
    .map((path) => ({ path, score: scoreRiverPathForDedup(path, flowAccum) }))
    .sort((a, b) => b.score - a.score);
  const mainCount = Math.min(longPaths.length, template.terrainType === "kanto_alluvial" ? 5 : 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.040 - (template.bigRiverChance > 0.56 ? 0.038 : 0) - (template.terrainType === "kanto_alluvial" ? 0.035 : 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 visibleRiverPaths) {
    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 kantoMain = main && template.terrainType === "kanto_alluvial";
      const boost = main ? (kantoMain ? 0.66 : 0.54) + downstream * (kantoMain ? 0.50 : 0.42) : 0.30 + downstream * 0.22;
      river[i] = clamp(Math.max(river[i], boost + flowAccum[i] * (main ? (kantoMain ? 0.88 : 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 nonMainPaths = sortedPaths.filter((p) => !mainSet.has(p));
  const tributaryRivers = nonMainPaths.slice(0, 24);
  const smallStreams = nonMainPaths.slice(24, 74);
  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.52 - e) * 1.70);
      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.12 + template.plainBias * 0.34 + riverGate * 0.84));
      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);
  }
}

function rectTerrainProfile(template) {
  const id = String(template?.terrainType || "auto");
  if (id.includes("oceanic")) return {
    base: 0.28, relief: 0.11, ridge: 0.20, ridgeWidth: 18, ridgeSpacing: 70, coast: 0.52, archipelago: 0.18,
    seaQuantile: Math.max(0.95, template.seaRatio ?? 0.965), plain: 0.12, moisture: 0.66, capStart: 0.56, capMax: 0.74,
  };
  if (id.includes("setouchi") || id.includes("archipelago")) return {
    base: 0.43, relief: 0.18, ridge: 0.34, ridgeWidth: 30, ridgeSpacing: 94, coast: 0.25, archipelago: 0.20,
    seaQuantile: Math.max(0.30, template.seaRatio ?? 0.36), plain: 0.34, moisture: 0.58, capStart: 0.72, capMax: 0.94,
  };
  if (id.includes("chubu") || id.includes("mountain")) return {
    base: 0.52, relief: 0.26, ridge: 0.62, ridgeWidth: 36, ridgeSpacing: 108, coast: 0.12, archipelago: 0.03,
    seaQuantile: Math.min(0.22, template.seaRatio ?? 0.18), plain: 0.15, moisture: 0.48, capStart: 0.90, capMax: 1.10,
  };
  if (id.includes("kanto") || id.includes("alluvial")) return {
    base: 0.48, relief: 0.12, ridge: 0.18, ridgeWidth: 42, ridgeSpacing: 130, coast: 0.16, archipelago: 0.04,
    seaQuantile: template.seaRatio ?? 0.13, plain: 0.66, moisture: 0.56, capStart: 0.84, capMax: 1.00,
  };
  if (id.includes("tohoku") || id.includes("spine")) return {
    base: 0.49, relief: 0.19, ridge: 0.46, ridgeWidth: 24, ridgeSpacing: 88, coast: 0.18, archipelago: 0.03,
    seaQuantile: template.seaRatio ?? 0.22, plain: 0.26, moisture: 0.52, capStart: 0.78, capMax: 0.98,
  };
  return {
    base: 0.45, relief: 0.18, ridge: 0.34, ridgeWidth: 32, ridgeSpacing: 100, coast: 0.18, archipelago: 0.06,
    seaQuantile: template.seaRatio ?? 0.20, plain: 0.30, moisture: 0.52, capStart: 0.86, capMax: 1.04,
  };
}

function rectSeed(seed, variant, salt) {
  let h = (seed >>> 0) ^ Math.imul((variant || 0) >>> 0, 0x9e3779b9) ^ (salt >>> 0);
  h ^= h >>> 16;
  h = Math.imul(h, 0x7feb352d) >>> 0;
  h ^= h >>> 15;
  h = Math.imul(h, 0x846ca68b) >>> 0;
  return (h ^ (h >>> 16)) >>> 0;
}

function periodicRidgeField(wx, wy, template, profile, seed) {
  const angle = template.mountainAngle || 0;
  const c = Math.cos(angle);
  const s = Math.sin(angle);
  const u = wx * c + wy * s;
  const v = -wx * s + wy * c;
  const spacing = Math.max(18, profile.ridgeSpacing);
  const shifted = v / spacing + valueNoise(wx, wy, seed ^ 0x654f6d23, 115) * 0.70;
  const nearest = Math.abs((shifted - Math.round(shifted)) * spacing);
  const ridgeCore = Math.exp(-Math.pow(nearest / Math.max(4, profile.ridgeWidth), 2.0));
  const along = valueNoise(u, v, seed ^ 0x27d4eb2f, 86);
  const cut = valueNoise(u, v, seed ^ 0x165667b1, 31);
  return clamp(ridgeCore * (0.62 + along * 0.62) * (0.74 + cut * 0.40));
}

function worldMarinePressure(wx, wy, template, profile, seed) {
  const angle = template.coastAngle || 0;
  const c = Math.cos(angle);
  const s = Math.sin(angle);
  const axis = wx * c + wy * s;
  const cross = -wx * s + wy * c;
  const period = template.coastStyle === "oceanic_archipelago" ? 160 : template.coastStyle === "inland_sea" ? 220 : 300;
  const broad = Math.sin((axis + valueNoise(wx, wy, seed ^ 0xc2b2ae35, 190) * 90) / period * Math.PI * 2);
  const channel = Math.exp(-Math.pow((cross + (valueNoise(wx, wy, seed ^ 0x85ebca6b, 130) - 0.5) * 80) / (profile.ridgeSpacing * 0.85), 2.0));
  const radial = valueNoise(wx, wy, seed ^ 0x9e3779b9, 260);
  let pressure = clamp((broad * 0.5 + 0.5) * profile.coast + channel * profile.coast * 0.62 + radial * profile.coast * 0.52);
  if (template.coastStyle === "oceanic_archipelago") {
    const gap = clamp((fbm(wx * 0.75 + 33, wy * 0.75 - 17, seed ^ 0x3c6ef372) - 0.42) * 2.2);
    pressure = clamp(pressure + gap * profile.archipelago);
  }
  if (template.coastStyle === "open_bay") pressure = clamp(pressure + channel * 0.12);
  return pressure;
}

function classifyRectWater(ctx, fields, seaLevel) {
  const { elevation, sea, ocean, lake } = fields;
  sea.fill(0); ocean.fill(0); lake.fill(0);
  const water = new Uint8Array(ctx.size);
  for (let i = 0; i < ctx.size; i++) water[i] = elevation[i] <= seaLevel ? 1 : 0;
  const seen = new Uint8Array(ctx.size);
  let oceanCells = 0;
  for (let i = 0; i < ctx.size; i++) {
    if (!water[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 % ctx.width;
      const y = Math.floor(cur / ctx.width);
      if (x === 0 || y === 0 || x === ctx.width - 1 || y === ctx.height - 1) touchesEdge = true;
      for (const [nx, ny] of rectNeighbors8(ctx, x, y)) {
        const ni = rectIndexOf(ctx, nx, ny);
        if (!water[ni] || seen[ni]) continue;
        seen[ni] = 1;
        queue.push(ni);
      }
    }
    const isOcean = touchesEdge || cells.length > Math.max(96, ctx.size * 0.018);
    if (isOcean || cells.length >= 20) {
      for (const ci of cells) {
        sea[ci] = 1;
        if (isOcean) ocean[ci] = 1;
        else lake[ci] = 1;
      }
      if (isOcean) oceanCells += cells.length;
    } else {
      for (const ci of cells) elevation[ci] = seaLevel + 0.012;
    }
  }
  return oceanCells;
}

function recomputeRectSlope(ctx, fields) {
  const { elevation, sea, slope } = fields;
  slope.fill(0);
  for (let y = 1; y < ctx.height - 1; y++) {
    for (let x = 1; x < ctx.width - 1; x++) {
      const i = rectIndexOf(ctx, x, y);
      if (sea[i]) continue;
      const gx = elevation[rectIndexOf(ctx, x + 1, y)] - elevation[rectIndexOf(ctx, x - 1, y)];
      const gy = elevation[rectIndexOf(ctx, x, y + 1)] - elevation[rectIndexOf(ctx, x, y - 1)];
      slope[i] = clamp(Math.hypot(gx, gy) * 8.2);
    }
  }
}

function priorityFloodRect(ctx, fields) {
  const { elevation, sea, flowTo } = fields;
  const filled = new Float32Array(elevation);
  const visited = new Uint8Array(ctx.size);
  const heap = new MinHeap();
  let seeds = 0;
  for (let i = 0; i < ctx.size; i++) {
    const x = i % ctx.width;
    const y = Math.floor(i / ctx.width);
    if (sea[i] || x === 0 || y === 0 || x === ctx.width - 1 || y === ctx.height - 1) {
      visited[i] = 1;
      heap.push({ i, f: filled[i] });
      seeds++;
    }
  }
  if (!seeds) return filled;
  while (heap.length) {
    const cur = heap.pop();
    if (!cur || cur.f > filled[cur.i] + 1e-5) continue;
    const x = cur.i % ctx.width;
    const y = Math.floor(cur.i / ctx.width);
    for (const [nx, ny] of rectNeighbors8(ctx, x, y)) {
      const ni = rectIndexOf(ctx, 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] });
    }
  }
  flowTo.fill(-1);
  for (let y = 0; y < ctx.height; y++) {
    for (let x = 0; x < ctx.width; x++) {
      const i = rectIndexOf(ctx, x, y);
      if (sea[i]) continue;
      let best = -1;
      let bestScore = filled[i];
      for (const [nx, ny] of rectNeighbors8(ctx, x, y)) {
        const ni = rectIndexOf(ctx, nx, ny);
        const stepPenalty = (nx !== x && ny !== y) ? 0.000015 : 0;
        const score = filled[ni] + stepPenalty + hash2(ctx.originX + nx, ctx.originY + ny, 9000) * 0.000002;
        if (score < bestScore - 0.000001 || sea[ni]) {
          bestScore = score;
          best = ni;
          if (sea[ni]) break;
        }
      }
      flowTo[i] = best;
    }
  }
  return filled;
}

function computeRectFlowAccumulation(ctx, fields, filled) {
  const { sea, flowTo, flowAccum } = fields;
  const area = new Float32Array(ctx.size);
  const order = [];
  for (let i = 0; i < ctx.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 < ctx.size; i++) if (!sea[i]) maxArea = Math.max(maxArea, area[i]);
  for (let i = 0; i < ctx.size; i++) flowAccum[i] = sea[i] ? 0 : clamp(Math.pow(area[i] / maxArea, 0.42));
}

function rectStableId(seed, wx, wy, salt) {
  const x = Math.floor(wx) | 0;
  const y = Math.floor(wy) | 0;
  let h = (seed >>> 0) ^ Math.imul(x, 0x9e3779b1) ^ Math.imul(y, 0x85ebca77) ^ (salt >>> 0);
  h ^= h >>> 16;
  h = Math.imul(h, 0x7feb352d) >>> 0;
  h ^= h >>> 15;
  h = Math.imul(h, 0x846ca68b) >>> 0;
  return (h ^ (h >>> 16)) & 0x7fffffff;
}

function traceRectSink(ctx, start, fields, maxSteps = 4096) {
  const { sea, flowTo } = fields;
  let i = start;
  let last = i;
  const seen = new Set();
  for (let step = 0; step < maxSteps; step++) {
    if (i < 0 || i >= ctx.size || seen.has(i)) break;
    seen.add(i);
    last = i;
    if (sea[i]) break;
    const next = flowTo[i];
    if (next < 0 || next === i) break;
    i = next;
  }
  return last;
}

function buildRectWatershedId(ctx, fields, seed) {
  const { sea, flowAccum, watershedId } = fields;
  if (!watershedId) return { watershedCount: 0 };
  watershedId.fill(-1);
  const sinkToId = new Map();
  let watershedCount = 0;
  for (let i = 0; i < ctx.size; i++) {
    if (sea[i]) continue;
    const sink = traceRectSink(ctx, i, fields);
    const sx = sink % ctx.width;
    const sy = Math.floor(sink / ctx.width);
    const wx = ctx.originX + sx;
    const wy = ctx.originY + sy;
    const coarseX = Math.round(wx / 12);
    const coarseY = Math.round(wy / 12);
    const key = `${coarseX},${coarseY}`;
    let id = sinkToId.get(key);
    if (!Number.isFinite(id)) {
      id = 50000000 + rectStableId(seed, coarseX, coarseY, 0x51ed270b) % 40000000;
      sinkToId.set(key, id);
      watershedCount++;
    }
    watershedId[i] = id;
  }
  // Merge tiny or noisy drainage islands into their strongest neighbor.
  for (let pass = 0; pass < 2; pass++) {
    const changes = [];
    for (let y = 1; y < ctx.height - 1; y++) {
      for (let x = 1; x < ctx.width - 1; x++) {
        const i = rectIndexOf(ctx, x, y);
        if (sea[i] || watershedId[i] < 0) continue;
        const counts = new Map();
        for (const [nx, ny] of rectNeighbors8(ctx, x, y)) {
          const ni = rectIndexOf(ctx, nx, ny);
          const id = watershedId[ni];
          if (id < 0) continue;
          counts.set(id, (counts.get(id) || 0) + 1 + (flowAccum[ni] || 0));
        }
        let best = watershedId[i];
        let bestScore = counts.get(best) || 0;
        for (const [id, score] of counts) if (score > bestScore) { best = id; bestScore = score; }
        if (best !== watershedId[i] && bestScore >= 5.5) changes.push([i, best]);
      }
    }
    for (const [i, id] of changes) watershedId[i] = id;
    if (!changes.length) break;
  }
  return { watershedCount };
}

function buildRectNaturalRegions(ctx, fields, seed, template) {
  const { sea, ridgeField, valleyField, basinField, flowAccum, naturalBarrierScore, watershedId, naturalCompartmentId, regionId } = fields;
  if (!naturalCompartmentId || !regionId) return { naturalCompartmentCount: 0, regionCount: 0 };
  naturalCompartmentId.fill(-1);
  regionId.fill(-1);
  const type = String(template?.terrainType || "auto");
  const spacing = type.includes("oceanic") ? 30 : type.includes("kanto") ? 44 : type.includes("chubu") ? 34 : 38;
  const coarseSpacing = spacing * 2.55;
  const seeds = [];
  const gx0 = Math.floor((ctx.originX - spacing) / spacing) - 1;
  const gx1 = Math.ceil((ctx.originX + ctx.width + spacing) / spacing) + 1;
  const gy0 = Math.floor((ctx.originY - spacing) / spacing) - 1;
  const gy1 = Math.ceil((ctx.originY + ctx.height + spacing) / spacing) + 1;
  for (let gy = gy0; gy <= gy1; gy++) {
    for (let gx = gx0; gx <= gx1; gx++) {
      const jitterX = (hash2(gx, gy, seed ^ 0x6a09e667) - 0.5) * spacing * 0.74;
      const jitterY = (hash2(gx, gy, seed ^ 0xbb67ae85) - 0.5) * spacing * 0.74;
      const wx = gx * spacing + spacing * 0.5 + jitterX;
      const wy = gy * spacing + spacing * 0.5 + jitterY;
      const lx = Math.round(wx - ctx.originX);
      const ly = Math.round(wy - ctx.originY);
      let viability = 0.8;
      if (rectInside(ctx, lx, ly)) {
        const i = rectIndexOf(ctx, lx, ly);
        viability += (basinField[i] || 0) * 0.25 + (valleyField[i] || 0) * 0.16 - (ridgeField[i] || 0) * 0.14;
        if (sea[i]) viability -= 1.2;
      }
      if (viability < 0.18 && hash2(gx, gy, seed ^ 0x3c6ef372) < 0.82) continue;
      seeds.push({
        wx,
        wy,
        id: 40000000 + rectStableId(seed, gx, gy, 0xb5c0fbcf) % 42000000,
        coarseId: 30000000 + rectStableId(seed, Math.floor((gx * spacing) / coarseSpacing), Math.floor((gy * spacing) / coarseSpacing), 0xc2b2ae35) % 42000000,
      });
    }
  }
  if (!seeds.length) return { naturalCompartmentCount: 0, regionCount: 0 };
  for (let y = 0; y < ctx.height; y++) {
    for (let x = 0; x < ctx.width; x++) {
      const i = rectIndexOf(ctx, x, y);
      if (sea[i]) continue;
      const wx = ctx.originX + x;
      const wy = ctx.originY + y;
      let best = seeds[0];
      let bestScore = Infinity;
      const barrier = (naturalBarrierScore[i] || 0) + (ridgeField[i] || 0) * 0.55 + (flowAccum[i] || 0) * 0.18;
      const basinBonus = (basinField[i] || 0) * 0.18 + (valleyField[i] || 0) * 0.10;
      for (const s of seeds) {
        const dx = (wx - s.wx) * 1.05;
        const dy = wy - s.wy;
        const d = Math.hypot(dx, dy);
        const tileNoise = (valueNoise(wx + s.wx * 0.13, wy + s.wy * 0.13, seed ^ 0xa54ff53a, 52) - 0.5) * spacing * 0.34;
        const watershedPenalty = watershedId?.[i] >= 0 ? ((watershedId[i] ^ s.id) & 7) * 0.16 : 0;
        const score = d + barrier * spacing * 0.42 - basinBonus * spacing * 0.32 + tileNoise + watershedPenalty;
        if (score < bestScore) { bestScore = score; best = s; }
      }
      naturalCompartmentId[i] = best.id;
      regionId[i] = best.coarseId;
    }
  }
  for (let pass = 0; pass < 2; pass++) {
    const changes = [];
    for (let y = 1; y < ctx.height - 1; y++) {
      for (let x = 1; x < ctx.width - 1; x++) {
        const i = rectIndexOf(ctx, x, y);
        if (sea[i]) continue;
        if ((ridgeField[i] || 0) > 0.78) continue;
        const counts = new Map();
        for (const [nx, ny] of rectNeighbors8(ctx, x, y)) {
          const ni = rectIndexOf(ctx, nx, ny);
          const id = naturalCompartmentId[ni];
          if (id < 0) continue;
          counts.set(id, (counts.get(id) || 0) + 1 + (basinField[ni] || 0) * 0.3);
        }
        let best = naturalCompartmentId[i];
        let bestScore = counts.get(best) || 0;
        for (const [id, score] of counts) if (score > bestScore) { best = id; bestScore = score; }
        if (best !== naturalCompartmentId[i] && bestScore >= 5.8) changes.push([i, best]);
      }
    }
    for (const [i, id] of changes) naturalCompartmentId[i] = id;
    if (!changes.length) break;
  }
  const nset = new Set();
  const rset = new Set();
  for (let i = 0; i < ctx.size; i++) {
    if (naturalCompartmentId[i] >= 0) nset.add(naturalCompartmentId[i]);
    if (regionId[i] >= 0) rset.add(regionId[i]);
  }
  return { naturalCompartmentCount: nset.size, regionCount: rset.size };
}

function traceRectFlowPath(start, ctx, fields, maxSteps = 1200) {
  const { sea, flowTo } = fields;
  let i = start;
  const path = [];
  const seen = new Set();
  for (let step = 0; step < maxSteps; step++) {
    if (i < 0 || i >= ctx.size || seen.has(i)) break;
    seen.add(i);
    const x = i % ctx.width;
    const y = Math.floor(i / ctx.width);
    path.push([ctx.originX + x, ctx.originY + y]);
    if (sea[i]) break;
    const next = flowTo[i];
    if (next < 0 || next === i) break;
    i = next;
  }
  return path;
}

function scoreRectRiverPath(path, ctx, fields) {
  let score = 0;
  for (const [wx, wy] of path) {
    const x = wx - ctx.originX;
    const y = wy - ctx.originY;
    if (!rectInside(ctx, x, y)) continue;
    const i = rectIndexOf(ctx, x, y);
    score += (fields.flowAccum[i] || 0) + (fields.river[i] || 0) * 0.7;
  }
  return score;
}

function buildRectRiverPaths(ctx, fields, seed, template) {
  const { sea, flowAccum, river, erosionField } = fields;
  const candidates = [];
  const threshold = template?.terrainType === "oceanic_archipelago" ? 0.52 : template?.terrainType === "kanto_alluvial" ? 0.46 : 0.50;
  for (let i = 0; i < ctx.size; i++) {
    if (sea[i] || flowAccum[i] < threshold) continue;
    const x = i % ctx.width;
    const y = Math.floor(i / ctx.width);
    let upstream = 0;
    for (const [nx, ny] of rectNeighbors8(ctx, x, y)) {
      const ni = rectIndexOf(ctx, nx, ny);
      if (fields.flowTo[ni] === i) upstream++;
    }
    const sourceBias = hash2(ctx.originX + x, ctx.originY + y, seed ^ 0x1f123bb5);
    if (upstream <= 1 || sourceBias > 0.78) candidates.push({ i, score: flowAccum[i] + sourceBias * 0.12 });
  }
  candidates.sort((a, b) => b.score - a.score);
  const accepted = [];
  const occupied = new Set();
  const desired = Math.min(72, Math.max(8, Math.floor(ctx.size / 900)));
  for (const c of candidates) {
    if (accepted.length >= desired) break;
    const path = traceRectFlowPath(c.i, ctx, fields);
    if (path.length < 8) continue;
    const keyHits = path.reduce((n, [wx, wy], k) => k % 3 === 0 && occupied.has(`${wx},${wy}`) ? n + 1 : n, 0);
    if (keyHits > Math.max(5, path.length * 0.18)) continue;
    const score = scoreRectRiverPath(path, ctx, fields);
    if (score < 4.2) continue;
    accepted.push({ path, score });
    for (const [wx, wy] of path) occupied.add(`${wx},${wy}`);
  }
  accepted.sort((a, b) => b.score - a.score);
  const mainRivers = accepted.slice(0, Math.max(1, Math.min(10, Math.round(accepted.length * 0.25)))).map((r) => r.path);
  const tributaryRivers = accepted.slice(mainRivers.length, mainRivers.length + 28).map((r) => r.path);
  const smallStreams = accepted.slice(mainRivers.length + 28, mainRivers.length + 56).map((r) => r.path);
  for (const group of [mainRivers, tributaryRivers, smallStreams]) {
    const boost = group === mainRivers ? 0.72 : group === tributaryRivers ? 0.48 : 0.28;
    for (const path of group) {
      for (const [wx, wy] of path) {
        const x = wx - ctx.originX;
        const y = wy - ctx.originY;
        if (!rectInside(ctx, x, y)) continue;
        const i = rectIndexOf(ctx, x, y);
        if (sea[i]) continue;
        river[i] = clamp(Math.max(river[i], boost + (flowAccum[i] || 0) * 0.42));
        if (erosionField) erosionField[i] = clamp((erosionField[i] || 0) + river[i] * 0.18);
      }
    }
  }
  return { riverPaths: accepted.map((r) => r.path), mainRivers, tributaryRivers, smallStreams };
}

function deriveRectTerrainFields(ctx, fields, seaLevel) {
  const {
    elevation, sea, river, flowAccum, floodplain, plain, agriculture, ridgeField, valleyField, basinField,
    coastalLowland, erosionField, depositionField, depositionalLowland, alluvialFanField, deltaField,
    naturalBarrierScore, portSuitability, crossingSuitability, passSuitability, slope, moisture,
  } = fields;
  for (let i = 0; i < ctx.size; i++) {
    if (sea[i]) {
      river[i] = 0; plain[i] = 0; agriculture[i] = 0; naturalBarrierScore[i] = 0;
      continue;
    }
    const low = clamp((0.48 - elevation[i]) * 2.2);
    const flat = clamp(1 - slope[i] * 2.3);
    const coast = clamp((elevation[i] - seaLevel) * 18);
    river[i] = flowAccum[i] > 0.58 ? clamp((flowAccum[i] - 0.52) * 2.1 + (0.22 - slope[i]) * 0.40) : 0;
    floodplain[i] = clamp(river[i] * 0.72 + low * flat * 0.24);
    plain[i] = clamp(flat * (low * 0.78 + basinField[i] * 0.38 + floodplain[i] * 0.35));
    agriculture[i] = clamp(plain[i] * 0.72 + moisture[i] * 0.22 - slope[i] * 0.22);
    coastalLowland[i] = clamp((1 - coast) * flat * 0.90);
    erosionField[i] = clamp(slope[i] * 0.55 + river[i] * 0.34 + ridgeField[i] * 0.22);
    depositionField[i] = clamp(floodplain[i] * 0.58 + coastalLowland[i] * 0.34 + plain[i] * 0.18);
    depositionalLowland[i] = clamp(depositionField[i] * flat);
    alluvialFanField[i] = clamp(river[i] * slope[i] * 1.8);
    deltaField[i] = clamp(river[i] * coastalLowland[i] * 1.2);
    naturalBarrierScore[i] = clamp(ridgeField[i] * 0.72 + slope[i] * 0.42 + river[i] * 0.24);
    crossingSuitability[i] = clamp(flat * (1 - river[i] * 0.65) + plain[i] * 0.24);
    passSuitability[i] = clamp((1 - ridgeField[i]) * 0.55 + valleyField[i] * 0.40 - slope[i] * 0.15);
    portSuitability[i] = clamp(coastalLowland[i] * 0.65 + plain[i] * 0.22 - slope[i] * 0.26);
  }
}

export function generateTerrainRect(options = {}) {
  const seed = Number.isFinite(options.seed) ? options.seed >>> 0 : 0;
  const variant = Number.isFinite(options.variant) ? Math.max(0, Math.floor(options.variant)) >>> 0 : 0;
  const ctx = options.rectContext || createRectContext(options);
  const rectSeedValue = rectSeed(seed, variant, 0x5489a1f3);
  const terrainTemplate = buildTerrainTemplate(rectSeedValue, options);
  const profile = rectTerrainProfile(terrainTemplate);
  const fields = createRectTerrainFields(ctx);
  const {
    elevation, moisture, ridgeField, valleyField, basinField, coastalLowland, arcSpineField, branchRidgeField,
    visibleRavineField, surfaceTextureField,
  } = fields;

  for (let y = 0; y < ctx.height; y++) {
    for (let x = 0; x < ctx.width; x++) {
      const i = rectIndexOf(ctx, x, y);
      const wx = ctx.originX + x;
      const wy = ctx.originY + y;
      const broad = (fbm(wx * 0.58, wy * 0.58, rectSeedValue ^ 0x9e3779b9) - 0.5) * profile.relief;
      const regional = (valueNoise(wx, wy, rectSeedValue ^ 0x85ebca6b, 58) - 0.5) * profile.relief * 0.72;
      const detail = (valueNoise(wx, wy, rectSeedValue ^ 0xc2b2ae35, 19) - 0.5) * profile.relief * 0.22;
      const ridge = periodicRidgeField(wx, wy, terrainTemplate, profile, rectSeedValue ^ 0x27d4eb2f);
      const marine = worldMarinePressure(wx, wy, terrainTemplate, profile, rectSeedValue ^ 0x165667b1);
      const basin = clamp((valueNoise(wx, wy, rectSeedValue ^ 0xd3a2646c, 120) - 0.36) * 1.65) * profile.plain;
      const valley = clamp((1 - ridge) * (valueNoise(wx, wy, rectSeedValue ^ 0xfd7046c5, 42) - 0.42) * 1.7);
      const archipelago = terrainTemplate.terrainType === "oceanic_archipelago" || terrainTemplate.coastStyle === "inland_sea"
        ? clamp((fbm(wx * 0.72 + 49, wy * 0.72 - 31, rectSeedValue ^ 0x94d049bb) - 0.44) * 2.2) * profile.archipelago
        : 0;
      let e = profile.base + broad + regional + detail + ridge * profile.ridge + archipelago - marine + basin * 0.10;
      if (terrainTemplate.terrainType === "kanto_alluvial") e -= basin * 0.075;
      if (terrainTemplate.terrainType === "oceanic_archipelago") e -= marine * 0.16;
      e = softCapElevation(e, profile.capStart, profile.capMax);
      elevation[i] = clamp(e, 0.025, profile.capMax);
      ridgeField[i] = clamp(ridge * (0.62 + profile.ridge));
      branchRidgeField[i] = clamp(ridge * 0.82 + detail * 0.60);
      arcSpineField[i] = clamp(ridge * 0.90);
      valleyField[i] = clamp(valley + (1 - ridge) * marine * 0.20);
      basinField[i] = clamp(basin + valley * 0.35);
      coastalLowland[i] = clamp(marine * 0.82 + basin * 0.25);
      moisture[i] = clamp(profile.moisture + marine * 0.22 + basin * 0.15 - elevation[i] * 0.22 + (fbm(wx * 0.85, wy * 0.85, rectSeedValue ^ 0xa0761d65) - 0.5) * 0.13);
      visibleRavineField[i] = clamp(Math.abs(detail) * ridge * 1.9 + valley * 0.25);
      surfaceTextureField[i] = clamp(Math.abs(broad) * 0.55 + Math.abs(detail) * 1.3 + ridge * 0.22);
    }
  }

  const seaLevel = clamp(
    Number.isFinite(options.seaLevel) ? options.seaLevel : rectQuantile(elevation, profile.seaQuantile),
    0.13,
    terrainTemplate.terrainType === "oceanic_archipelago" ? 0.72 : 0.50
  );
  const oceanCells = classifyRectWater(ctx, fields, seaLevel);
  recomputeRectSlope(ctx, fields);
  const filled = priorityFloodRect(ctx, fields);
  computeRectFlowAccumulation(ctx, fields, filled);
  const watershedDebug = buildRectWatershedId(ctx, fields, rectSeedValue ^ 0x51ed270b);
  deriveRectTerrainFields(ctx, fields, seaLevel);
  const riverNetwork = buildRectRiverPaths(ctx, fields, rectSeedValue ^ 0x1f123bb5, terrainTemplate);
  const naturalDebug = buildRectNaturalRegions(ctx, fields, rectSeedValue ^ 0xb5c0fbcf, terrainTemplate);

  let landCount = 0;
  let mountainCount = 0;
  let plainCount = 0;
  for (let i = 0; i < ctx.size; i++) {
    if (fields.sea[i]) continue;
    landCount++;
    if (fields.elevation[i] > 0.56 || fields.ridgeField[i] > 0.52) mountainCount++;
    if (fields.plain[i] > 0.36) plainCount++;
  }

  return {
    rectContext: ctx,
    originX: ctx.originX,
    originY: ctx.originY,
    width: ctx.width,
    height: ctx.height,
    size: ctx.size,
    terrainTemplate,
    seaLevel,
    ...fields,
    terrainDebug: {
      terrainType: terrainTemplate.terrainType,
      terrainTypeLabel: terrainTemplate.terrainTypeLabel,
      coastStyle: terrainTemplate.coastStyle,
      rectNative: true,
      originX: ctx.originX,
      originY: ctx.originY,
      width: ctx.width,
      height: ctx.height,
      variant,
      seaRatio: fields.sea.reduce((sum, value) => sum + value, 0) / Math.max(1, ctx.size),
      landCount,
      oceanCells,
      mountainRatio: mountainCount / Math.max(1, landCount),
      plainRatio: plainCount / Math.max(1, landCount),
      watershedCount: watershedDebug.watershedCount,
      naturalCompartmentCount: naturalDebug.naturalCompartmentCount,
      regionCount: naturalDebug.regionCount,
      mainRiverCount: riverNetwork.mainRivers.length,
      tributaryRiverCount: riverNetwork.tributaryRivers.length,
      smallStreamCount: riverNetwork.smallStreams.length,
    },
    ...riverNetwork,
  };
}

export function finalizeRectTerrainForFixedMap(seed, terrain, options = {}) {
  if (!terrain || terrain.width !== MAP_W || terrain.height !== MAP_H || terrain.size !== SIZE) {
    throw new Error(`finalizeRectTerrainForFixedMap requires ${MAP_W}x${MAP_H} terrain, got ${terrain?.width}x${terrain?.height}`);
  }
  const {
    elevation, slope, sea, river, ridgeField, valleyField, basinField, coastalLowland, flowAccum,
    plain, agriculture, watershedId, landMask: existingLandMask, prefectureMask: existingPrefectureMask,
  } = terrain;
  const prefectureMask = existingPrefectureMask || makePrefectureMask(seed, sea, elevation, slope, river);
  const landMask = existingLandMask || new Uint8Array(SIZE);
  if (!existingLandMask) {
    for (let i = 0; i < SIZE; i++) landMask[i] = sea[i] ? 0 : 1;
  }
  const zeroDensity = new Float32Array(SIZE);
  const zeroLanduse = new Int8Array(SIZE);
  const landCount = landMask.reduce((sum, value, i) => sum + (value && !sea[i] ? 1 : 0), 0);
  const natural = buildNaturalCompartments(
    landMask, sea, elevation, slope, river, ridgeField, valleyField, basinField, coastalLowland, flowAccum,
    null, plain, agriculture, zeroDensity, zeroLanduse,
    {
      seed: (seed + 17003) >>> 0,
      watershedId,
      targetCompartmentCount: clamp(Math.round(landCount / 45), 70, 360),
    }
  );
  const prefectureBorder = extractMaskBorder(prefectureMask, sea);
  const terrainDebug = {
    ...(terrain.terrainDebug || {}),
    rectNativeInitialTerrain: true,
    rectInitialOriginX: terrain.originX || 0,
    rectInitialOriginY: terrain.originY || 0,
    sharedNaturalCompartmentLayer: true,
    naturalCompartmentCount: natural.compartments?.filter?.((unit) => unit && unit.area > 0).length || 0,
  };
  return {
    ...terrain,
    prefectureMask,
    landMask,
    prefectureBorder,
    naturalBarrierScore: natural.naturalBarrierScore || terrain.naturalBarrierScore,
    naturalCompartmentId: natural.compartmentId,
    naturalCompartments: natural.compartments,
    terrainDebug,
  };
}

export function generateInitialTerrainRect(seed, options = {}) {
  const variant = Number.isFinite(options.initialVariant) ? Math.max(0, Math.floor(options.initialVariant)) : 0;
  const terrain = generateTerrainRect({
    ...options,
    seed,
    variant,
    originX: 0,
    originY: 0,
    width: MAP_W,
    height: MAP_H,
    name: "initial-full-map",
  });
  return finalizeRectTerrainForFixedMap(seed, terrain, options);
}


export function generateTerrainAndRivers(seed, options = {}) {
  const generationContext = options.generationContext || {};
  const originX = Math.floor(Number.isFinite(options.originX) ? options.originX : (Number.isFinite(generationContext.originX) ? generationContext.originX : 0));
  const originY = Math.floor(Number.isFinite(options.originY) ? options.originY : (Number.isFinite(generationContext.originY) ? generationContext.originY : 0));
  const variant = Math.max(0, Math.floor(Number.isFinite(options.variant) ? options.variant : (Number.isFinite(generationContext.variant) ? generationContext.variant : 0))) >>> 0;
  const worldNative = options.worldNative === true || generationContext.worldNative === true;
  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, options);
  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 wx = originX + x;
      const wy = originY + y;
      const { px, py } = normalizeCoord(x, y);
      const terrainLarge = (fbm(wx * 0.65, wy * 0.65, seed + 1) - 0.5) * 0.23;
      const terrainRegional = (valueNoise(wx * 0.8, wy * 0.8, seed + 2, 42) - 0.5) * 0.16;
      const { pressure: coastPressure } = computeCoastLower(px, py, terrainTemplate, seed);
      let e = 0.42 + terrainLarge + terrainRegional - coastPressure * (terrainTemplate.coastStrength ?? 0.96) * (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(1.70, 1.16, system.massifness)) * system.height * lerp(0.22, 0.34, 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));
      }
      let macro;
      let scratch;
      if (terrainTemplate.terrainType === "tohoku_spine") {
        const dx = (px - 0.5) * ASPECT;
        const dy = py - 0.5;
        const { u, v } = rotate(dx, dy, terrainTemplate.mountainAngle);
        const warp = (valueNoise(wx * 0.50, wy * 0.50, seed + 504, 28) - 0.5) * 18;
        macro = (fbm(u * 520 + warp, v * 1850 - warp * 0.35, seed + 500) - 0.5) * 2;
        scratch = (fbm(u * 920 + warp * 0.8, v * 2300 + warp * 0.25, seed + 502) - 0.5) * 2;
        const passBreak = clamp((valueNoise(u * 760 + 17, v * 1800 - 11, seed + 505, 18) - 0.62) * 2.6);
        const lateralBranch = clamp((valueNoise(v * 1700 - 3, u * 540 + 5, seed + 506, 16) - 0.54) * 2.1);
        e += lateralBranch * mountainMaskMax * 0.022;
        e -= passBreak * mountainMaskMax * 0.052;
      } else {
        macro = (fbm(wx * terrainTemplate.macroNoiseScale * 48, wy * terrainTemplate.macroNoiseScale * 48, seed + 500) - 0.5) * 2;
        scratch = (fbm(wx * 2.2, wy * 2.2, seed + 502) - 0.5) * 2;
      }
      const global = (valueNoise(wx * 0.23, wy * 0.23, seed + 501, 38) - 0.5) * 2;
      e += macro * terrainTemplate.macroNoiseStrength * (0.38 + mountainMaskMax * 0.80);
      e += global * 0.020;
      e += scratch * terrainTemplate.scratchNoiseStrength * mountainMaskMax;
      if (terrainTemplate.terrainType === "tohoku_spine") {
        // Keep the broad Ou/backbone footprint, but compress peak height so the
        // range reads as a long Japanese spine rather than an alpine wall.
        const high = Math.max(0, e - 0.48);
        e -= high * clamp(0.18 + mountainMaskMax * 0.22, 0.18, 0.40);
      }
      if (terrainTemplate.terrainType === "setouchi_inland_sea") {
        // Setouchi should read as sea-dominant, with many compact wooded island
        // backbones rather than broad continental ridges.  The small-massif term
        // is band-limited so it forms believable islands, not one-cell speckle.
        const lowHillNoise = clamp((fbm(wx * 0.95 + 17, wy * 0.95 - 23, seed + 571) - 0.39) * 2.7);
        const islandMassif = clamp((fbm(wx * 1.85 - 31, wy * 1.85 + 19, seed + 573) - 0.50) * 3.4);
        const islandBackbone = clamp((valueNoise(wx * 2.8 + 7, wy * 2.8 - 11, seed + 574, 9) - 0.54) * 3.2);
        const coastalIslandBias = clamp(coastPressure * 0.64 + mountainMaskMax * 0.46 + lowHillNoise * 0.24);
        e += lowHillNoise * 0.090;
        e += islandMassif * coastalIslandBias * 0.105;
        e += islandBackbone * coastalIslandBias * 0.045;
        e -= clamp((coastPressure - 0.34) * 1.55) * 0.044;
        const high = Math.max(0, e - 0.60);
        e -= high * 0.48;
      }
      if (terrainTemplate.terrainType === "oceanic_archipelago") {
        // 海洋型は大きな山塊ではなく、島列を読むための低い起伏を点在させる。
        const islandCore = clamp((fbm(wx * 0.82 + 41, wy * 0.82 - 29, seed + 572) - 0.46) * 2.7);
        const islandChain = clamp(mountainMaskMax * 0.92 + islandCore * 0.54 - coastPressure * 0.24);
        e += islandChain * 0.135;
        e -= clamp((coastPressure - 0.38) * 1.55) * 0.040;
        const high = Math.max(0, e - 0.56);
        e -= high * 0.52;
      }
      const softCapStart = terrainTemplate.terrainType === "tohoku_spine" ? 0.78 : terrainTemplate.terrainType === "setouchi_inland_sea" ? 0.72 : terrainTemplate.terrainType === "oceanic_archipelago" ? 0.62 : terrainTemplate.terrainType === "chubu_mountain" ? 0.88 : 0.91;
      const softCapMax = terrainTemplate.terrainType === "tohoku_spine" ? 0.96 : terrainTemplate.terrainType === "setouchi_inland_sea" ? 0.92 : terrainTemplate.terrainType === "oceanic_archipelago" ? 0.84 : 1.08;
      elevation[i] = clamp(softCapElevation(e, softCapStart, softCapMax), 0.025, softCapMax);
      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(wx * 1.1, wy * 1.1, seed + 503) - 0.5) * 0.16);
    }
  }

  let seaLevel = quantile(elevation, terrainTemplate.seaRatio);
  seaLevel = clamp(seaLevel, 0.14, terrainTemplate.terrainType === "oceanic_archipelago" ? 0.72 : 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 watershedId = buildWatershedId(sea, flowTo, 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 landMask = new Uint8Array(SIZE);
  for (let i = 0; i < SIZE; i++) landMask[i] = sea[i] ? 0 : 1;
  const zeroDensity = new Float32Array(SIZE);
  const zeroLanduse = new Int8Array(SIZE);
  const natural = buildNaturalCompartments(
    landMask, sea, elevation, slope, river, ridgeField, valleyField, basinField, coastalLowland, flowAccum,
    null, plain, agriculture, zeroDensity, zeroLanduse,
    { seed: seed + 17003, watershedId, targetCompartmentCount: clamp(Math.round((SIZE - sea.reduce((sum, value) => sum + value, 0)) / 45), 70, 360) }
  );
  const sharedNaturalBarrierScore = natural.naturalBarrierScore || naturalBarrierScore;
  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.56 || ridgeField[i] > 0.52) mountainCount++;
    if (plain[i] > 0.36) plainCount++;
    if (arcSpineField[i] > 0.55) { primarySpineStrength += arcSpineField[i]; spineSamples++; }
  }
  primarySpineStrength /= Math.max(1, spineSamples);
  const terrainDebug = {
    terrainType: terrainTemplate.terrainType,
    terrainTypeLabel: terrainTemplate.terrainTypeLabel,
    coastStyle: terrainTemplate.coastStyle,
    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,
    originX,
    originY,
    width: MAP_W,
    height: MAP_H,
    variant,
    worldNative,
  };

  return {
    originX,
    originY,
    width: MAP_W,
    height: MAP_H,
    generationContext: { ...generationContext, originX, originY, width: MAP_W, height: MAP_H, variant, worldNative },
    terrainTemplate,
    seaLevel,
    elevation,
    moisture,
    slope,
    sea,
    ocean,
    lake,
    river,
    floodplain,
    plain,
    agriculture,
    ridgeField,
    valleyField,
    visibleRavineField,
    surfaceTextureField,
    basinField,
    coastalLowland,
    flowAccum,
    watershedId,
    erosionField,
    depositionField,
    arcSpineField,
    branchRidgeField,
    depositionalLowland,
    alluvialFanField,
    deltaField,
    naturalBarrierScore: sharedNaturalBarrierScore,
    portSuitability,
    crossingSuitability,
    passSuitability,
    prefectureMask,
    landMask,
    prefectureBorder,
    naturalCompartmentId: natural.compartmentId,
    naturalCompartments: natural.compartments,
    terrainDebug,
    riverPaths,
    mainRivers,
    tributaryRivers,
    smallStreams,
  };
}