City-Lifeline-Cockpit/tools/city-massing/segment-imagery-footprints.mjs

// Segment the local high-resolution satellite imagery inside the confirmed residential masks into
// building footprints, so generated city massing sits on satellite-visible roofs instead of being
// inferred from roads.
//
// Deterministic: no randomness, fixed iteration order, rounded numbers.
// Usage: node tools/city-massing/segment-imagery-footprints.mjs [maskId]
//   maskId is an optional debug filter that limits the run to one residential mask.
import { readFileSync, writeFileSync, existsSync } from 'node:fs'
import { resolve } from 'node:path'
import { PNG } from 'pngjs'
import jpeg from 'jpeg-js'
import * as turf from '@turf/turf'
import { isWithinPrimaryArea } from './primary-area.mjs'

const IMAGERY_DIR = process.env.GIS_DIR ? `${process.env.GIS_DIR}/yuanlingPic3` : 'D:/GIS/GIS/yuanlingPic3'
const CITY_DIR = resolve('tools/city-massing/city')
const OUT_GEOJSON = resolve(CITY_DIR, 'imagery-footprints.geojson')
const OUT_REPORT = resolve(CITY_DIR, 'imagery-footprints-report.json')
const OUT_AUDIT = resolve(CITY_DIR, 'imagery-footprints-audit.png')

const TILE = 256
const Z = 19 // segmentation zoom, about 0.26 m/px at this latitude
const Z_AUDIT = 18 // audit backdrop zoom, about 0.53 m/px (exactly half of Z)
const AUDIT_MAX_SIDE = 2400 // panel cap: below this a 10 m building is still several pixels wide
const MASK_FILTER = process.argv[2] ? Number(process.argv[2]) : null

// Every threshold the run uses, echoed verbatim into the report.
const THRESHOLDS = Object.freeze({
    greenMin: 8, // vegetation index g - max(r, b) must exceed this ...
    greenG: 55, // ... and g must exceed this
    waterBrightness: 90, // very dark ...
    waterBlueGreen: 14, // ... strongly blue-green ...
    waterTexture: 9, // ... and low-texture (5x5 brightness std below this)
    darkLightSplit: 95, // brightness separating the dark-roof and light-roof families
    edgeBarrier: 14, // Sobel gradient on blurred brightness above this blocks connectivity
    erodeRadius: 1, // px eroded so thin shadow gaps between neighbours split them
    minAreaSqm: 40,
    maxAreaSqm: 4000,
    minFillRatio: 0.55, // component area / minimum-area rotated rectangle area
    maxLongSideM: 90,
    minShortSideM: 4,
    ringRadiusPx: 3,
    maxRingVegetation: 0.5, // reject when the dilated ring is dominated by vegetation
})

const SOURCE = 'yuanlingPic3 z19 TMS local imagery, per-pixel classified then 8-connected segmented'

// ---------------------------------------------------------------- projections
const scale = 2 ** Z
const lonToPx = (lon) => ((lon + 180) / 360) * scale * TILE
const latToPy = (lat) => ((1 - Math.asinh(Math.tan((lat * Math.PI) / 180)) / Math.PI) / 2) * scale * TILE
const pxToLon = (x) => (x / (scale * TILE)) * 360 - 180
const pyToLat = (y) => (Math.atan(Math.sinh(Math.PI * (1 - (2 * y) / (scale * TILE)))) * 180) / Math.PI
/** Ground resolution in metres per pixel at a latitude. */
const metresPerPixel = (lat) => (156543.03392 * Math.cos((lat * Math.PI) / 180)) / scale

/** z/x/googleY tile range covering a WGS84 bbox at a zoom. */
function tileRange(bbox, zoom) {
    const n = 2 ** zoom
    const xOf = (lon) => ((lon + 180) / 360) * n
    const yOf = (lat) => ((1 - Math.asinh(Math.tan((lat * Math.PI) / 180)) / Math.PI) / 2) * n
    return { x0: Math.floor(xOf(bbox[0])), x1: Math.floor(xOf(bbox[2])), y0: Math.floor(yOf(bbox[3])), y1: Math.floor(yOf(bbox[1])) }
}

/** Read one TMS tile at `{z}/{x}/{tmsY}.png`; JPEG bytes when the first byte is 0xFF. */
function readTile(zoom, x, googleY) {
    const tmsY = 2 ** zoom - 1 - googleY
    const path = `${IMAGERY_DIR}/${zoom}/${x}/${tmsY}.png`
    if (!existsSync(path)) return null
    const bytes = readFileSync(path)
    return { image: bytes[0] === 0xff ? jpeg.decode(bytes) : PNG.sync.read(bytes), ref: [zoom, x, tmsY] }
}

const round = (value, digits) => Number(value.toFixed(digits))

// ------------------------------------------------------- classification raster
/** Copy one decoded tile's RGB into a window-local plane, 3 bytes per pixel. */
function blitTile(image, plane, winW, winH, ox, oy) {
    for (let y = 0; y < TILE; y++) {
        const wy = oy + y
        if (wy < 0 || wy >= winH) continue
        for (let x = 0; x < TILE; x++) {
            const wx = ox + x
            if (wx < 0 || wx >= winW) continue
            const k = (y * TILE + x) * 4
            const d = (wy * winW + wx) * 3
            const a = image.data[k + 3]
            plane[d] = a === 0 ? 0 : image.data[k]
            plane[d + 1] = a === 0 ? 0 : image.data[k + 1]
            plane[d + 2] = a === 0 ? 0 : image.data[k + 2]
        }
    }
}

/** Separable 3x3 box blur, borders clamped so window edges are not artificial features. */
function boxBlur3(src, winW, winH) {
    const temp = new Uint8Array(winW * winH)
    for (let y = 0; y < winH; y++) {
        const row = y * winW
        for (let x = 0; x < winW; x++) {
            const a = src[row + (x > 0 ? x - 1 : 0)]
            const b = src[row + x]
            const c = src[row + (x < winW - 1 ? x + 1 : winW - 1)]
            temp[row + x] = (a + b + c) / 3
        }
    }
    const out = new Uint8Array(winW * winH)
    for (let y = 0; y < winH; y++) {
        const row = y * winW
        const up = (y > 0 ? y - 1 : 0) * winW
        const down = (y < winH - 1 ? y + 1 : winH - 1) * winW
        for (let x = 0; x < winW; x++) out[row + x] = (temp[up + x] + temp[row + x] + temp[down + x]) / 3
    }
    return out
}

/**
 * Per-pixel classification into 0 barrier, 1 vegetation, 2 water, 3 dark roof, 4 light roof.
 * The barrier class is an edge test, so two roofs separated by a bright gutter or a dark shadow
 * line do not fuse into one component even when both sides fall in the same brightness family.
 */
function classify(rgb, blur, winW, winH) {
    const value = new Uint8Array(winW * winH)
    for (let y = 0; y < winH; y++) {
        for (let x = 0; x < winW; x++) {
            const i = y * winW + x
            const k = i * 3
            const r = rgb[k]
            const g = rgb[k + 1]
            const b = rgb[k + 2]
            const v = blur[i]
            if (g - Math.max(r, b) > THRESHOLDS.greenMin && g > THRESHOLDS.greenG) {
                value[i] = 1
            } else if (v < THRESHOLDS.waterBrightness && (g + b) / 2 - r > THRESHOLDS.waterBlueGreen && std5(blur, winW, winH, x, y) < THRESHOLDS.waterTexture) {
                value[i] = 2
            } else if (sobelMagnitude(blur, winW, winH, x, y) > THRESHOLDS.edgeBarrier) {
                value[i] = 0
            } else {
                value[i] = v < THRESHOLDS.darkLightSplit ? 3 : 4
            }
        }
    }
    return value
}

/** 5x5 standard deviation of brightness; water is smooth, roofs and vegetation are not. */
function std5(plane, winW, winH, x, y) {
    let sum = 0
    let sumSq = 0
    let n = 0
    for (let dy = -2; dy <= 2; dy++) {
        const yy = y + dy
        if (yy < 0 || yy >= winH) continue
        for (let dx = -2; dx <= 2; dx++) {
            const xx = x + dx
            if (xx < 0 || xx >= winW) continue
            const v = plane[yy * winW + xx]
            sum += v
            sumSq += v * v
            n++
        }
    }
    const mean = sum / n
    return Math.sqrt(Math.max(0, sumSq / n - mean * mean))
}

/** Sobel gradient magnitude of the blurred brightness. */
function sobelMagnitude(plane, winW, winH, x, y) {
    const at = (xx, yy) => plane[(yy < 0 ? 0 : yy >= winH ? winH - 1 : yy) * winW + (xx < 0 ? 0 : xx >= winW ? winW - 1 : xx)]
    const tl = at(x - 1, y - 1)
    const tc = at(x, y - 1)
    const tr = at(x + 1, y - 1)
    const ml = at(x - 1, y)
    const mr = at(x + 1, y)
    const bl = at(x - 1, y + 1)
    const bc = at(x, y + 1)
    const br = at(x + 1, y + 1)
    const gx = tr + 2 * mr + br - tl - 2 * ml - bl
    const gy = bl + 2 * bc + br - tl - 2 * tc - tr
    return Math.hypot(gx, gy) / 4
}

// ------------------------------------------------------------- polygon raster
/** Even-odd scanline fill of a WGS84 polygon with holes onto the window pixel grid. */
function rasterizePolygon(rings, px0, py0, winW, winH) {
    const raster = new Uint8Array(winW * winH)
    const edges = []
    for (const ring of rings) {
        for (let i = 1; i < ring.length; i++) {
            edges.push([lonToPx(ring[i - 1][0]) - px0, latToPy(ring[i - 1][1]) - py0, lonToPx(ring[i][0]) - px0, latToPy(ring[i][1]) - py0])
        }
    }
    const xs = []
    for (let y = 0; y < winH; y++) {
        const cy = y + 0.5
        xs.length = 0
        for (const [ax, ay, bx, by] of edges) {
            if (ay === by) continue
            if (cy < Math.min(ay, by) || cy >= Math.max(ay, by)) continue
            xs.push(ax + ((cy - ay) / (by - ay)) * (bx - ax))
        }
        if (xs.length < 2) continue
        xs.sort((a, b) => a - b)
        for (let i = 0; i + 1 < xs.length; i += 2) {
            const from = Math.max(0, Math.ceil(xs[i] - 0.5))
            const to = Math.min(winW - 1, Math.floor(xs[i + 1] - 0.5))
            for (let x = from; x <= to; x++) raster[y * winW + x] = 1
        }
    }
    return raster
}

// ------------------------------------------------------- connected components
/** Shrink each roof family by `radius` px so thin shadow gaps split neighbouring buildings. */
function erode(value, winW, winH, radius) {
    const out = new Uint8Array(winW * winH)
    for (let y = radius; y < winH - radius; y++) {
        for (let x = radius; x < winW - radius; x++) {
            const i = y * winW + x
            const c = value[i]
            if (c < 3) continue
            let keep = 1
            for (let dy = -radius; dy <= radius && keep; dy++) {
                const row = (y + dy) * winW
                for (let dx = -radius; dx <= radius; dx++) {
                    if (value[row + x + dx] !== c) { keep = 0; break }
                }
            }
            out[i] = keep ? c : 0
        }
    }
    return out
}

/** Two-pass 8-connected labelling; returns one label per pixel (0 = background). */
function labelComponents(mask, winW, winH) {
    const labels = new Int32Array(winW * winH)
    const parent = [0]
    const find = (v) => {
        while (parent[v] !== v) { parent[v] = parent[parent[v]]; v = parent[v] }
        return v
    }
    const union = (a, b) => { a = find(a); b = find(b); if (a !== b) parent[Math.max(a, b)] = Math.min(a, b) }
    let next = 1
    for (let y = 0; y < winH; y++) {
        for (let x = 0; x < winW; x++) {
            const i = y * winW + x
            if (!mask[i]) continue
            let best = 0
            if (y > 0) {
                const up = i - winW
                for (let dx = -1; dx <= 1; dx++) {
                    const xx = x + dx
                    if (xx < 0 || xx >= winW) continue
                    const l = labels[up + dx]
                    if (l && (best === 0 || l < best)) best = l
                }
            }
            if (x > 0 && labels[i - 1] && (best === 0 || labels[i - 1] < best)) best = labels[i - 1]
            if (best === 0) { best = next; parent[next] = next; next++ }
            labels[i] = best
            if (y > 0) {
                const up = i - winW
                for (let dx = -1; dx <= 1; dx++) {
                    const xx = x + dx
                    if (xx < 0 || xx >= winW) continue
                    const l = labels[up + dx]
                    if (l && l !== best) union(l, best)
                }
            }
            if (x > 0 && labels[i - 1] && labels[i - 1] !== best) union(labels[i - 1], best)
        }
    }
    const remap = new Int32Array(next).fill(-1)
    let count = 0
    for (let i = 0; i < labels.length; i++) {
        if (!labels[i]) continue
        const root = find(labels[i])
        if (remap[root] < 0) remap[root] = ++count
        labels[i] = remap[root]
    }
    return { labels, count }
}

/**
 * Grow each label back out by `radius` px, but only onto pixels of its own family, so the erode
 * step does not permanently shrink footprints and neighbours cannot fuse back together.
 */
function dilateLabels(labels, value, winW, winH, radius) {
    const out = Int32Array.from(labels)
    for (let y = 0; y < winH; y++) {
        for (let x = 0; x < winW; x++) {
            const i = y * winW + x
            if (labels[i] || value[i] < 3) continue
            let best = 0
            for (let dy = -radius; dy <= radius; dy++) {
                const yy = y + dy
                if (yy < 0 || yy >= winH) continue
                const row = yy * winW
                for (let dx = -radius; dx <= radius; dx++) {
                    const xx = x + dx
                    if (xx < 0 || xx >= winW) continue
                    const l = labels[row + xx]
                    if (l && value[row + xx] === value[i] && (best === 0 || l < best)) best = l
                }
            }
            out[i] = best
        }
    }
    return out
}

// ------------------------------------------------------------- geometry maths
/** Convex hull (monotone chain) of pixel centres. */
function convexHull(points) {
    const pts = points.slice().sort((a, b) => a[0] - b[0] || a[1] - b[1])
    if (pts.length < 3) return pts
    const cross = (o, a, b) => (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0])
    const lower = []
    for (const p of pts) {
        while (lower.length >= 2 && cross(lower[lower.length - 2], lower[lower.length - 1], p) <= 0) lower.pop()
        lower.push(p)
    }
    const upper = []
    for (let i = pts.length - 1; i >= 0; i--) {
        const p = pts[i]
        while (upper.length >= 2 && cross(upper[upper.length - 2], upper[upper.length - 1], p) <= 0) upper.pop()
        upper.push(p)
    }
    lower.pop()
    upper.pop()
    return lower.concat(upper)
}

/** Minimum-area rotated rectangle around a pixel set, by rotating calipers over hull edges. */
function minAreaRect(points) {
    const hull = convexHull(points)
    if (hull.length < 3) return null
    let best = null
    for (let i = 0; i < hull.length; i++) {
        const a = hull[i]
        const b = hull[(i + 1) % hull.length]
        const dx = b[0] - a[0]
        const dy = b[1] - a[1]
        const len = Math.hypot(dx, dy)
        if (len === 0) continue
        const ux = dx / len
        const uy = dy / len
        let minU = Infinity
        let maxU = -Infinity
        let minV = Infinity
        let maxV = -Infinity
        for (const p of hull) {
            const ox = p[0] - a[0]
            const oy = p[1] - a[1]
            const u = ox * ux + oy * uy
            const v = -ox * uy + oy * ux
            if (u < minU) minU = u
            if (u > maxU) maxU = u
            if (v < minV) minV = v
            if (v > maxV) maxV = v
        }
        const area = (maxU - minU) * (maxV - minV)
        if (!best || area < best.area) {
            best = { area, ux, uy, ax: a[0], ay: a[1], minU, maxU, minV, maxV, lenU: maxU - minU, lenV: maxV - minV }
        }
    }
    return best
}

/** The four corners of a rotated rectangle, in pixel space, counter-clockwise. */
const rectCorners = (rect) =>
    [[rect.minU, rect.minV], [rect.maxU, rect.minV], [rect.maxU, rect.maxV], [rect.minU, rect.maxV]].map(([u, v]) => [
        rect.ax + u * rect.ux - v * rect.uy,
        rect.ay + u * rect.uy + v * rect.ux,
    ])

// ------------------------------------------------------------- exclusion layers
function loadLayer(name) {
    return JSON.parse(readFileSync(resolve(CITY_DIR, name), 'utf8')).features.map((feature) => ({ feature, bbox: turf.bbox(feature) }))
}
const bboxOverlaps = (a, b) => !(a[2] < b[0] || a[0] > b[2] || a[3] < b[1] || a[1] > b[3])
/** First layer entry the polygon touches, or null. Layers are bbox-prefiltered for speed. */
function firstIntersection(polygon, bbox, layer) {
    for (const entry of layer) {
        if (!bboxOverlaps(bbox, entry.bbox)) continue
        if (turf.booleanIntersects(polygon, entry.feature)) return entry
    }
    return null
}

// ------------------------------------------------------------------- pipeline
const allMasks = JSON.parse(readFileSync(resolve(CITY_DIR, 'residential-clipped.geojson'), 'utf8')).features
const masks = allMasks
    .filter((mask) => MASK_FILTER === null || mask.properties.ORIG_FID === MASK_FILTER)
    .sort((a, b) => a.properties.ORIG_FID - b.properties.ORIG_FID)
if (!masks.length) throw new Error(`No residential mask matched ${MASK_FILTER}`)
const roads = loadLayer('roads-clipped.geojson')
const waters = loadLayer('water-clipped.geojson')
const osmBuildings = loadLayer('osm-buildings-clipped.geojson')
const regularRoofs = loadLayer('regular-estate-roofs.geojson')

const REJECTION_ORDER = [
    'areaBelowMin',
    'areaAboveMax',
    'fillRatioBelowMin',
    'longSideAboveMax',
    'shortSideBelowMin',
    'ringDominatedByVegetation',
    'outsideMask',
    'outsidePrimaryArea',
    'intersectsRoad',
    'intersectsWater',
    'intersectsOsmBuilding',
    'overlapsRegularEstateRoof',
]

const features = []
const auditPanels = []
const rejectionTotals = Object.fromEntries(REJECTION_ORDER.map((reason) => [reason, 0]))
const perMask = []
let tilesConsidered = 0
let tilesRead = 0
let tilesMissing = 0

/**
 * Segment one residential mask. Kept as a function so every window-sized buffer goes out of scope
 * before the next mask starts; the largest masks need several gigabytes at z19.
 */
function segmentMask(mask) {
    const maskId = mask.properties.ORIG_FID
    const rings = mask.geometry.type === 'Polygon' ? mask.geometry.coordinates : mask.geometry.coordinates.flat()
    const bbox = turf.bbox(mask)
    const range = tileRange(bbox, Z)
    // The window is the mask bbox snapped out to whole tiles, so the ring test has context off-mask.
    const px0 = range.x0 * TILE
    const py0 = range.y0 * TILE
    const winW = (range.x1 - range.x0 + 1) * TILE
    const winH = (range.y1 - range.y0 + 1) * TILE
    const pixelCount = winW * winH
    const mpp = metresPerPixel((bbox[1] + bbox[3]) / 2)

    const considered = (range.x1 - range.x0 + 1) * (range.y1 - range.y0 + 1)
    tilesConsidered += considered
    let read = 0
    let missing = 0

    const rgb = new Uint8Array(pixelCount * 3)
    for (let ty = range.y0; ty <= range.y1; ty++) {
        for (let tx = range.x0; tx <= range.x1; tx++) {
            const tile = readTile(Z, tx, ty)
            if (!tile) { missing++; continue }
            read++
            blitTile(tile.image, rgb, winW, winH, tx * TILE - px0, ty * TILE - py0)
        }
    }
    tilesRead += read
    tilesMissing += missing

    const brightness = new Uint8Array(pixelCount)
    for (let i = 0; i < pixelCount; i++) brightness[i] = (rgb[i * 3] + rgb[i * 3 + 1] + rgb[i * 3 + 2]) / 3
    const blur = boxBlur3(brightness, winW, winH)
    const value = classify(rgb, blur, winW, winH)

    // Candidates only exist inside the mask; off-mask pixels stay as context for the ring test.
    const inMask = rasterizePolygon(rings, px0, py0, winW, winH)
    for (let i = 0; i < pixelCount; i++) if (value[i] >= 3 && !inMask[i]) value[i] = 0

    const eroded = erode(value, winW, winH, THRESHOLDS.erodeRadius)
    const { labels, count } = labelComponents(eroded, winW, winH)
    const grown = dilateLabels(labels, value, winW, winH, THRESHOLDS.erodeRadius)

    // Bucket pixels by component so each component is one contiguous slice of `index`.
    const sizes = new Int32Array(count + 2)
    for (let i = 0; i < pixelCount; i++) if (grown[i]) sizes[grown[i]]++
    const start = new Int32Array(count + 2)
    for (let l = 1; l <= count; l++) start[l + 1] = start[l] + sizes[l]
    const cursor = start.slice()
    const index = new Int32Array(start[count + 1])
    for (let i = 0; i < pixelCount; i++) {
        const l = grown[i]
        if (l) index[cursor[l]++] = i
    }

    const maskFeatures = []
    const maskRejections = Object.fromEntries(REJECTION_ORDER.map((reason) => [reason, 0]))
    const rejectedRects = []
    const member = new Uint8Array(pixelCount)
    const ringFlag = new Uint8Array(pixelCount)
    const touched = []

    for (let label = 1; label <= count; label++) {
        const from = start[label]
        const to = start[label + 1]
        if (from === to) continue
        const points = []
        let brightnessSum = 0
        let family = 0
        for (let i = from; i < to; i++) {
            const p = index[i]
            member[p] = 1
            points.push([(p % winW) + 0.5, ((p / winW) | 0) + 0.5])
            brightnessSum += blur[p]
            family |= value[p] === 3 ? 1 : 2
        }

        const areaSqm = points.length * mpp * mpp
        let reason = null
        if (areaSqm < THRESHOLDS.minAreaSqm) reason = 'areaBelowMin'
        else if (areaSqm > THRESHOLDS.maxAreaSqm) reason = 'areaAboveMax'

        let rect = null
        let fillRatio = 0
        let longSideM = 0
        let shortSideM = 0
        if (!reason) {
            rect = minAreaRect(points)
            const rectArea = rect ? rect.area : 0
            fillRatio = rectArea > 0 ? points.length / rectArea : 0
            longSideM = rect ? Math.max(rect.lenU, rect.lenV) * mpp : 0
            shortSideM = rect ? Math.min(rect.lenU, rect.lenV) * mpp : 0
            if (fillRatio < THRESHOLDS.minFillRatio) reason = 'fillRatioBelowMin'
            else if (longSideM > THRESHOLDS.maxLongSideM) reason = 'longSideAboveMax'
            else if (shortSideM < THRESHOLDS.minShortSideM) reason = 'shortSideBelowMin'
        }

        // Ring test: a candidate ringed by vegetation is a field or a courtyard, not a roof.
        if (!reason) {
            let ringPixels = 0
            let vegetation = 0
            const R = THRESHOLDS.ringRadiusPx
            touched.length = 0
            for (let i = from; i < to; i++) {
                const p = index[i]
                const x = p % winW
                const y = (p / winW) | 0
                for (let dy = -R; dy <= R; dy++) {
                    const yy = y + dy
                    if (yy < 0 || yy >= winH) continue
                    const row = yy * winW
                    for (let dx = -R; dx <= R; dx++) {
                        const xx = x + dx
                        if (xx < 0 || xx >= winW) continue
                        const q = row + xx
                        if (member[q] || ringFlag[q]) continue
                        ringFlag[q] = 1
                        touched.push(q)
                        ringPixels++
                        if (value[q] === 1) vegetation++
                    }
                }
            }
            for (const q of touched) ringFlag[q] = 0
            if (ringPixels > 0 && vegetation / ringPixels > THRESHOLDS.maxRingVegetation) reason = 'ringDominatedByVegetation'
        }

        let ring = null
        if (!reason) {
            const windowCorners = rectCorners(rect)
            ring = windowCorners.map(([x, y]) => [pxToLon(px0 + x), pyToLat(py0 + y)])
            ring.push(ring[0])
            const polygon = turf.polygon([ring])
            const cheapInside = windowCorners.every(([x, y]) => {
                const xi = Math.round(x - 0.5)
                const yi = Math.round(y - 0.5)
                return xi >= 0 && yi >= 0 && xi < winW && yi < winH && inMask[yi * winW + xi] === 1
            })
            const polygonBbox = turf.bbox(polygon)
            if (!cheapInside || !turf.booleanWithin(polygon, mask)) reason = 'outsideMask'
            else if (!isWithinPrimaryArea(polygon)) reason = 'outsidePrimaryArea'
            else if (firstIntersection(polygon, polygonBbox, roads)) reason = 'intersectsRoad'
            else if (firstIntersection(polygon, polygonBbox, waters)) reason = 'intersectsWater'
            else if (firstIntersection(polygon, polygonBbox, osmBuildings)) reason = 'intersectsOsmBuilding'
            else if (firstIntersection(polygon, polygonBbox, regularRoofs)) reason = 'overlapsRegularEstateRoof'
        }

        if (reason) {
            maskRejections[reason]++
            rejectionTotals[reason]++
            // Only near misses are drawn on the audit: components too small to have a rectangle are
            // noise and would bury the accepted outlines.
            if (rect) rejectedRects.push(rectCorners(rect).map(([x, y]) => [px0 + x, py0 + y]))
        } else {
            // Heading of the long side, degrees clockwise from north, normalised to [0, 180).
            const width = Math.hypot(ring[1][0] - ring[0][0], ring[1][1] - ring[0][1])
            const height = Math.hypot(ring[2][0] - ring[1][0], ring[2][1] - ring[1][1])
            const longIsWidth = width >= height
            const a = longIsWidth ? ring[0] : ring[1]
            const b = longIsWidth ? ring[1] : ring[2]
            const east = (b[0] - a[0]) * Math.cos((((a[1] + b[1]) / 2) * Math.PI) / 180)
            const north = -(b[1] - a[1])
            const headingDeg = ((((Math.atan2(east, north) * 180) / Math.PI) % 180) + 180) % 180
            const tileRefs = new Set()
            for (let i = from; i < to; i++) {
                const p = index[i]
                tileRefs.add(`${Z}/${Math.floor((px0 + (p % winW)) / TILE)}/${Math.floor((py0 + ((p / winW) | 0)) / TILE)}`)
            }
            maskFeatures.push({
                type: 'Feature',
                properties: {
                    id: `img-${maskId}-${maskFeatures.length + 1}`,
                    maskId,
                    source: 'imagery-segmentation',
                    areaSqm: round(areaSqm, 1),
                    longSideM: round(longSideM, 1),
                    shortSideM: round(shortSideM, 1),
                    headingDeg: round(headingDeg, 1),
                    fillRatio: round(fillRatio, 3),
                    meanBrightness: round(brightnessSum / points.length, 1),
                    tileRefs: [...tileRefs].sort().map((ref) => ref.split('/').map(Number)),
                    roofFamily: family === 1 ? 'dark' : family === 2 ? 'light' : 'mixed',
                },
                geometry: { type: 'Polygon', coordinates: [ring] },
            })
        }
        for (let i = from; i < to; i++) member[index[i]] = 0
    }

    console.log(`mask ${maskId}: ${count} candidate components, ${maskFeatures.length} accepted (${read} tiles read, ${missing} missing)`)
    return {
        maskFeatures,
        auditPanel: { maskId, bbox, features: maskFeatures, rejected: rejectedRects },
        summary: {
            maskId,
            maskAreaSqm: Math.round(turf.area(mask)),
            tileRange: { z: Z, x0: range.x0, x1: range.x1, googleY0: range.y0, googleY1: range.y1 },
            tilesConsidered: considered,
            tilesRead: read,
            tilesMissing: missing,
            windowPx: [winW, winH],
            metresPerPixel: round(mpp, 4),
            candidateComponents: count,
            accepted: maskFeatures.length,
            rejected: REJECTION_ORDER.reduce((sum, key) => sum + maskRejections[key], 0),
            rejections: maskRejections,
        },
    }
}

for (const mask of masks) {
    const { maskFeatures, auditPanel, summary } = segmentMask(mask)
    features.push(...maskFeatures)
    auditPanels.push(auditPanel)
    perMask.push(summary)
}

// --------------------------------------------------------------- assertions
// Re-test every accepted footprint against all six exclusions, independently of the loop above.
const violations = []
for (const feature of features) {
    const polygon = turf.feature(feature.geometry)
    const polygonBbox = turf.bbox(polygon)
    const mask = masks.find((m) => m.properties.ORIG_FID === feature.properties.maskId)
    if (!turf.booleanWithin(polygon, mask)) violations.push([feature.properties.id, 'outsideMask'])
    if (!isWithinPrimaryArea(polygon)) violations.push([feature.properties.id, 'outsidePrimaryArea'])
    for (const [name, layer] of [['road', roads], ['water', waters], ['osmBuilding', osmBuildings], ['regularEstateRoof', regularRoofs]]) {
        if (firstIntersection(polygon, polygonBbox, layer)) violations.push([feature.properties.id, `intersects${name[0].toUpperCase()}${name.slice(1)}`])
    }
}
console.log(`exclusion assertions on ${features.length} accepted footprints: ${violations.length} violations`)
if (violations.length) throw new Error(`Accepted footprints violate exclusions: ${JSON.stringify(violations.slice(0, 10))}`)

// ------------------------------------------------------------------- outputs
const sortedReasons = Object.entries(rejectionTotals).sort(
    (a, b) => b[1] - a[1] || REJECTION_ORDER.indexOf(a[0]) - REJECTION_ORDER.indexOf(b[0])
)
writeFileSync(OUT_GEOJSON, JSON.stringify({ type: 'FeatureCollection', features }, null, 2))

const audit = renderAudit(auditPanels)
writeFileSync(OUT_AUDIT, PNG.sync.write(audit.image, { deflateLevel: 6 }))

writeFileSync(
    OUT_REPORT,
    JSON.stringify(
        {
            source: SOURCE,
            imageryDir: IMAGERY_DIR,
            command: 'node tools/city-massing/segment-imagery-footprints.mjs',
            thresholds: THRESHOLDS,
            tiles: {
                z19: { considered: tilesConsidered, read: tilesRead, missing: tilesMissing },
                z18Audit: { read: audit.tilesRead, missing: audit.tilesMissing },
            },
            tileRefFormat: 'tileRefs entries are [z, x, tmsY]; tmsY = 2**z - 1 - googleY matches the on-disk path',
            perMask,
            totals: {
                accepted: features.length,
                candidateComponents: perMask.reduce((sum, m) => sum + m.candidateComponents, 0),
                rejected: perMask.reduce((sum, m) => sum + m.rejected, 0),
            },
            rejectionReasons: Object.fromEntries(sortedReasons),
            exclusionAssertions: {
                checked: features.length,
                violations: violations.length,
                tests: ['fullyInsideResidentialMask', 'fullyInsidePrimaryArea', 'noRoadIntersection', 'noWaterIntersection', 'noOsmBuildingIntersection', 'noRegularEstateRoofOverlap'],
            },
            auditImage: {
                path: 'tools/city-massing/city/imagery-footprints-audit.png',
                layout: `per-mask z18 crops stacked vertically; accepted footprints 2 px bright yellow, rejected components 1 px red; panels top to bottom: ${auditPanels.map((p) => `mask ${p.maskId}`).join(', ')}`,
                size: [audit.image.width, audit.image.height],
            },
            limitations: [
                'This is imagery segmentation, not survey geometry. Boundaries follow visible roof pixels at about 0.26 m/px.',
                'Two touching buildings with no shadow or brightness edge between them merge into a single footprint.',
                'Very dark roofs on dark ground are missed or merged with the ground, because brightness cannot separate them.',
                'Dark blue roofs can be classified as water, because they are dark, blue-dominant and smooth.',
                'The emitted polygon is a minimum-area rotated rectangle: it does not follow L-shapes, courtyards or roof ridges.',
                'Bare bright ground (construction spoil, concrete yards) can be classified as light roof; the fill, side and ring filters remove most but not all of it.',
                'Positional uncertainty is on the order of one to two pixels, roughly 0.3-0.6 m.',
                'No height, storey count or roof material is derived; only the footprint outline.',
            ],
        },
        null,
        2
    )
)

console.log(
    JSON.stringify({
        command: 'node tools/city-massing/segment-imagery-footprints.mjs',
        accepted: features.length,
        perMask: perMask.map((m) => [m.maskId, m.accepted]),
        rejectionReasons: Object.fromEntries(sortedReasons),
        tilesRead,
        tilesMissing,
        audit: [audit.image.width, audit.image.height],
    })
)

// ---------------------------------------------------------------- audit image
function renderAudit(panels) {
    const rendered = panels.map(renderPanel)
    const width = Math.max(...rendered.map((panel) => panel.width))
    const height = rendered.reduce((sum, panel) => sum + panel.height, 0)
    const image = new PNG({ width, height })
    image.data.fill(255)
    let offset = 0
    for (const panel of rendered) {
        for (let y = 0; y < panel.height; y++) {
            const from = y * panel.width * 4
            image.data.set(panel.data.subarray(from, from + panel.width * 4), (offset + y) * width * 4)
        }
        offset += panel.height
    }
    return {
        image,
        tilesRead: rendered.reduce((sum, panel) => sum + panel.tilesRead, 0),
        tilesMissing: rendered.reduce((sum, panel) => sum + panel.tilesMissing, 0),
    }
}

/** One mask panel: z18 imagery backdrop, accepted outlines bright yellow, rejected components red. */
function renderPanel(panel) {
    const range = tileRange(panel.bbox, Z_AUDIT)
    const x0 = range.x0 * TILE
    const y0 = range.y0 * TILE
    const cropW = (range.x1 - range.x0 + 1) * TILE
    const cropH = (range.y1 - range.y0 + 1) * TILE
    const factor = Math.min(1, AUDIT_MAX_SIDE / cropW, AUDIT_MAX_SIDE / cropH)
    const width = Math.max(1, Math.round(cropW * factor))
    const height = Math.max(1, Math.round(cropH * factor))
    const image = new PNG({ width, height })
    image.data.fill(255)
    const cache = new Map()
    let tilesRead = 0
    let tilesMissing = 0
    for (let py = 0; py < height; py++) {
        const cy = Math.min(cropH - 1, Math.floor(py / factor))
        for (let px = 0; px < width; px++) {
            const cx = Math.min(cropW - 1, Math.floor(px / factor))
            const gx = x0 + cx
            const gy = y0 + cy
            const tx = Math.floor(gx / TILE)
            const ty = Math.floor(gy / TILE)
            const key = tx * 1e6 + ty
            let tile = cache.get(key)
            if (tile === undefined) {
                tile = readTile(Z_AUDIT, tx, ty)
                if (tile) tilesRead++
                else tilesMissing++
                cache.set(key, tile)
            }
            const d = (py * width + px) * 4
            if (tile) {
                const s = ((gy - ty * TILE) * TILE + (gx - tx * TILE)) * 4
                image.data[d] = tile.image.data[s]
                image.data[d + 1] = tile.image.data[s + 1]
                image.data[d + 2] = tile.image.data[s + 2]
            }
            image.data[d + 3] = 255
        }
    }
    const project = (x, y) => [(x - x0) * factor, (y - y0) * factor]
    const stroke = (a, b, color, half) => {
        const steps = Math.max(1, Math.ceil(Math.hypot(b[0] - a[0], b[1] - a[1])))
        for (let i = 0; i <= steps; i++) {
            const x = Math.round(a[0] + ((b[0] - a[0]) * i) / steps)
            const y = Math.round(a[1] + ((b[1] - a[1]) * i) / steps)
            for (let dy = -half; dy <= half; dy++) {
                for (let dx = -half; dx <= half; dx++) {
                    const xx = x + dx
                    const yy = y + dy
                    if (xx < 0 || yy < 0 || xx >= width || yy >= height) continue
                    const k = (yy * width + xx) * 4
                    image.data[k] = color[0]
                    image.data[k + 1] = color[1]
                    image.data[k + 2] = color[2]
                }
            }
        }
    }
    // Rectangles and polygon corners are held in z19 window pixels; the backdrop is at Z_AUDIT.
    const zoomShift = 2 ** (Z - Z_AUDIT)
    for (const corners of panel.rejected) {
        const pts = corners.map(([x, y]) => project(x / zoomShift, y / zoomShift))
        for (let i = 0; i < 4; i++) stroke(pts[i], pts[(i + 1) % 4], [255, 0, 0], 0)
    }
    for (const feature of panel.features) {
        const pts = feature.geometry.coordinates[0].slice(0, 4).map(([lon, lat]) => project(lonToPx(lon) / zoomShift, latToPy(lat) / zoomShift))
        for (let i = 0; i < 4; i++) stroke(pts[i], pts[(i + 1) % 4], [255, 255, 0], 1)
    }
    return { width, height, data: image.data, tilesRead, tilesMissing }
}