diff --git a/masque/file/oasis.py b/masque/file/oasis.py index 95a2039..b5d0cd8 100644 --- a/masque/file/oasis.py +++ b/masque/file/oasis.py @@ -562,7 +562,9 @@ def _shapes_to_elements( xy = rint_cast(shape.offset + shape.vertices[0] + rep_offset) deltas = rint_cast(numpy.diff(shape.vertices, axis=0)) half_width = rint_cast(shape.width / 2) - path_type = next(k for k, v in path_cap_map.items() if v == shape.cap) # reverse lookup + path_type = next((k for k, v in path_cap_map.items() if v == shape.cap), None) # reverse lookup + if path_type is None: + raise PatternError(f'OASIS writer does not support path cap {shape.cap}') extension_start = (path_type, shape.cap_extensions[0] if shape.cap_extensions is not None else None) extension_end = (path_type, shape.cap_extensions[1] if shape.cap_extensions is not None else None) path = fatrec.Path( diff --git a/masque/shapes/arc.py b/masque/shapes/arc.py index 4befb6d..ab2e47b 100644 --- a/masque/shapes/arc.py +++ b/masque/shapes/arc.py @@ -268,7 +268,7 @@ class Arc(PositionableImpl, Shape): """ Figure out the parameter values at which we should place vertices to meet the arclength constraint""" dr = -wh if inner else wh - n_pts = numpy.ceil(2 * pi * max(self.radii + dr) / max_arclen).astype(int) + n_pts = max(2, int(numpy.ceil(2 * pi * max(self.radii + dr) / max_arclen))) arc_lengths, thetas = get_arclens(n_pts, *a_ranges[0 if inner else 1], dr=dr) keep = [0] @@ -313,77 +313,48 @@ class Arc(PositionableImpl, Shape): return [poly] def get_bounds_single(self) -> NDArray[numpy.float64]: - """ - Equation for rotated ellipse is - `x = x0 + a * cos(t) * cos(rot) - b * sin(t) * sin(phi)` - `y = y0 + a * cos(t) * sin(rot) + b * sin(t) * cos(rot)` - where `t` is our parameter. - - Differentiating and solving for 0 slope wrt. `t`, we find - `tan(t) = -+ b/a cot(phi)` - where -+ is for x, y cases, so that's where the extrema are. - - If the extrema are innaccessible due to arc constraints, check the arc endpoints instead. - """ a_ranges = cast('_array2x2_t', self._angles_to_parameters()) + sin_r = numpy.sin(self.rotation) + cos_r = numpy.cos(self.rotation) - mins = [] - maxs = [] + def point(rx: float, ry: float, tt: float) -> NDArray[numpy.float64]: + return numpy.array(( + rx * numpy.cos(tt) * cos_r - ry * numpy.sin(tt) * sin_r, + rx * numpy.cos(tt) * sin_r + ry * numpy.sin(tt) * cos_r, + )) + + def points_in_interval(rx: float, ry: float, a0: float, a1: float) -> list[NDArray[numpy.float64]]: + candidates = [a0, a1] + if rx != 0 and ry != 0: + tx = numpy.arctan2(-ry * sin_r, rx * cos_r) + ty = numpy.arctan2(ry * cos_r, rx * sin_r) + candidates.extend((tx, tx + pi, ty, ty + pi)) + + lo = min(a0, a1) + hi = max(a0, a1) + pts = [] + for base in candidates: + k_min = int(numpy.floor((lo - base) / (2 * pi))) - 1 + k_max = int(numpy.ceil((hi - base) / (2 * pi))) + 1 + for kk in range(k_min, k_max + 1): + tt = base + kk * 2 * pi + if lo <= tt <= hi: + pts.append(point(rx, ry, tt)) + return pts + + pts = [] for aa, sgn in zip(a_ranges, (-1, +1), strict=True): wh = sgn * self.width / 2 rx = self.radius_x + wh ry = self.radius_y + wh - if rx == 0 or ry == 0: - # Single point, at origin - mins.append([0, 0]) - maxs.append([0, 0]) + pts.append(numpy.zeros(2)) continue + pts.extend(points_in_interval(rx, ry, aa[0], aa[1])) - a0, a1 = aa - a0_offset = a0 - (a0 % (2 * pi)) - - sin_r = numpy.sin(self.rotation) - cos_r = numpy.cos(self.rotation) - sin_a = numpy.sin(aa) - cos_a = numpy.cos(aa) - - # Cutoff angles - xpt = (-self.rotation) % (2 * pi) + a0_offset - ypt = (pi / 2 - self.rotation) % (2 * pi) + a0_offset - xnt = (xpt - pi) % (2 * pi) + a0_offset - ynt = (ypt - pi) % (2 * pi) + a0_offset - - # Points along coordinate axes - rx2_inv = 1 / (rx * rx) - ry2_inv = 1 / (ry * ry) - xr = numpy.abs(cos_r * cos_r * rx2_inv + sin_r * sin_r * ry2_inv) ** -0.5 - yr = numpy.abs(-sin_r * -sin_r * rx2_inv + cos_r * cos_r * ry2_inv) ** -0.5 - - # Arc endpoints - xn, xp = sorted(rx * cos_r * cos_a - ry * sin_r * sin_a) - yn, yp = sorted(rx * sin_r * cos_a + ry * cos_r * sin_a) - - # If our arc subtends a coordinate axis, use the extremum along that axis - if abs(a1 - a0) >= 2 * pi: - xn, xp, yn, yp = -xr, xr, -yr, yr - else: - if a0 <= xpt <= a1 or a0 <= xpt + 2 * pi <= a1: - xp = xr - - if a0 <= xnt <= a1 or a0 <= xnt + 2 * pi <= a1: - xn = -xr - - if a0 <= ypt <= a1 or a0 <= ypt + 2 * pi <= a1: - yp = yr - - if a0 <= ynt <= a1 or a0 <= ynt + 2 * pi <= a1: - yn = -yr - - mins.append([xn, yn]) - maxs.append([xp, yp]) - return numpy.vstack((numpy.min(mins, axis=0) + self.offset, - numpy.max(maxs, axis=0) + self.offset)) + all_pts = numpy.asarray(pts) + self.offset + return numpy.vstack((numpy.min(all_pts, axis=0), + numpy.max(all_pts, axis=0))) def rotate(self, theta: float) -> 'Arc': self.rotation += theta diff --git a/masque/shapes/circle.py b/masque/shapes/circle.py index 8dad165..a6f7af1 100644 --- a/masque/shapes/circle.py +++ b/masque/shapes/circle.py @@ -108,7 +108,7 @@ class Circle(PositionableImpl, Shape): n += [num_vertices] if max_arclen is not None: n += [2 * pi * self.radius / max_arclen] - num_vertices = int(round(max(n))) + num_vertices = max(3, int(round(max(n)))) thetas = numpy.linspace(2 * pi, 0, num_vertices, endpoint=False) xs = numpy.cos(thetas) * self.radius ys = numpy.sin(thetas) * self.radius diff --git a/masque/shapes/ellipse.py b/masque/shapes/ellipse.py index 985d252..55ce9fd 100644 --- a/masque/shapes/ellipse.py +++ b/masque/shapes/ellipse.py @@ -168,7 +168,7 @@ class Ellipse(PositionableImpl, Shape): n += [num_vertices] if max_arclen is not None: n += [perimeter / max_arclen] - num_vertices = int(round(max(n))) + num_vertices = max(3, int(round(max(n)))) thetas = numpy.linspace(2 * pi, 0, num_vertices, endpoint=False) sin_th, cos_th = (numpy.sin(thetas), numpy.cos(thetas)) @@ -180,9 +180,13 @@ class Ellipse(PositionableImpl, Shape): return [poly] def get_bounds_single(self) -> NDArray[numpy.float64]: - rot_radii = numpy.dot(rotation_matrix_2d(self.rotation), self.radii) - return numpy.vstack((self.offset - rot_radii[0], - self.offset + rot_radii[1])) + cos_r = numpy.cos(self.rotation) + sin_r = numpy.sin(self.rotation) + x_extent = numpy.sqrt((self.radius_x * cos_r) ** 2 + (self.radius_y * sin_r) ** 2) + y_extent = numpy.sqrt((self.radius_x * sin_r) ** 2 + (self.radius_y * cos_r) ** 2) + extents = numpy.array((x_extent, y_extent)) + return numpy.vstack((self.offset - extents, + self.offset + extents)) def rotate(self, theta: float) -> Self: self.rotation += theta diff --git a/masque/test/test_oasis.py b/masque/test/test_oasis.py index ad9bbf9..f549db7 100644 --- a/masque/test/test_oasis.py +++ b/masque/test/test_oasis.py @@ -1,9 +1,12 @@ +import io from pathlib import Path import pytest from numpy.testing import assert_equal +from ..error import PatternError from ..pattern import Pattern from ..library import Library +from ..shapes import Path as MPath def test_oasis_roundtrip(tmp_path: Path) -> None: @@ -42,3 +45,16 @@ def test_oasis_properties_to_annotations_merges_repeated_keys() -> None: ) assert annotations == {"k": [1, 2, 3]} + + +def test_oasis_write_rejects_circle_path_caps() -> None: + pytest.importorskip("fatamorgana") + from ..file import oasis + + lib = Library() + pat = Pattern() + pat.path((1, 0), vertices=[[0, 0], [10, 0]], width=2, cap=MPath.Cap.Circle) + lib["cell1"] = pat + + with pytest.raises(PatternError, match="does not support path cap"): + oasis.write(lib, io.BytesIO(), units_per_micron=1000) diff --git a/masque/test/test_shape_advanced.py b/masque/test/test_shape_advanced.py index 046df1a..fc76fcb 100644 --- a/masque/test/test_shape_advanced.py +++ b/masque/test/test_shape_advanced.py @@ -97,6 +97,31 @@ def test_arc_edge_cases() -> None: assert_allclose(bounds, [[-11, -11], [11, 11]], atol=1e-10) +def test_rotated_ellipse_bounds_match_polygonized_geometry() -> None: + ellipse = Ellipse(radii=(10, 20), rotation=pi / 4, offset=(100, 200)) + bounds = ellipse.get_bounds_single() + poly_bounds = ellipse.to_polygons(num_vertices=8192)[0].get_bounds_single() + assert_allclose(bounds, poly_bounds, atol=1e-3) + + +def test_rotated_arc_bounds_match_polygonized_geometry() -> None: + arc = Arc(radii=(10, 20), angles=(0, pi), width=2, rotation=pi / 4, offset=(100, 200)) + bounds = arc.get_bounds_single() + poly_bounds = arc.to_polygons(num_vertices=8192)[0].get_bounds_single() + assert_allclose(bounds, poly_bounds, atol=1e-3) + + +def test_curve_polygonizers_clamp_large_max_arclen() -> None: + for shape in ( + Circle(radius=10), + Ellipse(radii=(10, 20)), + Arc(radii=(10, 20), angles=(0, 1), width=2), + ): + polys = shape.to_polygons(num_vertices=None, max_arclen=1e9) + assert len(polys) == 1 + assert len(polys[0].vertices) >= 3 + + def test_path_edge_cases() -> None: # Zero-length segments p = MPath(vertices=[[0, 0], [0, 0], [10, 0]], width=2)