diff --git a/manim/mobject/opengl/opengl_vectorized_mobject.py b/manim/mobject/opengl/opengl_vectorized_mobject.py index d668e3a2b1..0b94427cce 100644 --- a/manim/mobject/opengl/opengl_vectorized_mobject.py +++ b/manim/mobject/opengl/opengl_vectorized_mobject.py @@ -15,13 +15,14 @@ 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, + subdivide_bezier, ) from manim.utils.color import BLACK, WHITE, ManimColor, ParsableManimColor from manim.utils.config_ops import _Data @@ -547,21 +548,23 @@ def is_closed(self): def subdivide_sharp_curves(self, angle_threshold=30 * DEGREES, recurse=True): vmobs = [vm for vm in self.get_family(recurse) if vm.has_points()] for vmob in vmobs: - new_points = [] - for tup in vmob.get_bezier_tuples(): - angle = angle_between_vectors(tup[1] - tup[0], tup[2] - tup[1]) - if angle > angle_threshold: - n = int(np.ceil(angle / angle_threshold)) - alphas = np.linspace(0, 1, n + 1) - new_points.extend( - [ - partial_quadratic_bezier_points(tup, a1, a2) - for a1, a2 in zip(alphas, alphas[1:]) - ], - ) - else: - new_points.append(tup) - vmob.set_points(np.vstack(new_points)) + beziers = vmob.get_bezier_tuples() + angles = np.empty(beziers.shape[0]) + for i, bez in enumerate(beziers): + angles[i] = angle_between_vectors(bez[1] - bez[0], bez[2] - bez[1]) + + num_parts_per_bezier = np.ceil(angles / angle_threshold).astype(int) + acc_parts = np.add.accumulate(num_parts_per_bezier) + end_indices = self.n_points_per_curve * acc_parts + num_final_points = end_indices[-1] + new_points = np.empty((num_final_points, vmob.points.shape[1])) + + start = 0 + for bez, end, n in zip(beziers, end_indices, num_parts_per_bezier): + new_points[start:end] = subdivide_bezier(bez, n) + start = end + + vmob.set_points(new_points) return self def add_points_as_corners(self, points): @@ -1255,53 +1258,30 @@ def insert_n_curves(self, n: int, recurse=True) -> OpenGLVMobject: return self def insert_n_curves_to_point_list(self, n: int, points: np.ndarray) -> np.ndarray: - """Given an array of k points defining a bezier curves - (anchors and handles), returns points defining exactly - k + n bezier curves. + """Given an array of ``points`` defining :math:`k` Bézier curves (anchors and handles), + returns an array of points defining exactly :math:`k + n` Bézier curves. Parameters ---------- n - Number of desired curves. + Number of desired curves to add. points Starting points. Returns ------- - np.ndarray + :class:`np.ndarray` Points generated. """ - nppc = self.n_points_per_curve if len(points) == 1: + nppc = self.n_points_per_curve 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) + new_number_of_curves = bezier_tuples.shape[0] + 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 +1334,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 +1342,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..486d34424b 100644 --- a/manim/mobject/types/vectorized_mobject.py +++ b/manim/mobject/types/vectorized_mobject.py @@ -30,7 +30,8 @@ ) from manim.utils.bezier import ( bezier, - get_smooth_handle_points, + bezier_remap, + get_smooth_cubic_bezier_handle_points, integer_interpolate, interpolate, partial_bezier_points, @@ -1011,7 +1012,7 @@ def change_anchor_mode(self, mode: Literal["jagged", "smooth"]) -> Self: # The append is needed as the last element is not reached when slicing with numpy. anchors = np.append(subpath[::nppcc], subpath[-1:], 0) if mode == "smooth": - h1, h2 = get_smooth_handle_points(anchors) + h1, h2 = get_smooth_cubic_bezier_handle_points(anchors) else: # mode == "jagged" # The following will make the handles aligned with the anchors, thus making the bezier curve a segment a1 = anchors[:-1] @@ -1675,58 +1676,30 @@ def insert_n_curves(self, n: int) -> Self: def insert_n_curves_to_point_list( self, n: int, points: Point3D_Array - ) -> npt.NDArray[BezierPoints]: - """Given an array of k points defining a bezier curves (anchors and handles), returns points defining exactly k + n bezier curves. + ) -> Point3D_Array: + """Given an array of ``points`` defining :math:`k` Bézier curves (anchors and handles), + returns an array of points defining exactly :math:`k + n` Bézier curves. Parameters ---------- n - Number of desired curves. + Number of desired curves to add. points Starting points. Returns ------- + :class:`Point3D_Array` Points generated. """ - 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) + new_number_of_curves = bezier_tuples.shape[0] + 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 3f2bac86bc..b9e4e08025 100644 --- a/manim/utils/bezier.py +++ b/manim/utils/bezier.py @@ -4,27 +4,24 @@ from manim.typing import ( BezierPoints, + BezierPoints_Array, ColVector, - MatrixMN, Point3D, Point3D_Array, - PointDType, - QuadraticBezierPoints, - QuadraticBezierPoints_Array, ) __all__ = [ "bezier", "partial_bezier_points", - "partial_quadratic_bezier_points", + "split_bezier", + "subdivide_bezier", + "bezier_remap", "interpolate", "integer_interpolate", "mid", "inverse_interpolate", "match_interpolate", - "get_smooth_handle_points", "get_smooth_cubic_bezier_handle_points", - "diag_to_matrix", "is_closed", "proportions_along_bezier_curve_for_point", "point_lies_on_bezier", @@ -37,61 +34,248 @@ import numpy as np import numpy.typing as npt -from scipy import linalg from ..utils.simple_functions import choose from ..utils.space_ops import cross2d, find_intersection +@overload +def bezier( + points: BezierPoints, +) -> Callable[[float | ColVector], Point3D | Point3D_Array]: ... + + +@overload def bezier( - points: Sequence[Point3D] | Point3D_Array, -) -> Callable[[float], Point3D]: - """Classic implementation of a bezier curve. + points: Sequence[Point3D_Array], +) -> Callable[[float | ColVector], Point3D_Array]: ... + + +def bezier(points): + """Classic implementation of a Bézier curve. Parameters ---------- points - points defining the desired bezier curve. + :math:`(d+1, 3)`-shaped array of :math:`d+1` control points defining a single Bézier + curve of degree :math:`d`. Alternatively, for vectorization purposes, ``points`` can + also be a :math:`(d+1, M, 3)`-shaped sequence of :math:`d+1` arrays of :math:`M` + control points each, which define `M` Bézier curves instead. Returns ------- - function describing the bezier curve. - You can pass a t value between 0 and 1 to get the corresponding point on the curve. + bezier_func : :class:`typing.Callable` [[:class:`float` | :class:`~.ColVector`], :class:`~.Point3D` | :class:`~.Point3D_Array`] + Function describing the Bézier curve. The behaviour of this function depends on + the shape of ``points``: + + * If ``points`` was a :math:`(d+1, 3)` array representing a single Bézier curve, + then ``bezier_func`` can receive either: + + * a :class:`float` ``t``, in which case it returns a + single :math:`(1, 3)`-shaped :class:`~.Point3D` representing the evaluation + of the Bézier at ``t``, or + + * an :math:`(n, 1)`-shaped :class:`~.ColVector` + containing :math:`n` values to evaluate the Bézier curve at, returning instead + an :math:`(n, 3)`-shaped :class:`~.Point3D_Array` containing the points + resulting from evaluating the Bézier at each of the :math:`n` values. + + .. warning:: + If passing a vector of :math:`t`-values to ``bezier_func``, it **must** + be a column vector/matrix of shape :math:`(n, 1)`. Passing an 1D array of + shape :math:`(n,)` is not supported and **will result in undefined behaviour**. + + * If ``points`` was a :math:`(d+1, M, 3)` array describing :math:`M` Bézier curves, + then ``bezier_func`` can receive either: + + * a :class:`float` ``t``, in which case it returns an + :math:`(M, 3)`-shaped :class:`~.Point3D_Array` representing the evaluation + of the :math:`M` Bézier curves at the same value ``t``, or + + * an :math:`(M, 1)`-shaped + :class:`~.ColVector` containing :math:`M` values, such that the :math:`i`-th + Bézier curve defined by ``points`` is evaluated at the corresponding :math:`i`-th + value in ``t``, returning again an :math:`(M, 3)`-shaped :class:`~.Point3D_Array` + containing those :math:`M` evaluations. + + .. warning:: + Unlike the previous case, if you pass a :class:`~.ColVector` to ``bezier_func``, + it **must** contain exactly :math:`M` values, each value for each of the :math:`M` + Bézier curves defined by ``points``. Any array of shape other than :math:`(M, 1)` + **will result in undefined behaviour**. """ - n = len(points) - 1 - # Cubic Bezier curve - if n == 3: - return lambda t: np.asarray( - (1 - t) ** 3 * points[0] - + 3 * t * (1 - t) ** 2 * points[1] - + 3 * (1 - t) * t**2 * points[2] - + t**3 * points[3], - dtype=PointDType, - ) - # Quadratic Bezier curve - if n == 2: - return lambda t: np.asarray( - (1 - t) ** 2 * points[0] + 2 * t * (1 - t) * points[1] + t**2 * points[2], - dtype=PointDType, - ) + P = np.asarray(points) + d = P.shape[0] - 1 - return lambda t: np.asarray( - np.asarray( - [ - (((1 - t) ** (n - k)) * (t**k) * choose(n, k) * point) - for k, point in enumerate(points) - ], - dtype=PointDType, - ).sum(axis=0) - ) + if d == 0: + def zero_bezier(t): + return np.ones_like(t) * P[0] -# !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 BezierPoints: + r"""Given an array of ``points`` which define a Bézier curve, and two numbers :math:`a, b` + such that :math:`0 \le a < b \le 1`, return an array of the same size, which describes the + portion of the original Bézier curve on the interval :math:`[a, b]`. + + .. seealso:: + See :func:`split_bezier` for an explanation on how to split Bézier curves. + + .. note:: + To find the portion of a Bézier curve with :math:`t` between :math:`a` and :math:`b`: + + 1. Split the curve at :math:`t = a` and extract its 2nd subcurve. + 2. We cannot evaluate the new subcurve at :math:`t = b` because its range of values for :math:`t` is different. + To find the correct value, we need to transform the interval :math:`[a, 1]` into :math:`[0, 1]` + by first subtracting :math:`a` to get :math:`[0, 1-a]` and then dividing by :math:`1-a`. Thus, our new + value must be :math:`t = \frac{b - a}{1 - a}`. Define :math:`u = \frac{b - a}{1 - a}`. + 3. Split the subcurve at :math:`t = u` and extract its 1st subcurve. + + The final portion is a linear combination of points, and thus the process can be + summarized as a linear transformation by some matrix in terms of :math:`a` and :math:`b`. + This matrix is given explicitly for Bézier curves up to degree 3, which are often used in Manim. + For higher degrees, the algorithm described previously is used. + + For the case of a quadratic Bézier curve: + + * Step 1: + + .. math:: + H'_1 + = + \begin{pmatrix} + (1-a)^2 & 2(1-a)a & a^2 \\ + 0 & (1-a) & a \\ + 0 & 0 & 1 + \end{pmatrix} + \begin{pmatrix} + p_0 \\ + p_1 \\ + p_2 + \end{pmatrix} + + * Step 2: + + .. math:: + H''_0 + &= + \begin{pmatrix} + 1 & 0 & 0 \\ + (1-u) & u & 0\\ + (1-u)^2 & 2(1-u)u & u^2 + \end{pmatrix} + H'_1 + \\ + & + \\ + &= + \begin{pmatrix} + 1 & 0 & 0 \\ + (1-u) & u & 0\\ + (1-u)^2 & 2(1-u)u & u^2 + \end{pmatrix} + \begin{pmatrix} + (1-a)^2 & 2(1-a)a & a^2 \\ + 0 & (1-a) & a \\ + 0 & 0 & 1 + \end{pmatrix} + \begin{pmatrix} + p_0 \\ + p_1 \\ + p_2 + \end{pmatrix} + \\ + & + \\ + &= + \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} + \begin{pmatrix} + p_0 \\ + p_1 \\ + p_2 + \end{pmatrix} + + from where one can define a :math:`(3, 3)` matrix :math:`P_2` which, when applied over + the array of ``points``, will return the desired partial quadratic Bézier curve: + + .. 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} + + Similarly, for the cubic Bézier curve case, one can define the following + :math:`(4, 4)` matrix :math:`P_3`: + + .. math:: + P_3 + = + \begin{pmatrix} + (1-a)^3 & 3(1-a)^2a & 3(1-a)a^2 & a^3 \\ + (1-a)^2(1-b) & 2(1-a)a(1-b) + (1-a)^2b & a^2(1-b) + 2(1-a)ab & a^2b \\ + (1-a)(1-b)^2 & a(1-b)^2 + 2(1-a)(1-b)b & 2a(1-b)b + (1-a)b^2 & ab^2 \\ + (1-b)^3 & 3(1-b)^2b & 3(1-b)b^2 & b^3 + \end{pmatrix} Parameters ---------- @@ -104,158 +288,710 @@ def partial_bezier_points(points: BezierPoints, a: float, b: float) -> 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) - - 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, - ) + arr = np.array(points) + arr[:] = arr[-1] + return arr + if b == 0: + arr = np.array(points) + arr[:] = arr[0] + return arr + points = np.asarray(points) + degree = points.shape[0] - 1 -# 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]]) + 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 - 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( + [ + [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 - # 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 == 2: + ma, mb = 1 - a, 1 - b + 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 -def split_quadratic_bezier(points: QuadraticBezierPoints, t: float) -> BezierPoints: - """Split a quadratic Bézier curve at argument ``t`` into two quadratic curves. + 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:: + + .. 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 -def subdivide_quadratic_bezier(points: QuadraticBezierPoints, n: int) -> BezierPoints: - """Subdivide a quadratic Bézier curve into ``n`` subcurves which have the same shape. + 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 _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 @@ -265,14 +1001,48 @@ def quadratic_bezier_remap( def interpolate(start: float, end: float, alpha: float) -> float: ... +@overload +def interpolate(start: float, end: float, alpha: ColVector) -> ColVector: ... + + @overload def interpolate(start: Point3D, end: Point3D, alpha: float) -> Point3D: ... -def interpolate( - start: int | float | Point3D, end: int | float | Point3D, alpha: float | Point3D -) -> float | Point3D: - return (1 - alpha) * start + alpha * end +@overload +def interpolate(start: Point3D, end: Point3D, alpha: ColVector) -> Point3D_Array: ... + + +def interpolate(start, end, alpha): + """Linearly interpolates between two values ``start`` and ``end``. + + Parameters + ---------- + start + The start of the range. + end + The end of the range. + alpha + A float between 0 and 1, or an :math:`(n, 1)` column vector containing + :math:`n` floats between 0 and 1 to interpolate in a vectorized fashion. + + Returns + ------- + :class:`float` | :class:`~.ColVector` | :class:`~.Point3D` | :class:`~.Point3D_Array` + The result of the linear interpolation. + + * If ``start`` and ``end`` are of type :class:`float`, and: + + * ``alpha`` is also a :class:`float`, the return is simply another :class:`float`. + * ``alpha`` is a :class:`~.ColVector`, the return is another :class:`~.ColVector`. + + * If ``start`` and ``end`` are of type :class:`~.Point3D`, and: + + * ``alpha`` is a :class:`float`, the return is another :class:`~.Point3D`. + * ``alpha`` is a :class:`~.ColVector`, the return is a :class:`~.Point3D_Array`. + + """ + return start + alpha * (end - start) def integer_interpolate( @@ -280,21 +1050,20 @@ def integer_interpolate( end: float, alpha: float, ) -> tuple[int, float]: - """ - This is a variant of interpolate that returns an integer and the residual + """Variant of :func:`interpolate` that returns an integer and the residual. Parameters ---------- start - The start of the range + The start of the range. end - The end of the range + The end of the range. alpha - a float between 0 and 1. + A float between 0 and 1. Returns ------- - tuple[int, float] + :class:`tuple` [:class:`int`, :class:`float`] This returns an integer between start and end (inclusive) representing appropriate interpolation between them, along with a "residue" representing a new proportion between the @@ -327,19 +1096,20 @@ def mid(start: float, end: float) -> float: ... def mid(start: Point3D, end: Point3D) -> Point3D: ... -def mid(start: float | Point3D, end: float | Point3D) -> float | Point3D: +def mid(start, end): """Returns the midpoint between two values. Parameters ---------- start - The first value + The first value. end - The second value + The second value. Returns ------- - The midpoint between the two values + :class:`float` | :class:`~.Point3D` + The midpoint between the two values. """ return (start + end) / 2.0 @@ -356,9 +1126,7 @@ def inverse_interpolate(start: float, end: float, value: Point3D) -> Point3D: .. def inverse_interpolate(start: Point3D, end: Point3D, value: Point3D) -> Point3D: ... -def inverse_interpolate( - start: float | Point3D, end: float | Point3D, value: float | Point3D -) -> float | Point3D: +def inverse_interpolate(start, end, value): """Perform inverse interpolation to determine the alpha values that would produce the specified ``value`` given the ``start`` and ``end`` values or points. @@ -375,6 +1143,7 @@ def inverse_interpolate( Returns ------- + :class:`float` | :class:`~.Point3D` The alpha values producing the given input when interpolating between ``start`` and ``end``. @@ -415,13 +1184,7 @@ def match_interpolate( ) -> Point3D: ... -def match_interpolate( - new_start: float, - new_end: float, - old_start: float, - old_end: float, - old_value: float | Point3D, -) -> float | Point3D: +def match_interpolate(new_start, new_end, old_start, old_end, old_value): """Interpolate a value from an old range to a new range. Parameters @@ -441,6 +1204,7 @@ def match_interpolate( Returns ------- + :class:`float` | :class:`~.Point3D` The interpolated value within the new range. Examples @@ -456,163 +1220,621 @@ def match_interpolate( ) +# Figuring out which Bézier curves most smoothly connect a sequence of points def get_smooth_cubic_bezier_handle_points( - points: Point3D_Array, -) -> tuple[BezierPoints, BezierPoints]: - points = np.asarray(points) - num_handles = len(points) - 1 - dim = points.shape[1] - if num_handles < 1: - return np.zeros((0, dim)), np.zeros((0, dim)) - # Must solve 2*num_handles equations to get the handles. - # l and u are the number of lower an upper diagonal rows - # in the matrix to solve. - l, u = 2, 1 - # diag is a representation of the matrix in diagonal form - # See https://www.particleincell.com/2012/bezier-splines/ - # for how to arrive at these equations - diag: MatrixMN = np.zeros((l + u + 1, 2 * num_handles)) - diag[0, 1::2] = -1 - diag[0, 2::2] = 1 - diag[1, 0::2] = 2 - diag[1, 1::2] = 1 - diag[2, 1:-2:2] = -2 - diag[3, 0:-3:2] = 1 - # last - diag[2, -2] = -1 - diag[1, -1] = 2 - # This is the b as in Ax = b, where we are solving for x, - # and A is represented using diag. However, think of entries - # to x and b as being points in space, not numbers - b: Point3D_Array = np.zeros((2 * num_handles, dim)) - b[1::2] = 2 * points[1:] - b[0] = points[0] - b[-1] = points[-1] - - def solve_func(b: ColVector) -> ColVector | MatrixMN: - return linalg.solve_banded((l, u), diag, b) # type: ignore - - use_closed_solve_function = is_closed(points) - if use_closed_solve_function: - # Get equations to relate first and last points - matrix = diag_to_matrix((l, u), diag) - # last row handles second derivative - matrix[-1, [0, 1, -2, -1]] = [2, -1, 1, -2] - # first row handles first derivative - matrix[0, :] = np.zeros(matrix.shape[1]) - matrix[0, [0, -1]] = [1, 1] - b[0] = 2 * points[0] - b[-1] = np.zeros(dim) - - def closed_curve_solve_func(b: ColVector) -> ColVector | MatrixMN: - return linalg.solve(matrix, b) # type: ignore - - handle_pairs = np.zeros((2 * num_handles, dim)) - for i in range(dim): - if use_closed_solve_function: - handle_pairs[:, i] = closed_curve_solve_func(b[:, i]) - else: - handle_pairs[:, i] = solve_func(b[:, i]) - return handle_pairs[0::2], handle_pairs[1::2] - - -def get_smooth_handle_points( - points: BezierPoints, -) -> tuple[BezierPoints, BezierPoints]: - """Given some anchors (points), compute handles so the resulting bezier curve is smooth. + anchors: Point3D_Array, +) -> tuple[Point3D_Array, Point3D_Array]: + """Given an array of anchors for a cubic spline (array of connected cubic + Bézier curves), compute the 1st and 2nd handle for every curve, so that + the resulting spline is smooth. Parameters ---------- - points - Anchors. + anchors + Anchors of a cubic spline. Returns ------- - typing.Tuple[np.ndarray, np.ndarray] - Computed handles. + :class:`tuple` [:class:`~.Point3D_Array`, :class:`~.Point3D_Array`] + A tuple of two arrays: one containing the 1st handle for every curve in + the cubic spline, and the other containing the 2nd handles. """ - # NOTE points here are anchors. - points = np.asarray(points) - num_handles = len(points) - 1 - dim = points.shape[1] - if num_handles < 1: + anchors = np.asarray(anchors) + n_anchors = anchors.shape[0] + + # If there's a single anchor, there's no Bézier curve. + # Return empty arrays. + if n_anchors == 1: + dim = anchors.shape[1] return np.zeros((0, dim)), np.zeros((0, dim)) - # Must solve 2*num_handles equations to get the handles. - # l and u are the number of lower an upper diagonal rows - # in the matrix to solve. - l, u = 2, 1 - # diag is a representation of the matrix in diagonal form - # See https://www.particleincell.com/2012/bezier-splines/ - # for how to arrive at these equations - diag: MatrixMN = np.zeros((l + u + 1, 2 * num_handles)) - diag[0, 1::2] = -1 - diag[0, 2::2] = 1 - diag[1, 0::2] = 2 - diag[1, 1::2] = 1 - diag[2, 1:-2:2] = -2 - diag[3, 0:-3:2] = 1 - # last - diag[2, -2] = -1 - diag[1, -1] = 2 - # This is the b as in Ax = b, where we are solving for x, - # and A is represented using diag. However, think of entries - # to x and b as being points in space, not numbers - b = np.zeros((2 * num_handles, dim)) - b[1::2] = 2 * points[1:] - b[0] = points[0] - b[-1] = points[-1] - - def solve_func(b: ColVector) -> ColVector | MatrixMN: - return linalg.solve_banded((l, u), diag, b) # type: ignore - - use_closed_solve_function = is_closed(points) - if use_closed_solve_function: - # Get equations to relate first and last points - matrix = diag_to_matrix((l, u), diag) - # last row handles second derivative - matrix[-1, [0, 1, -2, -1]] = [2, -1, 1, -2] - # first row handles first derivative - matrix[0, :] = np.zeros(matrix.shape[1]) - matrix[0, [0, -1]] = [1, 1] - b[0] = 2 * points[0] - b[-1] = np.zeros(dim) - - def closed_curve_solve_func(b: ColVector) -> ColVector | MatrixMN: - return linalg.solve(matrix, b) # type: ignore - - handle_pairs = np.zeros((2 * num_handles, dim)) - for i in range(dim): - if use_closed_solve_function: - handle_pairs[:, i] = closed_curve_solve_func(b[:, i]) - else: - handle_pairs[:, i] = solve_func(b[:, i]) - return handle_pairs[0::2], handle_pairs[1::2] + # If there are only two anchors (thus only one pair of handles), + # they can only be an interpolation of these two anchors with alphas + # 1/3 and 2/3, which will draw a straight line between the anchors. + if n_anchors == 2: + return interpolate(anchors[0], anchors[1], np.array([[1 / 3], [2 / 3]])) + + # Handle different cases depending on whether the points form a closed + # curve or not + curve_is_closed = is_closed(anchors) + if curve_is_closed: + return get_smooth_closed_cubic_bezier_handle_points(anchors) + else: + return get_smooth_open_cubic_bezier_handle_points(anchors) + + +CP_CLOSED_MEMO = np.array([1 / 3]) +UP_CLOSED_MEMO = np.array([1 / 3]) + + +def get_smooth_closed_cubic_bezier_handle_points( + anchors: Point3D_Array, +) -> tuple[Point3D_Array, Point3D_Array]: + r"""Special case of :func:`get_smooth_cubic_bezier_handle_points`, + when the ``anchors`` form a closed loop. + + .. note:: + A system of equations must be solved to get the first handles of + every Bézier curve (referred to as :math:`H_1`). + Then :math:`H_2` (the second handles) can be obtained separately. + + .. seealso:: + The equations were obtained from: + + * `Conditions on control points for continuous curvature. (2016). Jaco Stuifbergen. `_ + + In general, if there are :math:`N+1` anchors, there will be :math:`N` Bézier curves + and thus :math:`N` pairs of handles to find. We must solve the following + system of equations for the 1st handles (example for :math:`N = 5`): + + .. math:: + \begin{pmatrix} + 4 & 1 & 0 & 0 & 1 \\ + 1 & 4 & 1 & 0 & 0 \\ + 0 & 1 & 4 & 1 & 0 \\ + 0 & 0 & 1 & 4 & 1 \\ + 1 & 0 & 0 & 1 & 4 + \end{pmatrix} + \begin{pmatrix} + H_{1,0} \\ + H_{1,1} \\ + H_{1,2} \\ + H_{1,3} \\ + H_{1,4} + \end{pmatrix} + = + \begin{pmatrix} + 4A_0 + 2A_1 \\ + 4A_1 + 2A_2 \\ + 4A_2 + 2A_3 \\ + 4A_3 + 2A_4 \\ + 4A_4 + 2A_5 + \end{pmatrix} + + which will be expressed as :math:`RH_1 = D`. + + :math:`R` is almost a tridiagonal matrix, so we could use Thomas' algorithm. + + .. seealso:: + `Tridiagonal matrix algorithm. Wikipedia. `_ + + However, :math:`R` has ones at the opposite corners. A solution to this is + the first decomposition proposed in the link below, with :math:`\alpha = 1`: + + .. seealso:: + `Tridiagonal matrix algorithm # Variants. Wikipedia. `_ + + .. math:: + R + = + \begin{pmatrix} + 4 & 1 & 0 & 0 & 1 \\ + 1 & 4 & 1 & 0 & 0 \\ + 0 & 1 & 4 & 1 & 0 \\ + 0 & 0 & 1 & 4 & 1 \\ + 1 & 0 & 0 & 1 & 4 + \end{pmatrix} + &= + \begin{pmatrix} + 3 & 1 & 0 & 0 & 0 \\ + 1 & 4 & 1 & 0 & 0 \\ + 0 & 1 & 4 & 1 & 0 \\ + 0 & 0 & 1 & 4 & 1 \\ + 0 & 0 & 0 & 1 & 3 + \end{pmatrix} + + + \begin{pmatrix} + 1 & 0 & 0 & 0 & 1 \\ + 0 & 0 & 0 & 0 & 0 \\ + 0 & 0 & 0 & 0 & 0 \\ + 0 & 0 & 0 & 0 & 0 \\ + 1 & 0 & 0 & 0 & 1 + \end{pmatrix} + \\ + & + \\ + &= + \begin{pmatrix} + 3 & 1 & 0 & 0 & 0 \\ + 1 & 4 & 1 & 0 & 0 \\ + 0 & 1 & 4 & 1 & 0 \\ + 0 & 0 & 1 & 4 & 1 \\ + 0 & 0 & 0 & 1 & 3 + \end{pmatrix} + + + \begin{pmatrix} + 1 \\ + 0 \\ + 0 \\ + 0 \\ + 0 + \end{pmatrix} + \begin{pmatrix} + 1 & 0 & 0 & 0 & 1 + \end{pmatrix} + \\ + & + \\ + &= + T + uv^t + + We decompose :math:`R = T + uv^t`, where :math:`T` is a tridiagonal matrix, and + :math:`u, v` are :math:`N`-D vectors such that :math:`u_0 = u_{N-1} = v_0 = v_{N-1} = 1`, + and :math:`u_i = v_i = 0, \forall i \in \{1, ..., N-2\}`. + + Thus: + + .. math:: + RH_1 &= D \\ + \Rightarrow (T + uv^t)H_1 &= D + + If we find a vector :math:`q` such that :math:`Tq = u`: + + .. math:: + \Rightarrow (T + Tqv^t)H_1 &= D \\ + \Rightarrow T(I + qv^t)H_1 &= D \\ + \Rightarrow H_1 &= (I + qv^t)^{-1} T^{-1} D + + According to Sherman-Morrison's formula: + + .. seealso:: + `Sherman-Morrison's formula. Wikipedia. `_ + + .. math:: + (I + qv^t)^{-1} = I - \frac{1}{1 + v^tq} qv^t + + If we find :math:`Y = T^{-1} D`, or in other words, if we solve for + :math:`Y` in :math:`TY = D`: + + .. math:: + H_1 &= (I + qv^t)^{-1} T^{-1} D \\ + &= (I + qv^t)^{-1} Y \\ + &= (I - \frac{1}{1 + v^tq} qv^t) Y \\ + &= Y - \frac{1}{1 + v^tq} qv^tY + + Therefore, we must solve for :math:`q` and :math:`Y` in :math:`Tq = u` and :math:`TY = D`. + As :math:`T` is now tridiagonal, we shall use Thomas' algorithm. + + Define: + + * :math:`a = [a_0, \ a_1, \ ..., \ a_{N-2}]` as :math:`T`'s lower diagonal of :math:`N-1` elements, + such that :math:`a_0 = a_1 = ... = a_{N-2} = 1`, so this diagonal is filled with ones; + * :math:`b = [b_0, \ b_1, \ ..., \ b_{N-2}, \ b_{N-1}]` as :math:`T`'s main diagonal of :math:`N` elements, + such that :math:`b_0 = b_{N-1} = 3`, and :math:`b_1 = b_2 = ... = b_{N-2} = 4`; + * :math:`c = [c_0, \ c_1, \ ..., \ c_{N-2}]` as :math:`T`'s upper diagonal of :math:`N-1` elements, + such that :math:`c_0 = c_1 = ... = c_{N-2} = 1`: this diagonal is also filled with ones. + + If, according to Thomas' algorithm, we define: + + .. math:: + c'_0 &= \frac{c_0}{b_0} & \\ + c'_i &= \frac{c_i}{b_i - a_{i-1} c'_{i-1}}, & \quad \forall i \in \{1, ..., N-2\} \\ + & & \\ + u'_0 &= \frac{u_0}{b_0} & \\ + u'_i &= \frac{u_i - a_{i-1} u'_{i-1}}{b_i - a_{i-1} c'_{i-1}}, & \quad \forall i \in \{1, ..., N-1\} \\ + & & \\ + D'_0 &= \frac{1}{b_0} D_0 & \\ + D'_i &= \frac{1}{b_i - a_{i-1} c'_{i-1}} (D_i - a_{i-1} D'_{i-1}), & \quad \forall i \in \{1, ..., N-1\} + + Then: + + .. math:: + c'_0 &= \frac{1}{3} & \\ + c'_i &= \frac{1}{4 - c'_{i-1}}, & \quad \forall i \in \{1, ..., N-2\} \\ + & & \\ + u'_0 &= \frac{1}{3} & \\ + u'_i &= \frac{-u'_{i-1}}{4 - c'_{i-1}} = -c'_i u'_{i-1}, & \quad \forall i \in \{1, ..., N-2\} \\ + u'_{N-1} &= \frac{1 - u'_{N-2}}{3 - c'_{N-2}} & \\ + & & \\ + D'_0 &= \frac{1}{3} (4A_0 + 2A_1) & \\ + D'_i &= \frac{1}{4 - c'_{i-1}} (4A_i + 2A_{i+1} - D'_{i-1}) & \\ + &= c_i (4A_i + 2A_{i+1} - D'_{i-1}), & \quad \forall i \in \{1, ..., N-2\} \\ + D'_{N-1} &= \frac{1}{3 - c'_{N-2}} (4A_{N-1} + 2A_N - D'_{N-2}) & + + Finally, we can do Backward Substitution to find :math:`q` and :math:`Y`: + + .. math:: + q_{N-1} &= u'_{N-1} & \\ + q_i &= u'_{i} - c'_i q_{i+1}, & \quad \forall i \in \{0, ..., N-2\} \\ + & & \\ + Y_{N-1} &= D'_{N-1} & \\ + Y_i &= D'_i - c'_i Y_{i+1}, & \quad \forall i \in \{0, ..., N-2\} + + With those values, we can finally calculate :math:`H_1 = Y - \frac{1}{1 + v^tq} qv^tY`. + Given that :math:`v_0 = v_{N-1} = 1`, and :math:`v_1 = v_2 = ... = v_{N-2} = 0`, its dot products + with :math:`q` and :math:`Y` are respectively :math:`v^tq = q_0 + q_{N-1}` and + :math:`v^tY = Y_0 + Y_{N-1}`. Thus: + + .. math:: + H_1 = Y - \frac{1}{1 + q_0 + q_{N-1}} q(Y_0 + Y_{N-1}) + + Once we have :math:`H_1`, we can get :math:`H_2` (the array of second handles) as follows: + + .. math:: + H_{2, i} &= 2A_{i+1} - H_{1, i+1}, & \quad \forall i \in \{0, ..., N-2\} \\ + H_{2, N-1} &= 2A_0 - H_{1, 0} & + + Because the matrix :math:`R` always follows the same pattern (and thus :math:`T, u, v` as well), + we can define a memo list for :math:`c'` and :math:`u'` to avoid recalculation. We cannot + memoize :math:`D` and :math:`Y`, however, because they are always different matrices. We + cannot make a memo for :math:`q` either, but we can calculate it faster because :math:`u'` + can be memoized. + + Parameters + ---------- + anchors + Anchors of a closed cubic spline. -def diag_to_matrix( - l_and_u: tuple[int, int], diag: npt.NDArray[Any] -) -> npt.NDArray[Any]: + Returns + ------- + :class:`tuple` [:class:`~.Point3D_Array`, :class:`~.Point3D_Array`] + A tuple of two arrays: one containing the 1st handle for every curve in + the closed cubic spline, and the other containing the 2nd handles. """ - Converts array whose rows represent diagonal - entries of a matrix into the matrix itself. - See scipy.linalg.solve_banded + global CP_CLOSED_MEMO + global UP_CLOSED_MEMO + + A = np.asarray(anchors) + N = A.shape[0] - 1 + dim = A.shape[1] + + # Calculate cp (c prime) and up (u prime) with help from + # CP_CLOSED_MEMO and UP_CLOSED_MEMO. + len_memo = CP_CLOSED_MEMO.size + if len_memo < N - 1: + cp = np.empty(N - 1) + up = np.empty(N - 1) + cp[:len_memo] = CP_CLOSED_MEMO + up[:len_memo] = UP_CLOSED_MEMO + # Forward Substitution 1 + # Calculate up (at the same time we calculate cp). + for i in range(len_memo, N - 1): + cp[i] = 1 / (4 - cp[i - 1]) + up[i] = -cp[i] * up[i - 1] + CP_CLOSED_MEMO = cp + UP_CLOSED_MEMO = up + else: + cp = CP_CLOSED_MEMO[: N - 1] + up = UP_CLOSED_MEMO[: N - 1] + + # The last element of u' is different + cp_last_division = 1 / (3 - cp[N - 2]) + up_last = cp_last_division * (1 - up[N - 2]) + + # Backward Substitution 1 + # Calculate q. + q = np.empty((N, dim)) + q[N - 1] = up_last + for i in range(N - 2, -1, -1): + q[i] = up[i] - cp[i] * q[i + 1] + + # Forward Substitution 2 + # Calculate Dp (D prime). + Dp = np.empty((N, dim)) + AUX = 4 * A[:N] + 2 * A[1:] # Vectorize the sum for efficiency. + Dp[0] = AUX[0] / 3 + for i in range(1, N - 1): + Dp[i] = cp[i] * (AUX[i] - Dp[i - 1]) + Dp[N - 1] = cp_last_division * (AUX[N - 1] - Dp[N - 2]) + + # Backward Substitution + # Calculate Y, which is defined as a view of Dp for efficiency + # and semantic convenience at the same time. + Y = Dp + # Y[N-1] = Dp[N-1] (redundant) + for i in range(N - 2, -1, -1): + Y[i] = Dp[i] - cp[i] * Y[i + 1] + + # Calculate H1. + H1 = Y - 1 / (1 + q[0] + q[N - 1]) * q * (Y[0] + Y[N - 1]) + + # Calculate H2. + H2 = np.empty((N, dim)) + H2[0 : N - 1] = 2 * A[1:N] - H1[1:N] + H2[N - 1] = 2 * A[N] - H1[0] + + return H1, H2 + + +CP_OPEN_MEMO = np.array([0.5]) + + +def get_smooth_open_cubic_bezier_handle_points( + anchors: Point3D_Array, +) -> tuple[Point3D_Array, Point3D_Array]: + r"""Special case of :func:`get_smooth_cubic_bezier_handle_points`, + when the ``anchors`` do not form a closed loop. + + .. note:: + A system of equations must be solved to get the first handles of + every Bèzier curve (referred to as :math:`H_1`). + Then :math:`H_2` (the second handles) can be obtained separately. + + .. seealso:: + The equations were obtained from: + + * `Smooth Bézier Spline Through Prescribed Points. (2012). Particle in Cell Consulting LLC. `_ + * `Conditions on control points for continuous curvature. (2016). Jaco Stuifbergen. `_ + + .. warning:: + The equations in the first webpage have some typos which were corrected in the comments. + + In general, if there are :math:`N+1` anchors, there will be :math:`N` Bézier curves + and thus :math:`N` pairs of handles to find. We must solve the following + system of equations for the 1st handles (example for :math:`N = 5`): + + .. math:: + \begin{pmatrix} + 2 & 1 & 0 & 0 & 0 \\ + 1 & 4 & 1 & 0 & 0 \\ + 0 & 1 & 4 & 1 & 0 \\ + 0 & 0 & 1 & 4 & 1 \\ + 0 & 0 & 0 & 2 & 7 + \end{pmatrix} + \begin{pmatrix} + H_{1,0} \\ + H_{1,1} \\ + H_{1,2} \\ + H_{1,3} \\ + H_{1,4} + \end{pmatrix} + = + \begin{pmatrix} + A_0 + 2A_1 \\ + 4A_1 + 2A_2 \\ + 4A_2 + 2A_3 \\ + 4A_3 + 2A_4 \\ + 8A_4 + A_5 + \end{pmatrix} + + which will be expressed as :math:`TH_1 = D`. + :math:`T` is a tridiagonal matrix, so the system can be solved in :math:`O(N)` + operations. Here we shall use Thomas' algorithm or the tridiagonal matrix + algorithm. + + .. seealso:: + `Tridiagonal matrix algorithm. Wikipedia. `_ + + Define: + + * :math:`a = [a_0, \ a_1, \ ..., \ a_{N-2}]` as :math:`T`'s lower diagonal of :math:`N-1` elements, + such that :math:`a_0 = a_1 = ... = a_{N-3} = 1`, and :math:`a_{N-2} = 2`; + * :math:`b = [b_0, \ b_1, \ ..., \ b_{N-2}, \ b_{N-1}]` as :math:`T`'s main diagonal of :math:`N` elements, + such that :math:`b_0 = 2`, :math:`b_1 = b_2 = ... = b_{N-2} = 4`, and :math:`b_{N-1} = 7`; + * :math:`c = [c_0, \ c_1, \ ..., \ c_{N-2}]` as :math:`T`'s upper diagonal of :math:`{N-1}` elements, + such that :math:`c_0 = c_1 = ... = c_{N-2} = 1`: this diagonal is filled with ones. + + If, according to Thomas' algorithm, we define: + + .. math:: + c'_0 &= \frac{c_0}{b_0} & \\ + c'_i &= \frac{c_i}{b_i - a_{i-1} c'_{i-1}}, & \quad \forall i \in \{1, ..., N-2\} \\ + & & \\ + D'_0 &= \frac{1}{b_0} D_0 & \\ + D'_i &= \frac{1}{b_i - a_{i-1} c'{i-1}} (D_i - a_{i-1} D'_{i-1}), & \quad \forall i \in \{1, ..., N-1\} + + Then: + + .. math:: + c'_0 &= 0.5 & \\ + c'_i &= \frac{1}{4 - c'_{i-1}}, & \quad \forall i \in \{1, ..., N-2\} \\ + & & \\ + D'_0 &= 0.5A_0 + A_1 & \\ + D'_i &= \frac{1}{4 - c'_{i-1}} (4A_i + 2A_{i+1} - D'_{i-1}) & \\ + &= c_i (4A_i + 2A_{i+1} - D'_{i-1}), & \quad \forall i \in \{1, ..., N-2\} \\ + D'_{N-1} &= \frac{1}{7 - 2c'_{N-2}} (8A_{N-1} + A_N - 2D'_{N-2}) & + + Finally, we can do Backward Substitution to find :math:`H_1`: + + .. math:: + H_{1, N-1} &= D'_{N-1} & \\ + H_{1, i} &= D'_i - c'_i H_{1, i+1}, & \quad \forall i \in \{0, ..., N-2\} + + Once we have :math:`H_1`, we can get :math:`H_2` (the array of second handles) as follows: + + .. math:: + H_{2, i} &= 2A_{i+1} - H_{1, i+1}, & \quad \forall i \in \{0, ..., N-2\} \\ + H_{2, N-1} &= 0.5A_N + 0.5H_{1, N-1} & + + As the matrix :math:`T` always follows the same pattern, we can define a memo list + for :math:`c'` to avoid recalculation. We cannot do the same for :math:`D`, however, + because it is always a different matrix. + + Parameters + ---------- + anchors + Anchors of an open cubic spline. + + Returns + ------- + :class:`tuple` [:class:`~.Point3D_Array`, :class:`~.Point3D_Array`] + A tuple of two arrays: one containing the 1st handle for every curve in + the open cubic spline, and the other containing the 2nd handles. """ - l, u = l_and_u - dim = diag.shape[1] - matrix = np.zeros((dim, dim)) - for i in range(l + u + 1): - np.fill_diagonal( - matrix[max(0, i - u) :, max(0, u - i) :], - diag[i, max(0, u - i) :], - ) - return matrix + global CP_OPEN_MEMO + + A = np.asarray(anchors) + N = A.shape[0] - 1 + dim = A.shape[1] + + # Calculate cp (c prime) with help from CP_OPEN_MEMO. + len_memo = CP_OPEN_MEMO.size + if len_memo < N - 1: + cp = np.empty(N - 1) + cp[:len_memo] = CP_OPEN_MEMO + for i in range(len_memo, N - 1): + cp[i] = 1 / (4 - cp[i - 1]) + CP_OPEN_MEMO = cp + else: + cp = CP_OPEN_MEMO[: N - 1] + + # Calculate Dp (D prime). + Dp = np.empty((N, dim)) + Dp[0] = 0.5 * A[0] + A[1] + AUX = 4 * A[1 : N - 1] + 2 * A[2:N] # Vectorize the sum for efficiency. + for i in range(1, N - 1): + Dp[i] = cp[i] * (AUX[i - 1] - Dp[i - 1]) + Dp[N - 1] = (1 / (7 - 2 * cp[N - 2])) * (8 * A[N - 1] + A[N] - 2 * Dp[N - 2]) + + # Backward Substitution. + # H1 (array of the first handles) is defined as a view of Dp for efficiency + # and semantic convenience at the same time. + H1 = Dp + # H1[N-1] = Dp[N-1] (redundant) + for i in range(N - 2, -1, -1): + H1[i] = Dp[i] - cp[i] * H1[i + 1] + + # Calculate H2. + H2 = np.empty((N, dim)) + H2[0 : N - 1] = 2 * A[1:N] - H1[1:N] + H2[N - 1] = 0.5 * (A[N] + H1[N - 1]) + + return H1, H2 + + +@overload +def get_quadratic_approximation_of_cubic( + a0: Point3D, + h0: Point3D, + h1: Point3D, + a1: Point3D, +) -> Point3D_Array: ... -# 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: Point3D_Array, + h0: Point3D_Array, + h1: Point3D_Array, + a1: Point3D_Array, +) -> Point3D_Array: ... + + +def get_quadratic_approximation_of_cubic(a0, h0, h1, a1): + r"""If :math:`a_0, h_0, h_1, a_1` are :math:`(3,)`-ndarrays representing control points + for a cubic Bézier curve, returns a :math:`(6, 3)`-ndarray of 6 control points + :math:`[a'_0, \ h', \ a'_1, \ a''_0, \ h'', \ a''_1]` for 2 quadratic Bézier curves + approximating it. + + If :math:`a_0, h_0, h_1, a_1` are :math:`(M, 3)`-ndarrays of :math:`M` control points + for :math:`M` cubic Bézier curves, returns instead a :math:`(6M, 3)`-ndarray of :math:`6M` + control points, where each one of the :math:`M` groups of 6 control points defines the 2 + quadratic Bézier curves approximating the respective cubic Bézier curve. + + .. note:: + The algorithm splits the cubic Bézier curve at an inflection point, if it contains + at least one. Otherwise, it simply splits it at :math:`t = 0.5`. Then, it computes + the two intersections between the tangent line at the split point and the tangents at + both ends of the original cubic Bézier: these points will be the handles for the two + new quadratic Bézier curves. + + .. seealso:: + `A Primer on Bézier Curves #21: Curve inflections. Pomax. `_ + + The inflection points in a cubic Bézier curve are those where the acceleration + (2nd derivative) is either zero or perpendicular to the velocity (1st derivative). + This can be expressed with a cross product in the following equation: + + .. math:: + C'(t) \times C''(t) = 0 + + The best way to solve this equation is by expressing :math:`C(t)` in its + polynomial form. If the control points :math:`a_0, h_0, h_1, a_1` allow us to + express it in its Bernstein form: + + .. math:: + C(t) = (1-t)^3 a_0 + 3(1-t)^2t h_0 + 3(1-t)t^2 h_1 + t^3 a_1 + + this can be rearranged into a polynomial form: + + .. math:: + C(t) = a_0 + 3t(h_0 - a_0) + 3t^2(h_1 - 2h_0 + a_0) + t^3(a_1 - 3h_1 + 3h_0 - a_0) + + where the following auxiliary points can be defined: + + .. math:: + p &= h_0 - a_0 \\ + q &= h_1 - 2h_0 + a_0 \\ + r &= a_1 - 3h_1 + 3h_0 - a_0 + + Thus, the velocity and acceleration can be easily derived: + + .. math:: + C(t) &= a_0 + 3tp + 3t^2q + t^3r \\ + C'(t) &= 3p + 6tq + 3t^2r \\ + C''(t) &= 6q + 6tr + + from where the original equation becomes: + + .. math:: + C'(t) \times C''(t) &= 0 \\ + (3p + 6tq + 3t^2r) \times (6q + 6tr) &= 0 \\ + 18(p \times q) + 18t(p \times r) + 36t(q \times q) + 36t^2(q \times r) + 18t^2(r \times q) + 18t^3(r \times r) &= 0 \\ + 18(t^2(q \times r) + t(p \times r) + (p \times q)) &= 0 \\ + t^2a + tb + c &= 0 + + which has a similar form to a quadratic equation, where one can define + :math:`a = q \times r, \ b = p \times r, \ c = p \times q`. + + * If the cubic Bézier curve is 2D, then :math:`a, b, c` are scalars, and thus + the 0, 1 or 2 values for :math:`t` can be derived from the quadratic equation. + One has to take care, however, that :math:`t \in [0, 1]`; otherwise, the + found solution must be discarded. + + * If the cubic Bézier curve is 3D, then :math:`a, b, c` are 3D vectors, + which implies solving 3 different quadratic equations. This case is more + complicated, as the solutions for the 3 equations should be the same + in order for the point located at those :math:`t` values to be actually + an inflection point. + + .. warning:: + The algorithm of this function assumes that the given cubic Bézier curve + is 2D, because it makes use of :func:`cross2D` and obtains :math:`a, b, c` + as scalars. For 3D curves the results might be wrong. + + + Parameters + ---------- + a0 + A :math:`(3,)` or :math:`(M, 3)`-ndarray of the start anchor(s) of the cubic Bézier curve(s). + h0 + A :math:`(3,)` or :math:`(M, 3)`-ndarray of the first handle(s) of the cubic Bézier curve(s). + h1 + A :math:`(3,)` or :math:`(M, 3)`-ndarray of the second handle(s) of the cubic Bézier curve(s). + a1 + A :math:`(3,)` or :math:`(M, 3)`-ndarray of the end anchor(s) of the cubic Bézier curve(s). + + Returns + ------- + :class:`~.Point3D_Array` + A :math:`(6M, 3)`-ndarray, where each one of the :math:`M` groups of + consecutive 6 points defines the 2 quadratic Bézier curves which + approximate the respective cubic Bézier curve. + """ + # If a0 is a Point3D, it's converted into a Point3D_Array of a single point: + # its shape is now (1, 3). + # If it was already a Point3D_Array of M points, it keeps its (M, 3) shape. + # Same with the other parameters. a0 = np.array(a0, ndmin=2) h0 = np.array(h0, ndmin=2) h1 = np.array(h1, ndmin=2) @@ -621,109 +1843,147 @@ def get_quadratic_approximation_of_cubic( 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) + # Search for inflection points. + # If no inflection points are found, use the midpoint as the split point instead. + # Based on https://pomax.github.io/bezierinfo/#inflections + t_split = np.full(a0.shape[0], 0.5) + # Let C(t) be a cubic Bézier curve defined by [a0, h0, h1, a1]. + # The below variables p, q, r allow for expressing the curve in a + # polynomial form convenient for derivatives: + # C(t) = a0 + 3tp + 3t²q + t³r p = h0 - a0 q = h1 - 2 * h0 + a0 r = a1 - 3 * h1 + 3 * h0 - a0 + # Velocity: C'(t) = 3p + 6tq + 3t²r + # Acceleration: C''(t) = 6q + 6tr + # C'(t) x C''(t) = 18 (t²(qxr) + t(pxr) + (pxq)) = 0 + # Define a = (qxr), b = (pxr) and c = (pxq). + # If C(t) is a 2D curve, then a, b and c are real numbers and + # this is a simple quadratic equation. + # However, if it's a 3D curve, then a, b and c are 3D vectors + # and this would require solving 3 quadratic equations. + # TODO: this simplifies by considering the 2D case. It might fail with 3D curves! 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 + # Case a == 0: degenerate 1st degree equation bt + c = 0 => t = -c/b + is_quadratic = a != 0 + is_linear = (~is_quadratic) & (b != 0) + t_split[is_linear] = -c[is_linear] / b[is_linear] + # Note: If a == 0 and b == 0, there are 0 or infinite solutions. + # Thus there are no inflection points. Just leave as is: t = 0.5 + + # Case a != 0: 2nd degree equation at² + bt + c = 0 + # => t = -(b/2a) +- sqrt((b/2a)² - c/a) + # Define u = b/2a and v = c/a, so that t = -u +- sqrt(u² - v). + u = 0.5 * b[is_quadratic] / a[is_quadratic] + v = c[is_quadratic] / a[is_quadratic] + radical = u * u - v + is_real = radical >= 0 + sqrt_radical = np.sqrt(radical[is_real]) + + t_minus = u[is_real] - sqrt_radical + t_plus = u[is_real] + sqrt_radical + is_t_minus_valid = (t_minus > 0) & (t_minus < 1) + is_t_plus_valid = (t_plus > 0) & (t_plus < 1) + + t_split[is_quadratic][is_real][is_t_minus_valid] = t_minus[is_t_minus_valid] + t_split[is_quadratic][is_real][is_t_plus_valid] = t_plus[is_t_plus_valid] + + # Compute Bézier point and tangent at the chosen value of t (these are vectorized) + t_split = t_split.reshape(-1, 1) + split_point = bezier([a0, h0, h1, a1])(t_split) # type: ignore + tangent_at_split = bezier([h0 - a0, h1 - h0, a1 - h1])(t_split) # 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) + i0 = find_intersection(a0, T0, split_point, tangent_at_split) + i1 = find_intersection(a1, T1, split_point, tangent_at_split) - m, n = np.shape(a0) - result = np.zeros((6 * m, n)) + M, dim = np.shape(a0) + result = np.empty((6 * M, dim)) result[0::6] = a0 result[1::6] = i0 - result[2::6] = mid - result[3::6] = mid + result[2::6] = split_point + result[3::6] = split_point result[4::6] = i1 result[5::6] = a1 return result def is_closed(points: Point3D_Array) -> bool: - return np.allclose(points[0], points[-1]) # type: ignore + """Returns ``True`` if the spline given by ``points`` is closed, by checking if its + first and last points are close to each other, or ``False`` otherwise. + + .. note:: + + This function reimplements :meth:`np.allclose` (without a relative + tolerance ``rtol``), because repeated calling of :meth:`np.allclose` + for only 2 points is inefficient. + + Parameters + ---------- + points + An array of points defining a spline. + + Returns + ------- + :class:`bool` + Whether the first and last points of the array are close enough or not + to be considered the same, thus considering the defined spline as closed. + """ + start, end = points[0], points[-1] + atol = 1e-8 + if abs(end[0] - start[0]) > atol: + return False + if abs(end[1] - start[1]) > atol: + return False + if abs(end[2] - start[2]) > atol: + return False + return True def proportions_along_bezier_curve_for_point( point: Point3D, control_points: BezierPoints, round_to: float = 1e-6, -) -> npt.NDArray[Any]: - """Obtains the proportion along the bezier curve corresponding to a given point - given the bezier curve's control points. +) -> Vector: + """Obtains the proportion along the Bézier curve corresponding to a given ``point``, + given the Bézier curve's ``control_points``. - The bezier polynomial is constructed using the coordinates of the given point - as well as the bezier curve's control points. On solving the polynomial for each dimension, - if there are roots common to every dimension, those roots give the proportion along the - curve the point is at. If there are no real roots, the point does not lie on the curve. + .. note:: + The Bézier polynomial is constructed using the coordinates of the given point, + as well as the Bézier curve's control points. On solving the polynomial for each dimension, + if there are roots common to every dimension, those roots give the proportion along the + curve the point is at. If there are no real roots, the point does not lie on the curve. Parameters ---------- point - The Cartesian Coordinates of the point whose parameter - should be obtained. + The cartesian coordinates of the point whose parameter should be obtained. control_points - The Cartesian Coordinates of the ordered control - points of the bezier curve on which the point may - or may not lie. + The cartesian coordinates of the ordered control points of the Bézier curve + on which the point may or may not lie. round_to - A float whose number of decimal places all values - such as coordinates of points will be rounded. + A float whose number of decimal places all values such as coordinates of + points will be rounded. Returns ------- - np.ndarray[float] - List containing possible parameters (the proportions along the bezier curve) - for the given point on the given bezier curve. - This usually only contains one or zero elements, but if the - point is, say, at the beginning/end of a closed loop, may return - a list with more than 1 value, corresponding to the beginning and - end etc. of the loop. + :class:`~.Vector` + Array containing possible parameters (the proportions along the Bézier curve) + for the given ``point`` on the given Bézier curve. + This array usually only contains one or zero elements, but if the + point is, say, at the beginning/end of a closed loop, this may contain more + than 1 value, corresponding to the beginning, end, etc. of the loop. Raises ------ :class:`ValueError` - When ``point`` and the control points have different shapes. + When ``point`` and the elements of ``control_points`` have different shapes. """ # Method taken from # http://polymathprogrammer.com/2012/04/03/does-point-lie-on-bezier-curve/ @@ -734,14 +1994,14 @@ def proportions_along_bezier_curve_for_point( ) control_points = np.array(control_points) - n = len(control_points) - 1 + d = len(control_points) - 1 roots = [] for dim, coord in enumerate(point): control_coords = control_points[:, dim] terms = [] - for term_power in range(n, -1, -1): - outercoeff = choose(n, term_power) + for term_power in range(d, -1, -1): + outercoeff = choose(d, term_power) term = [] sign = 1 for subterm_num in range(term_power, -1, -1): @@ -753,7 +2013,7 @@ def proportions_along_bezier_curve_for_point( sign *= -1 terms.append(outercoeff * sum(np.array(term))) if all(term == 0 for term in terms): - # Then both Bezier curve and Point lie on the same plane. + # Then both Bézier curve and point lie on the same plane. # Roots will be none, but in this specific instance, we don't need to consider that. continue bezier_polynom = np.polynomial.Polynomial(terms[::-1]) @@ -775,27 +2035,27 @@ def point_lies_on_bezier( control_points: BezierPoints, round_to: float = 1e-6, ) -> bool: - """Checks if a given point lies on the bezier curves with the given control points. + """Checks if a given ``point`` lies on the Bézier curve defined by ``control_points``. - This is done by solving the bezier polynomial with the point as the constant term; if - any real roots exist, the point lies on the bezier curve. + .. note:: + This is done by solving the Bézier polynomial with the point as the constant term. + If any real roots exist, then the point lies on the Bézier curve. Parameters ---------- point - The Cartesian Coordinates of the point to check. + The cartesian coordinates of the point to check. control_points - The Cartesian Coordinates of the ordered control - points of the bezier curve on which the point may - or may not lie. + The cartesian coordinates of the ordered control points of the Bézier curve on + which the point may or may not lie. round_to - A float whose number of decimal places all values - such as coordinates of points will be rounded. + A float whose number of decimal places all values such as coordinates of points + will be rounded. Returns ------- - bool - Whether the point lies on the curve. + :class:`bool` + Whether ``point`` lies on the Bézier curve or not. """ roots = proportions_along_bezier_curve_for_point(point, control_points, round_to) 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, + )