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])