geom: port the pure-geometry half of sb-geom

Vectors, GEO and MEMBER records, member placement, sleeve and cavity
paths, exact polyline distance, corner-radius derivation, and the
monotone solver. Direct translation of legacy/openscad/lib/sb-geom.scad
at rev 8.0.0.

Angles stay in degrees, matching OpenSCAD, so every expression reads the
same as its source line. The 44 solver iterations, the 0.999 and 0.98
scale factors, the 0.05/179.95 cutoffs and the 1e9 sentinel are
reproduced exactly: they shaped the frozen oracle.

Region operations are not included -- they need a 2D boolean kernel and
follow with the Shapely layer.

49 unit tests, none of which touch the oracle. Harness proven by
mutation: radians for degrees, a shortened solver, a dropped scale
factor, a skipped crossing test and a flipped offset sign are each
caught. Oracle acceptance still skips; 236 unchanged.
This commit is contained in:
2026-08-19 02:58:35 -05:00
parent b67cc1290e
commit dfd02a4fd8
4 changed files with 1081 additions and 1 deletions
+49 -1
View File
@@ -1 +1,49 @@
"""Mechanical Compiler - geom."""
"""Geometry primitives ported from ``legacy/openscad/lib/sb-geom.scad``."""
from .primitives import ( # noqa: F401
SB_EPS,
SB_FACE_BOTH_OUT,
SB_FACE_MINUS_IN,
SB_FACE_PLUS_IN,
SB_SQRT3,
Path,
Point,
atan2_d,
cos_d,
sb_ccw,
sb_centroid,
sb_corner_radii,
sb_cross2,
sb_dist,
sb_line_isect,
sb_mid,
sb_path_gap,
sb_path_max_round,
sb_paths_cross,
sb_pt_path_dist,
sb_pt_seg_dist,
sb_region_min_gap,
sb_rot2,
sb_segs_cross,
sb_signed_area,
sb_solvable,
sb_solve,
sin_d,
tan_d,
vector_angle,
)
from .records import ( # noqa: F401
Geo,
Member,
cavity_path,
face_toward,
local_rect,
member_on_edge,
member_radial,
place,
sleeve_path,
sleeve_span,
sleeve_to_line,
strap_layer_paths,
strap_path,
)
+290
View File
@@ -0,0 +1,290 @@
"""
Port of ``legacy/openscad/lib/sb-geom.scad`` -- the pure-geometry half.
Everything here is a direct translation of the rev-8.0.0 reference. Nothing in
this module knows how many straps a profile has, so it is reused unchanged by
the 3x, 4x and any later N-strap arrangement.
Two conventions are carried over from the reference deliberately:
**Angles are degrees.** OpenSCAD's ``cos``/``sin``/``tan``/``atan2`` and BOSL2's
``vector_angle`` all work in degrees. Converting to radians at the boundary would
make each formula differ from its source line, so the degree helpers below are
used instead and every expression reads the same as the ``.scad``.
**The arbitrary constants are contract, not taste.** The 44 solver iterations,
the 0.999 and 0.98 scale factors, the 0.05/179.95 degree cutoffs and the 1e9
sentinel all shaped the frozen oracle. They are reproduced exactly.
Region operations -- area, connected parts, cleaning -- are not here. They need
a 2D boolean kernel and live in the Shapely-backed module.
"""
from __future__ import annotations
import math
from typing import Callable, Optional, Sequence, Tuple
Point = Tuple[float, float]
Path = Sequence[Point]
SB_SQRT3 = math.sqrt(3.0)
SB_EPS = 1e-7
# Face modes -- how a member's two broad faces are walled.
SB_FACE_PLUS_IN = 1 # local +Y faces an enclosed interior -> inside wall
SB_FACE_MINUS_IN = -1 # local -Y faces an enclosed interior -> inside wall
SB_FACE_BOTH_OUT = 0 # neither face encloses anything -> outside wall both
# ----------------------------------------------------------------------------
# Degree-based trigonometry
# ----------------------------------------------------------------------------
def cos_d(a: float) -> float:
return math.cos(math.radians(a))
def sin_d(a: float) -> float:
return math.sin(math.radians(a))
def tan_d(a: float) -> float:
return math.tan(math.radians(a))
def atan2_d(y: float, x: float) -> float:
return math.degrees(math.atan2(y, x))
def vector_angle(prev: Point, here: Point, nxt: Point) -> float:
"""
Angle at ``here`` between the segments to ``prev`` and ``nxt``, in degrees.
Matches BOSL2's three-point ``vector_angle``: always in [0, 180], never
signed. Degenerate zero-length arms return 0, which the callers' flat-vertex
cutoffs then treat as unroundable.
"""
ux, uy = prev[0] - here[0], prev[1] - here[1]
vx, vy = nxt[0] - here[0], nxt[1] - here[1]
nu = math.hypot(ux, uy)
nv = math.hypot(vx, vy)
if nu < SB_EPS or nv < SB_EPS:
return 0.0
c = (ux * vx + uy * vy) / (nu * nv)
return math.degrees(math.acos(max(-1.0, min(1.0, c))))
# ----------------------------------------------------------------------------
# Small vector helpers
# ----------------------------------------------------------------------------
def sb_rot2(p: Point, a: float) -> Point:
return (p[0] * cos_d(a) - p[1] * sin_d(a),
p[0] * sin_d(a) + p[1] * cos_d(a))
def sb_mid(a: Point, b: Point) -> Point:
return ((a[0] + b[0]) / 2.0, (a[1] + b[1]) / 2.0)
def sb_dist(a: Point, b: Point) -> float:
return math.hypot(b[0] - a[0], b[1] - a[1])
def sb_cross2(a: Point, b: Point) -> float:
return a[0] * b[1] - a[1] * b[0]
def sb_centroid(pts: Path) -> Point:
n = len(pts)
return (sum(p[0] for p in pts) / n, sum(p[1] for p in pts) / n)
def sb_signed_area(path: Path) -> float:
"""Signed area; positive means counter-clockwise."""
n = len(path)
total = 0.0
for i in range(n):
a = path[i]
b = path[(i + 1) % n]
total += a[0] * b[1] - b[0] * a[1]
return total / 2.0
def sb_ccw(path: Path) -> list[Point]:
return list(path) if sb_signed_area(path) >= 0 else list(reversed(path))
def sb_line_isect(p1: Point, d1: Point, p2: Point, d2: Point) -> Optional[Point]:
"""Intersection of line (p1,d1) with line (p2,d2). None if parallel."""
denom = sb_cross2(d1, d2)
if abs(denom) < SB_EPS:
return None
delta = (p2[0] - p1[0], p2[1] - p1[1])
t = sb_cross2(delta, d2) / denom
return (p1[0] + t * d1[0], p1[1] + t * d1[1])
# ----------------------------------------------------------------------------
# Measurement -- exact distance between two closed polylines
# ----------------------------------------------------------------------------
#
# The minimum distance between two disjoint polygons is always attained at a
# vertex of one of them, so sampling every vertex against every segment of the
# other (both ways round) is exact, not an approximation.
def sb_pt_seg_dist(p: Point, a: Point, b: Point) -> float:
abx, aby = b[0] - a[0], b[1] - a[1]
l2 = abx * abx + aby * aby
if l2 < SB_EPS:
return sb_dist(p, a)
t = ((p[0] - a[0]) * abx + (p[1] - a[1]) * aby) / l2
t = max(0.0, min(1.0, t))
return sb_dist(p, (a[0] + t * abx, a[1] + t * aby))
def sb_pt_path_dist(p: Point, path: Path) -> float:
n = len(path)
return min(sb_pt_seg_dist(p, path[i], path[(i + 1) % n]) for i in range(n))
def sb_segs_cross(a1: Point, a2: Point, b1: Point, b2: Point) -> bool:
"""Do two segments properly cross or touch?"""
d1 = (a2[0] - a1[0], a2[1] - a1[1])
d2 = (b2[0] - b1[0], b2[1] - b1[1])
den = sb_cross2(d1, d2)
if abs(den) < SB_EPS:
return False
w = (b1[0] - a1[0], b1[1] - a1[1])
t = sb_cross2(w, d2) / den
u = sb_cross2(w, d1) / den
return 0.0 <= t <= 1.0 and 0.0 <= u <= 1.0
def sb_paths_cross(p: Path, q: Path) -> bool:
lp, lq = len(p), len(q)
for i in range(lp):
for j in range(lq):
if sb_segs_cross(p[i], p[(i + 1) % lp], q[j], q[(j + 1) % lq]):
return True
return False
def sb_path_gap(p: Path, q: Path) -> float:
"""
Minimum distance between two closed paths.
Two polygons that CROSS may have no vertex near the other's boundary at all,
and the naive vertex test then reports a comfortable clearance across an
outright overlap -- exactly the false pass that lets a solver settle on a
degenerate arrangement. Crossing is therefore tested first and reported as
zero.
Nesting is deliberately not treated as overlap: a hole inside an outer
boundary is the normal case, and the distance between them is the wall
thickness this whole library exists to measure.
"""
if sb_paths_cross(p, q):
return 0.0
return min(min(sb_pt_path_dist(v, q) for v in p),
min(sb_pt_path_dist(v, p) for v in q))
def sb_region_min_gap(rgn: Sequence[Path]) -> float:
"""
Minimum distance between any two paths in a region.
For a finished cross-section this is the thinnest surviving piece of PLA+.
Fewer than two paths yields the reference's 1e9 sentinel rather than
infinity, because that value reaches the report unchanged.
"""
n = len(rgn)
if n < 2:
return 1e9
return min(sb_path_gap(rgn[i], rgn[j])
for i in range(n - 1) for j in range(i + 1, n))
# ----------------------------------------------------------------------------
# Corner rounding limits
# ----------------------------------------------------------------------------
def sb_corner_radii(path: Path, r: float) -> list[float]:
"""
Per-vertex roundover radius: ``r`` at reflex corners, 0 elsewhere, clamped
so the tangent points stay on their own edges.
Probing this by trial is not an option because BOSL2's rounding routine
raises a library-level error rather than returning a flag, so it is derived
up front.
"""
n = len(path)
cw = sb_signed_area(path) < 0
out: list[float] = []
for i in range(n):
prev = path[(i + n - 1) % n]
here = path[i]
nxt = path[(i + 1) % n]
turn = sb_cross2((here[0] - prev[0], here[1] - prev[1]),
(nxt[0] - here[0], nxt[1] - here[1]))
reflex = (turn > SB_EPS) if cw else (turn < -SB_EPS)
ang = vector_angle(prev, here, nxt)
if ang <= 0.05 or ang >= 179.95:
fits = 0.0
else:
fits = (0.98 * min(sb_dist(prev, here), sb_dist(here, nxt))
/ 2.0 * tan_d(ang / 2.0))
out.append(min(r, fits) if reflex else 0.0)
return out
def sb_path_max_round(path: Path) -> float:
"""Largest uniform corner radius the path can physically accept."""
n = len(path)
if n < 3:
return 0.0
terms: list[float] = []
for i in range(n):
prev = path[(i + n - 1) % n]
here = path[i]
nxt = path[(i + 1) % n]
l1 = sb_dist(prev, here)
l2 = sb_dist(here, nxt)
ang = vector_angle(prev, here, nxt)
if ang <= 0.05 or ang >= 179.95:
terms.append(1e9)
else:
terms.append(min(l1, l2) / 2.0 * tan_d(ang / 2.0))
return 0.999 * min(terms)
# ----------------------------------------------------------------------------
# Monotone solver
# ----------------------------------------------------------------------------
def sb_solve(f: Callable[[float], float], lo: float, hi: float,
target: float, iters: int = 44) -> float:
"""
Bisection for "place this member so that the resulting web is exactly W".
Rather than deriving a closed form per profile -- the source of most of the
wrong-by-a-cosine errors in earlier revisions -- the real measured quantity
is solved numerically. ``f`` must be non-decreasing on [lo, hi].
Fixed iteration count with no convergence test, matching the reference: the
number of evaluations is part of what produced the frozen values.
"""
while iters > 0:
mid = (lo + hi) / 2.0
if f(mid) < target:
lo = mid
else:
hi = mid
iters -= 1
return (lo + hi) / 2.0
def sb_solvable(f: Callable[[float], float], hi: float, target: float) -> bool:
"""True when f(hi) actually reaches the target, i.e. the solve is feasible."""
return f(hi) >= target
+308
View File
@@ -0,0 +1,308 @@
"""
The GEO and MEMBER records from ``sb-geom.scad``, and the cross-section paths
for a single member.
GEO
Everything about a single strap bundle and the PLA+ that wraps it. Built
once per build and threaded through every call.
MEMBER
One strap bundle's cross-section placement: centre, angle, and which broad
face (if either) looks into an enclosed interior.
Coordinate convention for a member
local +X = along the strap's WIDTH (the 15.875 mm direction)
local +Y = along the strap's THICKNESS (the 0.508 mm direction)
The member's angle rotates local +X onto the global direction given.
"normal+" is local +Y expressed globally.
The reference stores both records as bare indexed lists. They are dataclasses
here because the field names carry the meaning and Python has no reason to
reproduce a positional layout that existed to work around OpenSCAD's lack of
structures. Field order is preserved regardless, so the two can be compared.
"""
from __future__ import annotations
from dataclasses import dataclass
from typing import List, Optional
from .primitives import (
Path,
Point,
SB_FACE_BOTH_OUT,
SB_FACE_MINUS_IN,
SB_FACE_PLUS_IN,
atan2_d,
cos_d,
sb_ccw,
sb_line_isect,
sb_mid,
sin_d,
)
# ----------------------------------------------------------------------------
# GEO record
# ----------------------------------------------------------------------------
@dataclass(frozen=True)
class Geo:
"""A strap bundle and its surrounding PLA+."""
width: float # nominal strap width
strap_t: float # one strap's thickness
count: int # straps per bundle
clearance: float # fit clearance, applied to every cavity face
wall_inside: float # inside wall
wall_outside: float # outside wall
wall_edge: float # edge wall (caps the strap's narrow edges)
min_wall: float # minimum acceptable PLA thickness anywhere
@property
def bundle_t(self) -> float:
return self.count * self.strap_t
# Cavity = strap bundle grown by the fit clearance on all four faces.
@property
def cavity_w(self) -> float:
return self.width + 2.0 * self.clearance
@property
def cavity_t(self) -> float:
return self.bundle_t + 2.0 * self.clearance
def web_to_strap_gap(self, web: float) -> float:
"""
A declared "web" is the PLA+ that must survive between two neighbouring
strap CAVITIES. Because every cavity is inflated by the clearance, the
corresponding gap between the physical STRAPS is larger. Callers state
the web they want; this converts to the strap-to-strap spacing that
produces it, so a declared 1.2 mm web really is 1.2 mm of plastic.
"""
return web + 2.0 * self.clearance
# Distance from a member centreline out to each of its four sleeve faces.
def reach_plus(self, face: int) -> float:
return self.cavity_t / 2.0 + (self.wall_inside if face > 0
else self.wall_outside)
def reach_minus(self, face: int) -> float:
return self.cavity_t / 2.0 + (self.wall_inside if face < 0
else self.wall_outside)
@property
def reach_inside(self) -> float:
"""
Distance from centreline to the enclosed-interior side of the sleeve.
Only meaningful when the member actually has an interior face.
"""
return self.cavity_t / 2.0 + self.wall_inside
# ----------------------------------------------------------------------------
# MEMBER record
# ----------------------------------------------------------------------------
@dataclass(frozen=True)
class Member:
cx: float
cy: float
angle: float
face: int = SB_FACE_BOTH_OUT
@property
def centre(self) -> Point:
return (self.cx, self.cy)
@property
def axis(self) -> Point:
"""Unit vector along the strap's width."""
return (cos_d(self.angle), sin_d(self.angle))
@property
def normal(self) -> Point:
"""Local +Y expressed globally."""
return (-sin_d(self.angle), cos_d(self.angle))
@property
def inside_dir(self) -> Optional[Point]:
"""
Unit vector pointing from the member towards the profile interior.
None for SB_FACE_BOTH_OUT, which has no interior.
"""
if self.face == 0:
return None
nx, ny = self.normal
return (self.face * nx, self.face * ny)
def inside_wall_pt(self, g: Geo) -> Optional[Point]:
"""A point on the interior-facing surface of the member's sleeve."""
d = self.inside_dir
if d is None:
return None
return (self.cx + d[0] * g.reach_inside,
self.cy + d[1] * g.reach_inside)
def offset_out(self, d: float) -> "Member":
"""Translate along the outward normal (away from the interior)."""
nx, ny = self.normal
s = 1.0 if self.face == 0 else -self.face
return Member(self.cx + s * d * nx, self.cy + s * d * ny,
self.angle, self.face)
def face_toward(centre: Point, angle: float,
interior_target: Optional[Point]) -> int:
"""
Decide the face mode from a target point that lies inside the profile.
Pass ``interior_target=None`` for members with no enclosed side, which keeps
the member symmetric and stops the walls from becoming chiral.
"""
if interior_target is None:
return SB_FACE_BOTH_OUT
nx, ny = -sin_d(angle), cos_d(angle)
vx = interior_target[0] - centre[0]
vy = interior_target[1] - centre[1]
return SB_FACE_PLUS_IN if (nx * vx + ny * vy) >= 0 else SB_FACE_MINUS_IN
def member_on_edge(a: Point, b: Point, interior_target: Optional[Point],
shift: float = 0.0) -> Member:
"""Member lying on the segment a->b, optionally slid along its own axis."""
mid = sb_mid(a, b)
angle = atan2_d(b[1] - a[1], b[0] - a[0])
c = (mid[0] + shift * cos_d(angle), mid[1] + shift * sin_d(angle))
return Member(c[0], c[1], angle, face_toward(c, angle, interior_target))
def member_radial(r: float, angle: float, face: int = SB_FACE_BOTH_OUT) -> Member:
"""
Member placed radially: centre at distance r from origin along ``angle``,
with its width axis pointing outward. Used by spoke profiles.
"""
return Member(r * cos_d(angle), r * sin_d(angle), angle, face)
# ----------------------------------------------------------------------------
# Cross-section paths for one member
# ----------------------------------------------------------------------------
def place(m: Member, path: Path) -> List[Point]:
"""Place a locally-defined path into the member's frame."""
ca, sa = cos_d(m.angle), sin_d(m.angle)
return [(p[0] * ca - p[1] * sa + m.cx,
p[0] * sa + p[1] * ca + m.cy) for p in path]
def local_rect(half_w_lead: float, half_w_trail: float,
up: float, down: float) -> List[Point]:
"""
Rectangle in member-local coordinates.
half_w_lead : extent along +X (towards the member's leading end)
half_w_trail : extent along -X
up / down : extents along +Y / -Y
"""
return [(half_w_lead, -down),
(half_w_lead, up),
(-half_w_trail, up),
(-half_w_trail, -down)]
def strap_path(m: Member, g: Geo) -> List[Point]:
"""The physical strap bundle, as one rectangle."""
return place(m, local_rect(g.width / 2.0, g.width / 2.0,
g.bundle_t / 2.0, g.bundle_t / 2.0))
def strap_layer_paths(m: Member, g: Geo) -> List[List[Point]]:
"""Individual strap laminae, for display when count > 1."""
out = []
for i in range(g.count):
y = (i - (g.count - 1) / 2.0) * g.strap_t
t = g.strap_t / 2.0
rect = local_rect(g.width / 2.0, g.width / 2.0, t, t)
out.append(place(m, [(p[0], p[1] + y) for p in rect]))
return out
def cavity_path(m: Member, g: Geo) -> List[Point]:
"""The void the strap slides through."""
return place(m, local_rect(g.cavity_w / 2.0, g.cavity_w / 2.0,
g.cavity_t / 2.0, g.cavity_t / 2.0))
def sleeve_path(m: Member, g: Geo, ext_lead: float = 0.0,
ext_trail: float = 0.0) -> List[Point]:
"""
The PLA+ sleeve around one member.
``ext_lead`` / ``ext_trail`` extend the sleeve along its own axis beyond the
default edge wall. Junction construction uses this to make neighbouring
sleeves genuinely overlap instead of merely touching at a corner.
"""
half = g.cavity_w / 2.0 + g.wall_edge
f = m.face
return place(m, local_rect(half + ext_lead, half + ext_trail,
g.reach_plus(f), g.reach_minus(f)))
def sleeve_to_line(m: Member, g: Geo, line_pt: Point, line_dir: Point,
ext_lead: float = 0.0) -> List[Point]:
"""
Sleeve whose trailing end is cut by an arbitrary line rather than by a face
perpendicular to the axis.
This produces a butt joint flush against a neighbouring member's outer face,
which is how junctions are made structural rather than decorative. The
leading end stays perpendicular as usual. Falls back to a plain sleeve if
the line is parallel to the axis.
"""
c = m.centre
u = m.axis
n = m.normal
f = m.face
up = g.reach_plus(f)
dn = g.reach_minus(f)
half = g.cavity_w / 2.0 + g.wall_edge
p_up = (c[0] + n[0] * up, c[1] + n[1] * up)
p_dn = (c[0] - n[0] * dn, c[1] - n[1] * dn)
t_up = sb_line_isect(p_up, u, line_pt, line_dir)
t_dn = sb_line_isect(p_dn, u, line_pt, line_dir)
if t_up is None or t_dn is None:
return sleeve_path(m, g, ext_lead, 0.0)
reach = half + ext_lead
lead_up = (p_up[0] + u[0] * reach, p_up[1] + u[1] * reach)
lead_dn = (p_dn[0] + u[0] * reach, p_dn[1] + u[1] * reach)
return sb_ccw([lead_dn, lead_up, t_up, t_dn])
def sleeve_span(m: Member, g: Geo, pt_a: Point, dir_a: Point,
pt_b: Point, dir_b: Point) -> List[Point]:
"""
Sleeve cut by a line at BOTH ends.
A member that spans between two neighbours -- a gable crossbar, a chord
across a polygon -- butts flush against each of them instead of stopping
short or poking through.
"""
c = m.centre
u = m.axis
n = m.normal
f = m.face
p_up = (c[0] + n[0] * g.reach_plus(f), c[1] + n[1] * g.reach_plus(f))
p_dn = (c[0] - n[0] * g.reach_minus(f), c[1] - n[1] * g.reach_minus(f))
a_up = sb_line_isect(p_up, u, pt_a, dir_a)
a_dn = sb_line_isect(p_dn, u, pt_a, dir_a)
b_up = sb_line_isect(p_up, u, pt_b, dir_b)
b_dn = sb_line_isect(p_dn, u, pt_b, dir_b)
if a_up is None or a_dn is None or b_up is None or b_dn is None:
return sleeve_path(m, g)
return sb_ccw([a_dn, a_up, b_up, b_dn])