inire/inire/geometry/components.py
2026-03-11 09:37:54 -07:00

582 lines
21 KiB
Python

from __future__ import annotations
from typing import Literal, cast
import numpy
import shapely
from shapely.geometry import Polygon, box
from shapely.ops import unary_union
from .primitives import Port
# Search Grid Snap (5.0 µm default)
SEARCH_GRID_SNAP_UM = 5.0
def snap_search_grid(value: float, snap_size: float = SEARCH_GRID_SNAP_UM) -> float:
"""
Snap a coordinate to the nearest search grid unit.
Args:
value: Value to snap.
snap_size: The grid size to snap to.
Returns:
Snapped value.
"""
if snap_size <= 0:
return value
return round(value / snap_size) * snap_size
class ComponentResult:
"""
The result of a component generation: geometry, final port, and physical length.
"""
__slots__ = ('geometry', 'dilated_geometry', 'proxy_geometry', 'actual_geometry', 'end_port', 'length', 'bounds', 'dilated_bounds', '_t_cache')
geometry: list[Polygon]
""" List of polygons representing the component geometry (could be proxy or arc) """
dilated_geometry: list[Polygon] | None
""" Optional list of pre-dilated polygons for collision optimization """
proxy_geometry: list[Polygon] | None
""" Simplified conservative proxy for tiered collision checks """
actual_geometry: list[Polygon] | None
""" High-fidelity 'actual' geometry for visualization (always the arc) """
end_port: Port
""" The final port after the component """
length: float
""" Physical length of the component path """
bounds: numpy.ndarray
""" Pre-calculated bounds for each polygon in geometry """
dilated_bounds: numpy.ndarray | None
""" Pre-calculated bounds for each polygon in dilated_geometry """
_t_cache: dict[tuple[float, float], ComponentResult]
""" Cache for translated versions of this result """
def __init__(
self,
geometry: list[Polygon],
end_port: Port,
length: float,
dilated_geometry: list[Polygon] | None = None,
proxy_geometry: list[Polygon] | None = None,
actual_geometry: list[Polygon] | None = None,
skip_bounds: bool = False,
) -> None:
self.geometry = geometry
self.dilated_geometry = dilated_geometry
self.proxy_geometry = proxy_geometry
self.actual_geometry = actual_geometry
self.end_port = end_port
self.length = length
self._t_cache = {}
if not skip_bounds:
# Vectorized bounds calculation
self.bounds = shapely.bounds(geometry)
self.dilated_bounds = shapely.bounds(dilated_geometry) if dilated_geometry is not None else None
def translate(self, dx: float, dy: float) -> ComponentResult:
"""
Create a new ComponentResult translated by (dx, dy).
"""
dxr, dyr = round(dx, 3), round(dy, 3)
if (dxr, dyr) == (0.0, 0.0):
return self
if (dxr, dyr) in self._t_cache:
return self._t_cache[(dxr, dyr)]
# Vectorized translation if possible, else list comp
geoms = list(self.geometry)
num_geom = len(self.geometry)
offsets = [num_geom]
if self.dilated_geometry is not None:
geoms.extend(self.dilated_geometry)
offsets.append(len(geoms))
if self.proxy_geometry is not None:
geoms.extend(self.proxy_geometry)
offsets.append(len(geoms))
if self.actual_geometry is not None:
geoms.extend(self.actual_geometry)
offsets.append(len(geoms))
from shapely.affinity import translate
translated = [translate(p, dx, dy) for p in geoms]
new_geom = translated[:offsets[0]]
new_dil = translated[offsets[0]:offsets[1]] if self.dilated_geometry is not None else None
new_proxy = translated[offsets[1]:offsets[2]] if self.proxy_geometry is not None else None
new_actual = translated[offsets[2]:offsets[3]] if self.actual_geometry is not None else None
new_port = Port(self.end_port.x + dx, self.end_port.y + dy, self.end_port.orientation)
res = ComponentResult(new_geom, new_port, self.length, new_dil, new_proxy, new_actual, skip_bounds=True)
# Optimize: reuse and translate bounds
res.bounds = self.bounds + [dx, dy, dx, dy]
if self.dilated_bounds is not None:
res.dilated_bounds = self.dilated_bounds + [dx, dy, dx, dy]
self._t_cache[(dxr, dyr)] = res
return res
class Straight:
"""
Move generator for straight waveguide segments.
"""
@staticmethod
def generate(
start_port: Port,
length: float,
width: float,
snap_to_grid: bool = True,
dilation: float = 0.0,
snap_size: float = SEARCH_GRID_SNAP_UM,
) -> ComponentResult:
"""
Generate a straight waveguide segment.
Args:
start_port: Port to start from.
length: Requested length.
width: Waveguide width.
snap_to_grid: Whether to snap the end port to the search grid.
dilation: Optional dilation distance for pre-calculating collision geometry.
snap_size: Grid size for snapping.
Returns:
A ComponentResult containing the straight segment.
"""
rad = numpy.radians(start_port.orientation)
cos_val = numpy.cos(rad)
sin_val = numpy.sin(rad)
ex = start_port.x + length * cos_val
ey = start_port.y + length * sin_val
if snap_to_grid:
ex = snap_search_grid(ex, snap_size)
ey = snap_search_grid(ey, snap_size)
end_port = Port(ex, ey, start_port.orientation)
actual_length = numpy.sqrt((end_port.x - start_port.x)**2 + (end_port.y - start_port.y)**2)
# Create polygons using vectorized points
half_w = width / 2.0
pts_raw = numpy.array([
[0, half_w],
[actual_length, half_w],
[actual_length, -half_w],
[0, -half_w]
])
# Rotation matrix (standard 2D rotation)
rot_matrix = numpy.array([[cos_val, -sin_val], [sin_val, cos_val]])
# Transform points: P' = R * P + T
poly_points = (pts_raw @ rot_matrix.T) + [start_port.x, start_port.y]
geom = [Polygon(poly_points)]
dilated_geom = None
if dilation > 0:
# Direct calculation of dilated rectangle instead of expensive buffer()
half_w_dil = half_w + dilation
pts_dil = numpy.array([
[-dilation, half_w_dil],
[actual_length + dilation, half_w_dil],
[actual_length + dilation, -half_w_dil],
[-dilation, -half_w_dil]
])
poly_points_dil = (pts_dil @ rot_matrix.T) + [start_port.x, start_port.y]
dilated_geom = [Polygon(poly_points_dil)]
# For straight segments, geom IS the actual geometry
return ComponentResult(geometry=geom, end_port=end_port, length=actual_length, dilated_geometry=dilated_geom, actual_geometry=geom)
def _get_num_segments(radius: float, angle_deg: float, sagitta: float = 0.01) -> int:
"""
Calculate number of segments for an arc to maintain a maximum sagitta.
Args:
radius: Arc radius.
angle_deg: Total angle turned.
sagitta: Maximum allowed deviation.
Returns:
Minimum number of segments needed.
"""
if radius <= 0:
return 1
ratio = max(0.0, min(1.0, 1.0 - sagitta / radius))
theta_max = 2.0 * numpy.arccos(ratio)
if theta_max < 1e-9:
return 16
num = int(numpy.ceil(numpy.radians(abs(angle_deg)) / theta_max))
return max(8, num)
def _get_arc_polygons(
cx: float,
cy: float,
radius: float,
width: float,
t_start: float,
t_end: float,
sagitta: float = 0.01,
dilation: float = 0.0,
) -> list[Polygon]:
"""
Helper to generate arc-shaped polygons using vectorized NumPy operations.
Args:
cx, cy: Center coordinates.
radius: Arc radius.
width: Waveguide width.
t_start, t_end: Start and end angles (radians).
sagitta: Geometric fidelity.
dilation: Optional dilation to apply directly to the arc.
Returns:
List containing the arc polygon.
"""
num_segments = _get_num_segments(radius, float(numpy.degrees(abs(t_end - t_start))), sagitta)
angles = numpy.linspace(t_start, t_end, num_segments + 1)
cos_a = numpy.cos(angles)
sin_a = numpy.sin(angles)
inner_radius = radius - width / 2.0 - dilation
outer_radius = radius + width / 2.0 + dilation
inner_points = numpy.stack([cx + inner_radius * cos_a, cy + inner_radius * sin_a], axis=1)
outer_points = numpy.stack([cx + outer_radius * cos_a[::-1], cy + outer_radius * sin_a[::-1]], axis=1)
# Concatenate inner and outer points to form the polygon ring
poly_points = numpy.concatenate([inner_points, outer_points])
return [Polygon(poly_points)]
def _clip_bbox(
bbox: Polygon,
cx: float,
cy: float,
radius: float,
width: float,
clip_margin: float,
arc_poly: Polygon,
t_start: float | None = None,
t_end: float | None = None,
) -> Polygon:
"""
Clips corners of a bounding box for better collision modeling.
"""
r_out_cut = radius + width / 2.0 + clip_margin
r_in_cut = radius - width / 2.0 - clip_margin
# Angular range of the arc
if t_start is not None and t_end is not None:
ts, te = t_start, t_end
if ts > te:
ts, te = te, ts
# Sweep could cross 2pi boundary
sweep = (te - ts) % (2 * numpy.pi)
ts_norm = ts % (2 * numpy.pi)
else:
# Fallback: assume 90 deg based on centroid quadrant
ac = arc_poly.centroid
mid_angle = numpy.arctan2(ac.y - cy, ac.x - cx)
ts_norm = (mid_angle - numpy.pi/4) % (2 * numpy.pi)
sweep = numpy.pi/2
minx, miny, maxx, maxy = bbox.bounds
verts = [
numpy.array([minx, miny]),
numpy.array([maxx, miny]),
numpy.array([maxx, maxy]),
numpy.array([minx, maxy])
]
new_verts = []
for p in verts:
dx, dy = p[0] - cx, p[1] - cy
dist = numpy.sqrt(dx**2 + dy**2)
angle = numpy.arctan2(dy, dx)
# Check if corner angle is within the arc's angular sweep
angle_rel = (angle - ts_norm) % (2 * numpy.pi)
is_in_sweep = angle_rel <= sweep + 1e-6
d_line = -1.0
if is_in_sweep:
# We can clip if outside R_out or inside R_in
if dist > radius + width/2.0 - 1e-6:
d_line = r_out_cut * numpy.sqrt(2)
elif r_in_cut > 1e-3 and dist < radius - width/2.0 + 1e-6:
d_line = r_in_cut
else:
# Corner is outside angular sweep.
if dist > radius + width/2.0 - 1e-6:
d_line = r_out_cut * numpy.sqrt(2)
elif r_in_cut > 1e-3 and dist < radius - width/2.0 + 1e-6:
d_line = r_in_cut
if d_line > 0:
sx = 1.0 if dx > 0 else -1.0
sy = 1.0 if dy > 0 else -1.0
try:
# Intersection of line sx*(x-cx) + sy*(y-cy) = d_line with box edges
p_edge_x = numpy.array([p[0], cy + (d_line - sx * (p[0] - cx)) / sy])
p_edge_y = numpy.array([cx + (d_line - sy * (p[1] - cy)) / sx, p[1]])
# Check if intersection points are on the box boundary
if (minx - 1e-6 <= p_edge_y[0] <= maxx + 1e-6 and
miny - 1e-6 <= p_edge_x[1] <= maxy + 1e-6):
new_verts.append(p_edge_x)
new_verts.append(p_edge_y)
else:
new_verts.append(p)
except ZeroDivisionError:
new_verts.append(p)
else:
new_verts.append(p)
return Polygon(new_verts).convex_hull
def _apply_collision_model(
arc_poly: Polygon,
collision_type: Literal["arc", "bbox", "clipped_bbox"] | Polygon,
radius: float,
width: float,
cx: float = 0.0,
cy: float = 0.0,
clip_margin: float = 10.0,
t_start: float | None = None,
t_end: float | None = None,
) -> list[Polygon]:
"""
Applies the specified collision model to an arc geometry.
Args:
arc_poly: High-fidelity arc.
collision_type: Model type or custom polygon.
radius: Arc radius.
width: Waveguide width.
cx, cy: Arc center.
clip_margin: Safety margin for clipping.
t_start, t_end: Arc angles.
Returns:
List of polygons representing the collision model.
"""
if isinstance(collision_type, Polygon):
return [collision_type]
if collision_type == "arc":
return [arc_poly]
# Get bounding box
minx, miny, maxx, maxy = arc_poly.bounds
bbox = box(minx, miny, maxx, maxy)
if collision_type == "bbox":
return [bbox]
if collision_type == "clipped_bbox":
return [_clip_bbox(bbox, cx, cy, radius, width, clip_margin, arc_poly, t_start, t_end)]
return [arc_poly]
class Bend90:
"""
Move generator for 90-degree bends.
"""
@staticmethod
def generate(
start_port: Port,
radius: float,
width: float,
direction: str = "CW",
sagitta: float = 0.01,
collision_type: Literal["arc", "bbox", "clipped_bbox"] | Polygon = "arc",
clip_margin: float = 10.0,
dilation: float = 0.0,
snap_size: float = SEARCH_GRID_SNAP_UM,
) -> ComponentResult:
"""
Generate a 90-degree bend.
"""
turn_angle = -90 if direction == "CW" else 90
rad_start = numpy.radians(start_port.orientation)
c_angle = rad_start + (numpy.pi / 2 if direction == "CCW" else -numpy.pi / 2)
# Initial guess for center
cx_init = start_port.x + radius * numpy.cos(c_angle)
cy_init = start_port.y + radius * numpy.sin(c_angle)
t_start_init = c_angle + numpy.pi
t_end_init = t_start_init + (numpy.pi / 2 if direction == "CCW" else -numpy.pi / 2)
# Snap the target point
ex = snap_search_grid(cx_init + radius * numpy.cos(t_end_init), snap_size)
ey = snap_search_grid(cy_init + radius * numpy.sin(t_end_init), snap_size)
end_port = Port(ex, ey, float((start_port.orientation + turn_angle) % 360))
# Adjust geometry to perfectly hit snapped port
dx = ex - start_port.x
dy = ey - start_port.y
dist = numpy.sqrt(dx**2 + dy**2)
# New radius for the right triangle connecting start to end with 90 deg
actual_radius = dist / numpy.sqrt(2)
# Vector from start to end
mid_x, mid_y = (start_port.x + ex)/2, (start_port.y + ey)/2
# Normal vector (orthogonal to start->end)
# Flip direction based on CW/CCW
dir_sign = 1 if direction == "CCW" else -1
cx = mid_x - dir_sign * (ey - start_port.y) / 2
cy = mid_y + dir_sign * (ex - start_port.x) / 2
# Update angles based on new center
t_start = numpy.arctan2(start_port.y - cy, start_port.x - cx)
t_end = numpy.arctan2(ey - cy, ex - cx)
# Maintain directionality and angular span near pi/2
if direction == "CCW":
while t_end < t_start: t_end += 2 * numpy.pi
else:
while t_end > t_start: t_end -= 2 * numpy.pi
arc_polys = _get_arc_polygons(cx, cy, actual_radius, width, t_start, t_end, sagitta)
collision_polys = _apply_collision_model(
arc_polys[0], collision_type, actual_radius, width, cx, cy, clip_margin, t_start, t_end
)
proxy_geom = None
if collision_type == "arc":
# Auto-generate a clipped_bbox proxy for tiered collision checks
proxy_geom = _apply_collision_model(
arc_polys[0], "clipped_bbox", actual_radius, width, cx, cy, clip_margin, t_start, t_end
)
dilated_geom = None
if dilation > 0:
if collision_type == "arc":
dilated_geom = _get_arc_polygons(cx, cy, actual_radius, width, t_start, t_end, sagitta, dilation=dilation)
else:
dilated_geom = [p.buffer(dilation) for p in collision_polys]
return ComponentResult(
geometry=collision_polys,
end_port=end_port,
length=actual_radius * numpy.abs(t_end - t_start),
dilated_geometry=dilated_geom,
proxy_geometry=proxy_geom,
actual_geometry=arc_polys
)
class SBend:
"""
Move generator for parametric S-bends.
"""
@staticmethod
def generate(
start_port: Port,
offset: float,
radius: float,
width: float,
sagitta: float = 0.01,
collision_type: Literal["arc", "bbox", "clipped_bbox"] | Polygon = "arc",
clip_margin: float = 10.0,
dilation: float = 0.0,
snap_size: float = SEARCH_GRID_SNAP_UM,
) -> ComponentResult:
"""
Generate a parametric S-bend (two tangent arcs).
"""
if abs(offset) >= 2 * radius:
raise ValueError(f"SBend offset {offset} must be less than 2*radius {2 * radius}")
theta_init = numpy.arccos(1 - abs(offset) / (2 * radius))
dx_init = 2 * radius * numpy.sin(theta_init)
rad_start = numpy.radians(start_port.orientation)
# Snap the target point
ex = snap_search_grid(start_port.x + dx_init * numpy.cos(rad_start) - offset * numpy.sin(rad_start), snap_size)
ey = snap_search_grid(start_port.y + dx_init * numpy.sin(rad_start) + offset * numpy.cos(rad_start), snap_size)
end_port = Port(ex, ey, start_port.orientation)
# Solve for theta and radius that hit (ex, ey) exactly
local_dx = (ex - start_port.x) * numpy.cos(rad_start) + (ey - start_port.y) * numpy.sin(rad_start)
local_dy = -(ex - start_port.x) * numpy.sin(rad_start) + (ey - start_port.y) * numpy.cos(rad_start)
# tan(theta / 2) = local_dy / local_dx
theta = 2 * numpy.arctan2(abs(local_dy), local_dx)
# Avoid division by zero if theta is 0 (though unlikely due to offset check)
actual_radius = abs(local_dy) / (2 * (1 - numpy.cos(theta))) if theta > 1e-9 else radius
direction = 1 if local_dy > 0 else -1
c1_angle = rad_start + direction * numpy.pi / 2
cx1 = start_port.x + actual_radius * numpy.cos(c1_angle)
cy1 = start_port.y + actual_radius * numpy.sin(c1_angle)
ts1, te1 = c1_angle + numpy.pi, c1_angle + numpy.pi + direction * theta
c2_angle = rad_start - direction * numpy.pi / 2
cx2 = ex + actual_radius * numpy.cos(c2_angle)
cy2 = ey + actual_radius * numpy.sin(c2_angle)
te2 = c2_angle + numpy.pi
ts2 = te2 + direction * theta
arc1 = _get_arc_polygons(cx1, cy1, actual_radius, width, ts1, te1, sagitta)[0]
arc2 = _get_arc_polygons(cx2, cy2, actual_radius, width, ts2, te2, sagitta)[0]
arc_polys = [arc1, arc2]
# Use the provided collision model for primary geometry
col1 = _apply_collision_model(arc1, collision_type, actual_radius, width, cx1, cy1, clip_margin, ts1, te1)[0]
col2 = _apply_collision_model(arc2, collision_type, actual_radius, width, cx2, cy2, clip_margin, ts2, te2)[0]
collision_polys = [col1, col2]
proxy_geom = None
if collision_type == "arc":
# Auto-generate proxies
p1 = _apply_collision_model(arc1, "clipped_bbox", actual_radius, width, cx1, cy1, clip_margin, ts1, te1)[0]
p2 = _apply_collision_model(arc2, "clipped_bbox", actual_radius, width, cx2, cy2, clip_margin, ts2, te2)[0]
proxy_geom = [p1, p2]
dilated_geom = None
if dilation > 0:
if collision_type == "arc":
d1 = _get_arc_polygons(cx1, cy1, actual_radius, width, ts1, te1, sagitta, dilation=dilation)[0]
d2 = _get_arc_polygons(cx2, cy2, actual_radius, width, ts2, te2, sagitta, dilation=dilation)[0]
dilated_geom = [d1, d2]
else:
dilated_geom = [p.buffer(dilation) for p in collision_polys]
return ComponentResult(
geometry=collision_polys,
end_port=end_port,
length=2 * actual_radius * theta,
dilated_geometry=dilated_geom,
proxy_geometry=proxy_geom,
actual_geometry=arc_polys
)