简介B样条曲线是计算机辅助设计、图像处理与有限元分析中常用的数学工具面向图形学学习者与工程开发人员这套MATLAB程序包覆盖了从基函数构造、控制点定义到曲线绘制与拟合的完整流程。用户只需传入控制点矩阵和参数范围即可生成所需曲线并通过示例脚本直观理解节点向量、局部支撑等关键概念解决曲线设计时灵活性与精确性难以兼顾的问题也适合计算几何课程与工业建模需求。压缩包共25个文件主体为23个.m源码文件涵盖B样条基函数计算、de Boor求值、曲线导数、近似与估计等模块另有README和license说明便于快速上手整体体积仅22KB。目前已有360人学习下载。随包提供多个示例与GUI交互工具可观察控制点调整对曲线形态的影响也能将核心函数迁移至CAD建模、数据拟合等场景大幅节省底层算法开发时间非常适合课程作业、毕业设计和科研项目参考。1. 为什么在 MATLAB 里从零写 B 样条而不是直接查库做轨迹规划、CAD 数据光顺或传感器点云拟合时我常遇到一类需求手里只有一批离散坐标点想要一条形状可控、局部修改不影响全局的光滑曲线并且每个采样点的坐标都能精确算出来。B 样条bspline就是为这个场景设计的它把曲线拆成控制点、节点矢量knot vector和基函数三部分比 Bezier 更适合表示长曲线也比滤波后的折线更贴近原始数据。如果你有 Curve Fitting Toolboxspmak、fnplt确实几行就能出图但遇到曲线点要参与机器人插补、导数要作为速度指令下发、控制点需要在线反算这类场景黑盒函数往往不够。这篇文章从节点矢量和基函数讲起给出完整的 MATLAB bspline 编程路径算曲线点、反算控制点、求导、用 matlab 画图核对最后附上验证和调参技巧。适合已经会基本 MATLAB 编程、想自己掌握 B 样条核心逻辑的开发者。2. 节点矢量和基函数bspline 编程前必须先写对的两个基础B 样条曲线的定义是一条求和式C(u) Σ N_{i,p}(u) P_i其中 P_i 是控制点N_{i,p}(u) 是 p 次基函数。所有编程难点几乎都集中在基函数和它依赖的节点矢量上。这一节先把这两个概念用 MATLAB 代码钉死后面所有曲线点计算都建立在它们之上所以这里写错后面全是 NaN 和畸变曲线。2.1 节点矢量怎么定clamped、unclamped 与重复节点节点矢量 U [u_0, u_1, ..., u_m] 是一个非递减序列。对 p 次曲线、n1 个控制点节点个数必须等于 np2。最常见的 clamped 形式让首尾各 p1 个节点相等曲线端点正好落在首末控制点上unclamped 形式让节点均匀铺开曲线不再经过端点多用于做曲线拼接或周期性样条。类型节点矢量例子p3, n5端点行为clamped[0 0 0 0 0.4 0.6 1 1 1 1]C(0)P0C(1)P5unclamped[0 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1 1.125]不经过首末控制点内部重复节点[0 0 0 0 0.5 0.5 1 1 1 1]在 u0.5 处连续性降到 C^0注意表里的第三行内部节点重复 r 次连续性降为 C^{p-r}。重复 p 次就会出现一个明显的角点重复 p1 次时基函数计算会遇到零分母这是编程时最容易踩的坑后面 2.2 节会给出处理办法。工程上 90% 的场景都用 clamped 节点矢量因为轨迹规划、曲线拟合都希望曲线从第一个控制点出发、在最后一个控制点收尾。2.2 Cox-de Boor 递推MATLAB 里算基函数的最小实现基函数从零次开始递推N_{i,0}(u) 在 u_i ≤ u u_{i1} 时取 1否则取 0。高一次的基函数由两个低次基函数加权得到N_{i,p}(u) (u-u_i)/(u_{ip}-u_i) · N_{i,p-1}(u) (u_{ip1}-u)/(u_{ip1}-u_{i1}) · N_{i1,p-1}(u)直接按这个公式写成递归函数在 MATLAB 里性能差而且容易栈溢出。工程上常用 NURBS Book 里的三角递推版本一次循环同时算出某个节点区间上 p1 个非零基函数function N bsplineBasis(i, u, p, U) % 计算节点区间 [U(i1), U(i2)) 上的全部非零基函数 % i: 0-based 节点区间编号由 findSpan 返回 % 返回值 N(1)..N(p1) 对应 N_{i-p}(u) ... N_i(u) N zeros(1, p1); N(1) 1; left zeros(1, p1); right zeros(1, p1); for j 1:p left(j1) u - U(i2-j); % u 到左边节点的距离 right(j1) U(ij1) - u; % 右边节点到 u 的距离 saved 0; for r 0:j-1 temp N(r1) / (right(r2) left(j-r1)); N(r1) saved right(r2) * temp; saved left(j-r1) * temp; end N(j1) saved; end end这段代码的核心是saved变量每一层 j 的循环里前一个低次基函数的右半部分会作为下一个基函数的左半部分被累积避免重复计算。left和right数组存的是参数 u 到两侧节点的距离分母right(r2) left(j-r1)本质上就是相邻节点区间的长度。参数含义注意i节点区间编号0-based必须由 findSpan 返回不要手动猜u参数值clamped 下一般取 [0,1]p曲线次数p3 表示三次样条U节点矢量必须非递减否则求值结果乱跳提示当 u 恰好落在重复节点上、且两侧节点值相同时right left会等于 0MATLAB 会给出 Inf 或 NaN。后面所有调用点都需要先做端点判断这是 bspline 编程里最常见的隐性 bug。2.3 findSpan 区间定位基函数求值的第一步算基函数之前必须先找到 u 落在哪个节点区间 [u_k, u_{k1}) 里。线性扫描复杂度 O(n)对频繁求值的场景太慢工程上一般用二分。MATLAB 里数组从 1 开始编号所以边界判断要整体平移一个下标function k findSpan(n, p, u, U) % n: 控制点个数减 1 % 返回值 k 满足 U(k1) u U(k2)k 是 0-based 区间编号 if u U(end) k n; return; end if u U(1) k p; return; end low p 1; high n 2; mid floor((low high) / 2); while (u U(mid1) || u U(mid2)) if u U(mid1) high mid; else low mid; end mid floor((low high) / 2); end k mid; end两个特殊分支很重要u 到达末端节点时直接返回 nu 在起点时返回 p。这样做的原因是clamped 节点矢量下 u0 或 u1 时曲线分别退化为第一个和最后一个控制点二分法在这种边界上会陷入死循环或返回错误区间。二分结束后拿到的 k 就是 2.2 节bsplineBasis需要的参数 i。基函数还有一个必须满足的性质叫单位分解性任意参数 u 处所有非零基函数之和恒等于 1。这个性质在 6.2 节会作为代码自检项。现在把findSpan和bsplineBasis封装好就可以开始算曲线点了。3. 用 de Boor 算法计算 B 样条曲线点并用 MATLAB 画图有了基函数曲线点有两种算法直接对基函数求和或者用 de Boor 递推。两者的理论结果完全一致但工程代码里我几乎不用直接求和原因在 3.1 节说清。这一章的产出是一个可以反复调用的bsplinePoint函数以及一张能看出控制多边形和曲线关系的图。3.1 直接基函数求和与 de Boor 递推怎么选直接求和就是把 2.2 节的基函数结果套进定义式代码很直观function C bsplinePointDirect(P, p, U, u) n size(P, 1) - 1; k findSpan(n, p, u, U); N bsplineBasis(k, u, p, U); C zeros(1, size(P, 2)); for i 0:p C C N(i1) * P(k-pi1, :); % 只有 p1 个非零基函数 end end注意循环里只累加 p1 项因为 B 样条基函数具有局部支撑性对落在区间 [u_k, u_{k1}) 内的参数 u只有 N_{k-p}(u) 到 N_k(u) 这 p1 个基函数非零。所以直接求和法其实也是 O(p^2)不算慢。但 de Boor 递推的优势在于它完全不依赖基函数的显式三角递推而是用控制点自身的线性插值来完成求值数值上对高阶曲线更稳定而且求导、曲线求交这类算法都建立在 de Boor 的结构上。既然复杂度相同我建议直接上 de Boor代码量并没有增加。3.2 计算 B 样条曲线点的 MATLAB 函数de Boor 算法的思路是在区间 [u_k, u_{k1}) 上先取 k-p 到 k 这 p1 个控制点然后做 p 轮线性插值每轮点数减一最后一轮剩下的那个点就是曲线点function C bsplinePoint(P, p, U, u) % 计算 B 样条曲线上参数 u 对应的点 % P: (n1)行 x dim列 的控制点矩阵 % p: 曲线次数 % U: 节点矢量 % u: 参数值 if u U(1) C P(1, :); return; end if u U(end) C P(end, :); return; end n size(P, 1) - 1; k findSpan(n, p, u, U); d P(k-p1 : k1, :); % 取参与计算的 p1 个控制点 for r 1:p for j p:-1:r alpha (u - U(jk-p1)) / (U(jk-r2) - U(jk-p1)); d(j1, :) (1-alpha) * d(j, :) alpha * d(j1, :); end end C d(p1, :); end内层循环里j必须从 p 递减到 r因为每层 r 都要用上一层 r-1 的结果如果正向覆盖会污染数据。alpha是插值权重分子是 u 到区间左端点的距离分母是区间长度这个区间会随着 r 增大而逐渐收缩最终收缩到 u 所在的节点区间。两个端点分支直接返回首末控制点绕开了 2.2 节提到的重复节点除零问题。参数含义边界条件Pn1 行 x dim 列控制点矩阵行数必须 ≥ p1p曲线次数p3 是工程默认值U节点矢量长度 np2元素非递减u参数值端点由函数内部特判3.3 用 MATLAB 画图核对控制多边形与曲线算曲线点通常要批量采样。用 linspace 在参数域上生成一串均匀参数逐点调用bsplinePoint然后一次性 plot 出来。控制多边形用带圆圈的折线画出曲线用实线画出方便对比曲线是否被控制点骨架合理约束% 控制点一个略微弯曲的四边形 P [0 0; 1 2; 3 2; 4 0]; p 3; U [0 0 0 0 1 1 1 1]; % 4 个控制点三次 clamped uk linspace(0, 1, 201); C zeros(length(uk), 2); for t 1:length(uk) C(t, :) bsplinePoint(P, p, U, uk(t)); end figure; plot(P(:,1), P(:,2), o-, LineWidth, 1.2); hold on; plot(C(:,1), C(:,2), r-, LineWidth, 1.8); xlabel(x); ylabel(y); legend(控制多边形, B样条曲线点, Location, best); grid on; axis equal;这个例子只有 4 个控制点、节点矢量两端各重复 4 次曲线会退化成一条经过首末控制点的三次 Bezier是验证代码正确性的最小用例。采样密度方面屏幕上显示 101201 个点足够但如果曲线点要用于数控插补或机器人路径应该按弦高误差或速度规划来定采样步长而不是固定点数。肉眼检查时主要看两点曲线是否落在控制多边形凸包内以及移动一个控制点后曲线是否只在附近几个区间发生变化。4. 从离散点反算 B 样条控制点拟合与插值的 MATLAB 实现前面是正向求值已知控制点算曲线点。实际工程里更常见的是反向问题只有一批测量点或关键路径点想要一条 B 样条曲线穿过或贴近它们。这个问题的完整链路是数据参数化、构造节点矢量、组线性方程组反解控制点。三条链路各自有坑这一节按顺序讲。4.1 数据点参数化均匀、弦长与向心三选一反算控制点之前必须给每个数据点 Q_j 分配一个参数值 u_j否则方程组无从建立。最简单的均匀参数化让 u_j j/m但数据点间距不均匀时曲线会在密集区抖动。工程默认是弦长参数化把每段距离累加后归一化如果点列曲率变化剧烈用向心参数化更稳它对弦长开根号抑制过冲function u paramPoints(Q, alpha) % Q: (m1)行 x dim列 的数据点 % alpha1: 弦长参数化; alpha0.5: 向心参数化 d diff(Q); d sqrt(sum(d.^2, 2)); % 逐段弦长 if alpha ~ 1 d d .^ alpha; end u cumsum([0; d]); u u / u(end); end参数化公式适用场景均匀u_j j/m数据点等间距采集弦长弦长累加后归一化一般 CAD 数据默认选择向心弦长的 0.5 次幂累加曲率变化大的点列数据点重合时弦长会为 0此时应向心参数化或先剔除重复点。参数化直接影响最终曲线形状同样的数据点换一种参数化结果能差出几个像素到几个毫米在机器人轨迹里就是实打实的位置误差。4.2 最小二乘拟合组矩阵、解控制点拟合的目标是让曲线尽量贴近数据点通常控制点个数少于数据点个数。把每个数据点的参数 u_j 代入基函数得到系数矩阵 A其中 A(j,i) N_{i,p}(u_j)然后解最小二乘问题function P bsplineFit(Q, p, nCtrl, alpha) % 最小二乘拟合用 nCtrl 个控制点逼近数据点 Q m size(Q, 1) - 1; % 数据点个数减 1 n nCtrl - 1; % 控制点个数减 1 u paramPoints(Q, alpha); U knotAveraging(u, p, n); A zeros(size(Q, 1), nCtrl); for j 1:size(Q, 1) if u(j) U(1) A(j, 1) 1; % 起点处只有 N_0 非零 elseif u(j) U(end) A(j, end) 1; % 终点处只有 N_n 非零 else k findSpan(n, p, u(j), U); N bsplineBasis(k, u(j), p, U); A(j, k-p1 : k1) N; end end P (A*A) \ (A*Q); end其中knotAveraging是工程界标准的节点矢量构造方法内部节点取相邻 p 个参数值的平均保证每个基函数支撑区间内至少有一个数据点参数方程组良态function U knotAveraging(u, p, n) % 由参数值构造 clamped 节点矢量n 为控制点个数减 1 U zeros(1, n p 2); U(1:p1) u(1); U(end-p:end) u(end); for j 1:n-p s 0; for i j:jp-1 s s u(i1); end U(jp1) s / p; end end求解用的是正规方程 (AA)(AQ)。当控制点数量接近数据点数量时AA 的条件数会变大求解结果对噪声敏感。我的习惯是先看cond(A*A)超过 1e8 就减少控制点数量或改用等距节点。装了 matlab优化工具箱的话需要约束控制点范围或加平滑惩罚时可以把目标写成 min ||A P - Q||用lsqlin直接带 lb/ub 求解效果比手工加正则项更可控。4.3 插值反求控制点当数据点个数等于控制点个数当数据点个数等于控制点个数时A 是方阵问题从拟合退化为插值曲线精确穿过每个数据点。求解只有一行function P bsplineInterp(Q, p) % 插值控制点个数 数据点个数曲线精确穿过 Q m size(Q, 1) - 1; n m; u paramPoints(Q, 1); % 插值一般用弦长参数化 U knotAveraging(u, p, n); A zeros(size(Q, 1), n1); for j 1:size(Q, 1) % 组装 A 的代码与 bsplineFit 完全相同省略重复部分 end P A \ Q; endclamped 节点矢量天然保证第一条和最后一行系数矩阵是单位向量所以插值结果必定经过首末数据点不需要额外处理端点条件。但插值对噪声极其敏感数据点只要带上毫米级测量误差反算出的控制点就可能剧烈摆动。工程上我一般只在数据点本身可信、且数量不多时用插值比如关键路径点插值点云或传感器数据一律走 4.2 节的最小二乘并手动调 nCtrl 观察残差在贴合度和光滑度之间取平衡。5. B 样条曲线求导、闭合曲线与节点矢量调参曲线点算通之后很多应用还需要一阶导甚至二阶导机器人路径里一阶导对应速度二阶导对应加速度CAD 光顺里要检查曲率。B 样条求导不需要数值差分直接对控制点做一次线性变换得到导数控制点精度和速度都好得多。另外闭合曲线和节点调参也是高频需求这一章一次讲完。5.1 一阶导曲线用导数控制点算速度p 次 B 样条的一阶导是 p-1 次 B 样条导数控制点 Q_i 由原控制点差分得到Q_i p · (P_{i1} - P_i) / (u_{ip1} - u_{i1})i 0..n-1同时节点矢量要掐头去尾去掉第一个和最后一个节点。MATLAB 实现function dC bsplineDeriv(P, p, U, u) % 计算 B 样条曲线在参数 u 处的一阶导向量 n size(P, 1) - 1; Q zeros(n, size(P, 2)); for i 1:n denom U(ip1) - U(i1); % u_{ip1} - u_i Q(i, :) p * (P(i1, :) - P(i, :)) / denom; end U2 U(2:end-1); % 去掉首尾节点 dC bsplinePoint(Q, p-1, U2, u); end这里有个索引细节要特别说原曲线有 n1 个控制点导数曲线有 n 个控制点所以循环里 i 从 1 到 n对应 0-based 的 0 到 n-1。分母中 U(ip1) 是节点矢量里 u_{ip1} 的位置U(i1) 是 u_i 的位置差值是相邻非零区间的跨度。U2的长度自动满足 p-1 次曲线对节点数量的要求不需要手动算。用的时候在采样循环里同时调用两个函数for t 1:length(uk) C(t, :) bsplinePoint(P, p, U, uk(t)); dC(t, :) bsplineDeriv(P, p, U, uk(t)); % 速度向量 end speed vecnorm(dC, 2, 2); % 每个采样点的速率二阶导只要对 Q 和 U2 再套一次同样的公式注意 U2 要继续掐头去尾次数降到 p-2。对首末控制点重合的闭合轨迹一阶导验证尤其重要因为闭合点处的速度必须连续。5.2 闭合 B 样条周期扩展控制点的做法只要把首末控制点设成同一个点曲线并不会自动在接缝处光滑因为接缝两侧的基函数定义域没有打通。严格的做法是把前 p 个控制点周期性地接到尾部同时使用均匀周期节点矢量求值时把参数折叠回有效区间function C closedBSplinePoint(P, p, u) % P: n 行 x dim 列n p1曲线首尾自动相接 n size(P, 1); Pext [P(end-p1:end, :); P; P(1:p, :)]; % 周期扩展 U 0 : (size(Pext, 1) p); % 等距节点 u p mod(u - p, n - p); % 参数折回 [p, n) C bsplinePoint(Pext, p, U, u); endPext把末尾 p 个控制点搬到前面、开头 p 个搬到后面让接缝两侧的基函数都能拿到完整的支撑区间。节点矢量用等距整数序列这是周期样条的标配。参数折叠用了一个 mod 技巧u 落在 [p, n) 之外时先平移再取模保证接缝处左右求值结果完全一致一阶导也连续。画图时可以放心地让参数从 0 扫到 n曲线首尾自然闭合。这个函数我一般配合 5.1 节一起用检查接缝处速度方向是否突变。5.3 节点重复与次数选择局部性怎么被调出来B 样条相对 Bezier 最大的优势是局部性移动第 i 个控制点只有 [u_i, u_{ip1}) 这段参数区间内的曲线会变化。这个局部范围由节点矢量和次数共同决定。内部节点重复次数 r连续性曲线表现1C^{p-1}默认光滑p3 时有 C^2pC^0曲线在该处出现角点p1基函数退化曲线被分离成独立段实际调参时记住两个经验一是想在某个位置强行制造尖角与其移动控制点不如把该处的内部节点重复 p 次效果更可控二是次数不要盲目提高p3 的 C^2 连续对绝大多数机械应用已经足够p5 以上数值稳定性明显变差控制点的微小扰动会被放大求解线性方程组时尤其明显。控制点数量同样影响局部性控制点越密单个控制点影响的区间越短曲线越贴数据但也越容易把噪声带进来。提示调节点矢量时改完一定要重新检查 U 是否非递减重复节点被插到错误位置会让 findSpan 直接返回错误区间表现是曲线点跳变排查起来比写代码烦得多。6. 验证 B 样条代码和曲线结果的几个实用技巧写 B 样条代码最怕的是看起来对实际错。基函数索引偏一位、节点矢量长度不对、端点除零这类 bug 在图上很难一眼发现。我会在交付前跑三组验证数值差分检验导数、基函数不变量自检、批量求值的性能与导出验证。这三个技巧分别针对算法正确性、数据结构正确性和工程可用性。6.1 用有限差分验证导数曲线bsplineDeriv的公式只要索引写错就会出错而且图上不一定看得出来。最简单的验证是取一个远离节点的内部参数用中心差分和解析导数对比h 1e-6; u0 0.37; % 避开节点值选一个普通内部参数 num (bsplinePoint(P, p, U, u0h) - bsplinePoint(P, p, U, u0-h)) / (2*h); ana bsplineDeriv(P, p, U, u0); disp(norm(num - ana)); % 期望在 1e-7 量级以下误差在 1e-6 数量级说明实现正确如果差到 1e-2 以上优先查bsplineDeriv里 U 的下标偏移。注意 u0 不能取在节点值上否则中心差分跨过不连续点有限差分本身就失真。6.2 基函数归一性与端点插值自检基函数单位分解性是一个强检查任意参数处p1 个非零基函数之和必须恒为 1。用随机参数跑一百次任何一次偏差超过 1e-12 都说明bsplineBasis或findSpan有 bugfor t 1:100 u 0.001 0.998 * rand; k findSpan(n, p, u, U); N bsplineBasis(k, u, p, U); if abs(sum(N) - 1) 1e-12 error(基函数不满足单位分解性); end end同时验证 clamped 曲线的端点插值性质bsplinePoint 在 u0 时严格等于第一个控制点u1 时严格等于最后一个控制点。这两条过了基本可以确认索引体系是自洽的。注意 u 的取值范围要避开 0 和 1 的精确端点因为那两处走的是特判分支自检不到三角递推本身。6.3 批量采样提速与曲线点导出对几千个采样点逐点调用findSpan和三角递推MATLAB 循环会比较慢。常见做法是保持函数不变把外层循环换成 parfor或者先把所有参数一次性传入、在函数内部分段向量化。我一般先确认点数少于 1 万点直接循环超过 1 万点才向量化因为向量化的下标管理容易引入偏移 bug收益却有限。曲线点算完总要落到文件里和其他工具对接writematrix(C, curve_points.csv); % 供 CAD 或仿真软件读取 save(curve_points.mat, C, P, U); % 留档 MATLAB 变量导出前用isequal(C(1,:), P(1,:))这类断言做最后一次检查能挡住大多数索引错误。整套流程跑下来自制 bspline 代码的正确性和工程可用性就有了明确保障。最后再提醒一个操作层面的细节凡是改过节点矢量或控制点数量务必重新跑一遍 6.2 的自检再进画图环节否则一张看似正常的曲线图很容易掩盖底层的下标错误。本文还有配套的精品资源点击获取