"""Geometry helpers for the nutrition vat (exec'd by nv_build.py).

Everything is laid out in pixels of the reference picture's middle view (the pod's right side, seen square on):
u right (= toward the pod's rear, +y), v down; x across the pod in the same pixels (+x is the side the picture shows).
1 px = 4.62 mm (the head's top, v 70, stands 2.00 m over the floor, v 503). P() turns that into Blender metres
(x right, y to the rear, z up; the front faces -y).
Surfaces are quads: rounded polygons extruded with a bevel profile and capped by Coons patches, lathes, swept tubes and
cell solids. Nothing is a boolean. Parts are built with any winding; MB.to_mesh() makes each part's winding consistent
and turns it to face out.
"""
import bpy, bmesh, math
import numpy as np
from mathutils import Vector, Matrix

S = 0.00462                      # metres per picture pixel (middle view)
U0, V0 = 486.0, 503.0            # middle view: the body's middle (y 0) and the floor (z 0)
RMIN = 0.25 * S                  # smallest fillet left at a corner when an inset eats its radius


def X(px):
    return px * S


def Y(u):
    return (u - U0) * S


def Z(v):
    return (V0 - v) * S


def P(x, u, v):
    return Vector((x * S, (u - U0) * S, (V0 - v) * S))


def UV_OF(co):
    """middle-view pixels (u, v) of a world point"""
    return co.y / S + U0, V0 - co.z / S


# ------------------------------------------------------------------ frames: a plane (a, b) with a normal n
class Frame:
    def __init__(s, o=(0, 0, 0), a=(0, 1, 0), b=(0, 0, 1), n=(1, 0, 0)):
        s.o, s.a, s.b, s.n = Vector(o), Vector(a), Vector(b), Vector(n)

    def __call__(s, p, f):
        return s.o + s.a * p[0] + s.b * p[1] + s.n * f


SIDE = Frame()                                          # outline in (Y(u), Z(v)), depth along x
FRONT = Frame(a=(1, 0, 0), b=(0, 0, 1), n=(0, -1, 0))  # outline in (x, z), depth toward the viewer (-y)
TOP = Frame(a=(1, 0, 0), b=(0, 1, 0), n=(0, 0, 1))    # outline in (x, y), depth up


def uv2(u, v):
    """a side-outline point from middle-view pixels"""
    return (Y(u), Z(v))


# ------------------------------------------------------------------ rounded polygons with offset rings
class RPoly:
    """A simple polygon in 2D (metres) with a fillet at each corner.
    pts: corners; r: fillet radius per corner (0 = sharp); k: arc segments per corner (0 for a sharp corner; a corner
    between two collinear edges with different weights is a 'jog' and takes k 1); n: segments per edge i (corner i to
    i+1); w: bevel weight per edge (0 = the edge is not bevelled: the cap runs out to it; used where two slabs meet).
    ring(d) is the outline offset by d (outward +, inward -) with the same vertex count for every d."""

    def __init__(s, pts, r, k, n, w=None):
        s.p = [Vector((float(a), float(b))) for a, b in pts]
        N = len(s.p)
        s.r, s.k, s.n = list(r), list(k), list(n)
        s.w = list(w) if w is not None else [1.0] * N
        assert len(s.r) == len(s.k) == len(s.n) == len(s.w) == N
        area = sum(s.p[i].x * s.p[(i + 1) % N].y - s.p[(i + 1) % N].x * s.p[i].y for i in range(N))
        s.sgn = 1.0 if area > 0 else -1.0
        s.t, s.nrm, s.c = [], [], []
        for i in range(N):
            t = (s.p[(i + 1) % N] - s.p[i]).normalized()
            nn = Vector((t.y, -t.x)) * s.sgn
            s.t.append(t)
            s.nrm.append(nn)
            s.c.append(nn.dot(s.p[i]))
        s.kind, s.cen = [], []
        for i in range(N):
            a, b = (i - 1) % N, i
            cr = s.t[a].x * s.t[b].y - s.t[a].y * s.t[b].x
            if abs(cr) < 1e-9:
                s.kind.append('jog')
                assert s.k[i] == 1 and s.r[i] == 0, ("jog corner needs k 1, r 0", i)
                s.cen.append(None)
                continue
            convex = cr * s.sgn > 0
            if s.r[i] <= 0:
                s.kind.append('sharp')
                assert s.k[i] == 0, ("sharp corner needs k 0", i)
                s.cen.append(None)
                continue
            assert abs(s.w[a] - s.w[b]) < 1e-9, ("a rounded corner needs equal weights on both edges", i)
            sg = -1.0 if convex else 1.0
            s.kind.append('convex' if convex else 'concave')
            s.cen.append(s._solve(s.nrm[a], s.c[a] + sg * s.r[i], s.nrm[b], s.c[b] + sg * s.r[i]))
        s.start = []
        acc = 0
        for i in range(N):
            s.start.append(acc)
            acc += s.k[i] + 1 + s.n[i] - 1
        s.count = acc

    @staticmethod
    def _solve(n1, c1, n2, c2):
        det = n1.x * n2.y - n1.y * n2.x
        return Vector(((c1 * n2.y - c2 * n1.y) / det, (n1.x * c2 - n2.x * c1) / det))

    def corner_pts(s, i, d):
        N = len(s.p)
        a, b = (i - 1) % N, i
        ca, cb = s.c[a] + d * s.w[a], s.c[b] + d * s.w[b]
        kd = s.kind[i]
        if kd == 'sharp':
            return [s._solve(s.nrm[a], ca, s.nrm[b], cb)]
        if kd == 'jog':
            P0 = s.p[i]
            return [P0 + s.nrm[a] * (ca - s.nrm[a].dot(P0)), P0 + s.nrm[b] * (cb - s.nrm[b].dot(P0))]
        # the fillet of the offset outline: radius r +- d (kept above RMIN so the vertex count holds), its centre from the
        # offset lines
        if kd == 'convex':
            rad = max(s.r[i] + d * s.w[b], RMIN)
            C = s._solve(s.nrm[a], ca - rad, s.nrm[b], cb - rad)
            v0, v1 = s.nrm[a], s.nrm[b]
        else:
            rad = max(s.r[i] - d * s.w[b], RMIN)
            C = s._solve(s.nrm[a], ca + rad, s.nrm[b], cb + rad)
            v0, v1 = -s.nrm[a], -s.nrm[b]
        a0, a1 = math.atan2(v0.y, v0.x), math.atan2(v1.y, v1.x)
        da = a1 - a0
        while da > math.pi:
            da -= 2 * math.pi
        while da < -math.pi:
            da += 2 * math.pi
        k = s.k[i]
        return [C + Vector((math.cos(a0 + da * j / k), math.sin(a0 + da * j / k))) * rad for j in range(k + 1)]

    def ring(s, d=0.0):
        N = len(s.p)
        cps = [s.corner_pts(i, d) for i in range(N)]
        out = []
        for i in range(N):
            out += cps[i]
            A, B = cps[i][-1], cps[(i + 1) % N][0]
            m = s.n[i]
            out += [A.lerp(B, j / m) for j in range(1, m)]
        assert len(out) == s.count
        return out

    def mid(s, i):
        """loop index of corner i's middle point (an arc's middle, a sharp corner itself)"""
        return s.start[i] + s.k[i] // 2

    def edge_pt(s, i, m):
        """loop index of the m-th point along edge i (0 = the end of corner i's arc)"""
        return s.start[i] + s.k[i] + m

    def sides(s, corners):
        n = s.count
        return [((corners[(j + 1) % 4] - corners[j]) % n) for j in range(4)]


def rpoly_px(pts_uv, r_px, k, n, w=None, frame='side'):
    """a rounded polygon from picture pixels: side outlines take (u, v), front and top outlines take (x, z) / (x, y)
    already in pixels from the pod's axis"""
    if frame == 'side':
        pts = [uv2(u, v) for u, v in pts_uv]
    else:
        pts = [(a * S, b * S) for a, b in pts_uv]
    return RPoly(pts, [x * S for x in r_px], k, n, w)


# ------------------------------------------------------------------ mesh builder
class MB:
    def __init__(s):
        s.co, s.faces, s.mats, s.part = [], [], [], []
        s.parts = []                 # per part: hint (Vector or None)
        s.cur = -1
        s.smooth = 400               # relaxation sweeps on Coons caps (0: plain Coons)
        s.begin()

    def begin(s, hint=None):
        """start a new part: winding is made consistent per part and the part is turned to face out (hint: the
        direction its frontmost face looks; None: by volume for closed parts, by the centroid otherwise)"""
        s.parts.append(Vector(hint) if hint is not None else None)
        s.cur = len(s.parts) - 1

    def v(s, co):
        s.co.append(Vector(co))
        return len(s.co) - 1

    def vs(s, pts):
        return [s.v(p) for p in pts]

    def f(s, idx, mat):
        uq = []
        for i in idx:
            if i not in uq:
                uq.append(i)
        if len(uq) >= 3:
            s.faces.append(uq)
            s.mats.append(mat)
            s.part.append(s.cur)

    def bridge(s, la, lb, mat, closed=True):
        n = len(la)
        assert n == len(lb), (n, len(lb))
        for i in range(n if closed else n - 1):
            j = (i + 1) % n
            s.f([la[i], la[j], lb[j], lb[i]], mat(i) if callable(mat) else mat)

    def coons(s, loop, corners, mat):
        """fill a closed loop with a quad grid; corners: four positions in the loop (in loop order) that split it into
        sides of equal counts two by two"""
        n = len(loop)
        c = list(corners)
        side = lambda a, b: [loop[(a + k) % n] for k in range((b - a) % n + 1)]
        top, right = side(c[0], c[1]), side(c[1], c[2])
        bot, left = side(c[2], c[3])[::-1], side(c[3], c[0])[::-1]
        nu, nv = len(top) - 1, len(right) - 1
        assert len(bot) - 1 == nu and len(left) - 1 == nv, ("coons sides", nu, nv, len(bot) - 1, len(left) - 1)
        P_ = lambda i: s.co[i]

        def chord(pts):
            d = [0.0]
            for a, b in zip(pts, pts[1:]):
                d.append(d[-1] + (P_(b) - P_(a)).length)
            return [x / d[-1] for x in d]
        su = [(a + b) / 2 for a, b in zip(chord(top), chord(bot))]
        tv = [(a + b) / 2 for a, b in zip(chord(left), chord(right))]
        grid = [[None] * (nv + 1) for _ in range(nu + 1)]
        for i in range(nu + 1):
            grid[i][0], grid[i][nv] = top[i], bot[i]
        for j in range(nv + 1):
            grid[0][j], grid[nu][j] = left[j], right[j]
        T00, T10, T01, T11 = P_(top[0]), P_(top[nu]), P_(bot[0]), P_(bot[nu])
        pts = {}
        for i in range(1, nu):
            for j in range(1, nv):
                a, b = su[i], tv[j]
                pts[i, j] = ((1 - b) * P_(top[i]) + b * P_(bot[i]) + (1 - a) * P_(left[j]) + a * P_(right[j])
                             - ((1 - a) * (1 - b) * T00 + a * (1 - b) * T10 + (1 - a) * b * T01 + a * b * T11))
        if nu > 1 and nv > 1 and s.smooth:
            # relax the inside toward the average of its neighbours (a Tutte embedding: no folds on a convex outline)
            G = np.zeros((nu + 1, nv + 1, 3))
            for i in range(nu + 1):
                for j in range(nv + 1):
                    G[i, j] = pts[i, j] if (i, j) in pts else P_(grid[i][j])
            for _ in range(s.smooth):
                G[1:-1, 1:-1] = 0.25 * (G[:-2, 1:-1] + G[2:, 1:-1] + G[1:-1, :-2] + G[1:-1, 2:])
            for (i, j) in pts:
                pts[i, j] = Vector(G[i, j])
        for (i, j), p in pts.items():
            grid[i][j] = s.v(p)
        for i in range(nu):
            for j in range(nv):
                s.f([grid[i][j], grid[i + 1][j], grid[i + 1][j + 1], grid[i][j + 1]], mat)
        return grid

    def to_mesh(s, name, mat_names):
        me = bpy.data.meshes.get(name) or bpy.data.meshes.new(name)
        bm = bmesh.new()
        bv = [bm.verts.new(c) for c in s.co]
        bm.verts.ensure_lookup_table()
        idx = {m: i for i, m in enumerate(mat_names)}
        bad = 0
        byp = {}
        part_of = {}
        for f, m, pt in zip(s.faces, s.mats, s.part):
            try:
                face = bm.faces.new([bv[i] for i in f])
                face.material_index = idx[m]
                byp.setdefault(pt, []).append(face)
                part_of[face] = pt
            except ValueError:
                bad += 1
        for pt in list(byp):                    # points that coincide inside a part (a jog at full size) are one
            vs = list({v for f in byp[pt] for v in f.verts})
            bmesh.ops.remove_doubles(bm, verts=vs, dist=1e-8)
        byp = {}
        for f in bm.faces:
            byp.setdefault(part_of[f], []).append(f)
        for pt, faces in byp.items():
            bmesh.ops.recalc_face_normals(bm, faces=faces)
            hint = s.parts[pt]
            closed = all(len(e.link_faces) == 2 for f in faces for e in f.edges)
            if hint is not None:
                front = max(faces, key=lambda f: f.calc_center_median().dot(hint) + 1e-4 * f.normal.dot(hint))
                flip = front.normal.dot(hint) < 0
            elif closed:
                vol = sum(f.calc_center_median().dot(f.normal) * f.calc_area() for f in faces)
                flip = vol < 0
            else:
                cen = sum((f.calc_center_median() for f in faces), Vector()) / len(faces)
                tot = sum((f.calc_center_median() - cen).dot(f.normal) * f.calc_area() for f in faces)
                flip = tot < 0
            if flip:
                for f in faces:
                    f.normal_flip()
        bm.normal_update()
        me.clear_geometry() if len(me.vertices) else None
        me.materials.clear()                 # before the faces: clearing the slots also zeroes the faces' indices
        for m in mat_names:
            me.materials.append(bpy.data.materials.get(m) or bpy.data.materials.new(m))
        bm.to_mesh(me)
        bm.free()
        return me, bad


# ------------------------------------------------------------------ slabs: a rounded polygon extruded with a profile
def prof_full(h, b, nb=3, h_back=None, b_back=None):
    """a slab from a cap at +h round a bevel of radius b to the outline's full size, along to -h_back and round its
    bevel to a cap at -h_back (pixels; returns metres)"""
    hb = h if h_back is None else h_back
    bb = b if b_back is None else b_back
    out = []
    for j in range(nb + 1):
        ph = math.pi / 2 * j / nb
        out.append(((-b + b * math.sin(ph)) * S, (h - b + b * math.cos(ph)) * S))
    for j in range(nb + 1):
        ph = math.pi / 2 * (1 - j / nb)
        out.append(((-bb + bb * math.sin(ph)) * S, -(hb - bb + bb * math.cos(ph)) * S))
    return out


def prof_plate(x_out, x_in, b, nb=3):
    """a one-sided plate: a cap at x_out (px, either sign), bevel b, side faces back to x_in (no cap there)"""
    sg = 1.0 if x_out >= x_in else -1.0
    out = []
    for j in range(nb + 1):
        ph = math.pi / 2 * j / nb
        out.append(((-b + b * math.sin(ph)) * S, (x_out - sg * (b - b * math.cos(ph))) * S))
    out.append((0.0, x_in * S))
    return out


def slab(mb, poly, fr, prof, mats, caps=(True, True), corners=None, cap_mats=None, hint=None, skip=()):
    """rings of poly offset by d at depth f for (d, f) in prof, bridged; caps filled by a Coons patch on the first and
    the last ring. mats: one name, or one per ring gap. corners: the Coons corners (loop indices), default the middles
    of corners 0..3. skip: edges of poly whose faces are left out (where the slab sits on another one)."""
    mb.begin(hint)
    rings = [[mb.v(fr(p, f)) for p in poly.ring(d)] for d, f in prof]
    gone = set()
    for i in skip:
        a0 = poly.start[i] + poly.k[i]
        gone |= set(range(a0, a0 + poly.n[i]))
    n = len(rings[0])
    for k in range(len(rings) - 1):
        m = mats if isinstance(mats, str) else mats[k]
        for i in range(n):
            if i not in gone:
                j = (i + 1) % n
                mb.f([rings[k][i], rings[k][j], rings[k + 1][j], rings[k + 1][i]], m)
    cs = corners or [poly.mid(i) for i in range(4)]
    cm = cap_mats or (mats if isinstance(mats, str) else mats[0], mats if isinstance(mats, str) else mats[-1])
    if caps[0]:
        mb.coons(rings[0], cs, cm[0])
    if caps[1]:
        mb.coons(rings[-1], cs, cm[1])
    return rings


# ------------------------------------------------------------------ lathes
def lathe(mb, prof, o, axis, nseg, mats, cap_start=False, cap_end=False, ref=None, hint=None, phase=0.0, squash=1.0,
          closed=False, cap_mat=None):
    """prof: [(r, h)] in pixels along axis from o (a Vector in metres); mats: one name or one per segment.
    A point with r 0 is a single vertex. Caps close a ring with a Coons patch (nseg divisible by 4).
    closed: the profile is a loop (a ring section); its last segment runs back to the first point."""
    mb.begin(hint)
    if closed:
        if (Vector(prof[0]) - Vector(prof[-1])).length < 1e-9:
            prof = prof[:-1]
    n = Vector(axis).normalized()
    a = Vector(ref) if ref is not None else (Vector((0, 0, 1)) if abs(n.z) < 0.9 else Vector((1, 0, 0)))
    a = (a - n * a.dot(n)).normalized()
    b = n.cross(a)
    rings = []
    for (r, h) in prof:
        c = o + n * (h * S)
        if r <= 1e-9:
            rings.append([mb.v(c)])
        else:
            rings.append([mb.v(c + (a * math.cos(t) + b * math.sin(t) * squash) * (r * S))
                          for t in (phase + 2 * math.pi * k / nseg for k in range(nseg))])
    for k in range(len(rings) - 1):
        m = mats if isinstance(mats, str) else mats[k]
        ra, rb = rings[k], rings[k + 1]
        if len(ra) == 1 and len(rb) == 1:
            continue
        if len(ra) == 1 or len(rb) == 1:
            ring, tip = (rb, ra[0]) if len(ra) == 1 else (ra, rb[0])
            for i in range(nseg):
                mb.f([ring[i], ring[(i + 1) % nseg], tip], m)
        else:
            mb.bridge(ra, rb, m)
    if closed:
        mb.bridge(rings[-1], rings[0], mats if isinstance(mats, str) else mats[-1])
    q = nseg // 4
    if cap_start and len(rings[0]) > 1:
        mb.coons(rings[0], [0, q, 2 * q, 3 * q], cap_mat or (mats if isinstance(mats, str) else mats[0]))
    if cap_end and len(rings[-1]) > 1:
        mb.coons(rings[-1], [0, q, 2 * q, 3 * q], cap_mat or (mats if isinstance(mats, str) else mats[-1]))
    return rings


# ------------------------------------------------------------------ swept tubes
def fillet_path(pts, rad, segs=6):
    """a polyline (Vectors, metres) with each inner corner replaced by an arc of radius rad (metres; one value or one
    per inner corner)"""
    out = [pts[0]]
    for i in range(1, len(pts) - 1):
        r = rad if not isinstance(rad, (list, tuple)) else rad[i - 1]
        a, b, c = pts[i - 1], pts[i], pts[i + 1]
        d1, d2 = (a - b).normalized(), (c - b).normalized()
        ang = d1.angle(d2)
        if r <= 0 or ang > math.pi - 1e-3:
            out.append(b)
            continue
        t = r / math.tan(ang / 2)
        t = min(t, (a - b).length * 0.49, (c - b).length * 0.49)
        r = t * math.tan(ang / 2)
        p0, p1 = b + d1 * t, b + d2 * t
        cen = b + (d1 + d2).normalized() * (r / math.sin(ang / 2))
        v0, v1 = p0 - cen, p1 - cen
        R = v0.length
        ax = v0.cross(v1)
        if ax.length < 1e-12:
            out += [p0, p1]
            continue
        ax.normalize()
        sweep = v0.angle(v1)
        for j in range(segs + 1):
            m = Matrix.Rotation(sweep * j / segs, 3, ax)
            out.append(cen + (m @ v0).normalized() * R)
    out.append(pts[-1])
    return out


def tube(mb, path, radius, nseg, mats, cap_start=False, cap_end=False, up=None, hint=None, twist=0.0):
    """a tube swept along path (Vectors, metres); radius in metres (one value or one per path point); rings carried
    along by parallel transport"""
    mb.begin(hint)
    n = len(path)
    tans = []
    for i in range(n):
        t = (path[min(i + 1, n - 1)] - path[max(i - 1, 0)])
        tans.append(t.normalized())
    u = Vector(up) if up is not None else Vector((0, 0, 1))
    if abs(u.dot(tans[0])) > 0.95:
        u = Vector((1, 0, 0))
    u = (u - tans[0] * u.dot(tans[0])).normalized()
    rings = []
    for i in range(n):
        if i > 0:
            q = tans[i - 1].rotation_difference(tans[i])
            u = q @ u
            u = (u - tans[i] * u.dot(tans[i])).normalized()
        w = tans[i].cross(u)
        r = radius if not isinstance(radius, (list, tuple)) else radius[i]
        rings.append([mb.v(path[i] + (u * math.cos(2 * math.pi * k / nseg + twist) + w * math.sin(2 * math.pi * k / nseg + twist)) * r)
                      for k in range(nseg)])
    for k in range(n - 1):
        mb.bridge(rings[k], rings[k + 1], mats if isinstance(mats, str) else mats[k])
    q4 = nseg // 4
    m0 = mats if isinstance(mats, str) else mats[0]
    m1 = mats if isinstance(mats, str) else mats[-1]
    if cap_start:
        mb.coons(rings[0], [0, q4, 2 * q4, 3 * q4], m0)
    if cap_end:
        mb.coons(rings[-1], [0, q4, 2 * q4, 3 * q4], m1)
    return rings


# ------------------------------------------------------------------ cell solids (prisms on a grid, shared faces left out)
def cell_solid(mb, ni, nj, cells, F, hint=None):
    """cells: {(i, j): mat} solid cells on a grid of ni x nj cells; F(i, j, side) -> Vector for grid corner (i, j) on
    side 0 (back) or 1 (front). Faces: front and back of every cell, walls where a neighbour is missing."""
    mb.begin(hint)
    vid = {}

    def V_(i, j, sd):
        key = (i, j, sd)
        if key not in vid:
            vid[key] = mb.v(F(i, j, sd))
        return vid[key]
    for (i, j), mat in cells.items():
        mb.f([V_(i, j, 1), V_(i + 1, j, 1), V_(i + 1, j + 1, 1), V_(i, j + 1, 1)], mat)
        mb.f([V_(i, j, 0), V_(i, j + 1, 0), V_(i + 1, j + 1, 0), V_(i + 1, j, 0)], mat)
        for (di, dj), (a, b) in (((0, -1), ((i, j), (i + 1, j))), ((1, 0), ((i + 1, j), (i + 1, j + 1))),
                                 ((0, 1), ((i + 1, j + 1), (i, j + 1))), ((-1, 0), ((i, j + 1), (i, j)))):
            if (i + di, j + dj) not in cells:
                mb.f([V_(*a, 1), V_(*b, 1), V_(*b, 0), V_(*a, 0)], mat)


# ------------------------------------------------------------------ audit
def audit(ob):
    me = ob.data
    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)
    ng = sum(1 for f in bm.faces if len(f.verts) > 4)
    open_e = sum(1 for e in bm.edges if len(e.link_faces) == 1)
    nonman = sum(1 for e in bm.edges if len(e.link_faces) > 2)
    tiny = sum(1 for f in bm.faces if f.calc_area() < 1e-9)
    tris = sum(len(f.verts) - 2 for f in bm.faces)
    bm.free()
    return dict(quads=q, tris=t, ngons=ng, open_edges=open_e, nonmanifold_edges=nonman, zero_area=tiny, triangles=tris,
                verts=len(me.vertices))


# ------------------------------------------------------------------ decals (T_NV_Decals.png, 2048 x 1024)
# kind 'side': on the +x side, rect in middle-view pixels (u0, v0, u1, v1); kind 'front': facing -y at u, rect (x0, z0,
# x1, z1) in pixels with z up from the floor. atlas: (x0, y0, x1, y1) in the texture's pixels, y down. snap: the decal's
# grid is laid onto the surface under it (rays along -n); flat decals float at their plane.
DECALS = [
    dict(name="label", kind="side", x=55.0, rect=(526.0, 126.5, 562.0, 137.0), atlas=(0, 0, 1024, 297), grid=(12, 4),
         snap=True),
    dict(name="warn", kind="side", x=55.0, rect=(541.0, 342.0, 560.0, 353.0), atlas=(1024, 0, 1543, 300), grid=(6, 4),
         snap=True),
    dict(name="keys", kind="front", u=410.0 - 1.96, rect=(-35.0, 59.0, -14.0, 80.0), atlas=(1536, 0, 2048, 512),
         grid=(1, 1), snap=False),
]
DECAL_ATLAS = (2048, 1024)


def decal_frame(d):
    """(origin, a, b, n) in metres for a decal: a along the rect's first axis, b along its second (both as the
    picture reads: u to the rear / x to the right, and down for side decals, up for front decals)"""
    if d["kind"] == "side":
        return Vector((d["x"] * S, 0, 0)), Vector((0, 1, 0)), Vector((0, 0, -1)), Vector((1, 0, 0))
    return Vector((0, Y(d["u"]), 0)), Vector((1, 0, 0)), Vector((0, 0, 1)), Vector((0, -1, 0))


def decal_point(d, s, t):
    """world point at fraction s, t of a decal's rect (before snapping)"""
    o, a, b, n = decal_frame(d)
    r0, r1, r2, r3 = d["rect"]
    if d["kind"] == "side":
        return Vector((d["x"] * S, Y(r0 + (r2 - r0) * s), Z(r1 + (r3 - r1) * t)))
    return Vector(((r0 + (r2 - r0) * s) * S, Y(d["u"]), (r1 + (r3 - r1) * t) * S))


def decal_uv(d, co):
    """the atlas UV of a world point on decal d"""
    r0, r1, r2, r3 = d["rect"]
    if d["kind"] == "side":
        s, t = ((co.y / S + U0) - r0) / (r2 - r0), ((V0 - co.z / S) - r1) / (r3 - r1)
    else:
        s, t = (co.x / S - r0) / (r2 - r0), (r3 - co.z / S) / (r3 - r1)
    x0, y0, x1, y1 = d["atlas"]
    W_, H_ = DECAL_ATLAS
    return ((x0 + (x1 - x0) * s) / W_, 1.0 - (y0 + (y1 - y0) * t) / H_)


def which_decal(co):
    """the decal a world point lies on: near its plane and inside its rect"""
    best, bd = None, 1e9
    for d in DECALS:
        o, a, b, n = decal_frame(d)
        dist = abs((co - o).dot(n))
        r0, r1, r2, r3 = d["rect"]
        if d["kind"] == "side":
            s, t = ((co.y / S + U0) - r0) / (r2 - r0), ((V0 - co.z / S) - r1) / (r3 - r1)
        else:
            s, t = (co.x / S - r0) / (r2 - r0), (r3 - co.z / S) / (r3 - r1)
        if -0.05 <= s <= 1.05 and -0.05 <= t <= 1.05 and dist < bd:
            best, bd = d, dist
    return best