City-Lifeline-Cockpit/tools/city-massing/assemble-cad-bim.py

"""Export a CAD/BIM seated plan as a spatial-LOD 3D Tiles tree.

Run with Blender: blender --background --factory-startup --python-exit-code 1 \
    --python tools/city-massing/assemble-cad-bim.py -- \
    --plan tools/city-massing/cad-bim/seated-plan.json \
    --output tools/city-massing/cad-bim/lod-preview --grid 32 --coarse-grid 8

Same output conventions as assemble-lod-district.py (GLB cells, `box()` bounding volumes,
tileset.json/manifest.json, root-anchor.glb), so verify-lod-district.mjs still gates it.
Difference: input footprints are surveyed polygons, so `kind == "massing"` buildings are
built here as extruded prisms instead of being copied from the house library.
"""
import json
import math
import struct
import argparse
import sys
from pathlib import Path

import bpy
import bmesh
from mathutils import Vector
from mathutils.geometry import tessellate_polygon

ROOT = Path(__file__).resolve().parents[2]
parser = argparse.ArgumentParser()
parser.add_argument('--plan', required=True)
parser.add_argument('--output', required=True)
parser.add_argument('--grid', type=int, default=32)
parser.add_argument('--coarse-grid', type=int, default=0)
parser.add_argument('--limit', type=int, default=0, help='export only the first N buildings')
parser.add_argument('--decimate-lod1', type=float, default=1.0,
                    help='Decimate ratio applied to library asset meshes at LOD1 (1.0 = off)')
parser.add_argument('--decimate-lod0', type=float, default=1.0,
                    help='Decimate ratio applied to library asset meshes at LOD0 (1.0 = off)')
parser.add_argument('--merge-buildings', action='store_true',
                    help='Bake each cell into ONE merged mesh (per-material primitives) instead of one node per building')
parser.add_argument('--lod1-footprint-only', action='store_true',
                    help='Use each surveyed GIS footprint for LOD1/coarse geometry, including asset buildings')
parser.add_argument('--near-architecture', action='store_true', help='Add bounded window reveals and storey bands at LOD0')
parser.add_argument('--massing-storey-bands', action='store_true', help='Experimental continuous storey bands on massing LOD0; requires a facade profile or --near-architecture')
parser.add_argument('--massing-facade-profile', choices=('off', 'compact', 'full'), default='off',
                    help='Opt-in facade quality for surveyed massing: compact windows/bands/eaves without deep reveals, or full near architecture')
parser.add_argument('--compact-window-texture', help='Opt-in estimated window-frame texture manifest; compact profile only')
parser.add_argument('--roof-textures', help='JSON mapping building IDs to sourced roof image paths and ENU UV bounds')
parser.add_argument('--draco', action='store_true', help='Enable Draco mesh compression in GLB export')
args = parser.parse_args(sys.argv[sys.argv.index('--') + 1:] if '--' in sys.argv else [])
for ratio in (args.decimate_lod1, args.decimate_lod0):
    if not (0.0 < ratio <= 1.0):
        raise RuntimeError('--decimate-lod* ratios must be in (0, 1]')
if args.grid < 2:
    raise RuntimeError('--grid must be at least 2')
if args.coarse_grid and (args.coarse_grid < 2 or args.coarse_grid >= args.grid):
    raise RuntimeError('--coarse-grid must be at least 2 and smaller than --grid')
PLAN = Path(args.plan).resolve()
OUT = Path(args.output).resolve()
LIBRARY_DIR = ROOT / 'tools/city-massing/house-library'
plan = json.loads(PLAN.read_text(encoding='utf-8'))
lib = json.loads((LIBRARY_DIR / 'library.json').read_text(encoding='utf-8'))
entries = lib['archetypes']
# Same asset-fit contract the placement generator used, so the two stages cannot disagree about
# what counts as an acceptable match.
ASSET_FIT = json.loads((ROOT / 'tools/city-massing/asset-fit.json').read_text(encoding='utf-8'))['variants']
# Scene variety without random re-rolls: a stable hash of the building id picks one wall tone and
# one roof tone, so neighbours differ while the whole extrusion stays coherent.
WALL_TONES = ('Plaster_Cream', 'Plaster_Ivory', 'Plaster_Sand', 'Concrete_Warm')
ROOF_TONES = ('Roof_Slate', 'Roof_Slate_Light', 'Roof_Terracotta', 'Roof_Terra_Light')
ROOF_TEXTURES = json.loads(Path(args.roof_textures).read_text(encoding='utf-8')) if args.roof_textures else {}
def facade_profile_flags(profile, near_architecture, storey_bands):
    """Resolve explicit massing quality without silently upgrading the default path."""
    if profile == 'full' and not near_architecture:
        raise RuntimeError('--massing-facade-profile full requires --near-architecture')
    if storey_bands and not (near_architecture or profile != 'off'):
        raise RuntimeError('--massing-storey-bands requires --near-architecture or a massing-facade-profile')
    return {
        'profile': profile,
        'architecture': near_architecture or profile in ('compact', 'full'),
        'fullReveals': near_architecture or profile == 'full',
        # Compact already gets the ordinary per-edge storey rhythm; the extra
        # boundary band remains an explicit opt-in for controlled experiments.
        'bands': storey_bands,
    }

PROFILE_FLAGS = facade_profile_flags(args.massing_facade_profile, args.near_architecture, args.massing_storey_bands)
FACADE_PROFILE = PROFILE_FLAGS['profile']
FULL_REVEALS = PROFILE_FLAGS['fullReveals']
ARCHITECTURE = PROFILE_FLAGS['architecture']
BANDS = PROFILE_FLAGS['bands']
WINDOW_TEXTURE = None
if args.compact_window_texture:
    if FACADE_PROFILE != 'compact' or args.near_architecture:
        raise RuntimeError('--compact-window-texture requires compact without --near-architecture')
    sys.path.insert(0, str(Path(__file__).resolve().parent))
    from window_texture import load_window_texture, window_quad_uv
    WINDOW_TEXTURE = load_window_texture(args.compact_window_texture)


def compact_scope_ids(plan, profile, full_reveals=False, bands=False):
    ids = plan.get('compactFacadeBuildingIds')
    if ids is None:
        return None
    if profile != 'compact' or full_reveals or bands:
        raise ValueError('compactFacadeBuildingIds requires compact without reveals or extra bands')
    if not isinstance(ids, list) or not ids or any(not isinstance(i, str) for i in ids) or len(ids) != len(set(ids)):
        raise ValueError('compactFacadeBuildingIds must be a nonempty unique ID list')
    records = {building['id']: building for building in plan['buildings']}
    if any(i not in records or records[i]['kind'] != 'massing' for i in ids):
        raise ValueError('compactFacadeBuildingIds must name existing massing records only')
    return frozenset(ids)


COMPACT_BUILDING_IDS = compact_scope_ids(plan, FACADE_PROFILE, FULL_REVEALS, BANDS)


def facade_operation(building, operation, *values, **options):
    """Use legacy facade rules for nonselected context; restore state even on failure."""
    global FACADE_PROFILE, ARCHITECTURE, FULL_REVEALS, BANDS, WINDOW_TEXTURE
    if COMPACT_BUILDING_IDS is None or building['id'] in COMPACT_BUILDING_IDS:
        return operation(*values, **options)
    previous = FACADE_PROFILE, ARCHITECTURE, FULL_REVEALS, BANDS, WINDOW_TEXTURE
    try:
        FACADE_PROFILE, ARCHITECTURE, FULL_REVEALS, BANDS, WINDOW_TEXTURE = 'off', False, False, False, None
        return operation(*values, **options)
    finally:
        FACADE_PROFILE, ARCHITECTURE, FULL_REVEALS, BANDS, WINDOW_TEXTURE = previous


def select_asset(building):
    family = building.get('archetype') or 'house'
    prefix = {'house': 'house', 'apartment': 'slab', 'midrise': 'midrise'}.get(family, 'house')
    # A building may name the exact approved variants it was planned for; otherwise any variant of
    # its family that fits the shared scale window is eligible.
    preferred = building.get('preferredAssets') or []
    family_key = family if family in ASSET_FIT else prefix
    limits = ASSET_FIT.get(family_key) or ASSET_FIT.get(prefix) or ASSET_FIT['house']
    min_x, min_y = limits['minScale']
    max_x, max_y = limits['maxScale']
    candidates = []
    for entry in entries:
        if preferred:
            if entry['id'] not in preferred:
                continue
        elif not entry['id'].startswith(prefix):
            continue
        width, depth = entry['nominal_footprint_m']
        sx, sy = building['width'] / width, building['depth'] / depth
        if not (min_x <= sx <= max_x and min_y <= sy <= max_y):
            continue
        preference = preferred.index(entry['id']) if entry['id'] in preferred else len(preferred) + 1
        score = preference * 0.25 + abs(entry['floors'] - building['floors']) * 0.8 + abs(math.log(sx)) + abs(math.log(sy))
        candidates.append((score, entry, sx, sy))
    return min(candidates, key=lambda item: item[0]) if candidates else None


# ---------------------------------------------------------------- polygon helpers
def ring_area(points):
    """Signed XY area; positive means counter-clockwise."""
    return 0.5 * sum(points[i][0] * points[(i + 1) % len(points)][1]
                     - points[(i + 1) % len(points)][0] * points[i][1] for i in range(len(points)))


def as_ccw(points):
    """Keep surveyed coordinates, removing consecutive duplicates at CAD centimetre precision."""
    ring = []
    for point in points:
        if not ring or math.hypot(point[0] - ring[-1][0], point[1] - ring[-1][1]) > 0.011:
            ring.append(point)
    if len(ring) > 1 and math.hypot(ring[0][0] - ring[-1][0], ring[0][1] - ring[-1][1]) <= 0.011:
        ring.pop()
    return ring if ring_area(ring) > 0 else list(reversed(ring))


def convex_hull(points):
    """Monotone chain; used as the cheap LOD1 outline for a massing prism."""
    ordered = sorted({(round(x, 6), round(y, 6)) for x, y in points})
    if len(ordered) < 3:
        return ordered

    def half(seq):
        out = []
        for point in seq:
            while len(out) >= 2 and ((out[-1][0] - out[-2][0]) * (point[1] - out[-2][1])
                                     - (out[-1][1] - out[-2][1]) * (point[0] - out[-2][0])) <= 0:
                out.pop()
            out.append(point)
        return out

    return half(ordered)[:-1] + half(reversed(ordered))[:-1]


def prism_triangles(n):
    """n-gon prism: two n-gon caps plus n wall quads."""
    return 4 * n - 4


def facade_window_count(edge_length, level):
    """Small, bounded facade detail count; LOD1 keeps one readable window per edge."""
    if level == 'lod1':
        return 1
    # Compact LOD0 uses physical-scale bays. The cap bounds the work on unusually
    # long survey edges; walls beyond 105.6m still require a separate review.
    if FACADE_PROFILE == 'compact':
        return 0 if edge_length < 2.4 else max(1, min(32, int(edge_length / 3.2)))
    # Full near architecture gets a denser, still bounded facade rhythm; coarse LOD stays sparse.
    return max(1, min(20 if FULL_REVEALS else 3, int(edge_length / (2.4 if FULL_REVEALS else 3.5))))


def facade_entrance_edge(footprint):
    """Compact ties use surveyed endpoint coordinates; legacy profiles keep their order."""
    if FACADE_PROFILE != 'compact':
        return max(range(len(footprint)), key=lambda i: math.dist(footprint[i], footprint[(i + 1) % len(footprint)]))
    return min(range(len(footprint)), key=lambda i: (
        -round(math.dist(footprint[i], footprint[(i + 1) % len(footprint)]), 6),
        tuple(sorted((tuple(footprint[i]), tuple(footprint[(i + 1) % len(footprint)]))))))


def facade_panel_halfwidth(edge_length, width_ratio=0.42, physical_window=False):
    """Shared panel sizing for generated geometry and entrance-clearance checks."""
    return (min(0.8, edge_length * 0.5) if physical_window else
            min(edge_length * 0.22, max(0.45, edge_length / 3.5 * width_ratio)))


def facade_window_positions(edge_length, level, reserve_entrance=False):
    """Keep upper-floor rhythm; reserve 15cm between compact ground windows and door."""
    count = facade_window_count(edge_length, level)
    positions = [(index + 1) / (count + 1) for index in range(count)]
    if not reserve_entrance or FACADE_PROFILE != 'compact' or not count:
        return positions
    span = edge_length / (count + 1) * 0.34
    window_half = facade_panel_halfwidth(2 * span, physical_window=level == 'lod0')
    door_half = facade_panel_halfwidth(min(edge_length, 3.5), width_ratio=0.72)
    clearance = door_half + window_half + 0.15
    adjusted = []
    for t in positions:
        offset = (t - 0.5) * edge_length
        if abs(offset) >= clearance:
            adjusted.append(t)
            continue
        if abs(offset) < 1e-9:
            continue  # a centered window occupies the entrance; keep the upper one
        center = edge_length * 0.5 + math.copysign(clearance, offset)
        # Do not solve a door conflict by making a wall-end or neighboring-window conflict.
        if center - window_half < 0.15 or center + window_half > edge_length - 0.15:
            continue
        if any(abs(center - other * edge_length) < 2 * window_half + 0.15
               for other in positions if other != t):
            continue
        adjusted.append(center / edge_length)
    return adjusted


def massing_triangles(footprint, floors, level, facade_scope='building', facade_mask=None):
    """Reference triangles for a massing prism plus lightweight facade details."""
    if facade_scope not in ('building', 'shell-only'):
        raise ValueError('unknown massing facade scope: ' + str(facade_scope))
    base = prism_triangles(len(footprint))
    if facade_scope == 'shell-only':
        return base
    mask = facade_mask or {}
    hidden_windows = mask.get('windows', set())
    hidden_bands = mask.get('bands', set())
    longest = facade_entrance_edge(footprint)
    windows = 0
    for i in range(len(footprint)):
        a, b = footprint[i], footprint[(i + 1) % len(footprint)]
        length = math.hypot(b[0] - a[0], b[1] - a[1])
        count = facade_window_count(length, level)
        # add_panel/add_window_reveal reject a sub-0.35m window span. Surveyed
        # cornices/notches have short edges; they remain in the solid, without windows.
        if length / (count + 1) * 0.68 >= 0.35:
            windows += 2 * sum(len(facade_window_positions(length, level,
                reserve_entrance=i == longest and floor == 0 and not mask.get('door')))
                for floor in range(floors) if (i, floor) not in hidden_windows)
    roof_edges = sum(math.hypot(footprint[(i + 1) % len(footprint)][0] - a[0], footprint[(i + 1) % len(footprint)][1] - a[1]) >= .35 for i, a in enumerate(footprint)) if FULL_REVEALS else 0
    storey_bands = sum((i, floor) not in hidden_bands for i in range(len(footprint))
                       for floor in range(1, floors))
    bands = storey_bands if BANDS and ARCHITECTURE and level == 'lod0' else 0
    architectural = ((windows * 4 if FULL_REVEALS else 0)
                     + 2 * storey_bands
                     + 2 * roof_edges + bands * 2) if ARCHITECTURE and level == 'lod0' else 0
    return base + windows + architectural + (0 if mask.get('door') else 12)


def poly_prism(footprint, bottom, top, wall_material, roof_material, anchor, floors, level,
               window_material, door_material, facade_scope='building', facade_mask=None):
    """Build a surveyed-outline building with bounded facade details.

    The solid uses the exact GIS ring. Windows are shallow, outward-facing quads repeated per
    floor and door is one first-floor quad on the longest wall. They add a readable facade without
    adding separate Blender objects or changing the building footprint.
    """
    if facade_scope not in ('building', 'shell-only'):
        raise ValueError('unknown massing facade scope: ' + str(facade_scope))
    mask = facade_mask or {}
    hidden_windows = mask.get('windows', set())
    hidden_bands = mask.get('bands', set())
    # Create the mesh directly: repeated scene-wide select/cube operators make city exports quadratic.
    # This is an exterior-only massing shell: no rooms, floors, ceilings, or interior props.
    mesh = bpy.data.meshes.new('Surveyed_building')
    obj = bpy.data.objects.new('Surveyed_building', mesh)
    n = len(footprint)
    ax, ay, az = anchor
    vertices = [(x - ax, y - ay, bottom - az) for x, y in footprint] + \
               [(x - ax, y - ay, top - az) for x, y in footprint]
    faces = [(i, (i + 1) % n, n + (i + 1) % n, n + i) for i in range(n)]
    materials = [0] * n
    # Explicit caps survive merging/export and also handle concave surveyed outlines.
    # Use the same footprint vertices; no hull, overhang or XY displacement is introduced.
    cap_points = [Vector((x - ax, y - ay, 0)) for x, y in footprint]
    cap_indices = {tuple(point): i for i, point in enumerate(cap_points)}
    roof_area = 0.0
    for triangle in tessellate_polygon([cap_points]):
        # Blender 5.2 returns point indices; earlier versions returned Vector triples.
        indices = [point if isinstance(point, int) else cap_indices[tuple(point)] for point in triangle]
        area = ring_area([footprint[index] for index in indices])
        if abs(area) < 1e-10:
            continue
        if area < 0:
            indices.reverse()
        faces.append(tuple(n + index for index in indices))  # roof points up
        materials.append(1)
        faces.append(tuple(reversed(indices)))  # buried bottom points down
        materials.append(0)
        roof_area += abs(area)
    # GIS vertices are rounded to centimetres. One real ring contains a 1cm backtrack and its
    # signed area differs by 0.08m² from the filled roof. Bound that tolerance by source precision,
    # instead of accepting missing caps or moving the whole surveyed footprint.
    perimeter = sum(math.hypot(footprint[(i + 1) % n][0] - p[0],
                               footprint[(i + 1) % n][1] - p[1]) for i, p in enumerate(footprint))
    tolerance = max(0.001, roof_area * 1e-6, perimeter * 0.01)
    if abs(roof_area - ring_area(footprint)) > tolerance:
        raise RuntimeError('Triangulated roof does not cover its GIS footprint')
    floor_base = max(bottom, az) if ARCHITECTURE else bottom
    floor_height = max(2.4, (top - floor_base) / max(1, floors))

    def edge_panel_points(a, b, z0, z1, inset, width_ratio=0.42, physical_window=False):
        dx, dy = b[0] - a[0], b[1] - a[1]
        length = math.hypot(dx, dy)
        if length < 0.35 or z1 <= z0:
            return None
        nx, ny = dy / length, -dx / length
        half = facade_panel_halfwidth(length, width_ratio, physical_window)
        cx, cy = (a[0] + b[0]) * 0.5, (a[1] + b[1]) * 0.5
        tx, ty = dx / length, dy / length
        return [(cx - tx * half + nx * inset, cy - ty * half + ny * inset),
                (cx + tx * half + nx * inset, cy + ty * half + ny * inset)]

    def add_panel(a, b, z0, z1, material_index, inset=0.035):
        points = edge_panel_points(a, b, z0, z1, inset,
                                   physical_window=FACADE_PROFILE == 'compact' and level == 'lod0')
        if not points:
            return
        base = len(vertices)
        vertices.extend([(x - ax, y - ay, z - az) for x, y in
                         [(points[0][0], points[0][1]), (points[1][0], points[1][1])]
                         for z in (z0, z1)])
        faces.append((base, base + 2, base + 3, base + 1))
        materials.append(material_index)

    def add_window_reveal(a, b, z0, z1):
        # Four solid-depth reveal strips around glass; no unsupported balcony/roof invention.
        inner = edge_panel_points(a, b, z0, z1, 0.035)
        outer = edge_panel_points(a, b, z0, z1, 0.16, 0.52)
        if not inner or not outer:
            return
        base = len(vertices)
        vertices.extend([(x - ax, y - ay, z - az) for points, low, high in
                         ((outer, z0 - 0.10, z1 + 0.10), (inner, z0, z1))
                         for x, y in points for z in (low, high)])
        faces.extend([(base, base + 2, base + 6, base + 4),
                      (base + 1, base + 5, base + 7, base + 3),
                      (base, base + 4, base + 5, base + 1),
                      (base + 2, base + 3, base + 7, base + 6)])
        materials.extend([0] * 4)

    def add_box_panel(a, b, z0, z1, material_index):
        """Closed shallow door box; keeps a door material through merged GLB export."""
        if ARCHITECTURE:
            # Entrance width is architectural scale, not a fraction of a 70m apartment facade.
            dx, dy = b[0] - a[0], b[1] - a[1]
            length = math.hypot(dx, dy) or 1
            cx, cy = (a[0] + b[0]) / 2, (a[1] + b[1]) / 2
            half = min(length, 3.5) / 2
            a = (cx - dx / length * half, cy - dy / length * half)
            b = (cx + dx / length * half, cy + dy / length * half)
        points = edge_panel_points(a, b, z0, z1, 0.055, 0.72)
        if not points:
            return
        dx, dy = b[0] - a[0], b[1] - a[1]
        length = math.hypot(dx, dy) or 1
        nx, ny = dy / length, -dx / length
        outer = 0.14
        inner = 0.045
        p = points
        base = len(vertices)
        for depth in (inner, outer):
            vertices.extend([(p[0][0] - ax + nx * depth, p[0][1] - ay + ny * depth, z - az)
                             for z in (z0, z1)])
            vertices.extend([(p[1][0] - ax + nx * depth, p[1][1] - ay + ny * depth, z - az)
                             for z in (z0, z1)])
        faces.extend([(base, base + 1, base + 3, base + 2),
                      (base + 4, base + 6, base + 7, base + 5),
                      (base, base + 4, base + 5, base + 1),
                      (base + 2, base + 3, base + 7, base + 6),
                      (base, base + 2, base + 6, base + 4),
                      (base + 1, base + 5, base + 7, base + 3)])
        materials.extend([material_index] * 6)

    longest = facade_entrance_edge(footprint)
    # Reviewed attachments retain closed GIS geometry without independent house details.
    for i in range(n if facade_scope == 'building' else 0):
        a, b = footprint[i], footprint[(i + 1) % n]
        count = facade_window_count(math.hypot(b[0] - a[0], b[1] - a[1]), level)
        dx, dy = b[0] - a[0], b[1] - a[1]
        for floor in range(max(1, floors)):
            if (i, floor) in hidden_windows:
                continue
            compact_window = FACADE_PROFILE == 'compact' and level == 'lod0'
            z0 = floor_base + floor * floor_height + floor_height * (0.30 if compact_window else 0.34)
            z1 = min(top - 0.08, z0 + floor_height * (0.45 if compact_window else 0.36))
            for t in facade_window_positions(math.hypot(dx, dy), level,
                    reserve_entrance=i == longest and floor == 0 and not mask.get('door')):
                center = (a[0] + dx * t, a[1] + dy * t)
                span = math.hypot(dx, dy) / (count + 1) * 0.34
                wa = (center[0] - dx / (math.hypot(dx, dy) or 1) * span,
                      center[1] - dy / (math.hypot(dx, dy) or 1) * span)
                wb = (center[0] + dx / (math.hypot(dx, dy) or 1) * span,
                      center[1] + dy / (math.hypot(dx, dy) or 1) * span)
                add_panel(wa, wb, z0, z1, 2)
                if FULL_REVEALS and level == 'lod0':
                    add_window_reveal(wa, wb, z0, z1)
        if ARCHITECTURE and level == 'lod0':
            for floor in range(1, max(1, floors)):
                if (i, floor) in hidden_bands:
                    continue
                z = floor_base + floor * floor_height
                # A thin full-edge storey band uses the darker door tone, not extra objects.
                length = math.hypot(dx, dy) or 1
                nx, ny = dy / length, -dx / length
                base = len(vertices)
                vertices.extend([(p[0] - ax + nx * 0.025, p[1] - ay + ny * 0.025, zz - az)
                                 for p in (a, b) for zz in (z - 0.045, z + 0.045)])
                faces.append((base, base + 2, base + 3, base + 1))
                materials.append(3)

    if facade_scope == 'building' and BANDS and ARCHITECTURE and level == 'lod0':
        for i in range(n):
            a, b = footprint[i], footprint[(i + 1) % n]
            dx, dy = b[0] - a[0], b[1] - a[1]
            length = math.hypot(dx, dy)
            if length < .35: continue
            nx, ny = dy / length, -dx / length
            for floor in range(1, max(1, floors)):
                if (i, floor) in hidden_bands:
                    continue
                z = floor_base + floor * floor_height
                inset = .028
                base = len(vertices)
                vertices.extend([(a[0] - ax + nx * inset, a[1] - ay + ny * inset, z - .035 - az),
                                 (b[0] - ax + nx * inset, b[1] - ay + ny * inset, z - .035 - az),
                                 (b[0] - ax + nx * inset, b[1] - ay + ny * inset, z + .035 - az),
                                 (a[0] - ax + nx * inset, a[1] - ay + ny * inset, z + .035 - az)])
                faces.append((base, base + 2, base + 3, base + 1)); materials.append(3)


    if facade_scope == 'building' and FULL_REVEALS and level == 'lod0':
        for i in range(n):
            a, b = footprint[i], footprint[(i + 1) % n]
            length = math.hypot(b[0] - a[0], b[1] - a[1])
            if length < 0.35: continue
            dx, dy = b[0] - a[0], b[1] - a[1]
            norm = math.hypot(dx, dy) or 1
            nx, ny = dy / norm, -dx / norm
            inset = 0.045
            low, high = top - 0.18, top + 0.08
            base = len(vertices)
            vertices.extend([(a[0] - ax + nx * inset, a[1] - ay + ny * inset, low - az),
                             (b[0] - ax + nx * inset, b[1] - ay + ny * inset, low - az),
                             (b[0] - ax + nx * inset, b[1] - ay + ny * inset, high - az),
                             (a[0] - ax + nx * inset, a[1] - ay + ny * inset, high - az)])
            faces.append((base, base + 1, base + 2, base + 3)); materials.append(0)

    # One door on the original longest wall, visible from at least one approach direction.
    edge_a = footprint[longest]
    edge_b = footprint[(longest + 1) % n]
    door_h = min(2.2, max(1.8, floor_height * 0.76))
    if facade_scope == 'building' and not mask.get('door'):
        add_box_panel(edge_a, edge_b, floor_base + 0.08, min(top - 0.05, floor_base + door_h), 3)
    mesh.from_pydata(vertices, [], faces)
    obj.location = anchor
    for material in (wall_material, roof_material, window_material, door_material):
        mesh.materials.append(material)
    mesh.polygons.foreach_set('material_index', materials)
    # This slot belongs to our authored window planes, never to library assets.
    obj['cad_authored_window_slot'] = 2
    edit = bmesh.new()
    edit.from_mesh(mesh)
    # An isolated plane has no volume from which bmesh can infer outside.
    # Keep authored windows in every profile, plus the existing compact planes.
    normal_faces = [face for face in edit.faces if not (
        (face.material_index == 2 or FACADE_PROFILE == 'compact')
        and all(edge.is_boundary for edge in face.edges))]
    bmesh.ops.recalc_face_normals(edit, faces=normal_faces)
    edit.to_mesh(mesh)
    edit.free()
    mesh.update()
    return obj


def prism_volume_sign(obj):
    """Volume of obj's closed surface read from Blender's own face normals.

    V = sum(centre . normal * area) / 3 is positive exactly when the normals point outward, which is
    the property the GLB export and every downstream renderer depend on. Local space is enough: the
    prisms are built unrotated and unscaled, and reading Blender's normals avoids re-deriving the
    winding sign by hand.
    """
    return sum(polygon.center.dot(polygon.normal) * polygon.area
               for polygon in obj.data.polygons) / 3


def box(bounds):
    low, high = bounds
    centre = [(a + b) / 2 for a, b in zip(low, high)]
    half = [(b - a) / 2 for a, b in zip(low, high)]
    return [centre[0], centre[1], centre[2], half[0], 0, 0, 0, half[1], 0, 0, 0, half[2]]


# ---------------------------------------------------------------- plan intake
def build_geometry(building):
    """Return (level -> (outline, triangles)) for one building, validating its outline."""
    outline = [(float(x), float(y)) for x, y in building['footprintENU']]
    if len(outline) < 3:
        raise RuntimeError(building['id'] + ' has fewer than 3 outline points')
    outline = as_ccw(outline)
    if building['kind'] == 'asset' and not args.lod1_footprint_only:
        # The library mesh is already authored in this ring; only the transform differs per level.
        return {level: (None, None) for level in ('lod0', 'lod1')}
    if ring_area(outline) <= 0:
        raise RuntimeError(building['id'] + ' outline is degenerate')
    if building['positionENU'][2] + building['heightM'] <= building['foundationBottom']:
        raise RuntimeError(building['id'] + ' massing prism has zero volume')
    # Both runtime levels use the surveyed ring. LOD1 reduces facade detail, not the GIS footprint;
    # replacing an irregular ring with its hull was the source of visibly misplaced road-covering mass.
    scope = building.get('massingFacadeScope', 'building')
    masks = PARTY_WALL_MASKS.get(building['id'], {})
    return {'lod0': (outline, massing_triangles(outline, building['floors'], 'lod0', scope, masks.get('lod0'))),
            'lod1': (outline, massing_triangles(outline, building['floors'], 'lod1', scope, masks.get('lod1')))}


placements, rejected, expected_triangles = [], [], {}
selected_input = plan['buildings'][:args.limit or None]
PARTY_WALL_MASKS = {}
if any(building.get('partyWallNeighbors') for building in selected_input):
    sys.path.insert(0, str(Path(__file__).resolve().parent))
    from party_wall_facades import party_wall_masks
    if FULL_REVEALS or BANDS:
        raise RuntimeError('Reviewed party-wall masks require compact without reveals or extra storey bands')
    PARTY_WALL_MASKS = party_wall_masks(selected_input, as_ccw, FACADE_PROFILE)
if any('massingFacadeScope' in building for building in selected_input):
    # Blender's script runner does not always add the script directory to sys.path.
    sys.path.insert(0, str(Path(__file__).resolve().parent))
    from building_composition import validate_composition_scopes
    validate_composition_scopes(selected_input, FACADE_PROFILE)
library_ids = {entry['id'] for entry in entries}
for building in selected_input:
    if building['kind'] == 'asset':
        # A named variant that is not in the library is a broken plan, not a building to skip.
        unknown = set(building.get('preferredAssets') or []).difference(library_ids)
        if unknown:
            raise RuntimeError(f"{building['id']} names missing library objects: " + ', '.join(sorted(unknown)))
        choice = select_asset(building)
        if not choice:
            rejected.append(building['id'])
            continue
        _, entry, sx, sy = choice
        expected_triangles[building['id']] = None  # filled from the library once templates load
        placements.append((building, entry, sx, sy))
    elif building['kind'] == 'massing':
        geometry = facade_operation(building, build_geometry, building)
        expected_triangles[building['id']] = {level: geometry[level][1] for level in geometry}
        placements.append((building, None, None, None))
    else:
        raise RuntimeError(building['id'] + ' has unknown kind ' + str(building['kind']))
if not placements or len({b['id'] for b, _, _, _ in placements}) != len(placements):
    raise RuntimeError('The seated plan does not yield unique, nonempty placements')
placed_ids = {building['id'] for building, _, _, _ in placements}
for building, _, _, _ in placements:
    if building.get('massingFacadeScope') == 'shell-only' and building['composition']['primaryId'] not in placed_ids:
        raise RuntimeError('Attachment primary did not pass asset placement: ' + building['id'])

# ---------------------------------------------------------------- scene setup
bpy.ops.wm.read_factory_settings(use_empty=True)
bpy.context.scene.unit_settings.system = 'METRIC'
bpy.context.scene.unit_settings.scale_length = 1
OUT.mkdir(parents=True, exist_ok=True)
for pattern in ('*-lod0.glb', '*-lod1.glb'):
    for stale in OUT.glob(pattern):
        stale.unlink()
for stale in (OUT / 'root-anchor.glb', OUT / 'tileset.json', OUT / 'manifest.json'):
    if stale.exists():
        stale.unlink()
stone = bpy.data.materials.new('Buried_foundation_stone')
stone.diffuse_color = (0.32, 0.33, 0.32, 1)
stone.use_nodes = True
stone.node_tree.nodes['Principled BSDF'].inputs['Base Color'].default_value = stone.diffuse_color
stone.node_tree.nodes['Principled BSDF'].inputs['Roughness'].default_value = 0.92

# The dimensions/height vary per foundation; scale a shared unit cube mesh per instance.
foundation_mesh = bpy.data.meshes.new('Seated_foundation_unit_cube')
foundation_mesh.from_pydata(
    [(-.5, -.5, -.5), (.5, -.5, -.5), (.5, .5, -.5), (-.5, .5, -.5),
     (-.5, -.5, .5), (.5, -.5, .5), (.5, .5, .5), (-.5, .5, .5)],
    [],
    [(3, 2, 1, 0), (4, 5, 6, 7), (0, 1, 5, 4),
     (1, 2, 6, 5), (2, 3, 7, 6), (3, 0, 4, 7)],
)
foundation_mesh.materials.append(stone)

# Load the placed library templates together so common materials retain their names. Massing
# buildings reuse those same material datablocks, so no extra palette has to be authored here.
required = {entry['joined_object'] for _, entry, _, _ in placements if entry}
required.update(entry['lod1']['joined_object'] for _, entry, _, _ in placements if entry)
with bpy.data.libraries.load(str(LIBRARY_DIR / 'house-library.blend'), link=False) as (src, dst):
    missing = required.difference(src.objects)
    if missing:
        raise RuntimeError('Missing Blender library templates: ' + ', '.join(sorted(missing)))
    dst.objects = sorted(required)
templates = {obj.name: obj for obj in dst.objects if obj is not None}
if len(templates) != len(required) or any(obj.type != 'MESH' or not obj.data.polygons
                                          for obj in templates.values()):
    raise RuntimeError('A required Blender template is empty')

# Asset-level reduction: a Decimate modifier on the TEMPLATE mesh collapses every placed instance
# at once (object copies share mesh data). Ratios come from the CLI; 1.0 (default) leaves meshes
# untouched. Collapsing an inset window frame to one sheet is acceptable at these viewing
# distances; what must survive is the building silhouette, its openings as dark insets, and the
# storey rhythm — the maintainer's acceptance bar for "looks like a building".
DECIMATE_RATIOS = {'lod0': args.decimate_lod0, 'lod1': args.decimate_lod1}
decimated_templates = {}
if any(ratio < 1.0 for ratio in DECIMATE_RATIOS.values()):
    for name, obj in templates.items():
        level = 'lod0' if name.endswith('_LOD0') else 'lod1' if name.endswith('_LOD1') else None
        ratio = DECIMATE_RATIOS.get(level, 1.0)
        if ratio < 1.0:
            modifier = obj.modifiers.new('AssetDecimate', 'DECIMATE')
            modifier.ratio = ratio
            modifier.use_collapse_triangulate = True
            decimated_templates[name] = modifier
    has_asset_placements = any(entry for _, entry, _, _ in placements)
    if not decimated_templates and has_asset_placements:
        raise RuntimeError('Decimate requested but no template name ends with _LOD0/_LOD1')
    # A massing-only batch (no asset placements, e.g. the unfitted-rectangle tail) has no templates
    # to decimate — its prisms are already minimal, so Draco-only export is the correct behavior.
if decimated_templates:
    # Library.json triangle counts no longer describe the decimated meshes: link the decimated
    # templates into the view layer, evaluate them, and measure the real triangulated counts.
    scratch = bpy.data.collections.new('Template_eval_scratch')
    bpy.context.scene.collection.children.link(scratch)
    for name in decimated_templates:
        scratch.objects.link(templates[name])
    bpy.context.view_layer.update()
    depsgraph = bpy.context.evaluated_depsgraph_get()
    evaluated_triangles = {}
    for name in decimated_templates:
        evaluated = templates[name].evaluated_get(depsgraph)
        mesh = evaluated.to_mesh()
        evaluated_triangles[name] = sum(len(polygon.vertices) - 2 for polygon in mesh.polygons)
        evaluated.to_mesh_clear()
    for building, entry, _, _ in placements:
        if entry:
            expected_triangles[building['id']] = {}
            for level in ('lod0', 'lod1'):
                if args.lod1_footprint_only and level == 'lod1':
                    expected_triangles[building['id']][level] = facade_operation(building, build_geometry, building)[level][1]
                    continue
                template_name = entry['joined_object'] if level == 'lod0' else entry['lod1']['joined_object']
                fallback = (entry['validation'] if level == 'lod0' else entry['lod1']['validation'])['triangles']
                expected_triangles[building['id']][level] = evaluated_triangles.get(template_name, fallback)
else:
    for building, entry, _, _ in placements:
        if entry:
            expected_triangles[building['id']] = {
                level: (facade_operation(building, build_geometry, building)[level][1] if args.lod1_footprint_only and level == 'lod1'
                        else (entry['validation'] if level == 'lod0' else entry['lod1']['validation'])['triangles'])
                for level in ('lod0', 'lod1')}

# Append the whole library palette on top of the templates: templates that happen to share no
# material would otherwise leave a massing tone name unavailable.
with bpy.data.libraries.load(str(LIBRARY_DIR / 'house-library.blend'), link=False) as (src, dst):
    dst.materials = list(src.materials)
missing_tones = set(WALL_TONES + ROOF_TONES + ('Glass_Dark', 'Door_Walnut')).difference(bpy.data.materials.keys())
if missing_tones:
    raise RuntimeError('Library is missing material tones: ' + ', '.join(sorted(missing_tones)))


def id_seed(building_id):
    value = 2166136261
    for character in str(building_id):
        value = ((value ^ ord(character)) * 16777619) % 4294967296
    return value


def massing_tone(building_id):
    seed = id_seed(building_id)
    return (bpy.data.materials[WALL_TONES[seed % len(WALL_TONES)]],
            bpy.data.materials[ROOF_TONES[(seed // len(WALL_TONES)) % len(ROOF_TONES)]],
            bpy.data.materials['Glass_Dark'], bpy.data.materials['Door_Walnut'])


def sourced_roof(building_id, fallback):
    record = ROOF_TEXTURES.get(building_id)
    if not record:
        return fallback
    name = 'SatelliteRoof_' + building_id
    material = bpy.data.materials.get(name)
    if material:
        return material
    material = bpy.data.materials.new(name)
    material.use_nodes = True
    shader = material.node_tree.nodes.get('Principled BSDF')
    shader.inputs['Roughness'].default_value = 0.85
    shader.inputs['Metallic'].default_value = 0.0
    texture = material.node_tree.nodes.new('ShaderNodeTexImage')
    texture.image = bpy.data.images.load(record['path'], check_existing=True)
    texture.extension = 'CLIP'
    material.node_tree.links.new(texture.outputs['Color'], shader.inputs['Base Color'])
    return material


def compact_window_material():
    name = 'Compact_Window_Frame_' + WINDOW_TEXTURE['sha256'][:12]
    material = bpy.data.materials.get(name)
    if material:
        return material
    material = bpy.data.materials.new(name)
    material.use_nodes = True
    shader = material.node_tree.nodes.get('Principled BSDF')
    shader.inputs['Base Color'].default_value = (1, 1, 1, 1)
    shader.inputs['Roughness'].default_value = 0.65
    shader.inputs['Metallic'].default_value = 0.0
    texture = material.node_tree.nodes.new('ShaderNodeTexImage')
    texture.image = bpy.data.images.load(WINDOW_TEXTURE['path'], check_existing=True)
    texture.extension = 'CLIP'
    material.node_tree.links.new(texture.outputs['Color'], shader.inputs['Base Color'])
    return material


def apply_window_uv(obj):
    faces = [face for face in obj.data.polygons if face.material_index == 2]
    if not faces:
        return
    uv = obj.data.uv_layers.active or obj.data.uv_layers.new(name='UVMap')
    for face in faces:
        points = [tuple(obj.data.vertices[obj.data.loops[index].vertex_index].co)
                  for index in face.loop_indices]
        for index, value in zip(face.loop_indices, window_quad_uv(points)):
            uv.data[index].uv = value


def apply_roof_uv(obj, building_id):
    record = ROOF_TEXTURES.get(building_id)
    if not record:
        return
    low, high = record['boundsENU']
    uv = obj.data.uv_layers.active or obj.data.uv_layers.new(name='UVMap')
    roof_loops = (obj.data.loops[index] for face in obj.data.polygons
                  if face.material_index == 1 for index in face.loop_indices)
    for loop in roof_loops:
        point = obj.data.vertices[loop.vertex_index].co + obj.location
        uv.data[loop.index].uv = ((point.x - low[0]) / (high[0] - low[0]),
                                  (point.y - low[1]) / (high[1] - low[1]))


def material_remap(building_id):
    """Asset palette rotation, keyed by the same stable hash as the massing tones."""
    seed = id_seed(building_id)
    roof_shift = (seed // 3) % len(ROOF_TONES)
    wall_shift = seed % len(WALL_TONES)
    remap = {name: ROOF_TONES[(index + roof_shift) % len(ROOF_TONES)]
             for index, name in enumerate(ROOF_TONES)}
    remap.update({name: WALL_TONES[(index + wall_shift) % len(WALL_TONES)]
                  for index, name in enumerate(WALL_TONES)})
    remap['Roof_Membrane'] = ROOF_TONES[roof_shift]
    return remap


def template(entry, level):
    return templates[entry['joined_object'] if level == 'lod0' else entry['lod1']['joined_object']]


def make_cell(group, level):
    collection = bpy.data.collections.new('Cad_' + level)
    bpy.context.scene.collection.children.link(collection)
    for building, entry, sx, sy in group:
        x, y, top = building['positionENU']
        bottom = building['foundationBottom']
        if entry and not (args.lod1_footprint_only and level == 'lod1'):
            obj = template(entry, level).copy()  # Object copy shares mesh data with siblings.
            obj.name = building['id']
            collection.objects.link(obj)
            obj.location = building['positionENU']
            obj.rotation_euler = (0, 0, building['heading'])
            obj.scale = (sx, sy, 1)
            obj['source_footprint'] = building.get('sourceId', building['id'])
            remap = material_remap(building['id'])
            for slot in obj.material_slots:
                if slot.material and slot.material.name in remap:
                    material = bpy.data.materials.get(remap[slot.material.name])
                    if material:
                        slot.link = 'OBJECT'
                        slot.material = material
        else:
            geometry = facade_operation(building, build_geometry, building)[level]
            outline = geometry[0]
            wall, roof, window, door = massing_tone(building['id'])
            textured_windows = (WINDOW_TEXTURE is not None and level == 'lod0'
                                and (COMPACT_BUILDING_IDS is None or building['id'] in COMPACT_BUILDING_IDS)
                                and building.get('massingFacadeScope', 'building') == 'building')
            if textured_windows:
                window = compact_window_material()
            if level == 'lod0':
                roof = sourced_roof(building['id'], roof)
            top_z = top + building['heightM']
            obj = facade_operation(building, poly_prism, outline, bottom, top_z, wall, roof, building['positionENU'],
                            max(1, int(building.get('floors') or round(building['heightM'] / 3.2))), level,
                            window, door, facade_scope=building.get('massingFacadeScope', 'building'),
                                 facade_mask=PARTY_WALL_MASKS.get(building['id'], {}).get(level))
            # Self-check the prism is a solid with outward normals. If the ring somehow came out
            # clockwise, rebuild once reversed instead of silently emitting an inside-out building.
            if prism_volume_sign(obj) <= 0:
                bpy.data.meshes.remove(obj.data, do_unlink=True)
                if building['id'] in PARTY_WALL_MASKS:
                    raise RuntimeError('Reviewed party-wall prism needs normal correction; do not reuse reversed edge masks: ' + building['id'])
                obj = facade_operation(building, poly_prism, list(reversed(outline)), bottom, top_z, wall, roof, building['positionENU'],
                                 max(1, int(building.get('floors') or round(building['heightM'] / 3.2))), level,
                                 window, door, facade_scope=building.get('massingFacadeScope', 'building'),
                                 facade_mask=PARTY_WALL_MASKS.get(building['id'], {}).get(level))
                if prism_volume_sign(obj) <= 0:
                    raise RuntimeError(building['id'] + ' massing prism is inside-out or flat')
            if textured_windows:
                apply_window_uv(obj)
            if level == 'lod0':
                apply_roof_uv(obj, building['id'])
            obj.name = building['id']
            for collection_link in list(obj.users_collection):
                collection_link.objects.unlink(obj)
            collection.objects.link(obj)
        # Buried foundation. An asset instance is a library mesh seated at floor level, so it needs
        # its own box down to foundationBottom. A massing prism already spans foundationBottom to the
        # roof, i.e. its buried foundation IS the prism; a rectangular box under an L-shaped or
        # rotated outline would poke outside the building and z-fight its walls. So one cell holds
        # two mesh nodes per asset placement and one per massing building.
        if entry and not (args.lod1_footprint_only and level == 'lod1'):
            foundation = bpy.data.objects.new(building['id'] + '_buried_foundation', foundation_mesh)
            collection.objects.link(foundation)
            foundation.location = (x, y, (bottom + top) / 2)
            foundation.rotation_euler.z = building['heading']
            foundation.scale = (building['width'], building['depth'], top - bottom)
    bpy.context.view_layer.update()
    objects = list(collection.objects)
    if args.merge_buildings and objects:
        objects = merge_cell_objects(objects, collection)
    corners = [obj.matrix_world @ Vector(corner) for obj in objects for corner in obj.bound_box]
    bounds = [[min(p[axis] for p in corners) for axis in range(3)],
              [max(p[axis] for p in corners) for axis in range(3)]]
    return collection, objects, bounds


def merge_cell_objects(objects, collection):
    """Bake every building/foundation in the cell into ONE mesh.

    The per-building export produces a separate node and material primitives per placement.
    Merging reduces that structural cost to one primitive per used material; actual frame-time
    improvement still requires a runtime measurement. The decimate modifier is baked via the
    depsgraph. Transforms are baked into ENU vertices; normals are recalculated with flat shading.
    Material slots are preserved. Roof/window appearance must still be visually reviewed.
    """
    bpy.context.view_layer.update()
    depsgraph = bpy.context.evaluated_depsgraph_get()
    vertices, faces, face_material_names, face_uvs = [], [], [], []
    material_names = {}
    preserve_planar_faces = []
    scope_ids = globals().get('COMPACT_BUILDING_IDS')
    has_uv = False
    for obj in objects:
        evaluated = obj.evaluated_get(depsgraph)
        mesh = evaluated.to_mesh()
        has_uv = has_uv or mesh.uv_layers.active is not None
        matrix = obj.matrix_world.copy()
        base = len(vertices)
        authored_window_slot = obj.get('cad_authored_window_slot')
        for vertex in mesh.vertices:
            vertices.append(tuple(matrix @ vertex.co))
        for polygon in mesh.polygons:
            corner_indices = tuple(base + index for index in polygon.vertices)
            # Reject invalid faces with fewer than three distinct indices before rebuilding.
            if len(set(corner_indices)) < 3:
                continue
            faces.append(corner_indices)
            compact_plane = FACADE_PROFILE == 'compact' and (scope_ids is None or getattr(obj, 'name', '') in scope_ids)
            preserve_planar_faces.append(compact_plane or polygon.material_index == authored_window_slot)
            uv_layer = mesh.uv_layers.active
            face_uvs.append([tuple(uv_layer.data[i].uv) if uv_layer else (0.0, 0.0)
                             for i in polygon.loop_indices])
            material = obj.material_slots[polygon.material_index].material if obj.material_slots else None
            key = material.name if material else ''
            if key not in material_names:
                material_names[key] = len(material_names)
            face_material_names.append(key)
        evaluated.to_mesh_clear()
    merged_mesh = bpy.data.meshes.new('cell_merged')
    merged_mesh.from_pydata(vertices, [], faces)
    if ROOF_TEXTURES or (WINDOW_TEXTURE is not None and has_uv):
        uv_layer = merged_mesh.uv_layers.new(name='UVMap')
        for polygon, values in zip(merged_mesh.polygons, face_uvs):
            for loop_index, value in zip(polygon.loop_indices, values):
                uv_layer.data[loop_index].uv = value
    for name in material_names:
        merged_mesh.materials.append(bpy.data.materials[name])
    merged_mesh.polygons.foreach_set('material_index', [material_names[n] for n in face_material_names])
    edit = bmesh.new()
    edit.from_mesh(merged_mesh)
    # Preserve authored windows in every profile and the existing selected compact
    # planes. Assets/foundations without that provenance keep their normal path.
    normal_faces = list(edit.faces)
    if len(normal_faces) != len(preserve_planar_faces):
        raise RuntimeError('Merged face provenance changed')
    normal_faces = [face for face, preserve in zip(normal_faces, preserve_planar_faces)
                    if not (preserve and all(edge.is_boundary for edge in face.edges))]
    bmesh.ops.recalc_face_normals(edit, faces=normal_faces)
    edit.to_mesh(merged_mesh)
    edit.free()
    merged_mesh.update()
    merged = bpy.data.objects.new('cell_merged', merged_mesh)
    collection.objects.link(merged)
    for obj in objects:
        bpy.data.objects.remove(obj, do_unlink=True)
    return [merged]


def export_cell(path, objects, group, level):
    bpy.ops.object.select_all(action='DESELECT')
    for obj in objects:
        obj.select_set(True)
    bpy.context.view_layer.objects.active = objects[0]
    export_options = dict(filepath=str(path), export_format='GLB', use_selection=True,
                          export_yup=True, export_materials='EXPORT', export_extras=True,
                          # Bake modifiers (the asset Decimate) into the exported mesh; without
                          # this the GLB would carry the full-detail template meshes again.
                          export_apply=True)
    if args.draco:
        export_options['export_draco_mesh_compression_enable'] = True
        if ARCHITECTURE:
            # Surveyed roofs include near-collinear centimetre edges: default quantization can invert slivers.
            export_options['export_draco_position_quantization'] = 24
    bpy.ops.export_scene.gltf(**export_options)
    with path.open('rb') as stream:
        header = stream.read(20)
        if len(header) != 20 or header[:4] != b'glTF':
            raise RuntimeError('Invalid GLB export ' + str(path))
        _, version, length, chunk_length, chunk_type = struct.unpack('<4s4I', header)
        if version != 2 or length != path.stat().st_size or chunk_type != 0x4E4F534A:
            raise RuntimeError('Invalid GLB header ' + str(path))
        doc = json.loads(stream.read(chunk_length))
    # Asset instances are one node plus their buried foundation; a massing prism is a single node.
    if args.merge_buildings:
        # Merged cell: exactly one mesh node carries the whole cell. Per-building node identity
        # and seated transforms are baked into world-space vertices, so only the triangle total
        # (below) and the scene-root shape remain meaningful here.
        names = [n.get('name') for n in doc['nodes'] if 'mesh' in n]
        if len(names) != 1:
            raise RuntimeError(f'{path.name}: merged cell must contain exactly one mesh node')
    else:
        expected_names = []
        for building, entry, _, _ in group:
            expected_names.append(building['id'])
            if entry and not (args.lod1_footprint_only and level == 'lod1'):
                expected_names.append(building['id'] + '_buried_foundation')
        names = [n.get('name') for n in doc['nodes'] if 'mesh' in n]
        # One node per massing prism, two per asset placement (instance + buried foundation).
        if len(names) != len(expected_names):
            raise RuntimeError(f'{path.name}: {len(names)} mesh nodes for {len(expected_names)} '
                               'expected building/foundation nodes')
        if sorted(names) != sorted(expected_names):
            raise RuntimeError('GLB mesh nodes do not match placed buildings and foundations: ' + str(path))
        nodes = {node['name']: node for node in doc['nodes']}
        for building, entry, sx, sy in group:
            x, y, top = building['positionENU']
            bottom = building['foundationBottom']
            use_asset = entry and not (args.lod1_footprint_only and level == 'lod1')
            checks = [(building['id'], top, (sx, 1, sy) if use_asset else None)]
            if use_asset:
                checks.append((building['id'] + '_buried_foundation', (bottom + top) / 2,
                               (building['width'], top - bottom, building['depth'])))
            for name, centre, scale in checks:
                node = nodes[name]
                expected_pos = (x, centre, -y)  # Blender Z-up to glTF Y-up, only once.
                if 'matrix' in node or any(abs(a - b) > .002 for a, b in zip(
                        node.get('translation', (0, 0, 0)), expected_pos)):
                    raise RuntimeError('Incorrect seated transform for ' + name + ' in ' + str(path))
                if scale and any(abs(a - b) > .002 for a, b in zip(node.get('scale', (1, 1, 1)), scale)):
                    raise RuntimeError('Incorrect seated scale for ' + name + ' in ' + str(path))
    if set(doc['scenes'][doc.get('scene', 0)]['nodes']) != set(range(len(doc['nodes']))):
        raise RuntimeError('Unexpected parent transform in ' + str(path))

    def triangles(mesh):
        return sum(doc['accessors'][part['indices']]['count'] // 3 if 'indices' in part
                   else doc['accessors'][part['attributes']['POSITION']]['count'] // 3
                   for part in mesh['primitives'])

    count = sum(triangles(doc['meshes'][node['mesh']]) for node in doc['nodes'] if 'mesh' in node)
    if args.merge_buildings:
        # Merged mode: the reference is the authored merged mesh itself (degenerate faces were
        # dropped at collection time, so export must reproduce it one-for-one). Per-building
        # seated transforms were verified by the identical make_cell code path in unmerged runs.
        reference = sum(len(polygon.vertices) - 2 for polygon in objects[0].data.polygons)
    else:
        # Asset placements also carry a 12-triangle buried foundation box; massing prisms do not.
        reference = sum(expected_triangles[b['id']][level] + (12 if entry and not (args.lod1_footprint_only and level == 'lod1') else 0)
                          for b, entry, _, _ in group)
    if count != reference:
        raise RuntimeError(f'{path.name}: {count} triangles, expected {reference}')
    meshes = {node['mesh'] for node in doc['nodes'] if 'mesh' in node}
    # Per-building material remapping is an object-level material link, so a cell whose instances
    # were recoloured exports one mesh variant per building; shared mesh data is not asserted here.
    # `uniqueMeshes` is reported so a regression towards N meshes for N nodes stays visible.
    return {'uri': path.name, 'bytes': path.stat().st_size, 'triangles': count,
            'meshNodes': len(names), 'uniqueMeshes': len(meshes),
            'merged': bool(args.merge_buildings)}


# ---------------------------------------------------------------- spatial cells
# Partition by positionENU into a --grid x --grid grid, then group cells into --coarse-grid parents,
# exactly like the rectangular exporter so the same verifier accepts the tree.
xs = [b['positionENU'][0] for b, _, _, _ in placements]
ys = [b['positionENU'][1] for b, _, _, _ in placements]
min_x, max_x = min(xs), max(xs)
min_y, max_y = min(ys), max(ys)
span_x = max(max_x - min_x, 1e-6)
span_y = max(max_y - min_y, 1e-6)
cells = {}
for placement in placements:
    x, y, _ = placement[0]['positionENU']
    ix = min(args.grid - 1, int((x - min_x) / span_x * args.grid))
    iy = min(args.grid - 1, int((y - min_y) / span_y * args.grid))
    cells.setdefault(f'cell-{iy:02d}-{ix:02d}', []).append(placement)
cells = {name: group for name, group in sorted(cells.items()) if group}
if len(cells) < 2:
    raise RuntimeError('Not enough geographically separate building groups to demonstrate spatial LOD')

hierarchical = bool(args.coarse_grid)
coarse_groups, fine_coarse = {}, {}
if hierarchical:
    for cell_id, group in cells.items():
        x, y, _ = group[0][0]['positionENU']
        ix = min(args.coarse_grid - 1, int((x - min_x) / span_x * args.coarse_grid))
        iy = min(args.coarse_grid - 1, int((y - min_y) / span_y * args.coarse_grid))
        coarse_id = f'coarse-{iy:02d}-{ix:02d}'
        fine_coarse[cell_id] = coarse_id
        coarse_groups.setdefault(coarse_id, []).extend(group)

fine_tiles, manifest_cells, all_bounds = {}, [], []
for cell_id, group in cells.items():
    exports, bounds_per_level = {}, []
    for level in ('lod1', 'lod0'):
        collection, objects, bounds = make_cell(group, level)
        path = OUT / f'{cell_id}-{level}.glb'
        exports[level] = export_cell(path, objects, group, level)
        bounds_per_level.append(bounds)
        bpy.data.collections.remove(collection, do_unlink=True)
        for obj in objects:
            bpy.data.objects.remove(obj, do_unlink=True)
    bounds = [[min(b[0][axis] for b in bounds_per_level) for axis in range(3)],
              [max(b[1][axis] for b in bounds_per_level) for axis in range(3)]]
    all_bounds.append(bounds)
    volume = {'box': box(bounds)}
    if hierarchical:
        fine_tiles[cell_id] = {'boundingVolume': volume, 'geometricError': 0,
                               'content': {'uri': exports['lod0']['uri']}}
    else:
        fine_tiles[cell_id] = {'boundingVolume': volume, 'geometricError': 8,
                               'refine': 'REPLACE', 'content': {'uri': exports['lod1']['uri']},
                               'children': [{'boundingVolume': volume, 'geometricError': 0,
                                             'content': {'uri': exports['lod0']['uri']}}]}
    manifest_cells.append({'id': cell_id, 'count': len(group),
                           'buildingIds': [b['id'] for b, _, _, _ in group],
                           'boundsENU': {'min': bounds[0], 'max': bounds[1]}, **exports})

if hierarchical:
    cell_tiles, manifest_coarse = [], []
    for coarse_id, group in coarse_groups.items():
        collection, objects, bounds = make_cell(group, 'lod1')
        path = OUT / f'{coarse_id}-lod1.glb'
        coarse_export = export_cell(path, objects, group, 'lod1')
        bpy.data.collections.remove(collection, do_unlink=True)
        for obj in objects:
            bpy.data.objects.remove(obj, do_unlink=True)
        children = [fine_tiles[cell_id] for cell_id, parent in fine_coarse.items() if parent == coarse_id]
        cell_tiles.append({'boundingVolume': {'box': box(bounds)}, 'geometricError': 16,
                           'refine': 'REPLACE', 'content': {'uri': coarse_export['uri']},
                           'children': children})
        manifest_coarse.append({'id': coarse_id, 'count': len(group), **coarse_export})
else:
    cell_tiles = list(fine_tiles.values())
    manifest_coarse = []

root_bounds = [[min(b[0][axis] for b in all_bounds) for axis in range(3)],
               [max(b[1][axis] for b in all_bounds) for axis in range(3)]]
# The explicit root must be renderable to keep Cesium traversing its spatial children.
anchor_json = json.dumps({'asset': {'version': '2.0'}, 'scene': 0, 'scenes': [{'nodes': []}]},
                         separators=(',', ':')).encode('utf-8')
anchor_json += b' ' * ((-len(anchor_json)) % 4)
anchor = b'glTF' + struct.pack('<II', 2, 20 + len(anchor_json))
anchor += struct.pack('<I', len(anchor_json)) + b'JSON' + anchor_json
(OUT / 'root-anchor.glb').write_bytes(anchor)
tileset = {'asset': {'version': '1.1', 'generator': 'Blender CAD/BIM massing LOD'},
           'geometricError': 800,
           'root': {'transform': plan['transform'], 'boundingVolume': {'box': box(root_bounds)},
                    'geometricError': 200, 'refine': 'REPLACE',
                    'content': {'uri': 'root-anchor.glb'}, 'children': cell_tiles}}
lod1_records = manifest_coarse if hierarchical else [cell['lod1'] for cell in manifest_cells]
lod0_records = [cell['lod0'] for cell in manifest_cells]
manifest = {'scope': plan.get('scope', 'surveyed CAD/BIM massing sample; not full-city reconstruction'),
            'origin': plan['origin'], 'notice': plan.get('notice'),
             'provenance': {'source': plan.get('source'), 'plan': PLAN.name, 'planPath': str(PLAN),
                             'nearArchitecture': ARCHITECTURE, 'massingFacadeProfile': FACADE_PROFILE,
                              'compactFacadeBuildingIds': sorted(COMPACT_BUILDING_IDS) if COMPACT_BUILDING_IDS is not None else None,
                              'massingFacadeLayout': 'physical-windows-entrance-v2' if FACADE_PROFILE == 'compact' else 'legacy-panels',
                              'compactWindowTexture': ({key: WINDOW_TEXTURE[key] for key in
                                  ('sha256', 'source', 'rasterBytes', 'rgba8FullMipEstimate')}
                                  if WINDOW_TEXTURE else None),
                              'shellOnlyAttachmentIds': [b['id'] for b, _, _, _ in placements
                                                         if b.get('massingFacadeScope') == 'shell-only'],
                              'reviewedPartyWallIds': sorted(PARTY_WALL_MASKS),
                              'partyWallReview': plan.get('partyWallReview'),
                              'compositionReview': plan.get('compositionReview'),
                              'estimatedMassingStoreyBands': BANDS,
                              'satelliteRoofTextures': len(ROOF_TEXTURES),
                             'facadeFidelity': 'procedural estimate; not facade photography or surveyed BIM',
                            'gisFootprints': True, 'lod1FootprintOnly': bool(args.lod1_footprint_only),
                            'facadeDetails': True, 'roadsProcedural': False, 'closedRoofCaps': True},
            'reduction': {'decimateLod0': args.decimate_lod0, 'decimateLod1': args.decimate_lod1,
                          'draco': bool(args.draco)},
            'count': len(placements), 'cellCount': len(cells), 'cellIds': list(cells),
            'hierarchy': {'enabled': hierarchical, 'coarseGrid': args.coarse_grid if hierarchical else None,
                          'coarseCellCount': len(manifest_coarse)},
            'excludedUnfittedIds': rejected,
            'boundsENU': {'min': root_bounds[0], 'max': root_bounds[1]},
            'totals': {level: {'bytes': sum(item['bytes'] for item in (lod1_records if level == 'lod1' else lod0_records)),
                                'triangles': sum(item['triangles'] for item in (lod1_records if level == 'lod1' else lod0_records))}
                         for level in ('lod1', 'lod0')},
            'coarseCells': manifest_coarse,
            'cells': manifest_cells}
(OUT / 'tileset.json').write_text(json.dumps(tileset, indent=2), encoding='utf-8')
(OUT / 'manifest.json').write_text(json.dumps(manifest, indent=2), encoding='utf-8')
print('CAD_BIM_LOD', json.dumps({'count': manifest['count'], 'cells': manifest['cellCount'],
                                 'coarse': len(manifest_coarse), 'totals': manifest['totals'],
                                 'output': str(OUT)}))