"""Geometry helpers for the convict body build (version 2).

Everything is made as a quad cage meant for a Subdivision Surface modifier: lofted shells with raised panels,
lathes, rounded boxes, hoses, and shaped armour plates that are laid onto a sculpted body through a radial map
of that body (MeshProf). The sculpted bodies themselves come from organic(): overlapping muscle volumes fused by a
voxel remesh, smoothed, and rebuilt as clean quads. No boolean modifiers are left in the result.
Every vertex is weighted to one bone (or blended between two) as it is made.
"""
import json
import bmesh
import bpy
import numpy as np
from bisect import bisect_right
from contextlib import contextmanager
from math import sin, cos, radians, degrees, pi, floor
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 S3(x, y, z): return Matrix.Diagonal((x, y, z, 1.0))


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: range across (degrees from the front on a shell, 0..1 on 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)


def lin(keys, x):
    """Piecewise linear through sorted (x, y) keys."""
    if x <= keys[0][0]: return keys[0][1]
    if x >= keys[-1][0]: return keys[-1][1]
    for (x0, y0), (x1, y1) in zip(keys, keys[1:]):
        if x <= x1: return y1 if x1 <= x0 else y0 + (y1 - y0) * (x - x0) / (x1 - x0)
    return keys[-1][1]


class Prof:
    """Analytic 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)
        rr = [(k[0] + k[1]) / 2 for k in self.v]; self.rref = max(sum(rr) / len(rr), 1e-4)

    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 surf(self, th, z, off=0.0):
        rx, ry, n, cx, cy = self.ev(z); r = sup(th, rx, ry, n) + off
        return Vector((cx + r * cos(th), cy + r * sin(th), z))

    def pt(self, phi, z, off=0.0): return self.surf(self.th(phi), z, off)

    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 MeshProf:
    """The same interface as Prof, but the radius is measured on a sculpted mesh: rays from the local Z axis of frame M
    out to the mesh surface, stored as a blurred (angle, z) map. Plates laid on it follow the body's shape."""

    def __init__(self, bvh, M, z0, z1, front=90.0, pdir=-1, nphi=120, dz=0.006, blur=1.3, axis=(0.0, 0.0), reach=0.6):
        self.front, self.pdir, self.z0, self.z1, self.ax = front, pdir, z0, z1, axis
        nz = max(2, int(round((z1 - z0) / dz)) + 1); self.nphi, self.nz = nphi, nz
        R = M.to_3x3(); A = np.zeros((nz, nphi), dtype=np.float64); ok = np.zeros((nz, nphi), dtype=bool)
        dirs = [R @ Vector((cos(2 * pi * i / nphi), sin(2 * pi * i / nphi), 0.0)) for i in range(nphi)]
        for j in range(nz):
            o = M @ Vector((axis[0], axis[1], z0 + (z1 - z0) * j / (nz - 1)))
            for i, d in enumerate(dirs):
                far = 0.0; start = o; travelled = 0.0
                for _ in range(8):                              # outermost surface along the ray
                    loc, nr, idx, dist = bvh.ray_cast(start, d, reach - travelled)
                    if loc is None: break
                    travelled += dist + 1e-4; far = travelled; start = loc + d * 1e-4
                    if travelled >= reach: break
                if far > 0: A[j, i] = far; ok[j, i] = True
        if ok.any():
            mean = float(A[ok].mean())
            for j in range(nz):                                 # holes: row average, else the nearest filled row
                if ok[j].any(): A[j, ~ok[j]] = A[j, ok[j]].mean()
            rows = [j for j in range(nz) if ok[j].any()]
            for j in range(nz):
                if not ok[j].any(): A[j, :] = A[min(rows, key=lambda r: abs(r - j)), :] if rows else mean
        if blur > 0:
            r = int(3 * blur) + 1; k = np.exp(-0.5 * (np.arange(-r, r + 1) / blur) ** 2); k /= k.sum()
            A = sum(w * np.roll(A, s, axis=1) for w, s in zip(k, range(-r, r + 1)))
            Ap = np.pad(A, ((r, r), (0, 0)), mode='edge'); A = sum(w * Ap[s:s + nz] for w, s in zip(k, range(2 * r + 1)))
        self.A = A; self.rref = max(float(A.mean()), 1e-4)

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

    def rho(self, th, z):
        fi = (th % (2 * pi)) / (2 * pi) * self.nphi; i0 = int(floor(fi)) % self.nphi; i1 = (i0 + 1) % self.nphi; a = fi - floor(fi)
        fj = min(max((z - self.z0) / (self.z1 - self.z0), 0.0), 1.0) * (self.nz - 1); j0 = int(floor(fj)); j1 = min(j0 + 1, self.nz - 1); b = fj - j0
        A = self.A
        return float((A[j0, i0] * (1 - a) + A[j0, i1] * a) * (1 - b) + (A[j1, i0] * (1 - a) + A[j1, i1] * a) * b)

    def surf(self, th, z, off=0.0):
        r = self.rho(th, z) + off
        return Vector((self.ax[0] + r * cos(th), self.ax[1] + r * sin(th), z))

    def pt(self, phi, z, off=0.0): return self.surf(self.th(phi), z, off)

    def nrm(self, phi, z):
        a = self.pt(phi + 1.5, z) - self.pt(phi - 1.5, z); b = self.pt(phi, z + 4e-3) - self.pt(phi, z - 4e-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 = []; self.crease = 0.0; self._creased = []     # crease: how hard the edges of what is being made stay under subdivision

    @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
        if self.crease: self._creased.append((f, self.crease))
        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 'blend'))
        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 add_mesh(self, me, name='base', uv_frame=None, uv_scale=0.08):
        """Copy a finished mesh (world coordinates, left side) into the module. Weights come from self.wfn.
        UVs are a cylinder unwrap around the local Z of uv_frame, continuous per face."""
        pid = len(self.parts); self.parts.append((name, 'base'))
        vs = [self._vert(v.co) for v in me.vertices]
        inv = uv_frame.inverted() if uv_frame is not None else Matrix.Identity(4)
        loc = [inv @ v.co for v in me.vertices]
        for p in me.polygons:
            idx = list(p.vertices); a0 = None; uvs = []
            for i in idx:
                q = loc[i]; a = np.arctan2(q.y, q.x)
                if a0 is None: a0 = a
                while a - a0 > pi: a -= 2 * pi
                while a - a0 < -pi: a += 2 * pi
                uvs.append((float(a) * uv_scale, q.z))
            self._face([vs[i] for i in idx], uvs, None, pid)
        return pid

    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 ell(self, c, r, rot=None, seg=16, rings=10, name='ell'):
        """Ellipsoid with half sizes r, optionally turned by rot (a 4x4)."""
        M = Matrix.Translation(Vector(c))
        if rot is not None: M = M @ rot
        with self.at(M @ S3(*r)):
            return self.lathe([(-cos(pi * k / rings), 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.0015, 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 _dedupe(xs, tol):
    out = []
    for x in sorted(xs):
        if not out or x - out[-1] > tol: out.append(x)
    return out


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 an analytic 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])
    rref = prof.rref
    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)
    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)
        p = prof.surf(prof.th(phi) + twist * t, z, off + (-thick if inner else h))
        if zwarp: p.z += zwarp(phi, t, u)
        return p
    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)
    return _solid(B, cols, zs, P3, H, matof, thick, name)


def _solid(B, cols, zs, P3, H, matof, thick, name, back=True):
    """Outer surface + edge wall (+ back) = one plate. back=False leaves the hidden back out: the wall then ends inside the body it sits on."""
    nc = len(cols)
    if not back:
        nz = len(zs)
        pid, vs = B.grid([[P3(u, z, H(u, z)) for u in cols] for z in zs], False,
                         cellmat=lambda j, i: matof((cols[i] + cols[i + 1]) / 2, (zs[j] + zs[j + 1]) / 2), name=name)
        ring = [(0, i) for i in range(nc)] + [(j, nc - 1) for j in range(1, nz)] + [(nz - 1, i) for i in range(nc - 2, -1, -1)] + [(j, 0) for j in range(nz - 2, 0, -1)]
        wall = [B._vert(P3(cols[i], zs[j], 0, True)) for j, i in ring]; n = len(ring)
        for k in range(n):
            (j0, i0), (j1, i1) = ring[k], ring[(k + 1) % n]
            B._face((vs[j0][i0], wall[k], wall[(k + 1) % n], vs[j1][i1]), ((0, 0), (0, thick), (0.01, thick), (0.01, 0)), None, pid)
        return pid, vs
    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 plate(B, prof, z0, z1, left, right, bottom=None, top=None, panels=(), seg=5, dz=0.035, thick=0.009, off=0.005, rim=0.003,
          eps=0.0025, sharp=0.0, name='plate', back=False):
    """An armour plate with a shaped outline, laid on prof (analytic or measured from a sculpted body).

    left/right: [(t, angle)] the two side edges, t = 0 at z0 .. 1 at z1, angles in degrees from the front.
    bottom/top: [(u, dz)] how far the z0 / z1 edge is pushed along z, u = 0 at the left edge .. 1 at the right.
    Kinks in those lists become corners of the outline (sharp > 0 adds holding loops beside them).
    panels work as on a shell, with the range across given as 0..1.
    """
    left = sorted(left); right = sorted(right); bottom = sorted(bottom or [(0, 0), (1, 0)]); top = sorted(top or [(0, 0), (1, 0)])
    Z0, Z1 = z0, z1; rref = prof.rref
    wavg = radians(max(sum(abs(lin(right, t) - lin(left, t)) for t in (0, .25, .5, .75, 1)) / 5, 1.0))

    def conv(e): return 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'] = (1 - a1, 1 - a0); pans.append(q)
    if rim > 0:
        pans.insert(0, dict(a=(None, None), z=(None, None), h=rim, mat=None, sym=False, eps=rim, _rim=True)); off -= rim
    must_u = {0.0, 1.0}; must_z = {Z0, Z1}; spec = []

    def addu(u):
        if 1e-5 < u < 1 - 1e-5: must_u.add(u)

    def addz(z):
        if Z0 + 1e-5 < z < Z1 - 1e-5: must_z.add(z)
    for u, _ in bottom + top:
        addu(u)
        if sharp: addu(u - conv(sharp)); addu(u + conv(sharp))
    for t, _ in left + right:
        z = Z0 + t * (Z1 - Z0); addz(z)
        if sharp: addz(z - sharp); addz(z + sharp)
    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)
        for a, sg in ((a0, -1), (a1, 1)):
            addu(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))
    mu = _dedupe(must_u, 1e-5); gap = 0.4 / seg
    us = sorted(mu + [b for b in (i / seg for i in range(seg + 1)) if all(abs(b - m) > gap for m in mu)])
    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.4 * dz for m in mz)])

    def inside(u, z, s, grow=0.0):
        a0, a1, pz0, pz1, h, m, eu, e = s
        if z < pz0 - grow * e - 1e-7 or z > pz1 + grow * e + 1e-7: return False
        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); L = lin(left, t); R = lin(right, t)
        zz = z + (1 - t) * lin(bottom, u) + t * lin(top, u)
        return prof.surf(prof.th(L + u * (R - L)), zz, off + (-thick if inner else h))
    cols = us if (prof.pdir > 0) == (lin(right, .5) > lin(left, .5)) else us[::-1]

    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
    return _solid(B, cols, zs, P3, H, matof, thick, name, back)


def organic(name, build_fn, voxel=0.006, smooth=10, lam=0.55, faces=1500):
    """Fuse overlapping volumes into one smooth body. build_fn(B) adds lofts and ellipsoids (world coordinates).
    Returns (bvh of the detailed surface, low quad mesh for the game)."""
    B = Builder(['root'], ['m']); build_fn(B); bm = B.bm
    bmesh.ops.recalc_face_normals(bm, faces=bm.faces[:])
    for lay in (B.l_b1, B.l_b2): bm.verts.layers.int.remove(lay)
    bm.verts.layers.float.remove(B.l_w)
    me = bpy.data.meshes.new(name + '_src'); bm.to_mesh(me); bm.free()
    ob = bpy.data.objects.new(name + '_src', me); bpy.context.scene.collection.objects.link(ob)
    md = ob.modifiers.new('remesh', 'REMESH'); md.mode = 'VOXEL'; md.voxel_size = voxel; md.adaptivity = 0.0
    sm = ob.modifiers.new('smooth', 'LAPLACIANSMOOTH'); sm.iterations = smooth; sm.lambda_factor = lam; sm.use_volume_preserve = True
    dg = bpy.context.evaluated_depsgraph_get(); hi = bpy.data.meshes.new_from_object(ob.evaluated_get(dg))
    bh = bmesh.new(); bh.from_mesh(hi)
    seen = set(); isl = []                                     # the remesh can leave shells inside the body: keep the outer one
    for f in bh.faces:
        if f.index in seen: continue
        comp = [f]; seen.add(f.index); k = 0
        while k < len(comp):
            for e in comp[k].edges:
                for g in e.link_faces:
                    if g.index not in seen: seen.add(g.index); comp.append(g)
            k += 1
        isl.append(comp)
    isl.sort(key=len, reverse=True)
    if len(isl) > 1: bmesh.ops.delete(bh, geom=[f for c in isl[1:] for f in c], context='FACES')
    for _ in range(3):                                         # the quad remesher refuses edges shorter than 0.1 mm
        short = [e for e in bh.edges if e.calc_length() < 4e-4]
        if not short: break
        bmesh.ops.collapse(bh, edges=short)
    bmesh.ops.dissolve_degenerate(bh, dist=1e-5, edges=bh.edges[:])
    bmesh.ops.recalc_face_normals(bh, faces=bh.faces[:])
    nme = sum(1 for e in bh.edges if not e.is_manifold); nmv = sum(1 for v in bh.verts if not v.is_manifold)
    flip = sum(1 for e in bh.edges if len(e.link_loops) == 2 and e.link_loops[0].vert == e.link_loops[1].vert)
    print(f'ORGANIC {name}: islands {len(isl)} {[len(c) for c in isl[:4]]}, non-manifold edges {nme} verts {nmv}, flipped pairs {flip}, faces {len(bh.faces)}')
    bh.to_mesh(hi); bvh = BVHTree.FromBMesh(bh); bh.free()
    ob.modifiers.clear(); ob.data = hi
    bpy.ops.object.select_all(action='DESELECT'); ob.select_set(True); bpy.context.view_layer.objects.active = ob
    ok = True
    try:
        r = bpy.ops.object.quadriflow_remesh(use_mesh_symmetry=False, use_preserve_sharp=False, use_preserve_boundary=False, smooth_normals=False,
                                             mode='FACES', target_faces=faces, seed=0)
        ok = 'FINISHED' in r and len(ob.data.polygons) > faces * 0.3 and all(len(p.vertices) == 4 for p in ob.data.polygons)
    except Exception as e:
        print('QUADRIFLOW EXC', name, e); ok = False
    if not ok:
        print('QUADRIFLOW FAILED', name, '-> decimate')
        ob.data = hi; d = ob.modifiers.new('dec', 'DECIMATE'); d.ratio = min(1.0, faces * 2.0 / max(len(hi.polygons) * 2, 1))
        ob.data = bpy.data.meshes.new_from_object(ob.evaluated_get(bpy.context.evaluated_depsgraph_get())); ob.modifiers.clear()
    low = ob.data; bpy.data.objects.remove(ob)
    print(f'ORGANIC {name}: detailed {len(hi.polygons)} faces -> {len(low.polygons)} faces, quads only: {ok}')
    return bvh, low


def finish(B, name, arm_obj, coll, materials, subsurf=True):
    """Turn the builder into a skinned mesh object. Returns (object, audit dict)."""
    bm = B.bm
    bmesh.ops.recalc_face_normals(bm, faces=bm.faces[:])
    if B._creased:
        cl = bm.edges.layers.float.get('crease_edge') or bm.edges.layers.float.new('crease_edge')
        for f, cv in B._creased:
            if f.is_valid:
                for e in f.edges: e[cl] = max(e[cl], cv)
    if not subsurf:
        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-11: a['zero_faces'] += 1
    a['open_edges'] = sum(1 for e in bm.edges if len(e.link_faces) == 1)      # plate walls end inside the body: those rims are meant to be open
    a['bad_edges'] = sum(1 for e in bm.edges if len(e.link_faces) not in (1, 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
    if subsurf:
        sd = ob.modifiers.new('Subdivision', 'SUBSURF'); sd.levels = 1; sd.render_levels = 2
    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). Sculpted base bodies are left out."""
    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():
        if names[pid][1] == 'base': continue
        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)