简介:这是一份面向C语言开发者的B样条曲线算法实现,完整覆盖二次与三次样条曲线计算,能够用于曲线拟合、数据平滑以及二维平面/三维空间曲线生成。核心逻辑被封装为两个函数:一个针对给定三点或四点计算出样条曲线上的平滑点,另一个负责对一系列散点生成整体平滑曲线;样条阶次与平滑程度均通过函数参数灵活配置,方便适配不同精度要求的图形学、CAD辅助设计、路径规划等任务。资源包共6个文件,包含3个cpp源码、2个h头文件以及1个bmp效果示意图,压缩包仅6KB,体量轻巧、结构清晰。源码全程带必要注释,并附有独立测试程序和一个完整使用案例,读者可循着示例快速理解二次/三次样条曲线的计算流程,也可直接将函数移植到自己的项目中使用。目前已有1112人学习下载,适合正在钻研样条曲线算法、需要平滑点列或曲线拟合实现参考的开发者和相关专业学生。 做嵌入式图形界面或者运动控制的时候,谁没被“把一堆离散点连成光滑曲线”这件事折磨过?鼠标手写笔画出来全是锯齿、传感器采集的数据抖动得像心电图、数控加工轨迹在拐角处突然顿一下……这些问题背后基本都能归结到同一个需求:曲线拟合与曲线平滑。我当时就是被手写板笔迹平滑逼着去啃B样条曲线的,最后用纯C语言写了一套样条曲线算法实现代码,覆盖二次样条曲线、三次样条曲线、曲线拟合和曲线平滑,今天把这套东西完整拆开讲一遍。
这套代码的适用场景很明确:单片机、嵌入式Linux、工控设备、上位机图形处理,凡是不能随便引第三方图形库、又需要高性能样条曲线计算的场合,都可以直接参考。如果你正在学C语言,或者要在项目里做轨迹规划、数据平滑,这篇文章就是冲着你写的。
1. 理解B样条:先搞懂这算法解决什么问题
1.1 从贝塞尔到B样条:为什么非得用这个算法
说到曲线拟合,很多人第一反应是贝塞尔曲线。贝塞尔确实好理解——用控制点拉出曲线,移动一个点整条曲线跟着变。但这个“全局性”恰恰是它在工程里的软肋:你只想调整一个局部的形状,结果整条曲线全变形了,这在手写笔画平滑、轨迹规划这种场景里非常致命。
B样条曲线正好解决了这个问题。它的核心思想是“分段拼接”:把参数域切成若干段,每一段由少数几个控制点决定,段与段之间保持光滑连接。一个控制点移动了,影响的只是附近的几个分段,远处基本不受影响。这个性质叫局部支撑性,也是B样条成为工程首选的最重要原因。
我记得第一次看懂B样条的局部支撑性时,脑子里冒出来一个画面:控制点像一排钉子,曲线像一根弹性绳在钉子之间穿行。挪动一颗钉子,只有挨着它的那截绳子会动,远处的绳子纹丝不动。这个类比虽然不够严谨,但用来理解算法的工程优势特别管用。
1.2 二次样条曲线和三次样条曲线的选型差异
B样条曲线按次数分,工程里最常用的是二次和三次。
二次样条曲线的基函数是二次多项式,曲线整体只保证一阶导数连续,也就是C1连续。它的优势是计算量小,在性能紧张的MCU上表现友好,而且对于只需要视觉平滑的应用完全够用。缺点是曲率变化不够柔和,在要求运动控制高平滑性的场景下,加速度会有跳变。
三次样条曲线的基函数是三次多项式,曲线保证二阶导数连续,也就是C2连续。这意味着不仅位置、速度连续,加速度也连续,用在运动控制里更稳,曲线本身的形态也更柔顺。代价是每个点计算需要遍历的基函数更多,代码逻辑也稍微复杂一点。
我的选型经验是:纯显示用途(绘图、UI笔迹)选二次就够了,涉及运动规划(CNC、机器人、伺服控制)直接上三次,别纠结。三次多出来的那点计算量,在现在的MCU上基本可以忽略。
提示:标题里有“二次样条曲线”和“三次样条曲线”两个方向,代码实现上建议两者都保留,通过一个阶数参数切换,因为实际项目里你永远不知道下一个需求要哪个。
2. C语言数据结构与算法核心实现
2.1 控制点与节点向量的存储设计
B样条曲线的数学表达是C(u) = Σ Ni,p(u) * Pi,即所有控制点Pi乘以对应的基函数Ni,p(u)求和。要落地成C代码,首先得设计一套清晰的数据结构。
我的做法是定义一个结构体,把曲线信息集中管理:
typedef struct { int dim; /* 空间维度:2D或3D */ int num_control; /* 控制点个数 */ int degree; /* 曲线次数:2或3 */ double *controls; /* 控制点坐标数组,连续存储 */ double *knots; /* 节点向量 */ } BSplineCurve;这里有两个关键设计点。第一,控制点用一维数组连续存储,而不是用二维数组或指针数组,这样在做指针遍历、传给SIMD指令时都方便,也能减少内存碎片。第二,节点向量长度不是拍脑袋定的,它等于控制点个数加次数加1,这个公式必须记死。比如5个控制点、3次样条,节点向量长度就是9。
节点向量是B样条里最难懂的概念,但它在代码里不过是一个单调不减的double数组。它把参数区间切分成若干小段,每个分段对应曲线的一段,分段的数量和位置直接决定曲线在哪些地方“拐弯”。我建议初学时把节点向量打印出来,逐个值看,再看曲线上对应的分段,比单纯背定义有效得多。
2.2 基函数核心:Cox-de Boor递推实现
B样条的所有魔法都集中在基函数的计算上。业界标准算法是Cox-de Boor递推公式,它把p次基函数拆成两个p-1次基函数的线性组合,层层递归下去,直到0次基函数——0次的情况很简单,参数落在哪个区间就是1,否则是0。
下面这段代码是我项目里实际在用的,做了几个必要的安全检查:
/* 计算第i个p次B样条基函数在参数u处的值 */ double basisF(int i, int p, double u, double *knots) { if (p == 0) { /* 0次基函数:u落在左闭右开区间内返回1 */ return (u >= knots[i] && u < knots[i + 1]) ? 1.0 : 0.0; } double left = 0.0, right = 0.0; double denom1 = knots[i + p] - knots[i]; double denom2 = knots[i + p + 1] - knots[i + 1]; /* 分母为零时该段贡献记0,避免浮点除零 */ if (denom1 > 1e-12) { left = (u - knots[i]) / denom1 * basisF(i, p - 1, u, knots); } if (denom2 > 1e-12) { right = (knots[i + p + 1] - u) / denom2 * basisF(i + 1, p - 1, u, knots); } return left + right; }这里最容易出错的点就是“左闭右开”区间。很多新手在0次基函数里用了u >= knots[i] && u <= knots[i+1],导致相邻分段在节点处重复计算,细小的数值误差就出现了。我测过两台不同的开发板,这种边界问题在浮点运算下确实会造成肉眼可见的曲线端点毛刺,所以判断区间时一定要用左闭右开。
用纯递归实现基函数在“能跑”层面没有问题,但性能上要留个心眼。每个参数u都要重算所有相关基函数,递归层数越深重复计算越多。实际项目里我建议升级成自底向上的表格法:开一个(p+1)*(p+1)的数组,从0次开始逐层向上填表。这一点如果读者用的是C语言版本,值得优先优化。
2.3 曲线点计算:从参数u到位坐标
有了基函数,计算曲线点就成了遍历控制点并加权求和。但这里有个性能关键点:根据局部支撑性,当参数u落在某个节点区间里时,只有附近p+1个基函数非零,其他基函数严格为0。所以完整遍历所有控制点是浪费——正确做法是先找到u在节点向量中的位置,再只遍历那p+1个控制点。
/* 计算曲线上参数u对应的坐标点 */ void evalBSpline(BSplineCurve *c, double u, double *out) { int i, k; int p = c->degree; int n = c->num_control; /* 夹紧u到有效范围 */ if (u <= c->knots[p]) u = c->knots[p]; if (u >= c->knots[n]) u = c->knots[n] - 1e-10; /* 二分查找u所在的节点区间,返回区间序号k */ k = findSpan(n, p, u, c->knots); /* 清零输出坐标 */ for (i = 0; i < c->dim; i++) out[i] = 0.0; /* 只有下标k-p到k的控制点参与计算 */ for (i = k - p; i <= k; i++) { if (i < 0 || i >= n) continue; double N = basisF(i, p, u, c->knots); for (int d = 0; d < c->dim; d++) { out[d] += c->controls[i * c->dim + d] * N; } } }这里findSpan是标准的二分查找函数,返回参数u在节点向量中的区间序号。实现上有个细节:当u恰好等于最大节点值的时候,要先用knots[n] - 1e-10把它稍微拉回来一点,否则查找函数会越界。这个“-1e-10”的小动作,我调试过程中排查了好久才发现。
3. 曲线拟合与平滑实操:两条经验路线
3.1 拟合:用B样条逼近一批离散数据点
B样条曲线有两种用法。第一种是“给控制点,算曲线”,这是正算;第二种是“给数据点,反推控制点”,这就是曲线拟合——也是工程里更常见的需求。
拟合的思路本质上是最小二乘:先给每个数据点分配一个参数u(这一步叫参数化),然后建立线性方程组,让曲线在参数处的值尽量逼近数据点。
参数化这块,常见的弦长参数法实测效果比均匀参数法好很多。具体做法是累积相邻数据点之间的距离,除以总长度,得到每个点在0到1之间的占比,这一步可以先写出代码。
/* 累积弦长参数化 */ void chordLengthParam(double *data, int m, int dim, double *params) { double *dist = calloc(m, sizeof(double)); double total = 0.0; for (int i = 1; i < m; i++) { double sum = 0.0; for (int d = 0; d < dim; d++) { double diff = data[i*dim+d] - data[(i-1)*dim+d]; sum += diff * diff; } dist[i] = sqrt(sum); total += dist[i]; } params[0] = 0.0; for (int i = 1; i < m; i++) { params[i] = params[i-1] + dist[i] / total; /* 等比移动到非零区间,避免基函数边界退化 */ if (params[i] <= params[i-1]) params[i] = params[i-1] + 1e-6; } free(dist); }参数化完成后,把每个数据点代入基函数矩阵,得到一个超定方程组N * P = D,解得控制点P。因为B样条基函数是稀疏的,这个矩阵是带状的,解方程用高斯消元或者Cholesky分解都行。矩阵规模不大时用哪种都行,规模大了优先考虑带状结构。
3.2 平滑:压制数据噪声的曲线平滑策略
拟合和曲线平滑在数学上是一对双胞胎,区别在于目标函数。拟合只求“贴近数据点”,平滑则额外加入一个“曲线不要抖得太厉害”的惩罚项。最常见的做法是加上二阶导的能量项,目标函数变成:
F = ||N * P - D||² + λ * P' * K * P
其中K是二阶导的刚度矩阵,λ是平滑系数。λ越大,曲线越平滑但越偏离数据点;λ越小,曲线越贴近数据点但噪声越放大。这个λ在实际工程里没有固定标准,我自己的经验是先从0.01试起,看效果再逐步调整,噪声大的数据可以到0.1。
在C语言里实现时,上述目标函数最终归结为对线性方程组做一点装配修改,代码量不过三四十行。核心思路是:在常规系数矩阵的对角线上加上λ相关项,再解方程。对于只做数据平滑的场合,还可以直接用移动平均结合三次样条,先滤波再拟合,效果也很不错,但工业级项目还是推荐正则化方法,理论边界清晰,参数可解释。
3.3 触摸屏笔迹和传感器曲线:两个真实场景
我自己实测的两个场景可以参考。一个是对触摸屏采集的笔迹做平滑,原始坐标抖动大约1-2像素,用三次B样条配合λ=0.05的正则化,输出曲线肉眼完全看不到锯齿,而且笔画端点没有拖尾。另一个场景是ADC采集的温度曲线,采样有高频噪声,用二次样条加λ=0.02平滑后,曲线保留了温度变化的趋势,又滤掉了尖刺。
这两个场景有个共同点:数据点个数远远多于最终需要的控制点个数。这也正是B样条相比插值算法的优势——可以用远少于数据点的控制点逼近整条曲线,天然具有压缩特性和噪声抑制能力。
4. 完整代码骨架:拿来就能改的最小实现
4.1 一个可直接运行的拟合示例
下面给出一段完整可编译的最小C代码骨架,它读取一组二维点,用三次B样条拟合出控制点,然后均匀采样输出曲线点:
#include <stdio.h> #include <stdlib.h> #include <math.h> /* 省略basisF和findSpan实现,见上文 */ int main() { /* 示例数据点:一个圆弧的离散采样 */ double data[9][2] = { {0.0, 0.0}, {0.5, 0.2}, {1.0, 0.7}, {1.5, 1.3}, {2.0, 1.8}, {2.5, 2.0}, {3.0, 1.8}, {3.5, 1.3}, {4.0, 0.5} }; int m = 9, dim = 2, p = 3, n = 5; /* 5个控制点 */ BSplineCurve bc = {0}; bc.dim = dim; bc.degree = p; bc.num_control = n; bc.controls = calloc(n * dim, sizeof(double)); /* 构造clamped节点向量 */ bc.knots = malloc((n + p + 1) * sizeof(double)); for (int i = 0; i <= n + p; i++) { if (i <= p) bc.knots[i] = 0.0; else if (i >= n) bc.knots[i] = 1.0; else bc.knots[i] = (double)(i - p) / (n - p); } double *params = malloc(m * sizeof(double)); chordLengthParam((double*)data, m, dim, params); /* 拟合求控制点:装配N矩阵后解方程 */ /* 这里省略线性代数实现,读者可选用高斯消元 */ /* 采样输出 */ for (int s = 0; s <= 30; s++) { double u = (double)s / 30.0; double out[2]; evalBSpline(&bc, u, out); printf("%.4f %.4f\n", out[0], out[1]); } return 0; }上面代码里的拟合求解段我省略了,工程量主要集中在基函数矩阵装配上。矩阵的行数等于数据点数,列数等于控制点数,第i行第j列就是第j个基函数在第i个参数值处的值。把矩阵搭好,后面就是标准最小二乘解法,想深挖的可以看数值分析教材的莱文贝格-马夸特法章节,工程上我用的是QR分解。
4.2 Clamped节点向量与边界条件处理
新手最容易踩的坑是曲线两端“够不着”首尾数据点。这个问题99%的根源在于节点向量不是“夹紧”的。
所谓夹紧,就是节点向量首尾各重复p+1个相同的值。比如3次样条、5个控制点,节点向量应该长这样:
[0, 0, 0, 0, 0.25, 0.5, 0.75, 1, 1, 1, 1]首尾各4个重复值。这样做的好处是曲线严格经过第一个和最后一个控制点,对应到拟合场景就是曲线起点终点和第一条最后一条数据点严格对齐。如果不加重复,曲线是“悬浮”在控制点之间的,两端总会有一段偏差。
4.3 数值稳定性:C语言版本必做的两个防护
浮点运算在嵌入式环境里会暴露很多意想不到的问题。我代码里强制加了两条防护规则,也建议读者采纳。
第一条是分母保护。Cox-de Boor递推里出现分母knots[i+p] - knots[i]时,如果节点向量有重复值,分母可能变成0。必须在除之前检查,绝对值小于1e-12就直接把分子项置0,而不是跳过去不处理,否则影响其他项的求和。
第二条是参数越界钳制。所有对参数u的输入,先判断范围,超出就钳到有效区间内。别以为上位机调好的逻辑挪到单片机上没问题,浮点精度不同会带来微小偏移,一旦越界,数组下标直接负值,段错误防不胜防。
5. 常见问题与排查实录
5.1 曲线过起点和终点,但回头检查数组越界?
这是我把代码移植到不同平台时最常遇到的坑。表面症状是曲线两端正常,但运行一段时间后偶发数组越界。问题通常出在findSpan的二分查找逻辑:当参数u等于最后一个节点区间的右端点(也就是值为1.0的夹紧端点)时,查找结果是最后一个区间序号n-1,但后续遍历控制点下标是k-p到k,当k等于n-1时,k+1已经在节点向量有效范围内,可控制点数组却可能越界。解决方式只有一条:在evalBSpline入口处预先对u做knots[n] - 1e-10处理,或者把k钳制到[p, n-1]之间。
5.2 曲线在节点附近出现“鼓包”或尖角
曲线明明有C2连续保证,为什么还会出现尖角和鼓包?这里要区分数学连续和视觉连续。当某个区间内数据点特别密集、节点向量对应取值跨度特别小的时候,基函数在该区间的值会出现“尖峰”,尽管导数连续,但视觉上看起来就是鼓包。
我的排查步骤是:先打印节点向量的相邻差值,看看是不是出现小于0.001的跨度;再看参数化后的数据点分布,是不是在某个局部过度密集。解决办法是调整节点向量分布,让节点尽量和数据点密度匹配,或者减少控制点数量。我的实测经验是,当控制点数量超过数据点数量的一半时,鼓包出现的概率会明显上升。
5.3 浮点误差导致的曲线抖动
这个问题在双精度上不明显,一旦用float就会暴露。症状是曲线在节点附近有幅度极小的锯齿抖动,肉眼不易发现,但运动控制系统里会表现为微小的高频振荡。
原因在于基函数在节点边界处是通过两个大数相减得到的,float精度不够时尾数误差被放大。两套方案:一是关键计算改用double,输出时再转float,ARM Cortex-M系列支持double指令,性能损失不算大;二是在边界点处做特殊处理,当u与节点的距离小于1e-7时直接取上一个区间的端点值。我在STM32F4上实测下来,第一种方案更简单可靠。
写在最后的经验之谈
如果你正在学C语言,或者第一次接触样条曲线算法,我强烈建议先不急着上线最小二乘拟合,把基函数和曲线计算跑通,再深入拟合。拟合涉及矩阵求解,排错复杂度和基函数不在一个量级。反过来,如果你只需要曲线平滑功能,一个二次样条加λ正则化足够,别因为看到三次样条的“更平滑”就盲目上高次,计算量、数值稳定性都会多出不少麻烦。
调试B样条代码有一个很实用的技巧:把所有控制点和曲线点都打印成文本文件,用Excel或Gnuplot先绘出来,可视化问题基本一眼就能定位,比在调试器里盯着数组强多了。做曲线拟合这个方向,永远是数学清晰了代码才清晰,先写下公式、再写代码,真的能省去大量无效调试时间。
本文还有配套的精品资源,点击获取