简介这套为 MATLAB 用户打造的 NURBS非均匀有理 B样条曲线处理工具箱主要面向几何建模、CAD 设计与工程计算场景帮助使用者完成离散数据点的曲线拟合、插值以及控制点反算等核心任务。压缩包共包含 6 个 M 文件整体仅有 8KB囊括 findspan、bspkntins、bspeval、bspderiv 等典型子程序对应节点区间查找、节点插入、曲线求值与导数计算等关键步骤便于直接调用或二次封装。通过调整控制点位置能够构造平滑连续的曲线有效避免简单线性插值带来的尖角与不连续性。工具箱体量精简、函数划分明确无需额外依赖即可运行适合熟悉 MATLAB 的工程师与学生快速上手并集成至自有拟合或建模系统中。目前已有 482 人学习下载对于需要精确构建和操纵复杂几何形状的开发者来说这套代码提供了低成本、易扩展的入门与参考实现。1. NURBS工具箱的价值从离散点到自由曲线做几何建模和逆向工程的人迟早会撞上同一堵墙拿一堆测量点想拉出一条光滑且可编辑的曲线用多项式插值会震荡用折线又太糙。几年前我在做叶片截面拟合时对比过三次样条和NURBS两种方案最终彻底转向后者因为NURBS把“逼近精度”和“形状控制”解耦了——你可以先固定曲线阶数和节点向量再去反算控制点改曲线形状时只动控制点不牵动全局。这就是nurbs_toolbox2这类工具箱存在的意义它把B样条基函数、节点插入、升阶、求导、控制点反算这些底层操作封装成MATLAB函数让非数值分析出身的工程师也能直接上手。本文要讲的是这套工具箱中六个核心函数怎么协同完成曲线拟合其中包括findspan、basisfun、bspeval、bspderiv、bspkntins和bspdegelev以及它们在插值和逼近场景下的真实用法。适合正在做CAD数据转换、轨迹规划或者曲线拟合的读者。2. 基函数与节点定位理解B样条的底层数学2.1 从B样条到NURBS的递进关系NURBS的全称是Non-Uniform Rational B-Splines中文叫非均匀有理B样条。它比普通B样条多了一个“有理”修饰意思是每个控制点带一个权重这使得它能精确表示圆弧、椭圆等圆锥曲线。但要把NURBS用起来第一步不是理解权重而是理解基函数。基函数定义在节点向量上输入参数u输出一组数量等于控制点数的数值这些数值是曲线形状的“系数组合”。在nurbs_toolbox2里basisfun.m就是用来计算基函数值的。它的输入通常是节点向量、当前参数位置u和阶数p输出是该u点上的非零基函数值。节点向量的分布方式决定了B样条曲线是否均匀均匀节点向量让曲线在参数空间均匀铺开但控制点分布会扭曲非均匀节点向量则允许在曲线曲率大的地方加密节点提高局部控制能力。理解这点很重要因为后面的拟合质量调优本质上是在调节点向量。2.2 findspan参数u落在哪个区间findspan.m是所有求值函数的前置步骤。给定一个参数u它返回节点向量中满足knots(i) u knots(i1)的区间索引。别小看这个二分查找它决定了后续基函数计算的范围。在曲线求值、求导、样条插值中如果不先定位区间直接遍历全部节点计算基函数当控制点上千时性能会急剧下降。% 在MATLAB中调用findspan的常见方式 n length(net) - p - 1; % 控制点数量减1 s findspan(n, p, u, knots); % s是返回的节点区间索引满足 k_{s} u k_{s1}这里的net是控制点矩阵p是曲线阶数。findspan内部采用二分搜索复杂度O(log n)而全区间遍历是O(n)。当曲线控制点从几十增加到几百时两种方式的耗时差异会非常明显。我在做点云截面线拟合时一次要计算几千个参数点用findspan可以节省约40%的时间。2.3 基函数的递推计算与边界处理B样条基函数的标准算法是Cox-de Boor递推basisfun.m就是按这个递推实现的。递推从0阶基函数开始1阶就是0阶的组合以此类推。关键细节在于重复节点重复度等于阶数时基函数可能不连续和端点节点通常重复度为p1。工具箱代码里对边界情况做了裁剪这在实际使用中很关键——当u恰好等于最后一个节点时不做处理会下标越界。% basisfun.m核心逻辑示意 % 输入i是区间索引u是参数p是阶数knots是节点向量 % 输出N是p1个非零基函数值 N zeros(p1, 1); N(1) 1.0; for j 1:p saved 0; for r 0:j-1 temp N(r1); den1 knots(i1r) - knots(i1j-p-1); % 左分母 den2 knots(i2j) - knots(i2r-p); % 右分母 if den1 ~ 0 N(r1) temp * (u - knots(i1j-p)) / den1; else N(r1) 0; end if den2 ~ 0 N(r1) N(r1) saved * (knots(i2j) - u) / den2; end saved temp; end end这段代码体现了递推的核心思想每一层都用上一层的基函数线性组合。分母为0意味着节点重复对应基函数为0必须显式处理否则会产生NaN。p阶曲线需要p1个非零基函数所以N的长度是p1。递推完成后这些值会乘以控制点并求和得到曲线上的点。3. 曲线求值与求导拆开bspeval和bspderiv3.1 bspeval曲线点的并行计算bspeval.m是NURBS曲线求值的主入口。它的输入是控制点矩阵net每行一个坐标或一个权重、节点向量knots、阶数p和一组参数值uu列向量输出是曲线上对应的点坐标P。这个函数内部先对每个参数点调用findspan得到区间再调用basisfun计算基函数最后做加权求和。% 给定控制点、节点向量、阶数求曲线上的点 ctrlpts [0 0; 1 2; 3 1; 5 4]; % 4个控制点二维 knots [0 0 0 0.5 1 1 1]; % 3次B样条3重端点 p 3; u linspace(0, 1, 20); % 取20个参数点 P bspeval(p, ctrlpts, knots, u); plot(P(:,1), P(:,2), b-, ctrlpts(:,1), ctrlpts(:,2), ro--);这里bspeval的第四个参数u必须按升序排列否则部分索引定位会出错。如果参数点中包含0或1工具箱内部能正确处理两端节点。返回的P矩阵行数与u相同列数等于控制点维度。我在实际项目中通常先把目标曲线离散成100到500个点用于后续误差检查。3.2 bspderiv导数曲线与切向量bspderiv.m用于计算B样条曲线的k阶导数。它返回的是控制点导函数曲线的控制点dc和节点向量dk。对于NURBS求导数比B样条复杂因为权重参与分子分母的微分。但bspderiv只处理非有理B样条对于NURBS工具箱通常先将有理权重连同坐标当作高维有理曲线用齐次坐标求导后再除以权重分量。% 计算一阶导数的控制点 [dc, dk] bspderiv(p, ctrlpts, knots); % dc是一个(p1)行的元胞数组第2个元素对应一阶导数控制点 % dk是节点向量比原节点向量缩短了p个节点导数控制点的数量少于原控制点原因是一个p阶曲线的导数是一个p-1阶曲线。dc的每一行对应一阶导数的一个控制点行数等于原控制点数减p。如果你想在曲线上某一点求切向量必须把导数曲线再求值一次。这里常见的坑是直接对dc调用bspeval时节点向量要用dk而不是原knots否则区间索引会错。3.3 求值顺序对数值稳定性的影响理论上B样条求值是稳定的但实际浮点计算中如果控制点坐标绝对值很大比如超过1e6而节点间距离很小基函数递推中分母会接近0导致数值误差放大。我的处理方法是先对控制点做归一化变换再求值最后变换回原坐标。bspeval内部不会做归一化所以大坐标场景下需要外部预处理。函数输入输出主要用途findspann, p, u, knots区间索引定位参数点所在节点区间basisfuni, u, p, knots基函数值向量计算B样条混合权重bspevalp, ctrlpts, knots, u曲线点坐标批量求值bspderivp, ctrlpts, knots导数控制点节点求切向量、曲率bspkntinsp, u, knots, ctrlpts, newknts新控制点新节点局部加密节点bspdegelevp, t, ctrlpts, knots, times升阶后控制点节点提高曲线的阶数这张表是我根据工具箱文件列表整理的对应关系其中u在findspan中是标量在bspeval中是列向量。理解每个函数的边界条件是避免写出小重构代码的关键。4. 节点插入与升阶曲线编辑的两把手术刀4.1 bspkntins不改变形状地加密控制bspkntins.m负责向节点向量中插入新节点同时重新计算控制点坐标使曲线几何形状保持不变。这个操作在工程中非常重要你需要增加控制点的自由度来拟合更细的局部特征但又不想改变当前曲线形态。插入一个节点对应控制点数量增加1原控制点会被局部更新。% 在节点值0.3处插入一次如果0.3不在节点向量中 newknots 0.3; [newCtrl, newKnots] bspkntins(p, 0, ctrlpts, knots, newknots); % 参数0表示插入次数实际工具箱按unique(newknots)处理重复插入节点需要注意的是newknots必须位于节点向量有效范围内而且不能超过最大重复度p1。如果你插入一个已存在的节点它的重复度会增加这会降低曲线在该位置的连续阶。比如插入一个重复度p的节点曲线在该点变成C0连续出现尖角。我有时候会用这个特性刻意制造折线角点比重新建模快得多。4.2 bspdegelev升阶的几何意义bspdegelev.m将B样条曲线的阶数升高比如从3次升到5次。升阶不会改变曲线形状但增加了控制点数量和自由度并且让曲线在新的阶数下拥有更多的形状调节空间。在拟合时如果低阶曲线无法满足精度升阶是首选而不是直接删除重拟合数据。times 2; % 升阶次数 [newCtrl, newKnots] bspdegelev(p, times, ctrlpts, knots);升阶后控制点的分布会更“平坦”因为原曲线在更高维空间里控制点有了更多余量。这里要注意升阶次数过多会导致控制点密集、数值稳定性下降。普通三次曲线升到五次通常够了再高没有实际意义。4.3 局部修改的完整流程实际编辑曲线时我经常组合使用节点插入和升阶。比如先定位一个高误差区域插入几个节点增加局部控制然后升阶提升整体平滑度。这个流程能应对大多数CAD曲线修改需求。值得注意的是bspkntins在内部也会调用basisfun来计算相关控制点所以底层函数的效率直接决定了编辑交互的流畅度。5. 从数据点到控制点拟合与反算的完整流程5.1 参数化数据点如何映射到u拟合的第一步是给每个数据点分配一个参数值u。最常用的是弦长参数化chord length它按照数据点间距的比例分配参数这样曲线沿数据点分布的“速度”更均匀。另一种是向心参数化它对数据点密集的局部更宽容。工具箱没有显式提供参数化函数但常见做法是% 均匀参数化 n length(x); u linspace(0, 1, n); % 弦长参数化 d sqrt(diff(x).^2 diff(y).^2); cumd cumsum([0; d]); u cumd / cumd(end);均匀参数化适用于数据点等距弦长参数化适用于不均匀分布。我用弦长参数化拟合一般曲线时平均误差能下降15%~30%因为曲线弧长和参数值更匹配。5.2 插值让曲线精确穿过每个数据点插值要求数据点都在曲线上这等价于求解一个线性方程组对于每个参数点u_i曲线方程C(u_i) sum N_{j,p}(u_i) * P_j。方程数量等于数据点数未知数是控制点P_j。通常令方程数量等于控制点数量并采用重建基函数矩阵的方式求解。% nodepoints是数据点坐标n是数据点个数 % 设定控制点数量 数据点数量阶数取3 nc n; p 3; knots zeros(nc p 1, 1); % 采用clamped节点向量节点两端重复p1次 % 内部节点均匀分布或者用到ur改造成UKnots求解控制点前的关键步骤是构造基函数矩阵A(i,j) N_{j,p}(u_i)然后P A \ D其中D是数据点坐标矩阵。数值上直接解可能病态我通常加一个微小的正则化项lambda * I让矩阵可逆。工具箱里没有直接提供A的构造函数但可以用basisfun针对每个u_i和区间循环组装。5.3 使用工具箱完成NURBS拟合完整的NURBS拟合流程比B样条复杂之处在于还要求解权重。如果数据点是二维或三维坐标且不需要精确表示圆弧可以暂时固定权重为1先做B样条拟合。若需要NURBS则先做B样条拟合再通过奇异值分解迭代调整权重或者用齐次坐标把每个点(x,y,w)当成三维坐标拟合后除以w得到(x,y)。下面是一个标准的流程示例使用工具箱函数完成曲线拟合并绘制结果% 原始数据点示例 D [0 0; 0.5 1.2; 1 1.8; 1.5 1.5; 2 0.8; 2.5 1.0; 3 1.6]; % 1. 参数化 d sqrt(sum(diff(D).^2, 2)); u [0; cumsum(d)] / sum(d); % 2. 阶数和节点向量设计 p 3; nc length(D); % 控制点数量 数据点数量用于插值 knots augknt(linspace(0,1,nc-p1), p1); % MATLAB自带的augknt可生成clamped节点 % 3. 组装基函数矩阵 A zeros(nc, nc); for i 1:nc ui u(i); si findspan(nc-1, p, ui, knots); N basisfun(si, ui, p, knots); % 基函数N对应区间[si-p:si]的控制点索引 idx si-p : si; A(i, idx1) N; % 注意MATLAB索引从1开始 end % 4. 求解控制点 P A \ D; % 5. 用bspeval求曲线上的点 uu linspace(0, 1, 200); C bspeval(p, P, knots, uu);逻辑说明第2步用augknt生成clamped节点向量保证曲线起点和终点经过首末控制点。第3步中findspan返回区间索引但basisfun实际上需要知道当前参数点对应的基函数在哪个控制点上非零因此通过si-p:si取得控制点范围。A矩阵第i行只在这一段上有非零值其余为零这是B样条局部支撑性的直接体现。最后一步bspeval能够批量求出曲线上的点用于可视化或误差分析。5.4 出现振荡或误差过大的调整方向拟合结果不理想时先检查节点向量是否太密或太疏。节点过多曲线会试图跟随数据点的微小噪声导致振荡节点过少则整体误差大。常见做法是增加节点数量然后观察误差分布把新节点添加到误差峰值附近。bspkntins在这里派上用场——先拟合再在误差大的位置插入节点重新拟合如此迭代。另一个方向是降低阶数。高次曲线容易在数据点之间摆动三次一般足够。如果数据点本身有噪声需要在拟合目标中加入正则化项最小化sum ||C(u_i)-D_i||^2 lambda * sum ||P_j - 2P_{j-1} P_{j-2}||^2第二项是二阶差分平滑。工具箱没有提供正则化求解器但求解线性方程组的地方换成稀疏矩阵linsolve即可。6. 精进技巧用findspan做快速误差评估与动态拟合6.1 基于findspan的误差分布可视化多数人评估拟合曲线时只计算最大误差和平均误差。但最大误差所在区间往往就是需要插入节点的位置而平均误差不能反映局部细节。我常用的技巧是用findspan把参数区间分段统计每一段的最大误差然后画出误差随参数u变化的曲线。这样能直观看到哪些区段拟合不足。% 已有拟合曲线bspeval得到C数据点D和参数u已知 err sqrt(sum((C - D).^2, 2)); % 每个数据点处的误差 % 对参数区间分成50段统计每段的最大误差 seg linspace(0, 1, 50); maxErr zeros(length(seg)-1, 1); for i 1:length(seg)-1 sel u seg(i) u seg(i1); if any(sel) maxErr(i) max(err(sel)); end end bar(midSeg, maxErr); % 绘制柱状图纵轴为误差横轴为参数段这段代码的价值在于柱状图的峰值段对应控制点的薄弱区域。配合节点插入我可以迭代地在最大误差段的中点添加一个节点通常两三轮后误差能下降一个数量级。6.2 半参数化拟合引入惩罚权重当数据点本身有测量噪声时精确插值会导致曲线产生锯齿。我更倾向于做最小二乘拟合即控制点数量少于数据点数量。此时A矩阵是欠定的需要求解min ||A P - D||^2用左除P A \ D即可自动得到最小二乘解。但这个解容易出现过拟合。改进方式是在目标函数中加入权重矩阵W让某些重要点位比如装配基准点的误差权重加大。% 定义权重向量基准点权重为10其他为1 w ones(size(D,1), 1); w([1 end]) 10; % 首末点更严格 W diag(w); P (A * W * A) \ (A * W * D);这里如果A*W*A的条件数很大建议用pinv代替\或者使用lsqlin添加边界约束。在参数化曲线时控制点数可以设定为数据点数的50%~70%既能保持形状又避免过拟合。6.3 把局部修改做成交互式小工具如果你要频繁调整曲线可以把上面的逻辑封装成一个MATLAB脚本点击曲线上的任意位置程序自动找到最近的参数u用findspan定位节点区间提示“插入节点”还是“移动控制点”。移动控制点时只需修改P的某一行然后用bspeval重新求值整个过程在毫秒级。这是NURBS工具箱最吸引我的地方——它的底层函数简单透明适合在此基础上构建上层应用。实际经验里用这个套路处理三维轨迹规划时控制点的数量从200降到60误差反而更均匀了。原因是节点向量经过误差驱动迭代加密后控制点分配更合理。如果你拿到nurbs_toolbox2建议先自己写一个拟合脚本把bspkntins和bspdegelev的结果分别画出来对比一下控制多边形和曲线的关系。这样做一次等于手工推导了一遍B样条的局部支撑性和升阶几何意义。本文还有配套的精品资源点击获取