codexbezier-subdivision-and-tangents.py
bezier-subdivision-and-tangents.py
"""Exact cubic subdivision and tangents; 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
import numpy as np
def 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 derivative(points, t):
t = np.asarray(t)
weights = np.stack(((1 - t) ** 2, 2 * t * (1 - t), t**2), axis=-1)
return weights @ (3 * np.diff(points, axis=0))
def main():
points = np.array([[0, 0], [1, 2], [3, 2], [4, 0]], dtype=float)
q = (points[:-1] + points[1:]) / 2
r = (q[:-1] + q[1:]) / 2
s = (r[0] + r[1]) / 2
left = np.array([points[0], q[0], r[0], s])
right = np.array([s, r[1], q[2], points[3]])
np.testing.assert_array_equal(left, [[0, 0], [0.5, 1], [1.25, 1.5], [2, 1.5]])
np.testing.assert_array_equal(right, [[2, 1.5], [2.75, 1.5], [3.5, 1], [4, 0]])
u = np.linspace(0, 1, 1001)
left_error = float(np.max(np.abs(cubic(left, u) - cubic(points, u / 2))))
right_error = float(np.max(np.abs(cubic(right, u) - cubic(points, (1 + u) / 2))))
assert left_error < 3e-15 and right_error < 3e-15
np.testing.assert_allclose(derivative(left, u), derivative(points, u / 2) / 2,
rtol=0, atol=3e-15)
np.testing.assert_allclose(derivative(right, u), derivative(points, (1 + u) / 2) / 2,
rtol=0, atol=3e-15)
np.testing.assert_array_equal(derivative(points, 0), [3, 6])
np.testing.assert_array_equal(derivative(points, 1), [3, -6])
np.testing.assert_array_equal(derivative(left, 1), [2.25, 0])
np.testing.assert_array_equal(derivative(right, 0), [2.25, 0])
speeds = np.linalg.norm(derivative(points, np.array([0, 0.5])), axis=1)
np.testing.assert_allclose(speeds, [np.sqrt(45), 4.5], rtol=0, atol=1e-14)
sample_t = np.linspace(0, 1, 5)
samples = cubic(points, sample_t)
chords = np.linalg.norm(np.diff(samples, axis=0), axis=1)
assert chords[0] > chords[1] + 0.2
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, dark = "#155b8b", "#ba570d", "#263c50"
ax = axes[0]
ax.plot(*points.T, "--", color="#c0c6cc", lw=1.2, label="Original polygon")
ax.plot(*left.T, "o--", color=blue, lw=1.5, ms=6, label="Left control polygon")
ax.plot(*right.T, "s--", color=orange, lw=1.5, ms=6, label="Right control polygon")
ax.plot(*cubic(left, u).T, color=blue, lw=3)
ax.plot(*cubic(right, u).T, color=orange, lw=3)
for label, point, offset, color in [
("$P_0$", left[0], (0, -15), blue),
("$Q_0$", left[1], (-18, 11), blue),
("$R_0$", left[2], (-8, 15), blue),
("$R_1$", right[1], (8, 15), orange),
("$Q_2$", right[2], (18, 11), orange),
("$P_3$", right[3], (0, -15), orange),
]:
ax.annotate(label, point, xytext=offset, textcoords="offset points",
ha="center", color=color)
ax.scatter(*s, s=60, color=dark, zorder=5)
ax.annotate("Shared endpoint\n$S=(2,1.5)$", s, xytext=(2, 0.6),
ha="center", fontsize=11, color=dark,
arrowprops={"arrowstyle": "->", "color": dark})
ax.set_title("Subdivision changes the controls, preserves the curve", fontsize=12, pad=16)
ax.legend(loc="lower center", bbox_to_anchor=(0.5, -0.39),
frameon=False, fontsize=10)
ax = axes[1]
ax.plot(*points.T, "o--", color="#a1aab3", lw=1.2, ms=5, label="Control polygon")
ax.plot(*cubic(points, u).T, color=dark, lw=2.6, label="Bézier curve")
for start, end, color in [
(points[0], points[0] + 0.18 * derivative(points, 0), blue),
(points[3] - 0.18 * derivative(points, 1), points[3], orange),
]:
ax.annotate("", end, xytext=start,
arrowprops={"arrowstyle": "->", "lw": 2.7, "color": color}, zorder=6)
ax.scatter(*samples.T, color=dark, s=36, zorder=5)
for value, point, offset in zip(sample_t, samples,
[(0, -15), (-20, 13), (0, 20), (20, 13), (0, -15)]):
ax.annotate(f"$t={value:g}$", point, xytext=offset,
textcoords="offset points", ha="center", fontsize=10)
ax.text(0.45, 2.55, "$C'(0)=(3,6)$", color=blue, ha="center")
ax.text(3.55, 2.55, "$C'(1)=(3,-6)$", color=orange, ha="center")
ax.set_title("Endpoint tangents and equal steps in t", pad=16)
ax.text(0.5, -0.28, "Tangent arrows show direction (lengths scaled).\n"
r"Speed: $\|C'(0)\|=\sqrt{45}$, $\|C'(1/2)\|=4.5$."
"\nEqual parameter steps need not give equal distances.",
transform=ax.transAxes, ha="center", va="top", fontsize=10, linespacing=1.5)
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; 1001-point subdivision errors "
f"left={left_error:.3g}, right={right_error:.3g}; "
f"endpoint tangents=(3,6),(3,-6); local join tangents=(2.25,0); "
f"speeds={speeds.tolist()}; sampled chord lengths={chords.tolist()}")
if __name__ == "__main__":
main()