#!/usr/bin/env python3 """ enrobe.py — Smoothly Enrobing Circles, single-file edition ============================================================ A self-contained implementation of implicit-field circle enrobing: given a set of circles and a list of which pairs should be bridged, computes a single smooth outline that wraps around all of them like a "coat" with rounded morphological-closing bridges between connected pairs. This file consolidates the enrobe/ package used throughout development into one script, with every bug found and fixed along the way folded in directly: 1. TAPERED TUBE RADIUS (fixes cap-drift / unevenness near a circle). Early version used tube_radius = max(Rc_i, Rc_j) for the whole bridge segment, which pulled the contour outward past the SMALLER circle's own coat boundary near the junction (off by tens of px). Switching to min(Rc_i, Rc_j) fixed the smaller circle but then under-shot near the LARGER circle's junction. The actual fix is to taper the tube radius linearly along the bridge segment, matching each circle's own coat radius exactly at its own end. See `build_field`, the tube-construction loop. 2. HALF-PLANE PHANTOM SELECTION (fixes inward-collapsing bridges for N>=3 circles). For exactly 2 circles, both possible bridge arcs are symmetric and both get drawn. For 3+, naively picking either phantom can carve the bridge into the middle of the cluster instead of around its outside. Fixed by checking, for each bridge, which side of the connecting line the *other* circles' centres sit on (averaged), and always choosing the phantom on the opposite side. See `build_field`, the `else` branch under `if n == 2`. 3. ARTEFACT-LOOP FILTERING (fixes spurious small contours showing up at full brightness in 3+ circle clusters). The tube fields from different bridges can cross near a cluster's centroid, producing a small extra zero-crossing loop distinct from the real outline. `contour_segs` can return multiple disconnected loops; only the largest is the real shape. `largest()` / `enrobe()` apply this filter by default; `preview()` now does too (it originally drew every loop at full glow intensity, which was misleading). 4. PINCH-OFF VALIDATION (catches a bad rho before wasting a full grid evaluation). A bridge with rho too small relative to the gap between circles will sever and read as just two separate hard circles touching the convex hull. `pinch_off_gap()` computes the safety margin analytically so this can be checked up front. 5. INPUT VALIDATION (clear errors instead of cryptic KeyErrors). `validate()` checks bridge indices are in range and not self-referential before any field math runs. 6. SVG EXPORT STEP SIZE (avoids absurdly long path strings). A contour from a fine grid (nx=1800-2400) has thousands of points; emitting one Q-command per point produces 30-40k character path strings. `step` subsampling (with an auto-scaled default in the CLI) keeps paths a manageable size without visible faceting, since the quadratic-Bezier midpoint spline still smooths between subsampled points. Usage as a library ------------------- from enrobe import Circle, Bridge, enrobe, enrobe_group_svg circles = [Circle(75, 400, 40), Circle(300, 400, 78)] bridges = [Bridge(0, 1, rho=420)] outline = enrobe(circles, bridges, m=18, s=16) svg_block = enrobe_group_svg("enrobe-kg", outline) Usage from the command line ----------------------------- python3 enrobe.py preview --circle 75,400,40 --circle 300,400,78 \\ --bridge 0,1,420 --m 18 --s 16 --out preview.png --check python3 enrobe.py svg --circle 75,400,40 --circle 300,400,78 \\ --bridge 0,1,420 --m 18 --s 16 --id enrobe-kg Geometry notes -------------- For unequal coats R_cA = R_A + m, R_cB = R_B + m, the phantom circle of radius rho that is externally tangent to both must satisfy: |P - c_A| = R_cA + rho |P - c_B| = R_cB + rho This is the intersection of two circles of radii (R_cA+rho) and (R_cB+rho) centred at c_A and c_B: d_A = R_cA + rho, d_B = R_cB + rho, D = |c_B - c_A| a = (d_A^2 - d_B^2 + D^2) / (2D) (foot along AB) h = sqrt(d_A^2 - a^2) (perpendicular offset) P = c_A + a*u_hat +/- h*n_hat where u_hat = (c_B - c_A)/D and n_hat = u_hat rotated 90 degrees. Pinch-off: the bridge severs when h = 0, i.e. when D >= d_A + d_B (circles too far apart for this rho), or D <= |d_A - d_B| (one coat would swallow the other) """ from __future__ import annotations import argparse import json import sys from dataclasses import dataclass import numpy as np import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt from matplotlib.path import Path from matplotlib.patches import PathPatch # ════════════════════════════════════════════════════════════════════ # 1. DATA TYPES # ════════════════════════════════════════════════════════════════════ @dataclass class Circle: """A circle: centre (x, y) and radius r, in scene units (px).""" x: float y: float r: float def __post_init__(self): if self.r <= 0: raise ValueError(f"Circle radius must be > 0, got {self.r}") @dataclass class Bridge: """ A connection between circles[a] and circles[b]. rho is the phantom-circle radius, the single free parameter controlling the bridge's look: large rho -> flat, taut connector small rho -> deep, curved-inward waist too small -> bridge pinches off entirely (see pinch_off_gap) """ a: int b: int rho: float def __post_init__(self): if self.rho <= 0: raise ValueError(f"Bridge rho must be > 0, got {self.rho}") def validate(circles: list[Circle], bridges: list[Bridge]) -> None: """Raise a clear error before any field math runs (BUG FIX #5).""" if not circles: raise ValueError("circles list is empty") n = len(circles) for i, br in enumerate(bridges): if not (0 <= br.a < n) or not (0 <= br.b < n): raise ValueError( f"bridge[{i}] references circle index {br.a} or {br.b}, " f"but only {n} circles were given (valid range 0..{n-1})" ) if br.a == br.b: raise ValueError(f"bridge[{i}] connects circle {br.a} to itself") def pinch_off_gap(c_a: Circle, c_b: Circle, rho: float, m: float) -> float: """ Safety margin (px) before this bridge pinches off (BUG FIX #4). Positive = safe. Negative = the bridge will sever / cusp. """ D = np.hypot(c_b.x - c_a.x, c_b.y - c_a.y) Rc_a, Rc_b = c_a.r + m, c_b.r + m d_a, d_b = Rc_a + rho, Rc_b + rho too_far_margin = (d_a + d_b) - D swallow_margin = D - abs(d_a - d_b) return min(too_far_margin, swallow_margin) # ════════════════════════════════════════════════════════════════════ # 2. SMOOTH MIN / MAX (quadratic, C1-continuous blend) # ════════════════════════════════════════════════════════════════════ def smin(a, b, k): if k <= 0: return np.minimum(a, b) h = np.maximum(1 - np.abs(a - b) / k, 0.0) return np.minimum(a, b) - 0.25 * k * h * h def smax(a, b, k): if k <= 0: return np.maximum(a, b) h = np.maximum(1 - np.abs(a - b) / k, 0.0) return np.maximum(a, b) + 0.25 * k * h * h # ════════════════════════════════════════════════════════════════════ # 3. DISTANCE PRIMITIVES # ════════════════════════════════════════════════════════════════════ def dist_pt(X, Y, cx, cy): return np.hypot(X - cx, Y - cy) def dist_seg(X, Y, a, b): """Vectorised point-to-segment distance. a, b are (x, y) tuples.""" ax, ay = a bx, by = b dx, dy = bx - ax, by - ay l2 = dx * dx + dy * dy if l2 == 0: return dist_pt(X, Y, ax, ay) t = np.clip(((X - ax) * dx + (Y - ay) * dy) / l2, 0, 1) return np.hypot(X - (ax + t * dx), Y - (ay + t * dy)) def phantom_centres(cA, cB, RcA, RcB, rho): """ Up to 2 phantom centres P with |P-cA| = RcA+rho and |P-cB| = RcB+rho (externally tangent to both coats). [] if geometry is infeasible. """ cA = np.array(cA, dtype=float) cB = np.array(cB, dtype=float) D = np.linalg.norm(cB - cA) if D < 1e-9: return [] d_a, d_b = RcA + rho, RcB + rho a = (d_a**2 - d_b**2 + D**2) / (2 * D) h2 = d_a**2 - a**2 if h2 < 0: return [] h = np.sqrt(h2) u = (cB - cA) / D n = np.array([-u[1], u[0]]) foot = cA + a * u return [tuple(foot + h * n), tuple(foot - h * n)] # ════════════════════════════════════════════════════════════════════ # 4. THE FIELD (this is where bug fixes #1 and #2 live) # ════════════════════════════════════════════════════════════════════ def build_field(X, Y, circles: list[Circle], bridges: list[Bridge], m: float, s: float): """ Scalar field whose zero level set is the enrobe contour. field < 0 inside the enrobed shape, > 0 outside. Construction: (a) hard union (np.minimum) of every circle's coat SDF (b) hard union with a TAPERED tube SDF per bridge — BUG FIX #1: radius interpolates linearly from Rc_i at circle i's end to Rc_j at circle j's end, so the tube matches each circle's own coat exactly at its own end instead of using a single max()/min() radius for the whole segment (which previously caused the contour to drift tens of px off the smaller or larger circle's true coat boundary near the junction). (c) smooth-max (smax) blend of a phantom-circle SDF per bridge, rounding the inward cusp where the bridge meets the coats. BUG FIX #2: for 3+ circles, the phantom side is chosen via a half-plane test against every other circle's centre, so the bridge always bulges outward (away from the cluster interior) rather than risking carving into the middle. """ validate(circles, bridges) centres = [(c.x, c.y) for c in circles] radii = [c.r for c in circles] n = len(circles) # (a) per-circle coat SDFs field = np.full(X.shape, 1e9) for (cx, cy), r in zip(centres, radii): field = np.minimum(field, dist_pt(X, Y, cx, cy) - (r + m)) # (b) tapered tube per bridge [BUG FIX #1] for br in bridges: i, j = br.a, br.b Rc_i, Rc_j = radii[i] + m, radii[j] + m ci, cj = centres[i], centres[j] dx, dy = cj[0] - ci[0], cj[1] - ci[1] L2 = dx * dx + dy * dy if L2 < 1e-9: continue t = np.clip(((X - ci[0]) * dx + (Y - ci[1]) * dy) / L2, 0.0, 1.0) tube_r = Rc_i * (1.0 - t) + Rc_j * t # <-- the taper proj_x, proj_y = ci[0] + t * dx, ci[1] + t * dy tube = np.hypot(X - proj_x, Y - proj_y) - tube_r field = np.minimum(field, tube) # (c) phantom carve per bridge [BUG FIX #2 in the n != 2 branch] for br in bridges: i, j, rho = br.a, br.b, br.rho cA, cB = centres[i], centres[j] RcA, RcB = radii[i] + m, radii[j] + m cands = phantom_centres(cA, cB, RcA, RcB, rho) if not cands: continue if n == 2: # symmetric: both arcs get carved (top + bottom bridge) for P in cands: field = smax(field, rho - dist_pt(X, Y, P[0], P[1]), s) else: # half-plane criterion: pick the phantom on the side # facing away from the rest of the cluster dx, dy = cB[0] - cA[0], cB[1] - cA[1] D = np.hypot(dx, dy) if D < 1e-9: continue nx, ny = -dy / D, dx / D mid = ((cA[0] + cB[0]) / 2, (cA[1] + cB[1]) / 2) sides = [ (centres[k][0] - mid[0]) * nx + (centres[k][1] - mid[1]) * ny for k in range(n) if k != i and k != j ] avg = sum(sides) / len(sides) if sides else 0.0 P = cands[1] if avg > 0 else cands[0] field = smax(field, rho - dist_pt(X, Y, P[0], P[1]), s) return field # ════════════════════════════════════════════════════════════════════ # 5. GRID + CONTOUR EXTRACTION (bug fix #3: artefact filtering) # ════════════════════════════════════════════════════════════════════ def make_grid(circles: list[Circle], pad: float = 40, nx: int = 1800): xmin = min(c.x - c.r for c in circles) - pad xmax = max(c.x + c.r for c in circles) + pad ymin = min(c.y - c.r for c in circles) - pad ymax = max(c.y + c.r for c in circles) + pad ny = max(2, int(nx * (ymax - ymin) / (xmax - xmin))) X, Y = np.meshgrid(np.linspace(xmin, xmax, nx), np.linspace(ymin, ymax, ny)) return X, Y def contour_segs(X, Y, Z, level: float = 0.0): """Marching-squares contour extraction -> list of (N,2) point arrays.""" fig0, ax0 = plt.subplots() cs = ax0.contour(X, Y, Z, levels=[level]) segs = [np.asarray(seg) for seg in cs.allsegs[0] if len(seg) > 2] plt.close(fig0) return segs def largest(segs: list[np.ndarray]): """Largest contour loop, discarding small artefact loops [BUG FIX #3].""" if not segs: return None return max(segs, key=len) def enrobe(circles: list[Circle], bridges: list[Bridge], m: float = 8, s: float = 13, pad: float = 40, nx: int = 1800): """Full pipeline: circles + bridges -> single outline polyline.""" X, Y = make_grid(circles, pad=pad, nx=nx) Z = build_field(X, Y, circles, bridges, m, s) segs = contour_segs(X, Y, Z) return largest(segs) # ════════════════════════════════════════════════════════════════════ # 6. SVG EXPORT (bug fix #6: auto-scaled step size) # ════════════════════════════════════════════════════════════════════ def path_from_polyline(pts: np.ndarray, step: int = 3, closed: bool = True) -> str: """ Quadratic-Bezier midpoint spline: smooth SVG path from a polyline. `step` subsamples points first; the spline still smooths between them, so this keeps paths short without visible faceting. """ pts = np.asarray(pts)[::step] n = len(pts) if n < 2: return "" mids = [ ((pts[i][0] + pts[(i + 1) % n][0]) / 2, (pts[i][1] + pts[(i + 1) % n][1]) / 2) for i in range(n) ] d = f"M {mids[-1][0]:.2f},{mids[-1][1]:.2f}" for i in range(n): d += f" Q {pts[i][0]:.2f},{pts[i][1]:.2f} {mids[i][0]:.2f},{mids[i][1]:.2f}" if closed: d += " Z" return d @dataclass class GlowStyle: color: str = "#ffffff" widths: tuple = (7.0, 3.5, 1.8) opacities: tuple = (0.04, 0.10, 0.80) def render(self, group_id: str, path_d: str) -> str: lines = [f'', f''] for w, o in zip(self.widths, self.opacities): lines.append(f' ') lines.append('') return "\n".join(lines) FOREGROUND = GlowStyle(widths=(7.0, 3.5, 1.8), opacities=(0.04, 0.10, 0.80)) BACKGROUND = GlowStyle(widths=(5.0, 2.5, 1.4), opacities=(0.035, 0.08, 0.65)) def enrobe_group_svg(group_id: str, pts: np.ndarray, style: GlowStyle = FOREGROUND, step: int = 3) -> str: path_d = path_from_polyline(pts, step=step) return style.render(group_id, path_d) def inject_before_svg_close(html: str, svg_snippet: str, count: int = 1) -> str: marker = "\n" if marker not in html: raise ValueError("Could not find '\\n' injection point.") return html.replace(marker, svg_snippet.rstrip() + "\n" + marker, count) # ════════════════════════════════════════════════════════════════════ # 7. RASTER PREVIEW (bug fix #3 applied here too: was drawing all # segments at full glow, including artefact loops) # ════════════════════════════════════════════════════════════════════ def _hero_bg(xmin, xmax, ymin, ymax, w=1100): h = int(w * (ymax - ymin) / (xmax - xmin)) yy = np.linspace(0, 1, h)[:, None, None] top = np.array([0.06, 0.05, 0.12]) bot = np.array([0.55, 0.18, 0.45]) base = top * (1 - yy) + bot * yy return np.clip(base[:, :, 0] * 0 + base, 0, 1) def _glow_line(ax, x, y, lw=1.4, glow=True): if glow: ax.plot(x, y, color="white", lw=lw * 5, alpha=0.04, solid_capstyle="round") ax.plot(x, y, color="white", lw=lw * 2.4, alpha=0.10, solid_capstyle="round") ax.plot(x, y, color="white", lw=lw, alpha=0.9, solid_capstyle="round") def preview(circles: list[Circle], bridges: list[Bridge], m: float = 8, s: float = 13, out: str = "preview.png", dpi: int = 190, lw: float = 1.4, glow: bool = True, pad: float = 60, nx: int = 1400, show_all_segments: bool = False) -> dict: """ Render a PNG preview. By default draws only the largest contour loop [BUG FIX #3] — pass show_all_segments=True to see artefacts. Returns {'segments': total found, 'drawn': how many rendered, 'path': out}. """ xs = [c.x for c in circles] ys = [c.y for c in circles] rs = [c.r for c in circles] xmin, xmax = min(xs[i] - rs[i] for i in range(len(circles))) - pad, \ max(xs[i] + rs[i] for i in range(len(circles))) + pad ymin, ymax = min(ys[i] - rs[i] for i in range(len(circles))) - pad, \ max(ys[i] + rs[i] for i in range(len(circles))) + pad X, Y = make_grid(circles, pad=pad, nx=nx) Z = build_field(X, Y, circles, bridges, m, s) all_segs = contour_segs(X, Y, Z) segs = all_segs if show_all_segments else ([max(all_segs, key=len)] if all_segs else []) fig, ax = plt.subplots(figsize=((xmax - xmin) / 200, (ymax - ymin) / 200)) ax.imshow(_hero_bg(xmin, xmax, ymin, ymax), extent=[xmin, xmax, ymax, ymin], aspect="auto", zorder=0) rng = np.random.default_rng(0) n_stars = 120 sx = rng.uniform(xmin, xmax, n_stars) sy = rng.uniform(ymin, ymax, n_stars) ax.scatter(sx, sy, s=rng.uniform(0.3, 4, n_stars), c="white", alpha=rng.uniform(0.10, 0.65, n_stars), zorder=1, linewidths=0) for c in circles: ax.add_patch(plt.Circle((c.x, c.y), c.r, fill=False, ec=(0.70, 0.66, 0.90), lw=0.7, alpha=0.4, zorder=2)) if segs: big = max(segs, key=len) ax.add_patch(PathPatch(Path(big), facecolor=(0.74, 0.69, 0.96), alpha=0.04, lw=0, zorder=2)) for seg in segs: _glow_line(ax, seg[:, 0], seg[:, 1], lw=lw, glow=glow) ax.set_xlim(xmin, xmax) ax.set_ylim(ymax, ymin) ax.axis("off") fig.savefig(out, dpi=dpi, bbox_inches="tight", pad_inches=0, facecolor=(0.06, 0.05, 0.12)) plt.close(fig) return {"segments": len(all_segs), "drawn": len(segs), "path": out} # ════════════════════════════════════════════════════════════════════ # 8. CLI # ════════════════════════════════════════════════════════════════════ def _parse_circles(args) -> list[Circle]: circles = [] for spec in (args.circle or []): try: x, y, r = (float(v) for v in spec.split(",")) except ValueError: sys.exit(f"--circle expects 'x,y,r', got: {spec!r}") circles.append(Circle(x, y, r)) if args.circles_json: try: raw = json.loads(args.circles_json) except json.JSONDecodeError as e: sys.exit(f"--circles-json is not valid JSON: {e}") for i, c in enumerate(raw): if not all(k in c for k in ("x", "y", "r")): sys.exit(f"--circles-json[{i}] missing one of x/y/r: {c}") circles.append(Circle(c["x"], c["y"], c["r"])) if not circles: sys.exit("No circles given. Use --circle x,y,r (repeatable) or --circles-json.\n" "Note: for negative coordinates use --circle=-10,20,30 (with '=').") return circles def _parse_bridges(args) -> list[Bridge]: bridges = [] for spec in (args.bridge or []): try: a, b, rho = spec.split(",") bridges.append(Bridge(int(a), int(b), float(rho))) except ValueError: sys.exit(f"--bridge expects 'a,b,rho' (a,b ints, rho float), got: {spec!r}") if args.bridges_json: try: raw = json.loads(args.bridges_json) except json.JSONDecodeError as e: sys.exit(f"--bridges-json is not valid JSON: {e}") for i, br in enumerate(raw): if not all(k in br for k in ("a", "b", "rho")): sys.exit(f"--bridges-json[{i}] missing one of a/b/rho: {br}") bridges.append(Bridge(br["a"], br["b"], br["rho"])) return bridges def _add_common_args(p: argparse.ArgumentParser): p.add_argument("--circle", action="append", metavar="x,y,r", help="one circle (repeatable). e.g. --circle 75,400,40 " "(use --circle=-10,20,30 for negative coords)") p.add_argument("--circles-json", metavar="JSON") p.add_argument("--bridge", action="append", metavar="a,b,rho", help="one bridge (repeatable). e.g. --bridge 0,1,420") p.add_argument("--bridges-json", metavar="JSON") p.add_argument("--m", type=float, default=8, help="coat margin (default: 8)") p.add_argument("--s", type=float, default=13, help="cusp smoothing width (default: 13)") p.add_argument("--pad", type=float, default=40) p.add_argument("--nx", type=int, default=1800) p.add_argument("--check", action="store_true", help="print pinch-off margin for each bridge before rendering") def _check_bridges(circles, bridges, m): for i, br in enumerate(bridges): margin = pinch_off_gap(circles[br.a], circles[br.b], br.rho, m) status = "OK" if margin > 5 else ("MARGINAL" if margin > 0 else "WILL PINCH OFF") print(f" bridge[{i}] ({br.a}->{br.b}, rho={br.rho}): " f"margin={margin:.1f}px [{status}]", file=sys.stderr) def cmd_preview(args): circles = _parse_circles(args) bridges = _parse_bridges(args) validate(circles, bridges) if args.check: _check_bridges(circles, bridges, args.m) result = preview(circles, bridges, m=args.m, s=args.s, out=args.out, nx=min(args.nx, 1400)) print(f"saved -> {result['path']} ({result['drawn']} of {result['segments']} " f"contour segment(s) drawn)") if result["segments"] > result["drawn"]: print(f" note: {result['segments'] - result['drawn']} small artefact loop(s) " f"excluded (normal for 3+ circle clusters)", file=sys.stderr) def cmd_svg(args): circles = _parse_circles(args) bridges = _parse_bridges(args) validate(circles, bridges) if args.check: _check_bridges(circles, bridges, args.m) outline = enrobe(circles, bridges, m=args.m, s=args.s, pad=args.pad, nx=args.nx) if outline is None: sys.exit("No contour found — check geometry (try --check)") style = BACKGROUND if args.background else FOREGROUND step = args.step if args.step is not None else max(1, len(outline) // 220) svg = enrobe_group_svg(args.id, outline, style=style, step=step) if len(svg) > 15000: print(f" note: path is {len(svg):,} chars (step={step}). " f"Pass --step {step*2} to shorten further.", file=sys.stderr) print(svg) def main(): p = argparse.ArgumentParser(prog="enrobe", description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) sub = p.add_subparsers(dest="command", required=True) p_preview = sub.add_parser("preview", help="render a PNG preview") _add_common_args(p_preview) p_preview.add_argument("--out", default="preview.png") p_preview.set_defaults(func=cmd_preview) p_svg = sub.add_parser("svg", help="print the SVG group") _add_common_args(p_svg) p_svg.add_argument("--id", default="enrobe-1") p_svg.add_argument("--step", type=int, default=None, help="default: auto-scaled from contour length") p_svg.add_argument("--background", action="store_true", help="use the toned-down background glow style") p_svg.set_defaults(func=cmd_svg) args = p.parse_args() args.func(args) if __name__ == "__main__": main()