refactor: rewrite triangulation using Triangle library (constrained Delaunay + quality refinement)
- Replace scipy.spatial.Delaunay with Shewchuk's Triangle (PSLG-based) - Boundary conforming: PSLG constrains edges along inset contour + hole keepout rings - Quality: min angle 25°, no slivers - Per-triangle density-based area refinement (s_min=20, s_max=80) - Clean boundary plotting (no more crooked v1 line resampling) - Triangulation plot shows inset contour (red dashed) + keepout rings (orange dashed) - Add sandbox2_brain_input.json geometry file
This commit is contained in:
@@ -37,8 +37,8 @@ DEFAULT_PARAMS = {
|
||||
'p': 2.0,
|
||||
'beta': 0.3,
|
||||
'R_edge': 15.0,
|
||||
's_min': 45.0,
|
||||
's_max': 55.0,
|
||||
's_min': 20.0,
|
||||
's_max': 80.0,
|
||||
't_min': 2.5,
|
||||
't_0': 3.0,
|
||||
'gamma': 1.0,
|
||||
|
||||
@@ -43,44 +43,21 @@ def _merge_params(geometry: Dict[str, Any], params_file: Path | None) -> Dict[st
|
||||
|
||||
|
||||
def _plot_boundary_polyline(geometry: Dict[str, Any], arc_pts: int = 64) -> np.ndarray:
|
||||
"""Get the sandbox boundary as a clean polyline for plotting.
|
||||
|
||||
For v2 typed segments: densify arcs, keep lines exact.
|
||||
For v1: just use the raw outer_boundary vertices — no resampling needed
|
||||
since these are the true polygon corners from NX.
|
||||
"""
|
||||
typed = geometry.get("outer_boundary_typed")
|
||||
if typed:
|
||||
pts = typed_segments_to_polyline(typed, arc_pts=arc_pts)
|
||||
if len(pts) >= 3:
|
||||
return np.asarray(pts, dtype=float)
|
||||
|
||||
# v1 fallback: use raw boundary vertices directly — they ARE the geometry.
|
||||
outer = np.asarray(geometry["outer_boundary"], dtype=float)
|
||||
if len(outer) < 4:
|
||||
return outer
|
||||
# v1 fallback: dense polyline boundaries may encode fillets.
|
||||
# Resample uniformly for smoother plotting while preserving the polygon path.
|
||||
if np.allclose(outer[0], outer[-1]):
|
||||
ring = outer
|
||||
else:
|
||||
ring = np.vstack([outer, outer[0]])
|
||||
seg = np.diff(ring, axis=0)
|
||||
seg_len = np.hypot(seg[:, 0], seg[:, 1])
|
||||
nonzero = seg_len[seg_len > 1e-9]
|
||||
if len(nonzero) == 0:
|
||||
return outer
|
||||
step = max(float(np.median(nonzero) * 0.5), 0.5)
|
||||
cum = np.r_[0.0, np.cumsum(seg_len)]
|
||||
total = float(cum[-1])
|
||||
if total <= step:
|
||||
return outer
|
||||
samples = np.arange(0.0, total, step, dtype=float)
|
||||
if samples[-1] < total:
|
||||
samples = np.r_[samples, total]
|
||||
out = []
|
||||
j = 0
|
||||
for s in samples:
|
||||
while j < len(cum) - 2 and cum[j + 1] < s:
|
||||
j += 1
|
||||
den = max(cum[j + 1] - cum[j], 1e-12)
|
||||
t = (s - cum[j]) / den
|
||||
p = ring[j] + t * (ring[j + 1] - ring[j])
|
||||
out.append([float(p[0]), float(p[1])])
|
||||
return np.asarray(out, dtype=float)
|
||||
return outer
|
||||
|
||||
|
||||
def _plot_density(geometry: Dict[str, Any], params: Dict[str, Any], out_path: Path, resolution: float = 3.0) -> None:
|
||||
@@ -114,65 +91,83 @@ def _plot_density(geometry: Dict[str, Any], params: Dict[str, Any], out_path: Pa
|
||||
|
||||
def _plot_triangulation(geometry: Dict[str, Any], triangulation: Dict[str, Any], out_path: Path, params: Dict[str, Any] = None) -> None:
|
||||
"""
|
||||
Plot the Delaunay triangulation clipped to the inner frame (w_frame inset).
|
||||
Green boundary, blue rib edges, blue hole circles.
|
||||
Plot the triangulation with:
|
||||
- Green: original sandbox boundary
|
||||
- Red dashed: inset contour (w_frame offset) — this is what the mesh fills
|
||||
- Blue: triangle edges
|
||||
- Blue circles: hole boundaries
|
||||
- Orange dashed: hole keepout rings
|
||||
"""
|
||||
from shapely.geometry import Polygon as ShapelyPolygon, LineString, Point
|
||||
from shapely.geometry import Polygon as ShapelyPolygon, Point as ShapelyPoint
|
||||
|
||||
verts = triangulation["vertices"]
|
||||
tris = triangulation["triangles"]
|
||||
|
||||
outer = _plot_boundary_polyline(geometry)
|
||||
plate_poly = ShapelyPolygon(outer)
|
||||
w_frame = (params or {}).get("w_frame", 2.0)
|
||||
w_frame = (params or {}).get("w_frame", 8.0)
|
||||
d_keep = (params or {}).get("d_keep", 1.5)
|
||||
inner_plate = plate_poly.buffer(-w_frame)
|
||||
if inner_plate.is_empty or not inner_plate.is_valid:
|
||||
inner_plate = plate_poly
|
||||
|
||||
fig, ax = plt.subplots(figsize=(8, 6), dpi=160)
|
||||
fig, ax = plt.subplots(figsize=(10, 8), dpi=160)
|
||||
|
||||
# Draw triangle edges clipped to inner plate
|
||||
# Draw all triangle edges
|
||||
drawn_edges = set()
|
||||
for tri in tris:
|
||||
centroid = Point(verts[tri].mean(axis=0))
|
||||
if not inner_plate.contains(centroid):
|
||||
continue
|
||||
for i in range(3):
|
||||
edge = tuple(sorted([tri[i], tri[(i + 1) % 3]]))
|
||||
if edge in drawn_edges:
|
||||
continue
|
||||
drawn_edges.add(edge)
|
||||
p1, p2 = verts[edge[0]], verts[edge[1]]
|
||||
line = LineString([p1, p2])
|
||||
clipped = inner_plate.intersection(line)
|
||||
if clipped.is_empty:
|
||||
continue
|
||||
if clipped.geom_type == "LineString":
|
||||
cx, cy = clipped.xy
|
||||
ax.plot(cx, cy, color="#1f77b4", lw=0.5, alpha=0.85)
|
||||
elif clipped.geom_type == "MultiLineString":
|
||||
for seg in clipped.geoms:
|
||||
cx, cy = seg.xy
|
||||
ax.plot(cx, cy, color="#1f77b4", lw=0.5, alpha=0.85)
|
||||
ax.plot([p1[0], p2[0]], [p1[1], p2[1]],
|
||||
color="#1f77b4", lw=0.5, alpha=0.85)
|
||||
|
||||
# Green sandbox boundary
|
||||
# Green: original sandbox boundary
|
||||
ax.plot(np.r_[outer[:, 0], outer[0, 0]], np.r_[outer[:, 1], outer[0, 1]],
|
||||
color="#228833", lw=1.8, zorder=5)
|
||||
color="#228833", lw=1.8, zorder=5, label="Sandbox boundary")
|
||||
|
||||
# Blue hole circles
|
||||
# Red dashed: inset contour (w_frame)
|
||||
if not inner_plate.is_empty:
|
||||
polys = [inner_plate] if inner_plate.geom_type == 'Polygon' else list(inner_plate.geoms)
|
||||
for poly in polys:
|
||||
ix, iy = poly.exterior.xy
|
||||
ax.plot(ix, iy, color="#cc3333", lw=1.2, ls="--", zorder=4,
|
||||
label=f"Inset contour (w_frame={w_frame}mm)")
|
||||
|
||||
# Blue hole circles + Orange keepout rings
|
||||
for hole in geometry.get("holes", []):
|
||||
if hole.get("is_circular", False) and "center" in hole and "diameter" in hole:
|
||||
cx, cy = hole["center"]
|
||||
r = float(hole["diameter"]) * 0.5
|
||||
a = np.linspace(0.0, 2.0 * np.pi, 90, endpoint=True)
|
||||
# Hole circle
|
||||
hb = np.column_stack([cx + r * np.cos(a), cy + r * np.sin(a)])
|
||||
ax.plot(np.r_[hb[:, 0], hb[0, 0]], np.r_[hb[:, 1], hb[0, 1]],
|
||||
"b-", lw=0.8, zorder=3)
|
||||
# Keepout ring
|
||||
kr = r * (1.0 + d_keep)
|
||||
kb = np.column_stack([cx + kr * np.cos(a), cy + kr * np.sin(a)])
|
||||
ax.plot(np.r_[kb[:, 0], kb[0, 0]], np.r_[kb[:, 1], kb[0, 1]],
|
||||
color="#ff8800", lw=0.8, ls="--", zorder=3)
|
||||
else:
|
||||
hb = np.asarray(hole["boundary"])
|
||||
ax.plot(np.r_[hb[:, 0], hb[0, 0]], np.r_[hb[:, 1], hb[0, 1]],
|
||||
"b-", lw=0.8, zorder=3)
|
||||
ax.plot(np.r_[hb[:, 0], hb[0, 0]], np.r_[hb[:, 1], hb[0, 1]],
|
||||
"b-", lw=0.8, zorder=3)
|
||||
|
||||
# De-duplicate legend entries
|
||||
handles, labels = ax.get_legend_handles_labels()
|
||||
seen = set()
|
||||
unique_handles, unique_labels = [], []
|
||||
for h, l in zip(handles, labels):
|
||||
if l not in seen:
|
||||
seen.add(l)
|
||||
unique_handles.append(h)
|
||||
unique_labels.append(l)
|
||||
ax.legend(unique_handles, unique_labels, loc="upper right", fontsize=8)
|
||||
|
||||
ax.set_aspect("equal", adjustable="box")
|
||||
ax.set_title("Constrained Delaunay Triangulation / Rib Pattern")
|
||||
ax.set_title(f"Isogrid Triangulation ({len(tris)} triangles, {len(drawn_edges)} edges)")
|
||||
ax.set_xlabel("x [mm]")
|
||||
ax.set_ylabel("y [mm]")
|
||||
fig.tight_layout()
|
||||
|
||||
@@ -1,397 +1,365 @@
|
||||
"""
|
||||
Isogrid triangulation — generates a proper isogrid rib pattern.
|
||||
|
||||
Strategy: lay a regular equilateral triangle grid, then adapt it:
|
||||
1. Generate hex-packed vertex grid at base spacing
|
||||
2. Scale spacing locally based on density field (growing away from features)
|
||||
3. Delaunay triangulate the point set
|
||||
4. Clip to plate boundary, exclude hole keepouts
|
||||
5. Result: well-proportioned triangles everywhere, denser near holes/edges
|
||||
Uses Shewchuk's Triangle library (constrained Delaunay + quality refinement)
|
||||
for boundary-conforming, clean triangulations.
|
||||
|
||||
This is NOT a mesh — it's an isogrid layout. Every triangle should be
|
||||
a viable machining pocket.
|
||||
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
|
||||
from scipy.spatial import Delaunay
|
||||
from shapely.geometry import Polygon, Point, LinearRing
|
||||
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, typed_segments_to_sparse
|
||||
from src.shared.arc_utils import inset_arc, typed_segments_to_polyline
|
||||
from .density_field import evaluate_density, density_to_spacing
|
||||
|
||||
|
||||
def _geometry_to_list(geom: BaseGeometry) -> list[BaseGeometry]:
|
||||
# ---------------------------------------------------------------------------
|
||||
# Helpers
|
||||
# ---------------------------------------------------------------------------
|
||||
|
||||
def _geometry_to_list(geom: BaseGeometry) -> list:
|
||||
if geom.is_empty:
|
||||
return []
|
||||
return [geom] if geom.geom_type == 'Polygon' else list(getattr(geom, 'geoms', []))
|
||||
if geom.geom_type == 'Polygon':
|
||||
return [geom]
|
||||
if geom.geom_type == 'MultiPolygon':
|
||||
return list(geom.geoms)
|
||||
return list(getattr(geom, 'geoms', []))
|
||||
|
||||
|
||||
def _add_boundary_layer_seed_points(points, geometry, params, inner_plate, keepout_union):
|
||||
"""Add a row of points just inside the inset boundary to align boundary triangles."""
|
||||
offset = float(params.get('boundary_layer_offset', 0.0) or 0.0)
|
||||
if offset <= 0.0:
|
||||
offset = max(params.get('s_min', 10.0) * 0.6, 0.5)
|
||||
|
||||
layer_geom = inner_plate.buffer(-offset)
|
||||
if layer_geom.is_empty:
|
||||
return points
|
||||
|
||||
boundary_pts = []
|
||||
for poly in _geometry_to_list(layer_geom):
|
||||
ring = poly.exterior
|
||||
length = float(ring.length)
|
||||
if length <= 1e-9:
|
||||
continue
|
||||
n_pts = max(int(np.ceil(length / max(params.get('s_min', 10.0), 1e-3))), 8)
|
||||
for i in range(n_pts):
|
||||
frac = i / n_pts
|
||||
p = ring.interpolate(frac, normalized=True)
|
||||
if not inner_plate.buffer(1e-6).covers(p):
|
||||
continue
|
||||
if not keepout_union.is_empty and keepout_union.buffer(1e-6).covers(p):
|
||||
continue
|
||||
boundary_pts.append([p.x, p.y])
|
||||
|
||||
if boundary_pts:
|
||||
return np.vstack([points, np.asarray(boundary_pts, dtype=np.float64)])
|
||||
return points
|
||||
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 _np(pt):
|
||||
return np.asarray([float(pt[0]), float(pt[1])], dtype=float)
|
||||
|
||||
|
||||
def _hole_keepout_seed_points(geometry, params):
|
||||
"""Return sparse keepout seed points: 3 for circular holes, coarse ring for others."""
|
||||
d_keep = float(params['d_keep'])
|
||||
pts = []
|
||||
for hole in geometry.get('holes', []):
|
||||
if hole.get('is_circular', False) and 'center' in hole:
|
||||
cx, cy = hole['center']
|
||||
diameter = float(hole.get('diameter', 10.0) or 10.0)
|
||||
hole_radius = diameter / 2.0
|
||||
keepout_radius = hole_radius * (1.0 + d_keep)
|
||||
for k in range(3):
|
||||
a = (2.0 * np.pi * k) / 3.0
|
||||
pts.append([cx + keepout_radius * np.cos(a), cy + keepout_radius * np.sin(a)])
|
||||
continue
|
||||
|
||||
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 keepout.is_empty:
|
||||
continue
|
||||
for geom in _geometry_to_list(keepout):
|
||||
ext = geom.exterior
|
||||
n = max(int(np.ceil(ext.length / max(params.get('s_min', 10.0), 1e-3))), 6)
|
||||
for i in range(n):
|
||||
p = ext.interpolate(i / n, normalized=True)
|
||||
pts.append([p.x, p.y])
|
||||
|
||||
if not pts:
|
||||
def _sample_ring(ring, spacing: float) -> np.ndarray:
|
||||
"""Sample points along a Shapely ring at given spacing, returning Nx2 array (not closed)."""
|
||||
length = float(ring.length)
|
||||
if length < 1e-9:
|
||||
return np.empty((0, 2), dtype=np.float64)
|
||||
return np.asarray(pts, 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)
|
||||
|
||||
|
||||
def _generate_hex_grid(bbox, base_spacing):
|
||||
"""
|
||||
Generate a regular hexagonal-packed point grid.
|
||||
|
||||
This creates the classic isogrid vertex layout: rows of points
|
||||
offset by half-spacing, producing equilateral triangles when
|
||||
Delaunay triangulated.
|
||||
|
||||
Returns ndarray of shape (N, 2).
|
||||
"""
|
||||
x_min, y_min, x_max, y_max = bbox
|
||||
# Padding
|
||||
pad = base_spacing
|
||||
x_min -= pad
|
||||
y_min -= pad
|
||||
x_max += pad
|
||||
y_max += pad
|
||||
# ---------------------------------------------------------------------------
|
||||
# Step 1 — Build the inset plate polygon (frame inner edge)
|
||||
# ---------------------------------------------------------------------------
|
||||
|
||||
row_height = base_spacing * np.sqrt(3) / 2.0
|
||||
points = []
|
||||
row = 0
|
||||
y = y_min
|
||||
while y <= y_max:
|
||||
# Alternate rows offset by half spacing
|
||||
x_offset = (base_spacing / 2.0) if (row % 2) else 0.0
|
||||
x = x_min + x_offset
|
||||
while x <= x_max:
|
||||
points.append([x, y])
|
||||
x += base_spacing
|
||||
y += row_height
|
||||
row += 1
|
||||
|
||||
return np.array(points, dtype=np.float64)
|
||||
|
||||
|
||||
def _adaptive_hex_grid(geometry, params):
|
||||
"""
|
||||
Generate density-adaptive hex grid.
|
||||
|
||||
Start with a coarse grid at s_max spacing, then iteratively add
|
||||
points in high-density regions until local spacing matches the
|
||||
density-driven target.
|
||||
"""
|
||||
outer = np.array(geometry['outer_boundary'])
|
||||
x_min, y_min = outer.min(axis=0)
|
||||
x_max, y_max = outer.max(axis=0)
|
||||
bbox = (x_min, y_min, x_max, y_max)
|
||||
|
||||
s_min = params['s_min']
|
||||
s_max = params['s_max']
|
||||
plate_poly = Polygon(geometry['outer_boundary'])
|
||||
|
||||
# Build hole keepout polygons
|
||||
d_keep = params['d_keep']
|
||||
keepout_polys = []
|
||||
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
|
||||
boss_radius = hole_radius + d_keep * hole_radius
|
||||
keepout = Point(cx, cy).buffer(boss_radius, resolution=16)
|
||||
keepout_polys.append(keepout)
|
||||
|
||||
keepout_union = unary_union(keepout_polys) if keepout_polys else Polygon()
|
||||
|
||||
# Valid region = plate minus keepouts
|
||||
valid_region = plate_poly.difference(keepout_union)
|
||||
|
||||
# Inset valid region by frame width to keep grid points inside the frame
|
||||
w_frame = params.get('w_frame', 8.0)
|
||||
inner_region = plate_poly.buffer(-w_frame)
|
||||
if inner_region.is_empty:
|
||||
inner_region = plate_poly
|
||||
|
||||
# Generate base grid at s_max (coarsest)
|
||||
base_pts = _generate_hex_grid(bbox, s_max)
|
||||
|
||||
# Filter to plate interior (outside keepouts, inside frame)
|
||||
valid_pts = []
|
||||
for pt in base_pts:
|
||||
p = Point(pt[0], pt[1])
|
||||
if inner_region.contains(p) and not keepout_union.contains(p):
|
||||
valid_pts.append(pt)
|
||||
|
||||
valid_pts = np.array(valid_pts, dtype=np.float64) if valid_pts else np.empty((0, 2))
|
||||
|
||||
if s_min >= s_max * 0.95:
|
||||
# Uniform mode — no refinement needed
|
||||
return valid_pts, valid_region, keepout_union
|
||||
|
||||
# Adaptive refinement: add denser points near high-density areas
|
||||
# Generate a fine grid at s_min
|
||||
fine_pts = _generate_hex_grid(bbox, s_min)
|
||||
|
||||
# For each fine point, check if local density warrants adding it
|
||||
extra_pts = []
|
||||
for pt in fine_pts:
|
||||
p = Point(pt[0], pt[1])
|
||||
if not inner_region.contains(p) or keepout_union.contains(p):
|
||||
continue
|
||||
|
||||
# What spacing does the density field want here?
|
||||
eta = evaluate_density(pt[0], pt[1], geometry, params)
|
||||
target_spacing = density_to_spacing(eta, params)
|
||||
|
||||
# If target spacing < s_max * 0.8, this is a refinement zone
|
||||
if target_spacing < s_max * 0.8:
|
||||
# Check distance to nearest existing point
|
||||
if len(valid_pts) > 0:
|
||||
dists = np.sqrt(np.sum((valid_pts - pt)**2, axis=1))
|
||||
min_dist = np.min(dists)
|
||||
# Add point if nearest neighbor is farther than target spacing
|
||||
if min_dist > target_spacing * 0.85:
|
||||
extra_pts.append(pt)
|
||||
# Also add to valid_pts for future distance checks
|
||||
valid_pts = np.vstack([valid_pts, pt])
|
||||
|
||||
if extra_pts:
|
||||
all_pts = np.vstack([valid_pts, np.array(extra_pts)])
|
||||
# Deduplicate (shouldn't be needed but safe)
|
||||
all_pts = np.unique(np.round(all_pts, 6), axis=0)
|
||||
else:
|
||||
all_pts = valid_pts
|
||||
|
||||
return all_pts, valid_region, keepout_union
|
||||
|
||||
|
||||
def _add_boundary_vertices(points, geometry, params, keepout_union):
|
||||
"""Add inset boundary vertices (arc-aware) and keepout boundary points."""
|
||||
s_min = params['s_min']
|
||||
w_frame = params.get('w_frame', 8.0)
|
||||
|
||||
new_pts = list(points)
|
||||
def _build_inner_plate(geometry, params) -> Polygon:
|
||||
"""Offset sandbox boundary inward by w_frame."""
|
||||
w_frame = float(params.get('w_frame', 8.0))
|
||||
plate_poly = Polygon(geometry['outer_boundary'])
|
||||
|
||||
typed_segments = geometry.get('outer_boundary_typed')
|
||||
if typed_segments:
|
||||
ring = LinearRing(geometry['outer_boundary'])
|
||||
is_ccw = bool(ring.is_ccw)
|
||||
inset_segments = []
|
||||
|
||||
inset_segments = []
|
||||
for seg in typed_segments:
|
||||
stype = seg.get('type', 'line')
|
||||
if stype == 'arc':
|
||||
center_inside = plate_poly.contains(Point(seg['center']))
|
||||
inset_segments.append(inset_arc({**seg, 'center_inside': center_inside}, w_frame))
|
||||
continue
|
||||
else:
|
||||
x1, y1 = seg['start']
|
||||
x2, y2 = seg['end']
|
||||
dx, dy = (x2 - x1), (y2 - y1)
|
||||
ln = np.hypot(dx, dy)
|
||||
if ln < 1e-12:
|
||||
continue
|
||||
nx_l, ny_l = (-dy / ln), (dx / ln)
|
||||
nx, ny = (nx_l, ny_l) if is_ccw else (-nx_l, -ny_l)
|
||||
inset_segments.append({
|
||||
'type': 'line',
|
||||
'start': [x1 + w_frame * nx, y1 + w_frame * ny],
|
||||
'end': [x2 + w_frame * nx, y2 + w_frame * ny],
|
||||
})
|
||||
|
||||
x1, y1 = seg['start']
|
||||
x2, y2 = seg['end']
|
||||
dx, dy = (x2 - x1), (y2 - y1)
|
||||
ln = np.hypot(dx, dy)
|
||||
if ln < 1e-12:
|
||||
continue
|
||||
nx_l, ny_l = (-dy / ln), (dx / ln)
|
||||
nx, ny = (nx_l, ny_l) if is_ccw else (-nx_l, -ny_l)
|
||||
inset_segments.append({
|
||||
'type': 'line',
|
||||
'start': [x1 + w_frame * nx, y1 + w_frame * ny],
|
||||
'end': [x2 + w_frame * nx, y2 + w_frame * ny],
|
||||
})
|
||||
|
||||
# Sparse points forced into Delaunay (3/arc, 2/line)
|
||||
sparse_pts = typed_segments_to_sparse(inset_segments)
|
||||
new_pts.extend(sparse_pts)
|
||||
|
||||
# Dense inset polygon for containment checks
|
||||
dense = typed_segments_to_polyline(inset_segments, arc_pts=16)
|
||||
dense = typed_segments_to_polyline(inset_segments, arc_pts=32)
|
||||
if len(dense) >= 3:
|
||||
inner_plate = Polygon(dense)
|
||||
if not inner_plate.is_valid:
|
||||
inner_plate = inner_plate.buffer(0)
|
||||
if inner_plate.is_empty:
|
||||
inner_plate = plate_poly.buffer(-w_frame)
|
||||
else:
|
||||
inner_plate = plate_poly.buffer(-w_frame)
|
||||
|
||||
# Even spacing along inset boundary (as before)
|
||||
if not inner_plate.is_empty and inner_plate.is_valid:
|
||||
for inner_poly in _geometry_to_list(inner_plate):
|
||||
ext = inner_poly.exterior
|
||||
length = ext.length
|
||||
n_pts = max(int(length / s_min), 4)
|
||||
for i in range(n_pts):
|
||||
frac = i / n_pts
|
||||
pt = ext.interpolate(frac, normalized=True)
|
||||
new_pts.append([pt.x, pt.y])
|
||||
else:
|
||||
# v1 geometry fallback: inset from polygon and force all inset ring vertices.
|
||||
inner_plate = plate_poly.buffer(-w_frame)
|
||||
if not inner_plate.is_empty:
|
||||
for inner_poly in _geometry_to_list(inner_plate):
|
||||
ext = inner_poly.exterior
|
||||
# Always keep true inset corners/vertices from the offset polygon.
|
||||
coords = list(ext.coords)[:-1]
|
||||
for cx, cy in coords:
|
||||
new_pts.append([cx, cy])
|
||||
# Plus evenly spaced boundary points for perimeter alignment.
|
||||
length = ext.length
|
||||
n_pts = max(int(length / s_min), 4)
|
||||
for i in range(n_pts):
|
||||
frac = i / n_pts
|
||||
pt = ext.interpolate(frac, normalized=True)
|
||||
new_pts.append([pt.x, pt.y])
|
||||
|
||||
# Add sparse points for hole keepouts (3 points per circular keepout).
|
||||
keepout_seeds = _hole_keepout_seed_points(geometry, params)
|
||||
if len(keepout_seeds) > 0:
|
||||
new_pts.extend(keepout_seeds.tolist())
|
||||
if not inner_plate.is_empty:
|
||||
return inner_plate
|
||||
|
||||
inner_plate = plate_poly.buffer(-w_frame)
|
||||
if inner_plate.is_empty or not inner_plate.is_valid:
|
||||
inner_plate = plate_poly
|
||||
return inner_plate
|
||||
|
||||
return np.array(new_pts, dtype=np.float64), 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: hex-packed grid + Delaunay + clip.
|
||||
|
||||
Returns dict with 'vertices' (ndarray Nx2) and 'triangles' (ndarray Mx3).
|
||||
Generate isogrid triangulation using Triangle library.
|
||||
|
||||
Returns dict with 'vertices' (Nx2) and 'triangles' (Mx3).
|
||||
"""
|
||||
# Generate adaptive point grid
|
||||
grid_pts, valid_region, keepout_union = _adaptive_hex_grid(geometry, params)
|
||||
s_min = float(params['s_min'])
|
||||
s_max = float(params['s_max'])
|
||||
|
||||
if len(grid_pts) < 3:
|
||||
return {'vertices': grid_pts, 'triangles': np.empty((0, 3), dtype=int)}
|
||||
|
||||
# Add boundary-conforming vertices and get inset plate polygon for clipping
|
||||
all_pts, inner_plate = _add_boundary_vertices(grid_pts, geometry, params, keepout_union)
|
||||
|
||||
# Add structured boundary-layer seed row along straight edges
|
||||
plate_poly = Polygon(geometry['outer_boundary'])
|
||||
all_pts = _add_boundary_layer_seed_points(all_pts, geometry, params, inner_plate, keepout_union)
|
||||
|
||||
# Deduplicate close points
|
||||
all_pts = np.unique(np.round(all_pts, 4), axis=0)
|
||||
|
||||
if len(all_pts) < 3:
|
||||
return {'vertices': all_pts, 'triangles': np.empty((0, 3), dtype=int)}
|
||||
|
||||
# Delaunay triangulation of the point set
|
||||
tri = Delaunay(all_pts)
|
||||
triangles = tri.simplices
|
||||
|
||||
# Filter: keep only triangles whose centroid is inside valid region
|
||||
plate_poly = Polygon(geometry['outer_boundary'])
|
||||
# Step 1: Build inner plate
|
||||
inner_plate = _build_inner_plate(geometry, params)
|
||||
if inner_plate.is_empty:
|
||||
inner_plate = plate_poly
|
||||
return {'vertices': np.empty((0, 2)), 'triangles': np.empty((0, 3), dtype=int)}
|
||||
|
||||
inner_valid_region = inner_plate
|
||||
if not keepout_union.is_empty:
|
||||
inner_valid_region = inner_plate.difference(keepout_union)
|
||||
# Step 2: Build keepouts
|
||||
keepouts = _build_keepouts(geometry, params)
|
||||
keepout_union = unary_union(keepouts) if keepouts else Polygon()
|
||||
|
||||
keep = []
|
||||
for i, t in enumerate(triangles):
|
||||
cx = np.mean(all_pts[t, 0])
|
||||
cy = np.mean(all_pts[t, 1])
|
||||
centroid = Point(cx, cy)
|
||||
# Keep if centroid is inside plate frame and outside keepouts,
|
||||
# AND all 3 vertices are inside the plate boundary (no crossovers)
|
||||
if not inner_plate.contains(centroid):
|
||||
continue
|
||||
if keepout_union.contains(centroid):
|
||||
continue
|
||||
all_inside = True
|
||||
for vi in t:
|
||||
if not plate_poly.contains(Point(all_pts[vi])):
|
||||
# Allow small tolerance (vertex on boundary is OK)
|
||||
if not plate_poly.buffer(0.5).contains(Point(all_pts[vi])):
|
||||
all_inside = False
|
||||
break
|
||||
if not all_inside:
|
||||
continue
|
||||
# Step 3: Build PSLG
|
||||
# Boundary sampling at intermediate spacing for clean boundary conformance
|
||||
boundary_spacing = max(s_min, min(s_max * 0.5, 30.0))
|
||||
pslg = _build_pslg(inner_plate, keepouts, boundary_spacing)
|
||||
|
||||
tri_poly = Polygon(all_pts[t])
|
||||
if tri_poly.is_empty or not tri_poly.is_valid:
|
||||
continue
|
||||
if not inner_valid_region.buffer(1e-6).covers(tri_poly):
|
||||
continue
|
||||
keep.append(i)
|
||||
if pslg is None or len(pslg['vertices']) < 3:
|
||||
return {'vertices': np.empty((0, 2)), 'triangles': np.empty((0, 3), dtype=int)}
|
||||
|
||||
triangles = triangles[keep] if keep else 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}'
|
||||
|
||||
# Filter degenerate / tiny triangles
|
||||
min_area = params.get('min_triangle_area', 20.0)
|
||||
if len(triangles) > 0:
|
||||
v0 = all_pts[triangles[:, 0]]
|
||||
v1 = all_pts[triangles[:, 1]]
|
||||
v2 = all_pts[triangles[:, 2]]
|
||||
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])
|
||||
)
|
||||
triangles = triangles[areas >= min_area]
|
||||
tris = tris[areas >= min_area_filter]
|
||||
|
||||
return {'vertices': all_pts, 'triangles': triangles}
|
||||
return {'vertices': verts, 'triangles': tris}
|
||||
|
||||
100
tools/adaptive-isogrid/tests/sandbox2_brain_input.json
Normal file
100
tools/adaptive-isogrid/tests/sandbox2_brain_input.json
Normal file
@@ -0,0 +1,100 @@
|
||||
{
|
||||
"plate_id": "sandbox_2",
|
||||
"units": "mm",
|
||||
"thickness": 12.7,
|
||||
"outer_boundary": [
|
||||
[
|
||||
0.0,
|
||||
0.0
|
||||
],
|
||||
[
|
||||
7.5,
|
||||
-7.5
|
||||
],
|
||||
[
|
||||
7.5,
|
||||
-22.6
|
||||
],
|
||||
[
|
||||
22.5,
|
||||
-22.6
|
||||
],
|
||||
[
|
||||
22.5,
|
||||
-13.496098
|
||||
],
|
||||
[
|
||||
74.5,
|
||||
-13.496098
|
||||
],
|
||||
[
|
||||
74.5,
|
||||
-22.6
|
||||
],
|
||||
[
|
||||
102.5,
|
||||
-22.6
|
||||
],
|
||||
[
|
||||
102.5,
|
||||
-7.5
|
||||
],
|
||||
[
|
||||
117.5,
|
||||
-7.5
|
||||
],
|
||||
[
|
||||
117.5,
|
||||
-22.6
|
||||
],
|
||||
[
|
||||
140.748693,
|
||||
-22.6
|
||||
],
|
||||
[
|
||||
140.748693,
|
||||
124.4
|
||||
],
|
||||
[
|
||||
117.5,
|
||||
124.4
|
||||
],
|
||||
[
|
||||
117.5,
|
||||
102.5
|
||||
],
|
||||
[
|
||||
102.5,
|
||||
102.5
|
||||
],
|
||||
[
|
||||
102.5,
|
||||
124.4
|
||||
],
|
||||
[
|
||||
7.5,
|
||||
124.4
|
||||
],
|
||||
[
|
||||
7.5,
|
||||
102.5
|
||||
],
|
||||
[
|
||||
0.0,
|
||||
95.0
|
||||
],
|
||||
[
|
||||
-13.5,
|
||||
95.0
|
||||
],
|
||||
[
|
||||
-13.5,
|
||||
0.0
|
||||
],
|
||||
[
|
||||
0.0,
|
||||
0.0
|
||||
]
|
||||
],
|
||||
"holes": []
|
||||
}
|
||||
Reference in New Issue
Block a user