"""Build a small, deterministic architectural asset library with Blender 5.x.

Run from the workspace (no add-ons, textures, downloads, or app server required):
  blender --background --factory-startup --python-exit-code 1 --python tools/city-massing/blender-house-library.py
  blender --background --factory-startup --python-exit-code 1 --python tools/city-massing/blender-house-library.py -- --samples 48

Only this script's sibling house-library directory is written by default.
Authoring coordinates: metres, Z up, +Y front, footprint centred, ground at Z=0.
Each runtime object is one mesh composed of individually closed solids. Detail
solids deliberately intersect; these are instancing assets, not boolean-unioned
construction CAD. Optional foundations are separate, strictly below ground.
"""

import argparse
import json
import math
import random
import struct
import sys
import zlib
from pathlib import Path

import bpy
import bmesh
from mathutils import Vector

SEED = 20260918
RNG = random.Random(SEED)
PALETTE = {}
MATERIALS = []
CUBE_FACES = [(0, 3, 2, 1), (0, 1, 5, 4), (1, 2, 6, 5),
              (2, 3, 7, 6), (3, 0, 4, 7), (4, 5, 6, 7)]


def arguments():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--output', type=Path, default=Path(__file__).resolve().parent / 'house-library')
    parser.add_argument('--samples', type=int, default=40)
    parser.add_argument('--skip-render', action='store_true')
    return parser.parse_args(sys.argv[sys.argv.index('--') + 1:] if '--' in sys.argv else [])


def linear(v):
    return v / 12.92 if v <= 0.04045 else ((v + 0.055) / 1.055) ** 2.4


def material(name, color, roughness, metallic=0.0):
    rgb = tuple(linear(int(color[i:i + 2], 16) / 255) for i in (0, 2, 4))
    mat = bpy.data.materials.new(name)
    mat.use_nodes = True
    mat.diffuse_color = (*rgb, 1.0)
    mat.use_backface_culling = True
    shader = mat.node_tree.nodes.get('Principled BSDF')
    shader.inputs['Base Color'].default_value = (*rgb, 1.0)
    shader.inputs['Roughness'].default_value = roughness
    shader.inputs['Metallic'].default_value = metallic
    mat['palette_hex_srgb'] = '#' + color
    PALETTE[name] = len(MATERIALS)
    MATERIALS.append(mat)


def make_materials():
    for entry in [
        ('Plaster_Cream', 'E4DAC4', 0.87),
        ('Plaster_Ivory', 'F2ECDD', 0.86),
        ('Plaster_Sand', 'C2AF96', 0.88),
        ('Reveal_Shadow', '8D8B80', 0.91),
        ('Stone_Base', 'A1A39C', 0.89),
        ('Concrete_Warm', 'C8C7BD', 0.90),
        ('Roof_Terracotta', '98634A', 0.84),
        ('Roof_Terra_Light', 'AC7758', 0.86),
        ('Roof_Slate', '596264', 0.82),
        ('Roof_Slate_Light', '697272', 0.83),
        ('Frame_Ivory', 'DEDCD0', 0.49, 0.15),
        ('Metal_Charcoal', '475052', 0.48, 0.48),
        ('Glass_Blue', '617E88', 0.24, 0.20),
        ('Glass_Soft', '91A3A3', 0.30, 0.12),
        ('Glass_Dark', '405F6D', 0.24, 0.20),
        ('Door_Walnut', '715747', 0.68),
        ('Roof_Membrane', '777F7D', 0.93),
    ]:
        material(*entry)


def cube_vertices(center, size):
    x, y, z = center
    a, b, c = (s / 2 for s in size)
    return [(x - a, y - b, z - c), (x + a, y - b, z - c),
            (x + a, y + b, z - c), (x - a, y + b, z - c),
            (x - a, y - b, z + c), (x + a, y - b, z + c),
            (x + a, y + b, z + c), (x - a, y + b, z + c)]


def face_point(side, u, depth, z, center=(0, 0)):
    x, y = center
    if side == 'front':
        return (x + u, y + depth, z)
    if side == 'right':
        return (x + depth, y - u, z)
    if side == 'back':
        return (x - u, y - depth, z)
    return (x - depth, y + u, z)


class Geometry:
    """Accumulate closed pieces, tagging dispensable small detail for LOD1."""

    def __init__(self):
        self.pieces = []

    def solid(self, vertices, faces, mat, tiny=False):
        indices = [PALETTE[m] for m in mat] if isinstance(mat, list) else [PALETTE[mat]] * len(faces)
        self.pieces.append((vertices, faces, indices, tiny))

    def box(self, center, size, mat, tiny=False):
        assert min(size) > 0, (center, size)
        self.solid(cube_vertices(center, size), CUBE_FACES, mat, tiny)

    def fbox(self, side, u, depth, z, width, thickness, height, mat, center=(0, 0), tiny=False):
        vertices = [face_point(side, *v, center) for v in cube_vertices((u, depth, z), (width, thickness, height))]
        self.solid(vertices, CUBE_FACES, mat, tiny)

    def ring(self, side, u, depth, z, width, height, border, thickness, mat, center=(0, 0)):
        vertices = []
        for d in (depth - thickness / 2, depth + thickness / 2):
            for w, h in ((width, height), (width - 2 * border, height - 2 * border)):
                vertices += [face_point(side, u + a * w / 2, d, z + b * h / 2, center)
                             for a, b in ((-1, -1), (1, -1), (1, 1), (-1, 1))]
        faces = []
        for i in range(4):
            j = (i + 1) % 4
            faces += [(i, j, 8 + j, 8 + i), (4 + i, 12 + i, 12 + j, 4 + j),
                      (i, 4 + i, 4 + j, j), (8 + i, 8 + j, 12 + j, 12 + i)]
        self.solid(vertices, faces, mat)

    def object(self, name, collection, lod=0):
        vertices, faces, indices = [], [], []
        for verts, polys, mats, tiny in self.pieces:
            if lod and tiny:
                continue
            offset = len(vertices)
            vertices.extend(verts)
            faces.extend(tuple(offset + i for i in poly) for poly in polys)
            indices.extend(mats)
        mesh = bpy.data.meshes.new(name + '_Mesh')
        mesh.from_pydata(vertices, [], faces)
        mesh.update()
        for mat in MATERIALS:
            mesh.materials.append(mat)
        for poly, index in zip(mesh.polygons, indices):
            poly.material_index = index
        bm = bmesh.new()
        bm.from_mesh(mesh)
        bmesh.ops.recalc_face_normals(bm, faces=list(bm.faces))
        # Coalesce the coplanar window-wall grid, preserving reveal geometry.
        bmesh.ops.dissolve_limit(bm, angle_limit=0.00001, verts=list(bm.verts),
                                 edges=list(bm.edges), delimit={'MATERIAL'}, use_dissolve_boundaries=False)
        bm.to_mesh(mesh)
        bm.free()
        mesh.update()
        obj = bpy.data.objects.new(name, mesh)
        collection.objects.link(obj)
        return obj


def pocket_wall(geo, side, width, depth, height, openings, mat, center=(0, 0), base=0):
    """A watertight wall with actual recessed window wells (no alpha planes)."""
    thickness, recess = 0.32, 0.18
    xs = sorted(set([-width / 2, width / 2] + [v for o in openings for v in (o['u'] - o['w'] / 2, o['u'] + o['w'] / 2)]))
    zs = sorted(set([base, base + height] + [v for o in openings for v in (o['z'] - o['h'] / 2, o['z'] + o['h'] / 2)]))
    assert xs[0] >= -width / 2 and xs[-1] <= width / 2
    assert zs[0] >= base and zs[-1] <= base + height
    vertices, faces, mats, cache = [], [], [], {}

    def vertex(u, d, z):
        key = (round(u, 8), round(d, 8), round(z, 8))
        if key not in cache:
            cache[key] = len(vertices)
            vertices.append(face_point(side, u, d, z, center))
        return cache[key]

    def quad(points, material_name=mat):
        faces.append(tuple(vertex(*p) for p in points))
        mats.append(material_name)

    depths = []
    for k in range(len(zs) - 1):
        row = []
        for j in range(len(xs) - 1):
            u, z = (xs[j] + xs[j + 1]) / 2, (zs[k] + zs[k + 1]) / 2
            inside = any(abs(u - o['u']) < o['w'] / 2 - 1e-8 and abs(z - o['z']) < o['h'] / 2 - 1e-8 for o in openings)
            row.append(depth - recess if inside else depth)
        depths.append(row)
    inner = depth - thickness
    for k, row in enumerate(depths):
        z0, z1 = zs[k:k + 2]
        for j, d in enumerate(row):
            u0, u1 = xs[j:j + 2]
            quad([(u0, d, z0), (u0, d, z1), (u1, d, z1), (u1, d, z0)], 'Reveal_Shadow' if d < depth else mat)
            quad([(u0, inner, z0), (u1, inner, z0), (u1, inner, z1), (u0, inner, z1)])
            if j == 0:
                quad([(u0, inner, z0), (u0, inner, z1), (u0, d, z1), (u0, d, z0)])
            if j == len(row) - 1:
                quad([(u1, inner, z0), (u1, d, z0), (u1, d, z1), (u1, inner, z1)])
            if k == 0:
                quad([(u0, inner, z0), (u0, d, z0), (u1, d, z0), (u1, inner, z0)])
            if k == len(depths) - 1:
                quad([(u0, inner, z1), (u1, inner, z1), (u1, d, z1), (u0, d, z1)])
            if j and d != row[j - 1]:
                prev = row[j - 1]
                quad([(u0, d, z0), (u0, prev, z0), (u0, prev, z1), (u0, d, z1)], 'Reveal_Shadow')
            if k and d != depths[k - 1][j]:
                prev = depths[k - 1][j]
                quad([(u0, d, z0), (u1, d, z0), (u1, prev, z0), (u0, prev, z0)], 'Reveal_Shadow')
    geo.solid(vertices, faces, mats)


def opening(u, z, w=1.8, h=1.65, door=False):
    return {'u': u, 'z': z, 'w': w, 'h': h, 'door': door}


def window(geo, side, depth, op, center=(0, 0)):
    u, z, w, h = (op[k] for k in ('u', 'z', 'w', 'h'))
    is_door = op.get('door', False)
    frame = 'Metal_Charcoal' if is_door else 'Frame_Ivory'
    geo.ring(side, u, depth + 0.008, z, w + 0.10, h + 0.10, 0.095, 0.11, frame, center)
    cols = 2 if w >= 1.35 else 1
    pane_w = (w - 0.12) / cols
    for c in range(cols):
        pane_u = u + (c - (cols - 1) / 2) * pane_w
        glass = RNG.choice(['Glass_Blue', 'Glass_Blue', 'Glass_Soft', 'Glass_Dark'])
        geo.fbox(side, pane_u, depth - 0.095, z, pane_w - 0.035, 0.028, h - 0.12, glass, center)
    if cols > 1:
        geo.fbox(side, u, depth - 0.016, z, 0.052, 0.068, h - 0.09, frame, center)
    geo.fbox(side, u, depth - 0.012, z + h * 0.20, w - 0.08, 0.071, 0.048, frame, center, tiny=True)
    if is_door:
        geo.fbox(side, u, depth - 0.047, z - h * 0.29, w - 0.12, 0.056, h * 0.30, 'Door_Walnut', center)
        for sign in (-1, 1):
            geo.fbox(side, u + sign * 0.10, depth + 0.055, z - 0.04, 0.027, 0.065, 0.27, 'Metal_Charcoal', center, tiny=True)
    else:
        geo.fbox(side, u, depth + 0.07, z - h / 2 - 0.08, w + 0.24, 0.29, 0.095, 'Concrete_Warm', center)
        geo.fbox(side, u, depth + 0.023, z + h / 2 + 0.085, w + 0.18, 0.14, 0.07, 'Plaster_Ivory', center, tiny=True)


def building(parts, width, depth, height, openings, center=(0, 0), base=0, mat='Plaster_Cream'):
    core, glazing = parts['Structure'], parts['Windows']
    core.box((*center, base + height / 2), (width - 0.64, depth - 0.64, height), mat)
    for side in ('front', 'back', 'right', 'left'):
        span, d = (width, depth / 2) if side in ('front', 'back') else (depth - 0.64, width / 2)
        ops = openings.get(side, [])
        pocket_wall(core, side, span, d, height, ops, mat, center, base)
        for op in ops:
            window(glazing, side, d, op, center)


def belt(geo, width, depth, z, height=0.12, mat='Plaster_Ivory', center=(0, 0), front_gaps=()):
    for side in ('front', 'back', 'left', 'right'):
        span, d = (width, depth / 2) if side in ('front', 'back') else (depth, width / 2)
        cursor = -(span + 0.10) / 2
        gaps = sorted(front_gaps) if side == 'front' else []
        for left, right in gaps + [((span + 0.10) / 2, (span + 0.10) / 2)]:
            if left > cursor:
                geo.fbox(side, (cursor + left) / 2, d + 0.025, z, left - cursor, 0.13, height, mat, center)
            cursor = right


def rail(geo, side, u, depth, floor, width, height=0.95, center=(0, 0), spacing=0.22):
    for z, thick in ((floor + 0.13, 0.048), (floor + height, 0.065)):
        geo.fbox(side, u, depth, z, width, 0.058, thick, 'Metal_Charcoal', center)
    count = max(2, math.ceil(width / spacing))
    for i in range(count + 1):
        geo.fbox(side, u - width / 2 + width * i / count, depth, floor + height / 2,
                 0.027, 0.034, height - 0.055, 'Metal_Charcoal', center, tiny=(i % 2 == 1))


def balcony(geo, side, u, depth, floor, width, reach=1.0, center=(0, 0), parapet=False, spacing=0.22):
    geo.fbox(side, u, depth + reach / 2 - 0.025, floor - 0.08, width, reach + 0.12, 0.16, 'Concrete_Warm', center)
    outer = depth + reach - 0.06
    if parapet:
        geo.fbox(side, u, outer, floor + 0.17, width, 0.10, 0.34, 'Plaster_Ivory', center)
    rail(geo, side, u, outer, floor, width - 0.10, center=center, spacing=spacing)
    # Side rails are also true, closed rectangular bars.
    for sign in (-1, 1):
        for z in (floor + 0.13, floor + 0.95):
            geo.fbox(side, u + sign * (width / 2 - 0.045), depth + reach / 2, z, 0.055, reach - 0.10, 0.055, 'Metal_Charcoal', center)
        n = max(2, math.ceil(reach / spacing))
        for i in range(n):
            geo.fbox(side, u + sign * (width / 2 - 0.045), depth + (i + 0.5) * reach / n,
                     floor + 0.48, 0.034, 0.027, 0.88, 'Metal_Charcoal', center, tiny=(i % 2 == 1))


def gable(geo, width, depth, eave, rise, mat, center=(0, 0), roof=False):
    """Extruded closed gable; roof uses a six-sided closed thick profile."""
    x, y = center
    if roof:
        t = 0.18
        profile = [(-depth / 2, eave - t), (0, eave + rise - t), (depth / 2, eave - t),
                   (depth / 2, eave), (0, eave + rise), (-depth / 2, eave)]
    else:
        profile = [(-depth / 2, eave), (depth / 2, eave), (0, eave + rise)]
    vertices = [(x + a, y + b, z) for a in (-width / 2, width / 2) for b, z in profile]
    n = len(profile)
    faces = [tuple(reversed(range(n))), tuple(range(n, n * 2))]
    faces += [(i, (i + 1) % n, (i + 1) % n + n, i + n) for i in range(n)]
    geo.solid(vertices, faces, mat)


def pitched_roof(parts, width, depth, eave, rise, slate=False, center=(0, 0)):
    geo = parts['Roof']
    mat = 'Roof_Slate' if slate else 'Roof_Terracotta'
    light = 'Roof_Slate_Light' if slate else 'Roof_Terra_Light'
    rw, rd = width + 0.74, depth + 0.72
    # Match the gable to the underside of the overhanging roof, then overlap by
    # 25 mm. A flat soffit closes the small eave gap above the vertical walls.
    gable_base = eave - 0.13 + rise * (1 - depth / rd) + 0.025
    gable(parts['Structure'], width, depth, gable_base, rise * depth / rd, 'Plaster_Ivory', center)
    parts['Structure'].box((*center, (eave - 0.03 + gable_base) / 2),
                           (width, depth, gable_base - eave + 0.07), 'Plaster_Ivory')
    gable(geo, rw, rd, eave + 0.05, rise, mat, center, roof=True)
    cx, cy = center
    # Low-profile ridge cap and tile courses; the roof shell itself has no holes.
    geo.box((cx, cy, eave + rise + 0.073), (rw + 0.08, 0.21, 0.13), light)
    slope = rise / (rd / 2)
    course_count = math.ceil((rd / 2) / 0.34)
    for sign in (-1, 1):
        for i in range(1, course_count):
            d = i * rd / 2 / course_count
            z = eave + 0.05 + rise - d * slope + 0.008
            geo.box((cx, cy + sign * d, z), (rw, 0.045, 0.025), light, tiny=True)
        # Raised tile ribs give directional texture without huge geometry counts.
        for k in range(math.ceil(rw / 0.48)):
            x = cx - rw / 2 + (k + 0.5) * rw / math.ceil(rw / 0.48)
            y0, y1 = 0.09 * sign, (rd / 2 - 0.06) * sign
            z0 = eave + 0.05 + rise - abs(y0) * slope + 0.024
            z1 = eave + 0.05 + rise - abs(y1) * slope + 0.024
            beam(geo, (x, cy + y0, z0), (x, cy + y1, z1), 0.018, 0.027, light, tiny=True)
        geo.box((cx, cy + sign * (rd / 2 - 0.035), eave - 0.01), (rw, 0.095, 0.11), 'Metal_Charcoal')


def beam(geo, start, end, width, height, mat, tiny=False):
    a, b = Vector(start), Vector(end)
    direction = (b - a).normalized()
    axis = Vector((1, 0, 0)) if abs(direction.x) < 0.8 else Vector((0, 1, 0))
    cross = direction.cross(axis).normalized()
    side = cross.cross(direction).normalized()
    center = (a + b) / 2
    vertices = [tuple(center + side * v[0] + cross * v[1] + direction * v[2])
                for v in cube_vertices((0, 0, 0), (width, height, (b - a).length))]
    geo.solid(vertices, CUBE_FACES, mat, tiny)


def canopy(geo, u, depth, z, width=2.4, reach=1.0, center=(0, 0)):
    geo.fbox('front', u, depth + reach / 2, z, width, reach + 0.12, 0.11, 'Metal_Charcoal', center)
    for sign in (-1, 1):
        beam(geo, (center[0] + u + sign * width * 0.36, center[1] + depth + 0.02, z - 0.55),
             (center[0] + u + sign * width * 0.36, center[1] + depth + reach * 0.8, z - 0.07),
             0.055, 0.055, 'Metal_Charcoal')


def ac_unit(geo, side, u, depth, z, center=(0, 0)):
    geo.fbox(side, u, depth + 0.19, z, 0.72, 0.36, 0.48, 'Plaster_Ivory', center, tiny=True)
    for k in range(4):
        geo.fbox(side, u - 0.11, depth + 0.38, z - 0.15 + k * 0.10, 0.41, 0.022, 0.022, 'Metal_Charcoal', center, tiny=True)


def new_parts():
    return {name: Geometry() for name in ('Structure', 'Windows', 'Balconies', 'Roof', 'Details')}


def make_house(floors=3):
    p = new_parts()
    height = floors * 3.05
    ops = {side: [] for side in ('front', 'back', 'right', 'left')}
    for floor in range(floors):
        base = floor * 3.05
        ops['front'] += [opening(-3.85, base + 1.66, 2.05),
                         opening(-0.45, base + (1.20 if floor == 0 else 1.69), 1.50, 2.30 if floor == 0 else 1.65, floor == 0),
                         opening(3.45, base + (1.67 if floor == 0 else 1.27), 2.15, 1.65 if floor == 0 else 2.24)]
        ops['back'] += [opening(u, base + 1.67, 1.65) for u in (-3.85, 0, 3.85)]
        for side in ('left', 'right'):
            ops[side] += [opening(u, base + 1.67, 1.75) for u in (-2.40, 2.40)]
    building(p, 12, 9, height, ops)
    for f in range(1, floors + 1):
        belt(p['Details'], 12, 9, min(f * 3.05, height - 0.02))
    belt(p['Details'], 12, 9, 0.20, 0.40, 'Stone_Base', front_gaps=[(-1.27, 0.37)])
    for x in (-5.93, 5.93):
        p['Details'].box((x, 4.535, height / 2 + 0.18), (0.20, 0.13, height - 0.42), 'Plaster_Ivory')
    for f in range(1, floors):
        balcony(p['Balconies'], 'front', 3.45, 4.50, f * 3.05 + 0.015, 3.05, 1.12)
    canopy(p['Details'], -0.45, 4.5, 2.65, 2.45, 1.15)
    pitched_roof(p, 12, 9, height, 2.05)
    for x in (-5.67, 5.67):
        p['Details'].box((x, -4.58, height / 2), (0.085, 0.085, height - 0.15), 'Metal_Charcoal', tiny=True)
    for f in range(1, floors):
        ac_unit(p['Details'], 'back', 1.40, 4.5, f * 3.05 + 0.65)
    # Two louvred attic vents on the gable ends.
    for side in ('left', 'right'):
        p['Details'].fbox(side, 0, 6.015, height + 0.65, 0.65, 0.07, 0.48, 'Metal_Charcoal')
        for k in range(4):
            p['Details'].fbox(side, 0, 6.06, height + 0.50 + k * 0.10, 0.66, 0.03, 0.034, 'Frame_Ivory', tiny=True)
    return p


def make_slab(floors=6):
    p = new_parts()
    height = floors * 3.05
    ops = {side: [] for side in ('front', 'back', 'right', 'left')}
    for floor in range(floors):
        base = floor * 3.05
        for u in (-12, -8, -2, 2, 8, 12):
            tall = floor > 0 and u in (-12, -2, 8)
            ops['front'].append(opening(u, base + (1.29 if tall else 1.70), 1.95, 2.25 if tall else 1.64))
            ops['back'].append(opening(-u, base + 1.68, 1.75, 1.62))
        for u in (-5, 5):
            ops['front'].append(opening(u, base + (1.23 if floor == 0 else 1.65), 1.18, 2.34 if floor == 0 else 1.97, floor == 0))
            ops['back'].append(opening(u, base + 1.70, 0.95, 1.70))
        for side in ('left', 'right'):
            ops[side] += [opening(u, base + 1.69, 1.75) for u in (-3.3, 0, 3.3)]
    building(p, 30, 12, height, ops, mat='Plaster_Ivory')
    belt(p['Details'], 30, 12, 0.28, 0.56, 'Stone_Base', front_gaps=[(-5.67, -4.33), (4.33, 5.67)])
    for f in range(1, floors + 1):
        belt(p['Details'], 30, 12, min(f * 3.05, height - 0.04), 0.14, 'Concrete_Warm')
    # Vertical stair strips stop and restart around actual staircase windows.
    for u in (-5, 5):
        for x in (u - 0.92, u + 0.92):
            p['Details'].fbox('front', x, 6.022, height / 2, 0.35, 0.12, height, 'Plaster_Sand')
        canopy(p['Details'], u, 6, 2.65, 2.2, 1.10)
    for f in range(1, floors):
        for u in (-10, 0, 10):
            balcony(p['Balconies'], 'front', u, 6, f * 3.05 + 0.01, 6.02, 1.02, parapet=True)
    pitched_roof(p, 30, 12, height, 1.8, slate=True)
    for u in (-10, 0, 10):
        p['Roof'].box((u, 0.9, height + 1.74), (0.68, 0.60, 0.85), 'Plaster_Sand')
        p['Roof'].box((u, 0.9, height + 2.20), (0.82, 0.74, 0.12), 'Roof_Slate')
    for f in range(floors):
        for u in (-13.5, -3.6, 6.4):
            ac_unit(p['Details'], 'back', u, 6, f * 3.05 + 0.65)
    for u in (-14.75, 14.75):
        p['Details'].fbox('back', u, 6.08, height / 2, 0.095, 0.095, height - 0.06, 'Metal_Charcoal', tiny=True)
    return p


def make_midrise(floors=10):
    p = new_parts()
    lower_floors = floors - 2
    setback_z, roof_z = lower_floors * 3.05, floors * 3.05
    ops = {side: [] for side in ('front', 'back', 'right', 'left')}
    for floor in range(lower_floors):
        base = floor * 3.05
        for u in (-7, -3.6, 0, 3.6, 7):
            stair = u == 0
            tall = floor > 0 and abs(u) == 7
            ops['front'].append(opening(u, base + (1.25 if tall or (stair and floor == 0) else 1.7),
                                        1.5 if stair else 2.10, 2.30 if tall or (stair and floor == 0) else 1.75, stair and floor == 0))
            ops['back'].append(opening(u, base + 1.7, 1.35 if stair else 1.95, 1.70))
        for side in ('left', 'right'):
            ops[side] += [opening(u, base + 1.7, 2.02, 1.75) for u in (-5.1, 0, 5.1)]
    building(p, 20, 16, setback_z, ops)
    upper = {side: [] for side in ('front', 'back', 'right', 'left')}
    for floor in range(2):
        base = setback_z + floor * 3.05
        for side in ('front', 'back'):
            upper[side] += [opening(u, base + 1.67, 2.05, 1.75) for u in (-5.8, -1.9, 1.9, 5.8)]
        for side in ('left', 'right'):
            upper[side] += [opening(u, base + 1.67, 2.00, 1.75) for u in (-4.5, 0, 4.5)]
    center = (-0.8, -0.8)
    building(p, 16.6, 14.4, 6.1, upper, center, setback_z, 'Plaster_Ivory')
    belt(p['Details'], 20, 16, 0.30, 0.60, 'Stone_Base', front_gaps=[(-0.83, 0.83)])
    for f in range(1, lower_floors + 1):
        belt(p['Details'], 20, 16, min(f * 3.05, setback_z - 0.04), 0.13, 'Concrete_Warm')
    for f in range(1, lower_floors):
        for u in (-5.3, 5.3):
            balcony(p['Balconies'], 'front', u, 8, f * 3.05 + 0.01, 5.45, 1.10, parapet=True)
    for u in (-9.9, 9.9):
        p['Details'].fbox('front', u, 8.03, setback_z / 2, 0.23, 0.16, setback_z, 'Plaster_Sand')
    for u in (-1.16, 1.16):
        p['Details'].fbox('front', u, 8.025, setback_z / 2, 0.21, 0.13, setback_z, 'Plaster_Sand')
    canopy(p['Details'], 0, 8, 2.70, 3.15, 1.22)
    # Usable setback terrace and modest rooftop plant room, all closed volumes.
    p['Roof'].box((0, 0, setback_z + 0.02), (20.15, 16.15, 0.16), 'Roof_Membrane')
    rail(p['Balconies'], 'front', 0, 8.015, setback_z + 0.10, 19.90)
    rail(p['Balconies'], 'right', 0, 10.015, setback_z + 0.10, 15.85)
    belt(p['Details'], 16.6, 14.4, setback_z + 3.05, 0.12, center=center)
    p['Roof'].box((*center, roof_z + 0.01), (16.8, 14.6, 0.20), 'Roof_Membrane')
    for side in ('front', 'back', 'left', 'right'):
        span, d = (16.85, 7.28) if side in ('front', 'back') else (14.45, 8.28)
        p['Roof'].fbox(side, 0, d, roof_z + 0.35, span, 0.17, 0.68, 'Plaster_Ivory', center)
        p['Roof'].fbox(side, 0, d, roof_z + 0.70, span + 0.08, 0.24, 0.07, 'Roof_Slate', center)
    p['Roof'].box((-2.1, -2.1, roof_z + 0.80), (5.0, 4.4, 1.55), 'Plaster_Sand')
    p['Roof'].box((-2.1, -2.1, roof_z + 1.62), (5.25, 4.65, 0.16), 'Roof_Slate')
    for f in range(lower_floors):
        for u in (-8.6, 1.4, 8.6):
            ac_unit(p['Details'], 'back', u, 8, f * 3.05 + 0.66)
    return p


def make_wide_midrise(floors=10):
    """28 x 24 m module: add facade bays, never stretch the 20 x 16 m windows."""
    p = new_parts()
    lower_floors = floors - 2
    setback_z, roof_z = lower_floors * 3.05, floors * 3.05
    front_bays = (-10.8, -7.2, -3.6, 0, 3.6, 7.2, 10.8)
    side_bays = (-8.4, -4.2, 0, 4.2, 8.4)
    ops = {side: [] for side in ('front', 'back', 'right', 'left')}
    for floor in range(lower_floors):
        base = floor * 3.05
        for u in front_bays:
            door = floor == 0 and u == 0
            balcony_door = floor > 0 and abs(u) == 7.2
            ops['front'].append(opening(u, base + (1.25 if door or balcony_door else 1.70),
                                        1.8 if door else 2.05, 2.30 if door or balcony_door else 1.72, door))
            ops['back'].append(opening(u, base + 1.70, 1.95, 1.72))
        for side in ('left', 'right'):
            ops[side] += [opening(u, base + 1.70, 2.05, 1.72) for u in side_bays]
    building(p, 28, 24, setback_z, ops)
    upper = {side: [] for side in ('front', 'back', 'right', 'left')}
    center = (-0.8, -0.8)
    for floor in range(2):
        base = setback_z + floor * 3.05
        for side in ('front', 'back'):
            upper[side] += [opening(u, base + 1.70, 2.05, 1.72) for u in (-9, -5.4, -1.8, 1.8, 5.4, 9)]
        for side in ('left', 'right'):
            upper[side] += [opening(u, base + 1.70, 2.05, 1.72) for u in side_bays]
    building(p, 24.6, 22.4, 6.1, upper, center, setback_z, 'Plaster_Ivory')
    # Keep broad-building budget focused on actual windows and modest balconies.
    p['Windows'].pieces = [piece for piece in p['Windows'].pieces if not piece[3]]
    belt(p['Details'], 28, 24, 0.30, 0.60, 'Stone_Base', front_gaps=[(-0.99, 0.99)])
    for floor in range(1, lower_floors + 1):
        belt(p['Details'], 28, 24, min(floor * 3.05, setback_z - 0.04), 0.12, 'Concrete_Warm')
    for floor in range(1, lower_floors):
        for u in (-7.2, 7.2):
            balcony(p['Balconies'], 'front', u, 12, floor * 3.05 + 0.01, 3.25, 1.02, parapet=True, spacing=0.38)
    for u in (-13.9, -1.25, 1.25, 13.9):
        p['Details'].fbox('front', u, 12.03, setback_z / 2, 0.20, 0.14, setback_z, 'Plaster_Sand')
    canopy(p['Details'], 0, 12, 2.70, 3.15, 1.20)
    p['Roof'].box((0, 0, setback_z + 0.02), (28.15, 24.15, 0.16), 'Roof_Membrane')
    rail(p['Balconies'], 'front', 0, 12.015, setback_z + 0.1, 27.90, spacing=0.38)
    rail(p['Balconies'], 'right', 0, 14.015, setback_z + 0.1, 23.85, spacing=0.38)
    belt(p['Details'], 24.6, 22.4, setback_z + 3.05, 0.12, center=center)
    p['Roof'].box((*center, roof_z + 0.01), (24.8, 22.6, 0.20), 'Roof_Membrane')
    for side in ('front', 'back', 'left', 'right'):
        span, d = (24.85, 11.28) if side in ('front', 'back') else (22.45, 12.28)
        p['Roof'].fbox(side, 0, d, roof_z + 0.35, span, 0.17, 0.68, 'Plaster_Ivory', center)
        p['Roof'].fbox(side, 0, d, roof_z + 0.70, span + 0.08, 0.24, 0.07, 'Roof_Slate', center)
    p['Roof'].box((-3, -2.1, roof_z + 0.80), (5, 4.4, 1.55), 'Plaster_Sand')
    p['Roof'].box((-3, -2.1, roof_z + 1.62), (5.25, 4.65, 0.16), 'Roof_Slate')
    return p


def collection(name, parent):
    col = bpy.data.collections.new(name)
    parent.children.link(col)
    return col


def join_meshes(objects, name, col):
    vertices, faces, indices = [], [], []
    for obj in objects:
        offset = len(vertices)
        vertices += [tuple(v.co) for v in obj.data.vertices]
        faces += [tuple(offset + i for i in p.vertices) for p in obj.data.polygons]
        indices += [p.material_index for p in obj.data.polygons]
    mesh = bpy.data.meshes.new(name + '_Mesh')
    mesh.from_pydata(vertices, [], faces)
    mesh.update()
    for mat in MATERIALS:
        mesh.materials.append(mat)
    for face, idx in zip(mesh.polygons, indices):
        face.material_index = idx
    obj = bpy.data.objects.new(name, mesh)
    col.objects.link(obj)
    return obj


def validate_mesh(obj):
    mesh = obj.data
    mesh.calc_loop_triangles()
    degenerate = sum(t.area <= 1e-10 for t in mesh.loop_triangles)
    bm = bmesh.new()
    bm.from_mesh(mesh)
    bm.faces.ensure_lookup_table()
    bad_edges = sum(not edge.is_manifold for edge in bm.edges)
    bad_vertices = sum(not vertex.is_manifold for vertex in bm.verts)
    inconsistent = sum(edge.is_manifold and not edge.is_contiguous for edge in bm.edges)
    zero_faces = sum(face.calc_area() <= 1e-10 for face in bm.faces)
    unseen = set(bm.faces)
    volumes = []
    while unseen:
        seed = unseen.pop()
        stack, shell = [seed], [seed]
        while stack:
            face = stack.pop()
            for edge in face.edges:
                for neighbor in edge.link_faces:
                    if neighbor in unseen:
                        unseen.remove(neighbor)
                        stack.append(neighbor)
                        shell.append(neighbor)
        volume = 0.0
        for face in shell:
            verts = [v.co for v in face.verts]
            a = verts[0]
            for i in range(1, len(verts) - 1):
                volume += a.dot(verts[i].cross(verts[i + 1])) / 6.0
        volumes.append(volume)
    bm.free()
    negative = sum(v < -1e-9 for v in volumes)
    zero_volume = sum(abs(v) <= 1e-9 for v in volumes)
    bounds = [[round(min(v.co[i] for v in mesh.vertices), 5) for i in range(3)],
              [round(max(v.co[i] for v in mesh.vertices), 5) for i in range(3)]]
    report = {'vertices': len(mesh.vertices), 'triangles': len(mesh.loop_triangles),
              'closed_shells': len(volumes), 'nonmanifold_edges': bad_edges, 'nonmanifold_vertices': bad_vertices,
              'inconsistent_winding_edges': inconsistent, 'negative_volume_shells': negative,
              'zero_volume_shells': zero_volume, 'zero_area_faces': zero_faces,
              'degenerate_triangles': degenerate, 'bounds_m': bounds,
              'dimensions_m': [round(bounds[1][i] - bounds[0][i], 5) for i in range(3)],
              'materials': sorted({mesh.materials[p.material_index].name for p in mesh.polygons})}
    errors = {k: report[k] for k in ('nonmanifold_edges', 'nonmanifold_vertices', 'inconsistent_winding_edges',
                                   'negative_volume_shells', 'zero_volume_shells', 'zero_area_faces', 'degenerate_triangles') if report[k]}
    if errors:
        raise RuntimeError(f'{obj.name}: invalid geometry: {errors}')
    return report


def glb_summary(path):
    data = path.read_bytes()
    magic, version, total = struct.unpack_from('<4sII', data)
    assert magic == b'glTF' and version == 2 and total == len(data)
    size, chunk = struct.unpack_from('<II', data, 12)
    assert chunk == 0x4E4F534A
    doc = json.loads(data[20:20 + size].decode('utf-8'))
    triangles = sum(doc['accessors'][p['indices']]['count'] // 3
                    for mesh in doc['meshes'] for p in mesh['primitives'])
    assert all(not m.get('doubleSided', False) for m in doc.get('materials', []))
    assert len(doc['meshes']) == 1
    assert all(p.get('mode', 4) == 4 for mesh in doc['meshes'] for p in mesh['primitives'])
    assert not doc.get('images') and all('uri' not in b for b in doc['buffers'])
    bin_size, bin_kind = struct.unpack_from('<II', data, 20 + size)
    assert bin_kind == 0x004E4942
    binary = data[28 + size:28 + size + bin_size]

    def accessor(index):
        acc = doc['accessors'][index]
        view = doc['bufferViews'][acc['bufferView']]
        components = {'SCALAR': 1, 'VEC3': 3}[acc['type']]
        fmt = '<' + {5123: 'H', 5125: 'I', 5126: 'f'}[acc['componentType']] * components
        stride = view.get('byteStride', struct.calcsize(fmt))
        start = view.get('byteOffset', 0) + acc.get('byteOffset', 0)
        return [struct.unpack_from(fmt, binary, start + i * stride) for i in range(acc['count'])]

    degenerate, wrong_normals = 0, 0
    for primitive in doc['meshes'][0]['primitives']:
        positions = [Vector(v) for v in accessor(primitive['attributes']['POSITION'])]
        normals = [Vector(v) for v in accessor(primitive['attributes']['NORMAL'])]
        indices = [v[0] for v in accessor(primitive['indices'])]
        assert len(indices) % 3 == 0 and all(0 <= i < len(positions) for i in indices)
        assert all(math.isfinite(v) for point in positions for v in point)
        assert all(abs(n.length - 1) < 1e-4 for n in normals)
        for k in range(0, len(indices), 3):
            a, b, c = indices[k:k + 3]
            cross = (positions[b] - positions[a]).cross(positions[c] - positions[a])
            if cross.length <= 2e-10:
                degenerate += 1
            elif any(cross.dot(normals[i]) <= 0 for i in (a, b, c)):
                wrong_normals += 1
    assert degenerate == 0 and wrong_normals == 0, (path.name, degenerate, wrong_normals)
    return {'file': path.name, 'bytes': len(data), 'mesh_count': len(doc['meshes']),
            'triangles': triangles, 'nodes': [n.get('name') for n in doc['nodes']],
            'materials': [m.get('name') for m in doc.get('materials', [])], 'all_double_sided_false': True,
            'decoded_buffer_checks': {'degenerate_triangles': degenerate, 'triangles_opposing_vertex_normals': wrong_normals,
                                      'indices_in_range': True, 'finite_positions': True, 'unit_normals': True}}


def export_glb(obj, path):
    bpy.ops.object.select_all(action='DESELECT')
    obj.hide_set(False)
    obj.select_set(True)
    bpy.context.view_layer.objects.active = obj
    bpy.ops.export_scene.gltf(filepath=str(path), export_format='GLB', use_selection=True,
                              export_yup=True, export_apply=True, export_extras=True,
                              export_cameras=False, export_lights=False, export_animations=False,
                              export_materials='EXPORT', export_texcoords=False, export_normals=True)
    obj.select_set(False)
    return glb_summary(path)


def light(scene, name, kind, location, energy, size):
    data = bpy.data.lights.new(name, kind)
    data.energy = energy
    if kind == 'AREA':
        data.shape = 'DISK'
        data.size = size
    else:
        data.angle = size
    obj = bpy.data.objects.new(name, data)
    scene.collection.objects.link(obj)
    obj.location = location
    obj.rotation_euler = (Vector((0, 0, 8)) - obj.location).to_track_quat('-Z', 'Y').to_euler()
    return obj


def camera_text(scene, camera, text, x, y, size, name, font):
    curve = bpy.data.curves.new(name, 'FONT')
    curve.body, curve.size = text, size
    curve.space_character = 1.10
    if font:
        curve.font = font
        if not font.packed_file:
            font.pack()
    obj = bpy.data.objects.new(name, curve)
    scene.collection.objects.link(obj)
    obj.parent = camera
    obj.location = (x, y, -2.0)
    mat = bpy.data.materials.get('Preview_Ink')
    if not mat:
        mat = bpy.data.materials.new('Preview_Ink')
        mat.use_nodes = True
        nodes = mat.node_tree.nodes
        nodes.clear()
        emission = nodes.new('ShaderNodeEmission')
        emission.inputs['Color'].default_value = (0.035, 0.055, 0.06, 1)
        output = nodes.new('ShaderNodeOutputMaterial')
        mat.node_tree.links.new(emission.outputs[0], output.inputs['Surface'])
    curve.materials.append(mat)
    return obj


def preview_scene(obj, entry, index, samples):
    scene = bpy.data.scenes.new('Preview_' + entry['id'])
    scene.unit_settings.system = 'METRIC'
    scene.unit_settings.scale_length = 1.0
    instance = bpy.data.objects.new('Preview_' + obj.name, obj.data)
    scene.collection.objects.link(instance)
    scene.render.engine = 'CYCLES'
    scene.cycles.device = 'CPU'
    scene.cycles.samples = samples
    scene.cycles.use_denoising = True
    scene.cycles.seed = SEED
    scene.render.resolution_x, scene.render.resolution_y = 840, 1080
    scene.render.resolution_percentage = 100
    scene.render.image_settings.file_format = 'PNG'
    scene.render.image_settings.color_mode = 'RGBA'
    scene.render.image_settings.color_depth = '8'
    scene.view_settings.view_transform = 'AgX'
    scene.world = bpy.data.worlds.new('Preview_World_' + entry['id'])
    scene.world.use_nodes = True
    scene.world.node_tree.nodes['Background'].inputs[0].default_value = (0.75, 0.82, 0.87, 1)
    scene.world.node_tree.nodes['Background'].inputs[1].default_value = 0.65
    ground = bpy.data.materials.get('Preview_Ground')
    if not ground:
        ground = bpy.data.materials.new('Preview_Ground')
        ground.use_nodes = True
        ground.node_tree.nodes['Principled BSDF'].inputs['Base Color'].default_value = (0.77, 0.79, 0.76, 1)
        ground.node_tree.nodes['Principled BSDF'].inputs['Roughness'].default_value = 0.95
    mesh = bpy.data.meshes.new('Preview_Ground_Mesh_' + entry['id'])
    mesh.from_pydata([(-200, -200, -0.012), (200, -200, -0.012), (200, 200, -0.012), (-200, 200, -0.012)], [], [(0, 1, 2, 3)])
    floor = bpy.data.objects.new('Preview_Ground_' + entry['id'], mesh)
    floor.data.materials.append(ground)
    scene.collection.objects.link(floor)
    light(scene, 'Preview_Sun_' + entry['id'], 'SUN', (-30, 45, 60), 2.6, 0.12)
    light(scene, 'Preview_Fill_' + entry['id'], 'AREA', (25, 15, 42), 1400, 30)
    camera_data = bpy.data.cameras.new('Preview_Camera_' + entry['id'])
    camera = bpy.data.objects.new(camera_data.name, camera_data)
    scene.collection.objects.link(camera)
    scene.camera = camera
    low, high = entry['validation']['bounds_m']
    target = Vector(((low[0] + high[0]) / 2, (low[1] + high[1]) / 2, high[2] * 0.50))
    camera.location = target + Vector((1.0, 1.7, 1.0)).normalized() * 120
    camera.rotation_euler = (target - camera.location).to_track_quat('-Z', 'Y').to_euler()
    camera_data.type = 'ORTHO'
    rotation = camera.rotation_euler.to_matrix().transposed()
    corners = [rotation @ (Vector((x, y, z)) - target) for x in (low[0], high[0]) for y in (low[1], high[1]) for z in (0, high[2])]
    aspect = scene.render.resolution_x / scene.render.resolution_y
    span_x = max(c.x for c in corners) - min(c.x for c in corners)
    span_y = max(c.y for c in corners) - min(c.y for c in corners)
    camera_data.ortho_scale = max(span_x / aspect * 1.18, span_y * 1.42)
    scale = camera_data.ortho_scale * aspect
    font_path = Path('C:/Windows/Fonts/segoeui.ttf')
    font = bpy.data.fonts.load(str(font_path), check_existing=True) if font_path.exists() else None
    left, top = -scale * 0.45, scale / aspect * 0.445
    camera_text(scene, camera, f'{index:02d}  {entry["display_name"].upper()}', left, top, scale / 35, 'Preview_Title_' + entry['id'], font)
    footprint = entry['nominal_footprint_m']
    subtitle = f'{footprint[0]} x {footprint[1]} m  /  {entry["floors"]} floors  /  {high[2]:.1f} m'
    camera_text(scene, camera, subtitle, left, top - scale * 0.052, scale / 53, 'Preview_Subtitle_' + entry['id'], font)
    camera_text(scene, camera, f'{entry["validation"]["triangles"]:,} triangles   |   LOD 0', left, -scale / aspect * 0.434, scale / 55, 'Preview_Footer_' + entry['id'], font)
    camera_text(scene, camera, 'HOUSE LIBRARY   /   PROTOTYPE 01', left, -scale / aspect * 0.466, scale / 77, 'Preview_Caption_' + entry['id'], font)
    return scene


def png_rows(path):
    """Read our 8-bit RGBA renders losslessly, without external image packages."""
    data = path.read_bytes()
    assert data[:8] == b'\x89PNG\r\n\x1a\n'
    pos, compressed = 8, bytearray()
    width = height = channels = None
    while pos < len(data):
        length = struct.unpack_from('>I', data, pos)[0]
        kind, body = data[pos + 4:pos + 8], data[pos + 8:pos + 8 + length]
        if kind == b'IHDR':
            width, height, bits, mode, _, _, interlace = struct.unpack('>IIBBBBB', body)
            assert bits == 8 and mode == 6 and interlace == 0
            channels = 4
        elif kind == b'IDAT':
            compressed.extend(body)
        pos += 12 + length
    raw, stride = zlib.decompress(compressed), width * channels
    rows, previous = [], bytearray(stride)
    for y in range(height):
        filter_type = raw[y * (stride + 1)]
        row = bytearray(raw[y * (stride + 1) + 1:(y + 1) * (stride + 1)])
        for i in range(stride):
            a = row[i - channels] if i >= channels else 0
            b, c = previous[i], previous[i - channels] if i >= channels else 0
            if filter_type == 1:
                prediction = a
            elif filter_type == 2:
                prediction = b
            elif filter_type == 3:
                prediction = (a + b) // 2
            elif filter_type == 4:
                p = a + b - c
                distances = (abs(p - a), abs(p - b), abs(p - c))
                prediction = (a, b, c)[distances.index(min(distances))]
            elif filter_type == 0:
                prediction = 0
            else:
                raise ValueError(f'Unknown PNG filter {filter_type}')
            row[i] = (row[i] + prediction) & 255
        rows.append(bytes(row))
        previous = row
    return width, height, rows


def contact_sheet(paths, output):
    panels = [png_rows(path) for path in paths]
    width, height = sum(p[0] for p in panels), panels[0][1]
    assert all(p[1] == height for p in panels)
    payload = b''.join(b'\x00' + b''.join(p[2][y] for p in panels) for y in range(height))

    def chunk(kind, body):
        return struct.pack('>I', len(body)) + kind + body + struct.pack('>I', zlib.crc32(kind + body) & 0xffffffff)

    output.write_bytes(b'\x89PNG\r\n\x1a\n' + chunk(b'IHDR', struct.pack('>IIBBBBB', width, height, 8, 6, 0, 0, 0))
                       + chunk(b'sRGB', b'\x00') + chunk(b'IDAT', zlib.compress(payload, 6)) + chunk(b'IEND', b''))


def main():
    args = arguments()
    assert args.samples > 0
    out = args.output.resolve()
    out.mkdir(parents=True, exist_ok=True)
    bpy.ops.wm.read_factory_settings(use_empty=True)
    source_scene = bpy.context.scene
    source_scene.name = 'Library_Sources'
    source_scene.unit_settings.system = 'METRIC'
    source_scene.unit_settings.scale_length = 1.0
    bpy.context.preferences.filepaths.save_version = 0
    make_materials()
    source_root = collection('HL_Source_Collections', source_scene.collection)
    joined_col = collection('HL_Joined_LOD0', source_scene.collection)
    lod_col = collection('HL_Joined_LOD1', source_scene.collection)
    foundation_col = collection('HL_Optional_Foundations_Below_Zero', source_scene.collection)
    specs = [('house', 'Gabled house', 'House', [12, 9], 3, make_house),
             ('slab', 'Apartment slab', 'Slab', [30, 12], 6, make_slab),
             ('midrise', 'Stepped midrise', 'Midrise', [20, 16], 10, make_midrise),
             ('house-2f', 'Gabled house 2F', 'House', [12, 9], 2, make_house),
             ('house-4f', 'Gabled house 4F', 'House', [12, 9], 4, make_house),
             ('slab-4f', 'Apartment slab 4F', 'Slab', [30, 12], 4, make_slab),
             ('midrise-8f', 'Stepped midrise 8F', 'Midrise', [20, 16], 8, make_midrise),
             ('midrise-wide', 'Wide midrise', 'Midrise', [28, 24], 10, make_wide_midrise)]
    entries, preview_scenes = [], []
    for index, (slug, title, kind, footprint, floors, build) in enumerate(specs, 1):
        stem = f'HL_{kind}_{footprint[0]}x{footprint[1]}_{floors}F'
        print(f'BUILD {slug}', flush=True)
        src = collection('HL_SOURCE_' + slug, source_root)
        RNG.seed(SEED + sum(ord(c) for c in slug))
        parts = build(floors)
        source_objects = [geometry.object('SRC_' + slug + '__' + name, src) for name, geometry in parts.items()]
        source_reports = {obj.name: validate_mesh(obj) for obj in source_objects}
        joined = join_meshes(source_objects, stem + '_LOD0', joined_col)
        low_parts = [geometry.object('TMP_' + slug + '_' + name, lod_col, lod=1) for name, geometry in parts.items()]
        low = join_meshes(low_parts, stem + '_LOD1', lod_col)
        for temp in low_parts:
            mesh = temp.data
            bpy.data.objects.remove(temp, do_unlink=True)
            bpy.data.meshes.remove(mesh)
        report, low_report = validate_mesh(joined), validate_mesh(low)
        assert report['triangles'] <= 35000, (slug, report['triangles'])
        assert low_report['triangles'] < report['triangles']
        assert abs(report['bounds_m'][0][2]) < 1e-6, report['bounds_m']
        for obj, lod in ((joined, 0), (low, 1)):
            obj['archetype'], obj['lod'], obj['floors'] = slug, lod, floors
            obj['nominal_width_m'], obj['nominal_depth_m'] = footprint
            obj['ground_z_m'], obj['authoring_up'], obj['authoring_front'] = 0.0, '+Z', '+Y'
            obj['seed'] = SEED
        glb = export_glb(joined, out / (slug + '.glb'))
        low_glb = export_glb(low, out / (slug + '-lod1.glb'))
        assert glb['triangles'] == report['triangles']
        assert low_glb['triangles'] == low_report['triangles']
        footing_geo = Geometry()
        footing_geo.box((0, 0, -0.275), (footprint[0] + 0.16, footprint[1] + 0.16, 0.50), 'Stone_Base')
        footing = footing_geo.object(stem + '_Foundation', foundation_col)
        footing_report = validate_mesh(footing)
        assert footing_report['bounds_m'][1][2] < 0
        footing_glb = export_glb(footing, out / (slug + '-foundation.glb'))
        entry = {'id': slug, 'display_name': title, 'floors': floors, 'nominal_footprint_m': footprint,
                 'ground_z_m': 0, 'joined_object': joined.name, 'blend_object_names': [joined.name, low.name, footing.name],
                 'source_collection': src.name, 'source_objects': list(source_reports),
                 'validation': report, 'source_validation': source_reports, 'glb': glb,
                 'lod1': {'joined_object': low.name, 'validation': low_report, 'glb': low_glb,
                          'omits': ['roof tile ribs and courses', 'alternate railing pickets', 'window transoms and lintels', 'AC units and downpipes', 'small door handles']},
                 'foundation': {'object': footing.name, 'included_in_main_glb': False, 'validation': footing_report, 'glb': footing_glb}}
        entries.append(entry)
        if index <= 3 or slug == 'midrise-wide':
            preview_scenes.append((entry, preview_scene(joined, entry, index, args.samples)))
        print(f'VALID {slug}: LOD0={report["triangles"]}, LOD1={low_report["triangles"]}, closed_shells={report["closed_shells"]}', flush=True)
    source_root.hide_viewport = True
    source_root.hide_render = True
    lod_col.hide_viewport = True
    lod_col.hide_render = True
    foundation_col.hide_viewport = True
    foundation_col.hide_render = True
    # All linkable source assets retain identity transforms at the local origin.
    for obj in joined_col.objects:
        obj.hide_set(obj.name != entries[0]['joined_object'])
    preview_paths = []
    if not args.skip_render:
        for entry, scene in preview_scenes:
            path = out / ('preview-' + entry['id'] + '.png')
            scene.render.filepath = str(path)
            bpy.context.window.scene = scene
            print('RENDER ' + entry['id'], flush=True)
            bpy.ops.render.render(write_still=True)
            preview_paths.append(path)
        contact_sheet(preview_paths[:3], out / 'preview.png')
    bpy.context.window.scene = preview_scenes[0][1]
    bpy.ops.wm.save_as_mainfile(filepath=str(out / 'house-library.blend'))
    doc = {'schema_version': 1, 'library': 'county-residential-prototype-01', 'seed': SEED,
           'blender_version': bpy.app.version_string, 'blend_file': 'house-library.blend',
           'authoring_coordinates': {'units': 'metres', 'unit_scale': 1, 'up': '+Z', 'front': '+Y',
                                     'origin': 'nominal footprint centre at ground Z=0', 'transforms': 'identity'},
           'gltf_coordinates': {'up': '+Y', 'front': '-Z', 'from_blender_xyz': ['x', 'z', '-y'],
                                'note': 'Standard glTF axis conversion is already applied. Blender import restores Z-up; do not rotate twice.'},
           'recommended_scale_xyz': {'x': [0.85, 1.15], 'y': [0.85, 1.15], 'z': [0.90, 1.12]},
           'assembly': 'Link joined_object from the blend, or import its GLB. Do not instance source collections together with joined objects. Foundation is opt-in below-ground geometry.',
           'validation_method': 'Shell-level Blender bmesh checks: all edges and vertices manifold, contiguous winding, positive signed volume per connected shell; zero-area faces and degenerate loop triangles rejected. Decoded GLB POSITION/NORMAL/index buffers checked for finite positions, unit normals, valid indices, zero degenerate triangles and agreement of face winding with vertex normals. Export also checked for one mesh, embedded data, triangle count, and backface culling. Components intersect intentionally; no boolean-union or inter-component intersection-free certification.',
           'materials': [{'name': m.name, 'base_color_srgb': m['palette_hex_srgb'],
                          'roughness': m.node_tree.nodes['Principled BSDF'].inputs['Roughness'].default_value,
                          'metallic': m.node_tree.nodes['Principled BSDF'].inputs['Metallic'].default_value,
                          'doubleSided': False} for m in MATERIALS],
           'preview': {'file': 'preview.png' if not args.skip_render else None,
                       'individual_files': [p.name for p in preview_paths], 'engine': 'Cycles', 'device': 'CPU',
                       'samples': args.samples, 'contact_sheet_size_px': [2520, 1080],
                       'note': 'Three independently framed actual camera renders; panel scales differ to show architecture.'},
           'limits': ['First reusable asset prototype; not a whole-city model or surveyed architecture.',
                      'Joined meshes contain intersecting but individually closed solids, not a boolean-unioned interior model.',
                      'Opaque tinted PBR glazing avoids alpha sorting and interior dependencies.',
                      'No textures or external downloads. Material names are shared across all archetypes and LODs.',
                      'No Cesium/runtime integration, terrain seating, or georeferencing is performed by this script.'],
           'archetypes': entries}
    (out / 'library.json').write_text(json.dumps(doc, ensure_ascii=False, indent=2) + '\n', encoding='utf-8')
    print('LIBRARY_COMPLETE ' + str(out), flush=True)


if __name__ == '__main__':
    main()