diff --git a/masque/file/oasis.py b/masque/file/oasis.py index b5d0cd8..95a2039 100644 --- a/masque/file/oasis.py +++ b/masque/file/oasis.py @@ -562,9 +562,7 @@ 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), None) # reverse lookup - if path_type is None: - raise PatternError(f'OASIS writer does not support path cap {shape.cap}') + path_type = next(k for k, v in path_cap_map.items() if v == shape.cap) # reverse lookup 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 ab2e47b..4befb6d 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 = max(2, int(numpy.ceil(2 * pi * max(self.radii + dr) / max_arclen))) + n_pts = numpy.ceil(2 * pi * max(self.radii + dr) / max_arclen).astype(int) arc_lengths, thetas = get_arclens(n_pts, *a_ranges[0 if inner else 1], dr=dr) keep = [0] @@ -313,48 +313,77 @@ 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) - 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 = [] + mins = [] + maxs = [] 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: - pts.append(numpy.zeros(2)) - continue - pts.extend(points_in_interval(rx, ry, aa[0], aa[1])) - all_pts = numpy.asarray(pts) + self.offset - return numpy.vstack((numpy.min(all_pts, axis=0), - numpy.max(all_pts, axis=0))) + if rx == 0 or ry == 0: + # Single point, at origin + mins.append([0, 0]) + maxs.append([0, 0]) + continue + + 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)) def rotate(self, theta: float) -> 'Arc': self.rotation += theta diff --git a/masque/shapes/circle.py b/masque/shapes/circle.py index a6f7af1..8dad165 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 = max(3, int(round(max(n)))) + num_vertices = 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 55ce9fd..985d252 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 = max(3, int(round(max(n)))) + num_vertices = 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,13 +180,9 @@ class Ellipse(PositionableImpl, Shape): return [poly] def get_bounds_single(self) -> NDArray[numpy.float64]: - 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)) + rot_radii = numpy.dot(rotation_matrix_2d(self.rotation), self.radii) + return numpy.vstack((self.offset - rot_radii[0], + self.offset + rot_radii[1])) def rotate(self, theta: float) -> Self: self.rotation += theta diff --git a/masque/test/test_oasis.py b/masque/test/test_oasis.py index f549db7..ad9bbf9 100644 --- a/masque/test/test_oasis.py +++ b/masque/test/test_oasis.py @@ -1,12 +1,9 @@ -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: @@ -45,16 +42,3 @@ 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 fc76fcb..046df1a 100644 --- a/masque/test/test_shape_advanced.py +++ b/masque/test/test_shape_advanced.py @@ -97,31 +97,6 @@ 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)