"""Geometry helpers for the convict body build.

Pure bmesh: lofted shells with raised panels, lathes, rounded boxes, plates with thickness, hoses.
Everything is quads except fan caps at poles and cylinder ends. No boolean modifiers anywhere.
Every vertex is weighted to one bone (or blended between two) as it is made.
"""
import json
import bmesh
import bpy
from bisect import bisect_right
from contextlib import contextmanager
from math import sin, cos, radians, degrees, pi
from mathutils import Vector, Matrix
from mathutils.bvhtree import BVHTree

V = Vector


def frame(o, z, xh):
    """Right-handed frame: local Z along z, local X as near xh as possible."""
    z = Vector(z).normalized(); xh = Vector(xh)
    x = (xh - xh.dot(z) * z).normalized(); y = z.cross(x)
    return Matrix(((x.x, y.x, z.x, o[0]), (x.y, y.y, z.y, o[1]), (x.z, y.z, z.z, o[2]), (0, 0, 0, 1)))


def T(x, y, z): return Matrix.Translation(Vector((x, y, z)))
def RX(d): return Matrix.Rotation(radians(d), 4, 'X')
def RY(d): return Matrix.Rotation(radians(d), 4, 'Y')
def RZ(d): return Matrix.Rotation(radians(d), 4, 'Z')


def K(z, rx, ry=None, n=2.0, cx=0.0, cy=0.0):
    """Profile key: superellipse half sizes rx, ry, exponent n, centre offset."""
    return (z, rx, rx if ry is None else ry, n, cx, cy)


def P(a0, a1, z0=None, z1=None, h=0.0, mat=None, sym=False, eps=None):
    """Panel on a shell: angle range (degrees from the front; or 0..1 across a plate), z range, height (additive)."""
    return dict(a=(a0, a1), z=(z0, z1), h=h, mat=mat, sym=sym, eps=eps)


def sup(t, rx, ry, n):
    c = abs(cos(t)) / max(rx, 1e-7); s = abs(sin(t)) / max(ry, 1e-7)
    return (c ** n + s ** n) ** (-1.0 / n)


class Prof:
    """Cross-section along local Z. front = angle of the 'front' direction in the local XY plane, pdir = turn sense."""

    def __init__(self, keys, front=90.0, pdir=-1, lin=False):
        ks = sorted(keys, key=lambda k: k[0])
        self.z = [k[0] for k in ks]; self.v = [tuple(k[1:]) for k in ks]
        self.front, self.pdir, self.lin = front, pdir, lin
        n = len(ks); self.m = []
        for i in range(n):
            mi = []
            for c in range(5):
                if n == 1: mi.append(0.0); continue
                if i == 0: s = (self.v[1][c] - self.v[0][c]) / (self.z[1] - self.z[0])
                elif i == n - 1: s = (self.v[-1][c] - self.v[-2][c]) / (self.z[-1] - self.z[-2])
                else:
                    a = (self.v[i][c] - self.v[i - 1][c]) / (self.z[i] - self.z[i - 1])
                    b = (self.v[i + 1][c] - self.v[i][c]) / (self.z[i + 1] - self.z[i])
                    s = 0.0 if a * b <= 0 else 2 * a * b / (a + b)
                mi.append(s)
            self.m.append(mi)

    def ev(self, z):
        zs = self.z
        if z <= zs[0]: return self.v[0]
        if z >= zs[-1]: return self.v[-1]
        i = bisect_right(zs, z) - 1
        h = zs[i + 1] - zs[i]; t = (z - zs[i]) / h
        if self.lin: return tuple(a + (b - a) * t for a, b in zip(self.v[i], self.v[i + 1]))
        t2 = t * t; t3 = t2 * t
        h00 = 2 * t3 - 3 * t2 + 1; h10 = t3 - 2 * t2 + t; h01 = -2 * t3 + 3 * t2; h11 = t3 - t2
        return tuple(h00 * self.v[i][c] + h10 * h * self.m[i][c] + h01 * self.v[i + 1][c] + h11 * h * self.m[i + 1][c]
                     for c in range(5))

    def th(self, phi): return radians(self.front + self.pdir * phi)

    def pt(self, phi, z, off=0.0):
        rx, ry, n, cx, cy = self.ev(z); t = self.th(phi)
        r = sup(t, rx, ry, n) + off
        return Vector((cx + r * cos(t), cy + r * sin(t), z))

    def nrm(self, phi, z):
        a = self.pt(phi + 0.5, z) - self.pt(phi - 0.5, z); b = self.pt(phi, z + 1e-3) - self.pt(phi, z - 1e-3)
        n = a.cross(b).normalized(); rad = self.pt(phi, z, 0.01) - self.pt(phi, z)
        return n if n.dot(rad) > 0 else -n


class Builder:
    def __init__(self, bones, mats):
        self.bm = bmesh.new(); bm = self.bm
        self.bone_names = list(bones); self.bones = {b: i for i, b in enumerate(bones)}
        self.mat_names = list(mats); self.mats = {m: i for i, m in enumerate(mats)}
        self.l_b1 = bm.verts.layers.int.new('b1'); self.l_b2 = bm.verts.layers.int.new('b2')
        self.l_w = bm.verts.layers.float.new('w2')
        self.l_part = bm.faces.layers.int.new('part'); self.l_uv = bm.loops.layers.uv.new('UVMap')
        self.G = Matrix.Identity(4); self.L = Matrix.Identity(4); self.sfx = '_l'
        self.bone = 'root'; self.mat = mats[0]; self.wfn = None
        self.parts = []

    @property
    def M(self): return self.G @ self.L

    @contextmanager
    def at(self, M):
        old = self.L; self.L = old @ M
        try: yield
        finally: self.L = old

    @contextmanager
    def side(self, s):
        og, osf, ol = self.G, self.sfx, self.L
        self.G = Matrix.Scale(-1, 4, (1, 0, 0)) if s < 0 else Matrix.Identity(4)
        self.sfx = '_r' if s < 0 else '_l'; self.L = Matrix.Identity(4)
        try: yield
        finally: self.G, self.sfx, self.L = og, osf, ol

    def _vert(self, p, w=None):
        v = self.bm.verts.new(self.M @ Vector(p))
        if w is not None: b1, b2, ww = w
        elif self.wfn: b1, b2, ww = self.wfn(v.co)
        else: b1, b2, ww = self.bone, self.bone, 0.0
        v[self.l_b1] = self.bones[b1]; v[self.l_b2] = self.bones[b2]; v[self.l_w] = ww
        return v

    def _face(self, vl, uvl, mat, pid):
        try: f = self.bm.faces.new(vl)
        except ValueError: return None
        f.material_index = self.mats[mat or self.mat]; f[self.l_part] = pid; f.smooth = True
        for lp, uv in zip(f.loops, uvl): lp[self.l_uv].uv = uv
        return f

    def grid(self, rows, closed=True, pole0=None, pole1=None, cap0=False, cap1=False, cellmat=None, name='part',
             wrows=None):
        """rows[j][i] local points. Quads between rows; optional pole point or centroid fan at either end."""
        pid = len(self.parts); self.parts.append((name, self.bone if not self.wfn else 'spine'))
        rows = [[Vector(p) for p in r] for r in rows]
        nr = len(rows); nc = len(rows[0]); ncell = nc if closed else nc - 1
        vs = [[self._vert(p, wrows[j] if wrows else None) for p in row] for j, row in enumerate(rows)]
        U = [0.0]
        for i in range(ncell):
            i2 = (i + 1) % nc
            U.append(U[-1] + sum((rows[j][i2] - rows[j][i]).length for j in range(nr)) / nr)
        W = [0.0]
        for j in range(nr - 1):
            W.append(W[-1] + sum((rows[j + 1][i] - rows[j][i]).length for i in range(nc)) / nc)
        for j in range(nr - 1):
            for i in range(ncell):
                i2 = (i + 1) % nc
                self._face((vs[j][i], vs[j][i2], vs[j + 1][i2], vs[j + 1][i]),
                           ((U[i], W[j]), (U[i + 1], W[j]), (U[i + 1], W[j + 1]), (U[i], W[j + 1])),
                           cellmat(j, i) if cellmat else None, pid)
        for end, pole, cap in ((0, pole0, cap0), (nr - 1, pole1, cap1)):
            if pole is None and not cap: continue
            c = Vector(pole) if pole is not None else sum(rows[end], Vector()) / nc
            cv = self._vert(c, wrows[end] if wrows else None)
            r = max((p - c).length for p in rows[end]); ring = vs[end]
            jm = 0 if end == 0 else max(nr - 2, 0)
            for i in range(ncell):
                i2 = (i + 1) % nc
                ua = (r * cos(2 * pi * i / ncell), r * sin(2 * pi * i / ncell))
                ub = (r * cos(2 * pi * (i + 1) / ncell), r * sin(2 * pi * (i + 1) / ncell))
                m = cellmat(jm, i) if (cellmat and nr > 1) else None
                if end == 0: self._face((cv, ring[i2], ring[i]), ((0, 0), ub, ua), m, pid)
                else: self._face((cv, ring[i], ring[i2]), ((0, 0), ua, ub), m, pid)
        return pid, vs

    def lathe(self, prof, seg=14, name='lathe'):
        """Revolve (z, r) around local Z. r = 0 at an end makes a pole."""
        prof = [(z, max(r, 0.0)) for z, r in prof]; pole0 = pole1 = None
        if prof[0][1] < 1e-6: pole0 = (0, 0, prof[0][0]); prof = prof[1:]
        if prof[-1][1] < 1e-6: pole1 = (0, 0, prof[-1][0]); prof = prof[:-1]
        rows = [[(r * cos(2 * pi * i / seg), r * sin(2 * pi * i / seg), z) for i in range(seg)] for z, r in prof]
        return self.grid(rows, True, pole0, pole1, cap0=pole0 is None, cap1=pole1 is None, name=name)

    def cyl(self, p0, p1, r, seg=12, ch=0.0, r1=None, name='cyl'):
        p0 = Vector(p0); p1 = Vector(p1); d = p1 - p0; L = d.length
        xh = Vector((1, 0, 0)) if abs(d.normalized().x) < 0.9 else Vector((0, 1, 0))
        r1 = r if r1 is None else r1
        with self.at(frame(p0, d, xh)):
            if ch > 0: prof = [(0, 0), (0, r - ch), (ch, r), (L - ch, r1), (L, r1 - ch), (L, 0)]
            else: prof = [(0, 0), (0, r), (L, r1), (L, 0)]
            return self.lathe(prof, seg, name)

    def tube(self, z0, z1, r, seg=14, ch=0.002, r1=None, name='tube'):
        return self.cyl((0, 0, z0), (0, 0, z1), r, seg, ch, r1, name)

    def ball(self, c, r, seg=16, rings=8, name='ball'):
        with self.at(Matrix.Translation(Vector(c))):
            return self.lathe([(-r * cos(pi * k / rings), r * sin(pi * k / rings)) for k in range(rings + 1)], seg, name)

    def rbox(self, c, size, n=6.0, seg=16, ch=0.003, t0=(1.0, 1.0), t1=(1.0, 1.0), rot=None, name='box'):
        """Rounded box lofted along local Z. t0/t1 scale the -Z / +Z end (wedges)."""
        hx, hy, hz = size[0] / 2, size[1] / 2, size[2] / 2
        ch = min(ch, hx * 0.45, hy * 0.45, hz * 0.45)
        M = Matrix.Translation(Vector(c))
        if rot is not None: M = M @ rot

        def ring(z, shrink):
            t = (z + hz) / (2 * hz); sx = t0[0] + (t1[0] - t0[0]) * t; sy = t0[1] + (t1[1] - t0[1]) * t
            ax = max(hx * sx - shrink, 1e-4); ay = max(hy * sy - shrink, 1e-4); out = []
            for i in range(seg):
                a = 2 * pi * i / seg; ca, sa = cos(a), sin(a)
                out.append((ax * (1 if ca >= 0 else -1) * abs(ca) ** (2 / n), ay * (1 if sa >= 0 else -1) * abs(sa) ** (2 / n), z))
            return out
        with self.at(M):
            return self.grid([ring(-hz, ch), ring(-hz + ch, 0), ring(hz - ch, 0), ring(hz, ch)], True, cap0=True,
                             cap1=True, name=name)

    def hose(self, pts, r, seg=8, bones=None, sub=6, rib=0.0, name='hose'):
        """Tube through the points (Catmull-Rom). bones: one bone per point, blended along each span."""
        pts = [Vector(p) for p in pts]; n = len(pts)

        def crp(i, t):
            p0 = pts[max(i - 1, 0)]; p1 = pts[i]; p2 = pts[min(i + 1, n - 1)]; p3 = pts[min(i + 2, n - 1)]
            return 0.5 * ((2 * p1) + (-p0 + p2) * t + (2 * p0 - 5 * p1 + 4 * p2 - p3) * t * t + (-p0 + 3 * p1 - 3 * p2 + p3) * t ** 3)
        cs = []; ws = []
        for i in range(n - 1):
            for k in range(sub): cs.append(crp(i, k / sub)); ws.append((i, k / sub))
        cs.append(pts[-1]); ws.append((n - 2, 1.0))
        tg = [(cs[min(k + 1, len(cs) - 1)] - cs[max(k - 1, 0)]).normalized() for k in range(len(cs))]
        nr = tg[0].orthogonal().normalized(); rows = []; wr = []
        for k, (c, t) in enumerate(zip(cs, tg)):
            nr = (nr - nr.dot(t) * t).normalized(); bn = t.cross(nr)
            rr = r * (1 + (rib if k % 2 else 0.0))
            rows.append([c + rr * (cos(2 * pi * i / seg) * nr + sin(2 * pi * i / seg) * bn) for i in range(seg)])
            if bones: i, tt = ws[k]; wr.append((bones[i], bones[i + 1], tt))
        return self.grid(rows, True, cap0=True, cap1=True, name=name, wrows=wr if bones else None)

    def bolt(self, p, n, r=0.005, h=0.003, seg=6, name='bolt'):
        p = Vector(p); n = Vector(n).normalized()
        return self.cyl(p - n * 0.001, p + n * h, r, seg, ch=min(r * 0.25, h * 0.4), name=name)


def _nrm180(a): return (a + 180.0) % 360.0 - 180.0
def _angd(a, b): return abs(_nrm180(a - b))


def shell(B, prof, panels=(), seg=24, dz=0.02, z0=None, z1=None, span=None, thick=0.0, off=0.0, rim=0.0,
          cap0=True, cap1=True, zwarp=None, twist=0.0, eps=0.0025, name='shell'):
    """Loft the profile along local Z as one quad grid; panels raise/sink regions with chamfered steps.

    span=None: closed shell all the way round (angles in degrees from the front).
    span=(centre, width[, centre1, width1]): a plate over that angle range (panel angles become 0..1 across it),
    made a closed solid with `thick`; `rim` chamfers its outline.
    """
    Z0 = prof.z[0] if z0 is None else z0; Z1 = prof.z[-1] if z1 is None else z1
    full = span is None
    if not full and len(span) == 2: span = (span[0], span[1], span[0], span[1])
    rr = [(k[0] + k[1]) / 2 for k in prof.v]; rref = max(sum(rr) / len(rr), 1e-4)
    wavg = None if full else radians((span[1] + span[3]) / 2)

    def conv(e): return degrees(e / rref) if full else e / max(rref * wavg, 1e-5)
    pans = []
    for p in panels:
        pans.append(p)
        if p['sym']:
            a0, a1 = p['a']; q = dict(p); q['a'] = (-a1, -a0) if full else (1 - a1, 1 - a0); pans.append(q)
    if rim > 0 and not full:
        pans.insert(0, dict(a=(None, None), z=(None, None), h=rim, mat=None, sym=False, eps=rim, _rim=True)); off -= rim
    must_u = set() if full else {0.0, 1.0}; must_z = {Z0, Z1}; spec = []

    def addz(z):
        if Z0 + 1e-5 < z < Z1 - 1e-5: must_z.add(z)
    for p in pans:
        e = eps if p['eps'] is None else p['eps']
        if p['h'] == 0: e = 0.0
        eu = conv(e); a0, a1 = p['a']; pz0, pz1 = p['z']
        if p.get('_rim'): a0, a1 = eu, 1 - eu; pz0, pz1 = Z0 + e, Z1 - e
        pz0 = Z0 if pz0 is None else max(pz0, Z0); pz1 = Z1 if pz1 is None else min(pz1, Z1)
        ring = full and (a1 - a0) >= 360 - 1e-6
        if not ring:
            for a, sg in ((a0, -1), (a1, 1)):
                if full:
                    must_u.add(_nrm180(a))
                    if eu > 0: must_u.add(_nrm180(a + sg * eu))
                else:
                    if 1e-6 < a < 1 - 1e-6: must_u.add(a)
                    b = a + sg * eu
                    if eu > 0 and -1e-9 <= b <= 1 + 1e-9: must_u.add(min(max(b, 0.0), 1.0))
        for z, sg in ((pz0, -1), (pz1, 1)):
            addz(z)
            if e > 0: addz(z + sg * e)
        spec.append((a0, a1, pz0, pz1, p['h'], p['mat'], eu, e, ring))
    if prof.lin:
        for z in prof.z: addz(z)

    def dedupe(xs, tol):
        out = []
        for x in sorted(xs):
            if not out or x - out[-1] > tol: out.append(x)
        return out
    mu = dedupe(must_u, 1e-5)
    if full:
        if len(mu) > 1 and (mu[0] + 360 - mu[-1]) < 1e-5: mu.pop()
        gap = 0.35 * 360 / seg
        us = mu + [b for b in (-180 + 360 * i / seg for i in range(seg)) if all(_angd(b, m) > gap for m in mu)]
    else:
        gap = 0.35 / seg
        us = mu + [b for b in (i / seg for i in range(seg + 1)) if all(abs(b - m) > gap for m in mu)]
    us = sorted(us)
    mz = dedupe(must_z, 1e-6); nz = max(1, round((Z1 - Z0) / dz))
    zs = sorted(mz + [b for b in (Z0 + (Z1 - Z0) * j / nz for j in range(nz + 1)) if all(abs(b - m) > 0.35 * dz for m in mz)])

    def inside(u, z, s, grow=0.0):
        a0, a1, pz0, pz1, h, m, eu, e, ring = s
        if z < pz0 - grow * e - 1e-7 or z > pz1 + grow * e + 1e-7: return False
        if ring: return True
        if full: return ((u - a0 + grow * eu + 1e-6) % 360.0) <= (a1 - a0) + 2 * grow * eu + 2e-6
        return a0 - grow * eu - 1e-7 <= u <= a1 + grow * eu + 1e-7

    def H(u, z): return sum(s[4] for s in spec if inside(u, z, s))

    def P3(u, z, h, inner=False):
        t = (z - Z0) / (Z1 - Z0) if Z1 > Z0 else 0.0
        if full: phi = u
        else: phi = span[0] + (span[2] - span[0]) * t + (u - 0.5) * (span[1] + (span[3] - span[1]) * t)
        th = prof.th(phi) + twist * t
        rx, ry, n, cx, cy = prof.ev(z)
        rho = sup(th, rx, ry, n) + off + (-thick if inner else h)
        return Vector((cx + rho * cos(th), cy + rho * sin(th), z + (zwarp(phi, t, u) if zwarp else 0.0)))
    cols = us if prof.pdir > 0 else us[::-1]; nc = len(cols)

    def matof(um, zm):
        m = None
        for s in spec:
            if s[5] and inside(um, zm, s, 1.0): m = s[5]
        return m
    if full:
        pole0 = pole1 = None
        if max(prof.ev(Z0)[:2]) < 3e-4: pole0 = P3(0, Z0, 0); pole0.x, pole0.y = prof.ev(Z0)[3:5]; zs = zs[1:]
        if max(prof.ev(Z1)[:2]) < 3e-4: pole1 = P3(0, Z1, 0); pole1.x, pole1.y = prof.ev(Z1)[3:5]; zs = zs[:-1]
        rows = [[P3(u, z, H(u, z)) for u in cols] for z in zs]

        def cm(j, i):
            a = cols[i]; b = cols[(i + 1) % nc]
            return matof(_nrm180(a + _nrm180(b - a) / 2), (zs[j] + zs[min(j + 1, len(zs) - 1)]) / 2)
        return B.grid(rows, True, pole0, pole1, cap0=cap0 and pole0 is None, cap1=cap1 and pole1 is None, cellmat=cm, name=name)
    rows = [[P3(u, z, H(u, z)) for u in cols] + [P3(u, z, 0, True) for u in cols[::-1]] for z in zs]

    def cm2(j, i):
        if i >= nc - 1: return None
        return matof((cols[i] + cols[i + 1]) / 2, (zs[j] + zs[min(j + 1, len(zs) - 1)]) / 2)
    pid, vs = B.grid(rows, True, cellmat=cm2, name=name)
    for j in (0, len(zs) - 1):
        for i in range(nc - 1):
            B._face((vs[j][i], vs[j][i + 1], vs[j][2 * nc - 2 - i], vs[j][2 * nc - 1 - i]),
                    ((0, 0), (0.01, 0), (0.01, thick), (0, thick)), None, pid)
    return pid, vs


def finish(B, name, arm_obj, coll, materials):
    """Turn the builder into a skinned mesh object. Returns (object, audit dict)."""
    bm = B.bm
    bmesh.ops.recalc_face_normals(bm, faces=bm.faces[:])
    for e in bm.edges:
        if len(e.link_faces) == 2 and e.calc_face_angle(0.0) > radians(33): e.smooth = False
    used_b = sorted({v[B.l_b1] for v in bm.verts} | {v[B.l_b2] for v in bm.verts if v[B.l_w] > 1e-4})
    bmap = {b: i for i, b in enumerate(used_b)}
    used_m = sorted({f.material_index for f in bm.faces}); mmap = {m: i for i, m in enumerate(used_m)}
    for f in bm.faces: f.material_index = mmap[f.material_index]
    me = bpy.data.meshes.new(name); ob = bpy.data.objects.new(name, me); coll.objects.link(ob)
    for m in used_m: me.materials.append(materials[B.mat_names[m]])
    for b in used_b: ob.vertex_groups.new(name=B.bone_names[b])
    dl = bm.verts.layers.deform.verify()
    for v in bm.verts:
        b1, b2, w = v[B.l_b1], v[B.l_b2], v[B.l_w]
        if b1 == b2 or w <= 1e-4: v[dl][bmap[b1]] = 1.0
        elif w >= 1 - 1e-4: v[dl][bmap[b2]] = 1.0
        else: v[dl][bmap[b1]] = 1.0 - w; v[dl][bmap[b2]] = w
    a = dict(verts=len(bm.verts), tris=0, quads=0, ngons=0, tri_count=0, open_edges=0, zero_faces=0,
             bones=[B.bone_names[b] for b in used_b], parts=len(B.parts))
    for f in bm.faces:
        k = len(f.verts)
        if k == 3: a['tris'] += 1
        elif k == 4: a['quads'] += 1
        else: a['ngons'] += 1
        a['tri_count'] += k - 2
        if f.calc_area() < 1e-10: a['zero_faces'] += 1
    a['open_edges'] = sum(1 for e in bm.edges if len(e.link_faces) != 2)
    for lay in (B.l_b1, B.l_b2): bm.verts.layers.int.remove(lay)
    bm.verts.layers.float.remove(B.l_w)
    bm.to_mesh(me); bm.free()
    ob['parts'] = json.dumps(B.parts)
    if arm_obj is not None:
        ob.parent = arm_obj
        md = ob.modifiers.new('Armature', 'ARMATURE'); md.object = arm_obj
    return ob, a


def part_trees(ob, dg=None):
    """One BVH per part of a module, in world space (posed if a depsgraph is given)."""
    names = json.loads(ob['parts'])
    src = ob.evaluated_get(dg) if dg else ob
    me = src.to_mesh() if dg else ob.data
    att = me.attributes['part'].data; mw = ob.matrix_world
    co = [mw @ v.co for v in me.vertices]; groups = {}
    for p in me.polygons: groups.setdefault(att[p.index].value, []).append(tuple(p.vertices))
    out = []
    for pid, polys in groups.items():
        idx = sorted({i for p in polys for i in p}); rm = {i: k for k, i in enumerate(idx)}
        vs = [co[i] for i in idx]
        lo = Vector((min(v.x for v in vs), min(v.y for v in vs), min(v.z for v in vs)))
        hi = Vector((max(v.x for v in vs), max(v.y for v in vs), max(v.z for v in vs)))
        out.append(dict(ob=ob.name, name=names[pid][0], bone=names[pid][1], lo=lo, hi=hi,
                        tree=BVHTree.FromPolygons(vs, [tuple(rm[i] for i in p) for p in polys])))
    if dg: src.to_mesh_clear()
    return out


def overlaps(parts_a, parts_b, allow=()):
    """Pairs of parts on different bones whose surfaces cross. allow: substrings of part names to skip."""
    res = []; same = parts_a is parts_b
    for i, a in enumerate(parts_a):
        for j, b in enumerate(parts_b):
            if (same and j <= i) or a['bone'] == b['bone']: continue
            if any(k in a['name'] or k in b['name'] for k in allow): continue
            if any(a['lo'][k] > b['hi'][k] or b['lo'][k] > a['hi'][k] for k in range(3)): continue
            n = len(a['tree'].overlap(b['tree']))
            if n: res.append((n, a['ob'], a['name'], a['bone'], b['ob'], b['name'], b['bone']))
    return sorted(res, reverse=True)