[Ellipse / Arc] improve bounds calculation

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 = []
+        def point(rx: float, ry: float, tt: float) -> NDArray[numpy.float64]:
-        maxs = []
+            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
+                pts.append(numpy.zeros(2))
-                mins.append([0, 0])
-                maxs.append([0, 0])
                 continue
+            pts.extend(points_in_interval(rx, ry, aa[0], aa[1]))
 
-            a0, a1 = aa
+        all_pts = numpy.asarray(pts) + self.offset
-            a0_offset = a0 - (a0 % (2 * pi))
+        return numpy.vstack((numpy.min(all_pts, axis=0),
-
+                             numpy.max(all_pts, axis=0)))
-            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

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)
+        cos_r = numpy.cos(self.rotation)
-        return numpy.vstack((self.offset - rot_radii[0],
+        sin_r = numpy.sin(self.rotation)
-                             self.offset + rot_radii[1]))
+        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

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)