pure Python Bézier curve implementation

Question

I came up with this recursive pure-Python implementation of De Casteljau's algorithm for computing points on a Bézier curve:

def bezier_curve(control_points, number_of_curve_points):
    return [
        bezier_point(control_points, t)
        for t in (
            i / (number_of_curve_points - 1) for i in range(number_of_curve_points)
        )
    ]


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

Assumptions about control_points

The elements of control_points represent the control points of the Bézier curve. They must be of the same type of mutually compatible types fulfilling the following rules:

• points shall all be of the same dimension.
• Multiplying a point by a scalar shall result in a point of the same dimension and with a value according to vector-scalar multiplication (i.e., multiply each of the point's Cartesian coordinates with the scalar)
• Adding two points shall result in a point of the same dimension and with a value according to vector addition (i.e., component-wise addition of the points' Cartesian coordinates)

Some examples that work as control_points:

• list of turtle.Vec2D
• list of complex
• list of numpy.array with shape (2,)
• numpy.array with shape (n, 2), where n is the number of control points

(tuples instead of lists work, too. Probably any sequencey container will work.)

Why pure Python?

Because I want this to be usable in a QGIS plugin, but SciPy, NumPy, etc. (usually) aren't available to QGIS plugins. As I'm unsure which Python libraries are available in QGIS (and the answer to that seems to be platform-dependent), I'd like to avoid external libraries (those that would have to be installed with pip or one of its alternatives) completely.

Using standard library functions should be fine, so if any part of the implementation could benefit from them, please point that out.

What I'd like to know in this review

• Could / should the readability and comprehensibility of this implementation be improved?
• Did I go foul of any performance (computation speed, memory usage, etc.) no gos? (It doesn't need to be super-fast, but it shouldn't be unnecessarily slow if I can avoid it.)
  • Performance for low degrees (e.g., degree 2, i.e. cubic Bézier with three control points per curve) will probably be more relevant than performance at high degrees (many control points per curve)
  • Performance for large outputs (large number_of_curve_points) might be relevant
• About the destructuring assignment result, = control_points to unpack the single point while at the same time making sure it really is exactly one point
  • Is this idiomatic in Python (i.e. "pythonic")?
  • Is this readable and comprehensible enough or too obscure?
  • Is there any good alternative that's an expression, i.e. that could be used directly in the return statement without going through an assignment? (control_points[0] is an expression but doesn't fail when there's more than one element in control_points.)

Answer 1

Unnecessary Generator

You've got an unnecessary generator expression here:

def bezier_curve(control_points, number_of_curve_points):
    return [
        bezier_point(control_points, t)
        for t in (
            i / (number_of_curve_points - 1) for i in range(number_of_curve_points)
        )
    ]

You don't need to generate i / (n-1); you could simply pass that argument to the bezier_point() function:

def bezier_curve(control_points, number_of_curve_points):
    return [ bezier_point(control_points, i / (number_of_curve_points - 1))
             for i in range(number_of_curve_points)
           ]

Slight optimization: instead of computing number_of_curve_points - 1 \$O(N)\$ times (pure Python will not cache the result), precompute it:

def bezier_curve(control_points, number_of_curve_points):
    last_point = number_of_curve_points - 1
    return [ bezier_point(control_points, i / last_point )
             for i in range(number_of_curve_points)
           ]

Tail Recursion

Python does not do Tail Call Optimization, so with M control points, you will recursively enter and exit M calls, for each of the N points along your curve. That is M*N unnecessary stack frame entry/exits. You should do the looping yourself:

def bezier_point(control_points, t):
    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]

Since we loop while len(control_points) > 1, it should be guaranteed that control_points will only have one point when the loop exits, so return control_points[0] is safe. The exception is if the function is called with zero control points, but then control_points[0] will properly fail with an IndexError.

Answer 2

About this code:

def bezier_point(control_points, t):
    if len(control_points) == 1:
        result, = control_points  # <-- here

you ask:

Is [the tuple-unpacking] idiom Pythonic?

Yes.

Is it comprehensible?

Yes.

Does the Python standard library offer any handy alternatives for this?

Technically there is operator.itemgetter but I do not recommend that you use that in this case. For one thing it would only provide an equivalent to control_points[0], without effectively asserting for length.

Is there a way that is itself a single expression, so that it can be used inline in other expressions (e.g., in lambdas or in list comprehensions)?

To put unpacking as an expression on the right-hand side of an assignment, no, this effectively can't be done without a really silly comprehension hack:

next(iter(cp for (cp,) in (control_points,)))

Please do not do this. Doing anything more complicated than what you have now (for instance defining your own "unpacking function") is not advisable.

The exception might be if you also want to do some of your own validation, i.e. wrapping an exception in your own:

def get_only_point(control_points: Iterable[float]) -> float:
    try:
        point, = control_points
    except ValueError as e:
        raise MyDataError('too many control points') from e
    return point

Answer 3

Going in quite a different direction: let's see exactly how much AJ's improvements make a difference, and how and why to vectorize. I know you claim that

SciPy, NumPy, etc. (usually) aren't available to QGIS plugins

but given these results, it would be worth doing a

try:
    import numpy as np
except ImportError:
    # sad face
    from .fallbacks import *

In other words, keep both vectorized and non-vectorized implementations, using the best one possible.

This (somewhat hacky) profiling code:

#!/usr/bin/env python3

from matplotlib import pyplot
from matplotlib.axes import Axes
from matplotlib.figure import Figure
from numpy.random._generator import default_rng
from seaborn import color_palette
from timeit import timeit
from typing import List, Sequence
import numpy as np


def original_curve(control_points, number_of_curve_points):
    return [
        original_point(control_points, t)
        for t in (
            i/(number_of_curve_points - 1) for i in range(number_of_curve_points)
        )
    ]


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


def aj_curve(control_points, number_of_curve_points):
    last_point = number_of_curve_points - 1
    return [
        aj_point(control_points, i / last_point)
        for i in range(number_of_curve_points)
    ]


def aj_point(control_points, t):
    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]


def vectorized_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] = vectorized_point(control_points, i / last_point)
    return result


def vectorized_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]


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

    controls = np.random.default_rng().random((3, 2), dtype=np.float64)
    n_points = 10_000

    expected: List[complex] = original_curve(controls, n_points)

    for alt in (aj_curve, vectorized_curve):
        actual = alt(controls, n_points)
        assert np.isclose(expected, actual).all()


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

    METHOD_NAMES = (
        'original',
        'aj',
        'vectorized',
    )
    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 control_figures(self, colours):
        control_indices = (
            0,
            len(self.control_counts) // 2,
            -1,
        )

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

        for ax, i_control in zip(axes, control_indices):
            n_control = self.control_counts[i_control]
            ax.set_title(f'nc={n_control}')
            if i_control == len(self.control_counts) // 2:
                ax.set_xlabel('Curve points')
            if i_control == 0:
                ax.set_ylabel('Time (s)')

            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):
                data = self.quantiles[:, i_control, :, i_method]
                ax.plot(
                    self.point_counts,
                    data[1, :],
                    label=method_name if i_control == 0 else '',
                    c=colours[i_method],
                )
                ax.fill_between(
                    self.point_counts,
                    data[0, :],
                    data[2, :],
                    facecolor=colours[i_method],
                    alpha=0.3,
                )
        fig.legend()

    def point_figures(self, colours):
        point_indices = (
            0,
            len(self.point_counts)//2,
            -1,
        )

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

        for ax, i_point in zip(axes, point_indices):
            n_points = self.point_counts[i_point]
            ax.set_title(f'np={n_points}')

            if i_point == len(self.point_counts) // 2:
                ax.set_xlabel('Control points')
            if i_point == 0:
                ax.set_ylabel('Time (s)')

            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):
                data = self.quantiles[:, :, i_point, i_method]
                ax.plot(
                    self.control_counts,
                    data[1, :],
                    label=method_name if i_point == 0 else '',
                    c=colours[i_method],
                )
                ax.fill_between(
                    self.control_counts,
                    data[0, :],
                    data[2, :],
                    facecolor=colours[i_method],
                    alpha=0.3,
                )
        fig.legend()

    def plot(self):
        colours = color_palette('husl', len(self.METHODS))
        self.control_figures(colours)
        self.point_figures(colours)
        pyplot.show()


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

produces these:

I didn't give this profiling a lot of CPU time so the results are a little bumpy (inter-quantile shading shown between 0.2 and 0.8), but quite clear. Vectorization is definitely worth doing, even if it can't always be done. Some efficiencies might be found on top of what I have shown because I am not a Numpy expert.