master
John Lauer ESC routing round: planes fixture (GND In1, +3V3 In2), Fable's grid router, offline DRC gate, live replay f564403 25d ago
#!/usr/bin/env python3
"""Board model, pad geometry and grid rasterisation shared by the ESC routing scripts.

Reads a KiCad 10 `.kicad_pcb` with the span-keeping tokenizer from
tools/make_placement_fixture.py and turns it into what a grid router needs:

  * every copper pad in BOARD coordinates (footprint position + rotation, KiCad's
    y-down convention: x' = x*cos + y*sin, y' = -x*sin + y*cos), with its outline as
    one or more polygons (rect/roundrect as the rectangle, oval as a stadium, circle as
    a 32-gon, custom pads as their anchor plus every gr_poly/gr_rect/gr_circle primitive)
  * the Edge.Cuts outline as one closed polygon (arcs sampled)
  * the net format (KiCad 10 boards carry `(net "NAME")` on pads and no net table)
  * a 0.1 mm grid with conservative rasterisers (a cell is copper if any part of the
    shape touches the cell square) and a disk dilation for clearance maps

Python 3 plus numpy. Nothing here writes a file.
"""
import math
import os
import sys

import numpy as np

_TOOLS = os.path.join(os.path.dirname(os.path.abspath(__file__)), "..", "..", "..", "tools")
sys.path.insert(0, os.path.abspath(_TOOLS))
from make_placement_fixture import parse, reference_of, num  # noqa: E402

RES = 0.1  # mm per grid cell
F_CU, B_CU = 0, 1
LAYER_NAMES = {F_CU: "F.Cu", B_CU: "B.Cu"}
LAYER_IDS = {"F.Cu": F_CU, "B.Cu": B_CU}


# ----------------------------------------------------------------------------- geometry

def rot_kicad(px, py, deg):
    """KiCad's board rotation of a local offset (y-down, positive angle = counterclockwise on screen)."""
    a = math.radians(deg)
    c, s = math.cos(a), math.sin(a)
    return px * c + py * s, -px * s + py * c


def circle_poly(cx, cy, r, n=32):
    return [(cx + r * math.cos(2 * math.pi * k / n), cy + r * math.sin(2 * math.pi * k / n)) for k in range(n)]


def stadium_poly(sx, sy, n=10):
    """Oval pad of size sx x sy centred at 0 (local frame)."""
    if sx >= sy:
        r = sy / 2.0
        d = (sx - sy) / 2.0
        pts = []
        for k in range(n + 1):
            a = -math.pi / 2 + math.pi * k / n
            pts.append((d + r * math.cos(a), r * math.sin(a)))
        for k in range(n + 1):
            a = math.pi / 2 + math.pi * k / n
            pts.append((-d + r * math.cos(a), r * math.sin(a)))
        return pts
    pts = stadium_poly(sy, sx, n)
    return [(y, x) for x, y in pts]


def arc_points(p0, pm, p1, max_step=0.35):
    """Sample a KiCad three-point arc (start, mid, end) into a polyline including both ends."""
    (x0, y0), (xm, ym), (x1, y1) = p0, pm, p1
    ax, ay = xm - x0, ym - y0
    bx, by = x1 - x0, y1 - y0
    d = 2.0 * (ax * by - ay * bx)
    if abs(d) < 1e-9:
        return [p0, p1]
    a2 = ax * ax + ay * ay
    b2 = bx * bx + by * by
    ux = (by * a2 - ay * b2) / d
    uy = (ax * b2 - bx * a2) / d
    cx, cy = x0 + ux, y0 + uy
    r = math.hypot(ux, uy)
    t0 = math.atan2(y0 - cy, x0 - cx)
    tm = math.atan2(ym - cy, xm - cx)
    t1 = math.atan2(y1 - cy, x1 - cx)

    def ccw(a, b):
        return (b - a) % (2 * math.pi)

    if ccw(t0, tm) <= ccw(t0, t1):
        sweep = ccw(t0, t1)
    else:
        sweep = -ccw(t1, t0)
    n = max(2, int(math.ceil(abs(sweep) * r / max_step)))
    return [(cx + r * math.cos(t0 + sweep * k / n), cy + r * math.sin(t0 + sweep * k / n)) for k in range(n + 1)]


def chain_polylines(pieces, tol=0.01):
    """Join open polylines end to end into one closed loop (the board outline)."""
    pieces = [list(p) for p in pieces if len(p) >= 2]
    if not pieces:
        return []
    loop = pieces.pop(0)
    while pieces:
        end = loop[-1]
        best = None
        for idx, p in enumerate(pieces):
            if math.dist(end, p[0]) < tol:
                best = (idx, False)
                break
            if math.dist(end, p[-1]) < tol:
                best = (idx, True)
                break
        if best is None:
            raise ValueError("Edge.Cuts outline is not one closed loop (gap at %s)" % (end,))
        idx, rev = best
        p = pieces.pop(idx)
        if rev:
            p.reverse()
        loop.extend(p[1:])
    if math.dist(loop[0], loop[-1]) < tol:
        loop.pop()
    return loop


def polygon_area(poly):
    a = 0.0
    for i in range(len(poly)):
        x0, y0 = poly[i]
        x1, y1 = poly[(i + 1) % len(poly)]
        a += x0 * y1 - x1 * y0
    return a / 2.0


# ----------------------------------------------------------------------------- pads

class Pad:
    __slots__ = ("ref", "name", "net", "x", "y", "angle", "shape", "size", "drill", "thru", "layers",
                 "polys", "key", "fp_x", "fp_y", "long_axis", "half_long", "half_short", "uuid", "polys_centre", "clearance")

    def __init__(self):
        self.polys = []

    def __repr__(self):
        return "Pad(%s %s at %.3f,%.3f)" % (self.key, self.net, self.x, self.y)

    @property
    def min_dim(self):
        return 2.0 * self.half_short

    def local_extent(self):
        """Extents of the pad copper in its own (rotated) frame: (half_long, half_short, axis unit vector)."""
        return self.half_long, self.half_short, self.long_axis


def _pad_local_polys(pad_node, shape, sx, sy):
    """Polygons of the pad shape in the pad-local unrotated frame, centred on the pad origin."""
    polys = []
    if shape in ("rect", "roundrect", "trapezoid", "chamfer"):
        polys.append([(-sx / 2, -sy / 2), (sx / 2, -sy / 2), (sx / 2, sy / 2), (-sx / 2, sy / 2)])
    elif shape == "oval":
        polys.append(stadium_poly(sx, sy))
    elif shape == "circle":
        polys.append(circle_poly(0.0, 0.0, sx / 2.0))
    elif shape == "custom":
        opts = pad_node.find("options")
        anchor = "rect"
        if opts is not None and opts.find("anchor") is not None:
            anchor = opts.find("anchor").value()
        if anchor == "circle":
            polys.append(circle_poly(0.0, 0.0, sx / 2.0))
        else:
            polys.append([(-sx / 2, -sy / 2), (sx / 2, -sy / 2), (sx / 2, sy / 2), (-sx / 2, sy / 2)])
        prims = pad_node.find("primitives")
        if prims is not None:
            for prim in prims.items:
                if not hasattr(prim, "tag"):
                    continue
                if prim.tag == "gr_poly":
                    pts = prim.find("pts")
                    if pts is not None:
                        poly = [(xy.number(1), xy.number(2)) for xy in pts.find_all("xy")]
                        w = prim.find("width")
                        wv = w.number(1) if w is not None else 0.0
                        if wv and wv > 0:
                            # stroked outline: widen by half the stroke (conservative bbox growth)
                            poly = _grow_poly(poly, wv / 2.0)
                        polys.append(poly)
                elif prim.tag == "gr_rect":
                    s, e = prim.find("start"), prim.find("end")
                    if s is not None and e is not None:
                        x0, y0, x1, y1 = s.number(1), s.number(2), e.number(1), e.number(2)
                        polys.append([(x0, y0), (x1, y0), (x1, y1), (x0, y1)])
                elif prim.tag == "gr_circle":
                    c, e = prim.find("center"), prim.find("end")
                    if c is not None and e is not None:
                        r = math.hypot(e.number(1) - c.number(1), e.number(2) - c.number(2))
                        w = prim.find("width")
                        r += (w.number(1) if w is not None else 0.0) / 2.0
                        polys.append(circle_poly(c.number(1), c.number(2), r))
                elif prim.tag == "gr_line":
                    s, e = prim.find("start"), prim.find("end")
                    w = prim.find("width")
                    if s is not None and e is not None:
                        hw = (w.number(1) if w is not None else 0.1) / 2.0
                        polys.append(_thick_line_poly((s.number(1), s.number(2)), (e.number(1), e.number(2)), hw))
    else:
        polys.append([(-sx / 2, -sy / 2), (sx / 2, -sy / 2), (sx / 2, sy / 2), (-sx / 2, sy / 2)])
    return polys


def _grow_poly(poly, d):
    """Cheap outward growth of a polygon: push every vertex away from the centroid by d."""
    cx = sum(p[0] for p in poly) / len(poly)
    cy = sum(p[1] for p in poly) / len(poly)
    out = []
    for x, y in poly:
        vx, vy = x - cx, y - cy
        n = math.hypot(vx, vy) or 1.0
        out.append((x + vx / n * d, y + vy / n * d))
    return out


def _thick_line_poly(a, b, hw):
    dx, dy = b[0] - a[0], b[1] - a[1]
    n = math.hypot(dx, dy) or 1.0
    nx, ny = -dy / n * hw, dx / n * hw
    return [(a[0] + nx, a[1] + ny), (b[0] + nx, b[1] + ny), (b[0] - nx, b[1] - ny), (a[0] - nx, a[1] - ny)]


def _oriented_extent(polys):
    xs = [p[0] for poly in polys for p in poly]
    ys = [p[1] for poly in polys for p in poly]
    return min(xs), min(ys), max(xs), max(ys)


# ----------------------------------------------------------------------------- board

class Board:
    def __init__(self, path):
        self.path = path
        with open(path, "r", encoding="utf-8", newline="") as f:
            self.text = f.read()
        self.root = parse(self.text)
        if self.root.tag != "kicad_pcb":
            raise ValueError("%s is not a kicad_pcb" % path)
        self.net_table = {}  # number -> name when the file has a (net N "NAME") table
        for n in self.root.find_all("net"):
            try:
                self.net_table[int(n.atom(1))] = n.atom(2) or ""
            except (TypeError, ValueError):
                pass
        self.net_format = "number" if self.net_table else "name"
        self.pads = []
        self.footprints = {}
        self._read_footprints()
        self.outline = self._read_outline()
        xs = [p[0] for p in self.outline]
        ys = [p[1] for p in self.outline]
        self.bbox = (min(xs), min(ys), max(xs), max(ys))
        self.nets = {}
        for p in self.pads:
            if p.net:
                self.nets.setdefault(p.net, []).append(p)
        self.copper_layers = [a.atom(1) for a in self.root.find("layers").items[1:] if a.atom(2) == "signal" or a.atom(2) == "power" or a.atom(2) == "mixed"]

    # -- footprints and pads
    def _read_footprints(self):
        for fp in self.root.find_all("footprint"):
            ref = reference_of(fp)
            at = fp.find("at")
            fx, fy = at.number(1), at.number(2)
            rot = at.number(3) or 0.0
            layer = fp.find("layer").value() if fp.find("layer") is not None else "F.Cu"
            fp_clr = fp.find("clearance")
            fp_clearance = fp_clr.number(1) if fp_clr is not None else 0.0
            self.footprints[ref] = {"x": fx, "y": fy, "rot": rot, "layer": layer, "lib": fp.atom(1), "node": fp}
            for pn in fp.find_all("pad"):
                name = pn.value()
                ptype = pn.atom(2)
                shape = pn.atom(3)
                layers_node = pn.find("layers")
                layers = [a.value for a in layers_node.items[1:]] if layers_node is not None else []
                copper = set()
                for l in layers:
                    if l == "*.Cu":
                        copper.update(("F.Cu", "B.Cu"))
                    elif l.endswith(".Cu"):
                        copper.add(l)
                if not copper or not name:
                    continue  # paste-only apertures and unnamed helper pads carry no copper
                pat = pn.find("at")
                px, py = pat.number(1), pat.number(2)
                pang = pat.number(3) if pat.atom(3) is not None else rot
                size = pn.find("size")
                sx, sy = size.number(1), size.number(2)
                dr = pn.find("drill")
                pad = Pad()
                pad.ref, pad.name, pad.key = ref, name, "%s.%s" % (ref, name)
                nn = pn.find("net")
                pad.net = nn.value() if nn is not None else ""
                if self.net_format == "number" and nn is not None:
                    try:
                        pad.net = self.net_table.get(int(nn.atom(1)), pad.net)
                    except (TypeError, ValueError):
                        pass
                ox, oy = rot_kicad(px, py, rot)
                pad.x, pad.y = fx + ox, fy + oy
                pad.angle = pang
                pad.shape, pad.size = shape, (sx, sy)
                pad.drill = dr.number(1) if dr is not None else None
                pad.thru = ptype == "thru_hole" or ptype == "np_thru_hole"
                pad.layers = frozenset(l for l in copper if l in ("F.Cu", "B.Cu"))
                pad.fp_x, pad.fp_y = fx, fy
                u = pn.find("uuid")
                pad.uuid = u.value() if u is not None else ""
                pclr = pn.find("clearance")
                # KiCad: the effective clearance between two items is the largest of the rule and both local overrides
                pad.clearance = max(fp_clearance, pclr.number(1) if pclr is not None else 0.0)
                local = _pad_local_polys(pn, shape, sx, sy)
                x0, y0, x1, y1 = _oriented_extent(local)
                lx, ly = (x1 - x0) / 2.0, (y1 - y0) / 2.0
                # centre of the copper in the local frame (custom pads are often offset from the origin)
                lcx, lcy = (x0 + x1) / 2.0, (y0 + y1) / 2.0
                if lx >= ly:
                    pad.half_long, pad.half_short = lx, ly
                    axis = rot_kicad(1.0, 0.0, pang)
                else:
                    pad.half_long, pad.half_short = ly, lx
                    axis = rot_kicad(0.0, 1.0, pang)
                pad.long_axis = axis
                for poly in local:
                    pad.polys.append([(pad.x + rot_kicad(x, y, pang)[0], pad.y + rot_kicad(x, y, pang)[1]) for x, y in poly])
                # keep the copper centre for escapes (custom pads)
                ccx, ccy = rot_kicad(lcx, lcy, pang)
                pad.polys_centre = (pad.x + ccx, pad.y + ccy)
                self.pads.append(pad)

    # -- outline
    def _read_outline(self):
        pieces = []
        for node in self.root.items:
            if not hasattr(node, "tag"):
                continue
            layer = node.find("layer") if node.tag and node.tag.startswith("gr_") else None
            if layer is None or layer.value() != "Edge.Cuts":
                continue
            if node.tag == "gr_line":
                s, e = node.find("start"), node.find("end")
                pieces.append([(s.number(1), s.number(2)), (e.number(1), e.number(2))])
            elif node.tag == "gr_arc":
                s, m, e = node.find("start"), node.find("mid"), node.find("end")
                pieces.append(arc_points((s.number(1), s.number(2)), (m.number(1), m.number(2)), (e.number(1), e.number(2))))
            elif node.tag == "gr_rect":
                s, e = node.find("start"), node.find("end")
                x0, y0, x1, y1 = s.number(1), s.number(2), e.number(1), e.number(2)
                return [(x0, y0), (x1, y0), (x1, y1), (x0, y1)]
            elif node.tag == "gr_circle":
                c, e = node.find("center"), node.find("end")
                r = math.hypot(e.number(1) - c.number(1), e.number(2) - c.number(2))
                return circle_poly(c.number(1), c.number(2), r, 64)
            elif node.tag == "gr_poly":
                pts = node.find("pts")
                return [(xy.number(1), xy.number(2)) for xy in pts.find_all("xy")]
        if not pieces:
            raise ValueError("no Edge.Cuts outline")
        return chain_polylines(pieces)

    def net_ref_sexpr(self, net_name):
        """The `(net ...)` line for copper on this board, matching the bridge's net_sexpr."""
        if self.net_format == "number":
            number = next((n for n, name in self.net_table.items() if name == net_name), None)
            if number is None:
                raise KeyError("net %r is not in the net table" % net_name)
            return "\t\t(net %d)\n" % number
        return "\t\t(net %s)\n" % _quote(net_name)

    def pad_by_key(self, key):
        for p in self.pads:
            if p.key == key:
                return p
        return None


def _quote(s):
    return '"' + s.replace("\\", "\\\\").replace('"', '\\"') + '"'


# ----------------------------------------------------------------------------- grid

class Grid:
    """A 0.1 mm cell grid covering the outline bbox. Cell (i, j) is centred on (x0 + i*RES, y0 + j*RES)."""

    def __init__(self, bbox, res=RES):
        self.res = res
        self.x0, self.y0 = bbox[0], bbox[1]
        self.W = int(round((bbox[2] - bbox[0]) / res)) + 1
        self.H = int(round((bbox[3] - bbox[1]) / res)) + 1

    def to_cell(self, x, y):
        # round half UP (not Python's half-to-even) so pads at a 0.5 mm pitch on .x5 coordinates
        # snap to cells a consistent 5 apart
        # 1e-3 cell (0.1 um) of tolerance: KiCad coordinates carry nanometre noise (123.849998)
        return int(math.floor((x - self.x0) / self.res + 0.5 + 1e-3)), int(math.floor((y - self.y0) / self.res + 0.5 + 1e-3))

    def to_xy(self, i, j):
        return round(self.x0 + i * self.res, 4), round(self.y0 + j * self.res, 4)

    def in_bounds(self, i, j):
        return 0 <= i < self.W and 0 <= j < self.H

    # -- rasterisers (all conservative: any touch of the cell square marks the cell)
    def poly_mask(self, poly, mask=None):
        if mask is None:
            mask = np.zeros((self.H, self.W), dtype=bool)
        xs = np.array([p[0] for p in poly])
        ys = np.array([p[1] for p in poly])
        i0 = max(0, int(math.floor((xs.min() - self.x0) / self.res)) - 1)
        i1 = min(self.W - 1, int(math.ceil((xs.max() - self.x0) / self.res)) + 1)
        j0 = max(0, int(math.floor((ys.min() - self.y0) / self.res)) - 1)
        j1 = min(self.H - 1, int(math.ceil((ys.max() - self.y0) / self.res)) + 1)
        if i1 < i0 or j1 < j0:
            return mask
        ci = np.arange(i0, i1 + 1)
        cj = np.arange(j0, j1 + 1)
        cx = self.x0 + ci * self.res
        cy = self.y0 + cj * self.res
        CX, CY = np.meshgrid(cx, cy)
        inside = np.zeros(CX.shape, dtype=bool)
        n = len(poly)
        for k in range(n):
            xa, ya = poly[k]
            xb, yb = poly[(k + 1) % n]
            if ya == yb:
                continue
            cond = (ya > CY) != (yb > CY)
            xint = xa + (CY - ya) * (xb - xa) / (yb - ya)
            inside ^= cond & (CX < xint)
        # boundary supercover: sample each edge finely and mark the containing cells
        for k in range(n):
            xa, ya = poly[k]
            xb, yb = poly[(k + 1) % n]
            L = math.hypot(xb - xa, yb - ya)
            steps = max(1, int(math.ceil(L / 0.025)))
            t = np.linspace(0.0, 1.0, steps + 1)
            ex = xa + (xb - xa) * t
            ey = ya + (yb - ya) * t
            ii = np.round((ex - self.x0) / self.res).astype(int)
            jj = np.round((ey - self.y0) / self.res).astype(int)
            ok = (ii >= i0) & (ii <= i1) & (jj >= j0) & (jj <= j1)
            inside[jj[ok] - j0, ii[ok] - i0] = True
        mask[j0:j1 + 1, i0:i1 + 1] |= inside
        return mask

    def segment_cells(self, a, b, half_width, slack=None):
        """Cells whose centre is within half_width + slack of segment a-b. Returns (jj, ii) index arrays."""
        if slack is None:
            slack = self.res / 2.0
        r = half_width + slack
        xa, ya = a
        xb, yb = b
        i0 = max(0, int(math.floor((min(xa, xb) - r - self.x0) / self.res)))
        i1 = min(self.W - 1, int(math.ceil((max(xa, xb) + r - self.x0) / self.res)))
        j0 = max(0, int(math.floor((min(ya, yb) - r - self.y0) / self.res)))
        j1 = min(self.H - 1, int(math.ceil((max(ya, yb) + r - self.y0) / self.res)))
        if i1 < i0 or j1 < j0:
            return np.array([], dtype=int), np.array([], dtype=int)
        ci = np.arange(i0, i1 + 1)
        cj = np.arange(j0, j1 + 1)
        CX, CY = np.meshgrid(self.x0 + ci * self.res, self.y0 + cj * self.res)
        dx, dy = xb - xa, yb - ya
        L2 = dx * dx + dy * dy
        if L2 < 1e-12:
            t = np.zeros_like(CX)
        else:
            t = np.clip(((CX - xa) * dx + (CY - ya) * dy) / L2, 0.0, 1.0)
        px = xa + t * dx
        py = ya + t * dy
        d = np.hypot(CX - px, CY - py)
        jj, ii = np.nonzero(d <= r + 1e-9)
        return jj + j0, ii + i0

    def disk_cells(self, c, radius, slack=None):
        return self.segment_cells(c, c, radius, slack)

    def dilate(self, mask, r_cells):
        """Binary dilation by a disk of radius r_cells (float, in cells)."""
        if r_cells <= 0:
            return mask.copy()
        R = int(math.floor(r_cells))
        H, W = mask.shape
        padded = np.zeros((H + 2 * R, W + 2 * R), dtype=bool)
        padded[R:R + H, R:R + W] = mask
        out = np.zeros_like(mask)
        r2 = r_cells * r_cells + 1e-9
        for dy in range(-R, R + 1):
            for dx in range(-R, R + 1):
                if dx * dx + dy * dy <= r2:
                    out |= padded[R + dy:R + dy + H, R + dx:R + dx + W]
        return out

    def outline_inside(self, outline, margin):
        """Cells whose centre is inside the outline polygon and at least `margin` from every edge."""
        inside, dmin = self.outline_distance(outline)
        return inside & (dmin >= margin)

    def outline_distance(self, outline):
        """(inside mask, distance of every cell centre to the nearest outline edge)."""
        inside = np.zeros((self.H, self.W), dtype=bool)
        ci = np.arange(self.W)
        cj = np.arange(self.H)
        CX, CY = np.meshgrid(self.x0 + ci * self.res, self.y0 + cj * self.res)
        n = len(outline)
        for k in range(n):
            xa, ya = outline[k]
            xb, yb = outline[(k + 1) % n]
            if ya == yb:
                continue
            cond = (ya > CY) != (yb > CY)
            xint = xa + (CY - ya) * (xb - xa) / (yb - ya)
            inside ^= cond & (CX < xint)
        dmin = np.full(CX.shape, np.inf)
        for k in range(n):
            xa, ya = outline[k]
            xb, yb = outline[(k + 1) % n]
            dx, dy = xb - xa, yb - ya
            L2 = dx * dx + dy * dy
            if L2 < 1e-12:
                continue
            t = np.clip(((CX - xa) * dx + (CY - ya) * dy) / L2, 0.0, 1.0)
            d = np.hypot(CX - (xa + t * dx), CY - (ya + t * dy))
            dmin = np.minimum(dmin, d)
        return inside, dmin


def fmt(v):
    return num(v)