City-Lifeline-Cockpit/tools/city-massing/build-town.mjs

/**
 * 沅陵县城「轻量体块」生成器 · 第二版（票 35，替代票 34 被否的第一版）。
 *
 * ⚠️ 示意体块，不是实测建筑：除 OSM 真实轮廓外，其余按真实路网/建成区程序化排布。
 *    产物 meta.json/README 均显著标注；不得用于量测、统计或应急决策。
 *
 * 对第一版三条否决理由的修复（票 34 Answer 记录的根因）：
 *   ① 漏顶      → 所有三角面做**程序化定向**（按「面重心 − 体中心」翻正绕序），
 *                 且材质 doubleSided —— 从机制上消灭"背面被剔除"整类漏面；
 *   ② 尺寸不搭边 → 单元尺寸**从 74 个真实轮廓的分位数推导**（长边中位 53m/p90 83m、
 *                 短边 27–45m），第一版的 9×7/14×10/26×11 全部过小；
 *   ③ 不贴合地形 → 高程网格 17×17 → **64×64**（~27m），并**逐栋取足印四角采样**：
 *                 基准 = min(四角) − 裙边 4m，楼体向下延伸埋入地下 ⇒ 数学上不悬空，
 *                 正确落位后地形自然遮挡（"模型压在 Cesium 之上"随之消失）。
 *   中近景可辨  → 楼层色带（每层一圈略深的墙）+ 平顶女儿墙 / 坡屋顶出檐。
 *
 * 地理配准由构造保证（同第一版）：经纬度 → 目标点 ENU → 3D Tiles transform。
 * 用法：node tools/city-massing/build-town.mjs
 */
import { readFileSync, writeFileSync, mkdirSync, readdirSync } from 'node:fs'
import { join } from 'node:path'

const GIS = 'D:/GIS/GIS'
const OSM = join(GIS, 'planet_110.059,28.258_110.744,28.625-shp/shpOutput')
const DATA_DIR = process.env.CITY_DATA_DIR || join(GIS, 'cityMassing')

const ORIGIN = { lon: 110.39028164, lat: 28.4436787 }
const WINDOW = { lon: 0.02, lat: 0.015 }
const CORE_RADIUS = 500
const CAP = 5000
const SKIRT = 4 // 每栋向下多埋 4m：斜坡上不悬空、不露缝

const DEG = Math.PI / 180
const M_PER_DEG_LAT = 110540
const M_PER_DEG_LON = 111320 * Math.cos(ORIGIN.lat * DEG)

const toLocal = (lon, lat) => ({
    x: (lon - ORIGIN.lon) * M_PER_DEG_LON,
    y: (lat - ORIGIN.lat) * M_PER_DEG_LAT,
})
const inWindow = (lon, lat) => Math.abs(lon - ORIGIN.lon) <= WINDOW.lon && Math.abs(lat - ORIGIN.lat) <= WINDOW.lat

function hash01(...nums) {
    let h = 2166136261
    for (const n of nums) {
        let v = Math.round(n * 1000) | 0
        for (let i = 0; i < 4; i++) {
            h ^= (v >> (i * 8)) & 0xff
            h = Math.imul(h, 16777619)
        }
    }
    return ((h >>> 0) % 100000) / 100000
}
const pick = (arr, r) => arr[Math.min(arr.length - 1, Math.floor(r * arr.length))]

// ---------------------------------------------------------------- 原料
const readJson = (p) => JSON.parse(readFileSync(p, 'utf8'))
const buildings = readJson(join(OSM, 'buildings.json')).features
const roads = readJson(join(OSM, 'roads.json')).features
const waterways = readJson(join(OSM, 'waterways.json')).features
const resa = readJson(join(GIS, 'output/resa.json')).features
const terrain = JSON.parse(readFileSync(join(DATA_DIR, 'terrain-grid.json'), 'utf8'))

/** 高程：64×64 网格双线性插值（缺失点回落到中心值）。 */
function elevationAt(lon, lat) {
    const gx = (lon - terrain.lon0) / terrain.dLon
    const gy = (lat - terrain.lat0) / terrain.dLat
    const n = terrain.n
    const cx = Math.max(0, Math.min(n - 1.001, gx))
    const cy = Math.max(0, Math.min(n - 1.001, gy))
    const x0 = Math.floor(cx), y0 = Math.floor(cy)
    const fx = cx - x0, fy = cy - y0
    const at = (x, y) => terrain.heights[y * n + x]
    const c = [at(x0, y0), at(x0 + 1, y0), at(x0, y0 + 1), at(x0 + 1, y0 + 1)].map((v) => (typeof v === 'number' ? v : 111.6))
    return (c[0] * (1 - fx) + c[1] * fx) * (1 - fy) + (c[2] * (1 - fx) + c[3] * fx) * fy
}

// 水系避让
const waterLines = []
for (const f of waterways) {
    const ring = f.geometry && f.geometry.type === 'LineString' ? f.geometry.coordinates : null
    if (!ring) continue
    const seg = []
    for (const [lon, lat] of ring) {
        if (inWindow(lon, lat)) seg.push(toLocal(lon, lat))
        else if (seg.length) { waterLines.push(seg); seg.length = 0 }
    }
    if (seg.length) waterLines.push(seg)
}
function distanceToWater(p) {
    let best = Infinity
    for (const line of waterLines) {
        for (let i = 1; i < line.length; i++) {
            const a = line[i - 1], b = line[i]
            const dx = b.x - a.x, dy = b.y - a.y
            const len2 = dx * dx + dy * dy || 1
            let t = Math.max(0, Math.min(1, ((p.x - a.x) * dx + (p.y - a.y) * dy) / len2))
            best = Math.min(best, Math.hypot(p.x - (a.x + t * dx), p.y - (a.y + t * dy)))
        }
    }
    return best
}

// 建成区
const resaPolys = []
for (const f of resa) {
    const ring = f.geometry && f.geometry.type === 'Polygon' ? f.geometry.coordinates[0] : null
    if (!ring || ring.length < 4) continue
    if (!ring.some(([lon, lat]) => inWindow(lon, lat))) continue
    resaPolys.push(ring.map(([lon, lat]) => toLocal(lon, lat)))
}
function inResa(p) {
    for (const ring of resaPolys) {
        let inside = false
        for (let i = 0, j = ring.length - 1; i < ring.length; j = i++) {
            const a = ring[i], b = ring[j]
            if ((a.y > p.y) !== (b.y > p.y) && p.x < ((b.x - a.x) * (p.y - a.y)) / (b.y - a.y || 1e-9) + a.x) inside = !inside
        }
        if (inside) return true
    }
    return false
}

const STREET_TYPES = new Set(['trunk', 'primary', 'secondary', 'tertiary', 'residential', 'unclassified', 'living_street'])
const streetSegments = []
for (const f of roads) {
    const type = (f.properties && f.properties.type) || ''
    if (!STREET_TYPES.has(type)) continue
    const ring = f.geometry && f.geometry.type === 'LineString' ? f.geometry.coordinates : null
    if (!ring) continue
    for (let i = 1; i < ring.length; i++) {
        const [lon1, lat1] = ring[i - 1]
        const [lon2, lat2] = ring[i]
        if (!inWindow(lon1, lat1) && !inWindow(lon2, lat2)) continue
        streetSegments.push({ a: toLocal(lon1, lat1), b: toLocal(lon2, lat2) })
    }
}

// ---------------------------------------------------------------- 可复用单元（尺寸来自真实轮廓分位数）
/**
 * 单元尺寸依据（tools 脚本对 74 个真实轮廓的实测分位数）：
 *   长边 中位 53m / p90 83m；短边 p25–p90 27–45m；面积 1000–3800㎡。
 * 据此（乘以 0.45 的体块系数 —— 体块只求轮廓感，不必占满整个真实地块）：
 *   house  长边 24–34m · block 长边 34–60m · slab 长边 60–84m
 */
const ARCHETYPES = {
    house: { long: [24, 34], short: [13, 18], floors: [2, 3], floorH: 3.1, pitched: true },
    block: { long: [34, 60], short: [16, 24], floors: [3, 6], floorH: 3.1, pitched: false },
    slab: { long: [60, 84], short: [22, 34], floors: [6, 11], floorH: 3.0, pitched: false },
}
// 墙色带：墙体主色 + 楼层带更深一档（中近景读出"层"）；屋色独立。
const WALLS = [
    { wall: 0xf1ede3, band: 0xd9d2c2, roof: 0x8c4a3a },
    { wall: 0xe8e4da, band: 0xcfc9bb, roof: 0x6f7478 },
    { wall: 0xdfe1e3, band: 0xc4c8cb, roof: 0x545a5f },
    { wall: 0xe3d9cd, band: 0xcabfae, roof: 0x7a5140 },
]

// ---------------------------------------------------------------- 几何发射（定向 pass 消灭漏面）
const buckets = new Map()
let placed = 0
const stats = { house: 0, block: 0, slab: 0, real: 0, skippedWater: 0 }
/** 回验样本：每栋楼的位置与基准高程（供 verify 脚本独立复核"逐栋贴地"）。 */
const verificationSample = []

function bucket(color) {
    if (!buckets.has(color)) buckets.set(color, { positions: [], normals: [], indices: [] })
    return buckets.get(color)
}

/** 把一个世界坐标三角形写入桶：**自动定向** —— 法线取"远离体中心"的方向，绕序随法线生成。
 *  第一版"漏顶"的根因是手写绕序出错 + 单面材质；这里从机制上保证每个面朝外。 */
function emitTri(b, p1, p2, p3, center) {
    // 面法线
    const ux = p2[0] - p1[0], uy = p2[1] - p1[1], uz = p2[2] - p1[2]
    const vx = p3[0] - p1[0], vy = p3[1] - p1[1], vz = p3[2] - p1[2]
    let nx = uy * vz - uz * vy
    let ny = uz * vx - ux * vz
    let nz = ux * vy - uy * vx
    const len = Math.hypot(nx, ny, nz) || 1
    nx /= len; ny /= len; nz /= len
    // 面重心指向"体中心反方向"即为朝外；若与法线同向则翻转（含绕序交换）
    const gx = (p1[0] + p2[0] + p3[0]) / 3 - center[0]
    const gy = (p1[1] + p2[1] + p3[1]) / 3 - center[1]
    const gz = (p1[2] + p2[2] + p3[2]) / 3 - center[2]
    const outward = nx * gx + ny * gy + nz * gz
    // outward = 法线 · (重心 − 体中心)：朝外的面这个点积应为**正**；
    // 第一版把翻转条件写反（>0 才翻），恰好一半面朝内 —— 这就是"漏顶"的机制根源。
    const flip = outward < 0
    const a = flip ? p3 : p1
    const bb = flip ? p2 : p2
    const c = flip ? p1 : p3
    const base = b.positions.length / 3
    for (const v of [a, bb, c]) b.positions.push(v[0], v[1], v[2])
    for (let i = 0; i < 3; i++) {
        const s = flip ? -1 : 1
        b.normals.push(nx * s, ny * s, nz * s)
    }
    b.indices.push(base, base + 1, base + 2)
}

function emitBuilding({ x, y, zBase, heading, type, floors, width, depth, palette }) {
    const spec = ARCHETYPES[type]
    const height = floors * spec.floorH
    const cos = Math.cos(heading), sin = Math.sin(heading)
    // 体中心（用于面定向；z 取半高）
    const center = [x, zBase + height / 2, -y]
    const put = (lx, ly, lz) => {
        const ex = lx * cos - ly * sin
        const ny = lx * sin + ly * cos
        return [x + ex, zBase + lz, -(y + ny)]
    }
    const hw = width / 2, hd = depth / 2
    const b = bucket(palette.wall)

    const corners = [put(-hw, -hd, 0), put(hw, -hd, 0), put(hw, hd, 0), put(-hw, hd, 0)]
    const top = [put(-hw, -hd, height), put(hw, -hd, height), put(hw, hd, height), put(-hw, hd, height)]

    // 四片墙
    for (let i = 0; i < 4; i++) {
        const j = (i + 1) % 4
        emitTri(b, corners[i], corners[j], top[j], center)
        emitTri(b, corners[i], top[j], top[i], center)
    }

    // 楼层色带：每层在墙面贴一圈更深的带（外凸 0.15m），中近景读出"层"
    const band = bucket(palette.band)
    for (let f = 1; f < floors; f++) {
        const z0 = f * spec.floorH - 0.35
        const z1 = f * spec.floorH - 0.05
        const b0 = [put(-hw - 0.15, -hd - 0.15, z0), put(hw + 0.15, -hd - 0.15, z0), put(hw + 0.15, hd + 0.15, z0), put(-hw - 0.15, hd + 0.15, z0)]
        const b1 = [put(-hw - 0.15, -hd - 0.15, z1), put(hw + 0.15, -hd - 0.15, z1), put(hw + 0.15, hd + 0.15, z1), put(-hw - 0.15, hd + 0.15, z1)]
        for (let i = 0; i < 4; i++) {
            const j = (i + 1) % 4
            emitTri(band, b0[i], b0[j], b1[j], center)
            emitTri(band, b0[i], b1[j], b1[i], center)
        }
    }

    // 屋顶
    const roof = bucket(palette.roof)
    if (spec.pitched) {
        const ridge = height + Math.min(3.2, width * 0.14)
        const r1 = put(0, -hd, ridge)
        const r2 = put(0, hd, ridge)
        // 出檐 0.5m
        const e = 0.5
        const eave = [put(-hw - e, -hd - e, height - 0.1), put(hw + e, -hd - e, height - 0.1), put(hw + e, hd + e, height - 0.1), put(-hw - e, hd + e, height - 0.1)]
        const slope = [0, 0, 0]
        slope[0] = 0
        slope[1] = Math.sin(Math.atan2(ridge - height, hw))
        slope[2] = -Math.cos(Math.atan2(ridge - height, hw))
        emitTri(roof, eave[0], eave[1], r1, center)
        emitTri(roof, eave[1], eave[2], r2, center)
        emitTri(roof, eave[2], eave[3], r2, center)
        emitTri(roof, eave[3], eave[0], r1, center)
        // 山墙两片
        emitTri(roof, eave[0], r1, r2, center)
        emitTri(roof, eave[2], r2, r1, center)
        void slope
    } else {
        // 平顶 + 女儿墙（高出屋面 0.6m、内收 0.4m），远景读出"楼"而不是"盒盖"
        emitQuadTop(roof, top, center)
        const p0 = [put(-hw + 0.4, -hd + 0.4, height + 0.6), put(hw - 0.4, -hd + 0.4, height + 0.6), put(hw - 0.4, hd + 0.4, height + 0.6), put(-hw + 0.4, hd + 0.4, height + 0.6)]
        const p1 = [put(-hw, -hd, height), put(hw, -hd, height), put(hw, hd, height), put(-hw, hd, height)]
        for (let i = 0; i < 4; i++) {
            const j = (i + 1) % 4
            emitTri(roof, p1[i], p1[j], p0[j], center)
            emitTri(roof, p1[i], p0[j], p0[i], center)
        }
    }
    placed += 1
    stats[type] += 1
}

/** 顶面四边形（绕序按定向 pass 自动处理）。 */
function emitQuadTop(b, top, center) {
    emitTri(b, top[0], top[1], top[2], center)
    emitTri(b, top[0], top[2], top[3], center)
}

// ---------------------------------------------------------------- 布点
function chooseType(x, y, r) {
    const dist = Math.hypot(x, y)
    if (dist < CORE_RADIUS) return r < 0.3 ? 'slab' : 'block'
    if (dist < CORE_RADIUS * 2.4) return r < 0.65 ? 'block' : 'house'
    return r < 0.2 ? 'block' : 'house'
}

function emitInstance(x, y, heading, type, r) {
    if (placed >= CAP) return
    const lon = ORIGIN.lon + x / M_PER_DEG_LON
    const lat = ORIGIN.lat + y / M_PER_DEG_LAT
    if (!inWindow(lon, lat)) return
    if (distanceToWater({ x, y }) < 28) { stats.skippedWater += 1; return }

    const spec = ARCHETYPES[type]
    const long = spec.long[0] + (spec.long[1] - spec.long[0]) * hash01(x, y, 7)
    const short = spec.short[0] + (spec.short[1] - spec.short[0]) * hash01(x, y, 11)
    const floors = spec.floors[0] + Math.floor(hash01(x, y, 13) * (spec.floors[1] - spec.floors[0] + 1))
    // 足印四角采地形，取最小值 - 裙边 ⇒ 不悬空（楼体下延埋入地下）
    const cornersLonLat = [
        [lon - long / 2 / M_PER_DEG_LON, lat - short / 2 / M_PER_DEG_LAT],
        [lon + long / 2 / M_PER_DEG_LON, lat - short / 2 / M_PER_DEG_LAT],
        [lon + long / 2 / M_PER_DEG_LON, lat + short / 2 / M_PER_DEG_LAT],
        [lon - long / 2 / M_PER_DEG_LON, lat + short / 2 / M_PER_DEG_LAT],
    ].map(([a, b]) => elevationAt(a, b))
    const zBase = Math.min(...cornersLonLat) - SKIRT
    verificationSample.push({ lon, lat, zBase: Math.round(zBase * 100) / 100, type })

    emitBuilding({
        x, y, zBase, heading,
        type, floors,
        width: heading % (Math.PI) < Math.PI / 2 ? long : short, // 长边沿街
        depth: heading % Math.PI < Math.PI / 2 ? short : long,
        palette: pick(WALLS, hash01(x, y, 3)),
    })
}

// 1) 真实 OSM 轮廓
for (const f of buildings) {
    const ring = f.geometry && f.geometry.type === 'Polygon' ? f.geometry.coordinates[0] : null
    if (!ring) continue
    const pts = ring.filter(([lon, lat]) => inWindow(lon, lat))
    if (pts.length < 4) continue
    const ex = pts.map(([lon]) => (lon - ORIGIN.lon) * M_PER_DEG_LON)
    const ny = pts.map(([, lat]) => (lat - ORIGIN.lat) * M_PER_DEG_LAT)
    const minX = Math.min(...ex), maxX = Math.max(...ex)
    const minY = Math.min(...ny), maxY = Math.max(...ny)
    const cx = (minX + maxX) / 2, cy = (minY + maxY) / 2
    const w = Math.max(12, maxX - minX), d = Math.max(10, maxY - minY)
    const lon = ORIGIN.lon + cx / M_PER_DEG_LON
    const lat = ORIGIN.lat + cy / M_PER_DEG_LAT
    if (distanceToWater({ x: cx, y: cy }) < 24) continue
    const r = hash01(cx, cy, 13)
    const area = w * d
    const type = area > 3200 ? 'slab' : area > 1200 ? 'block' : 'house'
    const spec = ARCHETYPES[type]
    const floors = 2 + Math.floor(r * 5)
    const cornerEls = [
        elevationAt(ORIGIN.lon + minX / M_PER_DEG_LON, ORIGIN.lat + minY / M_PER_DEG_LAT),
        elevationAt(ORIGIN.lon + maxX / M_PER_DEG_LON, ORIGIN.lat + minY / M_PER_DEG_LAT),
        elevationAt(ORIGIN.lon + maxX / M_PER_DEG_LON, ORIGIN.lat + maxY / M_PER_DEG_LAT),
        elevationAt(ORIGIN.lon + minX / M_PER_DEG_LON, ORIGIN.lat + maxY / M_PER_DEG_LAT),
    ]
    emitBuilding({
        x: cx, y: cy, zBase: Math.min(...cornerEls) - SKIRT,
        heading: 0, type, floors, width: w, depth: d,
        palette: pick(WALLS, r),
    })
    verificationSample.push({ lon, lat, zBase: Math.round((Math.min(...cornerEls) - SKIRT) * 100) / 100, type })
    placed += 1
    stats.real += 1
}

// 2) 沿街排布（长边沿街）
const SPACING = 34
for (const seg of streetSegments) {
    const dx = seg.b.x - seg.a.x, dy = seg.b.y - seg.a.y
    const len = Math.hypot(dx, dy)
    if (len < 14) continue
    const heading = Math.atan2(dy, dx)
    const steps = Math.floor(len / SPACING)
    for (let s = 0; s <= steps; s++) {
        const t = steps === 0 ? 0.5 : s / steps
        const px = seg.a.x + dx * t
        const py = seg.a.y + dy * t
        for (const side of [-1, 1]) {
            const off = 16 + 10 * hash01(px, py, side)
            const ox = px - Math.sin(heading) * off * side
            const oy = py + Math.cos(heading) * off * side
            const r = hash01(ox, oy, 19)
            emitInstance(ox, oy, heading, chooseType(ox, oy, hash01(ox, oy, 23)), r)
        }
    }
}

// 3) 建成区补白（朝向跟最近的路向难以一一对应，这里给小角度扰动）
const FILL = 52
for (const ring of resaPolys) {
    const xs = ring.map((p) => p.x), ys = ring.map((p) => p.y)
    for (let gx = Math.min(...xs); gx <= Math.max(...xs); gx += FILL) {
        for (let gy = Math.min(...ys); gy <= Math.max(...ys); gy += FILL) {
            const jx = gx + (hash01(gx, gy, 29) - 0.5) * FILL * 0.7
            const jy = gy + (hash01(gx, gy, 31) - 0.5) * FILL * 0.7
            if (!inResa({ x: jx, y: jy })) continue
            const r = hash01(jx, jy, 37)
            if (r > 0.7) continue
            emitInstance(jx, jy, (hash01(jx, jy, 41) - 0.5) * 0.7, chooseType(jx, jy, hash01(jx, jy, 43)), r)
        }
    }
}

// ---------------------------------------------------------------- 写 GLB（含定向复查）
function writeGlb(path, bucketList) {
    const json = {
        asset: { version: '2.0', generator: 'cockpit-city-massing-v2' },
        scene: 0,
        scenes: [{ nodes: [0] }],
        nodes: [{ mesh: 0, name: 'yuanling-town-massing' }],
        meshes: [{ name: 'town', primitives: [] }],
        materials: [],
        accessors: [],
        bufferViews: [],
        buffers: [],
    }
    const bin = []
    let byteOffset = 0
    const pushView = (data, target) => {
        const buf = Buffer.from(data.buffer, data.byteOffset, data.byteLength)
        json.bufferViews.push({ buffer: 0, byteOffset, byteLength: buf.length, target })
        bin.push(buf)
        byteOffset += buf.length + ((4 - (buf.length % 4)) % 4)
        return json.bufferViews.length - 1
    }
    for (const [color, b] of bucketList) {
        if (!b.indices.length) continue
        const positions = new Float32Array(b.positions)
        const normals = new Float32Array(b.normals)
        const indices = new Uint32Array(b.indices)
        const min = [Infinity, Infinity, Infinity]
        const max = [-Infinity, -Infinity, -Infinity]
        for (let i = 0; i < positions.length; i += 3) {
            for (let k = 0; k < 3; k++) {
                min[k] = Math.min(min[k], positions[i + k])
                max[k] = Math.max(max[k], positions[i + k])
            }
        }
        const posView = pushView(positions, 34962)
        const norView = pushView(normals, 34962)
        const idxView = pushView(indices, 34963)
        const posAcc = json.accessors.push({ bufferView: posView, componentType: 5126, count: positions.length / 3, type: 'VEC3', min, max }) - 1
        const norAcc = json.accessors.push({ bufferView: norView, componentType: 5126, count: normals.length / 3, type: 'VEC3' }) - 1
        const idxAcc = json.accessors.push({ bufferView: idxView, componentType: 5125, count: indices.length, type: 'SCALAR' }) - 1
        const material = json.materials.push({
            pbrMetallicRoughness: {
                baseColorFactor: [((color >> 16) & 255) / 255, ((color >> 8) & 255) / 255, (color & 255) / 255, 1],
                metallicFactor: 0,
                roughnessFactor: 0.85,
            },
            doubleSided: true, // 票 35：配合定向 pass，杜绝任何"漏面"
        }) - 1
        json.meshes[0].primitives.push({ attributes: { POSITION: posAcc, NORMAL: norAcc }, indices: idxAcc, material })
    }
    const binBuffer = Buffer.concat(bin)
    json.buffers.push({ byteLength: binBuffer.length })
    const jsonText = Buffer.from(JSON.stringify(json), 'utf8')
    const jsonChunk = Buffer.concat([jsonText, Buffer.alloc((4 - (jsonText.length % 4)) % 4, 0x20)])
    const binChunk = Buffer.concat([binBuffer, Buffer.alloc((4 - (binBuffer.length % 4)) % 4, 0)])
    const total = 12 + 8 + jsonChunk.length + 8 + binChunk.length
    const header = Buffer.alloc(12)
    header.write('glTF', 0, 'ascii')
    header.writeUInt32LE(2, 4)
    header.writeUInt32LE(total, 8)
    const jsonHeader = Buffer.alloc(8)
    jsonHeader.writeUInt32LE(jsonChunk.length, 0)
    jsonHeader.write('JSON', 4, 'ascii')
    const binHeader = Buffer.alloc(8)
    binHeader.writeUInt32LE(binChunk.length, 0)
    binHeader.write('BIN\0', 4, 'ascii')
    writeFileSync(path, Buffer.concat([header, jsonHeader, jsonChunk, binHeader, binChunk]))
    return { bytes: total, primitives: json.meshes[0].primitives.length, materials: json.materials.length }
}

// ---------------------------------------------------------------- tileset + meta
const R = 6378137
const F = 1 / 298.257223563
const E2 = 2 * F - F * F
function enuTransform(lon, lat, height) {
    const lonR = lon * DEG, latR = lat * DEG
    const sinLat = Math.sin(latR), cosLat = Math.cos(latR)
    const sinLon = Math.sin(lonR), cosLon = Math.cos(lonR)
    const N = R / Math.sqrt(1 - E2 * sinLat * sinLat)
    const ox = (N + height) * cosLat * cosLon
    const oy = (N + height) * cosLat * sinLon
    const oz = (N * (1 - E2) + height) * sinLat
    return [-sinLon, cosLon, 0, 0, -sinLat * cosLon, -sinLat * sinLon, cosLat, 0, cosLat * cosLon, cosLat * sinLon, sinLat, 0, ox, oy, oz, 1]
}

mkdirSync(DATA_DIR, { recursive: true })
const bucketList = [...buckets.entries()]
const glbInfo = writeGlb(join(DATA_DIR, 'town.glb'), bucketList)

// 包围盒
let bx = [-Infinity, Infinity], by = [-Infinity, Infinity], bz = [-Infinity, Infinity]
for (const [, b] of bucketList) {
    for (let i = 0; i < b.positions.length; i += 3) {
        const east = b.positions[i], up = b.positions[i + 1], north = -b.positions[i + 2]
        bx = [Math.max(bx[0], east), Math.min(bx[1], east)]
        by = [Math.max(by[0], north), Math.min(by[1], north)]
        bz = [Math.max(bz[0], up), Math.min(bz[1], up)]
    }
}
const box = [
    (bx[0] + bx[1]) / 2, (by[0] + by[1]) / 2, (bz[0] + bz[1]) / 2,
    (bx[0] - bx[1]) / 2, 0, 0,
    0, (by[0] - by[1]) / 2, 0,
    0, 0, (bz[0] - bz[1]) / 2,
]
const tileset = {
    asset: { version: '1.1', generator: 'cockpit-city-massing-v2' },
    geometricError: 800, // 必须 > 0：票 34 实测 0 时 Cesium 不请求根瓦片内容
    root: {
        transform: enuTransform(ORIGIN.lon, ORIGIN.lat, 0),
        boundingVolume: { box },
        geometricError: 0,
        refine: 'ADD',
        content: { uri: 'town.glb' },
    },
}
writeFileSync(join(DATA_DIR, 'tileset.json'), JSON.stringify(tileset, null, 2), 'utf8')

const meta = {
    _disclaimer: '示意体块，不是实测建筑。除 OSM 少量真实轮廓外，其余体块按真实路网/建成区程序化排布；不得用于量测、统计或应急决策。',
    _v2_changes: {
        尺寸依据: '74 个真实轮廓分位数（长边中位 53m/p90 83m，短边 27-45m），单元=真实尺寸×0.45 体块系数',
        贴地: '高程网格 64×64 + 逐栋足印四角采样，基准=min-4m 裙边（楼体向下延伸，不悬空）',
        漏顶: '所有三角面程序化定向（法线取远离体中心方向）+ 材质 doubleSided',
        中近景: '楼层色带（每层外凸 0.15m 深色带）+ 平顶女儿墙 / 坡屋顶出檐',
    },
    origin: ORIGIN,
    window: WINDOW,
    counts: { total: placed, ...stats, cap: CAP },
    archetypes: Object.fromEntries(Object.entries(ARCHETYPES).map(([k, v]) => [k, { long: v.long, short: v.short, floors: v.floors }])),
    glbBytes: glbInfo.bytes,
    localBox: { x: [bx[1], bx[0]], y: [by[1], by[0]], z: [bz[1], bz[0]] },
    // 回验样本：让 verify 脚本能独立复核"逐栋贴地"（重新按足印四角对高程网格插值）
    verificationSample: verificationSample.slice(0, 300),
    generatedAt: new Date().toISOString(),
}
writeFileSync(join(DATA_DIR, 'meta.json'), JSON.stringify(meta, null, 2), 'utf8')

console.log(`生成完成：${placed} 栋（真实轮廓 ${stats.real} / 民居 ${stats.house} / 条楼 ${stats.block} / 板楼 ${stats.slab}，避让水系 ${stats.skippedWater}）`)
console.log(`GLB：${(glbInfo.bytes / 1024 / 1024).toFixed(2)} MB，${glbInfo.primitives} 个图元 / ${glbInfo.materials} 种材质`)
console.log(`局部范围(米)：x ${bx[1].toFixed(0)}~${bx[0].toFixed(0)}，y ${by[1].toFixed(0)}~${by[0].toFixed(0)}，z ${bz[1].toFixed(0)}~${bz[0].toFixed(0)}`)