City-Lifeline-Cockpit/tools/city-massing/seat-cad-bim.mjs

// Seat every CAD BIM building on the original DEM.
// Flat ground is seated directly; a gently sloped footprint gets a bounded stone plinth; steeper
// footprints are refused rather than hidden behind a fabricated platform.
import { readFileSync, writeFileSync, existsSync } from 'node:fs'
import { resolve } from 'node:path'
import { Cartesian3, Matrix4 } from 'cesium'
import { localSeatingLevels } from './terrain-height-frame.mjs'

// Existing published assets keep their original coordinates. New candidates opt
// into the corrected frame; never silently reseat the accepted city in place.
const useLocalHeightFrame = process.argv.includes('--local-height-frame')

const dir = 'tools/city-massing/city/cad-bim/'
const option = name => { const at = process.argv.indexOf(name); return at < 0 ? null : process.argv[at + 1] }
const output = option('--output') || dir + 'seated-plan.json'
if (useLocalHeightFrame && (!option('--plan') || !option('--samples') || !option('--output') ||
    resolve(output) === resolve(dir + 'seated-plan.json') || existsSync(output))) {
    throw new Error('local height frame requires explicit candidate --plan, --samples and separate --output')
}
const plan = JSON.parse(readFileSync(option('--plan') || dir + 'plan.json'))
const sampled = JSON.parse(readFileSync(option('--samples') || dir + 'terrain-samples.json'))
if (plan.samples.length !== sampled.heights.length) throw new Error('DEM sample count does not match the plan')

const FLAT_RELIEF = 2.5
const MAX_PLINTH = 3.2
const PLINTH_HEIGHT_RATIO = 0.4
const FLOOR_HEIGHT_M = 3.2
const MAX_HEIGHT_M = 30

const inverse = Matrix4.inverseTransformation(Matrix4.fromArray(plan.transform), new Matrix4())
const local = (lon, lat, h) => Matrix4.multiplyByPoint(inverse, Cartesian3.fromDegrees(lon, lat, h), new Cartesian3())

// Survey outlines occasionally fold back on themselves (a "bow tie"), repeat a vertex, or enclose a
// negligible area. Those cannot be extruded into a solid, so they are rejected here — with a reason —
// instead of reaching the exporter and aborting a full-city run.
const MIN_MASSING_AREA = 4
const DUPLICATE_VERTEX_M = 0.05

function cleanRing(points) {
    const out = []
    for (const point of points) {
        const previous = out[out.length - 1]
        if (previous && Math.hypot(point[0] - previous[0], point[1] - previous[1]) < DUPLICATE_VERTEX_M) continue
        out.push(point)
    }
    if (out.length > 2) {
        const first = out[0], last = out[out.length - 1]
        if (Math.hypot(first[0] - last[0], first[1] - last[1]) < DUPLICATE_VERTEX_M) out.pop()
    }
    return out
}

function signedArea(points) {
    let total = 0
    for (let i = 0; i < points.length; i++) {
        const [x1, y1] = points[i]
        const [x2, y2] = points[(i + 1) % points.length]
        total += x1 * y2 - x2 * y1
    }
    return total / 2
}

function segmentsCross(a, b) {
    const cross = (o, p, q) => (p[0] - o[0]) * (q[1] - o[1]) - (p[1] - o[1]) * (q[0] - o[0])
    const d1 = cross(a[0], a[1], b[0]), d2 = cross(a[0], a[1], b[1])
    const d3 = cross(b[0], b[1], a[0]), d4 = cross(b[0], b[1], a[1])
    return ((d1 > 0) !== (d2 > 0)) && ((d3 > 0) !== (d4 > 0))
}

/** True when no two non-adjacent edges of the ring intersect. */
function isSimpleRing(points) {
    const n = points.length
    for (let i = 0; i < n; i++) {
        const a = [points[i], points[(i + 1) % n]]
        for (let j = i + 1; j < n; j++) {
            if (j === i || (j + 1) % n === i || j === (i + 1) % n) continue
            const b = [points[j], points[(j + 1) % n]]
            if (segmentsCross(a, b)) return false
        }
    }
    return true
}

const buildings = []
const excluded = []
const construction = { flat: 0, plinth: 0 }
let missing = 0
let invalidGeometry = 0

for (const building of plan.buildings) {
    const heights = building.sampleIndices.map((i) => sampled.heights[i])
    if (heights.some((h) => !Number.isFinite(h) || h <= 1)) {
        missing++
        excluded.push({ id: building.id, reason: 'missing or invalid DEM sample' })
        continue
    }
    let ring = building.footprintENU
    if (building.kind === 'massing') {
        ring = cleanRing(ring)
        if (ring.length < 3) {
            invalidGeometry++
            excluded.push({ id: building.id, tier: building.tier, reason: 'outline collapses to fewer than 3 distinct points' })
            continue
        }
        if (!isSimpleRing(ring)) {
            invalidGeometry++
            excluded.push({ id: building.id, tier: building.tier, reason: 'self-intersecting survey outline cannot be extruded' })
            continue
        }
        if (Math.abs(signedArea(ring)) < MIN_MASSING_AREA) {
            invalidGeometry++
            excluded.push({ id: building.id, tier: building.tier, reason: 'outline encloses negligible area' })
            continue
        }
    }
    const levels = useLocalHeightFrame ? localSeatingLevels(plan.transform,
        building.sampleIndices.map(i => plan.samples[i]), heights) : null
    const min = levels ? levels.minimum : Math.min(...heights)
    const max = levels ? levels.maximum : Math.max(...heights)
    const relief = max - min
    const heightM = Math.min(MAX_HEIGHT_M, Math.max(FLOOR_HEIGHT_M, building.floors * FLOOR_HEIGHT_M))
    if (relief > FLAT_RELIEF) {
        if (relief > MAX_PLINTH || relief > heightM * PLINTH_HEIGHT_RATIO) {
            excluded.push({
                id: building.id, tier: building.tier,
                reason: 'terrain relief exceeds the allowed plinth; needs stepped construction',
                relief: Math.round(relief * 100) / 100,
            })
            continue
        }
        construction.plinth++
    } else {
        construction.flat++
    }
    // Floor sits just below the highest measured ground so the building never floats.
    const floorZ = max - 0.25
    const centroid = building.sampleIndices[building.sampleIndices.length - 1]
    const anchor = plan.samples[centroid]
    const anchorHeight = levels ? levels.ellipsoidMaximum - 0.25 : floorZ
    buildings.push({
        id: building.id,
        kind: building.kind,
        tier: building.tier,
        maskId: building.maskId,
        footprintENU: ring,
        width: building.width,
        depth: building.depth,
        heading: building.heading,
        archetype: building.archetype,
        preferredAssets: building.preferredAssets,
        floors: building.floors,
        heightM: Math.round(heightM * 100) / 100,
        heightSource: building.heightSource,
        positionENU: [
            Math.round(local(anchor.lon, anchor.lat, anchorHeight).x * 100) / 100,
            Math.round(local(anchor.lon, anchor.lat, anchorHeight).y * 100) / 100,
            Math.round(floorZ * 100) / 100,
        ],
        ...(levels ? { heightFrame: levels.heightFrame,
            terrainEllipsoidMin: levels.ellipsoidMinimum, terrainEllipsoidMax: levels.ellipsoidMaximum } : {}),
        foundationBottom: Math.round((min - 3) * 100) / 100,
        terrainMin: Math.round(min * 100) / 100,
        terrainMax: Math.round(max * 100) / 100,
        terrainRelief: Math.round(relief * 100) / 100,
        construction: relief > FLAT_RELIEF ? 'plinth' : 'flat',
    })
}

const out = {
    ...plan,
    buildings,
    excluded: [...(plan.excluded || []), ...excluded],
    construction,
    notice: plan.disclaimer + ' Seated on the original quantized-mesh DEM; steep footprints stay excluded.',
}
writeFileSync(output, JSON.stringify(out, null, 2))
console.log(JSON.stringify({
    planned: plan.buildings.length,
    seated: buildings.length,
    byKind: {
        asset: buildings.filter((b) => b.kind === 'asset').length,
        massing: buildings.filter((b) => b.kind === 'massing').length,
    },
    byConstruction: construction,
    missingSamples: missing,
    invalidGeometry,
    excludedRelief: excluded.filter((e) => e.relief !== undefined).length,
    urbanSeated: buildings.filter((b) => b.tier === 'urban').length,
    ruralSeated: buildings.filter((b) => b.tier === 'rural').length,
}))