贝塞尔曲线的 C 语言实现：用有限差分避免内层乘法

问题背景

在单片机上做轨迹插值、路径规划或简单图形渲染时，经常需要按等参数步长计算三次贝塞尔曲线上的若干点。最直接的写法是每次迭代都做一次三次多项式求值，涉及多次乘法和浮点运算。在 200 MHz 左右、没有硬件 FPU 的 MCU 上，如果点数较多（比如上千点），纯乘法方案会明显拖慢实时性。

本文采用的思路是有限差分（差值表）：把三次多项式展开后，利用差分递推，让内层循环只剩加法，从而大幅降低单步耗时。

原理简述

三次贝塞尔曲线

$B(t) = (1-t)^3 P_0 + 3(1-t)^2 t\, P_1 + 3(1-t) t^2 P_2 + t^3 P_3$

展开后是关于 $t$ 的三次多项式：

$B(t) = a + b\,t + c\,t^2 + d\,t^3$

其中（以 $x$ 分量为例，$y$ 同理）：

• $a = x_0$
• $b = 3(x_1 - x_0)$
• $c = 3(x_0 - 2x_1 + x_2)$
• $d = -x_0 + 3x_1 - 3x_2 + x_3$

取步长 $h = 1/n$（$n$ 为输出点数），定义一阶差分 $\Delta f(t) = f(t+h) - f(t)$，可以推出：

递推规则（每步只做加法）：

f     += df
df    += ddf
ddf   += dddf   // 二阶差分每步递增一个常数 dddf
dddf  += 0      // 三阶差分恒定，无需更新（此行仅为示意）

这样内层循环没有乘法，只有几次浮点加法。

代码实现

void cubic_bezier(int pointNum, float x2, float y2, float x3, float y3, float* out_data_x, float* out_data_y) {
	// 曲线点的个数

	float t = 1.0 / pointNum;
	float t2 = t * t;
	float t3 = t * t2;

	float x1 = 0, y1 = 0;
	float x4 = 1, y4 = 1;

	// 预计算
	float pre_3t = 3 * t;
	float pre_3t2 = 3 * t2;
	float pre_6t2 = 6 * t2;
	float pre_6t3 = 6 * t3;

	// (p1 - 2p2 + p3)
	float tmp1x = x1 - x2 * 2.0 + x3;
	float tmp1y = y1 - y2 * 2.0 + y3;

	// (-p1 + 3p2 - 3p3 + p4)
	float tmp2x = (x2 - x3) * 3.0 - x1 + x4;
	float tmp2y = (y2 - y3) * 3.0 - y1 + y4;

	// p1
	float fx = x1;
	float fy = y1;

	// (p2 - p1) * 3t +
	// (p1 - 2p2 + p3) * 3t2  +
	// (-p1 + 3p2 - 3p3 + p4) * t3
	float dfx = (x2 - x1) * pre_3t + tmp1x * pre_3t2 + tmp2x * t3;
	float dfy = (y2 - y1) * pre_3t + tmp1y * pre_3t2 + tmp2y * t3;

	// (p1 - 2p2 + p3) * 6t2  +
	// (-p1 + 3p2 - 3p3 + p4) * 6t3
	float ddfx = tmp1x * pre_6t2 + tmp2x * pre_6t3;
	float ddfy = tmp1y * pre_6t2 + tmp2y * pre_6t3;

	// (-p1 + 3p2 - 3p3 + p4) * 6t3
	float dddfx = tmp2x * pre_6t3;
	float dddfy = tmp2y * pre_6t3;

	// fx, fy 就是t值对应计算出的曲线值
	while (pointNum--) {
		fx += dfx;
		fy += dfy;
		dfx += ddfx;
		dfy += ddfy;
		ddfx += dddfx;
		ddfy += dddfy;

		*out_data_x = fx;
		*out_data_y = fy;
		out_data_x++;
		out_data_y++;
	}
}

代码要点

• 端点固定：$P_0 = (0, 0)$，$P_3 = (1, 1)$ 被硬编码在函数内部。两个控制点 $P_1 = (x_2, y_2)$、$P_2 = (x_3, y_3)$ 通过参数传入。这意味着该函数生成的曲线起点和终点固定为 $(0,0)$ 和 $(1,1)$；如果需要任意端点，需将 $x_1, y_1, x_4, y_4$ 也改为参数并调整预计算部分。
• 输出起点：循环第一次迭代输出的是 $B(h)$，即 $t = 1/n$ 处的点。$t = 0$ 处的起点 $(0,0)$ 不在输出数组中。如果调用方需要包含起点，需自行在数组头部补一个 $(0,0)$。
• 输出终点：最后一次迭代对应 $t = n \cdot h = 1$，理论上输出 $(1,1)$。由于浮点累加误差，实际值可能与 $1.0$ 有微小偏差。
• 预计算：$3h$、$3h^2$、$6h^2$、$6h^3$ 以及 $c$、$d$ 的系数在循环外一次性算好，循环体内只有浮点加法、存储和指针递增。

适用场景与限制

• 适用：主频不高、无硬件 FPU（或 FPU 性能有限）的 MCU；需要一次性生成大量等参数间隔的曲线点；曲线端点可归一化到 $(0,0)$ 和 $(1,1)$。
• 不适用 / 需注意：
  • 端点不是 $(0,0)$ 和 $(1,1)$ 时，需要修改函数签名并调整预计算部分。
  • pointNum 必须大于 0。若为 0，1.0 / pointNum 在 IEEE 754 下产生 +inf，后续预计算将得到 inf/NaN（C 标准将浮点除零归为实现定义行为，可能触发 FE_DIVBYZERO 信号）；若为负数，步长为负，输出无意义。调用前应做参数校验。
  • out_data_x、out_data_y 指向的缓冲区必须至少有 pointNum 个 float 的空间，函数内部不做越界检查。
  • 使用 float（32 位）累加，点数很大时（例如数万点）累积误差可能不可忽略。如果对精度敏感，可考虑用 double（前提是 MCU 支持）或在关键节点用完整多项式重新计算一次。
  • 该实现是单线程、无锁的。如果多个任务同时调用，需自行加锁或分配独立缓冲区。

性能参考

本文最初于 2023/8/30 发布在 CSDN。