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('--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 {}
ARCHITECTURE = args.near_architecture


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
    # Near architecture gets a denser, still bounded facade rhythm; coarse LOD stays intentionally sparse.
    return max(1, min(16 if ARCHITECTURE else 3, int(edge_length / (2.8 if ARCHITECTURE else 3.5))))


def massing_triangles(footprint, floors, level):
    """Reference triangles for a massing prism plus lightweight facade details."""
    base = prism_triangles(len(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 * count * floors
    architectural = (windows * 4 + 2 * len(footprint) * max(0, floors - 1)) if ARCHITECTURE and level == 'lod0' else 0
    return base + windows + architectural + 12  # bounded reveals/bands and one closed entrance


def poly_prism(footprint, bottom, top, wall_material, roof_material, anchor, floors, level,
               window_material, door_material):
    """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.
    """
    # 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):
        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 = min(length * 0.22, max(0.45, length / 3.5 * width_ratio))
        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)
        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)

    # One or a few windows per wall and floor, with a dark blue material that cannot read as a flat cube.
    for i in range(n):
        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)):
            z0 = floor_base + floor * floor_height + floor_height * 0.34
            z1 = min(top - 0.08, z0 + floor_height * 0.36)
            for window in range(count):
                t = (window + 1) / (count + 1)
                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 ARCHITECTURE and level == 'lod0':
                    add_window_reveal(wa, wb, z0, z1)
        if ARCHITECTURE and level == 'lod0':
            for floor in range(1, max(1, floors)):
                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)

    # One door on the longest wall, visible from at least one approach direction.
    longest = max(range(n), key=lambda i: math.hypot(footprint[(i + 1) % n][0] - footprint[i][0],
                                                      footprint[(i + 1) % n][1] - footprint[i][1]))
    edge_a = footprint[longest]
    edge_b = footprint[(longest + 1) % n]
    door_h = min(2.2, max(1.8, floor_height * 0.76))
    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
    edit = bmesh.new()
    edit.from_mesh(mesh)
    bmesh.ops.recalc_face_normals(edit, faces=edit.faces)
    edit.to_mesh(mesh)
    edit.free()
    mesh.update()
    for material in (wall_material, roof_material, window_material, door_material):
        mesh.materials.append(material)
    mesh.polygons.foreach_set('material_index', materials)
    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.
    return {'lod0': (outline, massing_triangles(outline, building['floors'], 'lod0')),
            'lod1': (outline, massing_triangles(outline, building['floors'], 'lod1'))}


placements, rejected, expected_triangles = [], [], {}
library_ids = {entry['id'] for entry in entries}
for building in plan['buildings'][:args.limit or None]:
    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 = 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')

# ---------------------------------------------------------------- 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] = 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: (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 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.new(name='UVMap')
    for loop in obj.data.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 = build_geometry(building)[level]
            outline = geometry[0]
            wall, roof, window, door = massing_tone(building['id'])
            if level == 'lod0':
                roof = sourced_roof(building['id'], roof)
            top_z = top + building['heightM']
            obj = 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)
            # 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)
                obj = 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)
                if prism_volume_sign(obj) <= 0:
                    raise RuntimeError(building['id'] + ' massing prism is inside-out or flat')
            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 = {}
    for obj in objects:
        evaluated = obj.evaluated_get(depsgraph)
        mesh = evaluated.to_mesh()
        matrix = obj.matrix_world.copy()
        base = len(vertices)
        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)
            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:
        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)
    bmesh.ops.recalc_face_normals(edit, faces=edit.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, '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)}))