"""Door robot arm: geometry library (models/door-robot-arm). Run inside Blender (imported by dra_build.py).

Metres, Z up. Everything is built from closed shells of quads: surfaces of revolution (lathe), boxes with real bevels and
recessed or raised panels (rbox), outlines pushed out with rounded rims (slab), tubes along paths, knurled rings, and
quad patches with real blind holes (holed_surface). No booleans. Triangles appear only as fans at the poles of turned
parts, at the two ends of a convex cap (column fill), on non-convex caps (triangle fill joined back into quads), and
where a loop cut meets a bevel.

A Link collects the pieces of one moving part (one object). Every piece gets a part id (face layer 'part') and the
list of parts it is allowed to sit in (its seats), so the audit can tell a cap sunk 1 mm into its housing (fine) from a
cable that runs through a wall (not fine).
Ideas and code taken from models/hardsurface-kit/build/kitlib.py and models/pressure-door/build/pd_lib.py.
"""
import math
import bmesh
import bpy
from mathutils import Vector, Matrix, Euler
from mathutils.bvhtree import BVHTree

K = 0.0022                     # metres per ref1 pixel (design.json scale)
TAU = 2.0 * math.pi
W, G, S, R, O, GL, LA = range(7)   # material slots, the same order in every link
MAT_KEYS = ['paint_white', 'gunmetal', 'steel', 'rubber', 'orange', 'glass_dark', 'lamp']
MAT_NAMES = ['M_DoorArm_PaintRed', 'M_DoorArm_Gunmetal', 'M_DoorArm_Steel', 'M_DoorArm_Rubber',
             'M_DoorArm_Marking', 'M_DoorArm_Glass', 'M_DoorArm_Lamp']   # red look 2026-10-19 (were PaintWhite, Orange)
MKEY = {k: i for i, k in enumerate(MAT_KEYS)}
AXES = {'X': (1, 0, 0), '-X': (-1, 0, 0), 'Y': (0, 1, 0), '-Y': (0, -1, 0), 'Z': (0, 0, 1), '-Z': (0, 0, -1)}
DECAL_LIFT = 0.0004            # markings float this far over their surface


def vec(a):
    return Vector(AXES[a]) if isinstance(a, str) else Vector(a).normalized()


def V(*a):
    return Vector(a)


def frame(at=(0, 0, 0), axis='Z', spin=0.0):
    """matrix taking a part built about local +Z into place, local +Z -> axis"""
    a = vec(axis)
    if (a - Vector((0, 0, -1))).length < 1e-9:
        Rm = Matrix.Rotation(math.pi, 4, 'X')
    else:
        Rm = Vector((0, 0, 1)).rotation_difference(a).to_matrix().to_4x4()
    return Matrix.Translation(Vector(at)) @ Rm @ Matrix.Rotation(spin, 4, 'Z')


def basis(at, w, u):
    """matrix: local +Z -> w, local +X -> u (made square to w), origin at"""
    w = Vector(w).normalized()
    u = Vector(u)
    u = (u - w * u.dot(w)).normalized()
    v = w.cross(u)
    M = Matrix((u, v, w)).transposed().to_4x4()
    M.translation = Vector(at)
    return M


# ------------------------------------------------------------------ 2D outlines
def _round(P, Rr, M, closed, n):
    out = []
    N = len(P)
    for i in range(N):
        r = Rr[i]
        if r <= 0 or (not closed and i in (0, N - 1)):
            out.append((P[i].copy(), M[i]))
            continue
        da, db = P[i - 1] - P[i], P[(i + 1) % N] - P[i]
        if da.length < 1e-9 or db.length < 1e-9:
            out.append((P[i].copy(), M[i]))
            continue
        r = min(r, da.length * 0.5, db.length * 0.5)          # tangent length of the round
        ua, ub = da.normalized(), db.normalized()
        if ua.dot(ub) < -0.9999 or ua.dot(ub) > 0.9999:
            out.append((P[i].copy(), M[i]))
            continue
        p0, p2 = P[i] + ua * r, P[i] + ub * r
        half = ua.angle(ub) / 2
        c = P[i] + (ua + ub).normalized() * (r / math.cos(half))
        a, b = p0 - c, p2 - c
        phi = a.angle(b)
        for k in range(n + 1):
            t = k / n
            if phi < 1e-6:
                q = p0.lerp(p2, t)
            else:
                q = c + (a * math.sin((1 - t) * phi) + b * math.sin(t * phi)) / math.sin(phi)
            out.append((q, M[i]))
    ded = []
    for p, m in out:
        if ded and (ded[-1][0] - p).length < 1e-6:
            continue
        ded.append((p, m))
    if closed and len(ded) > 2 and (ded[0][0] - ded[-1][0]).length < 1e-6:
        ded.pop()
    return ded


def fillet(pts, closed=False, n=3):
    """pts: (x, y[, round[, mat]]); 'round' is how far from the corner the rounding starts. -> [(Vector2, mat)]"""
    P = [Vector((p[0], p[1])) for p in pts]
    Rr = [p[2] if len(p) > 2 else 0.0 for p in pts]
    M = [p[3] if len(p) > 3 else 0 for p in pts]
    return _round(P, Rr, M, closed, n)


def fillet3d(pts, rad, n=5):
    P = [Vector(p) for p in pts]
    return [p for p, _ in _round(P, [rad] * len(P), [0] * len(P), False, n)]


def area2(P):
    return sum(P[i - 1].x * P[i].y - P[i].x * P[i - 1].y for i in range(len(P))) / 2


def ccw(P):
    P = [Vector((p[0], p[1])) for p in P]
    return P if area2(P) > 0 else list(reversed(P))


def is_convex(P, eps=1e-9):
    n = len(P)
    for i in range(n):
        a, b, c = P[i - 1], P[i], P[(i + 1) % n]
        if (b - a).x * (c - b).y - (b - a).y * (c - b).x < -eps:
            return False
    return True


def inset_poly(P, d):
    """offset a closed CCW loop inward by d (miter, limited)"""
    n = len(P)
    out = []
    for i in range(n):
        p0, p1, p2 = P[i - 1], P[i], P[(i + 1) % n]
        e1, e2 = (p1 - p0), (p2 - p1)
        if e1.length < 1e-12 or e2.length < 1e-12:
            out.append(p1.copy())
            continue
        e1.normalize()
        e2.normalize()
        n1, n2 = Vector((-e1.y, e1.x)), Vector((-e2.y, e2.x))
        m = n1 + n2
        if m.length < 1e-9:
            m = n1
        m.normalize()
        c = max(0.35, m.dot(n1))
        out.append(p1 + m * (d / c))
    return out


def clip_half(P, nrm, d):
    """Sutherland-Hodgman: keep the part of polygon P with nrm . p <= d"""
    out = []
    n = len(P)
    for i in range(n):
        a, b = P[i - 1], P[i]
        da, db = nrm.dot(a) - d, nrm.dot(b) - d
        if db <= 0:
            if da > 0:
                out.append(a + (b - a) * (da / (da - db)))
            out.append(b.copy())
        elif da <= 0:
            out.append(a + (b - a) * (da / (da - db)))
    ded = []
    for p in out:
        if not ded or (ded[-1] - p).length > 1e-7:
            ded.append(p)
    if len(ded) > 2 and (ded[0] - ded[-1]).length < 1e-7:
        ded.pop()
    return ded


def split_poly(P, p0, p1, gap):
    """cut polygon P along the line p0 -> p1, leaving a gap; returns (left piece, right piece)"""
    p0, p1 = Vector(p0), Vector(p1)
    t = (p1 - p0).normalized()
    nrm = Vector((t.y, -t.x))                    # right of the line
    d = nrm.dot(p0)
    right = clip_half([q * -1 for q in P], nrm, -(d + gap / 2))   # keep nrm.p >= d + gap/2
    right = [q * -1 for q in right]
    left = clip_half(P, nrm, d - gap / 2)
    return left, right


def circle2(c, r, n, a0=0.0, a1=TAU, endpoint=False):
    c = Vector(c)
    m = n + 1 if endpoint else n
    return [c + Vector((math.cos(a0 + (a1 - a0) * k / n), math.sin(a0 + (a1 - a0) * k / n))) * r for k in range(m)]


def hull2(points):
    """convex hull (CCW) of 2D points"""
    pts = sorted(set((round(p[0], 7), round(p[1], 7)) for p in points))
    if len(pts) < 3:
        return [Vector(p) for p in pts]

    def cross(o, a, b):
        return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0])
    lower, upper = [], []
    for p in pts:
        while len(lower) >= 2 and cross(lower[-2], lower[-1], p) <= 0:
            lower.pop()
        lower.append(p)
    for p in reversed(pts):
        while len(upper) >= 2 and cross(upper[-2], upper[-1], p) <= 0:
            upper.pop()
        upper.append(p)
    return [Vector(p) for p in lower[:-1] + upper[:-1]]


def resample(P, step):
    """add points along the edges of a closed loop so no edge is longer than step"""
    out = []
    n = len(P)
    for i in range(n):
        a, b = P[i], P[(i + 1) % n]
        k = max(1, int(math.ceil((b - a).length / step)))
        for j in range(k):
            out.append(a.lerp(b, j / k))
    return out


# ------------------------------------------------------------------ cap filling
def _columns(P, eps=1e-6):
    """convex CCW polygon -> (new loop with matching column points on both sides, faces as index tuples)"""
    n = len(P)
    c = sum(P, Vector((0, 0))) / n
    sxx = sum((p.x - c.x) ** 2 for p in P)
    syy = sum((p.y - c.y) ** 2 for p in P)
    sxy = sum((p.x - c.x) * (p.y - c.y) for p in P)
    ang = 0.5 * math.atan2(2 * sxy, sxx - syy)
    d = Vector((math.cos(ang), math.sin(ang)))
    t = [(p - c).dot(d) for p in P]
    i0 = min(range(n), key=lambda i: (t[i], i))
    i1 = max(range(n), key=lambda i: (t[i], -i))
    chA = [i0]
    i = i0
    while i != i1:
        i = (i + 1) % n
        chA.append(i)
    chB = [i0]
    i = i0
    while i != i1:
        i = (i - 1) % n
        chB.append(i)
    T = sorted(t)
    Tu = []
    for v in T:
        if not Tu or v - Tu[-1] > eps:
            Tu.append(v)

    def expand(ch):
        res = []                                   # list of (point, t, is_original_index)
        for k in range(len(ch) - 1):
            a, b = ch[k], ch[k + 1]
            ta, tb = t[a], t[b]
            res.append((P[a], ta))
            lo, hi = (ta, tb) if ta < tb else (tb, ta)
            for v in Tu:
                if lo + eps < v < hi - eps:
                    f = (v - ta) / (tb - ta)
                    res.append((P[a].lerp(P[b], f), v))
        res.append((P[ch[-1]], t[ch[-1]]))
        return res
    A, B = expand(chA), expand(chB)

    def by_t(lst):
        cols = []
        for p, v in lst:
            if cols and abs(cols[-1][0] - v) <= eps:
                cols[-1][1].append(p)
            else:
                cols.append([v, [p]])
        return cols
    cA, cB = by_t(A), by_t(B)
    if len(cA) != len(cB):
        return None
    # the new loop: A forward, then B backward (without its shared ends)
    loop = [p for p, _ in A] + [p for p, _ in reversed(B[1:-1])]
    idx = {}
    for k, p in enumerate(loop):
        idx.setdefault((round(p.x, 9), round(p.y, 9)), k)

    def key(p):
        return idx[(round(p.x, 9), round(p.y, 9))]
    faces = []
    for k in range(len(cA) - 1):
        a0, a1 = cA[k][1][-1], cA[k + 1][1][0]
        b0, b1 = cB[k][1][-1], cB[k + 1][1][0]
        vs = []
        for p in (a0, a1, b1, b0):
            q = key(p)
            if q not in vs:
                vs.append(q)
        if len(vs) >= 3:
            faces.append(tuple(vs))
    return loop, faces


def fill_cap(bm, ring, out, mat=0, holes=(), faces_hint=None):
    """fill a planar loop of BMVerts (with optional hole loops) facing 'out'. Uses the column faces if given."""
    made = []
    if faces_hint is not None and not holes:
        for f in faces_hint:
            ff = face(bm, [ring[i] for i in f], mat, out=out)
            if ff:
                made.append(ff)
        return made
    edges = []
    for loop in [ring] + list(holes):
        for i in range(len(loop)):
            a, b = loop[i], loop[(i + 1) % len(loop)]
            e = bm.edges.get((a, b)) or bm.edges.new((a, b))
            edges.append(e)
    res = bmesh.ops.triangle_fill(bm, use_beauty=True, use_dissolve=False, edges=edges, normal=Vector(out))
    tri = [g for g in res['geom'] if isinstance(g, bmesh.types.BMFace)]
    if holes:
        nrm = Vector(out).normalized()
        ref = Vector((1, 0, 0)) if abs(nrm.x) < 0.9 else Vector((0, 1, 0))
        ua = (ref - nrm * ref.dot(nrm)).normalized()
        va = nrm.cross(ua)
        polys = [[Vector((v.co.dot(ua), v.co.dot(va))) for v in loop] for loop in holes]

        def inside(p, P):
            c = False
            j = len(P) - 1
            for i in range(len(P)):
                if (P[i].y > p.y) != (P[j].y > p.y) and p.x < (P[j].x - P[i].x) * (p.y - P[i].y) / (P[j].y - P[i].y) + P[i].x:
                    c = not c
                j = i
            return c
        kill = []
        for f in tri:
            cc = f.calc_center_median()
            q = Vector((cc.dot(ua), cc.dot(va)))
            if any(inside(q, P) for P in polys):
                kill.append(f)
        if kill:
            bmesh.ops.delete(bm, geom=kill, context='FACES_ONLY')
            tri = [f for f in tri if f.is_valid]
    for f in tri:
        f.material_index = mat
        f.normal_update()
        if f.normal.dot(Vector(out)) < 0:
            f.normal_flip()
    if tri:
        jr = bmesh.ops.join_triangles(bm, faces=tri, angle_face_threshold=math.radians(5),
                                      angle_shape_threshold=math.radians(70))
        made = [f for f in bm.faces if f.is_valid and f in set(tri)] if False else []
    return made


def face(bm, verts, mat=0, out=None, ref=None):
    uq = []
    for v in verts:
        if v not in uq:
            uq.append(v)
    if len(uq) < 3:
        return None
    try:
        f = bm.faces.new(uq)
    except ValueError:
        return None
    f.material_index = mat
    f.normal_update()
    if out is None and ref is not None:
        out = f.calc_center_median() - Vector(ref)
    if out is not None and f.normal.dot(Vector(out)) < 0:
        f.normal_flip()
    return f


# ------------------------------------------------------------------ primitives (each returns a new closed bmesh)
def lathe(prof, segs=32, closed=False, n=2, phase=0.0):
    """surface of revolution about +Z; prof = (r, z[, round[, mat]]); points with r = 0 sit on the axis (fans)"""
    pts = fillet(prof, closed, n)
    bm = bmesh.new()
    angs = [phase + TAU * j / segs for j in range(segs)]
    rings = []
    for p, _ in pts:
        if p.x < 1e-6:
            rings.append([bm.verts.new((0.0, 0.0, p.y))])
        else:
            rings.append([bm.verts.new((p.x * math.cos(a), p.x * math.sin(a), p.y)) for a in angs])
    L = len(rings)
    for i in range(L if closed else L - 1):
        A, B, m = rings[i], rings[(i + 1) % L], pts[i][1]
        if len(A) == 1 and len(B) == 1:
            continue
        for j in range(segs):
            j2 = (j + 1) % segs
            if len(A) == 1:
                vs = (A[0], B[j], B[j2])
            elif len(B) == 1:
                vs = (A[j], A[j2], B[0])
            else:
                vs = (A[j], A[j2], B[j2], B[j])
            try:
                bm.faces.new(vs).material_index = m
            except ValueError:
                pass
    bmesh.ops.remove_doubles(bm, verts=bm.verts[:], dist=1e-7)
    bmesh.ops.recalc_face_normals(bm, faces=bm.faces[:])
    return bm


def stack(rings, segs=32, f=0.0012, base=0.0, n=1):
    """a stack of rings sitting on a face (local z = 0) and growing along +Z: [(diameter, height, mat)], each on top of
    the last. base: how far the first ring sinks into its seat."""
    m0 = rings[0][2]
    prof = [(0.0, -base, 0, MKEY[m0] if isinstance(m0, str) else m0)]
    z = -base
    for k, (dia, h, m) in enumerate(rings):
        r = dia / 2
        mi = MKEY[m] if isinstance(m, str) else m
        prof.append((r, z, f if k else 0.0, mi))
        z += h + (base if k == 0 else 0.0)
        prof.append((r, z, min(f, h * 0.3), mi))
    prof.append((0.0, z, 0, prof[-1][3]))
    return lathe(prof, segs, n=n)


def rbox(size, r=0.003, ops=(), rbev=0.0015, segs=2, rv=None, vax=2, segs_v=3, rim_segs=1):
    """box centred on the origin with real bevels. rv: a larger radius on the 4 edges parallel to axis vax.
    ops: (side, u0, u1, v0, v1, depth) on side '+z' (or '-x', ...): the faces inside the rectangle (u, v = the side's
    other two axes in x, y, z order) are pushed by depth: < 0 a recess, > 0 a raised panel."""
    OTHER = {0: (1, 2), 1: (0, 2), 2: (0, 1)}
    half = [s / 2 for s in size]
    bm = bmesh.new()
    bmesh.ops.create_cube(bm, size=1.0)
    for v in bm.verts:
        v.co = Vector((v.co.x * size[0], v.co.y * size[1], v.co.z * size[2]))
    for (side, u0, u1, v0, v1, d) in ops:
        ai = 'xyz'.index(side[-1].lower())
        ui, vi = OTHER[ai]
        for c, ax in ((u0, ui), (u1, ui), (v0, vi), (v1, vi)):
            if abs(c) >= half[ax] - 1e-4:
                continue
            no, co = [0, 0, 0], [0, 0, 0]
            no[ax], co[ax] = 1.0, c
            bmesh.ops.bisect_plane(bm, geom=bm.verts[:] + bm.edges[:] + bm.faces[:], dist=1e-6, plane_co=co, plane_no=no)
    for (side, u0, u1, v0, v1, d) in ops:
        ai = 'xyz'.index(side[-1].lower())
        sgn = -1.0 if side[0] == '-' else 1.0
        ui, vi = OTHER[ai]
        dv = Vector((0, 0, 0))
        dv[ai] = sgn
        bm.normal_update()
        cand = []
        for f in bm.faces:
            if f.normal.dot(dv) < 0.999:
                continue
            c = f.calc_center_median()
            if u0 + 1e-5 < c[ui] < u1 - 1e-5 and v0 + 1e-5 < c[vi] < v1 - 1e-5:
                cand.append((c[ai] * sgn, f))
        if not cand:
            continue
        top = max(h for h, _ in cand)
        faces = [f for h, f in cand if h > top - 1e-5]
        ret = bmesh.ops.extrude_face_region(bm, geom=faces)
        nv = [g for g in ret["geom"] if isinstance(g, bmesh.types.BMVert)]
        bmesh.ops.translate(bm, vec=dv * d, verts=nv)
        bmesh.ops.delete(bm, geom=faces, context='FACES')
    bmesh.ops.recalc_face_normals(bm, faces=bm.faces[:])

    def planes(e):
        m = (e.verts[0].co + e.verts[1].co) / 2
        return [a for a in range(3) if abs(abs(m[a]) - half[a]) < 1e-5]

    def sharp(ang):
        return [e for e in bm.edges if len(e.link_faces) == 2 and e.calc_face_angle(0.0) > ang]
    if ops and rbev > 0:
        rims = [e for e in sharp(1.0) if len(planes(e)) < 2]
        if rims:
            bmesh.ops.bevel(bm, geom=rims, offset=rbev, offset_type='OFFSET', segments=rim_segs, profile=0.5,
                            affect='EDGES', clamp_overlap=True)
    if rv and rv > 0:
        vert = []
        for e in sharp(1.0):
            pl = planes(e)
            dvec = e.verts[1].co - e.verts[0].co
            if len(pl) >= 2 and vax not in pl and abs(dvec.normalized()[vax]) > 0.99:
                vert.append(e)
        if vert:
            bmesh.ops.bevel(bm, geom=vert, offset=rv, offset_type='OFFSET', segments=segs_v, profile=0.5,
                            affect='EDGES', clamp_overlap=True)
    if r > 0:
        outer = []
        for e in sharp(1.0):
            pl = planes(e)
            if len(pl) >= 2 or (rv and vax in pl):
                outer.append(e)
        if outer:
            bmesh.ops.bevel(bm, geom=outer, offset=r, offset_type='OFFSET', segments=segs, profile=0.5,
                            affect='EDGES', clamp_overlap=True)
    ngons = [f for f in bm.faces if len(f.verts) > 4]
    if ngons:
        bm.normal_update()
        bmesh.ops.triangulate(bm, faces=ngons, quad_method='BEAUTY', ngon_method='BEAUTY')
        bmesh.ops.join_triangles(bm, faces=[f for f in bm.faces if len(f.verts) == 3],
                                 angle_face_threshold=math.radians(2), angle_shape_threshold=math.radians(60))
    bmesh.ops.recalc_face_normals(bm, faces=bm.faces[:])
    return bm


def slab(outline, z0, z1, bev=0.003, n=3, steps=2, holes=(), bev0=None, mat=0, side_mat=None, step=None):
    """a closed outline (x, y[, round]) pushed from z0 to z1; the top rim rounded with bev, the bottom rim with bev0
    (default bev). holes: closed loops (x, y[, round]) cut through the slab. Caps: convex without holes -> quads in
    columns; otherwise triangle fill joined into quads."""
    bev0 = bev if bev0 is None else bev0
    pts = [p for p, _ in fillet(outline, closed=True, n=n)]
    pts = ccw(pts)
    if step:
        pts = resample(pts, step)
    hl = []
    for h in holes:
        hp = ccw([p for p, _ in fillet(h, closed=True, n=n)])
        hl.append(list(reversed(hp)))          # holes run clockwise
    cols = None
    if not hl and is_convex(pts, 1e-7):
        cc = _columns(pts)
        if cc:
            pts, cols = cc
    sm = mat if side_mat is None else side_mat
    bm = bmesh.new()
    th = z1 - z0
    b0 = min(bev0, th * 0.45)
    b1 = min(bev, th * 0.45)

    def levels(b, zbase, up):
        lv = []
        if b < 1e-6:
            return [(0.0, zbase)]
        for k in range(steps + 1):
            a = (math.pi / 2) * k / steps
            if up:      # from the bottom face up the rim
                lv.append((b * (1 - math.sin(a)), zbase + b * (1 - math.cos(a))))
            else:
                lv.append((b * (1 - math.cos(a)), zbase - b + b * math.sin(a)))
        return lv
    L = levels(b0, z0, True) + levels(b1, z1, False)

    def rings_for(loop, sign):
        rr = []
        for ins, z in L:
            Q = inset_poly(loop, ins * sign) if ins > 1e-9 else loop
            rr.append([bm.verts.new((q.x, q.y, z)) for q in Q])
        return rr
    outer = rings_for(pts, 1.0)

    def wall(rings, inward):
        for i in range(len(rings) - 1):
            A, B = rings[i], rings[i + 1]
            for k in range(len(A)):
                k2 = (k + 1) % len(A)
                e = Vector(A[k2].co - A[k].co)
                o = Vector((e.y, -e.x, 0.0))
                if inward:
                    o = -o
                face(bm, [A[k], A[k2], B[k2], B[k]], sm, out=o)
    wall(outer, False)
    hring = []
    for h in hl:
        rr = rings_for(h, 1.0)            # h is clockwise: inset of a CW loop grows it into the material
        wall(rr, False)
        hring.append(rr)
    fill_cap(bm, outer[-1], (0, 0, 1), mat, holes=[r[-1] for r in hring], faces_hint=cols)
    fill_cap(bm, outer[0], (0, 0, -1), mat, holes=[r[0] for r in hring], faces_hint=cols)
    bmesh.ops.remove_doubles(bm, verts=bm.verts[:], dist=1e-7)
    return bm


def tube(points, r, segs=12, rad=0.02, caps=True, n=5, mat=0):
    """a round bar or cable swept along a polyline whose corners are bent with radius rad"""
    P = fillet3d(points, rad, n) if rad > 0 and len(points) > 2 else [Vector(p) for p in points]
    T = []
    for i in range(len(P)):
        a, b = P[max(i - 1, 0)], P[min(i + 1, len(P) - 1)]
        T.append((b - a).normalized())
    ref = Vector((0, 0, 1)) if abs(T[0].z) < 0.9 else Vector((1, 0, 0))
    nrm = (ref - T[0] * ref.dot(T[0])).normalized()
    bm = bmesh.new()
    rings = []
    for i, p in enumerate(P):
        nrm = (nrm - T[i] * nrm.dot(T[i])).normalized()
        bi = T[i].cross(nrm)
        rings.append([bm.verts.new(p + r * (math.cos(TAU * j / segs) * nrm + math.sin(TAU * j / segs) * bi)) for j in range(segs)])
    for A, B in zip(rings, rings[1:]):
        for j in range(segs):
            j2 = (j + 1) % segs
            f = bm.faces.new((A[j], A[j2], B[j2], B[j]))
            f.material_index = mat
    if caps:
        for ring, p in ((rings[0], P[0]), (rings[-1], P[-1])):
            c = bm.verts.new(p)
            for j in range(segs):
                bm.faces.new((ring[j], ring[(j + 1) % segs], c)).material_index = mat
    bmesh.ops.recalc_face_normals(bm, faces=bm.faces[:])
    return bm


def path_frames(P):
    """parallel-transport frames (T, N, B) along a polyline"""
    T = []
    for i in range(len(P)):
        a, b = P[max(i - 1, 0)], P[min(i + 1, len(P) - 1)]
        T.append((b - a).normalized())
    ref = Vector((0, 0, 1)) if abs(T[0].z) < 0.9 else Vector((1, 0, 0))
    nrm = (ref - T[0] * ref.dot(T[0])).normalized()
    out = []
    for i in range(len(P)):
        nrm = (nrm - T[i] * nrm.dot(T[i])).normalized()
        out.append((T[i], nrm.copy(), T[i].cross(nrm)))
    return out


def bundle(points, r_cable, offsets, rad=0.05, segs=10, mat=R, n=4):
    """several cables side by side along one path (offsets in the path's N, B plane)"""
    P = fillet3d(points, rad, n)
    F = path_frames(P)
    bm = bmesh.new()
    for (ox, oy) in offsets:
        path = [p + Fn * ox + Fb * oy for p, (Ft, Fn, Fb) in zip(P, F)]
        merge(bm, tube(path, r_cable, segs, 0.0, mat=mat))
    return bm


def knurl(r_root, r_tip, length, ribs, ch=0.003, mat=S):
    """a knurled (ribbed) ring along +Z from 0 to length, closed with flat ends (fans)"""
    bm = bmesh.new()
    angs = []
    for k in range(ribs):
        a = TAU * k / ribs
        w = TAU / ribs
        angs += [(a, r_root), (a + w * 0.18, r_tip), (a + w * 0.50, r_tip), (a + w * 0.68, r_root)]
    zs = [(0.0, 0.0), (ch, 1.0), (length - ch, 1.0), (length, 0.0)]
    rings = []
    for z, f in zs:
        ring = []
        for a, rr in angs:
            r = r_root - ch + (rr - r_root + ch) * f if f < 1 else rr
            r = min(r, rr)
            ring.append(bm.verts.new((r * math.cos(a), r * math.sin(a), z)))
        rings.append(ring)
    n = len(angs)
    for A, B in zip(rings, rings[1:]):
        for j in range(n):
            j2 = (j + 1) % n
            bm.faces.new((A[j], A[j2], B[j2], B[j])).material_index = mat
    for ring, z in ((rings[0], 0.0), (rings[-1], length)):
        c = bm.verts.new((0, 0, z))
        for j in range(n):
            bm.faces.new((ring[j], ring[(j + 1) % n], c)).material_index = mat
    bmesh.ops.recalc_face_normals(bm, faces=bm.faces[:])
    return bm


def star_slab(outline, z0, z1, hub_r=0.0, bev=0.0015, mat=0):
    """a star-shaped outline (points in angular order round the origin) pushed from z0 to z1; caps filled with radial
    quads to a hub ring of hub_r, the hub closed with a fan. Rims chamfered by bev."""
    P = ccw(outline)
    bm = bmesh.new()
    th = z1 - z0
    b = min(bev, th * 0.4)
    L = [(b, z0), (0.0, z0 + b), (0.0, z1 - b), (b, z1)] if b > 1e-6 else [(0.0, z0), (0.0, z1)]
    rings = []
    for ins, z in L:
        Q = inset_poly(P, ins) if ins > 0 else P
        rings.append([bm.verts.new((q.x, q.y, z)) for q in Q])
    n = len(P)
    for A, B in zip(rings, rings[1:]):
        for k in range(n):
            k2 = (k + 1) % n
            e = P[k2] - P[k]
            face(bm, [A[k], A[k2], B[k2], B[k]], mat, out=(e.y, -e.x, 0))
    for ring, z, sg in ((rings[0], z0, -1), (rings[-1], z1, 1)):
        hub = []
        for v in ring:
            d = Vector((v.co.x, v.co.y))
            dd = d.normalized() * hub_r if hub_r > 0 else Vector((0, 0))
            hub.append(bm.verts.new((dd.x, dd.y, z)))
        if hub_r > 0:
            for k in range(n):
                k2 = (k + 1) % n
                face(bm, [ring[k], ring[k2], hub[k2], hub[k]], mat, out=(0, 0, sg))
            c = bm.verts.new((0, 0, z))
            for k in range(n):
                face(bm, [hub[k], hub[(k + 1) % n], c], mat, out=(0, 0, sg))
        else:
            for v in hub:
                bm.verts.remove(v)
            c = bm.verts.new((0, 0, z))
            for k in range(n):
                face(bm, [ring[k], ring[(k + 1) % n], c], mat, out=(0, 0, sg))
    bmesh.ops.remove_doubles(bm, verts=bm.verts[:], dist=1e-7)
    return bm


def holed_surface(bm, mapf, us, vs, holes, hole_r, depth, nrmf, periodic_u=False, mat=0, hole_mat=None, nh=8,
                  bottom_mat=None, ret_lattice=False):
    """a quad patch over the cells (us[i]..us[i+1]) x (vs[j]..vs[j+1]); the cells in 'holes' get a round blind hole of
    radius hole_r (metres) and depth. mapf(u, v) -> point; nrmf(u, v) -> outward normal.
    Every cell edge is split in two so cells with and without holes share their edges (no T-junctions).
    Returns the lattice rings along u at the first and last v (for joining) as lists of BMVerts."""
    hm = mat if hole_mat is None else hole_mat
    bmat = hm if bottom_mat is None else bottom_mat
    nu = len(us) - 1
    nv = len(vs) - 1
    U = []
    for i in range(nu):
        U += [us[i], (us[i] + us[i + 1]) / 2]
    if not periodic_u:
        U.append(us[-1])
    Vv = []
    for j in range(nv):
        Vv += [vs[j], (vs[j] + vs[j + 1]) / 2]
    Vv.append(vs[-1])
    lat = [[bm.verts.new(mapf(u, v)) for u in U] for v in Vv]
    NU = len(U)

    def L(iu, jv):
        return lat[jv][iu % NU] if periodic_u else lat[jv][iu]
    for i in range(nu):
        for j in range(nv):
            i0, j0 = 2 * i, 2 * j
            uc, vc = (us[i] + us[i + 1]) / 2, (vs[j] + vs[j + 1]) / 2
            n_out = nrmf(uc, vc)
            if (i, j) not in holes:
                for (a, b) in ((0, 0), (1, 0), (0, 1), (1, 1)):
                    face(bm, [L(i0 + a, j0 + b), L(i0 + a + 1, j0 + b), L(i0 + a + 1, j0 + b + 1), L(i0 + a, j0 + b + 1)],
                         mat, out=n_out)
                continue
            # boundary of the cell, counter-clockwise from the right middle
            bnd = [L(i0 + 2, j0 + 1), L(i0 + 2, j0 + 2), L(i0 + 1, j0 + 2), L(i0, j0 + 2),
                   L(i0, j0 + 1), L(i0, j0), L(i0 + 1, j0), L(i0 + 2, j0)]
            pc = mapf(uc, vc)
            du = (mapf(uc + (us[i + 1] - us[i]) * 0.01, vc) - pc)
            dv = (mapf(uc, vc + (vs[j + 1] - vs[j]) * 0.01) - pc)
            eu, ev = du.normalized(), dv.normalized()
            ev = (ev - eu * ev.dot(eu)).normalized()
            ring0, ring1 = [], []
            for k in range(8):
                a = TAU * k / 8
                off = eu * (math.cos(a) * hole_r) + ev * (math.sin(a) * hole_r)
                ring0.append(bm.verts.new(pc + off))
                ring1.append(bm.verts.new(pc + off - n_out * depth))
            for k in range(8):
                k2 = (k + 1) % 8
                face(bm, [bnd[k], bnd[k2], ring0[k2], ring0[k]], mat, out=n_out)
                face(bm, [ring0[k], ring0[k2], ring1[k2], ring1[k]], hm, ref=None,
                     out=(pc - (ring0[k].co + ring0[k2].co) / 2))
            cb = bm.verts.new(pc - n_out * depth)
            for k in range(8):
                face(bm, [ring1[k], ring1[(k + 1) % 8], cb], bmat, out=n_out)
    if ret_lattice:
        return lat
    first = [L(iu, 0) for iu in range(NU)]
    last = [L(iu, len(Vv) - 1) for iu in range(NU)]
    return first, last


def bridge(bm, A, B, mat=0, ref=None, out=None, closed=True):
    n = len(A)
    for k in range(n if closed else n - 1):
        k2 = (k + 1) % n
        face(bm, [A[k], A[k2], B[k2], B[k]], mat, out=out, ref=ref)


def merge(dst, src, mat=None):
    src.verts.index_update()
    vm = [dst.verts.new(v.co) for v in src.verts]
    lay_s = src.faces.layers.int.get('part')
    lay_d = dst.faces.layers.int.get('part')
    for f in src.faces:
        try:
            nf = dst.faces.new([vm[v.index] for v in f.verts])
        except ValueError:
            continue
        nf.material_index = f.material_index if mat is None else mat
        nf.smooth = f.smooth
        if lay_s and lay_d:
            nf[lay_d] = f[lay_s]


def xform(bm, M):
    bmesh.ops.transform(bm, matrix=M, verts=bm.verts[:])
    if M.determinant() < 0:
        bmesh.ops.reverse_faces(bm, faces=bm.faces[:])
    return bm


def mirror_y(bm):
    return xform(bm, Matrix.Scale(-1, 4, Vector((0, 1, 0))))


# ------------------------------------------------------------------ markings (thin floating geometry)
def decal_strip(pts2, width, closed=True):
    """an outline drawn as a strip of quads of the given width, in the local XY plane"""
    P = [Vector(p) for p in pts2]
    if closed:
        P = ccw(P)
        Pin = inset_poly(P, width / 2)
        Pout = inset_poly(P, -width / 2)
    else:
        Pin, Pout = [], []
        for i, p in enumerate(P):
            a, b = P[max(i - 1, 0)], P[min(i + 1, len(P) - 1)]
            t = (b - a).normalized()
            nn = Vector((-t.y, t.x))
            Pin.append(p + nn * width / 2)
            Pout.append(p - nn * width / 2)
    bm = bmesh.new()
    A = [bm.verts.new((p.x, p.y, 0)) for p in Pout]
    B = [bm.verts.new((p.x, p.y, 0)) for p in Pin]
    n = len(P)
    for k in range(n if closed else n - 1):
        k2 = (k + 1) % n
        face(bm, [A[k], A[k2], B[k2], B[k]], O, out=(0, 0, 1))
    return bm


def decal_poly(pts2):
    P = ccw([Vector(p) for p in pts2])
    bm = bmesh.new()
    vs = [bm.verts.new((p.x, p.y, 0)) for p in P]
    if len(vs) <= 4:
        face(bm, vs, O, out=(0, 0, 1))
    else:
        c = sum(P, Vector((0, 0))) / len(P)
        cv = bm.verts.new((c.x, c.y, 0))
        for k in range(len(vs)):
            face(bm, [vs[k], vs[(k + 1) % len(vs)], cv], O, out=(0, 0, 1))
    return bm


def decal_dot(c, r, n=10):
    return decal_poly(circle2(c, r, n))


# ------------------------------------------------------------------ the link builder
class Link:
    """the pieces of one object; every piece keeps its part id (for the audit) and its seats"""

    def __init__(self, name):
        self.name = name
        self.bm = bmesh.new()
        self.lay = self.bm.faces.layers.int.new('part')
        self.parts = []                 # [{'id', 'seats', 'decal'}]
        self.index = {}
        self._bvh = None

    def pid(self, pid, seats=(), decal=False):
        if pid not in self.index:
            self.index[pid] = len(self.parts)
            self.parts.append({'id': pid, 'seats': set(seats), 'decal': decal})
        else:
            self.parts[self.index[pid]]['seats'] |= set(seats)
        return self.index[pid]

    def put(self, piece, M=None, mat=None, pid='?', seats=(), decal=False):
        if M is not None:
            xform(piece, M)
        k = self.pid(pid, seats, decal)
        lay = piece.faces.layers.int.get('part') or piece.faces.layers.int.new('part')
        for f in piece.faces:
            f[lay] = k
        merge(self.bm, piece, mat)
        piece.free()
        self._bvh = None
        return self

    def bvh(self):
        if self._bvh is None:
            self._bvh = BVHTree.FromBMesh(self.bm)
        return self._bvh

    def snap(self, origin, direction, maxd=5.0):
        """first hit from origin along direction on what is built so far -> (point, normal) or (None, None)"""
        hit = self.bvh().ray_cast(Vector(origin), Vector(direction).normalized(), maxd)
        if hit[0] is None:
            return None, None
        n = Vector(hit[1])
        if n.length < 0.5:                      # a degenerate face: take the ray's own direction
            n = -Vector(direction).normalized()
        if n.dot(Vector(direction)) > 0:
            n = -n
        return hit[0], n.normalized()

    def on_surface(self, origin, direction, up=(0, 0, 1)):
        """a frame sitting on the surface hit from origin along direction: local Z = surface normal"""
        p, n = self.snap(origin, direction)
        if p is None:
            return None
        u = Vector(up)
        if abs(u.normalized().dot(n)) > 0.95:
            u = Vector((1, 0, 0)) if abs(n.x) < 0.9 else Vector((0, 1, 0))
        # local X along 'up' projected (so text-free markings keep their orientation)
        M = basis(p, n, u.cross(n) if False else (u - n * u.dot(n)))
        return M

    def decal(self, piece, M, pid):
        xform(piece, Matrix.Translation((0, 0, DECAL_LIFT)))
        return self.put(piece, M, None, pid, decal=True)


def finish(link, name, origin=(0, 0, 0), mats=None, sharp_deg=32.0, collection=None):
    """turn a Link into a mesh datablock 'name' with its origin at 'origin' (link coordinates); returns the mesh"""
    bm = link.bm
    bmesh.ops.remove_doubles(bm, verts=bm.verts[:], dist=1e-7)
    bmesh.ops.dissolve_degenerate(bm, dist=2e-6, edges=bm.edges[:])
    bad = [f for f in bm.faces if f.calc_area() < 1e-10]
    if bad:
        bmesh.ops.delete(bm, geom=bad, context='FACES')
    ngons = [f for f in bm.faces if len(f.verts) > 4]
    if ngons:
        bm.normal_update()
        bmesh.ops.triangulate(bm, faces=ngons, quad_method='BEAUTY', ngon_method='BEAUTY')
    me = bpy.data.meshes.get(name)
    if me is None:
        me = bpy.data.meshes.new(name)
    else:
        me.clear_geometry()
        me.materials.clear()
    bm.to_mesh(me)
    me.transform(Matrix.Translation(-Vector(origin)))
    used = sorted({p.material_index for p in me.polygons})
    remap = {old: new for new, old in enumerate(used)}
    for p in me.polygons:
        p.material_index = remap[p.material_index]
    for old in used:
        me.materials.append(bpy.data.materials[MAT_NAMES[old]] if mats is None else mats[old])
    me.shade_smooth()
    me.set_sharp_from_angle(angle=math.radians(sharp_deg))
    me["parts"] = [p['id'] for p in link.parts]
    me["seats"] = [";".join(sorted(p['seats'])) for p in link.parts]
    me["decals"] = [1 if p['decal'] else 0 for p in link.parts]
    me["origin_in_link"] = list(origin)
    return me


def uv_tile(me, scale=1.0):
    """UV_Tile: box projection in metres (by each face's main normal axis)"""
    if "UV_Tile" in me.uv_layers:
        me.uv_layers.remove(me.uv_layers["UV_Tile"])
    lay = me.uv_layers.new(name="UV_Tile")
    co = [v.co for v in me.vertices]
    for p in me.polygons:
        n = p.normal
        ax = max(range(3), key=lambda a: abs(n[a]))
        a1, a2 = [(1, 2), (0, 2), (0, 1)][ax]
        sg = 1.0 if n[ax] >= 0 else -1.0
        for li in p.loop_indices:
            c = co[me.loops[li].vertex_index]
            u = c[a1] * (sg if ax != 1 else -sg)
            lay.data[li].uv = (u / scale, c[a2] / scale)
    return lay


# ------------------------------------------------------------------ outlines from signed distances (marching squares)
def sd_circle(c, r):
    c = Vector(c)
    return lambda p: (p - c).length - r


def sd_convex(pts):
    """signed distance (inside < 0) to a convex polygon (any order): the max over its edges"""
    P = hull2(pts)
    E = []
    for i in range(len(P)):
        a, b = P[i], P[(i + 1) % len(P)]
        t = (b - a).normalized()
        n = Vector((t.y, -t.x))
        E.append((n, n.dot(a)))
    return lambda p: max(n.dot(p) - d for n, d in E)


def sd_hull(circles=(), points=(), n=48):
    pts = [Vector(p) for p in points]
    for c, r in circles:
        pts += circle2(c, r, n)
    return sd_convex(pts)


def sd_union(*fs):
    return lambda p: min(f(p) for f in fs)


def sd_inter(*fs):
    return lambda p: max(f(p) for f in fs)


def sd_minus(a, b):
    return lambda p: max(a(p), -b(p))


def sd_half(p0, nrm):
    """the half plane on the side of -nrm (inside where nrm . (p - p0) < 0)"""
    p0, nrm = Vector(p0), Vector(nrm).normalized()
    return lambda p: nrm.dot(p - p0)


def _dp(pts, tol):
    if len(pts) < 3:
        return pts
    a, b = pts[0], pts[-1]
    ab = b - a
    L = ab.length
    best, bi = -1.0, 0
    for i in range(1, len(pts) - 1):
        d = abs(ab.x * (pts[i] - a).y - ab.y * (pts[i] - a).x) / L if L > 1e-12 else (pts[i] - a).length
        if d > best:
            best, bi = d, i
    if best > tol:
        return _dp(pts[:bi + 1], tol)[:-1] + _dp(pts[bi:], tol)
    return [a, b]


def contour(f, bbox, cell=0.003, tol=0.0006):
    """the outline (largest loop, counter-clockwise) of the region f < 0 inside bbox (x0, y0, x1, y1)"""
    x0, y0, x1, y1 = bbox
    nx = int(math.ceil((x1 - x0) / cell)) + 3
    ny = int(math.ceil((y1 - y0) / cell)) + 3
    X = [x0 + (i - 1) * cell for i in range(nx)]
    Y = [y0 + (j - 1) * cell for j in range(ny)]
    F = [[f(Vector((x, y))) for x in X] for y in Y]
    for j in range(ny):
        F[j][0] = F[j][-1] = 1.0
    for i in range(nx):
        F[0][i] = F[-1][i] = 1.0
    segs = {}

    def ip(pa, pb, fa, fb):
        t = fa / (fa - fb)
        return (pa[0] + (pb[0] - pa[0]) * t, pa[1] + (pb[1] - pa[1]) * t)
    for j in range(ny - 1):
        for i in range(nx - 1):
            c = [(X[i], Y[j], F[j][i]), (X[i + 1], Y[j], F[j][i + 1]), (X[i + 1], Y[j + 1], F[j + 1][i + 1]), (X[i], Y[j + 1], F[j + 1][i])]
            idx = sum((1 << k) for k in range(4) if c[k][2] < 0)
            if idx in (0, 15):
                continue
            # edge points: e0 bottom (0-1), e1 right (1-2), e2 top (2-3), e3 left (3-0)
            def E(k):
                a, b = c[k], c[(k + 1) % 4]
                return ip(a[:2], b[:2], a[2], b[2])
            table = {1: [(3, 0)], 2: [(0, 1)], 3: [(3, 1)], 4: [(1, 2)], 6: [(0, 2)], 7: [(3, 2)], 8: [(2, 3)],
                     9: [(2, 0)], 11: [(2, 1)], 12: [(1, 3)], 13: [(1, 0)], 14: [(0, 3)]}
            if idx in (5, 10):
                centre = sum(q[2] for q in c) / 4
                if idx == 5:
                    pairs = [(3, 2), (1, 0)] if centre < 0 else [(3, 0), (1, 2)]
                else:
                    pairs = [(0, 3), (2, 1)] if centre < 0 else [(0, 1), (2, 3)]
            else:
                pairs = table[idx]
            for a, b in pairs:
                pa, pb = E(a), E(b)
                ka = (round(pa[0], 7), round(pa[1], 7))
                kb = (round(pb[0], 7), round(pb[1], 7))
                segs[ka] = kb
    loops = []
    used = set()
    for start in list(segs.keys()):
        if start in used:
            continue
        loop = [start]
        used.add(start)
        cur = segs.get(start)
        while cur is not None and cur != start and cur not in used:
            loop.append(cur)
            used.add(cur)
            cur = segs.get(cur)
        if len(loop) > 8:
            loops.append(loop)
    if not loops:
        return []
    best = max(loops, key=lambda L: abs(area2([Vector(p) for p in L])))
    P = [Vector(p) for p in best]
    if area2(P) < 0:
        P.reverse()
    # simplify (split the loop at two far points)
    k = max(range(len(P)), key=lambda i: (P[i] - P[0]).length)
    A = _dp(P[:k + 1], tol)
    B = _dp(P[k:] + [P[0]], tol)
    out = A[:-1] + B[:-1]
    return out