2
\$\begingroup\$

Recently I posted an answer on a question about Bézier curve calculations. As a micro-synopsis: there are three implementations of De Casteljau's algorithm here, including the original poster's, AJ Neufeld's, and my own. I've also added some non-recursive implementations based on a couple of different forms of the Bézier/Bernstein series, as well as a call into scipy's Bernstein implementation.

As an important attribution notice: those third-party implementations are included here because this code profiles them, but due to Code Review policy the first two should not be reviewed:

# de Casteljau's algorithm implementations
# This code courtesy @das-g


def dc_orig_curve(control_points, number_of_curve_points: int):
    return [
        dc_orig_point(control_points, t)
        for t in (
            i/(number_of_curve_points - 1) for i in range(number_of_curve_points)
        )
    ]


def dc_orig_point(control_points, t: float):
    if len(control_points) == 1:
        result, = control_points
        return result
    control_linestring = zip(control_points[:-1], control_points[1:])
    return dc_orig_point([(1 - t)*p1 + t*p2 for p1, p2 in control_linestring], t)


# This code courtesy @AJNeufeld


def dc_aj_curve(control_points, number_of_curve_points: int):
    last_point = number_of_curve_points - 1
    return [
        dc_aj_point(control_points, i / last_point)
        for i in range(number_of_curve_points)
    ]


def dc_aj_point(control_points, t: float):
    while len(control_points) > 1:
        control_linestring = zip(control_points[:-1], control_points[1:])
        control_points = [(1 - t) * p1 + t * p2 for p1, p2 in control_linestring]
    return control_points[0]

The remaining code is mine and can be reviewed:

# I adapted this code from @AJNeufeld's solution above


def dc_vec_curve(control_points, number_of_curve_points: int):
    last_point = number_of_curve_points - 1
    result = np.empty((number_of_curve_points, control_points.shape[1]))
    for i in range(number_of_curve_points):
        result[i] = dc_vec_point(control_points, i / last_point)
    return result


def dc_vec_point(control_points, t: float):
    while len(control_points) > 1:
        p1 = control_points[:-1]
        p2 = control_points[1:]
        control_points = (1 - t)*p1 + t*p2
    return control_points[0]


# The remaining code is not an adaptation.

def explicit_curve(control_points, n_curve_points: int):
    # https://en.wikipedia.org/wiki/B%C3%A9zier_curve#Explicit_definition
    n = len(control_points)
    # 0 <= t <= 1
    # B is the output
    # P_i are the control points

    binoms = [1]
    b = 1
    for k in range(1, (n + 1)//2):
        b = b*(n - k)//k
        binoms.append(b)

    mid = n//2 - 1
    binoms += binoms[mid::-1]

    t_delta = 1/(n_curve_points - 1)
    t = t_delta
    output = [None]*n_curve_points
    output[0] = control_points[0]
    output[-1] = control_points[-1]

    for ti in range(1, n_curve_points-1):
        B = 0
        tm = t/(1 - t)
        u = (1 - t)**(n - 1)

        for p, b in zip(control_points, binoms):
            B += b*u*p
            u *= tm

        output[ti] = B
        t += t_delta

    return output


def poly_curve(control_points, n_curve_points: int):
    # https://en.wikipedia.org/wiki/B%C3%A9zier_curve#Polynomial_form
    # B is the output
    n = len(control_points)
    # P is the array of control points
    # C is an array of coefficients
    # In the PI form,
    #    j is the index of the C coefficient
    #    m is the index of the factorial product
    #    i is the inner index of the sum

    C = [None]*n
    for j in range(n):
        product = 1
        for m in range(1, j + 1):
            product *= n - m

        total = 0
        for i, P in enumerate(control_points[:j+1]):
            addend = (1 if (i+j)&1 == 0 else -1)*P
            for f in range(2, i+1):
                addend /= f
            for f in range(2, j-i+1):
                addend /= f
            total += addend
        C[j] = product*total

    t_delta = 1/(n_curve_points - 1)
    t = t_delta
    output = [None]*n_curve_points
    output[0] = control_points[0]
    output[-1] = control_points[-1]

    for ti in range(1, n_curve_points-1):
        B = 0
        u = 1
        for c in C:
            B += u*c
            u *= t
        output[ti] = B
        t += t_delta

    return output


def bernvec_curve(control_points, n_curve_points: int):
    # scipy.interpolate.BPoly
    # This most closely resembles the "explicit" method because it calls comb()
    #     k: polynomial degree - len(control_points)
    #     m: number of breakpoints := 1
    #     n: coordinate dimensions := 2
    #     a: index of sum
    #     i: index into x
    #     x: m+1 array of "polynomial breakpoints"
    #     c: k * m * n array of polynomial coefficients
    #         - equal to control_points with an added dimension

    poly = BPoly(
        c=control_points[:, np.newaxis, :],
        x=[0, 1],
        extrapolate=False,
    )
    return poly(x=np.linspace(0, 1, n_curve_points))


def bernfun_curve(control_points, n_curve_points: int):
    # This does the same thing as bernvec, but bypasses the initialization and
    # validation layer
    out = np.empty((n_curve_points, 2), dtype=np.float64)
    evaluate_bernstein(
        c=control_points[:, np.newaxis, :],
        x=np.array([0., 1.]),
        xp=np.linspace(0, 1, n_curve_points),
        nu=0,
        extrapolate=False,
        out=out,
    )
    return out


# scipy.interpolate.BSpline
# https://github.com/numpy/numpy/issues/3845#issuecomment-227574158
# Bernstein polynomials are now available in SciPy as b-splines on a single interval.


def test():
    # degree 2, i.e. cubic Bézier with three control points per curve)
    # for large outputs (large number_of_curve_points)

    rng = np.random.default_rng(seed=0).random

    for n_controls in (2, 3, 6, 9):
        controls = rng((n_controls, 2), dtype=np.float64)
        for n_points in (2, 20, 2_000):
            expected = np.array(dc_orig_curve(controls, n_points), dtype=np.float64)

            for alt in (dc_aj_curve, dc_vec_curve, explicit_curve, poly_curve, bernvec_curve, bernfun_curve):
                actual = np.array(alt(controls, n_points), dtype=np.float64)
                assert actual.shape == expected.shape
                err = np.max(np.abs(expected - actual))
                print(f'nc={n_controls} np={n_points:4} method={alt.__name__:15} err={err:.1e}')
                assert err < 1e-12


class Profiler:
    MAX_CONTROLS = 10  # exclusive
    DECADES = 3
    PER_DECADE = 3
    N_ITERS = 30

    METHOD_NAMES = (
        'dc_orig',
        'dc_aj',
        'dc_vec',
        'explicit',
        'poly',
        'bernvec',
        'bernfun',
    )
    METHODS = {
        name: globals()[f'{name}_curve']
        for name in METHOD_NAMES
    }

    def __init__(self):
        self.all_control_points = default_rng().random((self.MAX_CONTROLS, 2), dtype=np.float64)
        self.control_counts = np.arange(2, self.MAX_CONTROLS, dtype=np.uint32)

        self.point_counts = np.logspace(
            0,
            self.DECADES,
            self.DECADES * self.PER_DECADE + 1,
            dtype=np.uint32,
        )

        self.quantiles = None

    def profile(self):
        times = np.empty(
            (
                len(self.control_counts),
                len(self.point_counts),
                len(self.METHODS),
                self.N_ITERS,
            ),
            dtype=np.float64,
        )

        times_vec = np.empty(self.N_ITERS, dtype=np.float64)

        for i, n_control in np.ndenumerate(self.control_counts):
            control_points = self.all_control_points[:n_control]
            for j, n_points in np.ndenumerate(self.point_counts):
                print(f'n_control={n_control} n_points={n_points})', end='\r')
                for k, method_name in enumerate(self.METHOD_NAMES):
                    method = lambda: self.METHODS[method_name](control_points, n_points)
                    for l in range(self.N_ITERS):
                        times_vec[l] = timeit(method, number=1)
                    times[i,j,k,:] = times_vec
        print()

        # Shape:
        #   Quantiles (3)
        #   Control counts
        #   Point counts
        #   Methods
        self.quantiles = np.quantile(times, (0.2, 0.5, 0.8), axis=3)

    def parametric_figure(
        self,
        x_series: np.ndarray,
        x_name: str,
        x_log: bool,
        z_series: np.ndarray,
        z_name: str,
        z_abbrev: str,
        colours: _ColorPalette,
    ):
        z_indices = (
            0,
            len(z_series)//2,
            -1,
        )

        fig: Figure
        axes: Sequence[Axes]
        fig, axes = pyplot.subplots(1, len(z_indices), sharey='all')
        fig.suptitle(f'Bézier curve calculation time, selected {z_name} counts')

        for ax, z in zip(axes, z_indices):
            ax.set_title(f'{z_abbrev}={z_series[z]}')

            if z == len(z_series) // 2:
                ax.set_xlabel(x_name)
            if z == 0:
                ax.set_ylabel('Time (s)')

            if x_log:
                ax.set_xscale('log')
            ax.set_yscale('log')
            ax.grid(axis='both', b=True, which='major', color='dimgray')
            ax.grid(axis='both', b=True, which='minor', color='whitesmoke')

            for i_method, method_name in enumerate(self.METHOD_NAMES):
                if z_abbrev == 'nc':
                    data = self.quantiles[:, z, :, i_method]
                elif z_abbrev == 'np':
                    data = self.quantiles[:, :, z, i_method]
                ax.plot(
                    x_series,
                    data[1, :],
                    label=method_name if z == 0 else '',
                    c=colours[i_method],
                )
                ax.fill_between(
                    x_series,
                    data[0, :],
                    data[2, :],
                    facecolor=colours[i_method],
                    alpha=0.3,
                )
        fig.legend()

    def plot(self):
        colours = color_palette('husl', len(self.METHODS))

        self.parametric_figure(
            x_series=self.point_counts,
            x_name='Point counts',
            x_log=True,
            z_series=self.control_counts,
            z_name='control',
            z_abbrev='nc',
            colours=colours,
        )
        self.parametric_figure(
            x_series=self.control_counts,
            x_name='Control counts',
            x_log=False,
            z_series=self.point_counts,
            z_name='point',
            z_abbrev='np',
            colours=colours,
        )

        pyplot.show()


if __name__ == '__main__':
    test()
    p = Profiler()
    p.profile()
    p.plot()

This produces two figures:

perf curves 1 perf curves 2

The approach to profiling in this case is to store method durations from timeit in an ndarray. timeit is being passed number=1 because I want to show variability in the graph as a shaded range between quantiles 0.2 and 0.8, with the line being on 0.5. The times array has the dimensions:

  • Control point count (8)
  • Curve output point count (10)
  • Method (3)
  • Iteration (30)

Specifically I am interested in whether I made appropriate and efficient use of matplotlib and numpy, and whether my own implementations of the curve functions make sense and are efficient.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.