diff --git a/manim/mobject/opengl/opengl_vectorized_mobject.py b/manim/mobject/opengl/opengl_vectorized_mobject.py index d668e3a2b1..41d6aac084 100644 --- a/manim/mobject/opengl/opengl_vectorized_mobject.py +++ b/manim/mobject/opengl/opengl_vectorized_mobject.py @@ -15,13 +15,13 @@ from manim.renderer.shader_wrapper import ShaderWrapper from manim.utils.bezier import ( bezier, + bezier_remap, get_quadratic_approximation_of_cubic, get_smooth_cubic_bezier_handle_points, integer_interpolate, interpolate, - partial_quadratic_bezier_points, + partial_bezier_points, proportions_along_bezier_curve_for_point, - quadratic_bezier_remap, ) from manim.utils.color import BLACK, WHITE, ManimColor, ParsableManimColor from manim.utils.config_ops import _Data @@ -555,7 +555,7 @@ def subdivide_sharp_curves(self, angle_threshold=30 * DEGREES, recurse=True): alphas = np.linspace(0, 1, n + 1) new_points.extend( [ - partial_quadratic_bezier_points(tup, a1, a2) + partial_bezier_points(tup, a1, a2) for a1, a2 in zip(alphas, alphas[1:]) ], ) @@ -1275,33 +1275,12 @@ def insert_n_curves_to_point_list(self, n: int, points: np.ndarray) -> np.ndarra if len(points) == 1: return np.repeat(points, nppc * n, 0) - bezier_groups = self.get_bezier_tuples_from_points(points) - norms = np.array([np.linalg.norm(bg[nppc - 1] - bg[0]) for bg in bezier_groups]) - total_norm = sum(norms) - # Calculate insertions per curve (ipc) - if total_norm < 1e-6: - ipc = [n] + [0] * (len(bezier_groups) - 1) - else: - ipc = np.round(n * norms / sum(norms)).astype(int) - - diff = n - sum(ipc) - for _ in range(diff): - ipc[np.argmin(ipc)] += 1 - for _ in range(-diff): - ipc[np.argmax(ipc)] -= 1 - - new_length = sum(x + 1 for x in ipc) - new_points = np.empty((new_length, nppc, 3)) - i = 0 - for group, n_inserts in zip(bezier_groups, ipc): - # What was once a single quadratic curve defined - # by "group" will now be broken into n_inserts + 1 - # smaller quadratic curves - alphas = np.linspace(0, 1, n_inserts + 2) - for a1, a2 in zip(alphas, alphas[1:]): - new_points[i] = partial_quadratic_bezier_points(group, a1, a2) - i = i + 1 - return np.vstack(new_points) + bezier_tuples = self.get_bezier_tuples_from_points(points) + current_number_of_curves = len(bezier_tuples) + new_number_of_curves = current_number_of_curves + n + new_bezier_tuples = bezier_remap(bezier_tuples, new_number_of_curves) + new_points = new_bezier_tuples.reshape(-1, 3) + return new_points def interpolate(self, mobject1, mobject2, alpha, *args, **kwargs): super().interpolate(mobject1, mobject2, alpha, *args, **kwargs) @@ -1354,7 +1333,7 @@ def pointwise_become_partial( return self if lower_index == upper_index: self.append_points( - partial_quadratic_bezier_points( + partial_bezier_points( bezier_triplets[lower_index], lower_residue, upper_residue, @@ -1362,24 +1341,18 @@ def pointwise_become_partial( ) else: self.append_points( - partial_quadratic_bezier_points( - bezier_triplets[lower_index], lower_residue, 1 - ), + partial_bezier_points(bezier_triplets[lower_index], lower_residue, 1), ) inner_points = bezier_triplets[lower_index + 1 : upper_index] if len(inner_points) > 0: if remap: - new_triplets = quadratic_bezier_remap( - inner_points, num_quadratics - 2 - ) + new_triplets = bezier_remap(inner_points, num_quadratics - 2) else: new_triplets = bezier_triplets self.append_points(np.asarray(new_triplets).reshape(-1, 3)) self.append_points( - partial_quadratic_bezier_points( - bezier_triplets[upper_index], 0, upper_residue - ), + partial_bezier_points(bezier_triplets[upper_index], 0, upper_residue), ) return self diff --git a/manim/mobject/types/vectorized_mobject.py b/manim/mobject/types/vectorized_mobject.py index cf7ce9a6b4..deb3690620 100644 --- a/manim/mobject/types/vectorized_mobject.py +++ b/manim/mobject/types/vectorized_mobject.py @@ -30,6 +30,7 @@ ) from manim.utils.bezier import ( bezier, + bezier_remap, get_smooth_handle_points, integer_interpolate, interpolate, @@ -1693,40 +1694,11 @@ def insert_n_curves_to_point_list( if len(points) == 1: nppcc = self.n_points_per_cubic_curve return np.repeat(points, nppcc * n, 0) - bezier_quads = self.get_cubic_bezier_tuples_from_points(points) - curr_num = len(bezier_quads) - target_num = curr_num + n - # This is an array with values ranging from 0 - # up to curr_num, with repeats such that - # it's total length is target_num. For example, - # with curr_num = 10, target_num = 15, this would - # be [0, 0, 1, 2, 2, 3, 4, 4, 5, 6, 6, 7, 8, 8, 9] - repeat_indices = (np.arange(target_num, dtype="i") * curr_num) // target_num - - # If the nth term of this list is k, it means - # that the nth curve of our path should be split - # into k pieces. - # In the above example our array had the following elements - # [0, 0, 1, 2, 2, 3, 4, 4, 5, 6, 6, 7, 8, 8, 9] - # We have two 0s, one 1, two 2s and so on. - # The split factors array would hence be: - # [2, 1, 2, 1, 2, 1, 2, 1, 2, 1] - split_factors = np.zeros(curr_num, dtype="i") - for val in repeat_indices: - split_factors[val] += 1 - - new_points = np.zeros((0, self.dim)) - for quad, sf in zip(bezier_quads, split_factors): - # What was once a single cubic curve defined - # by "quad" will now be broken into sf - # smaller cubic curves - alphas = np.linspace(0, 1, sf + 1) - for a1, a2 in zip(alphas, alphas[1:]): - new_points = np.append( - new_points, - partial_bezier_points(quad, a1, a2), - axis=0, - ) + bezier_tuples = self.get_cubic_bezier_tuples_from_points(points) + current_number_of_curves = len(bezier_tuples) + new_number_of_curves = current_number_of_curves + n + new_bezier_tuples = bezier_remap(bezier_tuples, new_number_of_curves) + new_points = new_bezier_tuples.reshape(-1, 3) return new_points def align_rgbas(self, vmobject: VMobject) -> Self: diff --git a/manim/utils/bezier.py b/manim/utils/bezier.py index 05e404248e..709f8b3bf4 100644 --- a/manim/utils/bezier.py +++ b/manim/utils/bezier.py @@ -16,7 +16,9 @@ __all__ = [ "bezier", "partial_bezier_points", - "partial_quadratic_bezier_points", + "split_bezier", + "subdivide_bezier", + "bezier_remap", "interpolate", "integer_interpolate", "mid", @@ -86,12 +88,120 @@ def bezier( ) -# !TODO: This function has still a weird implementation with the overlapping points def partial_bezier_points(points: BezierPoints, a: float, b: float) -> BezierPoints: - """Given an array of points which define bezier curve, and two numbers 0<=a BezierPoi Returns ------- - np.ndarray - Set of points defining the partial bezier curve. + :class:`~.BezierPoints` + An array containing the control points defining the partial Bézier curve. """ - _len = len(points) + # Border cases if a == 1: - return np.asarray([points[-1]] * _len, dtype=PointDType) + arr = np.array(points) + arr[:] = arr[-1] + return arr + if b == 0: + arr = np.array(points) + arr[:] = arr[0] + return arr - a_to_1 = np.asarray( - [bezier(points[i:])(a) for i in range(_len)], - dtype=PointDType, - ) - end_prop = (b - a) / (1.0 - a) - return np.asarray( - [bezier(a_to_1[: i + 1])(end_prop) for i in range(_len)], - dtype=PointDType, - ) + points = np.asarray(points) + degree = points.shape[0] - 1 + if degree == 3: + ma, mb = 1 - a, 1 - b + a2, b2, ma2, mb2 = a * a, b * b, ma * ma, mb * mb + a3, b3, ma3, mb3 = a2 * a, b2 * b, ma2 * ma, mb2 * mb -# Shortened version of partial_bezier_points just for quadratics, -# since this is called a fair amount -def partial_quadratic_bezier_points( - points: QuadraticBezierPoints, a: float, b: float -) -> QuadraticBezierPoints: - if a == 1: - return np.asarray(3 * [points[-1]]) + portion_matrix = np.array( + [ + [ma3, 3 * ma2 * a, 3 * ma * a2, a3], + [ma2 * mb, 2 * ma * a * mb + ma2 * b, a2 * mb + 2 * ma * a * b, a2 * b], + [ma * mb2, a * mb2 + 2 * ma * mb * b, 2 * a * mb * b + ma * b2, a * b2], + [mb3, 3 * mb2 * b, 3 * mb * b2, b3], + ] + ) + return portion_matrix @ points + + if degree == 2: + ma, mb = 1 - a, 1 - b - def curve(t: float) -> Point3D: - return np.asarray( - points[0] * (1 - t) * (1 - t) - + 2 * points[1] * t * (1 - t) - + points[2] * t * t + portion_matrix = np.array( + [ + [ma * ma, 2 * a * ma, a * a], + [ma * mb, a * mb + ma * b, a * b], + [mb * mb, 2 * b * mb, b * b], + ] ) + return portion_matrix @ points - # bezier(points) - h0 = curve(a) if a > 0 else points[0] - h2 = curve(b) if b < 1 else points[2] - h1_prime = (1 - a) * points[1] + a * points[2] - end_prop = (b - a) / (1.0 - a) - h1 = (1 - end_prop) * h0 + end_prop * h1_prime - return np.asarray((h0, h1, h2)) + if degree == 1: + direction = points[1] - points[0] + return np.array( + [ + points[0] + a * direction, + points[0] + b * direction, + ] + ) + + if degree == 0: + return points + + # Fallback case for nth degree Béziers + # It is convenient that np.array copies points + arr = np.array(points, dtype=float) + N = arr.shape[0] + + # Current state for an example Bézier curve C0 = [P0, P1, P2, P3]: + # arr = [P0, P1, P2, P3] + if a != 0: + for i in range(1, N): + # 1st iter: arr = [L0(a), L1(a), L2(a), P3] + # 2nd iter: arr = [Q0(a), Q1(a), L2(a), P3] + # 3rd iter: arr = [C0(a), Q1(a), L2(a), P3] + arr[: N - i] += a * (arr[1 : N - i + 1] - arr[: N - i]) + + # For faster calculations we shall define mu = 1 - u = (1 - b) / (1 - a). + # This is because: + # L0'(u) = P0' + u(P1' - P0') + # = (1-u)P0' + uP1' + # = muP0' + (1-mu)P1' + # = P1' + mu(P0' - P1) + # In this way, one can do something similar to the first loop. + # + # Current state: + # arr = [C0(a), Q1(a), L2(a), P3] + # = [P0', P1', P2', P3'] + if b != 1: + mu = (1 - b) / (1 - a) + for i in range(1, N): + # 1st iter: arr = [P0', L0'(u), L1'(u), L2'(u)] + # 2nd iter: arr = [P0', L0'(u), Q0'(u), Q1'(u)] + # 3rd iter: arr = [P0', L0'(u), Q0'(u), C0'(u)] + arr[i:] += mu * (arr[i - 1 : -1] - arr[i:]) + + return arr + + +def split_bezier(points: BezierPoints, t: float) -> Point3D_Array: + r"""Split a Bézier curve at argument ``t`` into two curves. + .. note:: -def split_quadratic_bezier(points: QuadraticBezierPoints, t: float) -> BezierPoints: - """Split a quadratic Bézier curve at argument ``t`` into two quadratic curves. + .. seealso:: + `A Primer on Bézier Curves #10: Splitting curves. Pomax. `_ + + As an example for a cubic Bézier curve, let :math:`p_0, p_1, p_2, p_3` be the points + needed for the curve :math:`C_0 = [p_0, \ p_1, \ p_2, \ p_3]`. + + Define the 3 linear Béziers :math:`L_0, L_1, L_2` as interpolations of :math:`p_0, p_1, p_2, p_3`: + + .. math:: + L_0(t) &= p_0 + t(p_1 - p_0) \\ + L_1(t) &= p_1 + t(p_2 - p_1) \\ + L_2(t) &= p_2 + t(p_3 - p_2) + + Define the 2 quadratic Béziers :math:`Q_0, Q_1` as interpolations of :math:`L_0, L_1, L_2`: + + .. math:: + Q_0(t) &= L_0(t) + t(L_1(t) - L_0(t)) \\ + Q_1(t) &= L_1(t) + t(L_2(t) - L_1(t)) + + Then :math:`C_0` is the following interpolation of :math:`Q_0` and :math:`Q_1`: + + .. math:: + C_0(t) = Q_0(t) + t(Q_1(t) - Q_0(t)) + + Evaluating :math:`C_0` at a value :math:`t=t'` splits :math:`C_0` into two cubic Béziers :math:`H_0` + and :math:`H_1`, defined by some of the points we calculated earlier: + + .. math:: + H_0 &= [p_0, &\ L_0(t'), &\ Q_0(t'), &\ C_0(t') &] \\ + H_1 &= [p_0(t'), &\ Q_1(t'), &\ L_2(t'), &\ p_3 &] + + As the resulting curves are obtained from linear combinations of ``points``, everything can + be encoded into a matrix for efficiency, which is done for Bézier curves of degree up to 3. + + .. seealso:: + `A Primer on Bézier Curves #11: Splitting curves using matrices. Pomax. `_ + + For the simpler case of a quadratic Bézier curve: + + .. math:: + H_0 + &= + \begin{pmatrix} + p_0 \\ + (1-t) p_0 + t p_1 \\ + (1-t)^2 p_0 + 2(1-t)t p_1 + t^2 p_2 \\ + \end{pmatrix} + &= + \begin{pmatrix} + 1 & 0 & 0 \\ + (1-t) & t & 0\\ + (1-t)^2 & 2(1-t)t & t^2 + \end{pmatrix} + \begin{pmatrix} + p_0 \\ + p_1 \\ + p_2 + \end{pmatrix} + \\ + & + \\ + H_1 + &= + \begin{pmatrix} + (1-t)^2 p_0 + 2(1-t)t p_1 + t^2 p_2 \\ + (1-t) p_1 + t p_2 \\ + p_2 + \end{pmatrix} + &= + \begin{pmatrix} + (1-t)^2 & 2(1-t)t & t^2 \\ + 0 & (1-t) & t \\ + 0 & 0 & 1 + \end{pmatrix} + \begin{pmatrix} + p_0 \\ + p_1 \\ + p_2 + \end{pmatrix} + + from where one can define a :math:`(6, 3)` split matrix :math:`S_2` which can multiply + the array of ``points`` to compute the return value: + + .. math:: + S_2 + &= + \begin{pmatrix} + 1 & 0 & 0 \\ + (1-t) & t & 0 \\ + (1-t)^2 & 2(1-t)t & t^2 \\ + (1-t)^2 & 2(1-t)t & t^2 \\ + 0 & (1-t) & t \\ + 0 & 0 & 1 + \end{pmatrix} + \\ + & + \\ + S_2 P + &= + \begin{pmatrix} + 1 & 0 & 0 \\ + (1-t) & t & 0 \\ + (1-t)^2 & 2(1-t)t & t^2 \\ + (1-t)^2 & 2(1-t)t & t^2 \\ + 0 & (1-t) & t \\ + 0 & 0 & 1 + \end{pmatrix} + \begin{pmatrix} + p_0 \\ + p_1 \\ + p_2 + \end{pmatrix} + = + \begin{pmatrix} + \vert \\ + H_0 \\ + \vert \\ + \vert \\ + H_1 \\ + \vert + \end{pmatrix} + + For the previous example with a cubic Bézier curve: + + .. math:: + H_0 + &= + \begin{pmatrix} + p_0 \\ + (1-t) p_0 + t p_1 \\ + (1-t)^2 p_0 + 2(1-t)t p_1 + t^2 p_2 \\ + (1-t)^3 p_0 + 3(1-t)^2 t p_1 + 3(1-t)t^2 p_2 + t^3 p_3 + \end{pmatrix} + &= + \begin{pmatrix} + 1 & 0 & 0 & 0 \\ + (1-t) & t & 0 & 0 \\ + (1-t)^2 & 2(1-t)t & t^2 & 0 \\ + (1-t)^3 & 3(1-t)^2 t & 3(1-t)t^2 & t^3 + \end{pmatrix} + \begin{pmatrix} + p_0 \\ + p_1 \\ + p_2 \\ + p_3 + \end{pmatrix} + \\ + & + \\ + H_1 + &= + \begin{pmatrix} + (1-t)^3 p_0 + 3(1-t)^2 t p_1 + 3(1-t)t^2 p_2 + t^3 p_3 \\ + (1-t)^2 p_1 + 2(1-t)t p_2 + t^2 p_3 \\ + (1-t) p_2 + t p_3 \\ + p_3 + \end{pmatrix} + &= + \begin{pmatrix} + (1-t)^3 & 3(1-t)^2 t & 3(1-t)t^2 & t^3 \\ + 0 & (1-t)^2 & 2(1-t)t & t^2 \\ + 0 & 0 & (1-t) & t \\ + 0 & 0 & 0 & 1 + \end{pmatrix} + \begin{pmatrix} + p_0 \\ + p_1 \\ + p_2 \\ + p_3 + \end{pmatrix} + + from where one can define a :math:`(8, 4)` split matrix :math:`S_3` which can multiply + the array of ``points`` to compute the return value: + + .. math:: + S_3 + &= + \begin{pmatrix} + 1 & 0 & 0 & 0 \\ + (1-t) & t & 0 & 0 \\ + (1-t)^2 & 2(1-t)t & t^2 & 0 \\ + (1-t)^3 & 3(1-t)^2 t & 3(1-t)t^2 & t^3 \\ + (1-t)^3 & 3(1-t)^2 t & 3(1-t)t^2 & t^3 \\ + 0 & (1-t)^2 & 2(1-t)t & t^2 \\ + 0 & 0 & (1-t) & t \\ + 0 & 0 & 0 & 1 + \end{pmatrix} + \\ + & + \\ + S_3 P + &= + \begin{pmatrix} + 1 & 0 & 0 & 0 \\ + (1-t) & t & 0 & 0 \\ + (1-t)^2 & 2(1-t)t & t^2 & 0 \\ + (1-t)^3 & 3(1-t)^2 t & 3(1-t)t^2 & t^3 \\ + (1-t)^3 & 3(1-t)^2 t & 3(1-t)t^2 & t^3 \\ + 0 & (1-t)^2 & 2(1-t)t & t^2 \\ + 0 & 0 & (1-t) & t \\ + 0 & 0 & 0 & 1 + \end{pmatrix} + \begin{pmatrix} + p_0 \\ + p_1 \\ + p_2 \\ + p_3 + \end{pmatrix} + = + \begin{pmatrix} + \vert \\ + H_0 \\ + \vert \\ + \vert \\ + H_1 \\ + \vert + \end{pmatrix} Parameters ---------- points - The control points of the bezier curve - has shape ``[a1, h1, b1]`` + The control points of the Bézier curve. t - The ``t``-value at which to split the Bézier curve + The ``t``-value at which to split the Bézier curve. Returns ------- - The two Bézier curves as a list of tuples, - has the shape ``[a1, h1, b1], [a2, h2, b2]`` + :class:`~.Point3D_Array` + An array containing the control points defining the two Bézier curves. """ - a1, h1, a2 = points - s1 = interpolate(a1, h1, t) - s2 = interpolate(h1, a2, t) - p = interpolate(s1, s2, t) - return np.array((a1, s1, p, p, s2, a2)) + points = np.asarray(points) + N, dim = points.shape + degree = N - 1 + + if degree == 3: + mt = 1 - t + mt2 = mt * mt + mt3 = mt2 * mt + t2 = t * t + t3 = t2 * t + two_mt_t = 2 * mt * t + three_mt2_t = 3 * mt2 * t + three_mt_t2 = 3 * mt * t2 + + # Split matrix S3 explained in the docstring + split_matrix = np.array( + [ + [1, 0, 0, 0], + [mt, t, 0, 0], + [mt2, two_mt_t, t2, 0], + [mt3, three_mt2_t, three_mt_t2, t3], + [mt3, three_mt2_t, three_mt_t2, t3], + [0, mt2, two_mt_t, t2], + [0, 0, mt, t], + [0, 0, 0, 1], + ] + ) + + return split_matrix @ points + + if degree == 2: + mt = 1 - t + mt2 = mt * mt + t2 = t * t + two_tmt = 2 * t * mt + + # Split matrix S2 explained in the docstring + split_matrix = np.array( + [ + [1, 0, 0], + [mt, t, 0], + [mt2, two_tmt, t2], + [mt2, two_tmt, t2], + [0, mt, t], + [0, 0, 1], + ] + ) + + return split_matrix @ points + + if degree == 1: + middle = points[0] + t * (points[1] - points[0]) + return np.array([points[0], middle, middle, points[1]]) + + if degree == 0: + return np.array([points[0], points[0]]) + + # Fallback case for nth degree Béziers + arr = np.empty((2, N, dim)) + arr[1] = points + arr[0, 0] = points[0] + + # Example for a cubic Bézier + # arr[0] = [P0 .. .. ..] + # arr[1] = [P0 P1 P2 P3] + for i in range(1, N): + # 1st iter: arr[1] = [L0 L1 L2 P3] + # 2nd iter: arr[1] = [Q0 Q1 L2 P3] + # 3rd iter: arr[1] = [C0 Q1 L2 P3] + arr[1, : N - i] += t * (arr[1, 1 : N - i + 1] - arr[1, : N - i]) + # 1st iter: arr[0] = [P0 L0 .. ..] + # 2nd iter: arr[0] = [P0 L0 Q0 ..] + # 3rd iter: arr[0] = [P0 L0 Q0 C0] + arr[0, i] = arr[1, 0] + + return arr.reshape(2 * N, dim) + +# Memos explained in subdivide_bezier docstring +SUBDIVISION_MATRICES = [{} for i in range(4)] -def subdivide_quadratic_bezier(points: QuadraticBezierPoints, n: int) -> BezierPoints: - """Subdivide a quadratic Bézier curve into ``n`` subcurves which have the same shape. + +def _get_subdivision_matrix(n_points: int, n_divisions: int) -> MatrixMN: + """Gets the matrix which subdivides a Bézier curve of + ``n_points`` control points into ``n_divisions`` parts. + + Auxiliary function for :func:`subdivide_bezier`. See its + docstrings for an explanation of the matrix build process. + + Parameters + ---------- + n_points + The number of control points of the Bézier curve to + subdivide. This function only handles up to 4 points. + n_divisions + The number of parts to subdivide the Bézier curve into. + + Returns + ------- + MatrixMN + The matrix which, upon multiplying the control points of the + Bézier curve, subdivides it into ``n_divisions`` parts. + """ + if n_points not in (1, 2, 3, 4): + raise NotImplementedError( + "This function does not support subdividing Bézier " + "curves with 0 or more than 4 control points." + ) + + subdivision_matrix = SUBDIVISION_MATRICES[n_points - 1].get(n_divisions, None) + if subdivision_matrix is not None: + return subdivision_matrix + + subdivision_matrix = np.empty((n_points * n_divisions, n_points)) + + # Cubic Bézier + if n_points == 4: + for i in range(n_divisions): + i2 = i * i + i3 = i2 * i + ip1 = i + 1 + ip12 = ip1 * ip1 + ip13 = ip12 * ip1 + nmi = n_divisions - i + nmi2 = nmi * nmi + nmi3 = nmi2 * nmi + nmim1 = nmi - 1 + nmim12 = nmim1 * nmim1 + nmim13 = nmim12 * nmim1 + + subdivision_matrix[4 * i : 4 * (i + 1)] = np.array( + [ + [ + nmi3, + 3 * nmi2 * i, + 3 * nmi * i2, + i3, + ], + [ + nmi2 * nmim1, + 2 * nmi * nmim1 * i + nmi2 * ip1, + nmim1 * i2 + 2 * nmi * i * ip1, + i2 * ip1, + ], + [ + nmi * nmim12, + nmim12 * i + 2 * nmi * nmim1 * ip1, + 2 * nmim1 * i * ip1 + nmi * ip12, + i * ip12, + ], + [ + nmim13, + 3 * nmim12 * ip1, + 3 * nmim1 * ip12, + ip13, + ], + ] + ) + subdivision_matrix /= n_divisions * n_divisions * n_divisions + + # Quadratic Bézier + elif n_points == 3: + for i in range(n_divisions): + ip1 = i + 1 + nmi = n_divisions - i + nmim1 = nmi - 1 + subdivision_matrix[3 * i : 3 * (i + 1)] = np.array( + [ + [nmi * nmi, 2 * i * nmi, i * i], + [nmi * nmim1, i * nmim1 + ip1 * nmi, i * ip1], + [nmim1 * nmim1, 2 * ip1 * nmim1, ip1 * ip1], + ] + ) + subdivision_matrix /= n_divisions * n_divisions + + # Linear Bézier (straight line) + elif n_points == 2: + aux_range = np.arange(n_divisions + 1) + subdivision_matrix[::2, 1] = aux_range[:-1] + subdivision_matrix[1::2, 1] = aux_range[1:] + subdivision_matrix[:, 0] = subdivision_matrix[::-1, 1] + subdivision_matrix /= n_divisions + + # Zero-degree Bézier (single point) + elif n_points == 1: + subdivision_matrix[:] = 1 + + SUBDIVISION_MATRICES[n_points - 1][n_divisions] = subdivision_matrix + return subdivision_matrix + + +def subdivide_bezier(points: BezierPoints, n_divisions: int) -> Point3D_Array: + r"""Subdivide a Bézier curve into :math:`n` subcurves which have the same shape. The points at which the curve is split are located at the - arguments :math:`t = i/n` for :math:`i = 1, ..., n-1`. + arguments :math:`t = \frac{i}{n}`, for :math:`i \in \{1, ..., n-1\}`. + + .. seealso:: + + * See :func:`split_bezier` for an explanation on how to split Bézier curves. + * See :func:`partial_bezier_points` for an extra understanding of this function. + + + .. note:: + The resulting subcurves can be expressed as linear combinations of + ``points``, which can be encoded in a single matrix that is precalculated + for 2nd and 3rd degree Bézier curves. + + As an example for a quadratic Bézier curve: taking inspiration from the + explanation in :func:`partial_bezier_points`, where the following matrix + :math:`P_2` was defined to extract the portion of a quadratic Bézier + curve for :math:`t \in [a, b]`: + + .. math:: + P_2 + = + \begin{pmatrix} + (1-a)^2 & 2(1-a)a & a^2 \\ + (1-a)(1-b) & a(1-b) + (1-a)b & ab \\ + (1-b)^2 & 2(1-b)b & b^2 + \end{pmatrix} + + the plan is to replace :math:`[a, b]` with + :math:`\left[ \frac{i-1}{n}, \frac{i}{n} \right], \ \forall i \in \{1, ..., n\}`. + + As an example for :math:`n = 2` divisions, construct :math:`P_1` for + the interval :math:`\left[ 0, \frac{1}{2} \right]`, and :math:`P_2` for the + interval :math:`\left[ \frac{1}{2}, 1 \right]`: + + .. math:: + P_1 + = + \begin{pmatrix} + 1 & 0 & 0 \\ + 0.5 & 0.5 & 0 \\ + 0.25 & 0.5 & 0.25 + \end{pmatrix} + , + \quad + P_2 + = + \begin{pmatrix} + 0.25 & 0.5 & 0.25 \\ + 0 & 0.5 & 0.5 \\ + 0 & 0 & 1 + \end{pmatrix} + + Therefore, the following :math:`(6, 3)` subdivision matrix :math:`D_2` can be + constructed, which will subdivide an array of ``points`` into 2 parts: + + .. math:: + D_2 + = + \begin{pmatrix} + M_1 \\ + M_2 + \end{pmatrix} + = + \begin{pmatrix} + 1 & 0 & 0 \\ + 0.5 & 0.5 & 0 \\ + 0.25 & 0.5 & 0.25 \\ + 0.25 & 0.5 & 0.25 \\ + 0 & 0.5 & 0.5 \\ + 0 & 0 & 1 + \end{pmatrix} + + For quadratic and cubic Bézier curves, the subdivision matrices are memoized for + efficiency. For higher degree curves, an iterative algorithm inspired by the + one from :func:`split_bezier` is used instead. + + .. image:: /_static/bezier_subdivision_example.png Parameters ---------- points - The control points of the Bézier curve in form ``[a1, h1, b1]`` + The control points of the Bézier curve. - n + n_divisions The number of curves to subdivide the Bézier curve into Returns ------- - The new points for the Bézier curve in the form ``[a1, h1, b1, a2, h2, b2, ...]`` - - .. image:: /_static/bezier_subdivision_example.png - + :class:`~.Point3D_Array` + An array containing the points defining the new :math:`n` subcurves. """ - beziers = np.empty((n, 3, 3)) - current = points - for j in range(0, n): - i = n - j - tmp = split_quadratic_bezier(current, 1 / i) - beziers[j] = tmp[:3] - current = tmp[3:] - return beziers.reshape(-1, 3) + if n_divisions == 1: + return points - -def quadratic_bezier_remap( - triplets: QuadraticBezierPoints_Array, new_number_of_curves: int -) -> QuadraticBezierPoints_Array: - """Remaps the number of curves to a higher amount by splitting bezier curves + points = np.asarray(points) + N, dim = points.shape + + if N <= 4: + subdivision_matrix = _get_subdivision_matrix(N, n_divisions) + return subdivision_matrix @ points + + # Fallback case for an nth degree Bézier: successive splitting + beziers = np.empty((n_divisions, N, dim)) + beziers[-1] = points + for curve_num in range(n_divisions - 1, 0, -1): + curr = beziers[curve_num] + prev = beziers[curve_num - 1] + prev[0] = curr[0] + a = (n_divisions - curve_num) / (n_divisions - curve_num + 1) + # Current state for an example cubic Bézier curve: + # prev = [P0 .. .. ..] + # curr = [P0 P1 P2 P3] + for i in range(1, N): + # 1st iter: curr = [L0 L1 L2 P3] + # 2nd iter: curr = [Q0 Q1 L2 P3] + # 3rd iter: curr = [C0 Q1 L2 P3] + curr[: N - i] += a * (curr[1 : N - i + 1] - curr[: N - i]) + # 1st iter: prev = [P0 L0 .. ..] + # 2nd iter: prev = [P0 L0 Q0 ..] + # 3rd iter: prev = [P0 L0 Q0 C0] + prev[i] = curr[0] + + return beziers.reshape(n_divisions * N, dim) + + +def bezier_remap( + bezier_tuples: BezierPoints_Array, + new_number_of_curves: int, +) -> BezierPoints_Array: + """Subdivides each curve in ``bezier_tuples`` into as many parts as necessary, until the final number of + curves reaches a desired amount, ``new_number_of_curves``. Parameters ---------- - triplets - The triplets of the quadratic bezier curves to be remapped shape(n, 3, 3) + bezier_tuples + An array of multiple Bézier curves of degree :math:`d` to be remapped. The shape of this array + must be ``(current_number_of_curves, nppc, dim)``, where: + * ``current_number_of_curves`` is the current amount of curves in the array ``bezier_tuples``, + * ``nppc`` is the amount of points per curve, such that their degree is ``nppc-1``, and + * ``dim`` is the dimension of the points, usually :math:`3`. new_number_of_curves The number of curves that the output will contain. This needs to be higher than the current number. Returns ------- - The new triplets for the quadratic bezier curves. + :class:`~.BezierPoints_Array` + The new array of shape ``(new_number_of_curves, nppc, dim)``, + containing the new Bézier curves after the remap. """ - difference = new_number_of_curves - len(triplets) - if difference <= 0: - return triplets - new_triplets = np.zeros((new_number_of_curves, 3, 3)) - idx = 0 - for triplet in triplets: - if difference > 0: - tmp_noc = int(np.ceil(difference / len(triplets))) + 1 - tmp = subdivide_quadratic_bezier(triplet, tmp_noc).reshape(-1, 3, 3) - for i in range(tmp_noc): - new_triplets[idx + i] = tmp[i] - difference -= tmp_noc - 1 - idx += tmp_noc - else: - new_triplets[idx] = triplet - idx += 1 - return new_triplets - - """ - This is an alternate version of the function just for documentation purposes - -------- + bezier_tuples = np.asarray(bezier_tuples) + current_number_of_curves, nppc, dim = bezier_tuples.shape + # This is an array with values ranging from 0 + # up to curr_num_curves, with repeats such that + # its total length is target_num_curves. For example, + # with curr_num_curves = 10, target_num_curves = 15, this + # would be [0, 0, 1, 2, 2, 3, 4, 4, 5, 6, 6, 7, 8, 8, 9]. + repeat_indices = ( + np.arange(new_number_of_curves, dtype="i") * current_number_of_curves + ) // new_number_of_curves + + # If the nth term of this list is k, it means + # that the nth curve of our path should be split + # into k pieces. + # In the above example our array had the following elements + # [0, 0, 1, 2, 2, 3, 4, 4, 5, 6, 6, 7, 8, 8, 9] + # We have two 0s, one 1, two 2s and so on. + # The split factors array would hence be: + # [2, 1, 2, 1, 2, 1, 2, 1, 2, 1] + split_factors = np.zeros(current_number_of_curves, dtype="i") + np.add.at(split_factors, repeat_indices, 1) + + new_tuples = np.empty((new_number_of_curves, nppc, dim)) + index = 0 + for curve, sf in zip(bezier_tuples, split_factors): + new_tuples[index : index + sf] = subdivide_bezier(curve, sf).reshape( + sf, nppc, dim + ) + index += sf - difference = new_number_of_curves - len(triplets) - if difference <= 0: - return triplets - new_triplets = [] - for triplet in triplets: - if difference > 0: - tmp_noc = int(np.ceil(difference / len(triplets))) + 1 - tmp = subdivide_quadratic_bezier(triplet, tmp_noc).reshape(-1, 3, 3) - for i in range(tmp_noc): - new_triplets.append(tmp[i]) - difference -= tmp_noc - 1 - else: - new_triplets.append(triplet) - return new_triplets - """ + return new_tuples # Linear interpolation variants diff --git a/tests/module/utils/_split_matrices.py b/tests/module/utils/_split_matrices.py new file mode 100644 index 0000000000..70c87497e6 --- /dev/null +++ b/tests/module/utils/_split_matrices.py @@ -0,0 +1,215 @@ +import numpy as np + +# Defined because pre-commit is inserting an unacceptable line-break +# between the "1" (or "2") and the "/ 3" +one_third = 1 / 3 +two_thirds = 2 / 3 + + +# Expected values for matrices in split_bezier +SPLIT_MATRICES = { + # For 0-degree Béziers + 0: { + 0: np.array([[1], [1]]), + one_third: np.array([[1], [1]]), + two_thirds: np.array([[1], [1]]), + 1: np.array([[1], [1]]), + }, + # For linear Béziers + 1: { + 0: np.array( + [ + [1, 0], + [1, 0], + [1, 0], + [0, 1], + ] + ), + one_third: np.array( + [ + [3, 0], + [2, 1], + [2, 1], + [0, 3], + ] + ) + / 3, + two_thirds: np.array( + [ + [3, 0], + [1, 2], + [1, 2], + [0, 3], + ] + ) + / 3, + 1: np.array( + [ + [1, 0], + [0, 1], + [0, 1], + [0, 1], + ] + ), + }, + # For quadratic Béziers + 2: { + 0: np.array( + [ + [1, 0, 0], + [1, 0, 0], + [1, 0, 0], + [1, 0, 0], + [0, 1, 0], + [0, 0, 1], + ] + ), + one_third: np.array( + [ + [9, 0, 0], + [6, 3, 0], + [4, 4, 1], + [4, 4, 1], + [0, 6, 3], + [0, 0, 9], + ] + ) + / 9, + two_thirds: np.array( + [ + [9, 0, 0], + [3, 6, 0], + [1, 4, 4], + [1, 4, 4], + [0, 3, 6], + [0, 0, 9], + ] + ) + / 9, + 1: np.array( + [ + [1, 0, 0], + [0, 1, 0], + [0, 0, 1], + [0, 0, 1], + [0, 0, 1], + [0, 0, 1], + ] + ), + }, + # For cubic Béziers + 3: { + 0: np.array( + [ + [1, 0, 0, 0], + [1, 0, 0, 0], + [1, 0, 0, 0], + [1, 0, 0, 0], + [1, 0, 0, 0], + [0, 1, 0, 0], + [0, 0, 1, 0], + [0, 0, 0, 1], + ] + ), + one_third: np.array( + [ + [27, 0, 0, 0], + [18, 9, 0, 0], + [12, 12, 3, 0], + [8, 12, 6, 1], + [8, 12, 6, 1], + [0, 12, 12, 3], + [0, 0, 18, 9], + [0, 0, 0, 27], + ] + ) + / 27, + two_thirds: np.array( + [ + [27, 0, 0, 0], + [9, 18, 0, 0], + [3, 12, 12, 0], + [1, 6, 12, 8], + [1, 6, 12, 8], + [0, 3, 12, 12], + [0, 0, 9, 18], + [0, 0, 0, 27], + ] + ) + / 27, + 1: np.array( + [ + [1, 0, 0, 0], + [0, 1, 0, 0], + [0, 0, 1, 0], + [0, 0, 0, 1], + [0, 0, 0, 1], + [0, 0, 0, 1], + [0, 0, 0, 1], + [0, 0, 0, 1], + ] + ), + }, + # Test case with a quartic Bézier + # to check if the fallback algorithms work + 4: { + 0: np.array( + [ + [1, 0, 0, 0, 0], + [1, 0, 0, 0, 0], + [1, 0, 0, 0, 0], + [1, 0, 0, 0, 0], + [1, 0, 0, 0, 0], + [1, 0, 0, 0, 0], + [0, 1, 0, 0, 0], + [0, 0, 1, 0, 0], + [0, 0, 0, 1, 0], + [0, 0, 0, 0, 1], + ] + ), + one_third: np.array( + [ + [81, 0, 0, 0, 0], + [54, 27, 0, 0, 0], + [36, 36, 9, 0, 0], + [24, 36, 18, 3, 0], + [16, 32, 24, 8, 1], + [16, 32, 24, 8, 1], + [0, 24, 36, 18, 3], + [0, 0, 36, 36, 9], + [0, 0, 0, 54, 27], + [0, 0, 0, 0, 81], + ] + ) + / 81, + two_thirds: np.array( + [ + [81, 0, 0, 0, 0], + [27, 54, 0, 0, 0], + [9, 36, 36, 0, 0], + [3, 18, 36, 24, 0], + [1, 8, 24, 32, 16], + [1, 8, 24, 32, 16], + [0, 3, 18, 36, 24], + [0, 0, 9, 36, 36], + [0, 0, 0, 27, 54], + [0, 0, 0, 0, 81], + ] + ) + / 81, + 1: np.array( + [ + [1, 0, 0, 0, 0], + [0, 1, 0, 0, 0], + [0, 0, 1, 0, 0], + [0, 0, 0, 1, 0], + [0, 0, 0, 0, 1], + [0, 0, 0, 0, 1], + [0, 0, 0, 0, 1], + [0, 0, 0, 0, 1], + [0, 0, 0, 0, 1], + [0, 0, 0, 0, 1], + ] + ), + }, +} diff --git a/tests/module/utils/_subdivision_matrices.py b/tests/module/utils/_subdivision_matrices.py new file mode 100644 index 0000000000..c020937dea --- /dev/null +++ b/tests/module/utils/_subdivision_matrices.py @@ -0,0 +1,167 @@ +import numpy as np + +# Expected values for matrices in subdivide_bezier and others +# Note that in bezier.py this is a list of dicts, +# not a dict of dicts! +SUBDIVISION_MATRICES = { + # For 0-degree Béziers + 0: { + 2: np.array([[1], [1]]), + 3: np.array([[1], [1], [1]]), + 4: np.array([[1], [1], [1], [1]]), + }, + # For linear Béziers + 1: { + 2: np.array( + [ + [2, 0], + [1, 1], + [1, 1], + [0, 2], + ] + ) + / 2, + 3: np.array( + [ + [3, 0], + [2, 1], + [2, 1], + [1, 2], + [1, 2], + [0, 3], + ] + ) + / 3, + 4: np.array( + [ + [4, 0], + [3, 1], + [3, 1], + [2, 2], + [2, 2], + [1, 3], + [1, 3], + [0, 4], + ] + ) + / 4, + }, + # For quadratic Béziers + 2: { + 2: np.array( + [ + [4, 0, 0], + [2, 2, 0], + [1, 2, 1], + [1, 2, 1], + [0, 2, 2], + [0, 0, 4], + ] + ) + / 4, + 3: np.array( + [ + [9, 0, 0], + [6, 3, 0], + [4, 4, 1], + [4, 4, 1], + [2, 5, 2], + [1, 4, 4], + [1, 4, 4], + [0, 3, 6], + [0, 0, 9], + ] + ) + / 9, + 4: np.array( + [ + [16, 0, 0], + [12, 4, 0], + [9, 6, 1], + [9, 6, 1], + [6, 8, 2], + [4, 8, 4], + [4, 8, 4], + [2, 8, 6], + [1, 6, 9], + [1, 6, 9], + [0, 4, 12], + [0, 0, 16], + ] + ) + / 16, + }, + # For cubic Béziers + 3: { + 2: np.array( + [ + [8, 0, 0, 0], + [4, 4, 0, 0], + [2, 4, 2, 0], + [1, 3, 3, 1], + [1, 3, 3, 1], + [0, 2, 4, 2], + [0, 0, 4, 4], + [0, 0, 0, 8], + ] + ) + / 8, + 3: np.array( + [ + [27, 0, 0, 0], + [18, 9, 0, 0], + [12, 12, 3, 0], + [8, 12, 6, 1], + [8, 12, 6, 1], + [4, 12, 9, 2], + [2, 9, 12, 4], + [1, 6, 12, 8], + [1, 6, 12, 8], + [0, 3, 12, 12], + [0, 0, 9, 18], + [0, 0, 0, 27], + ] + ) + / 27, + 4: np.array( + [ + [64, 0, 0, 0], + [48, 16, 0, 0], + [36, 24, 4, 0], + [27, 27, 9, 1], + [27, 27, 9, 1], + [18, 30, 14, 2], + [12, 28, 20, 4], + [8, 24, 24, 8], + [8, 24, 24, 8], + [4, 20, 28, 12], + [2, 14, 30, 18], + [1, 9, 27, 27], + [1, 9, 27, 27], + [0, 4, 24, 36], + [0, 0, 16, 48], + [0, 0, 0, 64], + ] + ) + / 64, + }, + # Test case with a quartic Bézier + # to check if the fallback algorithms work + 4: { + 2: np.array( + [ + [16, 0, 0, 0, 0], + [8, 8, 0, 0, 0], + [4, 8, 4, 0, 0], + [2, 6, 6, 2, 0], + [1, 4, 6, 4, 1], + [1, 4, 6, 4, 1], + [0, 2, 6, 6, 2], + [0, 0, 4, 8, 4], + [0, 0, 0, 8, 8], + [0, 0, 0, 0, 16], + ] + ) + / 16, + }, +} diff --git a/tests/module/utils/test_bezier.py b/tests/module/utils/test_bezier.py new file mode 100644 index 0000000000..7e1351c961 --- /dev/null +++ b/tests/module/utils/test_bezier.py @@ -0,0 +1,97 @@ +from __future__ import annotations + +import numpy as np +import numpy.testing as nt +from _split_matrices import SPLIT_MATRICES +from _subdivision_matrices import SUBDIVISION_MATRICES + +from manim.utils.bezier import ( + _get_subdivision_matrix, + partial_bezier_points, + split_bezier, + subdivide_bezier, +) + +QUARTIC_BEZIER = np.array( + [ + [-1, -1, 0], + [-1, 0, 0], + [0, 1, 0], + [1, 0, 0], + [1, -1, 0], + ], + dtype=float, +) + + +def test_partial_bezier_points() -> None: + """Test that :func:`partial_bezierpoints`, both in the + portion-matrix-building algorithm (degrees up to 3) and the + fallback algorithm (degree 4), works correctly. + """ + for degree, degree_dict in SUBDIVISION_MATRICES.items(): + n_points = degree + 1 + points = QUARTIC_BEZIER[:n_points] + for n_divisions, subdivision_matrix in degree_dict.items(): + for i in range(n_divisions): + a = i / n_divisions + b = (i + 1) / n_divisions + portion_matrix = subdivision_matrix[n_points * i : n_points * (i + 1)] + nt.assert_allclose( + partial_bezier_points(points, a, b), + portion_matrix @ points, + atol=1e-15, # Needed because of floating-point errors + ) + + +def test_split_bezier() -> None: + """Test that :func:`split_bezier`, both in the + split-matrix-building algorithm (degrees up to 3) and the + fallback algorithm (degree 4), works correctly. + """ + for degree, degree_dict in SPLIT_MATRICES.items(): + n_points = degree + 1 + points = QUARTIC_BEZIER[:n_points] + for t, split_matrix in degree_dict.items(): + nt.assert_allclose( + split_bezier(points, t), split_matrix @ points, atol=1e-15 + ) + + for degree, degree_dict in SUBDIVISION_MATRICES.items(): + n_points = degree + 1 + points = QUARTIC_BEZIER[:n_points] + # Split in half + split_matrix = degree_dict[2] + nt.assert_allclose( + split_bezier(points, 0.5), + split_matrix @ points, + ) + + +def test_get_subdivision_matrix() -> None: + """Test that the memos in .:meth:`_get_subdivision_matrix` + are being correctly generated. + """ + # Only for degrees up to 3! + for degree in range(4): + degree_dict = SUBDIVISION_MATRICES[degree] + for n_divisions, subdivision_matrix in degree_dict.items(): + nt.assert_allclose( + _get_subdivision_matrix(degree + 1, n_divisions), + subdivision_matrix, + ) + + +def test_subdivide_bezier() -> None: + """Test that :func:`subdivide_bezier`, both in the memoized cases + (degrees up to 3) and the fallback algorithm (degree 4), works + correctly. + """ + for degree, degree_dict in SUBDIVISION_MATRICES.items(): + n_points = degree + 1 + points = QUARTIC_BEZIER[:n_points] + for n_divisions, subdivision_matrix in degree_dict.items(): + nt.assert_allclose( + subdivide_bezier(points, n_divisions), + subdivision_matrix @ points, + )