import { INF, MAP_H, MAP_W, SIZE, clamp, fbm, hash2, indexOf, inside, pickEntities, rand, valueNoise } from "./mapUtils.js";
import { distanceToNearest, influenceFromPaths, influenceFromPoints, samplePath } from "./mapGeneratorHelpers.js";
import { LANDUSE } from "./landuseCodes.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;

  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 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 spine = (arcSpineField?.[i] || 0) * 0.58 + (branchRidgeField?.[i] || 0) * 0.38;
      confluenceField[i] = confluence;

      developable[i] = 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
      );
      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
      );
      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
      );
      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.50 + valleySettlement[i] * 0.26 + coastalSettlement[i] * 0.18 + agriculture[i] * 0.34 + plain[i] * 0.12) * clusterNoise);
      ruralSuitability[i] = clamp(
        agriculture[i] * 0.54 +
        developable[i] * 0.30 +
        valleySettlement[i] * 0.24 +
        coastalSettlement[i] * 0.15 +
        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.24 +
        coastalSettlement[i] * 0.20 +
        confluence * 0.34 +
        basinField[i] * 0.16 +
        plain[i] * 0.18 +
        settlementCluster[i] * 0.22 -
        slope[i] * 0.34 -
        ridgeField[i] * 0.17 -
        spine * 0.10
      );
      settlementScore[i] = clamp(ruralSuitability[i] * 0.58 + townSuitability[i] * 0.34 + confluence * 0.10);
      const naturalBarrier = naturalBarrierScore?.[i] || 0;
      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,
        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++;
      st.developableSum += developable[i];
      if (developable[i] > 0.16) st.developableCells++;
      if (valleySettlement[i] > 0.24) st.valleyCells++;
      if (coastalSettlement[i] > 0.25) st.coastCells++;
      if (townSuitability[i] > 0.28) st.townCells++;
      if (plain[i] > 0.24) st.plainCells++;
      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" }));

  const villageScore = new Float32Array(SIZE);
  for (let i = 0; i < SIZE; i++) {
    if (sea[i]) continue;
    villageScore[i] = clamp(ruralSuitability[i] * 0.60 + agriculture[i] * 0.30 + plain[i] * 0.15 + valleySettlement[i] * 0.24 + coastalSettlement[i] * 0.16 + settlementCluster[i] * 0.16);
  }
  const villages = pickRegionalPoints(villageScore, {
    stride: 2,
    threshold: 0.18 + rand(seed, 1031) * 0.030,
    totalMax: 280,
    minDistance: 4,
    seedOffset: 1030,
    kind: "Village",
    quotaForRegion: (regionId, st) => {
      if (!st || st.developableCells < 10) return 0;
      const vf = visibilityFactor(regionId, st);
      const raw = (st.developableCells / 28 + st.plainCells / 42 + st.valleyCells / 34 + st.coastCells / 46 + 3.2) * vf;
      const min = st.area > 2600 ? 20 : st.area > 1400 ? 12 : st.area > 520 ? 5 : st.area > 220 ? 2 : 0;
      const max = st.area > 3600 ? 72 : st.area > 2200 ? 50 : st.area > 900 ? 25 : 10;
      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 = valleySettlement[i] > 0.42 ? "Valley Village" : coastalSettlement[i] > 0.43 ? "Coastal Village" : "Village";
    const population = Math.round((300 + Math.pow(rand(seed, 18000 + n * 17 + p.x * 3 + p.y), 1.85) * 4700 + ruralSuitability[i] * 2600) / 100) * 100;
    return { ...p, kind, population };
  });

  const 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 featurePull = Math.max(
        distanceToNearest(ports, x, y) < 8 ? 0.10 : 0,
        distanceToNearest(crossings, x, y) < 6 ? 0.06 : 0,
        confluenceField[i] * 0.16
      );
      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.58 +
        agriculture[i] * 0.16 +
        plain[i] * 0.14 +
        villageInfluence[i] * 0.40 +
        featurePull +
        valleyMouth +
        basinField[i] * 0.12 +
        plain[i] * 0.14 +
        coastalLowland[i] * 0.08 -
        slope[i] * 0.18 -
        ridgeField[i] * 0.08
      );
    }
  }

  const markets = pickRegionalPoints(marketScore, {
    stride: 2,
    threshold: 0.245 + rand(seed, 1041) * 0.035,
    totalMax: 110,
    minDistance: 6,
    seedOffset: 1040,
    kind: "Market Town",
    quotaForRegion: (regionId, st) => {
      if (!st || st.townCells < 8) return 0;
      const vf = visibilityFactor(regionId, st);
      const raw = (st.developableCells / 105 + st.plainCells / 125 + st.valleyCells / 105 + st.coastCells / 130 + 2.1) * vf;
      const min = st.area > 2600 ? 9 : st.area > 1200 ? 5 : st.area > 520 ? 2 : 0;
      const max = st.area > 3600 ? 30 : st.area > 2200 ? 22 : st.area > 800 ? 10 : 5;
      return Math.round(clamp(raw + rand(seed, 1043 + regionId * 23) * 0.8, min, max));
    },
    extraScore: (x, y, i) => (distanceToNearest(commercialPorts, x, y) < 8 ? 0.07 : 0) + confluenceField[i] * 0.08,
  }).map((p, n) => {
    const i = indexOf(p.x, p.y);
    const kind = coastalSettlement[i] > 0.45 && distanceToNearest(ports, p.x, p.y) < 9 ? "Port Town" : valleySettlement[i] > 0.42 ? "Valley Market Town" : "Market Town";
    const population = Math.round((4000 + Math.pow(rand(seed, 18100 + n * 19 + p.x * 5 + p.y), 1.50) * 24000 + marketScore[i] * 13000) / 1000) * 1000;
    return { ...p, kind, population };
  });

  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 = developable[i];
        if (dev < 0.04) continue;
        const radial = clamp(1 - d / Math.max(1, radius));
        const terrainMultiplier = clamp(0.60 + plain[i] * 0.28 + basinField[i] * 0.18 + coastalLowland[i] * 0.16 + valleyField[i] * 0.13 - slope[i] * 0.42 - ridgeField[i] * 0.18, 0.26, 1.24);
        capacity += dev * (900 + 6200 * Math.pow(radial, 1.25)) * terrainMultiplier * densityBias;
      }
    }
    return Math.max(26000, Math.round(capacity / 1000) * 1000);
  }

  const urbanCandidates = [
    ...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 ? 28 : st && st.area > 900 ? 24 : 20;
    const capacity = estimateUrbanCapacity(p, cityRadius, 1.0);
    const score =
      Math.log10(capacity + 1) * 0.72 +
      townSuitability[i] * 1.40 +
      developable[i] * 1.05 +
      confluenceField[i] * 0.22 +
      (p.candidateKind === "port" ? 0.48 : 0) +
      (p.candidateKind === "castleTown" ? 0.22 : 0) +
      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 / 720 + 0.9) * vf + rand(seed, 12100 + regionId * 17) * 1.2),
      st.area > 1600 ? 1 : 0,
      st.area > 3600 ? 6 : st.area > 2200 ? 4 : st.area > 900 ? 3 : 1
    );
    const selected = pickEntities(list, {
      max: maxCities,
      minDistance: 17,
      threshold: 0,
      seed: seed + 12110 + regionId * 313,
      jitter: 0.02,
    });
    for (const p of selected) {
      if (usedCitySites.some((q) => Math.hypot(q.x - p.x, q.y - p.y) < 12)) 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);
  for (const [rank, city] of modernCities.entries()) {
    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 isRegionalCapital = isFirstInRegion;
    const rawPop = isRegionalCapital
      ? 150000 + rand(seed, 12201 + city.regionId * 17) * 520000
      : 32000 + Math.pow(rand(seed, 12202 + rank * 19 + city.x), 0.7) * 260000;
    const capMultiplier = isRegionalCapital ? 1.10 : 1.0;
    const population = Math.round(Math.min(rawPop, city.capacity * capMultiplier) / 1000) * 1000;
    city.population = Math.max(isRegionalCapital ? 90000 : 24000, population);
    city.isPrefecturalCapital = isPrefecturalCapital;
    city.isRegionalCapital = isRegionalCapital;
    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 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. Lightweight corridors -------------------------------------------
  function routeLight(a, b, snapRadius = 3) {
    if (!a || !b) return [];
    const steps = Math.max(2, Math.ceil(Math.hypot(a.x - b.x, a.y - b.y) * 1.15));
    const out = [];
    let lastKey = "";
    for (let s = 0; s <= steps; s++) {
      const t = s / steps;
      const fx = a.x + (b.x - a.x) * t;
      const fy = a.y + (b.y - a.y) * t;
      let best = null;
      let bestCost = INF;
      const radius = snapRadius + (s > 0 && s < steps ? 1 : 0);
      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]) continue;
          const lineDist = Math.hypot(x - fx, y - fy);
          const cost = lineDist * 0.72 + corridorCost[i] * 0.62 - valleySettlement[i] * 0.34 - developable[i] * 0.18 + hash2(x, y, seed + 13000 + s) * 0.05;
          if (cost < bestCost) {
            bestCost = cost;
            best = [x, y];
          }
        }
      }
      if (!best) best = [Math.round(fx), Math.round(fy)];
      const key = `${best[0]},${best[1]}`;
      if (key !== lastKey) {
        out.push(best);
        lastKey = key;
      }
    }
    return out;
  }

  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);
  }

  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 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);
    }
  }

  for (const regionId of [...regionStats.keys()].sort((a, b) => a - b)) {
    const nodes = importantNodesForRegion(regionId);
    if (nodes.length < 2) continue;
    const connected = [nodes[0]];
    const remaining = nodes.slice(1);
    const maxEdges = (regionStats.get(regionId)?.area || 0) > 2200 ? Math.min(13, nodes.length + 3) : Math.min(7, nodes.length + 1);
    while (remaining.length && nationalRoads.length < 48) {
      let best = null;
      let bestScore = INF;
      for (const a of connected) {
        for (const b of remaining) {
          const d = Math.hypot(a.x - b.x, a.y - b.y);
          const score = d - (a.nodeWeight + b.nodeWeight) * 0.9;
          if (score < bestScore) { bestScore = score; best = { a, b }; }
        }
      }
      if (!best) break;
      const path = routeLight(best.a, best.b, 3);
      if (path.length > 2) nationalRoads.push(path);
      connected.push(best.b);
      remaining.splice(remaining.indexOf(best.b), 1);
      if (connected.length - 1 >= maxEdges) break;
    }

    // A few k-nearest shortcuts for urbanized regions.
    const urbanNodes = nodes.filter((p) => p.population || p.portClass).slice(0, (regionStats.get(regionId)?.area || 0) > 2200 ? 7 : 4);
    for (let i = 0; i < urbanNodes.length; i++) {
      const a = urbanNodes[i];
      const b = urbanNodes.slice(i + 1).sort((p, q) => Math.hypot(a.x - p.x, a.y - p.y) - Math.hypot(a.x - q.x, a.y - q.y))[0];
      if (!b || Math.hypot(a.x - b.x, a.y - b.y) > 48) continue;
      const path = routeLight(a, b, 3);
      if (path.length > 2) nationalRoads.push(path);
    }

    // Railways: only high-order cities/ports, as a lightweight placeholder.
    const railNodes = nodes.filter((p) => (p.population || 0) > 80000 || p.portClass === "major" || p.portClass === "regional").slice(0, (regionStats.get(regionId)?.area || 0) > 2200 ? 6 : 4);
    railNodes.sort((a, b) => a.x - b.x || a.y - b.y);
    for (let i = 1; i < railNodes.length; i++) {
      const path = routeLight(railNodes[i - 1], railNodes[i], 4);
      if (path.length > 4) railways.push(path);
    }
  }

  // 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);
        if (path.length > 2) externalRoads.push(path);
      }
    }
  }

  // Approximate expressways as a very small subset of top inter-city links.
  const topCities = modernCities.slice().sort((a, b) => (b.population || 0) - (a.population || 0)).slice(0, 6);
  for (let i = 1; i < topCities.length && expressways.length < 4; i++) {
    const a = topCities[i - 1];
    const b = topCities[i];
    if (Math.hypot(a.x - b.x, a.y - b.y) < 85) {
      const path = routeLight(a, b, 5);
      if (path.length > 5) expressways.push(path);
    }
  }

  // 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 = influenceFromPaths([...nationalRoads, ...ringRoads, ...externalRoads], 2.25);
  const roadInfluence = influenceFromPaths([...nationalRoads, ...ringRoads, ...externalRoads], 5.0);
  const roadDensityInfluence = influenceFromPaths([...nationalRoads, ...ringRoads, ...externalRoads], 9.0);
  const railInfluence2 = influenceFromPaths([...railways, ...branchRailways, ...externalRailways], 4);

  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;
    const key = `${x},${y}`;
    if (usedStationKeys.has(key)) return;
    usedStationKeys.add(key);
    stations.push({ x, y, kind, score, regionId: regionIdAt(x, y) });
  }
  for (const city of modernCities) addStation(city.x, city.y, city.isPrefecturalCapital || city.isRegionalCapital ? "Major Station" : "Station", 1.5);
  for (const path of railways) for (const p of samplePath(path, 14)) addStation(p.x, p.y, "Station", 0.8);
  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 interchanges = [];
  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);
  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, populationDensity[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),
    nationalRoadPopulationCoverage: 0,
    nationalRoadUncoveredPopulation: 0,
  };

  return {
    ports,
    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,
  };
}