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

/**
 * 沅陵轻量体块的独立回验（票 34）。
 *
 * 这套几何是**程序化生成**的，所以"看起来对"不够 —— 这里验的是三件可证伪的事：
 *   1. **地理配准**：取一个真实 OSM 轮廓的质心，按生成时的同一条投影链推出它在 GLB 里的应有位置，
 *      再在 GLB 顶点里找得到它（容差 1m）→ 证明"经纬度 → 局部 ENU → glTF Y-up"这条链没写反；
 *   2. **tileset transform 自洽**：把 GLB 包围盒中心经 transform 送到经纬度，应落在目标城区窗口内；
 *   3. **GLB 结构合法**：header/JSON/BIN 长度自洽、accessor 与 bufferView 对得上、无 NaN。
 *
 * 用法：node tools/city-massing/verify-town.mjs
 */
import { readFileSync } from 'node:fs'
import { join } from 'node:path'

const DATA_DIR = process.env.CITY_DATA_DIR || 'D:/GIS/GIS/cityMassing'
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 OSM = 'D:/GIS/GIS/planet_110.059,28.258_110.744,28.625-shp/shpOutput'

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 magic = buf.toString('ascii', 0, 4)
const version = buf.readUInt32LE(4)
const totalLen = buf.readUInt32LE(8)
check('GLB-header', magic === 'glTF' && version === 2 && totalLen === buf.length, `${magic} v${version} len=${totalLen}/${buf.length}`)

const jsonLen = buf.readUInt32LE(12)
const jsonType = buf.toString('ascii', 16, 20)
const gltf = JSON.parse(buf.toString('utf8', 20, 20 + jsonLen))
check('GLB-json-chunk', jsonType === 'JSON' && !!gltf.meshes, `${jsonType}, primitives=${gltf.meshes[0].primitives.length}`)

const binStart = 20 + jsonLen + 8
const binChunkLen = buf.readUInt32LE(20 + jsonLen)
check('GLB-bin-chunk', binStart + binChunkLen <= buf.length, `bin=${binChunkLen}`)

function readAccessor(index) {
    const acc = gltf.accessors[index]
    const view = gltf.bufferViews[acc.bufferView]
    const offset = binStart + (view.byteOffset || 0)
    const out = []
    if (acc.componentType === 5126) {
        for (let i = 0; i < acc.count * 3; i++) out.push(buf.readFloatLE(offset + i * 4))
    } else if (acc.componentType === 5125) {
        for (let i = 0; i < acc.count; i++) out.push(buf.readUInt32LE(offset + i * 4))
    }
    return out
}

// ---------------------------------------------------------------- 3. 结构 + NaN
let vertexCount = 0, triCount = 0, nan = 0
for (const prim of gltf.meshes[0].primitives) {
    const pos = readAccessor(prim.attributes.POSITION)
    const idx = readAccessor(prim.indices)
    vertexCount += pos.length / 3
    triCount += idx.length / 3
    for (const v of pos) if (!Number.isFinite(v)) nan += 1
    const acc = gltf.accessors[prim.attributes.POSITION]
    check(
        `POSITION-accessor-minmax-${prim.material}`,
        Array.isArray(acc.min) && Array.isArray(acc.max) && acc.min.every((v, i) => v <= acc.max[i]),
        `min=[${acc.min.map((v) => v.toFixed(0))}] max=[${acc.max.map((v) => v.toFixed(0))}]`,
    )
}
check('geometry-sane', vertexCount > 1000 && triCount > 1000 && nan === 0, `顶点 ${vertexCount}，三角 ${triCount}，NaN ${nan}`)

// ---------------------------------------------------------------- 1. 真实轮廓回验
// 生成器把体块**居中**放在轮廓外接盒中心上，所以顶点落在**外接盒的四角**（真实轮廓不做旋转）。
// 因此判据是「GLB 里存在一个顶点落在该轮廓的外接盒角点上」——这条能证伪整条坐标链
// （经纬度 → 局部 ENU → glTF Y-up），写反任何一步都对不上。
// 同时按生成器的规则剔掉贴水的轮廓（那些本来就被避让掉了）。
const buildings = JSON.parse(readFileSync(join(OSM, 'buildings.json'), 'utf8')).features
const waterways = JSON.parse(readFileSync(join(OSM, 'waterways.json'), 'utf8')).features
const inWindow = (lon, lat) => Math.abs(lon - ORIGIN.lon) <= 0.02 && Math.abs(lat - ORIGIN.lat) <= 0.015
const toLocal = (lon, lat) => ({ x: (lon - ORIGIN.lon) * M_PER_DEG_LON, y: (lat - ORIGIN.lat) * M_PER_DEG_LAT })

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 = ((p.x - a.x) * dx + (p.y - a.y) * dy) / len2
            t = Math.max(0, Math.min(1, t))
            best = Math.min(best, Math.hypot(p.x - (a.x + t * dx), p.y - (a.y + t * dy)))
        }
    }
    return best
}

const targets = []
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
    if (distanceToWater({ x: cx, y: cy }) < 22) continue // 生成器会避让贴水的
    targets.push({ corners: [[minX, minY], [maxX, minY], [minX, maxY], [maxX, maxY]], area: (maxX - minX) * (maxY - minY) })
}
targets.sort((a, b) => b.area - a.area)
const sample = targets.slice(0, 8)

// 把所有顶点按 ENU 还原（glTF Y-up → ENU：east=x, up=y, north=-z）
const allVerts = []
for (const prim of gltf.meshes[0].primitives) {
    const pos = readAccessor(prim.attributes.POSITION)
    for (let i = 0; i < pos.length; i += 3) {
        allVerts.push({ x: pos[i], z: -pos[i + 2] }) // east, north（忽略 up）
    }
}
let hit = 0
const details = []
for (const t of sample) {
    const found = t.corners.some(([cx2, cy2]) => allVerts.some((v) => Math.abs(v.x - cx2) <= 2 && Math.abs(v.z - cy2) <= 2))
    if (found) hit += 1
    details.push(`${t.area.toFixed(0)}m²@(${t.corners[0][0].toFixed(0)},${t.corners[0][1].toFixed(0)})${found ? '✓' : '✗'}`)
}
check(
    'real-footprints-land-on-true-coords',
    sample.length >= 5 && hit === sample.length,
    `${hit}/${sample.length} 个最大真实轮廓的外接盒角点在 GLB 里 2m 内命中：${details.join(' ')}`,
)

// ---------------------------------------------------------------- 2. transform 自洽
const T = JSON.parse(readFileSync(join(DATA_DIR, 'tileset.json'), 'utf8')).root.transform
const box = JSON.parse(readFileSync(join(DATA_DIR, 'tileset.json'), 'utf8')).root.boundingVolume.box
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 }
}
// 原点自身（局部 0,0,0）应回到目标坐标
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（${originLL.lon.toFixed(6)}, ${originLL.lat.toFixed(6)}）`)

// 包围盒中心经 transform 后应落在窗口内
const corner = ecefToLonLat(applyTransform([box[0], box[1], box[2]]))
check(
    'bbox-center-in-window',
    Math.abs(corner.lon - ORIGIN.lon) <= 0.025 && Math.abs(corner.lat - ORIGIN.lat) <= 0.02,
    `包围盒中心 → ${corner.lon.toFixed(5)}, ${corner.lat.toFixed(5)}`,
)

const meta = JSON.parse(readFileSync(join(DATA_DIR, 'meta.json'), 'utf8'))
check('disclaimer-present', typeof meta._disclaimer === 'string' && meta._disclaimer.includes('不是实测建筑'), meta._disclaimer.slice(0, 40) + '…')

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