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\$

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.