import { INF, MAP_H, MAP_W, SIZE, MinHeap, clamp, fbm, hash2, indexOf, inside, pickEntities, rand, valueNoise, xyOf } from "./mapUtils.js";
import { distanceToNearest, influenceFromPaths, influenceFromPoints, samplePath } from "./mapGeneratorHelpers.js";
import { LANDUSE } from "./landuseCodes.js";
import { buildDensityFlowRoadTransportSystem, createPathInfluenceCache, packDebugField, pathAverageField, pathLengthCells, routeQualityAcceptable } from "./mapTransport.js";
import { buildUnifiedRailODNetwork } from "./mapTransportOD.js";

// Lightweight Human Geography V2
// --------------------------------
// This replaces the heavy iterative human stage with a sparse skeleton + raster
// synthesis model:
//   1. build terrain-derived human context once
//   2. place villages/towns/cities by region quotas
//   3. make sparse approximate transport paths without full-resolution A*
//   4. synthesize population and land-use fields in one raster pass

export function generateMapFeatures(seed, terrain) {
  const {
    elevation,
    moisture,
    slope,
    sea,
    river,
    floodplain,
    plain,
    agriculture,
    ridgeField,
    valleyField,
    basinField,
    coastalLowland,
    flowAccum,
    arcSpineField,
    branchRidgeField,
    depositionalLowland,
    alluvialFanField,
    deltaField,
    portSuitability,
    crossingSuitability,
    passSuitability,
    prefectureMask,
    prefectureRegionId,
    naturalBarrierScore,
  } = terrain;

  const geography = terrain.geography || {};
  const geoHabitability = geography.habitability || null;
  const geoAccessibility = geography.accessibility || null;
  const geoNaturalCentrality = geography.naturalCentrality || geography.centrality || null;
  const geoLowlandCapacity = geography.lowlandCapacity || null;
  const geoValleyAccess = geography.valleyAccess || null;
  const geoCoastalAccess = geography.coastalAccess || null;
  const geoBarrier = geography.geographicBarrier || naturalBarrierScore || null;
  const geoCorridorSuitability = geography.corridorSuitability || null;

  function fieldValue(field, i, fallback = 0) {
    const v = field?.[i];
    return Number.isFinite(v) ? v : fallback;
  }

  function regionIdAt(x, y) {
    if (!inside(x, y)) return -1;
    const i = indexOf(x, y);
    if (sea[i]) return -1;
    if (prefectureMask?.[i]) return 0;
    if (!prefectureRegionId) return 0;
    const id = prefectureRegionId?.[i];
    return id !== undefined && id >= 0 ? id : -1;
  }

  function inFocusedPrefecture(p) {
    return Boolean(p && inside(p.x, p.y) && prefectureMask[indexOf(p.x, p.y)] && !sea[indexOf(p.x, p.y)]);
  }

  function localConfluenceScore(x, y) {
    let arms = 0;
    let strong = 0;
    for (const [dx, dy] of [[1,0],[-1,0],[0,1],[0,-1],[1,1],[-1,1],[1,-1],[-1,-1]]) {
      const nx = x + dx;
      const ny = y + dy;
      if (!inside(nx, ny)) continue;
      const rv = river[indexOf(nx, ny)];
      if (rv > 0.18) arms++;
      if (rv > 0.34) strong++;
    }
    return clamp((arms >= 3 ? 0.22 : arms === 2 ? 0.09 : 0) + strong * 0.04);
  }

  // --- 1. Human context: one full raster pass -----------------------------
  const developable = new Float32Array(SIZE);
  const ruralSuitability = new Float32Array(SIZE);
  const townSuitability = new Float32Array(SIZE);
  const valleySettlement = new Float32Array(SIZE);
  const coastalSettlement = new Float32Array(SIZE);
  const confluenceField = new Float32Array(SIZE);
  const barrierCost = new Float32Array(SIZE);
  const corridorCost = new Float32Array(SIZE);
  const settlementCluster = new Float32Array(SIZE);
  const settlementScore = new Float32Array(SIZE);

  for (let y = 0; y < MAP_H; y++) {
    for (let x = 0; x < MAP_W; x++) {
      const i = indexOf(x, y);
      if (sea[i]) {
        barrierCost[i] = INF;
        corridorCost[i] = INF;
        continue;
      }
      const naturalBarrier = naturalBarrierScore?.[i] || 0;
      const depositional = (depositionalLowland?.[i] || 0) + (alluvialFanField?.[i] || 0) * 0.62 + (deltaField?.[i] || 0) * 0.90;
      const highPenalty = Math.max(0, elevation[i] - 0.56);
      const lowSlope = clamp(1 - slope[i] * 2.3);
      const confluence = x > 0 && y > 0 && x < MAP_W - 1 && y < MAP_H - 1 ? localConfluenceScore(x, y) : 0;
      const openPlainPotential = clamp(plain[i] * 0.68 + agriculture[i] * 0.54 + basinField[i] * 0.30 + lowSlope * 0.22 - river[i] * 0.18 - valleyField[i] * 0.08 - ridgeField[i] * 0.22 - slope[i] * 0.26);
      const spine = (arcSpineField?.[i] || 0) * 0.58 + (branchRidgeField?.[i] || 0) * 0.38;
      confluenceField[i] = confluence;

      const geoH = fieldValue(geoHabitability, i, 0);
      const geoLow = fieldValue(geoLowlandCapacity, i, 0);
      const geoValley = fieldValue(geoValleyAccess, i, 0);
      const geoCoast = fieldValue(geoCoastalAccess, i, 0);
      const geoB = fieldValue(geoBarrier, i, naturalBarrier);
      const localDevelopable = clamp(
        plain[i] * 0.34 +
        agriculture[i] * 0.24 +
        basinField[i] * 0.24 +
        valleyField[i] * 0.24 +
        coastalLowland[i] * 0.18 +
        depositional * 0.22 +
        lowSlope * 0.10 -
        slope[i] * 0.82 -
        ridgeField[i] * 0.52 -
        spine * 0.24 -
        highPenalty * 1.14 -
        floodplain[i] * 0.03
      );
      developable[i] = clamp(localDevelopable * 0.68 + geoH * 0.34 + geoLow * 0.16 - geoB * 0.05);
      valleySettlement[i] = clamp((
        valleyField[i] * 0.52 +
        river[i] * 0.08 +
        confluence * 0.38 +
        depositional * 0.20 +
        basinField[i] * 0.16 +
        plain[i] * 0.08 +
        lowSlope * 0.12 -
        slope[i] * 0.54 -
        ridgeField[i] * 0.30 -
        spine * 0.16 -
        highPenalty * 0.70 -
        floodplain[i] * 0.10
      ) * 0.74 + geoValley * 0.30 + geoH * 0.08 - geoB * 0.04);
      coastalSettlement[i] = clamp((
        coastalLowland[i] * 0.50 +
        (portSuitability?.[i] || 0) * 0.30 +
        (deltaField?.[i] || 0) * 0.20 +
        plain[i] * 0.10 -
        slope[i] * 0.52 -
        ridgeField[i] * 0.24 -
        spine * 0.12
      ) * 0.76 + geoCoast * 0.32 + geoH * 0.06 - geoB * 0.04);
      const clusterNoise = 0.72 + fbm(x * 0.34 + 13, y * 0.34 - 31, seed + 7001) * 0.46 + valueNoise(x, y, seed + 7002, 8) * 0.16;
      settlementCluster[i] = clamp((developable[i] * 0.42 + valleySettlement[i] * 0.12 + coastalSettlement[i] * 0.22 + agriculture[i] * 0.48 + plain[i] * 0.32 + openPlainPotential * 0.46) * clusterNoise);
      ruralSuitability[i] = clamp(
        agriculture[i] * 0.54 +
        developable[i] * 0.30 +
        valleySettlement[i] * 0.18 +
        coastalSettlement[i] * 0.20 +
        openPlainPotential * 0.34 +
        settlementCluster[i] * 0.30 -
        Math.max(0, elevation[i] - 0.64) * 0.56
      );
      townSuitability[i] = clamp(
        developable[i] * 0.38 +
        agriculture[i] * 0.18 +
        valleySettlement[i] * 0.16 +
        coastalSettlement[i] * 0.30 +
        confluence * 0.20 +
        basinField[i] * 0.18 +
        plain[i] * 0.26 +
        openPlainPotential * 0.44 +
        settlementCluster[i] * 0.22 -
        slope[i] * 0.34 -
        ridgeField[i] * 0.17 -
        spine * 0.10
      );
      settlementScore[i] = clamp(
        ruralSuitability[i] * 0.48 +
        townSuitability[i] * 0.30 +
        confluence * 0.08 +
        fieldValue(geoHabitability, i, developable[i]) * 0.18 +
        fieldValue(geoNaturalCentrality, i, 0) * 0.12 -
        fieldValue(geoBarrier, i, 0) * 0.06
      );
      barrierCost[i] = 1 + slope[i] * 6.4 + ridgeField[i] * 3.2 + spine * 2.2 + highPenalty * 5.8 + river[i] * 0.25 + naturalBarrier * 1.2 - valleyField[i] * 0.55 - plain[i] * 0.30 - coastalLowland[i] * 0.14;
      corridorCost[i] = Math.max(0.25, barrierCost[i] - developable[i] * 0.42 - valleySettlement[i] * 0.36 - coastalSettlement[i] * 0.12 + hash2(x, y, seed + 7011) * 0.05);
    }
  }

  // --- region statistics ---------------------------------------------------
  const regionStats = new Map();
  function ensureRegion(regionId) {
    let st = regionStats.get(regionId);
    if (!st) {
      st = {
        id: regionId,
        area: 0,
        developableCells: 0,
        developableSum: 0,
        valleyCells: 0,
        coastCells: 0,
        townCells: 0,
        plainCells: 0,
        highCentralityCells: 0,
        habitabilitySum: 0,
        accessibilitySum: 0,
        centralitySum: 0,
        lowlandCapacitySum: 0,
        minX: MAP_W,
        minY: MAP_H,
        maxX: 0,
        maxY: 0,
      };
      regionStats.set(regionId, st);
    }
    return st;
  }
  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 regionId = regionIdAt(x, y);
      if (regionId < 0) continue;
      const st = ensureRegion(regionId);
      st.area++;
      const gHabit = fieldValue(geoHabitability, i, developable[i]);
      const gAccess = fieldValue(geoAccessibility, i, 0);
      const gCentral = fieldValue(geoNaturalCentrality, i, townSuitability[i]);
      const gLow = fieldValue(geoLowlandCapacity, i, plain[i]);
      st.developableSum += developable[i];
      st.habitabilitySum += gHabit;
      st.accessibilitySum += gAccess;
      st.centralitySum += gCentral;
      st.lowlandCapacitySum += gLow;
      if (developable[i] > 0.16 || gHabit > 0.24) st.developableCells++;
      if (valleySettlement[i] > 0.24 || fieldValue(geoValleyAccess, i, 0) > 0.25) st.valleyCells++;
      if (coastalSettlement[i] > 0.25 || fieldValue(geoCoastalAccess, i, 0) > 0.24) st.coastCells++;
      if (townSuitability[i] > 0.28 || gCentral > 0.31) st.townCells++;
      if (plain[i] > 0.24 || gLow > 0.26) st.plainCells++;
      if (gCentral > 0.36 && gHabit > 0.18) st.highCentralityCells++;
      st.minX = Math.min(st.minX, x);
      st.minY = Math.min(st.minY, y);
      st.maxX = Math.max(st.maxX, x);
      st.maxY = Math.max(st.maxY, y);
    }
  }

  function visibilityFactor(regionId, st) {
    if (!st || st.area <= 0) return 0;
    // Treat the focused prefecture and neighboring prefectures with the same
    // density curve.  Only genuinely clipped map-edge slivers are downscaled.
    return clamp(Math.sqrt(st.area / 1900), 0.32, 1.05);
  }

  function pickRegionalPoints(scoreArray, {
    stride = 1,
    threshold = 0.25,
    minDistance = 6,
    totalMax = 100,
    seedOffset = 0,
    quotaForRegion,
    predicate = () => true,
    kind = "Point",
    extraScore = () => 0,
  }) {
    const byRegion = new Map();
    for (let y = 2; y < MAP_H - 2; y += stride) {
      for (let x = 2; x < MAP_W - 2; x += stride) {
        const i = indexOf(x, y);
        if (sea[i] || !predicate(x, y, i)) continue;
        const regionId = regionIdAt(x, y);
        if (regionId < 0) continue;
        const score = scoreArray[i] + extraScore(x, y, i) + hash2(x, y, seed + seedOffset) * 0.055;
        if (score < threshold) continue;
        if (!byRegion.has(regionId)) byRegion.set(regionId, []);
        byRegion.get(regionId).push({ x, y, score, kind, regionId });
      }
    }
    const out = [];
    for (const [regionId, candidates] of [...byRegion.entries()].sort((a, b) => a[0] - b[0])) {
      const st = regionStats.get(regionId);
      const quota = quotaForRegion ? quotaForRegion(regionId, st) : 0;
      if (quota <= 0) continue;
      out.push(...pickEntities(candidates, {
        max: quota,
        minDistance,
        threshold,
        seed: seed + seedOffset + regionId * 1009,
        jitter: 0.04,
      }));
    }
    return out.sort((a, b) => b.score - a.score).slice(0, totalMax);
  }

  function pickGlobalPoints(scoreArray, { threshold, max, minDistance, seedOffset = 0, predicate = () => true, stride = 1 }) {
    const candidates = [];
    for (let y = 2; y < MAP_H - 2; y += stride) {
      for (let x = 2; x < MAP_W - 2; x += stride) {
        const i = indexOf(x, y);
        if (sea[i] || !predicate(x, y, i)) continue;
        const score = scoreArray[i] + hash2(x, y, seed + seedOffset) * 0.07;
        if (score >= threshold) candidates.push({ x, y, score, regionId: regionIdAt(x, y) });
      }
    }
    return pickEntities(candidates, { max, minDistance, threshold, seed: seed + seedOffset });
  }

  // --- 2. Sparse points ----------------------------------------------------
  let ports = pickGlobalPoints(portSuitability || coastalSettlement, {
    threshold: 0.30 + rand(seed, 1001) * 0.08,
    max: 10,
    minDistance: 13,
    seedOffset: 1000,
    predicate: (x, y, i) => coastalSettlement[i] > 0.14 || (portSuitability?.[i] || 0) > 0.25,
  }).map((p, n) => {
    const i = indexOf(p.x, p.y);
    const harborPotential = (portSuitability?.[i] || 0) + coastalLowland[i] * 0.22 + (deltaField?.[i] || 0) * 0.08 - slope[i] * 0.18;
    const portClass = n === 0 ? "major" : n < 3 && harborPotential > 0.34 ? "regional" : harborPotential > 0.24 ? "fishing" : "lake";
    const kind = portClass === "major" ? "Major Port" : portClass === "regional" ? "Regional Port" : portClass === "lake" ? "Lake Port" : "Fishing Port";
    return { ...p, harborPotential, portClass, kind, score: harborPotential };
  }).sort((a, b) => b.harborPotential - a.harborPotential);
  if (ports.length && !ports.some((p) => p.portClass === "major")) {
    ports[0].portClass = "major";
    ports[0].kind = "Major Port";
  }
  const commercialPorts = ports.filter((p) => p.portClass === "major" || p.portClass === "regional");

  const crossings = pickGlobalPoints(crossingSuitability || confluenceField, {
    threshold: 0.30 + rand(seed, 1011) * 0.06,
    max: 18,
    minDistance: 9,
    seedOffset: 1010,
    predicate: (x, y, i) => river[i] > 0.12 || confluenceField[i] > 0.09,
  }).map((p) => ({ ...p, kind: "River Crossing" }));

  const passes = pickGlobalPoints(passSuitability || valleySettlement, {
    threshold: 0.18 + rand(seed, 1021) * 0.06,
    max: 12,
    minDistance: 11,
    seedOffset: 1020,
    predicate: (x, y, i) => elevation[i] > 0.42 && !sea[i],
  }).map((p) => ({ ...p, kind: "Pass" }));

  // Phase 2: provisional upper-tier settlement anchors are selected directly
  // from the unified geography fields before lower-tier villages and market
  // towns are placed.  These anchors are not rendered as separate settlements;
  // they guide city selection and lower-tier spacing.
  const geographicUrbanAnchorScore = new Float32Array(SIZE);
  for (let i = 0; i < SIZE; i++) {
    if (sea[i]) continue;
    const gHabit = fieldValue(geoHabitability, i, developable[i]);
    const gAccess = fieldValue(geoAccessibility, i, 0);
    const gCentral = fieldValue(geoNaturalCentrality, i, townSuitability[i]);
    const gLow = fieldValue(geoLowlandCapacity, i, plain[i]);
    const gValley = fieldValue(geoValleyAccess, i, valleySettlement[i]);
    const gCoast = fieldValue(geoCoastalAccess, i, coastalSettlement[i]);
    const gBarrier = fieldValue(geoBarrier, i, naturalBarrierScore?.[i] || 0);
    geographicUrbanAnchorScore[i] = clamp(
      gCentral * 0.58 +
      gHabit * 0.34 +
      gAccess * 0.28 +
      gLow * 0.20 +
      gValley * 0.08 +
      gCoast * 0.10 +
      (portSuitability?.[i] || 0) * 0.12 +
      (crossingSuitability?.[i] || 0) * 0.08 +
      confluenceField[i] * 0.06 -
      gBarrier * 0.26 -
      slope[i] * 0.10
    );
  }
  const geographicUrbanAnchors = pickRegionalPoints(geographicUrbanAnchorScore, {
    stride: 2,
    threshold: 0.33 + rand(seed, 1026) * 0.025,
    totalMax: 18,
    minDistance: 24,
    seedOffset: 1025,
    kind: "Geographic Urban Anchor",
    predicate: (x, y, i) => fieldValue(geoHabitability, i, developable[i]) > 0.15 && fieldValue(geoBarrier, i, 0) < 0.58 && slope[i] < 0.42,
    quotaForRegion: (regionId, st) => {
      if (!st || st.developableCells < 40) return 0;
      const vf = visibilityFactor(regionId, st);
      const centralCells = st.highCentralityCells || 0;
      const raw = (centralCells / 180 + st.developableCells / 980 + 0.85) * vf;
      const min = st.area > 1800 || st.developableCells > 260 ? 1 : 0;
      const max = st.area > 4200 ? 4 : st.area > 2400 ? 3 : st.area > 900 ? 2 : 1;
      return Math.round(clamp(raw + rand(seed, 1027 + regionId * 17) * 0.7, min, max));
    },
    extraScore: (x, y, i) => fieldValue(geoNaturalCentrality, i, 0) * 0.18 + fieldValue(geoAccessibility, i, 0) * 0.10,
  }).map((p) => ({ ...p, candidateKind: "geographicAnchor", anchorScore: p.score }));
  const geographicAnchorInfluence = influenceFromPoints(geographicUrbanAnchors, 20, (p) => clamp((p.anchorScore || p.score || 0.4) * 1.35, 0.45, 1.35));

  const villageScore = new Float32Array(SIZE);
  for (let i = 0; i < SIZE; i++) {
    if (sea[i]) continue;
    villageScore[i] = clamp(
      ruralSuitability[i] * 0.42 +
      agriculture[i] * 0.42 +
      plain[i] * 0.34 +
      fieldValue(geoHabitability, i, developable[i]) * 0.22 +
      fieldValue(geoLowlandCapacity, i, plain[i]) * 0.18 +
      Math.max(0, plain[i] * 1.12 + agriculture[i] * 0.78 + basinField[i] * 0.30 - river[i] * 0.36 - valleyField[i] * 0.14 - flowAccum[i] * 0.10) * 0.48 +
      valleySettlement[i] * 0.06 +
      coastalSettlement[i] * 0.28 +
      settlementCluster[i] * 0.18 -
      geographicAnchorInfluence[i] * 0.08 -
      fieldValue(geoBarrier, i, 0) * 0.10 -
      river[i] * 0.10 -
      flowAccum[i] * 0.04
    );
  }
  let villages = pickRegionalPoints(villageScore, {
    stride: 2,
    threshold: 0.18 + rand(seed, 1031) * 0.030,
    totalMax: 190,
    minDistance: 6,
    seedOffset: 1030,
    kind: "Village",
    quotaForRegion: (regionId, st) => {
      if (!st || st.developableCells < 10) return 0;
      const vf = visibilityFactor(regionId, st);
      const raw = (st.developableCells / 40 + st.plainCells / 58 + st.valleyCells / 48 + st.coastCells / 46 + 2.1) * vf;
      const min = st.area > 2600 ? 12 : st.area > 1400 ? 7 : st.area > 520 ? 3 : st.area > 220 ? 1 : 0;
      const max = st.area > 3600 ? 46 : st.area > 2200 ? 32 : st.area > 900 ? 16 : 7;
      return Math.round(clamp(raw + rand(seed, 1033 + regionId * 19) * 1.5, min, max));
    },
  }).map((p, n) => {
    const i = indexOf(p.x, p.y);
    const kind = coastalSettlement[i] > 0.36 ? "Coastal Village" : valleySettlement[i] > 0.46 ? "Valley Village" : "Village";
    const population = Math.round((900 + Math.pow(rand(seed, 18000 + n * 17 + p.x * 3 + p.y), 1.30) * 9800 + ruralSuitability[i] * 4700 + agriculture[i] * 3600) / 100) * 100;
    return { ...p, kind, population };
  });

  // Supplemental open-plain villages: broad Japanese-style farmland should not be empty
  // just because it lacks a river/confluence anchor.
  const openPlainVillageScore = new Float32Array(SIZE);
  for (let i = 0; i < SIZE; i++) {
    if (sea[i]) continue;
    const open = Math.max(0, plain[i] * 0.92 + agriculture[i] * 0.72 + basinField[i] * 0.26 + (depositionalLowland?.[i] || 0) * 0.20 - river[i] * 0.22 - valleyField[i] * 0.10 - flowAccum[i] * 0.08 - slope[i] * 0.24 - ridgeField[i] * 0.14);
    openPlainVillageScore[i] = clamp(open * 0.82 + settlementCluster[i] * 0.14 + ruralSuitability[i] * 0.16 + fieldValue(geoHabitability, i, developable[i]) * 0.16 + fieldValue(geoLowlandCapacity, i, plain[i]) * 0.12 - geographicAnchorInfluence[i] * 0.05);
  }
  const supplementalPlainVillages = pickRegionalPoints(openPlainVillageScore, {
    stride: 2,
    threshold: 0.235 + rand(seed, 1036) * 0.020,
    totalMax: 60,
    minDistance: 7,
    seedOffset: 1035,
    kind: "Plain Village",
    predicate: (x, y, i) => plain[i] > 0.20 && agriculture[i] > 0.16 && river[i] < 0.30 && valleyField[i] < 0.52 && slope[i] < 0.32,
    quotaForRegion: (regionId, st) => {
      if (!st || st.plainCells < 24) return 0;
      const vf = visibilityFactor(regionId, st);
      const raw = (st.plainCells / 92 + st.developableCells / 260 + 0.9) * vf;
      const min = st.plainCells > 360 ? 3 : st.plainCells > 160 ? 1 : st.plainCells > 90 ? 1 : 0;
      const max = st.plainCells > 720 ? 12 : st.plainCells > 360 ? 8 : st.plainCells > 140 ? 4 : 2;
      return Math.round(clamp(raw + rand(seed, 1037 + regionId * 29) * 1.4, min, max));
    },
    extraScore: (x, y, i) => Math.max(0, plain[i] * 0.34 + agriculture[i] * 0.26 - river[i] * 0.20 - valleyField[i] * 0.12),
  }).filter((p) => distanceToNearest(villages, p.x, p.y) >= 6.5)
    .map((p, n) => {
      const i = indexOf(p.x, p.y);
      const population = Math.round((1100 + Math.pow(rand(seed, 18220 + n * 31 + p.x * 7 + p.y), 1.12) * 7600 + agriculture[i] * 3900 + plain[i] * 2200) / 100) * 100;
      return { ...p, kind: "Plain Village", population };
    });
  villages = [...villages, ...supplementalPlainVillages];

  let villageInfluence = influenceFromPoints(villages, 6, (v) => clamp((v.population || 1800) / 4200, 0.35, 1.2));

  const marketScore = new Float32Array(SIZE);
  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 openPlainMarket = Math.max(0, plain[i] * 0.68 + agriculture[i] * 0.52 + basinField[i] * 0.24 + (depositionalLowland?.[i] || 0) * 0.18 - river[i] * 0.16 - valleyField[i] * 0.06 - slope[i] * 0.18);
      const featurePull = Math.max(
        distanceToNearest(ports, x, y) < 10 ? 0.16 : 0,
        distanceToNearest(crossings, x, y) < 6 ? 0.035 : 0,
        confluenceField[i] * 0.08
      );
      const valleyMouth = valleyField[i] > 0.22 && (plain[i] > 0.22 || basinField[i] > 0.18 || coastalLowland[i] > 0.18) ? 0.14 : 0;
      marketScore[i] = clamp(
        townSuitability[i] * 0.38 +
        fieldValue(geoNaturalCentrality, i, townSuitability[i]) * 0.36 +
        fieldValue(geoAccessibility, i, 0) * 0.18 +
        fieldValue(geoHabitability, i, developable[i]) * 0.18 +
        agriculture[i] * 0.24 +
        plain[i] * 0.22 +
        openPlainMarket * 0.58 +
        coastalSettlement[i] * 0.18 +
        villageInfluence[i] * 0.24 +
        geographicAnchorInfluence[i] * 0.08 +
        featurePull +
        valleyMouth +
        basinField[i] * 0.12 +
        plain[i] * 0.16 +
        coastalLowland[i] * 0.07 -
        fieldValue(geoBarrier, i, 0) * 0.12 -
        slope[i] * 0.16 -
        ridgeField[i] * 0.07 -
        river[i] * 0.08 -
        flowAccum[i] * 0.035
      );
    }
  }

  let markets = pickRegionalPoints(marketScore, {
    stride: 2,
    threshold: 0.245 + rand(seed, 1041) * 0.035,
    totalMax: 70,
    minDistance: 9,
    seedOffset: 1040,
    kind: "Market Town",
    quotaForRegion: (regionId, st) => {
      if (!st || st.townCells < 8) return 0;
      const vf = visibilityFactor(regionId, st);
      const raw = (st.developableCells / 132 + st.plainCells / 148 + st.valleyCells / 128 + st.coastCells / 104 + 1.8) * vf;
      const min = st.area > 2600 ? 5 : st.area > 1200 ? 3 : st.area > 520 ? 1 : 0;
      const max = st.area > 3600 ? 20 : st.area > 2200 ? 14 : st.area > 800 ? 7 : 3;
      return Math.round(clamp(raw + rand(seed, 1043 + regionId * 23) * 0.8, min, max));
    },
    extraScore: (x, y, i) => (distanceToNearest(commercialPorts, x, y) < 10 ? 0.12 : 0) + coastalSettlement[i] * 0.08 + Math.max(0, plain[i] * 0.32 + agriculture[i] * 0.20 - river[i] * 0.16) * 0.09 + confluenceField[i] * 0.035,
  }).map((p, n) => {
    const i = indexOf(p.x, p.y);
    const kind = coastalSettlement[i] > 0.38 && distanceToNearest(ports, p.x, p.y) < 11 ? "Port Town" : valleySettlement[i] > 0.48 ? "Valley Market Town" : "Market Town";
    const population = Math.round((9000 + Math.pow(rand(seed, 18100 + n * 19 + p.x * 5 + p.y), 1.02) * 70000 + marketScore[i] * 40000 + villageInfluence[i] * 7800 + Math.max(0, plain[i] * 0.48 + agriculture[i] * 0.30 + basinField[i] * 0.18 - river[i] * 0.14) * 22000 + coastalSettlement[i] * 12000) / 1000) * 1000;
    return { ...p, kind, population };
  });

  const openPlainMarketScore = new Float32Array(SIZE);
  for (let i = 0; i < SIZE; i++) {
    if (sea[i]) continue;
    const open = Math.max(0, plain[i] * 0.86 + agriculture[i] * 0.64 + basinField[i] * 0.28 + (depositionalLowland?.[i] || 0) * 0.22 - river[i] * 0.20 - valleyField[i] * 0.10 - flowAccum[i] * 0.08 - slope[i] * 0.24 - ridgeField[i] * 0.16);
    openPlainMarketScore[i] = clamp(open * 0.78 + townSuitability[i] * 0.16 + fieldValue(geoNaturalCentrality, i, townSuitability[i]) * 0.18 + fieldValue(geoHabitability, i, developable[i]) * 0.12 + villageInfluence[i] * 0.16 + settlementCluster[i] * 0.10 + geographicAnchorInfluence[i] * 0.05);
  }
  const supplementalPlainMarkets = pickRegionalPoints(openPlainMarketScore, {
    stride: 2,
    threshold: 0.335 + rand(seed, 1046) * 0.025,
    totalMax: 26,
    minDistance: 12,
    seedOffset: 1045,
    kind: "Plain Market Town",
    predicate: (x, y, i) => plain[i] > 0.22 && agriculture[i] > 0.18 && river[i] < 0.28 && valleyField[i] < 0.50 && slope[i] < 0.30,
    quotaForRegion: (regionId, st) => {
      if (!st || st.plainCells < 60) return 0;
      const vf = visibilityFactor(regionId, st);
      const raw = (st.plainCells / 380 + st.developableCells / 720 + 0.25) * vf;
      const min = st.plainCells > 520 ? 1 : st.plainCells > 260 ? 1 : 0;
      const max = st.plainCells > 900 ? 4 : st.plainCells > 420 ? 3 : st.plainCells > 160 ? 1 : 1;
      return Math.round(clamp(raw + rand(seed, 1047 + regionId * 31) * 0.9, min, max));
    },
    extraScore: (x, y, i) => Math.max(0, plain[i] * 0.24 + agriculture[i] * 0.18 - river[i] * 0.14 - valleyField[i] * 0.08),
  }).filter((p) => distanceToNearest(markets, p.x, p.y) >= 10.5 && distanceToNearest(villages, p.x, p.y) >= 4.5)
    .map((p, n) => {
      const i = indexOf(p.x, p.y);
      const population = Math.round((10000 + Math.pow(rand(seed, 18340 + n * 37 + p.x * 11 + p.y), 1.02) * 52000 + openPlainMarketScore[i] * 26000 + agriculture[i] * 9000 + plain[i] * 8000) / 1000) * 1000;
      return { ...p, kind: "Plain Market Town", population };
    });
  markets = [...markets, ...supplementalPlainMarkets];

  // Sparse-area towns: when a developable basin/plain/coast has few nearby towns,
  // add a small market town candidate. This avoids large inhabited regions being
  // empty while still keeping minimum spacing from existing settlements.
  const existingTownInfluenceForSparseFill = influenceFromPoints([...markets, ...commercialPorts], 18, (p) => clamp((p.population || 9000) / 32000, 0.35, 1.25));
  const existingSettlementInfluenceForSparseFill = influenceFromPoints([...markets, ...villages, ...ports], 12, (p) => clamp((p.population || 1800) / 16000, 0.18, 1.0));
  const sparseTownScore = new Float32Array(SIZE);
  for (let i = 0; i < SIZE; i++) {
    if (sea[i]) continue;
    const remoteness = clamp((0.30 - existingTownInfluenceForSparseFill[i]) / 0.30);
    const settlementGap = clamp((0.42 - existingSettlementInfluenceForSparseFill[i]) / 0.42);
    const livable = clamp(townSuitability[i] * 0.34 + fieldValue(geoNaturalCentrality, i, townSuitability[i]) * 0.22 + fieldValue(geoHabitability, i, developable[i]) * 0.20 + fieldValue(geoAccessibility, i, 0) * 0.12 + developable[i] * 0.22 + plain[i] * 0.24 + agriculture[i] * 0.18 + basinField[i] * 0.16 + coastalSettlement[i] * 0.14 + valleySettlement[i] * 0.12 - fieldValue(geoBarrier, i, 0) * 0.12 - slope[i] * 0.22 - ridgeField[i] * 0.10);
    sparseTownScore[i] = clamp(livable * (0.42 + remoteness * 0.88) + settlementGap * 0.16);
  }
  const sparseMarkets = pickRegionalPoints(sparseTownScore, {
    stride: 2,
    threshold: 0.315 + rand(seed, 1049) * 0.025,
    totalMax: 24,
    minDistance: 15,
    seedOffset: 1048,
    kind: "Sparse Market Town",
    predicate: (x, y, i) => existingTownInfluenceForSparseFill[i] < 0.34 && developable[i] > 0.12 && slope[i] < 0.36 && ridgeField[i] < 0.55,
    quotaForRegion: (regionId, st) => {
      if (!st || st.developableCells < 90) return 0;
      const vf = visibilityFactor(regionId, st);
      const underServed = clamp(1.0 - ((markets.filter((m) => m.regionId === regionId).length || 0) / Math.max(1, st.area / 850)));
      const raw = (st.developableCells / 900 + st.plainCells / 720 + st.coastCells / 560 + 0.55) * vf * (0.55 + underServed * 0.75);
      return Math.round(clamp(raw + rand(seed, 1050 + regionId * 37) * 0.45, 0, st.area > 2000 ? 2 : 1));
    },
    extraScore: (x, y, i) => clamp((0.34 - existingTownInfluenceForSparseFill[i]) * 0.38 + plain[i] * 0.10 + agriculture[i] * 0.08 + coastalSettlement[i] * 0.06),
  })
    .filter((p) => distanceToNearest(markets, p.x, p.y) >= 12 && distanceToNearest(villages, p.x, p.y) >= 4.5)
    .map((p, n) => {
      const i = indexOf(p.x, p.y);
      const population = Math.round((8000 + Math.pow(rand(seed, 18480 + n * 41 + p.x * 13 + p.y), 1.08) * 36000 + sparseTownScore[i] * 26000) / 1000) * 1000;
      return { ...p, kind: "Sparse Market Town", population };
    });
  markets = [...markets, ...sparseMarkets];

  const defenseScore = new Float32Array(SIZE);
  for (let i = 0; i < SIZE; i++) {
    if (sea[i]) continue;
    defenseScore[i] = clamp(
      confluenceField[i] * 0.38 +
      townSuitability[i] * 0.16 +
      ridgeField[i] * clamp(1 - Math.abs(elevation[i] - 0.52) / 0.25) * 0.42 +
      plain[i] * 0.08 -
      floodplain[i] * 0.36 -
      coastalLowland[i] * 0.08
    );
  }
  const castles = pickGlobalPoints(defenseScore, {
    threshold: 0.34 + rand(seed, 1051) * 0.06,
    max: 5,
    minDistance: 16,
    seedOffset: 1050,
  }).map((p) => ({
    ...p,
    kind: elevation[indexOf(p.x, p.y)] > 0.55 ? "Mountain Castle" : elevation[indexOf(p.x, p.y)] > 0.38 ? "Hilltop Castle" : "Flatland Castle",
  }));

  const castleTowns = castles.map((c, n) => {
    const near = markets.slice().sort((a, b) => Math.hypot(a.x - c.x, a.y - c.y) - Math.hypot(b.x - c.x, b.y - c.y))[0];
    const x = near && Math.hypot(near.x - c.x, near.y - c.y) < 10 ? near.x : c.x;
    const y = near && Math.hypot(near.x - c.x, near.y - c.y) < 10 ? near.y : c.y;
    return { x, y, score: c.score, kind: "Castle Town", population: 12000 + Math.round(rand(seed, 1060 + n) * 22000 / 1000) * 1000, regionId: regionIdAt(x, y) };
  });

  // --- 3. Cities by region, without detailed urban flood-fill --------------
  function estimateUrbanCapacity(p, radius = 22, densityBias = 1.0) {
    if (!p || !inside(p.x, p.y)) return 0;
    const centerRegion = regionIdAt(p.x, p.y);
    let capacity = 0;
    const r = Math.ceil(radius);
    for (let dy = -r; dy <= r; dy++) {
      for (let dx = -r; dx <= r; dx++) {
        const x = p.x + dx;
        const y = p.y + dy;
        if (!inside(x, y)) continue;
        const i = indexOf(x, y);
        if (sea[i]) continue;
        if (centerRegion >= 0 && regionIdAt(x, y) !== centerRegion) continue;
        const d = Math.hypot(dx, dy);
        if (d > radius) continue;
        const dev = clamp(developable[i] * 0.66 + fieldValue(geoHabitability, i, developable[i]) * 0.34);
        if (dev < 0.04) continue;
        const radial = clamp(1 - d / Math.max(1, radius));
        const terrainMultiplier = clamp(0.58 + plain[i] * 0.40 + agriculture[i] * 0.18 + basinField[i] * 0.24 + coastalLowland[i] * 0.16 + valleyField[i] * 0.06 + fieldValue(geoLowlandCapacity, i, plain[i]) * 0.24 + fieldValue(geoAccessibility, i, 0) * 0.12 - fieldValue(geoBarrier, i, 0) * 0.18 - slope[i] * 0.34 - ridgeField[i] * 0.14, 0.24, 1.46);
        capacity += dev * (900 + 6200 * Math.pow(radial, 1.25)) * terrainMultiplier * densityBias;
      }
    }
    return Math.max(26000, Math.round(capacity / 1000) * 1000);
  }

  const urbanCandidates = [
    ...geographicUrbanAnchors.map((p) => ({ ...p, candidateKind: "geographicAnchor" })),
    ...markets.map((p) => ({ ...p, candidateKind: "town" })),
    ...castleTowns.map((p) => ({ ...p, candidateKind: "castleTown" })),
    ...commercialPorts.map((p) => ({ ...p, candidateKind: "port" })),
    ...crossings.filter((p) => confluenceField[indexOf(p.x, p.y)] > 0.12).map((p) => ({ ...p, candidateKind: "crossing" })),
  ];

  const cityCandidateByRegion = new Map();
  for (const p of urbanCandidates) {
    const i = indexOf(p.x, p.y);
    const regionId = regionIdAt(p.x, p.y);
    if (regionId < 0) continue;
    const st = regionStats.get(regionId);
    const cityRadius = st && st.area > 2400 ? 30 : st && st.area > 900 ? 26 : 22;
    const capacity = estimateUrbanCapacity(p, cityRadius, p.candidateKind === "town" ? 1.14 : p.candidateKind === "port" ? 1.10 : 1.06);
    const score =
      Math.log10(capacity + 1) * 0.66 +
      townSuitability[i] * 0.86 +
      developable[i] * 0.68 +
      fieldValue(geoNaturalCentrality, i, townSuitability[i]) * 1.34 +
      fieldValue(geoHabitability, i, developable[i]) * 0.58 +
      fieldValue(geoAccessibility, i, 0) * 0.44 +
      confluenceField[i] * 0.18 +
      (p.candidateKind === "geographicAnchor" ? 0.42 : 0) +
      (p.candidateKind === "port" ? 0.48 : 0) +
      (p.candidateKind === "castleTown" ? 0.22 : 0) -
      fieldValue(geoBarrier, i, 0) * 0.36 +
      hash2(p.x, p.y, seed + 12000) * 0.16;
    if (!cityCandidateByRegion.has(regionId)) cityCandidateByRegion.set(regionId, []);
    cityCandidateByRegion.get(regionId).push({ ...p, score, capacity, regionId });
  }

  const modernCities = [];
  const usedCitySites = [];
  for (const [regionId, list] of [...cityCandidateByRegion.entries()].sort((a, b) => a[0] - b[0])) {
    const st = regionStats.get(regionId);
    if (!st || st.developableCells < 30) continue;
    const vf = visibilityFactor(regionId, st);
    const maxCities = clamp(
      Math.round((st.developableCells / 680 + (st.highCentralityCells || 0) / 360 + 0.95) * vf + rand(seed, 12100 + regionId * 17) * 1.1),
      (st.area > 900 || st.developableCells > 180) ? 1 : 0,
      st.area > 3600 ? 6 : st.area > 2200 ? 4 : st.area > 900 ? 3 : 1
    );
    const selected = pickEntities(list, {
      max: maxCities,
      minDistance: 24,
      threshold: 0,
      seed: seed + 12110 + regionId * 313,
      jitter: 0.02,
    });
    for (const p of selected) {
      const minD = p.capacity >= 420000 ? 34 : p.capacity >= 220000 ? 28 : p.capacity >= 110000 ? 23 : 20;
      if (usedCitySites.some((q) => {
        const qMinD = q.capacity >= 420000 ? 34 : q.capacity >= 220000 ? 28 : q.capacity >= 110000 ? 23 : 20;
        return Math.hypot(q.x - p.x, q.y - p.y) < Math.max(minD, qMinD) * 0.82;
      })) continue;
      usedCitySites.push(p);
      modernCities.push(p);
    }
  }

  // No focused-prefecture fallback: all prefecture regions use the same city
  // selection rules, so the highlighted region is not overwritten after the
  // regional pass.

  modernCities.sort((a, b) => b.capacity - a.capacity || b.score - a.score);
  const regionalCapitalSlots = Math.max(1, Math.min(4, Math.round(Math.sqrt(Math.max(1, modernCities.length)))));
  for (const [rank, city] of modernCities.entries()) {
    const i = indexOf(city.x, city.y);
    const st = regionStats.get(city.regionId);
    const isFirstInRegion = !modernCities.slice(0, rank).some((c) => c.regionId === city.regionId);
    const isPrefecturalCapital = inFocusedPrefecture(city) && !modernCities.slice(0, rank).some((c) => c.isPrefecturalCapital);
    const geoTierScore = clamp(
      fieldValue(geoNaturalCentrality, i, 0) * 0.48 +
      fieldValue(geoHabitability, i, 0) * 0.24 +
      fieldValue(geoAccessibility, i, 0) * 0.18 +
      Math.log10((city.capacity || 26000) + 1) / 7 * 0.26
    );
    const isTopCenter = rank < regionalCapitalSlots || geoTierScore > 0.58;
    const isRegionalCapital = isFirstInRegion && (isTopCenter || (city.capacity || 0) > 210000 || (st?.highCentralityCells || 0) > 220);
    const rawPop = isRegionalCapital
      ? 180000 + rand(seed, 12201 + city.regionId * 17) * (isTopCenter ? 760000 : 360000)
      : 52000 + Math.pow(rand(seed, 12202 + rank * 19 + city.x), 0.68) * 360000;
    const capMultiplier = isRegionalCapital ? (isTopCenter ? 1.66 : 1.42) : 1.20;
    const population = Math.round(Math.min(rawPop, city.capacity * capMultiplier) / 1000) * 1000;
    const floor = isRegionalCapital ? (isTopCenter ? 210000 : 120000) : 42000;
    city.population = Math.max(floor, population);
    city.isPrefecturalCapital = isPrefecturalCapital;
    city.isRegionalCapital = isRegionalCapital;
    city.geographicTierScore = Math.round(geoTierScore * 1000) / 1000;
    city.rank = isPrefecturalCapital ? "Prefectural Capital" : isRegionalCapital ? "Regional Capital" : city.population >= 200000 ? "Regional City" : "Local City";
    city.kind = city.rank;
    city.urbanRadius = clamp(5.8 + Math.sqrt(city.population) / 72, 8, isRegionalCapital ? 34 : 24);
    city.coreRadius = clamp(1.8 + Math.sqrt(city.population) / 380, 2.2, isRegionalCapital ? 6.5 : 5.6);
    city.sprawlRadius = clamp(city.urbanRadius * (isRegionalCapital ? 1.45 : city.population >= 120000 ? 1.28 : 1.15), city.urbanRadius + 2, isRegionalCapital ? 44 : 30);
    city.urbanWeight = clamp(0.95 + Math.log10(Math.max(10000, city.population)) * 0.29, 1.08, 2.5);
  }

  function urbanSettlementExclusionRadius(city, tier = "market") {
    const pop = city?.population || 0;
    const base = pop >= 500000 ? 9.5 : pop >= 240000 ? 7.6 : pop >= 110000 ? 6.0 : 4.8;
    return tier === "village" ? base + 2.0 : base;
  }
  const marketsBeforeHierarchyFilter = markets.length;
  const villagesBeforeHierarchyFilter = villages.length;
  markets = markets.filter((m) => {
    const nearestCity = modernCities.reduce((best, city) => {
      const d = Math.hypot(m.x - city.x, m.y - city.y);
      return d < best.d ? { city, d } : best;
    }, { city: null, d: Infinity });
    if (!nearestCity.city) return true;
    return nearestCity.d >= urbanSettlementExclusionRadius(nearestCity.city, "market");
  });
  villages = villages.filter((v) => {
    const nearestCity = modernCities.reduce((best, city) => {
      const d = Math.hypot(v.x - city.x, v.y - city.y);
      return d < best.d ? { city, d } : best;
    }, { city: null, d: Infinity });
    if (nearestCity.city && nearestCity.d < urbanSettlementExclusionRadius(nearestCity.city, "village")) return false;
    return distanceToNearest(markets, v.x, v.y) >= 3.4;
  });
  villageInfluence = influenceFromPoints(villages, 6, (v) => clamp((v.population || 1800) / 4200, 0.35, 1.2));
  const settlementHierarchyDebug = {
    version: "phase2-unified-settlement-hierarchy",
    geographicUrbanAnchors: geographicUrbanAnchors.length,
    marketsBeforeHierarchyFilter,
    marketsAfterHierarchyFilter: markets.length,
    villagesBeforeHierarchyFilter,
    villagesAfterHierarchyFilter: villages.length,
    regionalCapitalSlots,
  };

  function cityPopulationCap(city) {
    const radius = city?.isRegionalCapital ? 29 : city?.population >= 350000 ? 26 : 18;
    const bias = city?.isRegionalCapital ? 1.12 : 1.0;
    return estimateUrbanCapacity(city, radius, bias);
  }

  // --- 4. Field-derived transport corridors -------------------------------
  const preliminaryUrbanInfluence = influenceFromPoints(modernCities, 18, (c) => clamp((c.population || 60000) / 260000, 0.55, 2.0));
  const preliminaryTownInfluence = influenceFromPoints([...markets, ...commercialPorts], 10, (p) => p.portClass === "major" ? 1.35 : clamp((p.population || 12000) / 36000, 0.42, 1.1));
  const preliminaryVillageInfluence = influenceFromPoints(villages, 7, (v) => clamp((v.population || 1800) / 5200, 0.22, 0.9));
  const settlementDemand = new Float32Array(SIZE);
  const urbanEdge = new Float32Array(SIZE);
  const logisticsPreSuitability = new Float32Array(SIZE);

  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 density = clamp(preliminaryUrbanInfluence[i] * 0.62 + preliminaryTownInfluence[i] * 0.34 + preliminaryVillageInfluence[i] * 0.16);
      settlementDemand[i] = density;
      urbanEdge[i] = clamp(1 - Math.abs(density - 0.46) / 0.32);
      logisticsPreSuitability[i] = clamp(
        agriculture[i] * 0.30 +
        plain[i] * 0.24 +
        basinField[i] * 0.14 +
        coastalLowland[i] * 0.12 +
        preliminaryTownInfluence[i] * 0.18 +
        urbanEdge[i] * 0.34 -
        preliminaryUrbanInfluence[i] * 0.20 -
        slope[i] * 0.50 -
        ridgeField[i] * 0.32
      );
    }
  }

  function buildTransportCostFields() {
    const expressway = new Float32Array(SIZE);
    const rail = new Float32Array(SIZE);
    const national = new Float32Array(SIZE);
    const local = new Float32Array(SIZE);
    const expresswayPotential = new Float32Array(SIZE);
    const railPotential = new Float32Array(SIZE);
    const nationalPotential = new Float32Array(SIZE);
    const localPotential = new Float32Array(SIZE);

    function seaAdjacency(x, y, radius = 1) {
      let sum = 0;
      let total = 0;
      for (let dy = -radius; dy <= radius; dy++) {
        for (let dx = -radius; dx <= radius; dx++) {
          if (!dx && !dy) continue;
          const nx = x + dx;
          const ny = y + dy;
          if (!inside(nx, ny)) continue;
          total++;
          if (sea[indexOf(nx, ny)]) sum += 1;
        }
      }
      return total > 0 ? sum / total : 0;
    }

    function highAltitudeTransportClosed(i) {
      // Above this contour the generator should treat mountains as no-road
      // terrain.  A strong mapped pass is the exception, so genuine saddle
      // crossings can still exist without roads drilling through entire ranges.
      return elevation[i] >= 0.70;
    }

    for (let y = 0; y < MAP_H; y++) {
      for (let x = 0; x < MAP_W; x++) {
        const i = indexOf(x, y);
        if (sea[i] || highAltitudeTransportClosed(i)) {
          expressway[i] = rail[i] = national[i] = local[i] = INF;
          expresswayPotential[i] = railPotential[i] = nationalPotential[i] = localPotential[i] = 0;
          continue;
        }
        const density = settlementDemand[i];
        const mediumDensity = clamp(1 - Math.abs(density - 0.42) / 0.30);
        const highDensity = clamp((density - 0.32) / 0.50);
        const lowland = clamp(plain[i] * 0.48 + basinField[i] * 0.28 + valleyField[i] * 0.24 + coastalLowland[i] * 0.26 + agriculture[i] * 0.16);
        const pass = passSuitability?.[i] || 0;
        const crossing = crossingSuitability?.[i] || 0;
        const waterCrossingPenalty = river[i] > 0.18 ? (1 - crossing) * (0.72 + river[i] * 1.35) : 0;
        const seaNear = seaAdjacency(x, y, 1);
        const seaBroad = seaAdjacency(x, y, 3);
        const seaWide = seaAdjacency(x, y, 5);
        // Roads should use coastal lowlands when there is a settlement/port reason,
        // but should not casually trace beaches or hop over small bays.
        const coastalTraversePenalty = clamp(seaBroad * 1.72 + seaWide * 0.82 - coastalLowland[i] * 0.48 - (portSuitability?.[i] || 0) * 0.30);
        const highMountain = clamp((elevation[i] - 0.52) * 3.6 + slope[i] * 0.95 + ridgeField[i] * 1.05 - pass * 0.55 - valleyField[i] * 0.12);
        const extremeMountain = clamp((elevation[i] - 0.64) * 4.8 + slope[i] * 1.55 + ridgeField[i] * 1.45 - pass * 0.80);
        const denseCorePenalty = clamp((density - 0.66) / 0.28);
        const openPlainParallelPenalty = clamp(plain[i] * 0.48 + agriculture[i] * 0.26 - density * 0.20 - valleyField[i] * 0.20 - coastalLowland[i] * 0.12);
        const boundaryRidgePenalty = Math.pow(naturalBarrierScore?.[i] || 0, 2);

        if (extremeMountain > 0.92 && pass < 0.34) {
          expressway[i] = rail[i] = INF;
          national[i] = elevation[i] >= 0.68 ? INF : 2.9 + extremeMountain * 2.8 + waterCrossingPenalty;
          local[i] = elevation[i] >= 0.69 ? INF : 2.2 + extremeMountain * 2.3 + waterCrossingPenalty * 0.55;
          expresswayPotential[i] = 0;
          railPotential[i] = 0;
          nationalPotential[i] = clamp(lowland * 0.18 + pass * 0.24 - extremeMountain * 0.55);
          localPotential[i] = clamp(valleySettlement[i] * 0.14 + pass * 0.18 - extremeMountain * 0.28);
          continue;
        }

        expresswayPotential[i] = clamp(
          mediumDensity * 0.62 +
          urbanEdge[i] * 0.38 +
          logisticsPreSuitability[i] * 0.54 +
          lowland * 0.36 +
          agriculture[i] * 0.16 -
          denseCorePenalty * 0.54 -
          slope[i] * 0.82 -
          highMountain * 0.92 -
          coastalTraversePenalty * 0.32 -
          river[i] * 0.14
        );
        railPotential[i] = clamp(
          highDensity * 0.80 +
          preliminaryTownInfluence[i] * 0.22 +
          lowland * 0.46 +
          valleyField[i] * 0.22 +
          coastalLowland[i] * 0.22 -
          slope[i] * 1.28 -
          highMountain * 1.10 -
          coastalTraversePenalty * 0.26 -
          ridgeField[i] * 0.34
        );
        nationalPotential[i] = clamp(
          density * 0.46 +
          preliminaryTownInfluence[i] * 0.32 +
          preliminaryVillageInfluence[i] * 0.20 +
          agriculture[i] * 0.24 +
          valleyField[i] * 0.28 +
          coastalLowland[i] * 0.26 +
          pass * 0.18 +
          crossing * 0.18 -
          slope[i] * 0.48 -
          highMountain * 0.24 -
          coastalTraversePenalty * 0.16 -
          ridgeField[i] * 0.18
        );
        localPotential[i] = clamp(
          preliminaryVillageInfluence[i] * 0.52 +
          agriculture[i] * 0.38 +
          coastalSettlement[i] * 0.30 +
          valleySettlement[i] * 0.30 +
          developable[i] * 0.18 -
          slope[i] * 0.34 -
          coastalTraversePenalty * 0.08 -
          ridgeField[i] * 0.10
        );

        expressway[i] = Math.max(0.18,
          1.62 - expresswayPotential[i] * 0.96 +
          denseCorePenalty * 1.30 +
          slope[i] * 5.4 +
          highMountain * 5.8 +
          extremeMountain * 4.2 +
          boundaryRidgePenalty * 4.2 +
          waterCrossingPenalty * 2.1 +
          coastalTraversePenalty * 2.65 +
          seaNear * 2.35 +
          seaWide * 1.10 +
          openPlainParallelPenalty * 0.12 +
          hash2(x, y, seed + 13301) * 0.04
        );
        rail[i] = Math.max(0.16,
          1.48 - railPotential[i] * 1.02 +
          slope[i] * 7.2 +
          highMountain * 7.0 +
          extremeMountain * 4.8 +
          boundaryRidgePenalty * 2.4 +
          waterCrossingPenalty * 1.7 +
          coastalTraversePenalty * 1.75 +
          seaNear * 1.50 +
          seaWide * 0.70 +
          hash2(x, y, seed + 13302) * 0.03
        );
        national[i] = Math.max(0.16,
          1.28 - nationalPotential[i] * 0.84 +
          slope[i] * 2.8 +
          ridgeField[i] * 1.18 +
          Math.max(0, elevation[i] - 0.62) * 3.0 +
          highMountain * 2.9 +
          boundaryRidgePenalty * 1.8 +
          waterCrossingPenalty * 1.25 -
          valleyField[i] * 0.18 -
          coastalLowland[i] * 0.08 +
          coastalTraversePenalty * 1.55 +
          seaNear * 0.84 +
          seaWide * 0.48 -
          pass * 0.42 +
          hash2(x, y, seed + 13303) * 0.05
        );
        local[i] = Math.max(0.14,
          1.12 - localPotential[i] * 0.86 +
          slope[i] * 1.72 +
          ridgeField[i] * 0.82 +
          Math.max(0, elevation[i] - 0.68) * 1.9 +
          highMountain * 1.24 +
          boundaryRidgePenalty * 0.72 +
          waterCrossingPenalty * 0.65 -
          valleyField[i] * 0.22 -
          coastalLowland[i] * 0.10 +
          coastalTraversePenalty * 0.82 +
          seaNear * 0.48 +
          seaWide * 0.26 +
          hash2(x, y, seed + 13304) * 0.07
        );
      }
    }
    return { expressway, rail, national, local, expresswayPotential, railPotential, nationalPotential, localPotential };
  }

  const transportFields = buildTransportCostFields();
  const cachedInfluenceFromPaths = createPathInfluenceCache(influenceFromPaths);

  const componentCityInfluence = influenceFromPoints([...modernCities, ...markets], 11, (p) => clamp((p.population || 8000) / 50000, 0.18, 8.0));
  const componentCapitalInfluence = influenceFromPoints(modernCities.filter((p) => p.isPrefecturalCapital), 16, () => 5.0);

  const corridorSkeleton = new Uint8Array(SIZE);
  for (let i = 0; i < SIZE; i++) {
    if (sea[i]) continue;
    corridorSkeleton[i] = Math.round(clamp(
      valleyField[i] * 0.34 +
      coastalLowland[i] * 0.24 +
      plain[i] * 0.16 +
      basinField[i] * 0.16 +
      logisticsPreSuitability[i] * 0.18 +
      (passSuitability?.[i] || 0) * 0.18 +
      (crossingSuitability?.[i] || 0) * 0.10 -
      ridgeField[i] * 0.22 -
      slope[i] * 0.28
    ) * 255);
  }

  function chooseCorridorSeeds(potentialField, spacing, maxCount, threshold, predicate = () => true, seedOffset = 0, mode = "national") {
    const candidates = [];
    for (let y = 3; y < MAP_H - 3; y += 2) {
      for (let x = 3; x < MAP_W - 3; x += 2) {
        const i = indexOf(x, y);
        if (sea[i] || !predicate(x, y, i)) continue;
        const skeletonWeight = mode === "rail" ? 0.24 : mode === "expressway" ? 0.18 : 0.28;
        const score = potentialField[i] + (corridorSkeleton[i] / 255) * skeletonWeight + hash2(x, y, seed + seedOffset) * 0.055;
        if (score >= threshold) candidates.push({ x, y, score, regionId: regionIdAt(x, y) });
      }
    }
    return pickEntities(candidates, { max: maxCount, minDistance: spacing, threshold, seed: seed + seedOffset, jitter: 0.03 });
  }

  function corridorAllowance(i) {
    return clamp(settlementDemand[i] * 0.58 + valleyField[i] * 0.42 + coastalLowland[i] * 0.36 - plain[i] * 0.18 - agriculture[i] * 0.14);
  }

  function endpointSupport(i) {
    if (i < 0 || i >= SIZE || sea[i]) return 0;
    return clamp(
      settlementDemand[i] * 0.66 +
      preliminaryTownInfluence[i] * 0.36 +
      preliminaryVillageInfluence[i] * 0.24 +
      valleySettlement[i] * 0.22 +
      coastalSettlement[i] * 0.18 +
      agriculture[i] * 0.10 +
      (passSuitability?.[i] || 0) * 0.24 +
      (crossingSuitability?.[i] || 0) * 0.16 -
      slope[i] * 0.22 -
      ridgeField[i] * 0.14
    );
  }

  function nearMapEdgePoint(x, y, margin = 4) {
    return x <= margin || y <= margin || x >= MAP_W - 1 - margin || y >= MAP_H - 1 - margin;
  }

  function naturalEndpoint(x, y, i, mode = "national") {
    if (nearMapEdgePoint(x, y, 4)) return true;
    const support = endpointSupport(i);
    if (mode === "expressway") return support > 0.24 && settlementDemand[i] > 0.06;
    if (mode === "rail") return support > 0.18 && settlementDemand[i] > 0.07;
    if (mode === "local") return support > 0.12 || valleySettlement[i] > 0.14 || coastalSettlement[i] > 0.14;
    return support > 0.14 || transportFields.nationalPotential[i] > 0.37;
  }

  function traceCorridorByCost(start, goalRegionPredicate, costField, penaltyField, options = {}) {
    if (!start || !inside(start.x, start.y)) return [];
    const startIndex = indexOf(start.x, start.y);
    if (sea[startIndex] || costField[startIndex] >= INF) return [];
    const score = new Float32Array(SIZE);
    const cameFrom = new Int32Array(SIZE);
    const closed = new Uint8Array(SIZE);
    score.fill(INF);
    cameFrom.fill(-1);
    const heap = new MinHeap();
    score[startIndex] = 0;
    heap.push({ i: startIndex, f: 0 });
    const curvePenalty = options.curvePenalty ?? 0.12;
    const penaltyStrength = options.penaltyStrength ?? 1.0;
    const sameRegion = options.regionId ?? regionIdAt(start.x, start.y);
    const minGoalDistance = options.minGoalDistance ?? 18;
    const maxExpanded = options.maxExpanded ?? SIZE * 2;
    const bounds = options.bounds || null;
    const goalHint = options.goalHint || null;
    const heuristicWeight = options.heuristicWeight ?? 0;
    let goalIndex = -1;
    let expanded = 0;

    while (heap.length && expanded++ < maxExpanded) {
      const current = heap.pop();
      if (!current || closed[current.i]) continue;
      closed[current.i] = 1;
      const [cx, cy] = xyOf(current.i);
      if (current.i !== startIndex && Math.hypot(cx - start.x, cy - start.y) >= minGoalDistance && goalRegionPredicate(cx, cy, current.i)) {
        goalIndex = current.i;
        break;
      }
      for (let dy = -1; dy <= 1; dy++) {
        for (let dx = -1; dx <= 1; dx++) {
          if (!dx && !dy) continue;
          const nx = cx + dx;
          const ny = cy + dy;
          if (!inside(nx, ny)) continue;
          if (bounds && (nx < bounds.minX || nx > bounds.maxX || ny < bounds.minY || ny > bounds.maxY)) continue;
          const ni = indexOf(nx, ny);
          if (closed[ni] || sea[ni] || costField[ni] >= INF) continue;
          if (sameRegion >= 0 && options.keepRegion !== false && regionIdAt(nx, ny) !== sameRegion) continue;
          const prev = cameFrom[current.i];
          let turn = 0;
          if (prev >= 0) {
            const [px, py] = xyOf(prev);
            const ax = cx - px;
            const ay = cy - py;
            turn = Math.abs(ax * dy - ay * dx) > 0 ? curvePenalty : 0;
          }
          const existing = penaltyField?.[ni] || 0;
          const antiConcentration = existing * penaltyStrength * (1 - corridorAllowance(ni) * 0.72);
          const terrainFlowBias = (options.terrainFlowBias ?? 0) * clamp(
            valleyField[ni] * 0.54 +
            coastalLowland[ni] * 0.28 +
            plain[ni] * 0.16 +
            (passSuitability?.[ni] || 0) * 0.34 -
            ridgeField[ni] * 0.24 -
            slope[ni] * 0.22
          );
          const surfaceGrain = (options.surfaceGrain ?? 0) * valueNoise(nx, ny, seed + 13941, 18);
          const nd = score[current.i] + Math.max(0.08, costField[ni] + antiConcentration + turn - terrainFlowBias + surfaceGrain) * Math.hypot(dx, dy);
          if (nd < score[ni]) {
            score[ni] = nd;
            cameFrom[ni] = current.i;
            const h = goalHint ? Math.hypot(nx - goalHint.x, ny - goalHint.y) * heuristicWeight : 0;
            heap.push({ i: ni, f: nd + h });
          }
        }
      }
    }
    if (goalIndex < 0) return [];
    const path = [];
    for (let p = goalIndex; p >= 0; p = cameFrom[p]) {
      path.push(xyOf(p));
      if (p === startIndex) break;
    }
    return path.reverse();
  }

  function addCorridorInfluencePenalty(penaltyField, corridor, radius = 7, strength = 0.35) {
    for (const [x, y] of corridor || []) {
      for (let dy = -radius; dy <= radius; dy++) {
        for (let dx = -radius; dx <= radius; dx++) {
          const nx = x + dx;
          const ny = y + dy;
          if (!inside(nx, ny)) continue;
          const d = Math.hypot(dx, dy);
          if (d > radius) continue;
          const i = indexOf(nx, ny);
          if (sea[i]) continue;
          const openPlain = clamp(plain[i] * 0.54 + agriculture[i] * 0.34 - settlementDemand[i] * 0.24 - valleyField[i] * 0.22 - coastalLowland[i] * 0.18);
          const allowParallel = corridorAllowance(i);
          penaltyField[i] = Math.max(penaltyField[i], strength * (1 - d / radius) * (0.48 + openPlain * 1.15 - allowParallel * 0.42));
        }
      }
    }
  }

  function densifyPathByCost(path, costField) {
    if (!path || path.length < 2) return path || [];
    const out = [];
    for (let k = 0; k < path.length - 1; k++) {
      const [x0, y0] = path[k];
      const [x1, y1] = path[k + 1];
      const steps = Math.max(1, Math.ceil(Math.hypot(x1 - x0, y1 - y0)));
      for (let s = 0; s <= steps; s++) {
        if (k > 0 && s === 0) continue;
        const t = s / steps;
        const x = Math.round(x0 + (x1 - x0) * t);
        const y = Math.round(y0 + (y1 - y0) * t);
        if (!inside(x, y)) return [];
        const i = indexOf(x, y);
        if (sea[i] || costField?.[i] >= INF) return [];
        if (!out.length || out[out.length - 1][0] !== x || out[out.length - 1][1] !== y) out.push([x, y]);
      }
    }
    return out.length >= 2 ? out : [];
  }

  function relaxRouteToTerrain(path, costField, options = {}) {
    if (!path || path.length < 5) return path || [];
    const radius = options.radius ?? 2;
    const iterations = options.iterations ?? 1;
    const lineWeight = options.lineWeight ?? 0.42;
    const costWeight = options.costWeight ?? 1.0;
    const grain = options.grain ?? 0.05;
    let out = path.map((p) => [p[0], p[1]]);
    for (let iter = 0; iter < iterations; iter++) {
      const src = out.map((p) => [p[0], p[1]]);
      for (let k = 1; k < src.length - 1; k++) {
        const prev = out[k - 1];
        const cur = src[k];
        const next = src[k + 1];
        let best = cur;
        let bestScore = INF;
        const vx = next[0] - prev[0];
        const vy = next[1] - prev[1];
        const vLen2 = Math.max(1e-6, vx * vx + vy * vy);
        for (let dy = -radius; dy <= radius; dy++) {
          for (let dx = -radius; dx <= radius; dx++) {
            if (dx * dx + dy * dy > radius * radius) continue;
            const x = cur[0] + dx;
            const y = cur[1] + dy;
            if (!inside(x, y)) continue;
            const i = indexOf(x, y);
            if (sea[i] || costField[i] >= INF) continue;
            const t = clamp(((x - prev[0]) * vx + (y - prev[1]) * vy) / vLen2);
            const projX = prev[0] + vx * t;
            const projY = prev[1] + vy * t;
            const lineDist = Math.hypot(x - projX, y - projY);
            const neighborDist = Math.abs(Math.hypot(x - prev[0], y - prev[1]) - Math.hypot(cur[0] - prev[0], cur[1] - prev[1])) * 0.08;
            const terrainBonus = valleyField[i] * 0.34 + coastalLowland[i] * 0.16 + plain[i] * 0.08 + (passSuitability?.[i] || 0) * 0.18;
            const terrainPenalty = slope[i] * 0.42 + ridgeField[i] * 0.32 + Math.max(0, elevation[i] - 0.62) * 0.45;
            const localGrain = grain * valueNoise(x, y, seed + 15123 + k * 17, 14);
            const score = costField[i] * costWeight + lineDist * lineWeight + neighborDist + terrainPenalty - terrainBonus + localGrain;
            if (score < bestScore) {
              bestScore = score;
              best = [x, y];
            }
          }
        }
        out[k] = best;
      }
    }
    const deduped = [];
    let last = "";
    for (const p of out) {
      const key = `${p[0]},${p[1]}`;
      if (key !== last) {
        deduped.push(p);
        last = key;
      }
    }
    if (deduped.length < 2) return path;
    const dense = densifyPathByCost(deduped, costField);
    return dense.length >= 2 ? dense : path;
  }

  function endpointFromPath(path) {
    const p = path?.[path.length - 1];
    return p ? { x: p[0], y: p[1], regionId: regionIdAt(p[0], p[1]) } : null;
  }

  function qualityLimitsForMode(mode, overrides = {}) {
    const base = mode === "rail"
      ? { maxCompactness: 2.45, maxSteepShare: 0.18, maxHighElevationShare: 0, minAvgPotential: 0.12 }
      : mode === "expressway"
        ? { maxCompactness: 2.75, maxSteepShare: 0.24, maxHighElevationShare: 0, minAvgPotential: 0.10 }
        : mode === "local"
          ? { maxCompactness: 3.6, maxSteepShare: 0.46, maxHighElevationShare: 0.18, minAvgPotential: 0.02 }
          : { maxCompactness: 3.3, maxSteepShare: 0.38, maxHighElevationShare: 0.04, minAvgPotential: 0.04 };
    return { ...base, ...overrides };
  }

  function pathWaterCrossingStats(path) {
    let seaCells = 0;
    let maxSeaRun = 0;
    let currentSeaRun = 0;
    let sampled = 0;
    for (let k = 1; k < (path?.length || 0); k++) {
      const a = path[k - 1];
      const b = path[k];
      if (!a || !b) continue;
      const steps = Math.max(1, Math.ceil(Math.hypot(a[0] - b[0], a[1] - b[1])));
      for (let s = 0; s <= steps; s++) {
        const t = s / steps;
        const x = Math.round(a[0] + (b[0] - a[0]) * t);
        const y = Math.round(a[1] + (b[1] - a[1]) * t);
        sampled++;
        const isSea = !inside(x, y) || sea[indexOf(x, y)];
        if (isSea) {
          seaCells++;
          currentSeaRun++;
          maxSeaRun = Math.max(maxSeaRun, currentSeaRun);
        } else {
          currentSeaRun = 0;
        }
      }
    }
    return { seaCells, maxSeaRun, seaShare: sampled ? seaCells / sampled : 0 };
  }

  function pathTunnelStats(path) {
    let tunnelCells = 0;
    let maxTunnelRun = 0;
    let currentTunnelRun = 0;
    let sampled = 0;
    for (let k = 1; k < (path?.length || 0); k++) {
      const a = path[k - 1];
      const b = path[k];
      if (!a || !b) continue;
      const steps = Math.max(1, Math.ceil(Math.hypot(a[0] - b[0], a[1] - b[1])));
      for (let s = 0; s <= steps; s++) {
        const t = s / steps;
        const x = Math.round(a[0] + (b[0] - a[0]) * t);
        const y = Math.round(a[1] + (b[1] - a[1]) * t);
        if (!inside(x, y)) continue;
        sampled++;
        const i = indexOf(x, y);
        const isTunnel = !sea[i] && ((elevation[i] >= 0.72 && ridgeField[i] >= 0.34) || naturalBarrierScore[i] >= 0.72);
        if (isTunnel) {
          tunnelCells++;
          currentTunnelRun++;
          maxTunnelRun = Math.max(maxTunnelRun, currentTunnelRun);
        } else {
          currentTunnelRun = 0;
        }
      }
    }
    return { tunnelCells, maxTunnelRun, tunnelShare: sampled ? tunnelCells / sampled : 0 };
  }

  function routePhysicalAcceptable(path, mode, overrides = {}) {
    if (!path?.length) return false;
    const water = pathWaterCrossingStats(path);
    const tunnel = pathTunnelStats(path);
    const bridgeLimit = overrides.bridgeLimit ?? (mode === "expressway" ? 20 : 10);
    const tunnelLimit = overrides.tunnelLimit ?? (mode === "expressway" ? 10 : 0);
    const maxSeaRun = overrides.maxSeaRun ?? bridgeLimit;
    const maxTunnelRun = overrides.maxTunnelRun ?? tunnelLimit;
    const maxSeaShare = overrides.maxSeaShare ?? (mode === "expressway" ? 0.22 : mode === "rail" ? 0.030 : mode === "national" ? 0.10 : 0.05);
    const maxTunnelShare = overrides.maxTunnelShare ?? (mode === "expressway" ? 0.12 : 0);
    if (water.maxSeaRun > maxSeaRun || water.seaShare > maxSeaShare) return false;
    if (tunnel.maxTunnelRun > maxTunnelRun || tunnel.tunnelShare > maxTunnelShare) return false;
    if (water.seaCells > 0) {
      const first = path[0];
      const last = path[path.length - 1];
      const ai = inside(first?.[0], first?.[1]) ? indexOf(first[0], first[1]) : -1;
      const bi = inside(last?.[0], last?.[1]) ? indexOf(last[0], last[1]) : -1;
      const demand = clamp((ai >= 0 ? settlementDemand[ai] || 0 : 0) + (bi >= 0 ? settlementDemand[bi] || 0 : 0));
      const threshold = mode === "local" ? 0.62 : mode === "rail" ? 0.54 : 0.42;
      if (demand < threshold && mode !== "national" && mode !== "expressway") return false;
    }
    return true;
  }

  function transportRouteAcceptable(path, mode, potentialField, penaltyField = null, overrides = {}) {
    if (!routePhysicalAcceptable(path, mode, overrides)) return false;
    return routeQualityAcceptable(path, {
      sea,
      elevation,
      slope,
      potential: potentialField,
      penalty: penaltyField,
      highElevationThreshold: 0.70,
      steepThreshold: mode === "rail" ? 0.34 : mode === "expressway" ? 0.40 : 0.48,
    }, qualityLimitsForMode(mode, overrides));
  }

  function generateCorridorsFromField({ mode = "national", potentialField, costField, spacing, maxCount, threshold, minLength, penaltyRadius, penaltyStrength, curvePenalty, terrainFlowBias = 0, surfaceGrain = 0, relaxRadius = 2, relaxLineWeight = 0.40, seedOffset, startPredicate, goalPredicate }) {
    const paths = [];
    const penaltyField = new Float32Array(SIZE);
    const seeds = chooseCorridorSeeds(potentialField, spacing, maxCount * 2, threshold, startPredicate, seedOffset, mode);
    const usedEndpoints = [];
    for (const start of seeds) {
      if (paths.length >= maxCount) break;
      if (distanceToNearest(usedEndpoints, start.x, start.y) < spacing * 0.55) continue;
      let path = traceCorridorByCost(
        start,
        (x, y, i) => goalPredicate(start, x, y, i, usedEndpoints),
        costField,
        penaltyField,
        { curvePenalty, penaltyStrength: penaltyStrength * 2.2, minGoalDistance: minLength, regionId: start.regionId, terrainFlowBias, surfaceGrain }
      );
      if (path.length < minLength) continue;
      const rawPath = path;
      path = relaxRouteToTerrain(rawPath, costField, { radius: relaxRadius, lineWeight: relaxLineWeight, grain: surfaceGrain, iterations: 1 });
      if (path.length < Math.max(2, rawPath.length * 0.55)) path = rawPath;
      if (!transportRouteAcceptable(path, mode, potentialField, penaltyField, { minLength, maxLength: mode === "expressway" ? 130 : mode === "rail" ? 112 : 150 })) continue;
      paths.push(path);
      usedEndpoints.push(start);
      const end = endpointFromPath(path);
      if (end) usedEndpoints.push(end);
      addCorridorInfluencePenalty(penaltyField, path, penaltyRadius, penaltyStrength);
    }
    return paths;
  }

  function rasterizeNetworkComponents(paths, mode, potentialField) {
    const occupied = new Uint8Array(SIZE);
    for (const path of paths) {
      for (const [x, y] of path || []) {
        if (inside(x, y) && !sea[indexOf(x, y)]) occupied[indexOf(x, y)] = 1;
      }
    }
    const componentId = new Int32Array(SIZE);
    componentId.fill(-1);
    const components = [];
    for (let i = 0; i < SIZE; i++) {
      if (!occupied[i] || componentId[i] >= 0) continue;
      const id = components.length;
      const queue = [i];
      const cells = [];
      const boundary = [];
      componentId[i] = id;
      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);
        let edge = false;
        for (let dy = -1; dy <= 1; dy++) {
          for (let dx = -1; dx <= 1; dx++) {
            if (!dx && !dy) continue;
            const nx = x + dx;
            const ny = y + dy;
            if (!inside(nx, ny)) { edge = true; continue; }
            const ni = indexOf(nx, ny);
            if (!occupied[ni]) {
              edge = true;
              continue;
            }
            if (componentId[ni] < 0) {
              componentId[ni] = id;
              queue.push(ni);
            }
          }
        }
        if (edge) boundary.push({ x, y, i: cur });
      }
      let lengthScore = 0;
      let densityScore = 0;
      let townScore = 0;
      let logisticsScore = 0;
      let capitalScore = 0;
      let sx = 0;
      let sy = 0;
      const sampleStride = Math.max(1, Math.floor(cells.length / 80));
      for (let c = 0; c < cells.length; c += sampleStride) {
        const ci = cells[c];
        const x = ci % MAP_W;
        const y = Math.floor(ci / MAP_W);
        lengthScore += sampleStride;
        sx += x * sampleStride;
        sy += y * sampleStride;
        densityScore += settlementDemand[ci] * sampleStride;
        logisticsScore += logisticsPreSuitability[ci] * sampleStride;
        townScore += componentCityInfluence[ci] * sampleStride;
        capitalScore += componentCapitalInfluence[ci] * sampleStride;
      }
      const importance =
        Math.sqrt(lengthScore) * 1.20 +
        densityScore * 0.38 +
        townScore * 0.55 +
        logisticsScore * 0.24 +
        capitalScore;
      components.push({
        id,
        mode,
        cells,
        boundary,
        cx: sx / Math.max(1, lengthScore),
        cy: sy / Math.max(1, lengthScore),
        importance,
        length: cells.length,
        repairCount: 0,
        potential: cells.reduce((sum, ci) => sum + (potentialField?.[ci] || 0), 0) / Math.max(1, cells.length),
      });
    }
    return { occupied, componentId, components };
  }

  function componentAnchor(component, target, costField, usedAnchors = []) {
    if (!component?.boundary?.length) return null;
    let best = null;
    let bestScore = INF;
    const stride = Math.max(1, Math.floor(component.boundary.length / 90));
    for (let k = 0; k < component.boundary.length; k += stride) {
      const p = component.boundary[k];
      if (costField[p.i] >= INF) continue;
      if (distanceToNearest(usedAnchors, p.x, p.y) < 8) continue;
      const d = target ? Math.hypot(p.x - target.x, p.y - target.y) : 0;
      const score = d + costField[p.i] * 3.5 - corridorAllowance(p.i) * 2.4 + hash2(p.x, p.y, seed + 13701) * 1.8;
      if (score < bestScore) {
        bestScore = score;
        best = { x: p.x, y: p.y, componentId: component.id, regionId: regionIdAt(p.x, p.y) };
      }
    }
    return best;
  }

  function componentBounds(aComp, bComp, pad = 18) {
    let minX = MAP_W - 1;
    let minY = MAP_H - 1;
    let maxX = 0;
    let maxY = 0;
    for (const comp of [aComp, bComp]) {
      for (const p of comp.boundary || []) {
        minX = Math.min(minX, p.x);
        minY = Math.min(minY, p.y);
        maxX = Math.max(maxX, p.x);
        maxY = Math.max(maxY, p.y);
      }
    }
    return {
      minX: Math.max(0, minX - pad),
      minY: Math.max(0, minY - pad),
      maxX: Math.min(MAP_W - 1, maxX + pad),
      maxY: Math.min(MAP_H - 1, maxY + pad),
    };
  }

  function traceCoarseRepair(start, targetComp, raster, costField, penaltyField, options = {}) {
    const scale = options.coarseScale ?? 3;
    if (scale <= 1) return null;
    const cw = Math.ceil(MAP_W / scale);
    const ch = Math.ceil(MAP_H / scale);
    const cSize = cw * ch;
    const bounds = options.bounds || { minX: 0, minY: 0, maxX: MAP_W - 1, maxY: MAP_H - 1 };
    const cbounds = {
      minX: Math.max(0, Math.floor(bounds.minX / scale) - 1),
      minY: Math.max(0, Math.floor(bounds.minY / scale) - 1),
      maxX: Math.min(cw - 1, Math.ceil(bounds.maxX / scale) + 1),
      maxY: Math.min(ch - 1, Math.ceil(bounds.maxY / scale) + 1),
    };
    const cIndex = (x, y) => y * cw + x;
    const cCost = new Float32Array(cSize);
    const cTarget = new Uint8Array(cSize);
    cCost.fill(INF);
    for (let cy = cbounds.minY; cy <= cbounds.maxY; cy++) {
      for (let cx = cbounds.minX; cx <= cbounds.maxX; cx++) {
        let best = INF;
        let target = 0;
        for (let dy = 0; dy < scale; dy++) {
          for (let dx = 0; dx < scale; dx++) {
            const x = cx * scale + dx;
            const y = cy * scale + dy;
            if (!inside(x, y)) continue;
            const i = indexOf(x, y);
            if (sea[i] || costField[i] >= INF) continue;
            best = Math.min(best, costField[i] + (penaltyField?.[i] || 0) * (options.penaltyStrength ?? 1));
            if (raster.componentId[i] === targetComp.id) target = 1;
          }
        }
        const ci = cIndex(cx, cy);
        cCost[ci] = best;
        cTarget[ci] = target;
      }
    }
    const sx = Math.floor(start.x / scale);
    const sy = Math.floor(start.y / scale);
    if (sx < cbounds.minX || sx > cbounds.maxX || sy < cbounds.minY || sy > cbounds.maxY) return null;
    const startCi = cIndex(sx, sy);
    if (cCost[startCi] >= INF) return null;

    const score = new Float32Array(cSize);
    const cameFrom = new Int32Array(cSize);
    const closed = new Uint8Array(cSize);
    score.fill(INF);
    cameFrom.fill(-1);
    const heap = new MinHeap();
    score[startCi] = 0;
    heap.push({ i: startCi, f: 0 });
    let goal = -1;
    let guard = 0;
    while (heap.length && guard++ < cSize * 2) {
      const cur = heap.pop();
      if (!cur || closed[cur.i]) continue;
      closed[cur.i] = 1;
      if (cur.i !== startCi && cTarget[cur.i]) {
        goal = cur.i;
        break;
      }
      const cx = cur.i % cw;
      const cy = Math.floor(cur.i / cw);
      for (let dy = -1; dy <= 1; dy++) {
        for (let dx = -1; dx <= 1; dx++) {
          if (!dx && !dy) continue;
          const nx = cx + dx;
          const ny = cy + dy;
          if (nx < cbounds.minX || nx > cbounds.maxX || ny < cbounds.minY || ny > cbounds.maxY) continue;
          const ni = cIndex(nx, ny);
          if (closed[ni] || cCost[ni] >= INF) continue;
          const nd = score[cur.i] + cCost[ni] * Math.hypot(dx, dy);
          if (nd < score[ni]) {
            score[ni] = nd;
            cameFrom[ni] = cur.i;
            const h = Math.hypot(nx - targetComp.cx / scale, ny - targetComp.cy / scale) * (options.heuristicWeight ?? 0.18);
            heap.push({ i: ni, f: nd + h });
          }
        }
      }
    }
    if (goal < 0) return null;
    const coarse = [];
    for (let p = goal; p >= 0; p = cameFrom[p]) {
      const cx = p % cw;
      const cy = Math.floor(p / cw);
      coarse.push([Math.min(MAP_W - 1, Math.round(cx * scale + scale * 0.5)), Math.min(MAP_H - 1, Math.round(cy * scale + scale * 0.5))]);
      if (p === startCi) break;
    }
    coarse.reverse();
    if (coarse.length < 2) return null;
    const full = [[start.x, start.y]];
    for (let k = 1; k < coarse.length; k++) {
      const a = full[full.length - 1];
      const b = coarse[k];
      const steps = Math.max(1, Math.ceil(Math.hypot(a[0] - b[0], a[1] - b[1])));
      for (let s = 1; s <= steps; s++) {
        const t = s / steps;
        const x = Math.round(a[0] + (b[0] - a[0]) * t);
        const y = Math.round(a[1] + (b[1] - a[1]) * t);
        if (inside(x, y) && !sea[indexOf(x, y)] && costField[indexOf(x, y)] < INF) full.push([x, y]);
      }
    }
    const tail = full[full.length - 1];
    let nearest = null;
    let nearestD = INF;
    const stride = Math.max(1, Math.floor((targetComp.boundary?.length || 1) / 80));
    for (let k = 0; k < (targetComp.boundary?.length || 0); k += stride) {
      const p = targetComp.boundary[k];
      const d = Math.hypot(p.x - tail[0], p.y - tail[1]);
      if (d < nearestD) {
        nearestD = d;
        nearest = p;
      }
    }
    if (nearest && nearestD <= scale * 3 + 4) {
      const a = full[full.length - 1];
      const steps = Math.max(1, Math.ceil(nearestD));
      for (let s = 1; s <= steps; s++) {
        const t = s / steps;
        const x = Math.round(a[0] + (nearest.x - a[0]) * t);
        const y = Math.round(a[1] + (nearest.y - a[1]) * t);
        if (inside(x, y) && !sea[indexOf(x, y)] && costField[indexOf(x, y)] < INF) full.push([x, y]);
      }
    }
    return full;
  }

  function sampledBarrierBetween(a, b) {
    const steps = Math.max(1, Math.ceil(Math.hypot(a.cx - b.cx, a.cy - b.cy)));
    let sum = 0;
    let n = 0;
    let seaHits = 0;
    for (let s = 0; s <= steps; s++) {
      const t = s / steps;
      const x = Math.round(a.cx + (b.cx - a.cx) * t);
      const y = Math.round(a.cy + (b.cy - a.cy) * t);
      if (!inside(x, y)) continue;
      const i = indexOf(x, y);
      if (sea[i]) { seaHits++; sum += 1.25; n++; continue; }
      sum += clamp((naturalBarrierScore?.[i] || 0) * 0.75 + ridgeField[i] * 0.28 + slope[i] * 0.25 + Math.max(0, elevation[i] - 0.58) * 0.45 - (passSuitability?.[i] || 0) * 0.35);
      n++;
    }
    return n ? clamp(sum / n + seaHits / n) : 1;
  }

  function canRepairConnection(a, b, mode, d) {
    const demand = clamp((a.potential || 0) * 0.55 + (b.potential || 0) * 0.55 + Math.sqrt(Math.max(0, a.importance + b.importance)) / 11 - d / 180);
    const barrier = sampledBarrierBetween(a, b);
    if (mode === "rail" && demand < 0.45) return false;
    if (mode === "expressway" && demand < 0.55) return false;
    if (mode === "national" && demand < 0.24 && d > 48) return false;
    if (barrier > 0.66 && demand < 0.72) return false;
    if (d > 90 && demand < 0.75) return false;
    return true;
  }

  function repairTransportConnectivity(paths, mode, costField, potentialField, options = {}) {
    const debug = { mode, components: [], repairs: [] };
    const raster = rasterizeNetworkComponents(paths, mode, potentialField);
    debug.components = raster.components.map((c) => ({
      id: c.id,
      mode,
      importance: Math.round(c.importance * 10) / 10,
      length: c.length,
      potential: Math.round(c.potential * 100) / 100,
      cx: Math.round(c.cx),
      cy: Math.round(c.cy),
      cells: c.cells.filter((_, k) => k % Math.max(1, Math.floor(c.cells.length / 40)) === 0).slice(0, 40).map((i) => xyOf(i)),
    }));
    const important = raster.components
      .filter((c) => c.length >= (options.minComponentCells ?? 18) && c.importance >= (options.minImportance ?? 8))
      .sort((a, b) => b.importance - a.importance)
      .slice(0, options.maxComponents ?? 8);
    if (important.length < 2) return debug;

    const networkPenalty = cachedInfluenceFromPaths(paths, options.penaltyRadius ?? 8, `${mode}:repair`);
    const usedAnchors = [];
    const maxRepairs = options.maxRepairs ?? 4;
    for (let r = 0; r < maxRepairs; r++) {
      let bestPair = null;
      let bestScore = INF;
      for (let a = 0; a < important.length; a++) {
        for (let b = a + 1; b < important.length; b++) {
          const ca = important[a];
          const cb = important[b];
          if (ca.repairCount >= 2 || cb.repairCount >= 2) continue;
          const d = Math.hypot(ca.cx - cb.cx, ca.cy - cb.cy);
          if (d < (options.minRepairDistance ?? 14) || d > (options.maxRepairDistance ?? 120)) continue;
          if (!canRepairConnection(ca, cb, mode, d)) continue;
          const barrier = sampledBarrierBetween(ca, cb);
          const score = d / Math.sqrt(ca.importance + cb.importance) + barrier * 22 + (ca.repairCount + cb.repairCount) * 18;
          if (score < bestScore) {
            bestScore = score;
            bestPair = [ca, cb];
          }
        }
      }
      if (!bestPair) break;
      const [aComp, bComp] = bestPair;
      const roughTarget = bComp.boundary[Math.floor(bComp.boundary.length / 2)];
      const start = componentAnchor(aComp, roughTarget, costField, usedAnchors);
      const goalTarget = start ? componentAnchor(bComp, start, costField, usedAnchors) : null;
      if (!start || !goalTarget) break;
      const bounds = componentBounds(aComp, bComp, options.searchPad ?? 20);
      let path = traceCoarseRepair(start, bComp, raster, costField, networkPenalty, {
        ...options,
        bounds,
        heuristicWeight: 0.22,
      });
      if (!path) {
        path = traceCorridorByCost(
          start,
          (x, y, i) => raster.componentId[i] === bComp.id || (potentialField[i] > (options.highPotentialThreshold ?? 0.42) && Math.hypot(x - goalTarget.x, y - goalTarget.y) < 5),
          costField,
          networkPenalty,
          {
            curvePenalty: options.curvePenalty ?? 0.14,
            penaltyStrength: options.penaltyStrength ?? 1.8,
            terrainFlowBias: options.terrainFlowBias ?? 0.12,
            surfaceGrain: options.surfaceGrain ?? 0.012,
            minGoalDistance: Math.min(12, Math.max(5, Math.hypot(start.x - goalTarget.x, start.y - goalTarget.y) * 0.35)),
            keepRegion: false,
            maxExpanded: Math.floor(SIZE * 0.45),
            bounds,
            goalHint: goalTarget,
            heuristicWeight: 0.18,
          }
        );
      }
      if (path.length < (options.minAddedLength ?? 6) || path.length > (options.maxAddedLength ?? 120)) {
        aComp.repairCount++;
        continue;
      }
      const rawRepairPath = path;
      path = relaxRouteToTerrain(rawRepairPath, costField, { radius: options.relaxRadius ?? 2, lineWeight: options.relaxLineWeight ?? 0.36, grain: options.surfaceGrain ?? 0.012, iterations: 1 });
      if (path.length < Math.max(2, rawRepairPath.length * 0.55)) path = rawRepairPath;
      if (path.length < (options.minAddedLength ?? 6) || path.length > (options.maxAddedLength ?? 120)) {
        aComp.repairCount++;
        continue;
      }
      if (!transportRouteAcceptable(path, mode, potentialField, networkPenalty, { minLength: options.minAddedLength ?? 6, maxLength: options.maxAddedLength ?? 120 })) {
        aComp.repairCount++;
        continue;
      }
      paths.push(path);
      debug.repairs.push({ mode, path, from: aComp.id, to: bComp.id });
      usedAnchors.push(start, goalTarget);
      aComp.repairCount++;
      bComp.repairCount++;
      addCorridorInfluencePenalty(networkPenalty, path, options.penaltyRadius ?? 8, options.addedPenalty ?? 0.38);
    }
    return debug;
  }

  function routeLight(a, b, snapRadius = 3, costField = transportFields.local) {
    if (!a || !b) return [];
    const start = {
      x: Math.round(a.x),
      y: Math.round(a.y),
      regionId: regionIdAt(Math.round(a.x), Math.round(a.y)),
    };
    const target = { x: Math.round(b.x), y: Math.round(b.y) };
    if (!inside(start.x, start.y) || !inside(target.x, target.y)) return [];
    const dist = Math.hypot(start.x - target.x, start.y - target.y);
    const routed = traceCorridorByCost(
      start,
      (x, y) => Math.hypot(x - target.x, y - target.y) <= snapRadius,
      costField,
      null,
      {
        curvePenalty: 0.025,
        penaltyStrength: 0,
        minGoalDistance: Math.min(5, Math.max(2, dist * 0.10)),
        keepRegion: false,
        maxExpanded: Math.min(SIZE, Math.max(1800, Math.floor(dist * dist * 5.5))),
        terrainFlowBias: 0.16,
        surfaceGrain: 0.025,
      }
    );
    if (routed.length >= 2) return relaxRouteToTerrain(routed, costField, { radius: 2, lineWeight: 0.34, grain: 0.035, iterations: 1 });

    // Fallback for rare isolated cells: still use a snapped line, but keep the
    // radius small and terrain-weighted so it does not become a long artificial
    // chord across mountains.
    const steps = Math.max(2, Math.ceil(dist * 1.35));
    const out = [];
    let lastKey = "";
    for (let s = 0; s <= steps; s++) {
      const t = s / steps;
      const fx = start.x + (target.x - start.x) * t;
      const fy = start.y + (target.y - start.y) * t;
      let best = null;
      let bestCost = INF;
      const radius = Math.max(1, Math.min(snapRadius, 2));
      for (let dy = -radius; dy <= radius; dy++) {
        for (let dx = -radius; dx <= radius; dx++) {
          const x = Math.round(fx + dx);
          const y = Math.round(fy + dy);
          if (!inside(x, y)) continue;
          const i = indexOf(x, y);
          if (sea[i] || costField[i] >= INF) continue;
          const lineDist = Math.hypot(x - fx, y - fy);
          const cost = lineDist * 0.92 + costField[i] * 1.10 - valleyField[i] * 0.22 - coastalLowland[i] * 0.10 + ridgeField[i] * 0.22 + slope[i] * 0.18;
          if (cost < bestCost) {
            bestCost = cost;
            best = [x, y];
          }
        }
      }
      if (!best) continue;
      const key = `${best[0]},${best[1]}`;
      if (key !== lastKey) {
        out.push(best);
        lastKey = key;
      }
    }
    return relaxRouteToTerrain(out, costField, { radius: 2, lineWeight: 0.36, grain: 0.030, iterations: 1 });
  }

  function importantNodesForRegion(regionId) {
    const inRegion = (p) => regionIdAt(p.x, p.y) === regionId;
    return [
      ...modernCities.filter(inRegion).map((p) => ({ ...p, nodeWeight: 8 + (p.population || 0) / 120000 })),
      ...markets.filter(inRegion).map((p) => ({ ...p, nodeWeight: 3.2 + (p.population || 0) / 25000 })),
      ...commercialPorts.filter(inRegion).map((p) => ({ ...p, nodeWeight: p.portClass === "major" ? 6.5 : 4.6 })),
      ...passes.filter(inRegion).map((p) => ({ ...p, nodeWeight: 2.2 })),
    ].sort((a, b) => b.nodeWeight - a.nodeWeight).slice(0, (regionStats.get(regionId)?.area || 0) > 2200 ? 16 : 10);
  }

  function dedupePointCandidates(points, minDistance = 5) {
    const out = [];
    for (const raw of points) {
      if (!raw) continue;
      const x = Math.round(raw.x);
      const y = Math.round(raw.y);
      if (!inside(x, y)) continue;
      const i = indexOf(x, y);
      if (sea[i]) continue;
      if (out.some((p) => Math.hypot(p.x - x, p.y - y) < minDistance)) continue;
      out.push({ ...raw, x, y, regionId: regionIdAt(x, y), score: raw.score ?? raw.nodeWeight ?? 0.5 });
    }
    return out;
  }

  function modeSettlementAnchors(mode, potentialField, max = 48) {
    const anchors = [];
    const add = (p, baseScore, role) => {
      if (!p || !inside(p.x, p.y)) return;
      const x = Math.round(p.x);
      const y = Math.round(p.y);
      const i = indexOf(x, y);
      if (sea[i] || regionIdAt(x, y) < 0) return;
      anchors.push({
        x, y, role, regionId: regionIdAt(x, y),
        score: baseScore + (potentialField?.[i] || 0) * 1.4 + settlementDemand[i] * 0.75 + valleyField[i] * 0.18 + coastalLowland[i] * 0.12,
      });
    };

    for (const c of modernCities) add(c, mode === "expressway" ? 4.2 : mode === "rail" ? 4.0 : 3.5, "city");
    for (const pnt of commercialPorts) add(pnt, mode === "rail" ? 3.4 : mode === "expressway" ? 3.2 : 2.9, "port");
    if (mode !== "expressway") {
      for (const m of markets) add(m, mode === "rail" ? 2.1 : 2.4, "market");
      for (const pss of passes) add(pss, mode === "rail" ? 0.4 : 1.5, "pass");
    }
    if (mode === "national" || mode === "local") {
      for (const v of villages) add(v, mode === "local" ? 1.7 : 0.9, "village");
    }
    return dedupePointCandidates(anchors.sort((a, b) => b.score - a.score), 4.5).slice(0, max);
  }

  function preferenceCellAnchors(mode, potentialField, max = 36, minDistance = 9) {
    const step = mode === "expressway" ? 5 : mode === "rail" ? 4 : 4;
    const candidates = [];
    for (let y = 2; y < MAP_H - 2; y += step) {
      for (let x = 2; x < MAP_W - 2; x += step) {
        const i = indexOf(x, y);
        if (sea[i] || regionIdAt(x, y) < 0) continue;
        const pass = passSuitability?.[i] || 0;
        const base = potentialField[i] || 0;
        let score = base * 2.2 + settlementDemand[i] * (mode === "expressway" ? 0.35 : 0.65) + valleyField[i] * 0.42 + coastalLowland[i] * 0.20 + plain[i] * 0.12 + pass * (mode === "rail" ? 0.08 : 0.22);
        if (mode === "expressway") score += logisticsPreSuitability[i] * 0.65 - settlementDemand[i] * 0.10;
        if (mode === "rail") score += preliminaryTownInfluence[i] * 0.42 - slope[i] * 0.55;
        if (mode === "national") score += preliminaryVillageInfluence[i] * 0.24 + crossingSuitability[i] * 0.20;
        score -= ridgeField[i] * (mode === "expressway" ? 0.50 : mode === "rail" ? 0.65 : 0.30);
        score -= Math.max(0, elevation[i] - 0.62) * (mode === "expressway" ? 1.4 : mode === "rail" ? 1.7 : 0.85);
        score += hash2(x, y, seed + 17100 + (mode === "rail" ? 17 : mode === "expressway" ? 31 : 0)) * 0.055;
        if (score > (mode === "expressway" ? 0.64 : mode === "rail" ? 0.58 : 0.52)) candidates.push({ x, y, score, regionId: regionIdAt(x, y), role: "preference-cell" });
      }
    }
    return pickEntities(candidates, { max, minDistance, seed: seed + 17200 + (mode === "rail" ? 23 : mode === "expressway" ? 41 : 0) });
  }

  function transportCandidatePoints(mode, potentialField, options = {}) {
    const maxCells = options.maxCells ?? (mode === "expressway" ? 18 : mode === "rail" ? 26 : 40);
    const maxSettlements = options.maxSettlements ?? (mode === "expressway" ? 20 : mode === "rail" ? 34 : 56);
    const minDistance = options.minDistance ?? (mode === "expressway" ? 13 : mode === "rail" ? 10 : 8);
    const cells = preferenceCellAnchors(mode, potentialField, maxCells, minDistance);
    const settlements = modeSettlementAnchors(mode, potentialField, maxSettlements);
    return dedupePointCandidates([...settlements, ...cells].sort((a, b) => b.score - a.score), Math.max(4, minDistance * 0.55))
      .slice(0, options.maxTotal ?? (mode === "expressway" ? 32 : mode === "rail" ? 48 : 74));
  }

  function routeBetweenTrafficCandidates(a, b, mode, costField, penaltyField, options = {}) {
    if (!a || !b) return [];
    const start = { x: Math.round(a.x), y: Math.round(a.y), regionId: regionIdAt(Math.round(a.x), Math.round(a.y)) };
    const target = { x: Math.round(b.x), y: Math.round(b.y) };
    if (!inside(start.x, start.y) || !inside(target.x, target.y)) return [];
    if (sea[indexOf(start.x, start.y)] || sea[indexOf(target.x, target.y)]) return [];
    const d = Math.hypot(start.x - target.x, start.y - target.y);
    const snap = options.snapRadius ?? (mode === "expressway" ? 3 : mode === "rail" ? 2.5 : 2.25);
    const searchPad = options.searchPad ?? Math.ceil(Math.max(18, Math.min(58, d * (mode === "expressway" ? 0.46 : mode === "rail" ? 0.42 : 0.36))));
    const bounds = options.bounds || {
      minX: Math.max(0, Math.min(start.x, target.x) - searchPad),
      maxX: Math.min(MAP_W - 1, Math.max(start.x, target.x) + searchPad),
      minY: Math.max(0, Math.min(start.y, target.y) - searchPad),
      maxY: Math.min(MAP_H - 1, Math.max(start.y, target.y) + searchPad),
    };
    const path = traceCorridorByCost(
      start,
      (x, y) => Math.hypot(x - target.x, y - target.y) <= snap,
      costField,
      penaltyField || null,
      {
        curvePenalty: options.curvePenalty ?? (mode === "expressway" ? 0.12 : mode === "rail" ? 0.15 : mode === "national" ? 0.070 : 0.040),
        penaltyStrength: options.penaltyStrength ?? (mode === "expressway" ? 1.9 : mode === "rail" ? 1.35 : mode === "national" ? 1.05 : 0.78),
        minGoalDistance: Math.min(8, Math.max(2, d * 0.08)),
        keepRegion: false,
        maxExpanded: Math.min(SIZE, Math.max(1800, Math.floor(d * d * (mode === "expressway" ? 3.8 : mode === "rail" ? 4.6 : 5.2)))),
        terrainFlowBias: options.terrainFlowBias ?? (mode === "expressway" ? 0.10 : mode === "rail" ? 0.12 : mode === "national" ? 0.24 : 0.30),
        surfaceGrain: options.surfaceGrain ?? (mode === "national" ? 0.030 : mode === "local" ? 0.042 : 0.010),
        bounds,
        goalHint: options.goalHint || target,
        heuristicWeight: options.heuristicWeight ?? (mode === "expressway" ? 0.66 : mode === "rail" ? 0.54 : mode === "national" ? 0.46 : 0.30),
      }
    );
    if (path.length < 4) return [];
    if (options.maxPathLength && pathLengthCells(path) > options.maxPathLength) return [];
    const relaxed = relaxRouteToTerrain(path, costField, {
      radius: options.relaxRadius ?? (mode === "national" ? 2 : 1),
      lineWeight: options.relaxLineWeight ?? (mode === "national" ? 0.30 : mode === "expressway" ? 0.24 : 0.48),
      grain: options.surfaceGrain ?? 0.020,
      iterations: 1,
    });
    const crossesSea = (candidate) => (candidate || []).some(([x, y]) => !inside(x, y) || sea[indexOf(x, y)] || costField[indexOf(x, y)] >= INF);
    const chosen = relaxed.length >= Math.max(3, path.length * 0.55) ? relaxed : path;
    if (crossesSea(chosen)) return crossesSea(path) ? [] : path;
    if (!routePhysicalAcceptable(chosen, mode, options)) return routePhysicalAcceptable(path, mode, options) ? path : [];
    return chosen;
  }

function pruneParallelSameMode(paths, mode, potentialField, options = {}) {
    if (!paths?.length) return { mode, pruned: 0, kept: 0 };
    const scored = paths.map((path, originalIndex) => {
      const len = pathLengthCells(path);
      return { path, originalIndex, len, score: pathAverageField(path, potentialField) * 12 + Math.log1p(len) + (len > 80 ? 0.30 : 0) };
    }).sort((a, b) => b.score - a.score);
    const accepted = new Uint8Array(SIZE);
    const kept = [];
    let pruned = 0;
    const radius = options.radius ?? (mode === "expressway" ? 3 : mode === "rail" ? 2 : 2);
    const threshold = options.threshold ?? (mode === "expressway" ? 0.54 : mode === "rail" ? 0.58 : 0.62);
    function mark(path) {
      for (const [px, py] of path) {
        for (let dy = -radius; dy <= radius; dy++) {
          for (let dx = -radius; dx <= radius; dx++) {
            if (dx * dx + dy * dy > radius * radius) continue;
            const x = px + dx, y = py + dy;
            if (inside(x, y)) accepted[indexOf(x, y)] = 1;
          }
        }
      }
    }
    function overlap(path) {
      let hit = 0, n = 0;
      for (const [x, y] of path) {
        if (!inside(x, y)) continue;
        n++;
        if (accepted[indexOf(x, y)]) hit++;
      }
      return n ? hit / n : 0;
    }
    for (const item of scored) {
      const ov = overlap(item.path);
      const shortStub = item.len < (options.shortLength ?? (mode === "expressway" ? 40 : mode === "rail" ? 30 : 22));
      if (kept.length >= (options.minKeep ?? 2) && ov > threshold && (shortStub || ov > threshold + 0.13)) {
        pruned++;
        continue;
      }
      kept.push(item);
      mark(item.path);
    }
    kept.sort((a, b) => a.originalIndex - b.originalIndex);
    paths.length = 0;
    paths.push(...kept.map((item) => item.path));
    return { mode, pruned, kept: kept.length };
  }

  function endpointList(paths) {
    const out = [];
    paths.forEach((path, pathIndex) => {
      if (!path || path.length < 2) return;
      const a = path[0];
      const b = path[path.length - 1];
      out.push({ x: a[0], y: a[1], pathIndex, atStart: true });
      out.push({ x: b[0], y: b[1], pathIndex, atStart: false });
    });
    return out;
  }

  function validTransportEndpoint(p, mode, baseInfluence, civicPoints = []) {
    if (!p || !inside(p.x, p.y)) return true;
    const i = indexOf(p.x, p.y);
    if (sea[i]) return true;
    if (p.x < 4 || p.y < 4 || p.x > MAP_W - 5 || p.y > MAP_H - 5) return true;
    const nearNetwork = (baseInfluence?.[i] || 0) > (mode === "expressway" ? 0.22 : 0.16);
    const civicRadius = mode === "expressway" ? 9 : mode === "rail" ? 7 : 6;
    const nearCivic = civicPoints.some((q) => Math.hypot(q.x - p.x, q.y - p.y) < civicRadius);
    const settlementLike = settlementDemand[i] > (mode === "expressway" ? 0.22 : 0.14) || preliminaryTownInfluence[i] > 0.16 || preliminaryVillageInfluence[i] > 0.22;
    const passTerminus = mode !== "expressway" && (passSuitability?.[i] || 0) > 0.58 && settlementDemand[i] > 0.08;
    return nearNetwork || nearCivic || (settlementLike && naturalEndpoint(p.x, p.y, i, mode)) || passTerminus;
  }

  function repairDanglingTransportEndpoints(paths, mode, costField, targetPaths, potentialField, options = {}) {
    const debug = { mode, added: [], checked: 0 };
    if (!paths?.length) return debug;
    const endpoints = endpointList(paths);
    const civic = dedupePointCandidates([
      ...modernCities, ...markets, ...commercialPorts, ...(mode === "national" || mode === "local" ? villages : []),
    ], 4);
    const baseInfluence = cachedInfluenceFromPaths(targetPaths || [], options.targetRadius ?? (mode === "expressway" ? 7 : 5), `${mode}:endpoint-base`);
    const candidates = transportCandidatePoints(mode === "local" ? "national" : mode, potentialField, {
      maxCells: mode === "expressway" ? 10 : 20,
      maxSettlements: mode === "expressway" ? 12 : 36,
      maxTotal: mode === "expressway" ? 20 : 52,
      minDistance: mode === "expressway" ? 14 : 7,
    });
    const maxAdded = options.maxAdded ?? (mode === "expressway" ? 3 : mode === "rail" ? 5 : 16);
    const penalty = cachedInfluenceFromPaths(paths, options.penaltyRadius ?? 5, `${mode}:endpoint-penalty`);
    for (const ep of endpoints) {
      if (debug.added.length >= maxAdded) break;
      debug.checked++;
      if (validTransportEndpoint(ep, mode, baseInfluence, civic)) continue;
      const target = [...civic, ...candidates]
        .filter((q) => Math.hypot(q.x - ep.x, q.y - ep.y) >= (options.minTargetDistance ?? 7))
        .map((q) => {
          const d = Math.hypot(q.x - ep.x, q.y - ep.y);
          return { q, d, score: d / Math.sqrt(Math.max(0.35, q.score || 0.7)) };
        })
        .filter((e) => e.d <= (options.maxTargetDistance ?? (mode === "expressway" ? 58 : mode === "rail" ? 48 : 42)))
        .sort((a, b) => a.score - b.score)[0]?.q;
      if (!target) continue;
      const path = routeBetweenTrafficCandidates(ep, target, mode, costField, penalty, {
        ...options,
        snapRadius: options.snapRadius ?? 3,
        maxPathLength: options.maxPathLength ?? (mode === "expressway" ? 76 : 58),
      });
      if (path.length < 4) continue;
      if (!transportRouteAcceptable(path, mode, potentialField, penalty, { minLength: 4, maxLength: options.maxPathLength ?? (mode === "expressway" ? 76 : 58) })) continue;
      paths.push(path);
      debug.added.push({ mode: `${mode}-endpoint-repair`, path, from: "dangling-end", to: target.role || target.kind || "candidate" });
      addCorridorInfluencePenalty(penalty, path, options.penaltyRadius ?? 5, options.addedPenalty ?? 0.20);
    }
    return debug;
  }

  function pruneDanglingTerminalSegments(paths, mode, targetPaths, options = {}) {
    const debug = { mode, pruned: 0, kept: 0 };
    if (!paths?.length) return debug;
    const civic = dedupePointCandidates([
      ...modernCities, ...markets, ...commercialPorts, ...(mode === "national" || mode === "local" ? villages : []),
    ], 4);
    const baseInfluence = cachedInfluenceFromPaths(targetPaths || [], options.targetRadius ?? (mode === "expressway" ? 7 : 5), `${mode}:terminal-base`);
    const oneInvalidMax = options.oneInvalidMax ?? (mode === "expressway" ? 34 : mode === "rail" ? 24 : mode === "local" ? 10 : 18);
    const bothInvalidMax = options.bothInvalidMax ?? (mode === "expressway" ? 54 : mode === "rail" ? 42 : mode === "local" ? 18 : 34);
    const minKeep = options.minKeep ?? 1;
    const kept = [];
    for (const path of paths) {
      if (!path || path.length < 2) continue;
      const first = path[0];
      const last = path[path.length - 1];
      const a = { x: first[0], y: first[1] };
      const b = { x: last[0], y: last[1] };
      const len = pathLengthCells(path);
      const invalidA = !validTransportEndpoint(a, mode, baseInfluence, civic);
      const invalidB = !validTransportEndpoint(b, mode, baseInfluence, civic);
      const prune = paths.length - debug.pruned > minKeep && ((invalidA && invalidB && len < bothInvalidMax) || ((invalidA || invalidB) && len < oneInvalidMax));
      if (prune) {
        debug.pruned++;
      } else {
        kept.push(path);
      }
    }
    paths.length = 0;
    paths.push(...kept);
    debug.kept = paths.length;
    return debug;
  }

const premodernRoads = [];
  const nationalRoads = [];
  const minorRoads = [];
  const railways = [];
  const branchRailways = [];
  const externalRoads = [];
  const externalRailways = [];
  const expressways = [];
  const ringRoads = [];
  const ringRailways = [];
  const ringExpressways = [];
  const externalExpressways = [];
  const icAccessRoads = [];
  const interchanges = [];
  const externalGateways = [];

  // Premodern roads connect castles/markets/ports sparsely.
  for (const c of castles) {
    const near = [...markets, ...ports, ...crossings].sort((a, b) => Math.hypot(a.x - c.x, a.y - c.y) - Math.hypot(b.x - c.x, b.y - c.y)).slice(0, 2);
    for (const n of near) {
      const path = routeLight(c, n, 2);
      if (path.length > 2) premodernRoads.push(path);
    }
  }

  // National roads, expressways, and rail are generated from unified OD demand.
  // Roads are built by the density-flow portal system below; rail waits until
  // external gateways exist so the node model can include outside-region demand.

  let railODDebug = null;

  // Expressways are generated later by the density-flow portal system.

  // External gateways at land edges; used by naming/UI and later transport work.
  for (const regionId of [...regionStats.keys()].sort((a, b) => a - b)) {
    const st = regionStats.get(regionId);
    if (!st || st.area < 140) continue;
    const edgeCandidates = [];
    for (let y = st.minY; y <= st.maxY; y += 3) {
      for (const x of [st.minX, st.maxX]) {
        if (!inside(x, y)) continue;
        const i = indexOf(x, y);
        if (!sea[i] && regionIdAt(x, y) === regionId && (x < 5 || y < 5 || x > MAP_W - 6 || y > MAP_H - 6)) edgeCandidates.push({ x, y, score: developable[i] + valleySettlement[i] });
      }
    }
    for (let x = st.minX; x <= st.maxX; x += 3) {
      for (const y of [st.minY, st.maxY]) {
        if (!inside(x, y)) continue;
        const i = indexOf(x, y);
        if (!sea[i] && regionIdAt(x, y) === regionId && (x < 5 || y < 5 || x > MAP_W - 6 || y > MAP_H - 6)) edgeCandidates.push({ x, y, score: developable[i] + valleySettlement[i] });
      }
    }
    const gateway = pickEntities(edgeCandidates, { max: (regionStats.get(regionId)?.area || 0) > 2200 ? 2 : 1, minDistance: 16, seed: seed + 13200 + regionId * 11 })[0];
    if (gateway) {
      gateway.kind = "External Gateway";
      gateway.regionId = regionId;
      externalGateways.push(gateway);
      const target = importantNodesForRegion(regionId)[0];
      if (target) {
        const path = routeLight(gateway, target, 3, transportFields.national);
        if (path.length > 2) externalRoads.push(path);
      }
    }
  }

  const railOD = buildUnifiedRailODNetwork({
    seed,
    sea, elevation, slope, ridgeField, valleyField, basinField, coastalLowland, plain, agriculture, naturalBarrierScore, passSuitability,
    transportFields, settlementDemand, preliminaryUrbanInfluence, preliminaryTownInfluence, preliminaryVillageInfluence,
    modernCities, markets, ports, commercialPorts, externalGateways, geographicUrbanAnchors,
    regionIdAt, routeBetweenTrafficCandidates, addCorridorInfluencePenalty, transportRouteAcceptable, pruneParallelSameMode, cachedInfluenceFromPaths,
  });
  railways.push(...railOD.railways);
  branchRailways.push(...railOD.branchRailways);
  railODDebug = railOD.debug;


  const { transportDebugLayers, runLocalAccessPass, stitchRasterNearContacts, stitchLongLocalBranches, sanitizeLocalRoads, downgradeShortNationalRoads, connectAllRoadNetworksFinal } = buildDensityFlowRoadTransportSystem({
    seed,
    sea, elevation, slope, ridgeField, valleyField, coastalLowland, naturalBarrierScore,
    agriculture, basinField, plain, passSuitability, crossingSuitability,
    settlementDemand, preliminaryVillageInfluence, preliminaryTownInfluence,
    logisticsPreSuitability, urbanEdge,
    transportFields, cachedInfluenceFromPaths,
    nationalRoads, minorRoads, railways, externalRoads, externalRailways,
    expressways, externalExpressways, icAccessRoads, interchanges, externalGateways,
    modernCities, markets, villages, ports, commercialPorts, passes, regionStats,
    regionIdAt, inFocusedPrefecture, importantNodesForRegion, dedupePointCandidates,
    routeBetweenTrafficCandidates, traceCorridorByCost, addCorridorInfluencePenalty,
    relaxRouteToTerrain, transportRouteAcceptable, repairTransportConnectivity,
    repairDanglingTransportEndpoints, pruneDanglingTerminalSegments, pruneParallelSameMode,
  });

  // Land-use road influence intentionally excludes expressways.  Expressways
  // are through-corridors here, not automatic suburbanization generators.
  // A narrow field controls land-use attachment, while a broader field raises
  // population density around trunk roads without painting a wide suburb band.
  const roadLanduseInfluence = cachedInfluenceFromPaths([...nationalRoads, ...ringRoads, ...externalRoads], 2.25, "road:landuse");
  const roadInfluence = cachedInfluenceFromPaths([...nationalRoads, ...ringRoads, ...externalRoads], 5.0, "road:influence");
  const roadDensityInfluence = cachedInfluenceFromPaths([...nationalRoads, ...ringRoads, ...externalRoads], 9.0, "road:density");
  const railInfluence2 = cachedInfluenceFromPaths([...railways, ...branchRailways, ...externalRailways], 4, "rail:influence");

  const stations = [];
  const usedStationKeys = new Set();
  function addStation(x, y, kind = "Station", score = 1) {
    x = Math.round(x); y = Math.round(y);
    if (!inside(x, y) || sea[indexOf(x, y)]) return false;
    const key = `${x},${y}`;
    if (usedStationKeys.has(key)) return false;
    usedStationKeys.add(key);
    stations.push({ x, y, kind, score, regionId: regionIdAt(x, y) });
    return true;
  }
  function stationIntervalForCell(i) {
    const urbanDensity = preliminaryUrbanInfluence[i] || 0;
    const ruralDensity = Math.max(preliminaryTownInfluence[i] || 0, preliminaryVillageInfluence[i] || 0);
    return urbanDensity > 0.50 ? 6 : ruralDensity > 0.24 ? 12 : 20;
  }
  function shouldPlaceRailStation(x, y, i) {
    const nearCity = (preliminaryUrbanInfluence[i] || 0) > 0.12 || distanceToNearest(modernCities, x, y) < 4.8;
    const nearTown = (preliminaryTownInfluence[i] || 0) > 0.13 || distanceToNearest(markets, x, y) < 4.2;
    const nearVillage = (preliminaryVillageInfluence[i] || 0) > 0.20 || distanceToNearest(villages, x, y) < 3.4;
    const lowland = plain[i] > 0.16 || coastalLowland[i] > 0.16 || basinField[i] > 0.20 || valleyField[i] > 0.22;
    const terrainOk = slope[i] < 0.42 && ridgeField[i] < 0.60 && elevation[i] < 0.76;
    return terrainOk && lowland && (nearCity || nearTown || nearVillage);
  }
  for (const city of modernCities) addStation(city.x, city.y, city.isPrefecturalCapital || city.isRegionalCapital ? "Major Station" : "Station", 1.5);
  for (const path of [...railways, ...branchRailways]) {
    let lastStation = null;
    for (const p of samplePath(path, 5)) {
      const x = Math.round(p.x);
      const y = Math.round(p.y);
      if (!inside(x, y)) continue;
      const i = indexOf(x, y);
      if (sea[i] || !shouldPlaceRailStation(x, y, i)) continue;
      const interval = stationIntervalForCell(i);
      if (lastStation && Math.hypot(x - lastStation.x, y - lastStation.y) < interval) continue;
      if (distanceToNearest(stations, x, y) < Math.max(4.5, interval * 0.38)) continue;
      if (addStation(x, y, "Station", 0.8)) lastStation = { x, y };
    }
  }
  const stationInfluence = influenceFromPoints(stations, 7, (s) => s.kind === "Major Station" ? 1.35 : 0.85);
  const stationDensityInfluence = influenceFromPoints(stations, 10, (s) => s.kind === "Major Station" ? 1.85 : 1.05);

  // --- 5. Approximate city/town influence and land-use ---------------------
  const cityInfluence = new Float32Array(SIZE);
  const coreInfluence = new Float32Array(SIZE);
  const oldTownInfluence = influenceFromPoints([...markets, ...castleTowns, ...ports], 7, (p) => p.kind === "Major Port" ? 1.2 : 0.9);
  const populationDensity = new Float32Array(SIZE);

  function addKernel(grid, p, radius, weight, exponent = 1.7, terrainWeighted = true, combine = "max") {
    const r = Math.ceil(radius);
    const angle = hash2(p.x, p.y, seed + 14901) * Math.PI * 2;
    const stretch = 1.35 + hash2(p.x, p.y, seed + 14902) * 0.85;
    const squeeze = 0.62 + hash2(p.x, p.y, seed + 14903) * 0.28;
    const ca = Math.cos(angle);
    const sa = Math.sin(angle);
    for (let dy = -r; dy <= r; dy++) {
      for (let dx = -r; dx <= r; dx++) {
        const x = p.x + dx;
        const y = p.y + dy;
        if (!inside(x, y)) continue;
        const i = indexOf(x, y);
        if (sea[i]) continue;
        const along = (dx * ca + dy * sa) / stretch;
        const across = (-dx * sa + dy * ca) / squeeze;
        const baseD = Math.hypot(along, across);
        const conduit = clamp(
          plain[i] * 0.18 +
          valleySettlement[i] * 0.26 +
          coastalSettlement[i] * 0.16 +
          roadDensityInfluence[i] * 0.34 +
          roadInfluence[i] * 0.22 +
          railInfluence2[i] * 0.26 +
          stationDensityInfluence[i] * 0.12
        );
        const barrier = clamp(
          slope[i] * 0.62 +
          ridgeField[i] * 0.74 +
          Math.max(0, elevation[i] - 0.58) * 0.82 +
          (river[i] > 0.68 ? 0.60 : river[i] > 0.34 ? 0.22 : 0)
        );
        const noise = 0.78 + hash2(x, y, seed + 14910 + Math.round((p.population || 0) / 1000)) * 0.46;
        const d = baseD * (1.10 - conduit * 0.42 + barrier * 0.62) * noise;
        if (d > radius) continue;
        const terrain = terrainWeighted ? clamp(0.06 + developable[i] * 1.08 + valleySettlement[i] * 0.24 + coastalSettlement[i] * 0.14 + conduit * 0.38 - barrier * 0.70, 0, 1.42) : 1;
        const v = weight * Math.pow(1 - d / Math.max(1, radius), exponent) * terrain;
        if (combine === "add") grid[i] = Math.min(3.4, grid[i] + v);
        else if (v > grid[i]) grid[i] = v;
      }
    }
  }

  for (const city of modernCities) {
    addKernel(cityInfluence, city, city.sprawlRadius || Math.round((city.urbanRadius || 10) * 1.4), (city.urbanWeight || 1.0) * (city.isRegionalCapital ? 0.38 : 0.30), 2.75, true, "add");
    addKernel(cityInfluence, city, city.urbanRadius || 10, city.urbanWeight || 1.0, 1.18, true, "add");
    addKernel(coreInfluence, city, city.coreRadius || 3, (city.urbanWeight || 1.0) * 1.10, 1.65, true, "max");
  }
  const townInfluence = influenceFromPoints(markets, 6, (m) => clamp((m.population || 10000) / 26000, 0.45, 1.25));

  // Industrial/logistics/new town placeholders remain lightweight.  They are
  // routed by land-use proximity rather than expensive search passes.
  const industrialZones = [];
  for (const p of [...commercialPorts, ...modernCities.slice(0, 5)]) {
    const candidates = [];
    for (let dy = -10; dy <= 10; dy++) {
      for (let dx = -10; dx <= 10; dx++) {
        const x = p.x + dx;
        const y = p.y + dy;
        if (!inside(x, y)) continue;
        const i = indexOf(x, y);
        if (sea[i]) continue;
        const d = Math.hypot(dx, dy);
        if (d < 3 || d > 10) continue;
        const score = coastalLowland[i] * 0.22 + developable[i] * 0.22 + roadInfluence[i] * 0.20 + plain[i] * 0.12 - slope[i] * 0.25 + hash2(x, y, seed + 14000) * 0.06;
        if (score > 0.22) candidates.push({ x, y, score, kind: "Industrial Zone", regionId: regionIdAt(x, y) });
      }
    }
    const z = pickEntities(candidates, { max: 1, minDistance: 6, seed: seed + 14010 + p.x * 3 + p.y })[0];
    if (z && industrialZones.every((q) => Math.hypot(q.x - z.x, q.y - z.y) > 13)) industrialZones.push(z);
    if (industrialZones.length >= 8) break;
  }
  const industrialInfluence = influenceFromPoints(industrialZones, 5, () => 1.0);

  const satelliteCities = [];
  const newTowns = [];
  const logisticsScore = new Float32Array(SIZE);
  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 flatAgriculturalCorridor = clamp(
        agriculture[i] * 0.34 +
        plain[i] * 0.24 +
        developable[i] * 0.20 +
        basinField[i] * 0.12 +
        coastalLowland[i] * 0.08 +
        roadInfluence[i] * 0.26 +
        railInfluence2[i] * 0.18 +
        stationInfluence[i] * 0.10 -
        slope[i] * 0.44 -
        ridgeField[i] * 0.32
      );
      const nearMajorCity = modernCities.some((c) => Math.hypot(c.x - x, c.y - y) < 4);
      logisticsScore[i] = nearMajorCity ? 0 : flatAgriculturalCorridor;
    }
  }
  const logisticsParks = pickGlobalPoints(logisticsScore, {
    threshold: 0.34,
    max: 18,
    minDistance: 12,
    seedOffset: 1450,
    predicate: (x, y, i) => logisticsScore[i] > 0.30 && (roadInfluence[i] > 0.10 || railInfluence2[i] > 0.08 || stationInfluence[i] > 0.08),
  }).map((p) => ({ ...p, kind: "Logistics Park", score: logisticsScore[indexOf(p.x, p.y)], population: 0 }));
  const logisticsInfluence = influenceFromPoints(logisticsParks, 4.8, () => 1.0);

  function addFinalLocalAccessForUnservedSettlements() {
    const accessInfluence = cachedInfluenceFromPaths([...nationalRoads, ...railways, ...externalRoads, ...minorRoads], 8, "final-local:access");
    const candidates = [
      ...markets.filter((p) => (p.population || 0) >= 1800),
      ...ports,
      ...logisticsParks,
      ...villages.filter((p) => (p.population || 0) >= 450),
    ]
      .filter((p) => inside(p.x, p.y) && !sea[indexOf(p.x, p.y)] && accessInfluence[indexOf(p.x, p.y)] < 0.30)
      .map((p) => {
        const i = indexOf(p.x, p.y);
        const remoteness = Math.max(0, 0.32 - accessInfluence[i]);
        const ruralValue = agriculture[i] * 0.42 + valleySettlement[i] * 0.34 + coastalSettlement[i] * 0.24 + ruralSuitability[i] * 0.34;
        return { ...p, score: (p.population || 1200) / 16000 + remoteness * 3.2 + ruralValue + (p.portClass ? 0.8 : 0) + (p.kind === "Logistics Park" ? 1.0 : 0) + transportFields.localPotential[i] };
      })
      .sort((a, b) => b.score - a.score)
      .slice(0, 130);
    const localPenalty = cachedInfluenceFromPaths(minorRoads, 4, "final-local:minor");
    runLocalAccessPass({ candidates, accessInfluence, localPenalty, maxAdded: 90, maxLength: 92, debugMode: "local-access", from: "unserved", to: "network" });
  }

  function addRuralRoadMeshConnectors() {
    const roadInfluenceNow = cachedInfluenceFromPaths([...nationalRoads, ...externalRoads, ...railways, ...minorRoads], 7, "municipal-local:access");
    const localPenalty = cachedInfluenceFromPaths(minorRoads, 3, "municipal-local:minor");
    const candidates = villages
      .filter((p) => inside(p.x, p.y) && !sea[indexOf(p.x, p.y)] && roadInfluenceNow[indexOf(p.x, p.y)] < 0.46)
      .map((p) => {
        const i = indexOf(p.x, p.y);
        const ruralScore = agriculture[i] * 0.55 + valleySettlement[i] * 0.42 + coastalSettlement[i] * 0.25 + ruralSuitability[i] * 0.42 + Math.max(0, 0.46 - roadInfluenceNow[i]) * 2.4;
        return { ...p, score: ruralScore + (p.population || 700) / 24000 + transportFields.localPotential[i] * 0.55 };
      })
      .filter((p) => p.score > 0.30)
      .sort((a, b) => b.score - a.score)
      .slice(0, 170);
    runLocalAccessPass({
      candidates,
      accessInfluence: roadInfluenceNow,
      localPenalty,
      maxAdded: 70,
      minSpacing: 3.0,
      maxLength: 76,
      debugMode: "rural-mesh",
      from: "rural-settlement",
      to: "local-network",
      targetPredicate: (x, y, i) => roadInfluenceNow[i] > 0.18 || localPenalty[i] > 0.07,
    });
  }
  addFinalLocalAccessForUnservedSettlements();
  addRuralRoadMeshConnectors();
  transportDebugLayers.contactStitches = stitchRasterNearContacts();
  transportDebugLayers.longLocalStitches = stitchLongLocalBranches();
  transportDebugLayers.shortNationalDowngradeFinal = downgradeShortNationalRoads(18, 10);
  transportDebugLayers.localSanitizationFinal = sanitizeLocalRoads();
  transportDebugLayers.contactStitchesAfterConnectivity = stitchRasterNearContacts();
  // Keep this as the final road topology operation.  Later sanitization can cut
  // the short connectors that intentionally merge isolated components.
  transportDebugLayers.finalRoadNetworkConnectivity = connectAllRoadNetworksFinal(96);

  function dedupeTransportPathSet(paths, options = {}) {
    const before = paths.length;
    const minLength = options.minLength ?? 0;
    const minPoints = options.minPoints ?? 2;
    const sampleStep = Math.max(1, options.sampleStep ?? 1);
    const seen = new Set();
    const kept = [];
    let removedDuplicates = 0;
    let removedTooShort = 0;
    for (const path of paths) {
      if (!path || path.length < minPoints) { removedTooShort++; continue; }
      const cleaned = [];
      for (const pt of path) {
        if (!pt || pt.length < 2) continue;
        const x = Math.round(pt[0]);
        const y = Math.round(pt[1]);
        if (!cleaned.length || cleaned[cleaned.length - 1][0] !== x || cleaned[cleaned.length - 1][1] !== y) cleaned.push([x, y]);
      }
      if (cleaned.length < minPoints || pathLengthCells(cleaned) < minLength) { removedTooShort++; continue; }
      const sample = (candidate) => candidate
        .map((pt, idx) => (idx % sampleStep === 0 || idx === candidate.length - 1) ? `${pt[0]},${pt[1]}` : '')
        .filter(Boolean)
        .join('|');
      const forward = sample(cleaned);
      const backward = sample([...cleaned].reverse());
      const sig = forward < backward ? forward : backward;
      if (seen.has(sig)) { removedDuplicates++; continue; }
      seen.add(sig);
      kept.push(cleaned);
    }
    paths.length = 0;
    paths.push(...kept);
    return { before, after: kept.length, removedDuplicates, removedTooShort };
  }

  function pruneTransportPathSet(paths, minLength = 0, minKeep = 0) {
    const ranked = (paths || [])
      .map((path) => ({ path, len: pathLengthCells(path) }))
      .filter((row) => row.path?.length >= 2)
      .sort((a, b) => b.len - a.len);
    const kept = [];
    let pruned = 0;
    for (const row of ranked) {
      if (row.len >= minLength || kept.length < minKeep) kept.push(row.path);
      else pruned++;
    }
    paths.length = 0;
    paths.push(...kept);
    return { before: ranked.length, after: kept.length, pruned, minLength, minKeep };
  }

  function downgradeBranchNationalSpurs(maxLength = 30, importantRadius = 6.5, junctionRadius = 2.6) {
    const importantNodes = [
      ...modernCities.filter((p) => p.isPrefecturalCapital || p.isRegionalCapital || (p.population || 0) >= 90000),
      ...ports.filter((p) => p.portClass === 'major' || p.portClass === 'regional'),
      ...externalGateways,
    ];
    const otherTrunks = [...externalRoads, ...expressways, ...externalExpressways];
    const kept = [];
    const downgraded = [];
    function endpointNearImportant(endpoint) {
      return importantNodes.some((node) => Math.hypot(node.x - endpoint[0], node.y - endpoint[1]) <= importantRadius);
    }
    function endpointTouchesOtherTrunk(endpoint, currentPath) {
      const [ex, ey] = endpoint;
      for (const path of [...nationalRoads, ...otherTrunks]) {
        if (path === currentPath) continue;
        for (const pt of path) {
          if (Math.hypot(pt[0] - ex, pt[1] - ey) <= junctionRadius) return true;
        }
      }
      return false;
    }
    for (const path of nationalRoads) {
      if (!path || path.length < 2) continue;
      const len = pathLengthCells(path);
      const a = path[0];
      const b = path[path.length - 1];
      const aImportant = endpointNearImportant(a);
      const bImportant = endpointNearImportant(b);
      const aTouch = endpointTouchesOtherTrunk(a, path);
      const bTouch = endpointTouchesOtherTrunk(b, path);
      const oneSidedBranch = (aTouch && !bTouch) || (!aTouch && bTouch);
      const importantEndpoints = Number(aImportant) + Number(bImportant);
      if (oneSidedBranch && len <= maxLength && importantEndpoints <= 1 && !(aImportant && bImportant)) downgraded.push(path);
      else kept.push(path);
    }
    nationalRoads.length = 0;
    nationalRoads.push(...kept);
    minorRoads.push(...downgraded);
    return { threshold: maxLength, downgraded: downgraded.length, kept: kept.length };
  }

  transportDebugLayers.postConnectivityDedup = {
    national: dedupeTransportPathSet(nationalRoads, { minLength: 0.95, sampleStep: 1 }),
    externalRoads: dedupeTransportPathSet(externalRoads, { minLength: 1.5, sampleStep: 1 }),
    minor: dedupeTransportPathSet(minorRoads, { minLength: 0.95, sampleStep: 2 }),
    expressways: dedupeTransportPathSet(expressways, { minLength: 4, sampleStep: 2 }),
    externalExpressways: dedupeTransportPathSet(externalExpressways, { minLength: 4, sampleStep: 2 }),
  };
  transportDebugLayers.postConnectivityShortNationalDowngrade = downgradeShortNationalRoads(32, 7);
  transportDebugLayers.postConnectivityBranchNationalDowngrade = downgradeBranchNationalSpurs(42);
  transportDebugLayers.postConnectivityNationalPrune = pruneTransportPathSet(nationalRoads, 18, 7);
  transportDebugLayers.postConnectivityMinorDedup = dedupeTransportPathSet(minorRoads, { minLength: 0.95, sampleStep: 2 });
  transportDebugLayers.postConnectivityLocalSanitization = sanitizeLocalRoads();

  function smoothRasterPath(path, passes = 1) {
    let current = (path || []).map(([x, y]) => [Math.round(x), Math.round(y)]);
    for (let pass = 0; pass < passes; pass++) {
      if (current.length < 3) break;
      const next = [current[0]];
      for (let i = 1; i < current.length - 1; i++) {
        const [ax, ay] = current[i - 1];
        const [bx, by] = current[i];
        const [cx, cy] = current[i + 1];
        const nx = Math.round((ax + bx * 2 + cx) / 4);
        const ny = Math.round((ay + by * 2 + cy) / 4);
        if (next[next.length - 1][0] !== nx || next[next.length - 1][1] !== ny) next.push([nx, ny]);
      }
      next.push(current[current.length - 1]);
      current = next;
    }
    return current;
  }

  function directBridgeTunnelConnector(a, b, maxSegment = 20) {
    if (!a || !b) return [];
    const d = Math.hypot(a.x - b.x, a.y - b.y);
    const steps = Math.max(2, Math.ceil(d));
    const path = [];
    for (let s = 0; s <= steps; s++) {
      const t = s / steps;
      const x = Math.round(a.x + (b.x - a.x) * t);
      const y = Math.round(a.y + (b.y - a.y) * t);
      if (!path.length || path[path.length - 1][0] !== x || path[path.length - 1][1] !== y) path.push([x, y]);
    }
    return routePhysicalAcceptable(path, 'expressway', { maxSeaRun: maxSegment, maxTunnelRun: Math.min(10, maxSegment), maxSeaShare: 0.70, maxTunnelShare: 0.18 }) ? path : [];
  }

  function ensureInterchangePoint(x, y, kind = 'Interchange', source = 'expressway-endpoint') {
    x = Math.round(x); y = Math.round(y);
    if (!inside(x, y) || sea[indexOf(x, y)]) return false;
    if ((interchanges || []).some((p) => Math.hypot(p.x - x, p.y - y) <= 2.5)) return false;
    interchanges.push({ x, y, kind, score: 1, source });
    return true;
  }

  function ensureExpresswayEndpointsHaveICs() {
    let added = 0;
    for (const path of [...expressways, ...externalExpressways]) {
      if (!path || path.length < 2) continue;
      const endpoints = [path[0], path[path.length - 1]];
      for (const [x, y] of endpoints) if (ensureInterchangePoint(x, y)) added++;
    }
    return { added, total: interchanges.length };
  }

  function ensureNationalRoadCoverageForTowns(minPopulation = 5000) {
    const nationalInfluence = cachedInfluenceFromPaths([...nationalRoads, ...externalRoads], 3.2, 'final-national:coverage');
    const towns = dedupePointCandidates([
      ...modernCities.filter((p) => (p.population || 0) >= minPopulation),
      ...markets.filter((p) => (p.population || 0) >= minPopulation),
      ...ports.filter((p) => (p.population || 0) >= minPopulation || p.portClass === 'regional' || p.portClass === 'major'),
    ], 3.5);
    const debug = { minPopulation, checked: towns.length, added: 0 };
    const nationalPenalty = cachedInfluenceFromPaths([...nationalRoads, ...externalRoads], 6, 'final-national:penalty');
    for (const town of towns) {
      const ti = indexOf(town.x, town.y);
      if ((nationalInfluence[ti] || 0) > 0.24) continue;
      const candidates = dedupePointCandidates([...modernCities, ...externalGateways, ...ports, ...markets], 6)
        .filter((q) => q !== town)
        .map((q) => ({ q, d: Math.hypot(q.x - town.x, q.y - town.y) }))
        .filter((row) => row.d >= 10 && row.d <= 90)
        .sort((a, b) => a.d - b.d);
      let addedPath = null;
      for (const cand of candidates.slice(0, 8)) {
        const path = routeBetweenTrafficCandidates(town, cand.q, 'national', transportFields.national, nationalPenalty, {
          curvePenalty: 0.055,
          penaltyStrength: 0.60,
          terrainFlowBias: 0.18,
          surfaceGrain: 0.018,
          relaxRadius: 2,
          relaxLineWeight: 0.35,
          maxPathLength: cand.d * 2.8 + 40,
          maxSeaRun: 10,
          maxTunnelRun: 10,
        });
        if (path.length >= 4 && transportRouteAcceptable(path, 'national', transportFields.national, nationalPenalty, { minLength: 4, maxLength: cand.d * 3.0 + 48, maxSeaRun: 10, maxTunnelRun: 10 })) {
          addedPath = path;
          break;
        }
      }
      if (addedPath) {
        nationalRoads.push(addedPath);
        debug.added++;
      }
    }
    return debug;
  }

  function ensureMajorCityExpresswayConnections(minPopulation = 100000) {
    const majorCities = modernCities
      .filter((c) => (c.population || 0) >= minPopulation)
      .sort((a, b) => (b.population || 0) - (a.population || 0));
    const debug = { minPopulation, cityCount: majorCities.length, added: 0, pairs: [] };
    if (majorCities.length < 2) return debug;

    const permissiveExpresswayCost = new Float32Array(SIZE);
    for (let i = 0; i < SIZE; i++) {
      const terrainBase = Number.isFinite(transportFields.expressway[i]) && transportFields.expressway[i] < INF
        ? transportFields.expressway[i]
        : 0.28 + Math.max(0, slope[i] - 0.18) * 0.9 + Math.max(0, elevation[i] - 0.58) * 1.6 + ridgeField[i] * 0.45;
      permissiveExpresswayCost[i] = sea[i]
        ? 0.72 + (naturalBarrierScore?.[i] || 0) * 0.16
        : terrainBase + Math.max(0, elevation[i] - 0.70) * 1.7;
    }

    const componentOfCities = () => {
      const parent = new Map();
      const keyOf = (city) => city.name || `${city.x},${city.y}`;
      function find(k) {
        const p = parent.get(k);
        if (p === k) return k;
        const r = find(p);
        parent.set(k, r);
        return r;
      }
      function union(a, b) {
        const ra = find(a); const rb = find(b);
        if (ra !== rb) parent.set(ra, rb);
      }
      for (const city of majorCities) parent.set(keyOf(city), keyOf(city));
      const paths = [...expressways, ...externalExpressways];
      for (const path of paths) {
        const near = majorCities.filter((city) => path.some(([x, y]) => Math.hypot(city.x - x, city.y - y) <= 7.5));
        if (near.length >= 2) {
          const k0 = keyOf(near[0]);
          for (let i = 1; i < near.length; i++) union(k0, keyOf(near[i]));
        }
      }
      return new Map(majorCities.map((city) => [keyOf(city), find(keyOf(city))]));
    };

    const pairPriority = (a, b) => {
      const d = Math.hypot(a.x - b.x, a.y - b.y);
      const pop = Math.sqrt((a.population || minPopulation) * (b.population || minPopulation));
      return d / Math.max(1, Math.log2(pop));
    };

    for (let iter = 0; iter < majorCities.length * 2; iter++) {
      const comps = componentOfCities();
      const reps = new Set(comps.values());
      if (reps.size <= 1) break;
      let best = null;
      for (const a of majorCities) {
        for (const b of majorCities) {
          if (a === b) continue;
          const ka = a.name || `${a.x},${a.y}`;
          const kb = b.name || `${b.x},${b.y}`;
          if (comps.get(ka) === comps.get(kb)) continue;
          const score = pairPriority(a, b);
          if (!best || score < best.score) best = { a, b, score, d: Math.hypot(a.x - b.x, a.y - b.y) };
        }
      }
      if (!best) break;
      let path = routeBetweenTrafficCandidates(best.a, best.b, 'expressway', permissiveExpresswayCost, null, {
        curvePenalty: 0.045,
        penaltyStrength: 0.18,
        terrainFlowBias: 0.03,
        surfaceGrain: 0.001,
        relaxRadius: 2,
        relaxLineWeight: 0.15,
        maxPathLength: best.d * 3.5 + 110,
        snapRadius: 4.0,
        searchPad: Math.ceil(Math.max(20, Math.min(96, best.d * 0.55 + 16))),
        heuristicWeight: 0.78,
        maxSeaRun: 20,
        maxTunnelRun: 20,
        maxSeaShare: 0.45,
        maxTunnelShare: 0.45,
      });
      if (!path.length) path = directBridgeTunnelConnector(best.a, best.b, 20);
      if (path.length >= 2) {
        path = smoothRasterPath(path, 2);
        expressways.push(path);
        debug.added++;
        debug.pairs.push({ from: best.a.name, to: best.b.name, distance: Math.round(best.d), length: Math.round(pathLengthCells(path)) });
      } else {
        break;
      }
    }
    debug.finalExpresswayCount = expressways.length;
    return debug;
  }

  transportDebugLayers.postConnectivityNationalCoverage = ensureNationalRoadCoverageForTowns(5000);
  transportDebugLayers.postConnectivityMajorCityExpresswayGuarantee = ensureMajorCityExpresswayConnections(100000);
  transportDebugLayers.postConnectivityExpresswayDedup = dedupeTransportPathSet(expressways, { minLength: 8, sampleStep: 2 });
  transportDebugLayers.postConnectivityExpresswayEndpointICs = ensureExpresswayEndpointsHaveICs();
  var landuse = new Uint8Array(SIZE);

  // Re-run land-use classification after landuse allocation.  The loop above is
  // intentionally inside a helper to keep all thresholds in one place.
  function classifyLanduse() {
    landuse.fill(LANDUSE.RURAL);
    let maxDensity = 0;
    const baseNoiseSeed = seed + 15000;
    const urbanCapacity = new Float32Array(SIZE);
    const ruralDensityFloor = new Float32Array(SIZE);
    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 transport = Math.max(roadLanduseInfluence[i] * 0.95, railInfluence2[i] * 0.95, stationInfluence[i] * 0.82);
        const densityTransport = Math.max(roadDensityInfluence[i] * 0.95, stationDensityInfluence[i] * 1.05, railInfluence2[i] * 0.85);
        const urban = cityInfluence[i] * 0.76 + stationInfluence[i] * 0.30 + roadInfluence[i] * 0.14 + railInfluence2[i] * 0.10;
        const core = coreInfluence[i];
        const oldTown = oldTownInfluence[i] * 0.70 + townInfluence[i] * 0.38;
        const rural = villageInfluence[i] * 0.24 + ruralSuitability[i] * 0.30;
        const riverUrban = clamp(river[i] * 0.12 + valleyField[i] * 0.10 + plain[i] * 0.08 + basinField[i] * 0.08 - floodplain[i] * 0.10);
        urbanCapacity[i] = clamp(
          developable[i] * 0.66 +
          plain[i] * 0.16 +
          basinField[i] * 0.16 +
          valleyField[i] * 0.16 +
          coastalLowland[i] * 0.12 +
          roadInfluence[i] * 0.14 + roadDensityInfluence[i] * 0.16 + transport * 0.08 +
          riverUrban * 0.14 -
          slope[i] * 0.18 -
          ridgeField[i] * 0.12 -
          floodplain[i] * 0.08
        );
        const highPenaltyDensity = Math.max(0, elevation[i] - 0.58);
        const agrarianDensity = clamp(
          agriculture[i] * 0.045 +
          ruralSuitability[i] * 0.035 +
          developable[i] * 0.028 +
          plain[i] * 0.018 +
          basinField[i] * 0.014 +
          valleySettlement[i] * 0.014 +
          coastalSettlement[i] * 0.012 +
          villageInfluence[i] * 0.040 +
          townInfluence[i] * 0.022 +
          roadDensityInfluence[i] * 0.038 +
          stationDensityInfluence[i] * 0.020 +
          railInfluence2[i] * 0.012 -
          slope[i] * 0.030 -
          ridgeField[i] * 0.020 -
          highPenaltyDensity * 0.058
        );
        const remoteWilderness = elevation[i] > 0.60 && slope[i] > 0.34 && ridgeField[i] > 0.38 && densityTransport < 0.035 && villageInfluence[i] < 0.025 && townInfluence[i] < 0.025 && cityInfluence[i] < 0.025;
        ruralDensityFloor[i] = remoteWilderness ? 0 : clamp(agrarianDensity, 0, elevation[i] > 0.62 || slope[i] > 0.42 ? 0.052 : 0.105);
        populationDensity[i] = clamp(
          urban * 0.66 +
          core * 0.46 +
          oldTown * 0.28 +
          townInfluence[i] * 0.16 +
          villageInfluence[i] * 0.14 +
          roadDensityInfluence[i] * 0.42 +
          stationDensityInfluence[i] * 0.34 +
          railInfluence2[i] * 0.12 +
          transport * 0.05 +
          agrarianDensity * 0.34
        );
        maxDensity = Math.max(maxDensity[i]);

        if (elevation[i] > 0.67 || (slope[i] > 0.56 && ridgeField[i] > 0.30) || ridgeField[i] > 0.70) {
          landuse[i] = LANDUSE.FOREST;
          continue;
        }
        if (industrialInfluence[i] > 0.22 && urbanCapacity[i] > 0.10) {
          landuse[i] = LANDUSE.INDUSTRIAL;
          continue;
        }
        if (logisticsInfluence[i] > 0.24 && urbanCapacity[i] > 0.08 && (roadInfluence[i] > 0.08 || railInfluence2[i] > 0.06)) {
          landuse[i] = LANDUSE.LOGISTICS;
          continue;
        }
        if (core > 0.38 && urbanCapacity[i] > 0.10) {
          landuse[i] = LANDUSE.CBD;
          continue;
        }
        if (oldTown > 0.18 && urbanCapacity[i] > 0.09) {
          landuse[i] = LANDUSE.OLD_URBAN;
          continue;
        }

        const suburbanity = urban * 0.88 + transport * 0.23 + stationInfluence[i] * 0.12 + townInfluence[i] * 0.08 + riverUrban * 0.08;
        const edgeTaper = clamp(cityInfluence[i] * 0.52 + stationInfluence[i] * 0.16 + roadLanduseInfluence[i] * 0.10 + 0.28);
        const sprawlBias = clamp(0.58 + hash2(x, y, baseNoiseSeed) * 0.42);
        const sprawlScore = suburbanity * sprawlBias * edgeTaper - core * 0.12;
        if (sprawlScore > 0.24 && urbanCapacity[i] > 0.10) {
          landuse[i] = LANDUSE.SUBURB;
        } else if (roadLanduseInfluence[i] > 0.30 && urbanCapacity[i] > 0.10 && (townInfluence[i] > 0.05 || cityInfluence[i] > 0.11)) {
          landuse[i] = LANDUSE.SUBURB;
        } else if (agriculture[i] > 0.16 || rural > 0.18 || (developable[i] > 0.13 && plain[i] > 0.13) || (basinField[i] > 0.18 && slope[i] < 0.34) || (coastalLowland[i] > 0.16 && slope[i] < 0.32)) {
          landuse[i] = LANDUSE.FARMLAND;
        } else {
          const usablePlain = slope[i] < 0.30 && (plain[i] > 0.18 || developable[i] > 0.20 || basinField[i] > 0.20 || coastalLowland[i] > 0.18);
          landuse[i] = elevation[i] > 0.58 || slope[i] > 0.38 ? LANDUSE.FOREST : usablePlain ? LANDUSE.FARMLAND : LANDUSE.RURAL;
        }
      }
    }

    const baseLanduse = landuse.slice();
    const isBuilt = (lu) => lu >= LANDUSE.OLD_URBAN && lu <= LANDUSE.ROADSIDE;
    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] || baseLanduse[i] === LANDUSE.FOREST) continue;
        const transport = Math.max(roadLanduseInfluence[i] * 0.95, railInfluence2[i] * 0.95, stationInfluence[i] * 0.82);
        const densityTransport = Math.max(roadDensityInfluence[i] * 0.95, stationDensityInfluence[i] * 1.05, railInfluence2[i] * 0.85);
        let urbanNeighbors = 0;
        let cbdNeighbors = 0;
        for (let dy = -1; dy <= 1; dy++) {
          for (let dx = -1; dx <= 1; dx++) {
            if (!dx && !dy) continue;
            const lu = baseLanduse[indexOf(x + dx, y + dy)];
            if (isBuilt(lu)) urbanNeighbors++;
            if (lu === LANDUSE.CBD) cbdNeighbors++;
          }
        }
        if (baseLanduse[i] === LANDUSE.OLD_URBAN && coreInfluence[i] > 0.31 && cbdNeighbors >= 3) {
          landuse[i] = LANDUSE.CBD;
          continue;
        }
        if ((baseLanduse[i] === LANDUSE.FARMLAND || baseLanduse[i] === LANDUSE.RURAL) && urbanCapacity[i] > 0.10) {
          const fringeChance = urbanNeighbors * 0.055 + cityInfluence[i] * 0.13 + transport * 0.12 + stationInfluence[i] * 0.08;
          const noise = 0.23 + hash2(x, y, seed + 15050) * 0.24;
          if (fringeChance > 0.34 + noise) {
            landuse[i] = LANDUSE.SUBURB;
          }
        }
        if ((river[i] > 0.10 || railInfluence2[i] > 0.16 || roadLanduseInfluence[i] > 0.22) && urbanNeighbors >= 3 && landuse[i] <= LANDUSE.FARMLAND && urbanCapacity[i] > 0.09) {
          landuse[i] = cbdNeighbors >= 1 || coreInfluence[i] > 0.15 ? LANDUSE.OLD_URBAN : LANDUSE.SUBURB;
        }
      }
    }

    for (const park of logisticsParks) {
      const r = 3;
      for (let dy = -r; dy <= r; dy++) {
        for (let dx = -r; dx <= r; dx++) {
          const x = park.x + dx;
          const y = park.y + dy;
          if (!inside(x, y)) continue;
          const i = indexOf(x, y);
          if (sea[i] || landuse[i] === LANDUSE.CBD || landuse[i] === LANDUSE.FOREST) continue;
          if (Math.hypot(dx, dy) <= r && (agriculture[i] > 0.18 || plain[i] > 0.15 || roadInfluence[i] > 0.06 || railInfluence2[i] > 0.05)) {
            landuse[i] = LANDUSE.LOGISTICS;
          }
        }
      }
    }

    if (maxDensity > 0) {
      for (let i = 0; i < SIZE; i++) {
        if (sea[i]) continue;
        const lu = landuse[i];
        let floor = ruralDensityFloor[i];
        if (lu === LANDUSE.FARMLAND) {
          floor = Math.max(floor, clamp(0.024 + agriculture[i] * 0.044 + ruralSuitability[i] * 0.024 + roadDensityInfluence[i] * 0.030 + stationDensityInfluence[i] * 0.026 + villageInfluence[i] * 0.018, 0, 0.110));
        } else if (lu === LANDUSE.LOGISTICS) {
          floor = Math.max(floor, clamp(0.018 + roadDensityInfluence[i] * 0.026 + railInfluence2[i] * 0.014 + logisticsInfluence[i] * 0.012, 0, 0.060));
        } else if (lu === LANDUSE.RURAL) {
          floor = Math.max(floor, clamp(0.012 + ruralSuitability[i] * 0.020 + developable[i] * 0.014 + roadDensityInfluence[i] * 0.024 + stationDensityInfluence[i] * 0.020, 0, 0.070));
        } else if (lu === LANDUSE.FOREST) {
          floor = Math.min(floor, (roadDensityInfluence[i] > 0.04 || villageInfluence[i] > 0.03) ? 0.026 : 0);
        }
        const normalized = populationDensity[i] / maxDensity;
        populationDensity[i] = clamp(Math.max(normalized, floor));
        if (lu === LANDUSE.FOREST && floor === 0 && populationDensity[i] < 0.012) populationDensity[i] = 0;
      }
    }
  }
  classifyLanduse();

  for (const city of modernCities) {
    let urbanFootprintCells = 0;
    let coreFootprintCells = 0;
    const r = Math.ceil((city.urbanRadius || 8) * 1.3);
    for (let dy = -r; dy <= r; dy++) {
      for (let dx = -r; dx <= r; dx++) {
        const x = city.x + dx;
        const y = city.y + dy;
        if (!inside(x, y)) continue;
        const i = indexOf(x, y);
        if (sea[i]) continue;
        if (Math.hypot(dx, dy) > r) continue;
        if (landuse[i] >= LANDUSE.OLD_URBAN && landuse[i] <= LANDUSE.ROADSIDE) urbanFootprintCells++;
        if (landuse[i] === LANDUSE.CBD) coreFootprintCells++;
      }
    }
    city.urbanFootprintCells = urbanFootprintCells;
    city.coreFootprintCells = coreFootprintCells;
  }

  const transportDebug = {
    humanStageVersion: "v2-sparse-raster",
    aStarRoutes: 0,
    regionalNodeCount: [...regionStats.keys()].reduce((sum, regionId) => sum + importantNodesForRegion(regionId).length, 0),
    fieldCorridorTransport: false,
    unifiedODTransport: true,
    settlementHierarchy: settlementHierarchyDebug,
    railODTransport: railODDebug,
    expresswayFieldCorridors: expressways.length,
    railFieldCorridors: railways.length,
    nationalRoadFieldCorridors: nationalRoads.length,
    localRoadFieldCorridors: minorRoads.length,
    layers: transportDebugLayers,
    nationalRoadPopulationCoverage: 0,
    nationalRoadUncoveredPopulation: 0,
  };

  return {
    ports,
    geographicUrbanAnchors,
    crossings,
    passes,
    settlementCluster,
    settlementScore,
    villages,
    markets,
    castles,
    castleTowns,
    premodernRoads,
    minorRoads,
    modernCities,
    populationDensity,
    railways,
    branchRailways,
    ringRailways,
    externalRailways,
    stations,
    industrialZones,
    nationalRoads,
    ringRoads,
    expressways,
    ringExpressways,
    icAccessRoads,
    externalRoads,
    externalExpressways,
    interchanges,
    logisticsParks,
    satelliteCities,
    newTowns,
    landuse,
    stationInfluence,
    roadInfluence,
    roadDensityInfluence,
    stationDensityInfluence,
    railInfluence2,
    villageInfluence,
    externalGateways,
    cityPopulationCap,
    transportDebug,
  };
}