diff --git a/src/mechcomp/geom/__init__.py b/src/mechcomp/geom/__init__.py index f3e02e7..ec243eb 100644 --- a/src/mechcomp/geom/__init__.py +++ b/src/mechcomp/geom/__init__.py @@ -47,3 +47,16 @@ from .records import ( # noqa: F401 strap_layer_paths, strap_path, ) +from .rounding import ( # noqa: F401 + DEFAULT_FN, + EPSILON, + RoundoverTooLarge, + approx, + approx_pt, + arc, + deduplicate, + is_collinear, + path_merge_collinear, + round_corners, + segs, +) diff --git a/src/mechcomp/geom/rounding.py b/src/mechcomp/geom/rounding.py new file mode 100644 index 0000000..3452ddc --- /dev/null +++ b/src/mechcomp/geom/rounding.py @@ -0,0 +1,311 @@ +""" +The BOSL2 operations that ``sb-geom`` and ``sb-join`` call but Shapely has no +equivalent for: corner rounding, and the path cleanup that follows it. + +Ported from BOSL2 at commit ``92d697c2856de2fed93a33e858068589cefc2898``, which +is the commit the frozen oracle records. Read from source rather than from +recollection; the call chain is +``round_corners -> _circlecorner -> arc -> segs``. + +**Why this is transcribed rather than reimplemented.** A circular arc becomes a +finite number of straight segments, and the segment count determines the +enclosed area. ``SECTION_AREA_MM2`` is compared against the oracle at 1e-3 mm2, +and the extruded solid is these segments -- so the discretisation is not a +quality setting that a smoother modern approach could improve on. It is the +definition of the surface. Any adaptive subdivision or tolerance-based +flattening would produce a better curve and a failed build. + +The generators set ``$fn = facets`` with ``facets = 48``, and all 123 oracle +cases run at that default. With ``$fn`` positive, ``segs()`` ignores the radius +entirely, so segment counts depend only on swept angle. That is reproduced here +via the module-level default; it is a parameter rather than a constant because +the generators expose it as one. + +``round_corners`` **raises** when the requested roundovers do not fit the path, +matching BOSL2, which asserts rather than clamping. Callers are expected to +derive safe radii up front -- that is exactly what ``sb_corner_radii`` is for. +Silently clamping here would let a build succeed where the reference aborted. +""" + +from __future__ import annotations + +import math +from typing import List, Optional, Sequence, Tuple + +Point = Tuple[float, float] +Path = Sequence[Point] + +# BOSL2's private library epsilon (math.scad). +EPSILON = 1e-9 + +# OpenSCAD facet count set by the generators. Overridable, but every frozen +# oracle case was produced at this value. +DEFAULT_FN = 48 + + +class RoundoverTooLarge(ValueError): + """ + Raised where BOSL2 asserts "Roundovers are too big for the path." + + A distinct type because this is a caller error -- radii that were never + derived against the path -- and not a rejected profile. + """ + + +# ---------------------------------------------------------------------------- +# Comparison helpers (comparisons.scad) +# ---------------------------------------------------------------------------- + +def approx(a: float, b: float, eps: float = EPSILON) -> bool: + return abs(a - b) <= eps + + +def approx_pt(a: Point, b: Point, eps: float = EPSILON) -> bool: + """Componentwise, as BOSL2's list branch of ``approx``.""" + return abs(a[0] - b[0]) <= eps and abs(a[1] - b[1]) <= eps + + +def deduplicate(path: Path, closed: bool = False, + eps: float = EPSILON) -> List[Point]: + """ + Drop each point that equals the next one. + + Closed paths compare the last point against the first, which is how a + roundover that fully consumes a segment gets collapsed. + """ + pts = list(path) + n = len(pts) + if n == 0: + return [] + end = n if closed else n - 1 + return [pts[i] for i in range(n) + if i == end or not approx_pt(pts[i], pts[(i + 1) % n], eps)] + + +def _dist2line(d: Point, n: Point) -> float: + dot = d[0] * n[0] + d[1] * n[1] + return math.hypot(d[0] - dot * n[0], d[1] - dot * n[1]) + + +def is_collinear(pts: Sequence[Point], eps: float = EPSILON) -> bool: + """ + BOSL2's ``is_collinear`` via ``_noncollinear_triple``. + + The tolerance is *relative*: the furthest point from the first defines the + chord, and collinearity holds when every point lies within ``eps`` times + that chord length of the line. An absolute tolerance would behave + differently at the scales this library works at. + """ + if len(pts) < 3: + return True + pa = pts[0] + b = max(range(len(pts)), key=lambda i: math.hypot(pts[i][0] - pa[0], + pts[i][1] - pa[1])) + pb = pts[b] + nrm = math.hypot(pb[0] - pa[0], pb[1] - pa[1]) + if nrm <= eps: + return True + n = ((pb[0] - pa[0]) / nrm, (pb[1] - pa[1]) / nrm) + distlist = [_dist2line((p[0] - pa[0], p[1] - pa[1]), n) for p in pts] + return max(distlist) < eps * nrm + + +def path_merge_collinear(path: Path, closed: bool = True, + eps: float = EPSILON) -> List[Point]: + """ + Remove vertices that lie on the line between their neighbours. + + Exact butt joints and zero-radius fillets leave collinear vertices. They are + harmless in 2D but leave zero-area triangles the tessellator cannot resolve, + so a section that measures perfectly can still fail to extrude. + """ + pts = deduplicate(path, closed=closed, eps=eps) + n = len(pts) + if n <= 2: + return pts + out: List[Point] = [] + if not closed: + out.append(pts[0]) + rng = range(1, n - 1) + else: + rng = range(n) + for i in rng: + triple = (pts[(i - 1) % n], pts[i], pts[(i + 1) % n]) + if not is_collinear(triple, eps=eps): + out.append(triple[1]) + if not closed: + out.append(pts[-1]) + return out + + +# ---------------------------------------------------------------------------- +# Segmentation (utility.scad: segs) +# ---------------------------------------------------------------------------- + +def segs(r: float, angle: Optional[float] = None, + fn: int = DEFAULT_FN, fa: float = 12.0, fs: float = 2.0) -> int: + """ + Number of sides OpenSCAD gives a circle, or an arc of ``angle`` degrees. + + The ``2e-15`` subtraction is BOSL2's, guarding an angle that is fractionally + over its true value through rounding. It is reproduced because dropping it + can add a segment at exactly 360-divisible angles. + """ + if angle is not None: + return math.ceil(segs(r, None, fn, fa, fs) * abs(angle) / 360.0 - 2e-15) + if fn > 0: + return fn if fn > 3 else 3 + rr = r if math.isfinite(r) else 0.0 + return math.ceil(max(5.0, min(360.0 / fa, abs(rr) * 2.0 * math.pi / fs))) + + +# ---------------------------------------------------------------------------- +# Arc through two points about a centre (drawing.scad: arc) +# ---------------------------------------------------------------------------- + +def _vector_angle3(a: Point, b: Point, c: Point) -> float: + """Angle at ``b``, degrees, in [0, 180].""" + ux, uy = a[0] - b[0], a[1] - b[1] + vx, vy = c[0] - b[0], c[1] - b[1] + nu, nv = math.hypot(ux, uy), math.hypot(vx, vy) + if nu == 0.0 or nv == 0.0: + return 0.0 + cosv = (ux * vx + uy * vy) / (nu * nv) + return math.degrees(math.acos(max(-1.0, min(1.0, cosv)))) + + +def arc(n: int, cp: Point, points: Tuple[Point, Point]) -> List[Point]: + """ + ``n`` points along the arc from ``points[0]`` to ``points[1]`` about ``cp``. + + Sweep direction follows the sign of the 2D cross product, taking the short + way round -- BOSL2's ``long``/``cw``/``ccw`` flags are never passed by the + call sites this port needs, so the minor arc is always the one drawn. + """ + start, end = points + angle = _vector_angle3(start, cp, end) + v1 = (start[0] - cp[0], start[1] - cp[1]) + v2 = (end[0] - cp[0], end[1] - cp[1]) + prelim = v1[0] * v2[1] - v1[1] * v2[0] + direction = 1.0 if prelim > 0 else (-1.0 if prelim < 0 else 1.0) + r = math.hypot(v1[0], v1[1]) + final_angle = direction * angle + sa = math.degrees(math.atan2(v1[1], v1[0])) + out: List[Point] = [] + for i in range(n): + theta = sa + i * final_angle / (n - 1) + out.append((r * math.cos(math.radians(theta)) + cp[0], + r * math.sin(math.radians(theta)) + cp[1])) + return out + + +def _circlecorner(points: Tuple[Point, Point, Point], + parm: Tuple[float, float], fn: int = DEFAULT_FN) -> List[Point]: + """ + One rounded corner: ``parm`` is ``(d, r)``, the tangent setback and radius. + + A straight vertex -- half-angle 90 degrees -- degenerates to the two tangent + points with no arc between them. + """ + prev_p, here, next_p = points + angle = _vector_angle3(prev_p, here, next_p) / 2.0 + d, r = parm + + pux, puy = prev_p[0] - here[0], prev_p[1] - here[1] + nux, nuy = next_p[0] - here[0], next_p[1] - here[1] + pn = math.hypot(pux, puy) + nn = math.hypot(nux, nuy) + prev_u = (pux / pn, puy / pn) + next_u = (nux / nn, nuy / nn) + + start = (here[0] + prev_u[0] * d, here[1] + prev_u[1] * d) + end = (here[0] + next_u[0] * d, here[1] + next_u[1] * d) + + if approx(angle, 90.0): + return [start, end] + + bx, by = prev_u[0] + next_u[0], prev_u[1] + next_u[1] + bn = math.hypot(bx, by) + scale = r / math.sin(math.radians(angle)) + center = (scale * bx / bn + here[0], scale * by / bn + here[1]) + + n = max(3, math.ceil((90.0 - angle) / 180.0 * segs(r, None, fn))) + return arc(n, center, (start, end)) + + +# ---------------------------------------------------------------------------- +# round_corners (rounding.scad), method="circle", measure="radius" +# ---------------------------------------------------------------------------- + +def round_corners(path: Path, radius, closed: bool = True, + fn: int = DEFAULT_FN) -> List[Point]: + """ + Round each corner of ``path`` to its own radius. + + ``radius`` is a scalar or one value per vertex. Zero leaves a vertex + untouched, which is how ``sb_fillet_concave`` rounds only reflex corners + while keeping every convex corner bit-exact. + + Raises ``RoundoverTooLarge`` when the setbacks overrun an edge, matching + BOSL2's assertion. The message carries the same scale factors BOSL2 reports, + since those say directly how much too large the request was. + """ + pts = list(path) + n = len(pts) + if n < 3: + raise ValueError("Path has length %d. Length must be 3 or more." % n) + + parm = [float(radius)] * n if isinstance(radius, (int, float)) \ + else [float(x) for x in radius] + if len(parm) != n: + raise ValueError("radius list length %d does not match path length %d" + % (len(parm), n)) + + dk: List[Tuple[float, ...]] = [] + for i in range(n): + bit = (pts[(i - 1) % n], pts[i], pts[(i + 1) % n]) + degenerate = approx_pt(bit[0], bit[1]) or approx_pt(bit[1], bit[2]) + angle = None if degenerate else _vector_angle3(*bit) / 2.0 + + if not closed and (i == 0 or i == n - 1): + dk.append((0.0,)) + continue + if parm[i] == 0: + dk.append((0.0,)) + continue + if angle is None: + raise ValueError("Repeated point in path at index %d with nonzero " + "rounding" % i) + if approx(angle, 0.0): + raise ValueError("Path turns back on itself at index %d with " + "nonzero rounding" % i) + dk.append((parm[i] / math.tan(math.radians(angle)), parm[i])) + + lengths = [math.hypot(pts[i % n][0] - pts[(i - 1) % n][0], + pts[i % n][1] - pts[(i - 1) % n][1]) + for i in range(n + 1)] + + scalefactors: List[float] = [] + for i in range(n): + if not (closed or (i != 0 and i != n - 1)): + continue + back = dk[(i - 1) % n][0] + dk[i][0] + fwd = dk[i][0] + dk[(i + 1) % n][0] + scalefactors.append(min( + math.inf if back == 0 else lengths[i] / back, + math.inf if fwd == 0 else lengths[i + 1] / fwd, + )) + + if scalefactors and min(scalefactors) < 1.0: + raise RoundoverTooLarge( + "Roundovers are too big for the path. If you multiply them by this " + "vector they should fit: %r" % (scalefactors,)) + + out: List[Point] = [] + for i in range(n): + if dk[i][0] == 0: + out.append(pts[i]) + else: + bit = (pts[(i - 1) % n], pts[i], pts[(i + 1) % n]) + out.extend(_circlecorner(bit, (dk[i][0], dk[i][1]), fn)) + return deduplicate(out, closed=False) diff --git a/tests/test_rounding.py b/tests/test_rounding.py new file mode 100644 index 0000000..dc0e782 --- /dev/null +++ b/tests/test_rounding.py @@ -0,0 +1,300 @@ +""" +Unit tests for the BOSL2 rounding port. + +The expected values here come from the pinned BOSL2 source read directly, not +from this implementation. Where a number looks arbitrary -- 12 points on a right +angle, a relative collinearity tolerance -- it is arbitrary, and that is the +point: it is what produced the frozen oracle, and a tidier choice would be a +different surface. + +Geometric sanity is checked separately from segmentation, so a failure says +which of the two broke. +""" + +from __future__ import annotations + +import math + +import pytest + +from mechcomp.geom.rounding import ( + DEFAULT_FN, + EPSILON, + RoundoverTooLarge, + approx, + approx_pt, + arc, + deduplicate, + is_collinear, + path_merge_collinear, + round_corners, + segs, +) + +SQUARE10 = [(0.0, 0.0), (10.0, 0.0), (10.0, 10.0), (0.0, 10.0)] + + +# ---------------------------------------------------------------------------- +# segs -- $fn overrides everything +# ---------------------------------------------------------------------------- + +def test_segs_ignores_radius_when_fn_is_set(): + """With $fn positive the adaptive $fa/$fs path is never taken.""" + assert segs(1.0) == DEFAULT_FN + assert segs(1000.0) == DEFAULT_FN + assert segs(0.001) == DEFAULT_FN + + +def test_segs_floor_of_three(): + assert segs(5.0, None, 2) == 3 + assert segs(5.0, None, 12) == 12 + + +def test_segs_for_an_arc_scales_with_swept_angle(): + assert segs(5.0, 360.0) == 48 + assert segs(5.0, 180.0) == 24 + assert segs(5.0, 90.0) == 12 + + +def test_segs_epsilon_guard_prevents_an_extra_segment(): + """ + The 2e-15 subtraction stops a fractionally-over angle rounding up. Without + it, an angle a hair above an exact divisor gains a whole segment. + """ + assert segs(5.0, 90.0 + 1e-14) == 12 + + +def test_segs_adaptive_path_when_fn_is_zero(): + assert segs(10.0, None, 0, 12.0, 2.0) == 30 + + +# ---------------------------------------------------------------------------- +# arc +# ---------------------------------------------------------------------------- + +def test_arc_returns_exactly_n_points_on_the_circle(): + cp = (0.0, 0.0) + pts = arc(12, cp, ((1.0, 0.0), (0.0, 1.0))) + assert len(pts) == 12 + for p in pts: + assert math.hypot(*p) == pytest.approx(1.0) + + +def test_arc_endpoints_are_the_requested_points(): + pts = arc(7, (0.0, 0.0), ((1.0, 0.0), (0.0, 1.0))) + assert pts[0] == pytest.approx((1.0, 0.0)) + assert pts[-1] == pytest.approx((0.0, 1.0), abs=1e-12) + + +def test_arc_takes_the_short_way_round(): + """A quarter turn, not the three-quarter turn the other way.""" + pts = arc(5, (0.0, 0.0), ((1.0, 0.0), (0.0, 1.0))) + assert all(p[0] >= -1e-12 and p[1] >= -1e-12 for p in pts) + + +def test_arc_direction_follows_the_cross_product_sign(): + cw = arc(5, (0.0, 0.0), ((1.0, 0.0), (0.0, -1.0))) + assert cw[-1] == pytest.approx((0.0, -1.0), abs=1e-12) + assert all(p[1] <= 1e-12 for p in cw) + + +# ---------------------------------------------------------------------------- +# Comparison and cleanup helpers +# ---------------------------------------------------------------------------- + +def test_approx_uses_the_library_epsilon(): + assert approx(1.0, 1.0 + EPSILON / 2) + assert not approx(1.0, 1.0 + EPSILON * 10) + + +def test_deduplicate_open_keeps_the_last_point(): + pts = [(0.0, 0.0), (0.0, 0.0), (1.0, 0.0), (1.0, 0.0)] + assert deduplicate(pts, closed=False) == [(0.0, 0.0), (1.0, 0.0)] + + +def test_deduplicate_closed_compares_last_against_first(): + """A closed path whose final point repeats its first loses the duplicate.""" + pts = [(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 0.0)] + assert deduplicate(pts, closed=True) == [(0.0, 0.0), (1.0, 0.0), (1.0, 1.0)] + + +def test_collinearity_tolerance_is_relative_not_absolute(): + """ + A 1e-6 deviation over a 1 mm chord is not collinear; the same deviation over + a 1e6 mm chord is. An absolute epsilon would call both the same. + """ + assert not is_collinear([(0.0, 0.0), (0.5, 1e-6), (1.0, 0.0)]) + assert is_collinear([(0.0, 0.0), (5e5, 1e-6), (1e6, 0.0)]) + + +def test_path_merge_collinear_drops_the_midpoint(): + pts = [(0.0, 0.0), (5.0, 0.0), (10.0, 0.0), (10.0, 10.0), (0.0, 10.0)] + merged = path_merge_collinear(pts, closed=True) + assert (5.0, 0.0) not in merged + assert len(merged) == 4 + + +def test_path_merge_collinear_keeps_a_genuine_corner(): + assert len(path_merge_collinear(SQUARE10, closed=True)) == 4 + + +# ---------------------------------------------------------------------------- +# round_corners -- segmentation +# ---------------------------------------------------------------------------- + +def test_right_angle_corner_yields_twelve_points(): + """ + Half-angle 45, so n = max(3, ceil((90-45)/180 * 48)) = 12. This is the + single number the whole area comparison rests on. + """ + rounded = round_corners(SQUARE10, 2.0) + assert len(rounded) == 4 * 12 + + +def test_segment_count_is_independent_of_radius(): + """$fn segmentation depends on swept angle only. Halving r changes nothing.""" + assert len(round_corners(SQUARE10, 2.0)) == len(round_corners(SQUARE10, 1.0)) + + +def test_zero_radius_leaves_the_vertex_untouched(): + assert round_corners(SQUARE10, [0.0, 0.0, 0.0, 0.0]) == SQUARE10 + + +def test_mixed_radii_round_only_the_named_corners(): + rounded = round_corners(SQUARE10, [2.0, 0.0, 0.0, 0.0]) + assert len(rounded) == 12 + 3 + for v in SQUARE10[1:]: + assert v in rounded + assert (0.0, 0.0) not in rounded + + +def test_facet_count_is_a_parameter_not_a_constant(): + assert len(round_corners(SQUARE10, 2.0, fn=96)) == 4 * 24 + + +# ---------------------------------------------------------------------------- +# round_corners -- geometry +# ---------------------------------------------------------------------------- + +def test_rounded_square_stays_inside_the_original(): + rounded = round_corners(SQUARE10, 2.0) + for x, y in rounded: + assert -1e-9 <= x <= 10.0 + 1e-9 + assert -1e-9 <= y <= 10.0 + 1e-9 + + +def test_rounding_removes_area_and_the_loss_is_bounded(): + """ + A square loses one square minus one inscribed circle across its four + corners: 4r^2 - pi*r^2. The polygonal arc removes slightly more, so the + measured loss sits just above the exact figure. + """ + def shoelace(path): + n = len(path) + return abs(sum(path[i][0] * path[(i + 1) % n][1] + - path[(i + 1) % n][0] * path[i][1] + for i in range(n))) / 2.0 + + r = 2.0 + lost = shoelace(SQUARE10) - shoelace(round_corners(SQUARE10, r)) + exact = 4 * r * r - math.pi * r * r + assert exact < lost < exact * 1.05 + + +def test_tangent_points_sit_on_the_original_edges(): + """The roundover must start and end on the edges it replaces, not inside.""" + rounded = round_corners(SQUARE10, 2.0) + on_edge = [p for p in rounded + if approx(p[0], 0.0) or approx(p[0], 10.0) + or approx(p[1], 0.0) or approx(p[1], 10.0)] + assert len(on_edge) == 8 + + +def test_rounding_is_mirror_symmetric(): + """ + Chirality was a real failure in earlier revisions -- the morphological + closing it replaced quietly made symmetric profiles handed. + """ + rounded = round_corners(SQUARE10, 2.0) + xs = sorted(round(p[0], 9) for p in rounded) + mirrored = sorted(round(10.0 - p[0], 9) for p in rounded) + assert xs == mirrored + + +def test_corner_segment_count_follows_the_half_angle(): + """ + A 60-degree corner has half-angle 30, so n = max(3, ceil(60/180 * 48)) = 16. + An equilateral triangle has three of them. + """ + tri = [(0.0, 0.0), (10.0, 0.0), (5.0, 10.0 * math.sqrt(3) / 2)] + assert len(round_corners(tri, 0.5)) == 3 * 16 + + +def test_a_sharper_corner_gets_more_segments_than_a_blunter_one(): + """Segment count rises as the corner closes up, since the arc sweeps more.""" + sharp = [(0.0, 0.0), (10.0, 0.0), (5.0, 1.0)] + blunt = [(0.0, 0.0), (10.0, 0.0), (5.0, 20.0)] + assert len(round_corners(sharp, 0.05)) > len(round_corners(blunt, 0.05)) + + +# ---------------------------------------------------------------------------- +# round_corners -- failure behaviour +# ---------------------------------------------------------------------------- + +def test_roundover_too_large_raises_rather_than_clamping(): + """ + BOSL2 asserts here. Clamping instead would let a build succeed where the + reference aborted, which is a silent divergence from the oracle. + """ + with pytest.raises(RoundoverTooLarge): + round_corners(SQUARE10, 20.0) + + +def test_the_error_reports_the_scale_factors(): + with pytest.raises(RoundoverTooLarge) as exc: + round_corners(SQUARE10, 20.0) + assert "multiply them by this vector" in str(exc.value) + + +def test_the_exact_fit_boundary_falls_just_short_in_floating_point(): + """ + r = 5 on a 10 mm square is the exact-fit case on paper: each edge carries + two setbacks of 5/tan(45). But tan(45) is 0.9999999999999999, not 1, so the + setbacks total fractionally over the edge and the scale factor lands at + 0.9999999999999998. BOSL2 asserts here, and so must this port -- OpenSCAD + converts degrees to radians the same way and gets the same last bit. + + Asserting that the exact-fit case *passes* would look reasonable and would + quietly diverge from the reference. A hair under fits. + """ + with pytest.raises(RoundoverTooLarge): + round_corners(SQUARE10, 5.0) + round_corners(SQUARE10, 4.999999999) + + +def test_repeated_point_with_nonzero_rounding_is_rejected(): + with pytest.raises(ValueError, match="Repeated point"): + round_corners([(0.0, 0.0), (0.0, 0.0), (10.0, 0.0), (5.0, 5.0)], 1.0) + + +def test_path_too_short_is_rejected(): + with pytest.raises(ValueError, match="Length must be 3 or more"): + round_corners([(0.0, 0.0), (1.0, 1.0)], 1.0) + + +def test_radius_list_length_must_match(): + with pytest.raises(ValueError, match="does not match path length"): + round_corners(SQUARE10, [1.0, 1.0]) + + +def test_a_very_blunt_corner_still_gets_the_three_point_floor(): + """ + Half-angle 85 gives ceil(5/180 * 48) = 2, below BOSL2's max(3, ...) floor. + Without the floor an arc would degenerate to a single chord, which reads as + a slightly wrong area rather than as an error -- so it needs pinning. + """ + a = math.radians(10.0) + path = [(-10.0, 0.0), (0.0, 0.0), + (10.0 * math.cos(a), 10.0 * math.sin(a)), (0.0, -20.0)] + rounded = round_corners(path, [0.0, 0.05, 0.0, 0.0]) + assert len(rounded) == 3 + 3