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

// Verify decoded GLB transforms/feet (not self-reported zBase) against original outline + direct DEM samples.
import assert from 'node:assert/strict'
import { readFileSync, writeFileSync } from 'node:fs'
import { resolve } from 'node:path'
import { Cartesian3, Cartographic, Matrix4, Quaternion, Math as CM } from 'cesium'
const out = resolve(process.env.CITY_PREVIEW_DIR || 'tools/city-massing/district/generated')
const file = readFileSync(resolve(out, 'district.glb'))
assert.equal(file.toString('ascii', 0, 4), 'glTF')
assert.equal(file.readUInt32LE(4), 2)
assert.equal(file.readUInt32LE(8), file.length)
const jsonLen = file.readUInt32LE(12),
    gltf = JSON.parse(file.toString('utf8', 20, 20 + jsonLen)),
    binStart = 28 + jsonLen
const plan = JSON.parse(readFileSync('tools/city-massing/district/seated-plan.json'))
const samples = JSON.parse(readFileSync('tools/city-massing/district/terrain-samples.json'))
const library = JSON.parse(readFileSync('tools/city-massing/house-library/library.json'))
const report = JSON.parse(readFileSync(resolve(out, 'placement-report.json')))
const tile = JSON.parse(readFileSync(resolve(out, 'tileset.json')))
const root = Matrix4.fromArray(tile.root.transform)
function nodeMatrices(doc) {
    const map = new Map()
    function visit(i, parent) {
        const n = doc.nodes[i]
        const local = n.matrix
            ? Matrix4.fromArray(n.matrix)
            : Matrix4.fromTranslationQuaternionRotationScale(
                  Cartesian3.fromArray(n.translation || [0, 0, 0]),
                  new Quaternion(...(n.rotation || [0, 0, 0, 1])),
                  Cartesian3.fromArray(n.scale || [1, 1, 1])
              )
        const world = Matrix4.multiply(parent, local, new Matrix4())
        map.set(n.name, { node: n, world })
        for (const c of n.children || []) visit(c, world)
    }
    for (const i of doc.scenes[doc.scene || 0].nodes) visit(i, Matrix4.IDENTITY)
    return map
}
function positions(primitive) {
    const a = gltf.accessors[primitive.attributes.POSITION],
        v = gltf.bufferViews[a.bufferView],
        out = []
    assert.equal(a.componentType, 5126)
    assert.equal(a.type, 'VEC3')
    for (let i = 0; i < a.count; i++) {
        const o = binStart + (v.byteOffset || 0) + (a.byteOffset || 0) + i * (v.byteStride || 12)
        out.push(new Cartesian3(file.readFloatLE(o), file.readFloatLE(o + 4), file.readFloatLE(o + 8)))
    }
    return out
}
function enu(p) {
    return new Cartesian3(p.x, -p.z, p.y)
}
function fromNode(point, m) {
    return enu(Matrix4.multiplyByPoint(m, point, new Cartesian3()))
}
function errors(doc) {
    const failures = [],
        mats = nodeMatrices(doc)
    const metrics = []
    for (const b of report.buildings) {
        const actual = mats.get(b.id),
            foot = mats.get(b.id + '_buried_foundation')
        if (!actual || !foot) {
            failures.push(b.id + ' missing node')
            continue
        }
        const src = plan.buildings.find((x) => x.id === b.id)
        const entry = library.archetypes.find((x) => x.id === b.asset)
        const zero = fromNode(new Cartesian3(0, 0, 0), actual.world)
        const world = Matrix4.multiplyByPoint(root, zero, new Cartesian3())
        const ll = Cartographic.fromCartesian(world)
        // Independent absolute ECEF control using WGS84 converter, not degree-to-metre approximations.
        const expected = Cartesian3.fromDegrees(src.center.lon, src.center.lat, src.floorHeight)
        const displacement = Cartesian3.distance(expected, world)
        if (displacement > 0.03) failures.push(b.id + ' centre/height shifted ' + displacement)
        const nominal = entry.nominal_footprint_m,
            w = nominal[0],
            d = nominal[1]
        const corners = [
            [-w / 2, 0, d / 2],
            [w / 2, 0, d / 2],
            [w / 2, 0, -d / 2],
            [-w / 2, 0, -d / 2],
        ].map((c) => fromNode(Cartesian3.fromArray(c), actual.world))
        const cornerError = Math.max(
            ...corners.map((p) => Math.min(...src.footprintENU.map((q) => Math.hypot(p.x - q[0], p.y - q[1]))))
        )
        if (cornerError > 0.05) failures.push(b.id + ' footprint size/orientation mismatch ' + cornerError)
        const points = doc.meshes[actual.node.mesh].primitives.flatMap(positions).map((p) => fromNode(p, actual.world))
        const bottoms = doc.meshes[foot.node.mesh].primitives.flatMap(positions).map((p) => fromNode(p, foot.world))
        const minFoot = Math.min(...bottoms.map((p) => p.z)),
            maxFoot = Math.max(...bottoms.map((p) => p.z)),
            bodyMin = Math.min(...points.map((p) => p.z)),
            roofMax = Math.max(...points.map((p) => p.z))
        // Recompute the original DEM controls instead of trusting seated-plan's groundControls/zBase.
        const inverse = Matrix4.inverseTransformation(root, new Matrix4())
        const controls = src.sampleIndices.map((i) => {
            const coordinate = samples.planPoints[i]
            const height = samples.heights[i]
            assert(Number.isFinite(height), 'Missing original DEM height')
            return Matrix4.multiplyByPoint(inverse, Cartesian3.fromDegrees(coordinate.lon, coordinate.lat, height), new Cartesian3())
        })
        const minGround = Math.min(...controls.map((p) => p.z)),
            maxGround = Math.max(...controls.map((p) => p.z))
        if (minFoot > minGround - 2.8 || maxFoot < maxGround - 0.4)
            failures.push(b.id + ' foundation no longer intersects sampled ground')
        if (bodyMin > maxGround + 0.1 || bodyMin < maxGround - 0.5) failures.push(b.id + ' body floor floated or sank')
        if (roofMax - maxGround < b.actualFloors * 2.95) failures.push(b.id + ' usable facade/roof buried')
        if (points.some((p) => !Number.isFinite(p.x + p.y + p.z))) failures.push(b.id + ' NaN geometry')
        metrics.push({
            id: b.id,
            sourceId: b.sourceId,
            displacement,
            cornerError,
            bodyMin,
            roofMax,
            minGround,
            maxGround,
            foundationMin: minFoot,
            foundationMax: maxFoot,
            visibleFoundation: bodyMin - minGround,
            actualFloors: b.actualFloors,
        })
    }
    return { failures, metrics }
}
const baseline = errors(gltf)
console.log({
    nodes: report.count,
    geoFailures: baseline.failures.slice(0, 5),
    maxCornerError: Math.max(...baseline.metrics.map((x) => x.cornerError)),
})
assert.equal(baseline.failures.length, 0, baseline.failures.join('\n'))
// Negative controls modify actual final GLB scene data, keep all source/metadata untouched.
const first = report.buildings[0].id
for (const mode of ['horizontal', 'vertical', 'scale']) {
    let d = structuredClone(gltf)
    let n = d.nodes.find((n) => n.name === first)
    n.translation ||= [0, 0, 0]
    n.scale ||= [1, 1, 1]
    if (mode === 'horizontal') n.translation[0] += 10
    if (mode === 'vertical') n.translation[1] += 5
    if (mode === 'scale') n.scale[0] *= 1.2
    const r = errors(d)
    assert(r.failures.length > 0, 'Mutated ' + mode + ' should fail')
    console.log('negative control detected:', mode)
}
writeFileSync(
    resolve(out, 'geometry-verification.json'),
    JSON.stringify(
        {
            passed: true,
            sourceFootprints: plan.buildings.length,
            placed: report.count,
            directTerrainSamples: samples.heights.length,
            metrics: baseline.metrics,
            negativeControls: ['horizontal shift 10m', 'height shift 5m', 'X scale 1.2'],
        },
        null,
        2
    )
)
console.log('PASS actual GLB anchors, dimensions, terrain intersection, facade clearance, and mutation controls')