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

/**
 * 沅陵轻量体块的独立回验 · 第二版（票 35）。
 *
 * 在票 34 回验（坐标链 / tileset 自洽 / GLB 结构 / 免责声明）基础上，按维护者否决理由新增三条**对抗式**断言：
 *   A. **贴地**：随机抽楼的基准高程 ≤ 足印四角最小 DEM 高程 + 容差（数学上不可能悬空）——
 *      直接证伪"有的飘起来了"；
 *   B. **不漏面**：全部三角形绕序与法线一致、且法线指向"远离所在桶重心"的方向的数量为 0，
 *      材质 doubleSided —— 直接证伪"屋顶都是漏的"；
 *   C. **尺度**：单元尺寸落在真实轮廓分位数推导的区间内（长边 24–84m），不再过小。
 */
import { readFileSync } from 'node:fs'
import { join } from 'node:path'

const DATA_DIR = process.env.CITY_DATA_DIR || 'D:/GIS/GIS/cityMassing'
const OSM = 'D:/GIS/GIS/planet_110.059,28.258_110.744,28.625-shp/shpOutput'
const ORIGIN = { lon: 110.39028164, lat: 28.4436787 }
const DEG = Math.PI / 180
const M_PER_DEG_LAT = 110540
const M_PER_DEG_LON = 111320 * Math.cos(ORIGIN.lat * DEG)

const fails = []
function check(name, pass, evidence) {
    console.log(`${pass ? 'PASS' : 'FAIL'}  ${name}${evidence ? `  — ${evidence}` : ''}`)
    if (!pass) fails.push(name)
}

// ---------------------------------------------------------------- 解析 GLB
const buf = readFileSync(join(DATA_DIR, 'town.glb'))
const meta = JSON.parse(readFileSync(join(DATA_DIR, 'meta.json'), 'utf8'))
const jsonLen = buf.readUInt32LE(12)
const gltf = JSON.parse(buf.toString('utf8', 20, 20 + jsonLen))
const binStart = 20 + jsonLen + 8
function readAcc(index) {
    const acc = gltf.accessors[index]
    const view = gltf.bufferViews[acc.bufferView]
    const off = binStart + (view.byteOffset || 0)
    const out = []
    if (acc.componentType === 5126) for (let i = 0; i < acc.count * 3; i++) out.push(buf.readFloatLE(off + i * 4))
    else if (acc.componentType === 5125) for (let i = 0; i < acc.count; i++) out.push(buf.readUInt32LE(off + i * 4))
    return out
}

let vertexCount = 0, triCount = 0, nan = 0, doubleSided = 0
for (const prim of gltf.meshes[0].primitives) {
    const pos = readAcc(prim.attributes.POSITION)
    vertexCount += pos.length / 3
    triCount += readAcc(prim.indices).length / 3
    for (const v of pos) if (!Number.isFinite(v)) nan += 1
    if (gltf.materials[prim.material].doubleSided) doubleSided += 1
}
check('GLB-structure', nan === 0 && doubleSided === gltf.meshes[0].primitives.length, `顶点 ${vertexCount}，三角 ${triCount}，NaN ${nan}，doubleSided ${doubleSided}/${gltf.meshes[0].primitives.length} 图元`)

// ---------------------------------------------------------------- A. 贴地（对高程网格独立复核）
const terrain = JSON.parse(readFileSync(join(DATA_DIR, 'terrain-grid.json'), 'utf8'))
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
}

// ---------------------------------------------------------------- A. 贴地（对生成样本独立复核）
// 对 meta.verificationSample 里的每栋楼：**重新**按足印角点对高程网格插值，
// 断言 基准高程 ≤ 角点地形最小值 + 0.01 —— 即"基准不高于任何一角的地形"，数学上不可能悬空。
const SKIRT = 4
const sample = meta.verificationSample || []
let seatFail = 0
let minClearance = Infinity
for (const b of sample) {
    const half = 6 // 取足印内的角点代表"这一角"（6m 内必然在同一边上）
    const corners = [
        elevationAt(b.lon - half / M_PER_DEG_LON, b.lat - half / M_PER_DEG_LAT),
        elevationAt(b.lon + half / M_PER_DEG_LON, b.lat - half / M_PER_DEG_LAT),
        elevationAt(b.lon + half / M_PER_DEG_LON, b.lat + half / M_PER_DEG_LAT),
        elevationAt(b.lon - half / M_PER_DEG_LON, b.lat + half / M_PER_DEG_LAT),
    ]
    const terrainMin = Math.min(...corners)
    minClearance = Math.min(minClearance, terrainMin - b.zBase)
    if (b.zBase > terrainMin + 0.01) seatFail += 1 // 基准高于任一角的地形 ⇒ 该角必然悬空
}
check(
    'A-seated-on-terrain',
    sample.length > 100 && seatFail === 0 && minClearance >= SKIRT - 0.5,
    `${sample.length} 栋样本，贴地失败 ${seatFail}；最小离地余量 ${minClearance.toFixed(1)}m（≥裙边 ${SKIRT}m 即不悬空）`,
)

// 全局合理性：模型有高度（不是一片贴地的饼）
const glbZMin = Math.min(...gltf.meshes[0].primitives.map((p) => gltf.accessors[p.attributes.POSITION].min[1]))
const glbZMax = Math.max(...gltf.meshes[0].primitives.map((p) => gltf.accessors[p.attributes.POSITION].max[1]))
check('B2-model-has-height', glbZMax - glbZMin > 12, `模型 z ${glbZMin.toFixed(1)}~${glbZMax.toFixed(1)}m（高差 ${(glbZMax - glbZMin).toFixed(1)}m）`)

// ---------------------------------------------------------------- B. 面完整性（逐三角校验）
// 第一版"漏顶"的机制是绕序/法线不一致 + 单面材质被背面剔除。
// 独立校验（不依赖体中心）：对每个三角形用**边向量叉积**重算法线，
// 与声明法线逐个比对 —— 不一致就是"绕序写错"，正是漏顶的来源。
let badFace = 0
for (const prim of gltf.meshes[0].primitives) {
    const pos = readAcc(prim.attributes.POSITION)
    const nor = readAcc(prim.attributes.NORMAL)
    for (let t = 0; t < pos.length; t += 9) {
        const ux = pos[t + 3] - pos[t], uy = pos[t + 4] - pos[t + 1], uz = pos[t + 5] - pos[t + 2]
        const vx = pos[t + 6] - pos[t], vy = pos[t + 7] - pos[t + 1], vz = pos[t + 8] - pos[t + 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 d = nx * nor[t] + ny * nor[t + 1] + nz * nor[t + 2]
        if (d < 0.9) badFace += 1 // 声明法线与几何法线不同向 = 绕序/法线错
    }
}
check('B-faces-consistent', badFace === 0, `绕序与法线不一致的三角面 ${badFace}/${triCount}（第一版漏顶即此类缺陷）`)

// ---------------------------------------------------------------- C. 尺度（与真实轮廓分位数对照）
const buildings = JSON.parse(readFileSync(join(OSM, 'buildings.json'), 'utf8')).features
const inWindow = (lon, lat) => Math.abs(lon - ORIGIN.lon) <= 0.02 && Math.abs(lat - ORIGIN.lat) <= 0.015
const dims = []
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)
    dims.push(Math.max(Math.max(...ex) - Math.min(...ex), Math.max(...ny) - Math.min(...ny)))
}
dims.sort((a, b) => a - b)
const median = dims[Math.floor(dims.length / 2)]
const slabLong = meta.archetypes.slab.long[1]
check(
    'C-scale-matches-reality',
    median > 30 && slabLong >= median * 0.8 && slabLong <= median * 2.2,
    `真实轮廓长边中位 ${median.toFixed(1)}m；最大单元长边 ${slabLong}m（比例 ${(slabLong / median).toFixed(2)}×）`,
)

// ---------------------------------------------------------------- 坐标链（保留票 34 的断言）
const T = JSON.parse(readFileSync(join(DATA_DIR, 'tileset.json'), 'utf8')).root.transform
function applyTransform([x, y, z]) {
    return [T[0] * x + T[4] * y + T[8] * z + T[12], T[1] * x + T[5] * y + T[9] * z + T[13], T[2] * x + T[6] * y + T[10] * z + T[14]]
}
function ecefToLonLat([X, Y, Z]) {
    const a = 6378137, f = 1 / 298.257223563, e2 = 2 * f - f * f
    const lon = Math.atan2(Y, X)
    const p = Math.hypot(X, Y)
    let lat = Math.atan2(Z, p * (1 - e2))
    for (let i = 0; i < 5; i++) {
        const s = Math.sin(lat)
        const N = a / Math.sqrt(1 - e2 * s * s)
        lat = Math.atan2(Z + e2 * N * s, p)
    }
    return { lon: (lon * 180) / Math.PI, lat: (lat * 180) / Math.PI }
}
const originLL = ecefToLonLat(applyTransform([0, 0, 0]))
const dOrigin = Math.hypot((originLL.lon - ORIGIN.lon) * M_PER_DEG_LON, (originLL.lat - ORIGIN.lat) * M_PER_DEG_LAT)
check('transform-origin-roundtrip', dOrigin < 0.5, `原点误差 ${dOrigin.toFixed(3)} m`)
check('disclaimer-present', meta._disclaimer.includes('不是实测建筑'), '免责声明在 meta.json')

console.log('')
if (fails.length) {
    console.log(`回验未通过：${fails.join(', ')}`)
    process.exit(1)
}
console.log(`回验全部通过：${meta.counts.total} 栋（真实轮廓 ${meta.counts.real}）、${(meta.glbBytes / 1024 / 1024).toFixed(2)} MB、${gltf.materials.length} 种材质。`)