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.
This commit is contained in:
2026-08-19 03:30:03 -05:00
parent 545eee7217
commit 38ea024fdc
3 changed files with 665 additions and 0 deletions
+20
View File
@@ -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,
)
+325
View File
@@ -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
+320
View File
@@ -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])