diff --git a/masque/shapes/arc.py b/masque/shapes/arc.py index 4befb6d..a068700 100644 --- a/masque/shapes/arc.py +++ b/masque/shapes/arc.py @@ -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/ellipse.py b/masque/shapes/ellipse.py index 985d252..8a5c287 100644 --- a/masque/shapes/ellipse.py +++ b/masque/shapes/ellipse.py @@ -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_shape_advanced.py b/masque/test/test_shape_advanced.py index 046df1a..ab42b4b 100644 --- a/masque/test/test_shape_advanced.py +++ b/masque/test/test_shape_advanced.py @@ -97,6 +97,18 @@ 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_path_edge_cases() -> None: # Zero-length segments p = MPath(vertices=[[0, 0], [0, 0], [10, 0]], width=2)