Fix VMobject dimensions and critical points using Bézier extrema

manim/mobject/types/vectorized_mobject.py

@@ -1796,6 +1796,202 @@ def get_points_defining_boundary(self) -> Point3D_Array:
             tuple(it.chain(*(sm.get_anchors() for sm in self.get_family())))
         )
 
+    def _get_bezier_family_bounding_box(self) -> Point3D_Array | None:
+        """Return the exact bounding box of the curves of the family.
+
+        Aggregates the boxes of every :class:`~Mobject.get_family` member
+        that carries points.  Members storing curves other than cubics fall
+        back to their raw point bounds.
+        """
+        members = [
+            (m.points, m.n_points_per_cubic_curve)
+            for m in self.get_family()
+            if len(m.points) > 0
+        ]
+        if not members:
+            return None
+
+        if all(n == 4 and len(p) % 4 == 0 for p, n in members):
+            pts = np.concatenate([p for p, _ in members])
+            lower = np.zeros(3)
+            upper = np.zeros(3)
+            for dim in range(3):
+                vals = self._get_curve_extrema(dim, pts, 4)
+                lower[dim] = np.min(vals[:, 0])
+                upper[dim] = np.max(vals[:, 1])
+            return np.array([lower, upper])
+
+        bbox: Point3D_Array | None = None
+        for p, n in members:
+            if n == 4 and len(p) % n == 0:
+                bb = np.zeros((2, 3))
+                for dim in range(3):
+                    vals = self._get_curve_extrema(dim, p, n)
+                    bb[0, dim] = np.min(vals[:, 0])
+                    bb[1, dim] = np.max(vals[:, 1])
+            else:
+                bb = np.array([p.min(axis=0), p.max(axis=0)])
+            if bbox is None:
+                bbox = bb
+            else:
+                bbox[0] = np.minimum(bbox[0], bb[0])
+                bbox[1] = np.maximum(bbox[1], bb[1])
+        return bbox
+
+    def get_bezier_bounding_box(self) -> Point3D_Array | None:
+        """Return the exact axis-aligned bounding box of this
+        :class:`VMobject`'s Bézier curves.
+
+        Unlike bounds computed from ``self.points``, this includes only
+        points on the curves, accounting for interior extrema.
+
+        Returns ``None`` if the VMobject has no points.
+        """
+        pts = self.points
+        if len(pts) == 0:
+            return None
+        nppcc = self.n_points_per_cubic_curve
+        if nppcc != 4 or len(pts) % nppcc != 0:
+            return np.array([pts.min(axis=0), pts.max(axis=0)])
+
+        lower = np.zeros(3)
+        upper = np.zeros(3)
+        for dim in range(3):
+            vals = self._get_curve_extrema(dim, pts, nppcc)
+            lower[dim] = np.min(vals[:, 0])
+            upper[dim] = np.max(vals[:, 1])
+        return np.array([lower, upper])
+
+    def _get_curve_extrema(
+        self, dim: int, pts: Point3D_Array, nppcc: int
+    ) -> npt.NDArray[np.float64]:
+        """Return, for every curve in ``pts``, the exact minimum and maximum
+        value of its ``dim``-th coordinate (anchors and interior extrema).
+
+        Returns an array of shape ``(n_curves, 2)``.  Callers must ensure
+        ``nppcc == 4`` and ``len(pts) % nppcc == 0`` (standard VMobjects
+        store cubic curves; SVG quadratic segments are degree-elevated
+        before storage).
+        """
+        n_curves = len(pts) // nppcc
+        # Cubic Bézier: derivative roots of P'(t) = A t^2 + B t + C.
+        p0, p1, p2, p3 = (pts[i::nppcc, dim] for i in range(nppcc))
+        u = p1 - p0
+        v = p2 - p1
+        w = p3 - p2
+        aa = u - 2 * v + w
+        bb = 2 * (v - u)
+        cc = u
+        # Tolerances are relative to the coefficient scale so that scaling
+        # the whole curve does not change which roots are classified as
+        # interior extrema.
+        m = np.maximum(np.maximum(np.abs(aa), np.abs(bb)), np.abs(cc))
+        rel = 1e-12 * m
+        disc = bb * bb - 4 * aa * cc
+        disc_pos = disc > rel * rel
+        aa_zero = np.abs(aa) < rel
+        t = np.full((n_curves, 2), np.nan)
+        with np.errstate(divide="ignore", invalid="ignore"):
+            sq = np.sqrt(np.maximum(disc, 0))
+            # A == 0: degenerate quadratic, single root -C/B (B != 0).
+            t[:, 0] = np.where(
+                aa_zero & ~np.isclose(bb, 0, atol=rel),
+                -cc / bb,
+                (-bb + sq) / (2 * aa),
+            )
+            t[:, 1] = np.where(aa_zero, t[:, 0], (-bb - sq) / (2 * aa))
+        start = p0
+        end = p3
+        mins = np.minimum(start, end).astype(np.float64)
+        maxs = np.maximum(start, end).astype(np.float64)
+        for j in range(2):
+            tv = t[:, j]
+            quadratic_root = aa_zero & ~np.isclose(bb, 0, atol=rel)
+            valid = (
+                (disc_pos & ~aa_zero | quadratic_root) & (tv > 1e-12) & (tv < 1 - 1e-12)
+            )
+            if np.any(valid):
+                tvv = tv[valid]
+                omt = 1 - tvv
+                ev = (
+                    omt**3 * p0[valid]
+                    + 3 * omt**2 * tvv * p1[valid]
+                    + 3 * omt * tvv**2 * p2[valid]
+                    + tvv**3 * p3[valid]
+                )
+                mins[valid] = np.minimum(mins[valid], ev)
+                maxs[valid] = np.maximum(maxs[valid], ev)
+        return np.column_stack([mins, maxs])
+
+    def get_extremum_along_dim(
+        self,
+        points: Point3DLike_Array | None = None,
+        dim: int = 0,
+        key: int = 0,
+    ) -> float:
+        """Extremum of the exact curve bounds of the family along ``dim``.
+
+        When ``points`` are passed explicitly, the request refers to those
+        points and is delegated to
+        :meth:`~Mobject.get_extremum_along_dim`.  Otherwise the exact
+        bounding box of the full family is used, so planning methods such
+        as ``get_coord``, ``set_coord`` and ``align_to`` are consistent
+        with :meth:`get_critical_point` and the width/height/depth setters.
+        """
+        if points is not None:
+            return super().get_extremum_along_dim(points, dim, key)
+        bbox = self._get_bezier_family_bounding_box()
+        if bbox is None:
+            return 0.0
+        if key < 0:
+            return bbox[0, dim]
+        if key > 0:
+            return bbox[1, dim]
+        return 0.5 * (bbox[0, dim] + bbox[1, dim])
+
+    def get_critical_point(self, direction: Vector3DLike) -> Point3D:
+        """Return one of the 9 'critical points' of the bounding box, along the
+        given direction.
+
+        Unlike :meth:`~.Mobject.get_critical_point`, the bounding box is the
+        exact box of the rendered Bézier curves (see
+        :meth:`get_bezier_bounding_box`), not the box of their control points.
+        ``get_left()``, ``get_right()``, ``get_top()``, ``get_bottom()``,
+        ``get_center()`` etc. therefore always agree with ``width``, ``height``
+        and ``depth``.
+
+        See :meth:`~.Mobject.get_critical_point` for details and examples.
+        """
+        result = np.zeros(self.dim)
+        bbox = self._get_bezier_family_bounding_box()
+        if bbox is None:
+            return result
+        for dim in range(self.dim):
+            key = direction[dim]
+            if key > 0:
+                result[dim] = bbox[1][dim]
+            elif key < 0:
+                result[dim] = bbox[0][dim]
+            else:
+                result[dim] = 0.5 * (bbox[0][dim] + bbox[1][dim])
+        return result
+
+    def length_over_dim(self, dim: int) -> float:
+        """Find the length of this :class:`VMobject` in a certain direction.
+
+        Like :meth:`~.Mobject.length_over_dim`, this covers every point in
+        this :class:`VMobject` and its submobjects, and the length is
+        computed from the actual Bézier curves (including their interior
+        extrema) rather than from the raw control points.  Dimensions such
+        as :attr:`~Mobject.width` and :attr:`~Mobject.height` therefore
+        correspond to the physical extent of the rendered shape, and the
+        critical points (:meth:`~Mobject.get_left` etc.) agree with them.
+        """
+        bbox = self._get_bezier_family_bounding_box()
+        if bbox is None:
+            return 0.0
+        return bbox[1][dim] - bbox[0][dim]
+
     def get_arc_length(self, sample_points_per_curve: int | None = None) -> float:
         """Return the approximated length of the whole curve.
 

tests/module/mobject/types/vectorized_mobject/test_vectorized_mobject.py

@@ -4,6 +4,7 @@
 import pytest
 
 from manim import (
+    Arc,
     Circle,
     CurvesAsSubmobjects,
     Line,
@@ -15,7 +16,7 @@
     VGroup,
     VMobject,
 )
-from manim.constants import PI
+from manim.constants import DEGREES, LEFT, PI, RIGHT, TAU
 
 
 def test_vmobject_add():
@@ -774,3 +775,150 @@ def test_pointwise_become_partial_where_vmobject_is_self():
         ]
     )
     np.testing.assert_allclose(sq.points, expected_points)
+
+
+def test_width_height_account_for_curve_interior_extrema():
+    """Handles (control points) must not inflate width/height (#3619)."""
+    c = Circle(radius=3).rotate(30 * DEGREES)
+    assert c.width == pytest.approx(6.0, abs=1e-3)
+    assert c.height == pytest.approx(6.0, abs=1e-3)
+
+
+def test_width_height_include_interior_extrema_of_wild_cubic():
+    # Handles at x=5 and x=-5 with anchors x=0 and x=1: the control-point
+    # bounding box would give width 10, but the curve itself peaks at
+    # x ~= 1.45299 and dips to x ~= -0.98473, so the exact width is ~2.4377.
+    vmob = VMobject().set_points(
+        np.array([[0.0, 0.0, 0.0], [5.0, 3.0, 0.0], [-5.0, 3.0, 0.0], [1.0, 0.0, 0.0]])
+    )
+    assert vmob.width == pytest.approx(2.4377191199218955, abs=1e-4)
+    assert vmob.height == pytest.approx(2.25, abs=1e-4)
+
+
+def test_arc_critical_points_agree_with_width_height():
+    # Choose an arc whose x-extremum lies between subdivision anchors.
+    # Control-point, anchor-only, and exact curve bounds are all different.
+    arc = Arc(radius=2, start_angle=PI / 5, angle=TAU * 0.7)
+    assert arc.width == pytest.approx(3.618034, abs=1e-4)
+    assert arc.height == pytest.approx(4.0, abs=1e-3)
+    assert arc.get_right()[0] - arc.get_left()[0] == pytest.approx(arc.width, abs=1e-6)
+    assert arc.get_top()[1] - arc.get_bottom()[1] == pytest.approx(arc.height, abs=1e-6)
+
+
+@pytest.mark.parametrize("dim", ["width", "height", "depth"])
+def test_width_height_depth_setters_round_trip(dim):
+    # Setting a dimension must result in exactly that dimension, regardless of
+    # which one is set.  Explicit VMobject so the test does not depend on the
+    # point layout chosen by a specific shape.
+    mob = VMobject().set_points(
+        np.array(
+            [
+                [0.0, 0.0, 0.0],
+                [3.0, 2.0, 1.0],
+                [-2.0, 1.0, -1.0],
+                [1.0, 0.0, 0.0],
+            ]
+        )
+    )
+    setattr(mob, dim, 5.0)
+    assert getattr(mob, dim) == pytest.approx(5.0)
+
+
+def test_width_height_of_single_point_vmobject():
+    vmob = VMobject().set_points(np.array([[3.0, 4.0, 0.0]]))
+    assert vmob.width == pytest.approx(0.0)
+    assert vmob.height == pytest.approx(0.0)
+
+
+def test_critical_points_agree_with_width_height():
+    # Critical points and dimensions must use the same exact curve bounds.
+    vmob = VMobject().set_points(
+        np.array([[0.0, 0.0, 0.0], [5.0, 3.0, 0.0], [-5.0, 3.0, 0.0], [1.0, 0.0, 0.0]])
+    )
+    assert vmob.width == pytest.approx(2.4377191199218955, abs=1e-4)
+    assert vmob.get_left()[0] == pytest.approx(-0.98472845, abs=1e-6)
+    assert vmob.get_right()[0] == pytest.approx(1.45299067, abs=1e-6)
+    assert vmob.get_bottom()[1] == pytest.approx(0.0, abs=1e-6)
+    assert vmob.get_top()[1] == pytest.approx(2.25, abs=1e-6)
+    assert vmob.get_right()[0] - vmob.get_left()[0] == pytest.approx(
+        vmob.width, abs=1e-6
+    )
+    assert vmob.get_top()[1] - vmob.get_bottom()[1] == pytest.approx(
+        vmob.height, abs=1e-6
+    )
+    assert vmob.get_center()[0] == pytest.approx(
+        (vmob.get_left()[0] + vmob.get_right()[0]) / 2, abs=1e-6
+    )
+
+
+def test_depth_and_critical_points_cover_curve_extrema():
+    # Control points along z span [-5, 5], but the rendered curve only
+    # reaches ~[-0.985, 1.453].  depth and the OUT/IN critical points must
+    # agree with each other and stay within the control hull.
+    vmob = VMobject().set_points(
+        np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 5.0], [0.0, 0.0, -5.0], [0.0, 0.0, 1.0]])
+    )
+    assert vmob.depth == pytest.approx(2.437719, abs=1e-4)
+    out = vmob.get_critical_point(np.array([0.0, 0.0, 1.0]))
+    inn = vmob.get_critical_point(np.array([0.0, 0.0, -1.0]))
+    assert out[2] == pytest.approx(1.45299067, abs=1e-6)
+    assert inn[2] == pytest.approx(-0.98472845, abs=1e-6)
+    assert out[2] - inn[2] == pytest.approx(vmob.depth, abs=1e-6)
+
+
+def _wild_cubic() -> VMobject:
+    return VMobject().set_points(
+        np.array([[0.0, 0.0, 0.0], [5.0, 3.0, 0.0], [-5.0, 3.0, 0.0], [1.0, 0.0, 0.0]])
+    )
+
+
+def test_coord_planning_uses_exact_extrema():
+    # get_x/get_coord/set_coord/align_to must be consistent with the exact
+    # critical points, otherwise set_x and align_to move wrong amounts.
+    vmob = _wild_cubic()
+    assert vmob.get_x() == pytest.approx(vmob.get_center()[0])
+    assert vmob.get_x(RIGHT) == pytest.approx(vmob.get_right()[0])
+    assert vmob.get_x(LEFT) == pytest.approx(vmob.get_left()[0])
+    assert vmob.get_extremum_along_dim(dim=0, key=1) == pytest.approx(
+        vmob.get_right()[0]
+    )
+    assert vmob.get_extremum_along_dim(dim=0, key=-1) == pytest.approx(
+        vmob.get_left()[0]
+    )
+
+    vmob.set_x(0)
+    assert vmob.get_center()[0] == pytest.approx(0, abs=1e-9)
+
+    vmob = _wild_cubic()
+    vmob.set_x(0, RIGHT)
+    assert vmob.get_right()[0] == pytest.approx(0, abs=1e-9)
+
+    rect = Square(side_length=2).shift(3 * RIGHT)
+    vmob = _wild_cubic()
+    vmob.align_to(rect, RIGHT)
+    assert vmob.get_right()[0] == pytest.approx(rect.get_right()[0], abs=1e-9)
+
+
+def test_extrema_scale_invariance():
+    # Scaling a curve must not change which roots are treated as interior
+    # extrema, so width scales exactly with the overall scale factor.
+    true_width = 2.4377191199218955
+    pts = np.array(
+        [[0.0, 0.0, 0.0], [5.0, 3.0, 0.0], [-5.0, 3.0, 0.0], [1.0, 0.0, 0.0]]
+    )
+    for scale in (1e-16, 1e-14, 1e-12, 1e-9, 1e-3, 1.0, 1e3, 1e9):
+        vmob = VMobject().set_points(pts * scale)
+        assert vmob.width == pytest.approx(true_width * scale, rel=1e-9, abs=1e-20)
+
+    vmob = VMobject().set_points(pts * 1e-14)
+    vmob.width = 1.0
+    assert vmob.width == pytest.approx(1.0, abs=1e-9)
+
+
+def test_family_fast_path_aggregates_exact_bounds():
+    left = _wild_cubic().shift(20 * LEFT)
+    right = _wild_cubic().shift(20 * RIGHT)
+    group = VGroup(left, right)
+    assert group.get_left()[0] == pytest.approx(-20.98472845, abs=1e-6)
+    assert group.get_right()[0] == pytest.approx(21.45299067, abs=1e-6)
+    assert group.width == pytest.approx(42.4377191, abs=1e-4)