codexbezier-de-casteljau-construction.py
bezier-de-casteljau-construction.py
"""Cubic Bézier hull and midpoint construction; output PNG in the working directory.
Uses the repository's NumPy and Matplotlib dependencies; tested with Python 3.14.4.
"""
from pathlib import Path
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon, Rectangle
import numpy as np
def bernstein_cubic(points, t):
t = np.asarray(t)
weights = np.stack(((1 - t) ** 3, 3 * (1 - t) ** 2 * t,
3 * (1 - t) * t**2, t**3), axis=-1)
return weights @ points
def construction(points, t):
levels = [np.asarray(points, dtype=float)]
while len(levels[-1]) > 1:
row = levels[-1]
levels.append((1 - t) * row[:-1] + t * row[1:])
return levels
def main():
points = np.array([[0, 0], [1, 2], [3, 2], [4, 0]], dtype=float)
p, q, r, s = construction(points, 0.5)
np.testing.assert_array_equal(q, [[0.5, 1], [2, 2], [3.5, 1]])
np.testing.assert_array_equal(r, [[1.25, 1.5], [2.75, 1.5]])
np.testing.assert_array_equal(s, [[2, 1.5]])
t = np.linspace(0, 1, 1001)
curve = bernstein_cubic(points, t)
expanded = np.column_stack((3 * t + 3 * t**2 - 2 * t**3, 6 * t * (1 - t)))
np.testing.assert_allclose(curve, expanded, rtol=0, atol=2e-15)
recursive = np.array([construction(points, float(u))[-1][0] for u in t])
evaluation_error = float(np.max(np.abs(curve - recursive)))
assert evaluation_error < 2e-15
# The trapezoidal hull is 0 <= y <= min(2, 2x, 8 - 2x).
x, y = curve.T
assert np.all(y >= 0) and np.all(y <= 1.5)
assert np.all(y <= 2 * x + 1e-14) and np.all(y <= 8 - 2 * x + 1e-14)
assert 2 > 2 * 0 # Bounding-box corner (0, 2) is outside that hull.
matrix = np.array([[2, 1], [0, 1]], dtype=float)
shift = np.array([0, 1], dtype=float)
transformed = points @ matrix.T + shift
np.testing.assert_array_equal(transformed, [[0, 1], [4, 3], [8, 3], [8, 1]])
affine_error = float(np.max(np.abs(
bernstein_cubic(transformed, t) - (curve @ matrix.T + shift))))
assert affine_error < 4e-15
np.testing.assert_array_equal(bernstein_cubic(transformed, 0.5), [5.5, 2.5])
plt.rcParams.update({"font.size": 11, "axes.titlesize": 13,
"axes.labelsize": 11, "figure.facecolor": "white",
"axes.facecolor": "white"})
fig, axes = plt.subplots(1, 2, figsize=(11, 5.25), dpi=120)
blue, orange, green, purple = "#155b8b", "#ba570d", "#16744b", "#79439b"
ax = axes[0]
ax.add_patch(Polygon(points, closed=True, facecolor="#edf0f3",
edgecolor="#98a2ac", lw=1, label="Convex hull"))
ax.add_patch(Rectangle((0, 0), 4, 2, fill=False, ls=":", lw=1.8,
edgecolor="#606c78", label="Bounding box"))
ax.plot(*points.T, "o--", color="#687582", lw=1.5, ms=6,
label="Control polygon", zorder=3)
ax.plot(*curve.T, color=blue, lw=2.8, label="Bézier curve", zorder=4)
ax.scatter([2], [1.5], s=45, color=blue, zorder=5)
for i, point in enumerate(points):
offset = [(0, -16), (-12, 15), (12, 15), (0, -16)][i]
ax.annotate(f"$P_{i}$ ({point[0]:g}, {point[1]:g})", point,
xytext=offset, textcoords="offset points", ha="center")
ax.annotate("Box corner lies\noutside the hull", (0, 2), xytext=(0.15, 2.65),
fontsize=10, arrowprops={"arrowstyle": "->", "color": "#606c78"})
ax.annotate("$C(1/2)=(2,1.5)$", (2, 1.5), xytext=(2, 0.65),
ha="center", color=blue,
arrowprops={"arrowstyle": "->", "color": blue})
ax.set_title("Curve, controls and two different enclosures", pad=16)
ax.legend(loc="lower center", bbox_to_anchor=(0.5, -0.38), ncol=2,
frameon=False, fontsize=10)
ax = axes[1]
ax.plot(*curve.T, color=blue, lw=2.4, label="Bézier curve", zorder=3)
ax.plot(*p.T, "o--", color="#a1aab3", lw=1.2, ms=5, label="Original controls")
ax.plot(*q.T, "o-", color=orange, lw=1.8, ms=6, label="Level 1: Q")
ax.plot(*r.T, "s-", color=green, lw=2.1, ms=6, label="Level 2: R")
ax.scatter(*s.T, color=purple, marker="D", s=70, zorder=7, label="Level 3: S")
for i, point in enumerate(p):
ax.annotate(f"$P_{i}$", point, xytext=(0, 13 if i in (1, 2) else -15),
textcoords="offset points", ha="center", color="#687582")
for i, point in enumerate(q):
offsets = [(-22, -24), (0, 22), (22, -24)]
ax.annotate(f"$Q_{i}$\n({point[0]:g}, {point[1]:g})", point,
xytext=offsets[i], textcoords="offset points", ha="center",
fontsize=10, color=orange)
for i, point in enumerate(r):
ax.annotate(f"$R_{i}$\n({point[0]:g}, {point[1]:g})", point,
xytext=(0, 9),
textcoords="offset points", ha="center", fontsize=10, color=green,
bbox={"facecolor": "white", "edgecolor": "none", "pad": 1})
ax.annotate("$S=(2,1.5)$", s[0], xytext=(0, -40),
textcoords="offset points", ha="center", color=purple,
arrowprops={"arrowstyle": "->", "color": purple})
ax.set_title("Three levels: take each neighboring midpoint", pad=16)
ax.legend(loc="lower center", bbox_to_anchor=(0.5, -0.38), ncol=2,
frameon=False, fontsize=10)
for ax in axes:
ax.set(xlim=(-0.45, 4.45), ylim=(-0.35, 2.9), xlabel="x", ylabel="y")
ax.set_aspect("equal")
ax.set_xticks(range(5))
ax.set_yticks([0, 1, 2])
ax.spines[["top", "right"]].set_visible(False)
fig.subplots_adjust(left=0.055, right=0.98, top=0.91, bottom=0.26, wspace=0.18)
output = Path.cwd() / (Path(__file__).stem + ".png")
fig.savefig(output, facecolor="white", transparent=False)
plt.close(fig)
print(f"{output.name}: 1320x630; midpoint={s[0].tolist()}; "
f"1001-point Bernstein/recursion error={evaluation_error:.3g}; "
f"affine error={affine_error:.3g}; hull inequalities passed")
if __name__ == "__main__":
main()