"""Geometry library for the machine parts kit (models/machine-parts-kit), 2026-10-14.

The user sent a sheet of 39 tan-painted mechanical parts seen from above and asked for each one to be built in one
Blender file through the MCP. The sheet has no scale: 1 sheet pixel = 3 mm (U below), so the big plate is 0.83 m long
and the small capsule 0.15 m.

Item space, metres: X = right on the sheet, Y = up the sheet, Z = up towards the camera; the floor at z = 0; the
item's origin at its centre on the floor. Everything is a closed shell built in one shared bmesh per item (static
kitbash shells may overlap each other; nothing is cut with booleans). Rounded edges are real bevels with 2 segments
(quads, quad corner patches). Caps of round parts are poked (triangle fans); any other n-gon left at the end is
triangulated, which is the only place triangles appear.
"""
import math
import bmesh
from mathutils import Vector, Matrix

U = 0.003                     # metres per sheet pixel
PAINT, DARK, STEEL = 0, 1, 2  # material slots of every item
AX = {'X': (1, 0, 0), '-X': (-1, 0, 0), 'Y': (0, 1, 0), '-Y': (0, -1, 0), 'Z': (0, 0, 1), '-Z': (0, 0, -1)}


def V(*a):
    return Vector(a[0] if len(a) == 1 else a)


def axis_mat(axis):
    """rotation taking local +Z onto axis"""
    d = Vector(AX[axis]) if isinstance(axis, str) else Vector(axis).normalized()
    return Vector((0, 0, 1)).rotation_difference(d).to_matrix().to_4x4()


def fillet(pts, r, n=4, closed=False):
    """round every corner of a polyline (2D or 3D) with radius r (scalar or per-point list)"""
    P = [Vector(p) for p in pts]
    N = len(P)
    R = r if isinstance(r, (list, tuple)) else [r] * N
    out = []
    for i in range(N):
        if R[i] <= 0 or (not closed and i in (0, N - 1)):
            out.append(P[i]); continue
        a, b = P[i - 1] - P[i], P[(i + 1) % N] - P[i]
        if a.length < 1e-9 or b.length < 1e-9 or abs(a.normalized().dot(b.normalized())) > 0.9999:
            out.append(P[i]); continue
        ri = min(R[i], a.length * 0.49, b.length * 0.49)
        p0, p2 = P[i] + a.normalized() * ri, P[i] + b.normalized() * ri
        for k in range(n + 1):           # quadratic Bezier through the corner: smooth and simple
            t = k / n
            out.append(p0 * (1 - t) ** 2 + P[i] * 2 * t * (1 - t) + p2 * t * t)
    return out


def rrect(w, h, r, n=3):
    """rounded rectangle outline, CCW, centred"""
    x, y = w / 2, h / 2
    return [p.to_2d() for p in fillet([(-x, -y, 0), (x, -y, 0), (x, y, 0), (-x, y, 0)], r, n, closed=True)]


class Kit:
    """one item's mesh under construction"""

    def __init__(self):
        self.bm = bmesh.new()
        self.boxes = []          # (centre, size) world AABBs of solid primitives, for collision
        self.sockets = {}

    # ------------------------------------------------------------- helpers
    def _mat(self, faces, mat):
        for f in faces:
            f.material_index = mat

    def _note(self, verts, solid=True):
        if not solid or not verts:
            return
        lo = Vector([min(v.co[i] for v in verts) for i in range(3)])
        hi = Vector([max(v.co[i] for v in verts) for i in range(3)])
        self.boxes.append(((lo + hi) / 2, hi - lo))

    def _merge(self, t, mat):
        """copy a finished temporary bmesh into the item's bmesh, all faces given material mat"""
        vm = {v: self.bm.verts.new(v.co) for v in t.verts}
        out = []
        for f in t.faces:
            nf = self.bm.faces.new([vm[v] for v in f.verts])
            nf.material_index = mat
            out.append(nf)
        t.free()
        return out

    def _bevel_in(self, t, edges, r, seg=2):
        if r <= 0 or not edges:
            return
        bmesh.ops.bevel(t, geom=list(edges), offset=r, offset_type='OFFSET', segments=seg, profile=0.5,
                        affect='EDGES', clamp_overlap=True)

    # ------------------------------------------------------------- boxes
    def box(self, c, s, r=0.0, mat=PAINT, rot=None, solid=True):
        """box centred at c, size s (metres), all edges bevelled by r, rot = optional 4x4 rotation"""
        M = Matrix.Translation(V(c)) @ (rot or Matrix()) @ Matrix.Diagonal((s[0], s[1], s[2], 1))
        t = bmesh.new()
        bmesh.ops.create_cube(t, size=1.0, matrix=M)
        self._bevel_in(t, t.edges, min(r, min(s) * 0.45))
        fs = self._merge(t, mat)
        self._note(list({v for f in fs for v in f.verts}), solid)
        return fs

    def slab(self, outline, t, z0=0.0, mat=PAINT, M=None, r=0.0):
        """2D outline (CCW, metres) extruded from z0 to z0+t, then M applied. r bevels the top and bottom rims."""
        b = bmesh.new()
        vb = [b.verts.new((p[0], p[1], z0)) for p in outline]
        vt = [b.verts.new((p[0], p[1], z0 + t)) for p in outline]
        n = len(outline)
        for i in range(n):
            b.faces.new((vb[i], vb[(i + 1) % n], vt[(i + 1) % n], vt[i]))
        top = b.faces.new(vt)
        bot = b.faces.new(list(reversed(vb)))
        if M is not None:
            bmesh.ops.transform(b, matrix=M, verts=vb + vt)
        self._bevel_in(b, list(top.edges) + list(bot.edges), min(r, t * 0.45))
        fs = self._merge(b, mat)
        self._note(list({v for f in fs for v in f.verts}))
        return fs

    # ------------------------------------------------------------- round things
    def lathe(self, prof, c, axis='Z', segs=32, mat=PAINT, spin=0.0, solid=True):
        """surface of revolution. prof = [(radius, height)...] from one end to the other; radius 0 = pole (fan)."""
        bm = self.bm
        M = Matrix.Translation(V(c)) @ axis_mat(axis) @ Matrix.Rotation(spin, 4, 'Z')
        rings = []
        for (r, z) in prof:
            if r < 1e-7:
                rings.append([bm.verts.new(M @ Vector((0, 0, z)))])
            else:
                rings.append([bm.verts.new(M @ Vector((r * math.cos(2 * math.pi * i / segs),
                                                        r * math.sin(2 * math.pi * i / segs), z))) for i in range(segs)])
        fs = []
        for a, b in zip(rings, rings[1:]):
            if len(a) == 1 and len(b) == 1:
                continue
            for i in range(segs):
                j = (i + 1) % segs
                if len(a) == 1:
                    fs.append(bm.faces.new((a[0], b[j], b[i])))
                elif len(b) == 1:
                    fs.append(bm.faces.new((a[i], a[j], b[0])))
                else:
                    fs.append(bm.faces.new((a[i], a[j], b[j], b[i])))
        # profiles that run from the pole outwards and back give outward normals with this winding when z rises;
        # make it right whatever the direction
        bmesh.ops.recalc_face_normals(bm, faces=fs)
        self._mat(fs, mat)
        self._note([v for r in rings for v in r], solid)
        return fs

    def cyl(self, c, r, h, axis='Z', segs=24, bev=0.0, mat=PAINT, spin=0.0, solid=True):
        """closed cylinder from c along axis for length h, rims rounded by bev (lathe profile, poked caps)"""
        b = min(bev, r * 0.45, h * 0.45)
        if b > 0:
            k = 0.29289 * b
            prof = [(0, 0), (r - b, 0), (r - k, k), (r, b), (r, h - b), (r - k, h - k), (r - b, h), (0, h)]
        else:
            prof = [(0, 0), (r, 0), (r, h), (0, h)]
        return self.lathe(prof, c, axis, segs, mat, spin, solid)

    def can(self, c, r, L, axis='Z', end=0.35, segs=32, grooves=(), mat=PAINT, ends=(True, True)):
        """a tank: cylinder with domed (elliptic) ends; grooves = distances from the start"""
        e = r * end
        prof = []
        if ends[0]:
            prof += [(0, 0)] + [(r * math.sin(a), e * (1 - math.cos(a))) for a in
                                (math.pi / 2 * i / 5 for i in range(1, 6))]
        else:
            prof += [(0, 0), (r * 0.97, 0), (r, r * 0.03)]
        for g in sorted(grooves):
            prof += [(r, g - 0.006), (r * 0.94, g - 0.003), (r * 0.94, g + 0.003), (r, g + 0.006)]
        if ends[1]:
            prof += [(r * math.sin(a), L - e * (1 - math.cos(a))) for a in
                     (math.pi / 2 * (5 - i) / 5 for i in range(0, 5))] + [(0, L)]
        else:
            prof += [(r, L - r * 0.03), (r * 0.97, L), (0, L)]
        return self.lathe(prof, c, axis, segs, mat)

    def ring(self, c, r, w, h, axis='Z', segs=32, mat=PAINT):
        """flat washer / band: outer radius r, radial width w, height h (a closed square-section torus)"""
        prof = [(r - w, 0), (r, 0), (r, h), (r - w, h), (r - w, 0)]
        bm = self.bm
        M = Matrix.Translation(V(c)) @ axis_mat(axis)
        rings = [[bm.verts.new(M @ Vector((rr * math.cos(2 * math.pi * i / segs), rr * math.sin(2 * math.pi * i / segs), z)))
                  for i in range(segs)] for rr, z in prof[:-1]]
        fs = []
        for k in range(4):
            a, b = rings[k], rings[(k + 1) % 4]
            for i in range(segs):
                j = (i + 1) % segs
                fs.append(bm.faces.new((a[i], a[j], b[j], b[i])))
        bmesh.ops.recalc_face_normals(bm, faces=fs)
        self._mat(fs, mat)
        return fs

    def bolt(self, c, r=0.006, h=0.005, axis='Z', mat=STEEL):
        """hex-head bolt standing on c along axis"""
        return self.cyl(c, r, h, axis, segs=6, bev=r * 0.18, mat=mat, spin=math.pi / 6, solid=False)

    def dot(self, c, r, axis='Z', h=0.002, mat=DARK, segs=16):
        """a dark port/disc, proud by h (a hole read from above)"""
        return self.cyl(c, r, h, axis, segs=segs, bev=0, mat=mat, solid=False)

    def socket(self, c, r, depth, axis='Z', rim=0.006, segs=32, mat=PAINT):
        """a round open socket rising from a face: rim ring outside, bore inside down to a dark floor"""
        R2 = r + rim
        k = rim * 0.3
        f1 = self.lathe([(0, 0), (R2, 0), (R2, depth - k), (R2 - k, depth), (r + k, depth), (r, depth - k), (r, 0.002),
                         (0, 0.002)], c, axis, segs, mat, solid=False)
        d = Vector(AX[axis]) if isinstance(axis, str) else Vector(axis).normalized()
        f2 = self.dot(V(c) + d * 0.002, r * 0.999, axis, h=0.0015, mat=DARK, segs=segs)
        return f1 + f2

    # ------------------------------------------------------------- swept tube
    def tube(self, pts, r, bend=None, segs=12, mat=PAINT, closed=False, caps=True):
        """pipe of radius r along a polyline (metres), corners bent with radius `bend` (default 2.5 r)"""
        bm = self.bm
        P = fillet(pts, (bend if bend is not None else r * 2.5), 6, closed) if len(pts) > 2 else [Vector(p) for p in pts]
        N = len(P)
        T = []
        for i in range(N):
            if closed:
                t = P[(i + 1) % N] - P[i - 1]
            else:
                t = (P[min(i + 1, N - 1)] - P[max(i - 1, 0)])
            T.append(t.normalized())
        # parallel transport
        up = Vector((0, 0, 1)) if abs(T[0].z) < 0.9 else Vector((1, 0, 0))
        n0 = (up - T[0] * up.dot(T[0])).normalized()
        Ns = [n0]
        for i in range(1, N):
            q = T[i - 1].rotation_difference(T[i])
            n = q @ Ns[-1]
            Ns.append((n - T[i] * n.dot(T[i])).normalized())
        rings = []
        for i in range(N):
            b = T[i].cross(Ns[i])
            rings.append([bm.verts.new(P[i] + (Ns[i] * math.cos(2 * math.pi * k / segs) + b * math.sin(2 * math.pi * k / segs)) * r)
                          for k in range(segs)])
        fs = []
        last = N if closed else N - 1
        for i in range(last):
            a, b2 = rings[i], rings[(i + 1) % N]
            for k in range(segs):
                j = (k + 1) % segs
                fs.append(bm.faces.new((a[k], a[j], b2[j], b2[k])))
        if caps and not closed:
            for ring_, sgn in ((rings[0], -1), (rings[-1], 1)):
                f = bm.faces.new(ring_)
                fs += bmesh.ops.poke(bm, faces=[f])['faces']
        bmesh.ops.recalc_face_normals(bm, faces=fs)
        self._mat(fs, mat)
        return fs

    def handle(self, a, b, up, h, r=0.006, mat=PAINT, feet=True):
        """U handle standing from point a to point b, lifted h along `up`"""
        a, b, u = V(a), V(b), Vector(AX[up]) if isinstance(up, str) else V(up)
        fs = self.tube([a, a + u * h, b + u * h, b], r, bend=min(h * 0.45, r * 3), mat=mat)
        if feet:
            for p in (a, b):
                self.cyl(p - u * 0.002, r * 1.7, r * 1.2, axis=tuple(u), segs=12, bev=r * 0.3, mat=mat, solid=False)
        return fs

    # ------------------------------------------------------------- plate with real holes
    def holed_plate(self, xs, ys, holes, t, z0=0.0, k=4, mat=PAINT, M=None, wall_mat=None):
        """flat plate from cell edges xs, ys (metres, ascending). holes = {(i, j): (a, b, n)}: a superellipse hole of
        half sizes a, b, exponent n (2 = ellipse, 4+ = rounded slot) centred in cell (i, j). Every cell side carries
        k segments, so solid cells are k x k quad grids and hole cells a ring of 4k quads: no booleans, all quads.
        The hole walls come from extruding the plate (they are boundary edges)."""
        bm = self.bm
        nx, ny = len(xs) - 1, len(ys) - 1
        gx = [xs[i] + (xs[i + 1] - xs[i]) * s / k for i in range(nx) for s in range(k)] + [xs[-1]]
        gy = [ys[j] + (ys[j + 1] - ys[j]) * s / k for j in range(ny) for s in range(k)] + [ys[-1]]
        vert = {}

        def gv(ix, iy):
            key = (ix, iy)
            if key not in vert:
                vert[key] = bm.verts.new((gx[ix], gy[iy], z0))
            return vert[key]

        faces = []
        for i in range(nx):
            for j in range(ny):
                ox, oy = i * k, j * k
                if (i, j) not in holes:
                    for a in range(k):
                        for b in range(k):
                            faces.append(bm.faces.new((gv(ox + a, oy + b), gv(ox + a + 1, oy + b),
                                                       gv(ox + a + 1, oy + b + 1), gv(ox + a, oy + b + 1))))
                    continue
                ha, hb, n = holes[(i, j)]
                cx, cy = (xs[i] + xs[i + 1]) / 2, (ys[j] + ys[j + 1]) / 2
                loop = [(ox + s, oy) for s in range(k)] + [(ox + k, oy + s) for s in range(k)] + \
                       [(ox + k - s, oy + k) for s in range(k)] + [(ox, oy + k - s) for s in range(k)]
                outer = [gv(*q) for q in loop]
                inner = []
                for v in outer:
                    d = Vector((v.co.x - cx, v.co.y - cy))
                    th = math.atan2(d.y, d.x)
                    # even angular spacing looks better than the projected square: blend half way
                    rho = (abs(math.cos(th) / ha) ** n + abs(math.sin(th) / hb) ** n) ** (-1.0 / n)
                    inner.append(bm.verts.new((cx + rho * math.cos(th), cy + rho * math.sin(th), z0)))
                L = len(outer)
                for s in range(L):
                    faces.append(bm.faces.new((outer[s], outer[(s + 1) % L], inner[(s + 1) % L], inner[s])))
        ext = bmesh.ops.extrude_face_region(bm, geom=faces)
        nv = [e for e in ext['geom'] if isinstance(e, bmesh.types.BMVert)]
        bmesh.ops.translate(bm, vec=(0, 0, t), verts=nv)
        newf = [e for e in ext['geom'] if isinstance(e, bmesh.types.BMFace)]
        allf = list({f for v in list(vert.values()) + nv for f in v.link_faces})
        bmesh.ops.recalc_face_normals(bm, faces=allf)
        if M is not None:
            bmesh.ops.transform(bm, matrix=M, verts=list({v for f in allf for v in f.verts}))
        self._mat(allf, mat)
        if wall_mat is not None:     # hole walls darker (they are in shade on the sheet)
            for f in allf:
                if abs(f.normal.z) < 0.5 and not self._outer_wall(f, xs, ys, M):
                    f.material_index = wall_mat
        self._note(list({v for f in allf for v in f.verts}))
        return allf

    @staticmethod
    def _outer_wall(f, xs, ys, M):
        c = f.calc_center_median()
        if M is not None:
            c = M.inverted() @ c
        e = 1e-5
        return abs(c.x - xs[0]) < e or abs(c.x - xs[-1]) < e or abs(c.y - ys[0]) < e or abs(c.y - ys[-1]) < e

    # ------------------------------------------------------------- hidden faces
    @staticmethod
    def cull_hidden(bm, floor=True):
        """delete faces nobody can see (memory: game-ready-model-standards, 'drop hidden contact faces'): a face is
        hidden when its centre, corners (pulled 15 % in) and edge middles, nudged 0.4 mm out along its normal, all lie
        inside one and the same OTHER closed shell of the same mesh (ray parity per shell, two ray directions agreeing), or when it faces down on the floor (z = 0).
        Faces only partly buried are kept, so no hole can open."""
        from mathutils.bvhtree import BVHTree
        bm.faces.ensure_lookup_table()
        shell = {}
        sid = 0
        for f in bm.faces:
            if f in shell:
                continue
            stack = [f]
            shell[f] = sid
            while stack:
                g = stack.pop()
                for e in g.edges:
                    for h in e.link_faces:
                        if h not in shell:
                            shell[h] = sid
                            stack.append(h)
            sid += 1
        fshell = [shell[f] for f in bm.faces]
        tree = BVHTree.FromBMesh(bm)
        dirs = (Vector((0.5773, 0.5779, 0.5767)).normalized(), Vector((-0.3412, 0.8123, -0.4731)).normalized())

        def odd_shells(p, own, d):
            odd = {}
            o = p.copy()
            for _ in range(64):
                loc, nrm, idx, dist = tree.ray_cast(o, d)
                if loc is None:
                    break
                s = fshell[idx]
                if s != own:
                    odd[s] = not odd.get(s, False)
                o = loc + d * 1e-6
            return {s for s, v in odd.items() if v}

        def inside(p, own):
            a = odd_shells(p, own, dirs[0])
            return a & odd_shells(p, own, dirs[1]) if a else a

        dead = []
        for f in bm.faces:
            n = f.normal
            if floor and n.z < -0.99 and all(v.co.z < 3e-4 for v in f.verts):
                dead.append(f)
                continue
            c = f.calc_center_median()
            vs = [v.co for v in f.verts]
            pts = [c] + [c.lerp(v, 0.85) for v in vs] + [c.lerp((a + b) / 2, 0.9) for a, b in zip(vs, vs[1:] + vs[:1])]
            own = shell[f]
            common = None
            for p in pts:
                s = inside(p + n * 4e-4, own)
                common = s if common is None else common & s
                if not common:
                    break
            if common:              # every sample inside one and the same other shell
                dead.append(f)
        if dead:
            bmesh.ops.delete(bm, geom=dead, context='FACES')
        return len(dead)

    # ------------------------------------------------------------- finish
    def finish(self, cull=True, floor=True):
        bm = self.bm
        # no global weld: shells that touch must stay separate (welding them made non-manifold edges)
        ng = [f for f in bm.faces if len(f.verts) > 4]
        if ng:
            bmesh.ops.triangulate(bm, faces=ng, quad_method='BEAUTY', ngon_method='BEAUTY')
        bmesh.ops.dissolve_degenerate(bm, dist=1e-7, edges=bm.edges)
        bm.normal_update()
        self.culled = self.cull_hidden(bm, floor) if cull else 0
        bm.normal_update()
        return bm


def audit(me):
    """quads / tris / ngons / open edges / non-manifold of a mesh"""
    bm = bmesh.new()
    bm.from_mesh(me)
    q = sum(1 for f in bm.faces if len(f.verts) == 4)
    t = sum(1 for f in bm.faces if len(f.verts) == 3)
    n = sum(1 for f in bm.faces if len(f.verts) > 4)
    op = sum(1 for e in bm.edges if e.is_boundary)
    nm = sum(1 for e in bm.edges if not e.is_manifold and not e.is_boundary)
    za = sum(1 for f in bm.faces if f.calc_area() < 1e-10)
    bm.free()
    return {"faces": len(me.polygons), "quads": q, "tris": t, "ngons": n, "open": op, "nonmanifold": nm, "zero": za}