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

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')
args = parser.parse_args(sys.argv[sys.argv.index('--') + 1:] if '--' in sys.argv else [])
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')


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):
    """Winding fix: surveyed rings arrive CW or CCW, prisms need one known orientation."""
    return list(points) if ring_area(points) > 0 else list(reversed(points))


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 (2 * (n - 2)) plus n wall quads (2n)."""
    return 4 * n - 4


def poly_prism(footprint, bottom, top, wall_material, roof_material, anchor):
    """Build one extruded prism; the ring must already be CCW.

    Vertices are stored relative to `anchor` (the plan's positionENU) and the object is moved to it.
    That keeps local coordinates small — which the volume self-check needs — and makes the exported
    node translation equal the planned position, as the cell verification expects.
    """
    bpy.ops.object.select_all(action='DESELECT')
    bpy.ops.mesh.primitive_cube_add(size=1)
    obj = bpy.context.active_object
    mesh = obj.data
    mesh.clear_geometry()
    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]
    # Exactly n wall quads plus two caps. Adding the walls a second time (in reverse) would make
    # every wall edge belong to four faces, so the mesh is non-manifold and recalc_face_normals
    # cannot orient it — which silently produced a negative volume for thin outlines.
    faces = [(i, (i + 1) % n, n + (i + 1) % n, n + i) for i in range(n)]
    faces += [tuple(range(n)), tuple(reversed(range(n, 2 * n)))]
    mesh.from_pydata(vertices, [], faces)
    obj.location = anchor
    # Blender settles the normals for this closed manifold surface; prism_volume_sign in make_cell
    # is the independent outward-normal check on the result.
    edit = bmesh.new()
    edit.from_mesh(mesh)
    bmesh.ops.recalc_face_normals(edit, faces=edit.faces)
    edit.to_mesh(mesh)
    edit.free()
    mesh.update()
    mesh.materials.append(wall_material)
    mesh.materials.append(roof_material)
    # Face order is n walls followed by the two caps; slot 0 is the wall tone, slot 1 the roof.
    mesh.polygons.foreach_set('material_index', [0] * n + [1, 1])
    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':
        # 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')
    hull = convex_hull(outline)
    if len(hull) < 3:
        raise RuntimeError(building['id'] + ' LOD1 hull collapsed')
    return {'lod0': (outline, prism_triangles(len(outline))),
            'lod1': (hull, prism_triangles(len(hull)))}


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')
for building, entry, _, _ in placements:
    if entry:
        expected_triangles[building['id']] = {
            level: (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).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)]])


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:
            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 = massing_tone(building['id'])
            top_z = top + building['heightM']
            obj = poly_prism(outline, bottom, top_z, wall, roof, building['positionENU'])
            # 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'])
                if prism_volume_sign(obj) <= 0:
                    raise RuntimeError(building['id'] + ' massing prism is inside-out or flat')
            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:
            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)
    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 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]
    bpy.ops.export_scene.gltf(filepath=str(path), export_format='GLB', use_selection=True,
                              export_yup=True, export_materials='EXPORT', export_extras=True)
    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.
    expected_names = []
    for building, entry, _, _ in group:
        expected_names.append(building['id'])
        if entry:
            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))
    if set(doc['scenes'][doc.get('scene', 0)]['nodes']) != set(range(len(doc['nodes']))):
        raise RuntimeError('Unexpected parent transform in ' + 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']
        checks = [(building['id'], top, (sx, 1, sy) if entry else None)]
        if entry:
            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))

    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)
    # Asset placements also carry a 12-triangle buried foundation box; massing prisms do not.
    reference = sum(expected_triangles[b['id']][level] + (12 if entry 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)}


# ---------------------------------------------------------------- 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'),
            '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)}))