# Profiling for Bézier curve calculations

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):
for f in range(2, j-i+1):
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
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,
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:

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.