City-Lifeline-Cockpit/tools/city-massing/verify-roof-caps.mjs

// Verify REAL decoded roof triangles against seated GIS footprints, not material names.
// node tools/city-massing/verify-roof-caps.mjs <batch...> [--report <new-report.json>]
import assert from 'node:assert/strict'
import { readFileSync, writeFileSync } from 'node:fs'
import { resolve, join } from 'node:path'
import draco3d from 'draco3d'

const args = process.argv.slice(2)
const reportAt = args.indexOf('--report')
const reportPath = reportAt < 0 ? null : args.splice(reportAt, 2)[1]
assert.ok(args.length, 'one or more batch directories are required')
const draco = await draco3d.createDecoderModule({})
const report = { passed: false, batches: [] }
const area = (ring) => Math.abs(ring.reduce((sum, p, i) => {
    const q = ring[(i + 1) % ring.length]
    return sum + p[0] * q[1] - q[0] * p[1]
}, 0) / 2)
const perimeter = (ring) => ring.reduce((sum, p, i) => {
    const q = ring[(i + 1) % ring.length]
    return sum + Math.hypot(p[0] - q[0], p[1] - q[1])
}, 0)
const inside = (x, y, ring) => {
    let yes = false
    for (let i = 0, j = ring.length - 1; i < ring.length; j = i++) {
        const a = ring[i], b = ring[j]
        if ((a[1] > y) !== (b[1] > y) && x < (b[0] - a[0]) * (y - a[1]) / (b[1] - a[1]) + a[0]) yes = !yes
    }
    return yes
}
function readGlb(file) {
    const data = readFileSync(file)
    assert.equal(data.toString('ascii', 0, 4), 'glTF', file)
    const jsonLength = data.readUInt32LE(12)
    const json = JSON.parse(data.toString('utf8', 20, 20 + jsonLength))
    const offset = 20 + jsonLength
    assert.equal(data.toString('ascii', offset + 4, offset + 8), 'BIN\0', file)
    return { json, bin: data.subarray(offset + 8) }
}
function roofTriangles(file) {
    const { json, bin } = readGlb(file)
    const result = []
    let reversed = 0
    const node = json.nodes.find((value) => value.mesh !== undefined)
    assert.ok(node && !node.matrix && !node.rotation && !node.scale, `${file}: unexpected model transform`)
    const translation = node.translation || [0, 0, 0]
    for (const part of json.meshes[node.mesh].primitives) {
        if (!json.materials[part.material]?.name?.startsWith('Roof_')) continue
        const ext = part.extensions.KHR_draco_mesh_compression
        const view = json.bufferViews[ext.bufferView]
        const encoded = bin.subarray(view.byteOffset || 0, (view.byteOffset || 0) + view.byteLength)
        const decoder = new draco.Decoder(), buffer = new draco.DecoderBuffer(), mesh = new draco.Mesh()
        const allocated = []
        let status
        try {
            buffer.Init(new Int8Array(encoded.buffer, encoded.byteOffset, encoded.byteLength), encoded.byteLength)
            status = decoder.DecodeBufferToMesh(buffer, mesh)
            assert.ok(status.ok(), `${file}: Draco decode failed`)
            const values = (semantic) => {
                const attribute = decoder.GetAttributeByUniqueId(mesh, ext.attributes[semantic])
                const bytes = mesh.num_points() * 3 * 4
                const ptr = draco._malloc(bytes); allocated.push(ptr)
                decoder.GetAttributeDataArrayForAllPoints(mesh, attribute, draco.DT_FLOAT32, bytes, ptr)
                return new Float32Array(draco.HEAPF32.buffer, ptr, mesh.num_points() * 3).slice()
            }
            const positions = values('POSITION'), normals = values('NORMAL')
            const bytes = mesh.num_faces() * 3 * 4
            const ptr = draco._malloc(bytes); allocated.push(ptr)
            decoder.GetTrianglesUInt32Array(mesh, bytes, ptr)
            const indices = new Uint32Array(draco.HEAPU32.buffer, ptr, mesh.num_faces() * 3)
            for (let i = 0; i < indices.length; i += 3) {
                const points = [indices[i], indices[i + 1], indices[i + 2]].map((index) => [
                    positions[index * 3] + translation[0], -(positions[index * 3 + 2] + translation[2]),
                    positions[index * 3 + 1] + translation[1]])
                const [a, b, c] = points
                const projected = ((b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0])) / 2
                // Up-facing roof area in source ENU; below-ground downward faces aren't a roof.
                if (projected <= 1e-7) { if (projected < -1e-7) reversed++; continue }
                const flat = Math.max(a[2], b[2], c[2]) - Math.min(a[2], b[2], c[2]) < 0.001
                // Roof material also covers library fascia. Check vertical normals only after
                // identifying a triangle as a GIS cap, rather than judging every library detail.
                result.push({ x: (a[0] + b[0] + c[0]) / 3, y: (a[1] + b[1] + c[1]) / 3,
                    z: (a[2] + b[2] + c[2]) / 3, area: projected, flat, points,
                    upNormal: Math.min(...[indices[i], indices[i + 1], indices[i + 2]].map((index) => normals[index * 3 + 1])) })
            }
        } finally {
            for (const ptr of allocated) draco._free(ptr)
            if (status) draco.destroy(status)
            draco.destroy(mesh); draco.destroy(buffer); draco.destroy(decoder)
        }
    }
    return { triangles: result, reversed }
}
function checkFile(file, buildings, level, exactFootprints) {
    const { triangles } = roofTriangles(file)
    const expected = buildings.filter((building) => exactFootprints || building.kind === 'massing')
    const coverage = new Map(expected.map((building) => [building.id, 0]))
    // Heights + bounds reject other buildings cheaply; point in surveyed ring identifies the cap.
    const targets = expected.map((building) => ({ building,
        top: building.positionENU[2] + building.heightM,
        minX: Math.min(...building.footprintENU.map((point) => point[0])),
        maxX: Math.max(...building.footprintENU.map((point) => point[0])),
        minY: Math.min(...building.footprintENU.map((point) => point[1])),
        maxY: Math.max(...building.footprintENU.map((point) => point[1])),
    }))
    for (const triangle of triangles) {
        if (!triangle.flat) continue
        const matches = targets.filter(({ building, top, minX, maxX, minY, maxY }) => Math.abs(triangle.z - top) < 0.3
            && triangle.x >= minX && triangle.x <= maxX && triangle.y >= minY && triangle.y <= maxY
            && inside(triangle.x, triangle.y, building.footprintENU))
        // Distinct polygons in the same city sometimes overlap; avoid attributing one face twice.
        if (matches.length) {
            matches.sort((a, b) => Math.abs(a.top - triangle.z) - Math.abs(b.top - triangle.z))
            const id = matches[0].building.id
            assert.ok(triangle.upNormal > 0.99, `${file}: ${id} cap normal is not upwards`)
            coverage.set(id, coverage.get(id) + triangle.area)
        }
    }
    const missing = expected.filter((building) => coverage.get(building.id) <= 1e-7)
    assert.equal(missing.length, 0, `${file}: ${missing.length} buildings have no upward roof, first ${missing[0]?.id}`)
    let checkedArea = 0
    for (const building of expected) {
        const wanted = area(building.footprintENU), actual = coverage.get(building.id)
        const tolerance = Math.max(0.02, wanted * 0.015, perimeter(building.footprintENU) * 0.15)
        assert.ok(Math.abs(actual - wanted) <= tolerance,
            `${file}: ${building.id} roof coverage ${actual.toFixed(3)} differs from GIS area ${wanted.toFixed(3)}`)
        checkedArea += actual
    }
    return { level, buildings: expected.length, roofTriangles: triangles.length, roofAreaM2: checkedArea }
}
for (const batch of args) {
    const dir = resolve(batch)
    const manifest = JSON.parse(readFileSync(join(dir, 'manifest.json'), 'utf8'))
    const tree = JSON.parse(readFileSync(join(dir, 'tileset.json'), 'utf8'))
    // Older manifests may lack plan identity; the active batch name still maps to a seated source.
    const folder = dir.split(/[\\/]/).at(-1)
    const planName = manifest.provenance?.plan?.startsWith('seated-') ? manifest.provenance.plan
        : 'seated-' + folder.replace(/^lod-/, '').replace(/-(aligned|roofed)$/, '') + '.json'
    const plan = JSON.parse(readFileSync(join(dir, '..', planName), 'utf8'))
    const byId = new Map(plan.buildings.map((building) => [building.id, building]))
    const byContent = new Map(manifest.cells.map((cell) => [cell.lod0.uri, cell]))
    const result = { batch, buildings: manifest.count, checkedFiles: 0, checks: { coarse: 0, lod1: 0, lod0: 0 }, roofAreaM2: 0 }
    for (const cell of manifest.cells) {
        const buildings = cell.buildingIds.map((id) => { assert.ok(byId.has(id), `unknown GIS building ${id}`); return byId.get(id) })
        for (const level of ['lod1', 'lod0']) {
            const checked = checkFile(join(dir, cell[level].uri), buildings, level, level === 'lod1')
            result.checks[level] += checked.buildings; result.roofAreaM2 += checked.roofAreaM2; result.checkedFiles++
        }
    }
    for (const coarse of tree.root.children) {
        const buildings = coarse.children.flatMap((child) => byContent.get(child.content.uri).buildingIds).map((id) => byId.get(id))
        const checked = checkFile(join(dir, coarse.content.uri), buildings, 'coarse', true)
        result.checks.coarse += checked.buildings; result.roofAreaM2 += checked.roofAreaM2; result.checkedFiles++
    }
    report.batches.push(result)
    console.log(JSON.stringify(result))
}
report.passed = true
if (reportPath) writeFileSync(reportPath, JSON.stringify(report, null, 2) + '\n')
console.log(JSON.stringify({ passed: true, batches: report.batches.length,
    files: report.batches.reduce((sum, batch) => sum + batch.checkedFiles, 0),
    roofChecks: report.batches.reduce((sum, batch) => sum + Object.values(batch.checks).reduce((a, b) => a + b, 0), 0) }))