From 38ea024fdccde85fa1dda1e3490b23221ce36617 Mon Sep 17 00:00:00 2001 From: TheRON Date: Wed, 19 Aug 2026 03:30:03 -0500 Subject: [PATCH] geom: Shapely-backed region layer Booleans, decomposition, area, simplicity, hull, bounds, mitred offset and vertex cleaning. Booleans go to GEOS, which is the reason Shapely was chosen. The decomposition does not. BOSL2 region_parts counts by nesting parity, not connectivity: a path takes a level from how many others contain the midpoint of its first edge, even levels are outer boundaries, their odd children are holes. SECTION_PARTS == 1 is an exact assertion, and Shapely agreeing with that count is a coincidence that holds for well-formed input and not otherwise, so the decomposition is transcribed and both the part count and the area derive from it. is_region_simple is treated as a manifold precondition rather than a diagnostic. An outline that touches itself measures perfectly and cannot be tessellated, so it must fail here and not at export. Developed against Shapely 2.1.2 / GEOS 3.13.1, matching CT 100. Boolean results on near-degenerate geometry can shift between GEOS releases; if the oracle ever disagrees by one part after an upgrade, look there first. 39 tests, all arithmetic on rectangles. Mutation run found a real gap: nothing distinguished on-boundary from outside in the nesting probe until a shared-edge case was added. Eight mutations now caught. Oracle acceptance still skips; 236 unchanged. --- src/mechcomp/geom/__init__.py | 20 +++ src/mechcomp/geom/region.py | 325 ++++++++++++++++++++++++++++++++++ tests/test_region.py | 320 +++++++++++++++++++++++++++++++++ 3 files changed, 665 insertions(+) create mode 100644 src/mechcomp/geom/region.py create mode 100644 tests/test_region.py diff --git a/src/mechcomp/geom/__init__.py b/src/mechcomp/geom/__init__.py index ec243eb..fe7c563 100644 --- a/src/mechcomp/geom/__init__.py +++ b/src/mechcomp/geom/__init__.py @@ -60,3 +60,23 @@ from .rounding import ( # noqa: F401 round_corners, segs, ) +from .region import ( # noqa: F401 + Region, + area, + as_region, + clean_region, + difference, + from_shapely, + hull_region, + intersection, + is_path_simple, + is_region_simple, + nparts, + offset_path, + point_in_polygon, + pointlist_bounds, + region_area, + region_parts, + to_shapely, + union, +) diff --git a/src/mechcomp/geom/region.py b/src/mechcomp/geom/region.py new file mode 100644 index 0000000..204c086 --- /dev/null +++ b/src/mechcomp/geom/region.py @@ -0,0 +1,325 @@ +""" +Region operations: the boolean and measurement layer ``sb-join`` sits on. + +A *region* is BOSL2's: a list of closed paths, where nesting determines solid +from void. That representation is kept rather than replaced by Shapely +geometries, because it is what the reference passes around and what the +generators' call sites expect. + +**Why the decomposition is transcribed rather than delegated.** BOSL2's +``region_parts`` counts by nesting *parity*, not connectivity: each path takes a +level from how many other paths contain the midpoint of its first edge, and +even-level paths are outer boundaries with their odd-level children as holes. +``SECTION_PARTS == 1`` is an exact assertion in the oracle, and Shapely's +component count agreeing with BOSL2's parity count is a coincidence that holds +for well-formed input and not otherwise. Reproducing the decomposition means the +part count and the area come from the same reading of the geometry. + +Booleans are delegated to GEOS, which is the reason for choosing Shapely. The +results are converted back to path lists so nothing downstream needs to know. + +Pinned against BOSL2 ``92d697c2856de2fed93a33e858068589cefc2898``. Developed +against Shapely 2.1.2 / GEOS 3.13.1, matching CT 100; boolean results on +near-degenerate geometry can differ between GEOS releases, so that pairing is +worth keeping in step. +""" + +from __future__ import annotations + +import math +from typing import List, Optional, Sequence, Tuple + +from shapely.geometry import MultiPoint, MultiPolygon, Polygon +from shapely.geometry.base import BaseGeometry +from shapely.ops import unary_union + +from .primitives import SB_EPS, sb_signed_area +from .rounding import EPSILON, path_merge_collinear + +Point = Tuple[float, float] +Path = Sequence[Point] +Region = List[List[Point]] + +# Large enough that a mitred offset is never silently bevelled at the angles +# these profiles use, which run to fairly sharp apexes. +_MITRE_LIMIT = 1e6 + + +# ---------------------------------------------------------------------------- +# Point-in-polygon (geometry.scad: point_in_polygon, winding number) +# ---------------------------------------------------------------------------- + +def point_in_polygon(pt: Point, poly: Path, eps: float = EPSILON) -> int: + """ + 1 inside, 0 on the boundary, -1 outside. + + On-boundary is its own answer rather than folded into inside or outside, + because ``region_parts`` treats a point on a boundary as contained and the + distinction changes nesting levels. + """ + n = len(poly) + for i in range(n): + a, b = poly[i], poly[(i + 1) % n] + abx, aby = b[0] - a[0], b[1] - a[1] + l2 = abx * abx + aby * aby + if l2 < SB_EPS: + if math.hypot(pt[0] - a[0], pt[1] - a[1]) <= eps: + return 0 + continue + t = ((pt[0] - a[0]) * abx + (pt[1] - a[1]) * aby) / l2 + t = max(0.0, min(1.0, t)) + if math.hypot(pt[0] - (a[0] + t * abx), pt[1] - (a[1] + t * aby)) <= eps: + return 0 + + wind = 0 + for i in range(n): + a, b = poly[i], poly[(i + 1) % n] + if a[1] <= pt[1] < b[1] or b[1] <= pt[1] < a[1]: + cross = ((b[0] - a[0]) * (pt[1] - a[1]) + - (b[1] - a[1]) * (pt[0] - a[0])) + if a[1] < b[1]: + wind += 1 if cross > 0 else 0 + else: + wind -= 1 if cross < 0 else 0 + return 1 if wind != 0 else -1 + + +# ---------------------------------------------------------------------------- +# Decomposition (regions.scad: region_parts) +# ---------------------------------------------------------------------------- + +def _cw(path: Path) -> List[Point]: + return list(path) if sb_signed_area(path) < 0 else list(reversed(path)) + + +def _ccw(path: Path) -> List[Point]: + return list(path) if sb_signed_area(path) >= 0 else list(reversed(path)) + + +def region_parts(rgn: Sequence[Path]) -> List[List[List[Point]]]: + """ + Split a region into connected pieces, each ``[outer, *holes]``. + + The outer boundary comes back clockwise and its holes counter-clockwise, + matching the reference, because ``region_area`` depends on that convention + to make holes subtract. + """ + paths = [list(p) for p in rgn] + n = len(paths) + if n == 0: + return [] + + inside = [] + for i in range(n): + pt = ((paths[i][0][0] + paths[i][1][0]) / 2.0, + (paths[i][0][1] + paths[i][1][1]) / 2.0) + inside.append([0 if i == j else + (1 if point_in_polygon(pt, paths[j]) >= 0 else 0) + for j in range(n)]) + + level = [sum(row) for row in inside] + + out: List[List[List[Point]]] = [] + for i in range(n): + if level[i] % 2 != 0: + continue + holes = [j for j in range(n) + if level[j] == level[i] + 1 and inside[j][i] == 1] + out.append([_cw(paths[i])] + [_ccw(paths[j]) for j in holes]) + return out + + +def region_area(rgn: Sequence[Path]) -> float: + """Total enclosed area, holes subtracted.""" + return -sum(sb_signed_area(poly) + for part in region_parts(rgn) for poly in part) + + +def nparts(rgn: Sequence[Path]) -> int: + """Number of connected solids. Empty reads as zero, as the reference does.""" + return 0 if len(rgn) == 0 else len(region_parts(rgn)) + + +def area(rgn: Sequence[Path]) -> float: + """Empty-safe area: a legitimately empty result reads as zero.""" + return 0.0 if len(rgn) == 0 else region_area(rgn) + + +def as_region(x) -> Region: + """ + Accept a bare path where a region is expected. + + BOSL2's boolean functions return a bare ``[]`` when a result is empty, which + is not a valid region; every measurement goes through here so an empty + result reads as zero rather than raising. + """ + if not x: + return [] + first = x[0] + if isinstance(first, (tuple, list)) and len(first) == 2 \ + and isinstance(first[0], (int, float)): + return [list(x)] + return [list(p) for p in x] + + +# ---------------------------------------------------------------------------- +# Shapely bridge +# ---------------------------------------------------------------------------- + +def to_shapely(rgn: Sequence[Path]) -> BaseGeometry: + """Build a Shapely geometry from the reference's own decomposition.""" + polys = [] + for part in region_parts(rgn): + shell = part[0] + holes = part[1:] + p = Polygon(shell, holes) + if not p.is_valid: + p = p.buffer(0) + if not p.is_empty: + polys.append(p) + if not polys: + return Polygon() + return polys[0] if len(polys) == 1 else MultiPolygon( + [g for p in polys for g in (p.geoms if isinstance(p, MultiPolygon) else [p])]) + + +def from_shapely(geom: BaseGeometry) -> Region: + """ + Flatten a Shapely geometry back to a list of closed paths. + + Shapely repeats the first coordinate to close a ring; the reference's paths + are implicitly closed, so the duplicate is dropped. + """ + if geom.is_empty: + return [] + geoms = geom.geoms if isinstance(geom, MultiPolygon) else [geom] + out: Region = [] + for g in geoms: + if not isinstance(g, Polygon) or g.is_empty: + continue + out.append([(x, y) for x, y in list(g.exterior.coords)[:-1]]) + for ring in g.interiors: + out.append([(x, y) for x, y in list(ring.coords)[:-1]]) + return out + + +# ---------------------------------------------------------------------------- +# Booleans +# ---------------------------------------------------------------------------- + +def union(regions: Sequence[Sequence[Path]]) -> Region: + shapes = [to_shapely(r) for r in regions if len(r) > 0] + if not shapes: + return [] + return from_shapely(unary_union(shapes)) + + +def difference(a: Sequence[Path], b: Sequence[Path]) -> Region: + if len(a) == 0: + return [] + if len(b) == 0: + return as_region(a) + return from_shapely(to_shapely(a).difference(to_shapely(b))) + + +def intersection(a: Sequence[Path], b: Sequence[Path]) -> Region: + if len(a) == 0 or len(b) == 0: + return [] + return from_shapely(to_shapely(a).intersection(to_shapely(b))) + + +# ---------------------------------------------------------------------------- +# Simplicity -- the manifold precondition +# ---------------------------------------------------------------------------- + +def is_path_simple(path: Path, eps: float = EPSILON) -> bool: + """A closed path that neither crosses nor touches itself away from its seam.""" + pts = list(path) + if len(pts) < 3: + return False + ring = Polygon(pts + [pts[0]]).exterior + return bool(ring.is_simple) + + +def is_region_simple(rgn: Sequence[Path], eps: float = EPSILON) -> bool: + """ + Every path simple, and no two paths meeting at all -- touching included. + + This is not a cosmetic check. A section that touches itself at a point is a + valid 2D outline and cannot be tessellated, so it measures perfectly and + then fails to extrude into a sealed solid. + """ + paths = [list(p) for p in rgn] + for p in paths: + if not is_path_simple(p, eps): + return False + rings = [Polygon(p + [p[0]]).exterior for p in paths] + for i in range(len(rings)): + for j in range(i + 1, len(rings)): + if rings[i].intersects(rings[j]): + return False + return True + + +# ---------------------------------------------------------------------------- +# Hull, bounds, offset, cleanup +# ---------------------------------------------------------------------------- + +def hull_region(rgn: Sequence[Path]) -> List[Point]: + """Convex hull of every point in the region.""" + pts = [tuple(p) for path in rgn for p in path] + if len(pts) < 3: + return list(pts) + hull = MultiPoint(pts).convex_hull + if isinstance(hull, Polygon): + return [(x, y) for x, y in list(hull.exterior.coords)[:-1]] + # Collinear input degenerates to a line or a point; return it unchanged + # rather than inventing an area the caller would then measure. + return list(pts) + + +def pointlist_bounds(pts: Sequence[Point]) -> List[Point]: + """``[[min_x, min_y], [max_x, max_y]]``, as the reference returns it.""" + xs = [p[0] for p in pts] + ys = [p[1] for p in pts] + return [(min(xs), min(ys)), (max(xs), max(ys))] + + +def offset_path(path: Path, delta: float, closed: bool = True) -> List[Point]: + """ + Mitred offset of a closed path. Positive ``delta`` grows it. + + Mitred rather than rounded: BOSL2's ``offset`` with ``delta`` keeps corners + sharp, and the ring envelope's corners are rounded afterwards by an explicit + ``round_corners`` call, not by the offset itself. A rounded join here would + round them twice and by the wrong rule. + """ + pts = _ccw(path) + poly = Polygon(pts) + if not poly.is_valid: + poly = poly.buffer(0) + grown = poly.buffer(delta, join_style="mitre", mitre_limit=_MITRE_LIMIT) + if grown.is_empty: + return [] + if isinstance(grown, MultiPolygon): + grown = max(grown.geoms, key=lambda g: g.area) + return [(x, y) for x, y in list(grown.exterior.coords)[:-1]] + + +def clean_region(rgn: Sequence[Path], eps: float = EPSILON) -> Region: + """ + Drop duplicate and collinear vertices from every path, discarding any path + left with fewer than three. + + Exact butt joints and zero-radius fillets produce coincident or perfectly + 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. Cleaning once, at the end, removes that entire class + of failure. + """ + out: Region = [] + for path in rgn: + merged = path_merge_collinear(path, closed=True, eps=eps) + if len(merged) >= 3: + out.append(merged) + return out diff --git a/tests/test_region.py b/tests/test_region.py new file mode 100644 index 0000000..7ad5ab0 --- /dev/null +++ b/tests/test_region.py @@ -0,0 +1,320 @@ +""" +Unit tests for the region layer. + +Every expected area and count here is arithmetic on rectangles, checkable by +hand. The nesting cases matter most: BOSL2 counts parts by nesting parity rather +than by connectivity, and ``SECTION_PARTS == 1`` is an exact assertion in the +oracle, so a decomposition that is merely usually-right is not good enough. + +``is_region_simple`` is tested as a manifold precondition rather than as a +diagnostic, because that is what it is: an outline that touches itself measures +perfectly and will not extrude into a sealed solid. +""" + +from __future__ import annotations + +import math + +import pytest + +from mechcomp.geom.region import ( + area, + as_region, + clean_region, + difference, + from_shapely, + hull_region, + intersection, + is_path_simple, + is_region_simple, + nparts, + offset_path, + point_in_polygon, + pointlist_bounds, + region_area, + region_parts, + to_shapely, + union, +) + + +def rect(x0, y0, x1, y1): + return [(x0, y0), (x1, y0), (x1, y1), (x0, y1)] + + +SQ = rect(0, 0, 10, 10) +HOLE = rect(2, 2, 8, 8) +INNER = rect(4, 4, 6, 6) +FAR = rect(20, 0, 30, 10) + + +# ---------------------------------------------------------------------------- +# point_in_polygon +# ---------------------------------------------------------------------------- + +def test_point_in_polygon_three_way_answer(): + assert point_in_polygon((5.0, 5.0), SQ) == 1 + assert point_in_polygon((50.0, 5.0), SQ) == -1 + assert point_in_polygon((0.0, 5.0), SQ) == 0 + + +def test_boundary_is_its_own_answer_not_folded_into_inside(): + """ + region_parts treats on-boundary as contained. Collapsing it into inside or + outside would shift nesting levels and change the part count. + """ + assert point_in_polygon((10.0, 10.0), SQ) == 0 + assert point_in_polygon((5.0, 0.0), SQ) == 0 + + +def test_point_in_polygon_is_winding_not_bounding_box(): + ell = [(0, 0), (10, 0), (10, 4), (4, 4), (4, 10), (0, 10)] + assert point_in_polygon((8.0, 8.0), ell) == -1 + assert point_in_polygon((2.0, 2.0), ell) == 1 + + +# ---------------------------------------------------------------------------- +# Decomposition and area +# ---------------------------------------------------------------------------- + +def test_plain_square(): + assert nparts([SQ]) == 1 + assert area([SQ]) == pytest.approx(100.0) + + +def test_hole_subtracts_and_stays_one_part(): + assert nparts([SQ, HOLE]) == 1 + assert area([SQ, HOLE]) == pytest.approx(64.0) + + +def test_disjoint_paths_are_separate_parts(): + assert nparts([SQ, FAR]) == 2 + assert area([SQ, FAR]) == pytest.approx(200.0) + + +def test_nesting_parity_not_connectivity(): + """ + Square, hole, and an island inside the hole. Levels 0, 1, 2 -- so two parts, + and the island's area is added back rather than subtracted. + """ + assert nparts([SQ, HOLE, INNER]) == 2 + assert area([SQ, HOLE, INNER]) == pytest.approx(68.0) + + +def test_part_ordering_convention_outer_clockwise_holes_ccw(): + """region_area depends on this winding to make holes subtract.""" + from mechcomp.geom.primitives import sb_signed_area + parts = region_parts([SQ, HOLE]) + assert len(parts) == 1 + outer, hole = parts[0] + assert sb_signed_area(outer) < 0 + assert sb_signed_area(hole) > 0 + + +def test_path_order_does_not_change_the_answer(): + """Nesting is derived, not assumed from input order.""" + assert area([HOLE, SQ]) == pytest.approx(64.0) + assert nparts([HOLE, SQ]) == 1 + + +def test_empty_region_reads_as_zero_not_an_error(): + """BOSL2's booleans return a bare [] for an empty result.""" + assert nparts([]) == 0 + assert area([]) == 0.0 + assert region_parts([]) == [] + + +# ---------------------------------------------------------------------------- +# as_region +# ---------------------------------------------------------------------------- + +def test_as_region_wraps_a_bare_path(): + assert as_region(SQ) == [SQ] + + +def test_as_region_passes_a_region_through(): + assert as_region([SQ, HOLE]) == [SQ, HOLE] + + +def test_as_region_of_empty_is_empty(): + assert as_region([]) == [] + + +# ---------------------------------------------------------------------------- +# Booleans +# ---------------------------------------------------------------------------- + +def test_union_of_disjoint_squares_keeps_both(): + u = union([[SQ], [FAR]]) + assert nparts(u) == 2 + assert area(u) == pytest.approx(200.0) + + +def test_union_of_overlapping_squares_merges_them(): + u = union([[SQ], [rect(5, 0, 15, 10)]]) + assert nparts(u) == 1 + assert area(u) == pytest.approx(150.0) + + +def test_difference_makes_a_hole(): + d = difference([SQ], [HOLE]) + assert nparts(d) == 1 + assert area(d) == pytest.approx(64.0) + assert len(d) == 2 + + +def test_difference_can_split_one_solid_into_two(): + """A cut straight across is how SECTION_PARTS stops being 1.""" + d = difference([SQ], [rect(4, -1, 6, 11)]) + assert nparts(d) == 2 + assert area(d) == pytest.approx(80.0) + + +def test_difference_with_nothing_returns_the_original(): + assert area(difference([SQ], [])) == pytest.approx(100.0) + + +def test_intersection_of_overlapping_squares(): + i = intersection([SQ], [rect(5, 5, 15, 15)]) + assert area(i) == pytest.approx(25.0) + + +def test_union_ignores_empty_operands(): + assert area(union([[SQ], []])) == pytest.approx(100.0) + + +def test_round_trip_through_shapely_preserves_area_and_parts(): + for rgn in ([SQ], [SQ, HOLE], [SQ, HOLE, INNER], [SQ, FAR]): + back = from_shapely(to_shapely(rgn)) + assert area(back) == pytest.approx(area(rgn)) + assert nparts(back) == nparts(rgn) + + +def test_from_shapely_drops_the_repeated_closing_point(): + """Shapely closes rings explicitly; reference paths are implicitly closed.""" + back = from_shapely(to_shapely([SQ])) + assert len(back[0]) == 4 + + +# ---------------------------------------------------------------------------- +# Simplicity -- manifold precondition +# ---------------------------------------------------------------------------- + +def test_simple_region_is_simple(): + assert is_region_simple([SQ, FAR]) + + +def test_self_intersecting_path_is_not_simple(): + bowtie = [(0, 0), (10, 10), (10, 0), (0, 10)] + assert not is_path_simple(bowtie) + assert not is_region_simple([bowtie]) + + +def test_paths_that_merely_touch_are_not_simple(): + """ + Touching at an edge or a point is the case that extrudes into something + untessellatable. It has to fail here, not at STL export. + """ + assert not is_region_simple([SQ, rect(10, 0, 20, 10)]) + assert not is_region_simple([SQ, rect(10, 10, 20, 20)]) + + +def test_a_hole_inside_an_outline_is_still_simple(): + """Nesting is the normal case and must not read as an intersection.""" + assert is_region_simple([SQ, HOLE]) + + +# ---------------------------------------------------------------------------- +# Hull, bounds, offset +# ---------------------------------------------------------------------------- + +def test_hull_of_a_convex_path_is_itself(): + assert len(hull_region([SQ])) == 4 + assert area([hull_region([SQ])]) == pytest.approx(100.0) + + +def test_hull_spans_disjoint_paths(): + assert area([hull_region([SQ, FAR])]) == pytest.approx(300.0) + + +def test_hull_of_a_concave_path_fills_the_notch(): + ell = [(0, 0), (10, 0), (10, 4), (4, 4), (4, 10), (0, 10)] + assert area([hull_region([ell])]) > area([ell]) + + +def test_pointlist_bounds(): + assert pointlist_bounds(SQ) == [(0.0, 0.0), (10.0, 10.0)] + + +def test_offset_grows_by_delta_on_every_side(): + grown = offset_path(SQ, 2.0) + assert area([grown]) == pytest.approx(196.0) + assert pointlist_bounds(grown) == [(-2.0, -2.0), (12.0, 12.0)] + + +def test_offset_is_mitred_not_rounded(): + """ + A rounded join would give 4 corner arcs and an area of 100 + 4*10*2 + pi*4. + Mitred keeps the corners sharp, so the envelope stays a 4-vertex square and + the explicit round_corners call afterwards is the only rounding applied. + """ + grown = offset_path(SQ, 2.0) + assert len(grown) == 4 + assert area([grown]) != pytest.approx(100.0 + 80.0 + math.pi * 4.0) + + +def test_negative_offset_shrinks(): + assert area([offset_path(SQ, -2.0)]) == pytest.approx(36.0) + + +def test_offset_preserves_a_triangle_shape(): + tri = [(0.0, 0.0), (10.0, 0.0), (5.0, 8.0)] + grown = offset_path(tri, 1.0) + assert len(grown) == 3 + assert area([grown]) > area([tri]) + + +# ---------------------------------------------------------------------------- +# clean_region +# ---------------------------------------------------------------------------- + +def test_clean_drops_collinear_vertices(): + with_flat = [(0, 0), (5, 0), (10, 0), (10, 10), (0, 10)] + cleaned = clean_region([with_flat]) + assert len(cleaned[0]) == 4 + assert area(cleaned) == pytest.approx(100.0) + + +def test_clean_drops_duplicate_vertices(): + dup = [(0, 0), (0, 0), (10, 0), (10, 10), (0, 10)] + assert len(clean_region([dup])[0]) == 4 + + +def test_clean_discards_a_path_left_too_short(): + """A degenerate sliver collapses to fewer than three points and is dropped.""" + sliver = [(0.0, 0.0), (5.0, 0.0), (10.0, 0.0)] + assert clean_region([sliver]) == [] + + +def test_clean_preserves_a_genuine_outline(): + assert area(clean_region([SQ, HOLE])) == pytest.approx(64.0) + + +def test_on_boundary_probe_counts_as_contained(): + """ + Two squares sharing an edge. Each path's nesting is probed at the midpoint + of its first edge, and for the upper square that point lies exactly on the + lower square's top edge. + + BOSL2 counts on-boundary as contained, so the upper square reads as nested + and becomes a hole: one part, and the areas cancel. Counting it as outside + would give two parts and 200. The result is geometrically odd either way -- + which is exactly why is_region_simple rejects touching paths before any of + this is measured. + """ + upper = rect(0, 0, 10, 10) + lower = rect(0, -10, 10, 0) + assert point_in_polygon((5.0, 0.0), lower) == 0 + assert nparts([upper, lower]) == 1 + assert area([upper, lower]) == pytest.approx(0.0) + assert not is_region_simple([upper, lower])