Arc length parameterization of a Bezier curve.

bezier.py

"""
A Bezier curve implementation.

Requires 
- numpy >= 1.26.4
"""
from dataclasses import dataclass
from math import comb

import numpy as np
from numpy.typing import NDArray


@dataclass
class BezierCurve:
    """
    A Bezier curve.

    Parameters
    ----------
    control_points : NDArray[np.float32]
        Array of control points of Bezier curve with shape `(N, 2)`.
    arc_length_approximation : int, default: 50
        Number of evaluations for arc length approximation.

    Attributes
    ----------
    arc_length : float
        Approximate length of Bezier curve.
    arc_length_approximation : int
        Number of evaluations for arc length approximation.
    arc_lengths : NDArray[np.float32]
        Approximate arc lengths along Bezier curve.
    coef : NDArray[np.float32]
        Binomial coefficients of Bezier curve.
    control_points : NDArray[np.float32]
        Array of control points of Bezier curve with shape `(N, 2)`.
    degree : int
        Degree of Bezier curve.

    Methods
    -------
    evaluate(t)
        Evaluate the Bezier curve at `t` (0 <= t <= 1).
    arc_length_proportion(p)
        Evaluate the Bezier curve at a proportion of its total arc length.
    """

    control_points: NDArray[np.float32]
    """Array of control points of Bezier curve with shape `(N, 2)`."""
    arc_length_approximation: int = 50
    """Number of evaluations for arc length approximation."""

    def __post_init__(self):
        if self.degree == -1:
            raise ValueError("There must be at least one control point.")

        self.coef: NDArray[np.float32] = np.array(
            [comb(self.degree, i) for i in range(self.degree + 1)], dtype=float
        )
        """Binomial coefficients of Bezier curve."""

        evaluated = self.evaluate(np.linspace(0, 1, self.arc_length_approximation))
        norms = np.linalg.norm(evaluated[1:] - evaluated[:-1], axis=-1)
        self.arc_lengths: NDArray[np.float32] = np.append(0, norms.cumsum())
        """Approximate arc lengths along Bezier curve."""

    @property
    def degree(self) -> int:
        """Degree of Bezier curve."""
        return len(self.control_points) - 1

    @property
    def arc_length(self) -> float:
        """Approximate length of Bezier curve."""
        return self.arc_lengths[-1]

    def evaluate(self, t: float | NDArray[np.float32]) -> NDArray[np.float32]:
        """Evaluate the Bezier curve at `t` (0 <= t <= 1)."""
        t = np.asarray(t)
        terms = np.logspace(0, self.degree, num=self.degree + 1, base=t).T
        terms *= np.logspace(self.degree, 0, num=self.degree + 1, base=1 - t).T
        terms *= self.coef
        return terms @ self.control_points

    def arc_length_proportion(self, p: float) -> NDArray[np.float32]:
        """Evaluate the Bezier curve at a proportion of its total arc length."""
        target_length = self.arc_length * p
        i = np.searchsorted(self.arc_lengths, target_length, "right") - 1

        previous_length = self.arc_lengths[i]
        if previous_length == target_length:
            return self.evaluate(i / self.arc_length_approximation)

        target_dif = target_length - previous_length
        arc_dif = self.arc_lengths[i + 1] - previous_length

        t = (i + target_dif / arc_dif) / self.arc_length_approximation
        return self.evaluate(t)


if __name__ == "__main__":
    import matplotlib.pyplot as plt

    b = BezierCurve(np.random.random((5, 2)))
    plt.plot(*b.evaluate(np.linspace(0, 1, 51)).T)
    equi = [b.arc_length_proportion(i / 10) for i in range(11)]
    xs = [x for x, _ in equi]
    ys = [y for _, y in equi]
    plt.scatter(xs, ys, alpha=0.75)
    plt.scatter(
        *b.evaluate(np.linspace(0, 1, 11)).T, color=np.array([1.0, 0.0, 0.0, 0.5])
    )
    plt.show()