Antoine
78f56a68b0
Root cause: typed segment offsetting created self-intersecting geometry at concave corners (notches). Triangle's PSLG boundary didn't match the plotted inset contour, allowing vertices 7+ mm outside. Changes: - _build_inner_plate: always use Shapely buffer(-w_frame) (robust at concavities) - _sample_ring: use simplified polygon vertices + interpolated points on long edges (preserves tight features without vertex clustering) - Plot uses same inner_plate from triangulation (no mismatch) - Post-process: snap any residual outside vertices to boundary - Result: 0 vertices outside inner plate (was 10, up to 7.45mm)
381 lines
14 KiB
Python
381 lines
14 KiB
Python
"""
|
|
Isogrid triangulation — generates a proper isogrid rib pattern.
|
|
|
|
Uses Shewchuk's Triangle library (constrained Delaunay + quality refinement)
|
|
for boundary-conforming, clean triangulations.
|
|
|
|
Strategy:
|
|
1. Build PSLG (Planar Straight Line Graph):
|
|
- Outer contour = inset boundary (w_frame offset inward)
|
|
- Holes = keepout circles (d_keep offset outward from holes)
|
|
2. Triangulate with quality constraints:
|
|
- Minimum angle ~25° (no slivers)
|
|
- Maximum area varies by density heatmap (smaller near holes, larger away)
|
|
3. Result: clean mesh that fills 100% of valid area, with density-graded
|
|
triangle sizes and edges parallel to boundary/hole contours.
|
|
"""
|
|
|
|
import numpy as np
|
|
import triangle as tr
|
|
from shapely.geometry import Polygon, Point, LinearRing, MultiPolygon
|
|
from shapely.geometry.base import BaseGeometry
|
|
from shapely.ops import unary_union
|
|
|
|
from src.shared.arc_utils import inset_arc, typed_segments_to_polyline
|
|
from .density_field import evaluate_density, density_to_spacing
|
|
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# Helpers
|
|
# ---------------------------------------------------------------------------
|
|
|
|
def _geometry_to_list(geom: BaseGeometry) -> list:
|
|
if geom.is_empty:
|
|
return []
|
|
if geom.geom_type == 'Polygon':
|
|
return [geom]
|
|
if geom.geom_type == 'MultiPolygon':
|
|
return list(geom.geoms)
|
|
return list(getattr(geom, 'geoms', []))
|
|
|
|
|
|
def _ring_to_segments(coords: np.ndarray, start_idx: int):
|
|
"""Convert a ring (Nx2 array, NOT closed) to vertex array + segment index pairs."""
|
|
n = len(coords)
|
|
segments = []
|
|
for i in range(n):
|
|
segments.append([start_idx + i, start_idx + (i + 1) % n])
|
|
return segments
|
|
|
|
|
|
def _sample_ring(ring, spacing: float) -> np.ndarray:
|
|
"""Sample points along a Shapely ring at given spacing.
|
|
|
|
Uses Shapely simplify() to reduce vertex count on curved buffer segments,
|
|
then adds vertices from the simplified ring plus interpolated points on
|
|
long edges. This preserves corners/notches while avoiding vertex clusters.
|
|
"""
|
|
# Simplify to remove closely-spaced buffer curve points, preserving shape
|
|
simplified = ring.simplify(spacing * 0.15, preserve_topology=True)
|
|
coords = np.array(simplified.coords)
|
|
if len(coords) > 1 and np.allclose(coords[0], coords[-1]):
|
|
coords = coords[:-1]
|
|
if len(coords) < 3:
|
|
# Fallback to uniform interpolation
|
|
length = float(ring.length)
|
|
if length < 1e-9:
|
|
return np.empty((0, 2), dtype=np.float64)
|
|
n = max(int(np.ceil(length / max(spacing, 1e-3))), 8)
|
|
pts = []
|
|
for i in range(n):
|
|
p = ring.interpolate(i / n, normalized=True)
|
|
pts.append([p.x, p.y])
|
|
return np.array(pts, dtype=np.float64)
|
|
|
|
pts = []
|
|
n = len(coords)
|
|
for i in range(n):
|
|
p1 = coords[i]
|
|
p2 = coords[(i + 1) % n]
|
|
pts.append(p1.tolist())
|
|
|
|
# Add interpolated points on long edges
|
|
edge_len = np.linalg.norm(p2 - p1)
|
|
if edge_len > spacing * 1.5:
|
|
n_sub = int(np.ceil(edge_len / spacing))
|
|
for j in range(1, n_sub):
|
|
t = j / n_sub
|
|
pts.append((p1 + t * (p2 - p1)).tolist())
|
|
|
|
return np.array(pts, dtype=np.float64)
|
|
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# Step 1 — Build the inset plate polygon (frame inner edge)
|
|
# ---------------------------------------------------------------------------
|
|
|
|
def _build_inner_plate(geometry, params) -> Polygon:
|
|
"""Offset sandbox boundary inward by w_frame.
|
|
|
|
Uses Shapely buffer (robust at concave corners, handles self-intersections).
|
|
The typed segment approach was producing self-intersecting polygons at
|
|
concave corners (notches, L-junctions), causing triangle edges to extend
|
|
beyond the intended boundary.
|
|
"""
|
|
w_frame = float(params.get('w_frame', 8.0))
|
|
plate_poly = Polygon(geometry['outer_boundary'])
|
|
|
|
inner_plate = plate_poly.buffer(-w_frame, resolution=16)
|
|
if inner_plate.is_empty or not inner_plate.is_valid:
|
|
inner_plate = plate_poly
|
|
return inner_plate
|
|
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# Step 2 — Build hole keepout polygons
|
|
# ---------------------------------------------------------------------------
|
|
|
|
def _build_keepouts(geometry, params) -> list:
|
|
"""Build list of keepout Polygons (one per hole)."""
|
|
d_keep = float(params['d_keep'])
|
|
keepouts = []
|
|
|
|
for hole in geometry.get('holes', []):
|
|
diameter = float(hole.get('diameter', 10.0) or 10.0)
|
|
if hole.get('is_circular', False) and 'center' in hole:
|
|
cx, cy = hole['center']
|
|
hole_radius = diameter / 2.0
|
|
keepout_radius = hole_radius + d_keep * hole_radius
|
|
keepout = Point(cx, cy).buffer(keepout_radius, resolution=32)
|
|
keepouts.append(keepout)
|
|
else:
|
|
hb = np.asarray(hole.get('boundary', []), dtype=float)
|
|
if len(hb) < 3:
|
|
continue
|
|
ring = LinearRing(hb)
|
|
keepout = Polygon(ring).buffer(max(d_keep, 0.0))
|
|
if not keepout.is_empty:
|
|
keepouts.append(keepout)
|
|
|
|
return keepouts
|
|
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# Step 3 — Build PSLG for Triangle library
|
|
# ---------------------------------------------------------------------------
|
|
|
|
def _build_pslg(inner_plate: Polygon, keepouts: list, boundary_spacing: float):
|
|
"""
|
|
Build a Planar Straight Line Graph (PSLG) for the Triangle library.
|
|
|
|
Returns dict with:
|
|
'vertices': Nx2 array
|
|
'segments': Mx2 array of vertex index pairs
|
|
'holes': Hx2 array of hole marker points
|
|
"""
|
|
vertices = []
|
|
segments = []
|
|
hole_markers = []
|
|
|
|
# Outer boundary of the inner plate
|
|
# Handle MultiPolygon (take largest piece)
|
|
polys = _geometry_to_list(inner_plate)
|
|
if not polys:
|
|
return None
|
|
main_poly = max(polys, key=lambda p: p.area)
|
|
|
|
# Sample the outer ring
|
|
outer_pts = _sample_ring(main_poly.exterior, boundary_spacing)
|
|
if len(outer_pts) < 3:
|
|
return None
|
|
|
|
start_idx = 0
|
|
vertices.extend(outer_pts.tolist())
|
|
segs = _ring_to_segments(outer_pts, start_idx)
|
|
segments.extend(segs)
|
|
start_idx += len(outer_pts)
|
|
|
|
# Handle any interior rings of the inner plate (shouldn't happen normally)
|
|
for interior_ring in main_poly.interiors:
|
|
ring_pts = _sample_ring(interior_ring, boundary_spacing)
|
|
if len(ring_pts) < 3:
|
|
continue
|
|
vertices.extend(ring_pts.tolist())
|
|
segs = _ring_to_segments(ring_pts, start_idx)
|
|
segments.extend(segs)
|
|
# Hole marker inside this interior ring
|
|
centroid = Polygon(interior_ring).centroid
|
|
hole_markers.append([centroid.x, centroid.y])
|
|
start_idx += len(ring_pts)
|
|
|
|
# Keepout holes — clip to inner plate
|
|
for keepout in keepouts:
|
|
clipped = inner_plate.intersection(keepout)
|
|
if clipped.is_empty:
|
|
continue
|
|
for kpoly in _geometry_to_list(clipped):
|
|
if kpoly.area < 1.0:
|
|
continue
|
|
ring_pts = _sample_ring(kpoly.exterior, boundary_spacing)
|
|
if len(ring_pts) < 3:
|
|
continue
|
|
vertices.extend(ring_pts.tolist())
|
|
segs = _ring_to_segments(ring_pts, start_idx)
|
|
segments.extend(segs)
|
|
# Hole marker inside the keepout
|
|
centroid = kpoly.centroid
|
|
hole_markers.append([centroid.x, centroid.y])
|
|
start_idx += len(ring_pts)
|
|
|
|
return {
|
|
'vertices': np.array(vertices, dtype=np.float64),
|
|
'segments': np.array(segments, dtype=np.int32),
|
|
'holes': np.array(hole_markers, dtype=np.float64) if hole_markers else np.empty((0, 2)),
|
|
}
|
|
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# Step 4 — Compute area constraint from density field
|
|
# ---------------------------------------------------------------------------
|
|
|
|
def _max_area_for_spacing(spacing: float) -> float:
|
|
"""Area of an equilateral triangle with given edge length."""
|
|
return (np.sqrt(3) / 4.0) * spacing ** 2
|
|
|
|
|
|
def _refine_with_density(tri_result, geometry, params):
|
|
"""
|
|
Iterative refinement: check each triangle's area against the density-based
|
|
max area at its centroid. If too large, re-triangulate with tighter constraint.
|
|
|
|
Uses Triangle's region-based area constraints via the 'a' switch.
|
|
"""
|
|
s_min = float(params['s_min'])
|
|
s_max = float(params['s_max'])
|
|
|
|
verts = tri_result['vertices']
|
|
tris = tri_result['triangles']
|
|
|
|
# Compute max allowed area for each triangle based on density at centroid
|
|
max_areas = np.empty(len(tris))
|
|
for i, t in enumerate(tris):
|
|
cx = np.mean(verts[t, 0])
|
|
cy = np.mean(verts[t, 1])
|
|
eta = evaluate_density(cx, cy, geometry, params)
|
|
target_spacing = density_to_spacing(eta, params)
|
|
max_areas[i] = _max_area_for_spacing(target_spacing)
|
|
|
|
# Use the global maximum area as a single constraint
|
|
# (Triangle library doesn't support per-triangle area easily in basic mode)
|
|
# Instead, use a moderate constraint and let the quality refinement handle it
|
|
global_max_area = _max_area_for_spacing(s_max)
|
|
return global_max_area
|
|
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# Main entry point
|
|
# ---------------------------------------------------------------------------
|
|
|
|
def generate_triangulation(geometry, params, max_refinement_passes=3):
|
|
"""
|
|
Generate isogrid triangulation using Triangle library.
|
|
|
|
Returns dict with 'vertices' (Nx2) and 'triangles' (Mx3).
|
|
"""
|
|
s_min = float(params['s_min'])
|
|
s_max = float(params['s_max'])
|
|
|
|
# Step 1: Build inner plate
|
|
inner_plate = _build_inner_plate(geometry, params)
|
|
if inner_plate.is_empty:
|
|
return {'vertices': np.empty((0, 2)), 'triangles': np.empty((0, 3), dtype=int)}
|
|
|
|
# Step 2: Build keepouts
|
|
keepouts = _build_keepouts(geometry, params)
|
|
keepout_union = unary_union(keepouts) if keepouts else Polygon()
|
|
|
|
# Step 3: Build PSLG
|
|
# _sample_ring now uses actual polygon vertices (preserving tight features)
|
|
# and only adds interpolated points on long straight edges.
|
|
boundary_spacing = max(s_min, min(s_max * 0.4, 25.0))
|
|
pslg = _build_pslg(inner_plate, keepouts, boundary_spacing)
|
|
|
|
if pslg is None or len(pslg['vertices']) < 3:
|
|
return {'vertices': np.empty((0, 2)), 'triangles': np.empty((0, 3), dtype=int)}
|
|
|
|
# Step 4: Initial triangulation with coarse area constraint
|
|
# 'p' = use PSLG (segments + holes)
|
|
# 'q25' = minimum angle 25° (no slivers)
|
|
# 'a{area}' = max triangle area (coarse — s_max)
|
|
max_area = _max_area_for_spacing(s_max)
|
|
tri_switches = f'pq25a{max_area:.1f}'
|
|
|
|
tri_input = {
|
|
'vertices': pslg['vertices'],
|
|
'segments': pslg['segments'],
|
|
}
|
|
if len(pslg['holes']) > 0:
|
|
tri_input['holes'] = pslg['holes']
|
|
|
|
try:
|
|
result = tr.triangulate(tri_input, tri_switches)
|
|
except Exception:
|
|
return {'vertices': np.empty((0, 2)), 'triangles': np.empty((0, 3), dtype=int)}
|
|
|
|
verts = result.get('vertices', np.empty((0, 2)))
|
|
tris = result.get('triangles', np.empty((0, 3), dtype=int))
|
|
|
|
if len(tris) == 0:
|
|
return {'vertices': verts, 'triangles': tris}
|
|
|
|
# Step 5: Density-based refinement using per-triangle area constraints
|
|
# Triangle's 'r' switch refines an existing triangulation.
|
|
# We set 'triangle_max_area' per triangle based on density at centroid.
|
|
for pass_num in range(max_refinement_passes):
|
|
# Compute per-triangle max area from density field
|
|
per_tri_max = np.empty(len(tris))
|
|
needs_refinement = False
|
|
|
|
for i, t in enumerate(tris):
|
|
cx = np.mean(verts[t, 0])
|
|
cy = np.mean(verts[t, 1])
|
|
eta = evaluate_density(cx, cy, geometry, params)
|
|
target_spacing = density_to_spacing(eta, params)
|
|
per_tri_max[i] = _max_area_for_spacing(target_spacing)
|
|
|
|
# Check if this triangle actually needs refinement
|
|
v0, v1, v2 = verts[t[0]], verts[t[1]], verts[t[2]]
|
|
actual_area = 0.5 * abs(
|
|
(v1[0] - v0[0]) * (v2[1] - v0[1]) -
|
|
(v2[0] - v0[0]) * (v1[1] - v0[1])
|
|
)
|
|
if actual_area > per_tri_max[i] * 1.3:
|
|
needs_refinement = True
|
|
|
|
if not needs_refinement:
|
|
break
|
|
|
|
# Build refinement input with per-triangle area constraints
|
|
refine_input = {
|
|
'vertices': verts,
|
|
'triangles': tris,
|
|
'segments': pslg['segments'] if len(pslg['segments']) <= len(verts) else np.empty((0, 2), dtype=np.int32),
|
|
'triangle_max_area': per_tri_max,
|
|
}
|
|
if len(pslg['holes']) > 0:
|
|
refine_input['holes'] = pslg['holes']
|
|
|
|
try:
|
|
result = tr.triangulate(refine_input, 'rpq25')
|
|
verts = result.get('vertices', verts)
|
|
tris = result.get('triangles', tris)
|
|
except Exception:
|
|
break
|
|
|
|
# Step 6: Filter degenerate / tiny triangles
|
|
min_area_filter = float(params.get('min_triangle_area', 20.0))
|
|
if len(tris) > 0:
|
|
v0 = verts[tris[:, 0]]
|
|
v1 = verts[tris[:, 1]]
|
|
v2 = verts[tris[:, 2]]
|
|
areas = 0.5 * np.abs(
|
|
(v1[:, 0] - v0[:, 0]) * (v2[:, 1] - v0[:, 1]) -
|
|
(v2[:, 0] - v0[:, 0]) * (v1[:, 1] - v0[:, 1])
|
|
)
|
|
tris = tris[areas >= min_area_filter]
|
|
|
|
# Step 7: Snap out-of-bounds vertices to nearest boundary point
|
|
# Only snap vertices that are clearly outside (> 0.1mm), not boundary vertices
|
|
snap_tol = 0.1 # mm — don't touch vertices within this distance of boundary
|
|
inner_buffered = inner_plate.buffer(snap_tol)
|
|
for i in range(len(verts)):
|
|
p = Point(verts[i, 0], verts[i, 1])
|
|
if not inner_buffered.contains(p):
|
|
nearest = inner_plate.exterior.interpolate(
|
|
inner_plate.exterior.project(p)
|
|
)
|
|
verts[i, 0] = nearest.x
|
|
verts[i, 1] = nearest.y
|
|
|
|
return {'vertices': verts, 'triangles': tris, 'inner_plate': inner_plate}
|