From 1df1709e49fb95ff7e4409f7b26076b1e0a6c846 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Francisco=20Manr=C3=ADquez=20Novoa?= <49853152+chopan050@users.noreply.github.com> Date: Sat, 13 Jul 2024 12:27:45 -0400 Subject: [PATCH] Rewrite `manim.utils.bezier.get_quadratic_approximation_of_cubic()` to produce curves which can be animated smoothly (#3829) * Rewrite get_quadratic_approximation_of_cubic() and add test * Move test_get... to end of file * Complete docstring for get_quadratic...() --- manim/utils/bezier.py | 205 ++++++++++++++++++++---------- tests/module/utils/test_bezier.py | 48 +++++++ 2 files changed, 187 insertions(+), 66 deletions(-) diff --git a/manim/utils/bezier.py b/manim/utils/bezier.py index c0bbcdf912..5a4049b654 100644 --- a/manim/utils/bezier.py +++ b/manim/utils/bezier.py @@ -28,7 +28,6 @@ from manim.typing import PointDType from manim.utils.simple_functions import choose -from manim.utils.space_ops import cross2d, find_intersection if TYPE_CHECKING: import numpy.typing as npt @@ -39,6 +38,7 @@ MatrixMN, Point3D, Point3D_Array, + QuadraticBezierPoints_Array, ) # l is a commonly used name in linear algebra @@ -1610,74 +1610,147 @@ def get_smooth_open_cubic_bezier_handle_points( return H1, H2 -# Given 4 control points for a cubic bezier curve (or arrays of such) -# return control points for 2 quadratics (or 2n quadratics) approximating them. +@overload def get_quadratic_approximation_of_cubic( a0: Point3D, h0: Point3D, h1: Point3D, a1: Point3D -) -> BezierPoints: - a0 = np.array(a0, ndmin=2) - h0 = np.array(h0, ndmin=2) - h1 = np.array(h1, ndmin=2) - a1 = np.array(a1, ndmin=2) - # Tangent vectors at the start and end. - T0 = h0 - a0 - T1 = a1 - h1 - - # Search for inflection points. If none are found, use the - # midpoint as a cut point. - # Based on http://www.caffeineowl.com/graphics/2d/vectorial/cubic-inflexion.html - has_infl = np.ones(len(a0), dtype=bool) - - p = h0 - a0 - q = h1 - 2 * h0 + a0 - r = a1 - 3 * h1 + 3 * h0 - a0 - - a = cross2d(q, r) - b = cross2d(p, r) - c = cross2d(p, q) - - disc = b * b - 4 * a * c - has_infl &= disc > 0 - sqrt_disc = np.sqrt(np.abs(disc)) - settings = np.seterr(all="ignore") - ti_bounds = [] - for sgn in [-1, +1]: - ti = (-b + sgn * sqrt_disc) / (2 * a) - ti[a == 0] = (-c / b)[a == 0] - ti[(a == 0) & (b == 0)] = 0 - ti_bounds.append(ti) - ti_min, ti_max = ti_bounds - np.seterr(**settings) - ti_min_in_range = has_infl & (0 < ti_min) & (ti_min < 1) - ti_max_in_range = has_infl & (0 < ti_max) & (ti_max < 1) - - # Choose a value of t which starts at 0.5, - # but is updated to one of the inflection points - # if they lie between 0 and 1 - - t_mid = 0.5 * np.ones(len(a0)) - t_mid[ti_min_in_range] = ti_min[ti_min_in_range] - t_mid[ti_max_in_range] = ti_max[ti_max_in_range] - - m, n = a0.shape - t_mid = t_mid.repeat(n).reshape((m, n)) - - # Compute bezier point and tangent at the chosen value of t (these are vectorized) - mid = bezier([a0, h0, h1, a1])(t_mid) # type: ignore - Tm = bezier([h0 - a0, h1 - h0, a1 - h1])(t_mid) # type: ignore - - # Intersection between tangent lines at end points - # and tangent in the middle - i0 = find_intersection(a0, T0, mid, Tm) - i1 = find_intersection(a1, T1, mid, Tm) - - m, n = np.shape(a0) - result = np.zeros((6 * m, n)) +) -> QuadraticBezierPoints_Array: ... + + +@overload +def get_quadratic_approximation_of_cubic( + a0: Point3D_Array, + h0: Point3D_Array, + h1: Point3D_Array, + a1: Point3D_Array, +) -> QuadraticBezierPoints_Array: ... + + +def get_quadratic_approximation_of_cubic(a0, h0, h1, a1): + r"""If ``a0``, ``h0``, ``h1`` and ``a1`` are the control points of a cubic + Bézier curve, approximate the curve with two quadratic Bézier curves and + return an array of 6 points, where the first 3 points represent the first + quadratic curve and the last 3 represent the second one. + + Otherwise, if ``a0``, ``h0``, ``h1`` and ``a1`` are _arrays_ of :math:`N` + points representing :math:`N` cubic Bézier curves, return an array of + :math:`6N` points where each group of :math:`6` consecutive points + approximates each of the :math:`N` curves in a similar way as above. + + .. note:: + If the cubic spline given by the original cubic Bézier curves is + smooth, this algorithm will generate a quadratic spline which is also + smooth. + + If a cubic Bézier is given by + + .. math:: + C(t) = (1-t)^3 A_0 + 3(1-t)^2 t H_0 + 3(1-t)t^2 H_1 + t^3 A_1 + + where :math:`A_0`, :math:`H_0`, :math:`H_1` and :math:`A_1` are its + control points, then this algorithm should generate two quadratic + Béziers given by + + .. math:: + Q_0(t) &= (1-t)^2 A_0 + 2(1-t)t M_0 + t^2 K \\ + Q_1(t) &= (1-t)^2 K + 2(1-t)t M_1 + t^2 A_1 + + where :math:`M_0` and :math:`M_1` are the respective handles to be + found for both curves, and :math:`K` is the end anchor of the 1st curve + and the start anchor of the 2nd, which must also be found. + + To solve for :math:`M_0`, :math:`M_1` and :math:`K`, three conditions + can be imposed: + + 1. :math:`Q_0'(0) = \frac{1}{2}C'(0)`. The derivative of the first + quadratic curve at :math:`t = 0` should be proportional to that of + the original cubic curve, also at :math:`t = 0`. Because the cubic + curve is split into two parts, it is necessary to divide this by + two: the speed of a point travelling through the curve should be + half of the original. This gives: + + .. math:: + Q_0'(0) &= \frac{1}{2}C'(0) \\ + 2(M_0 - A_0) &= \frac{3}{2}(H_0 - A_0) \\ + 2M_0 - 2A_0 &= \frac{3}{2}H_0 - \frac{3}{2}A_0 \\ + 2M_0 &= \frac{3}{2}H_0 + \frac{1}{2}A_0 \\ + M_0 &= \frac{1}{4}(3H_0 + A_0) + + 2. :math:`Q_1'(1) = \frac{1}{2}C'(1)`. The derivative of the second + quadratic curve at :math:`t = 1` should be half of that of the + original cubic curve for the same reasons as above, also at + :math:`t = 1`. This gives: + + .. math:: + Q_1'(1) &= \frac{1}{2}C'(1) \\ + 2(A_1 - M_1) &= \frac{3}{2}(A_1 - H_1) \\ + 2A_1 - 2M_1 &= \frac{3}{2}A_1 - \frac{3}{2}H_1 \\ + -2M_1 &= -\frac{1}{2}A_1 - \frac{3}{2}H_1 \\ + M_1 &= \frac{1}{4}(3H_1 + A_1) + + 3. :math:`Q_0'(1) = Q_1'(0)`. The derivatives of both quadratic curves + should match at the point :math:`K`, in order for the final spline + to be smooth. This gives: + + .. math:: + Q_0'(1) &= Q_1'(0) \\ + 2(K - M_0) &= 2(M_1 - K) \\ + 2K - 2M_0 &= 2M_1 - 2K \\ + 4K &= 2M_0 + 2M_1 \\ + K &= \frac{1}{2}(M_0 + M_1) + + This is sufficient to find proper control points for the quadratic + Bézier curves. + + Parameters + ---------- + a0 + The start anchor of a single cubic Bézier curve, or an array of + :math:`N` start anchors for :math:`N` curves. + h0 + The first handle of a single cubic Bézier curve, or an array of + :math:`N` first handles for :math:`N` curves. + h1 + The second handle of a single cubic Bézier curve, or an array of + :math:`N` second handles for :math:`N` curves. + a1 + The end anchor of a single cubic Bézier curve, or an array of + :math:`N` end anchors for :math:`N` curves. + + Returns + ------- + result + An array containing either 6 points for 2 quadratic Bézier curves + approximating the original cubic curve, or :math:`6N` points for + :math:`2N` quadratic curves approximating :math:`N` cubic curves. + + Raises + ------ + ValueError + If ``a0``, ``h0``, ``h1`` and ``a1`` have different dimensions, or + if their number of dimensions is not 1 or 2. + """ + a0 = np.asarray(a0) + h0 = np.asarray(h0) + h1 = np.asarray(h1) + a1 = np.asarray(a1) + + if all(arr.ndim == 1 for arr in (a0, h0, h1, a1)): + num_curves, dim = 1, a0.shape[0] + elif all(arr.ndim == 2 for arr in (a0, h0, h1, a1)): + num_curves, dim = a0.shape + else: + raise ValueError("All arguments must be Point3D or Point3D_Array.") + + m0 = 0.25 * (3 * h0 + a0) + m1 = 0.25 * (3 * h1 + a1) + k = 0.5 * (m0 + m1) + + result = np.empty((6 * num_curves, dim)) result[0::6] = a0 - result[1::6] = i0 - result[2::6] = mid - result[3::6] = mid - result[4::6] = i1 + result[1::6] = m0 + result[2::6] = k + result[3::6] = k + result[4::6] = m1 result[5::6] = a1 return result diff --git a/tests/module/utils/test_bezier.py b/tests/module/utils/test_bezier.py index dca4be193e..e7e02a89ee 100644 --- a/tests/module/utils/test_bezier.py +++ b/tests/module/utils/test_bezier.py @@ -8,6 +8,7 @@ from manim.typing import ManimFloat from manim.utils.bezier import ( _get_subdivision_matrix, + get_quadratic_approximation_of_cubic, get_smooth_cubic_bezier_handle_points, partial_bezier_points, split_bezier, @@ -167,3 +168,50 @@ def test_get_smooth_cubic_bezier_handle_points() -> None: ] ), ) + + +def test_get_quadratic_approximation_of_cubic() -> None: + C = np.array( + [ + [-5, 2, 0], + [-4, 2, 0], + [-3, 2, 0], + [-2, 2, 0], + [-2, 2, 0], + [-7 / 3, 4 / 3, 0], + [-8 / 3, 2 / 3, 0], + [-3, 0, 0], + [-3, 0, 0], + [-1 / 3, -1, 0], + [7 / 3, -2, 0], + [5, -3, 0], + ] + ) + a0, h0, h1, a1 = C[::4], C[1::4], C[2::4], C[3::4] + + Q = get_quadratic_approximation_of_cubic(a0, h0, h1, a1) + assert np.allclose( + Q, + np.array( + [ + [-5, 2, 0], + [-17 / 4, 2, 0], + [-7 / 2, 2, 0], + [-7 / 2, 2, 0], + [-11 / 4, 2, 0], + [-2, 2, 0], + [-2, 2, 0], + [-9 / 4, 3 / 2, 0], + [-5 / 2, 1, 0], + [-5 / 2, 1, 0], + [-11 / 4, 1 / 2, 0], + [-3, 0, 0], + [-3, 0, 0], + [-1, -3 / 4, 0], + [1, -3 / 2, 0], + [1, -3 / 2, 0], + [3, -9 / 4, 0], + [5, -3, 0], + ] + ), + )