MATLAB NURBS工具箱：曲线曲面建模与拟合实践

简介：这是MATLAB工具箱集锦压缩包，面向科研人员、工程师与学生，将杂散于各领域的实用工具箱汇总到一起，省去逐个寻找安装包的麻烦。压缩包含57个文件，m脚本和函数约20个、png示意图34张，另有pdf说明、README与license文件，整体仅3.8MB，下载和离线使用都很方便。内容覆盖Nurbs曲面建模、粒子群与遗传算法、医学图像处理、SVM模式识别、马尔可夫决策、机器人工具箱、地震数据处理、凸优化CVX、时间序列hctsa等众多方向，还包含数学建模、聚类分析、语音与音频处理、相机标定等配套脚本。每个工具箱大多附带示例图片与函数源码，便于对照运行、理解输入输出和二次改造。目前已有1489人学习或下载，尤其适合做课程设计、算法对比和项目预研，也适合需要跨领域MATLAB算法工具、希望快速搭建实验环境的MATLAB使用者。

1. 从曲线失真说起：为什么需要 NURBS 工具箱

做 CAD 数据交换或者逆向工程时，你一定遇到过这种场景：从 STEP/IGES 文件里读入一条样条曲线，用 MATLAB 打开一看，控制点对得上，但插值出来的形状和原设计差了十万八千里。原因多半是底层用三次多项式或普通 Bezier 做了近似，而不是用 NURBS 的完整数学描述。NURBS（Non-Uniform Rational B-Spline，非均匀有理 B 样条）能精确表示圆锥曲线、球面和自由曲面，是 CAD/CAM 领域的工业标准。这套 Nurbs-surface 工具箱提供了一套完整的 MATLAB 实现，覆盖曲线、曲面建模、节点插入、升阶与可视化，适合做几何算法验证、逆向工程预处理和机器人路径规划前的曲面建模。下面我从数学原理讲到代码复现，把节点向量、控制点权重和曲面求值这几个容易翻车的点逐一拆开。

2. NURBS 数学模型与工具箱的代码骨架

2.1 从 B 样条到 NURBS：权重的作用

NURBS 与普通 B 样条的唯一区别，就是每个控制点带了一个权重因子 w。一条 p 次 NURBS 曲线的表达式为：

C(u) = Σ Ni,p(u) · wi · Pi / Σ Ni,p(u) · wi

其中 Ni,p(u) 是定义在节点向量 U 上的 B 样条基函数。分子是加权控制点求和，分母做归一化，这保证了几何形状不会因为权重整体缩放而改变。当所有 wi = 1 时，NURBS 退化为标准 B 样条。工具箱里的曲线脚本正是基于这个递推公式实现的，没有用符号计算，全部是数值迭代，速度上能满足实时交互。

2.1.1 基函数递推的 MATLAB 实现

function N = bspline_basis(i, p, u, U) % B样条基函数递推，Cox-de Boor公式 % 输入: i - 控制点下标, p - 次数, u - 参数值, U - 节点向量 % 输出: 基函数值 N if p == 0 if u >= U(i+1) && u < U(i+2) N = 1; else N = 0; end return; end % 计算左项和右项的系数，分母为零时该项视为0 left_num = u - U(i+1); left_den = U(i+p+1) - U(i+1); if left_den ~= 0 left = left_num / left_den * bspline_basis(i, p-1, u, U); else left = 0; end right_num = U(i+p+2) - u; right_den = U(i+p+2) - U(i+1+1); if right_den ~= 0 right = right_num / right_den * bspline_basis(i+1, p-1, u, U); else right = 0; end N = left + right; end

这段代码是 Cox-de Boor 递推的直接翻译，需要注意 MATLAB 数组下标从 1 开始，而节点向量在数学文献中通常从 0 开始，所以代码里做了 U(i+1) 的偏移。递归写法清晰但效率一般，工具箱里的生产版本用的是迭代写法，本质上结果一致。当分母为零时，意味着节点重复度达到了 p+1，此时基函数在该区间没有支撑，直接取零即可。

2.2 工具箱目录结构与工程分工

解压后你会看到 src 文件夹下分成 Nurbs curves 和 Nurbs Surfaces 两个子目录，另外根目录有 AddPath.m 和 EAN_rapport.pdf 文档。这种组织方式对 MATLAB 工具箱来说是比较标准的：曲线相关的脚本处理 1D 参数域问题，曲面脚本处理 2D 张量积问题，两者共用一套节点向量工具函数。

2.2.1 AddPath.m 的执行逻辑

# 在MATLAB命令行执行，将工具箱全部子目录加入搜索路径 run('AddPath.m')

AddPath.m 内部一般就是 addpath(genpath(pwd)) 这一行，genpath 会递归获取当前目录下所有子目录，再用 addpath 加入 MATLAB 搜索路径。执行之后可以在 MATLAB 里输入which nrbmak验证路径是否生效，返回的文件路径包含了 Nurbs-surface 目录就说明加载成功。我一般会在 startup.m 里加一行这个脚本，这样每次启动 MATLAB 都自动加载，省去手动运行的步骤。

3. NURBS 曲线的实现：控制点、节点向量与可视化

3.1 构造一条 2D NURBS 曲线

工具箱的核心数据结构是nrbmak生成的 struct，包含 form（曲线/曲面标记）、degree、control points（齐次坐标形式）和 knots 四个字段。下面用一条抛物线来演示完整流程，抛物线是二次曲线，必须用 NURBS 才能精确表示，普通多项式插值做不到。

% 构造一条二次NURBS抛物线：y = x^2 在 [0, 1] 区间 % 控制点取(0,0), (0.5,0.5), (1,1)，权重取 [1, 0.5, 1] coefs = [0 0.5 1; % x坐标 0 0.5 1; % y坐标 0 0 0; % z坐标（二维曲线z=0） 1 0.5 1]; % 权重w，第四行 knots = [0 0 0 1 1 1]; % 二次曲线需要6个节点值，两端各重复3次 curve = nrbmak(coefs, knots); % 在u=0到u=1之间均匀取50个点求值 u = linspace(0, 1, 50); p = nrbeval(curve, u); % 可视化 nrbplot(curve, 50); hold on; plot(coefs(1,:), coefs(2,:), 'ro--'); % 绘制控制多边形 legend('NURBS曲线', '控制多边形');

代码里 coefs 是一个 4×n 矩阵，前三行是 xyz 坐标，第四行是权重。nrbeval返回的 p 是一个 3×m 矩阵，每一列对应参数 u 处的曲线点坐标。把权重从 0.5 改成 2 再运行一次，你会看到曲线明显被拉向中间那个控制点，这就是权重对曲线形状的直观影响。二次 NURBS 曲线的节点向量长度遵循公式：n + p + 2，其中 n 是控制点数（3），p 是次数（2），所以长度为 7 的节点向量里内部只有一个非零节点值 1，两端重复 p+1=3 次保证曲线经过首末控制点。

3.2 节点插入与曲线细分

节点插入是 NURBS 最常用的操作之一，用于在不改变曲线形状的前提下增加控制点密度。工具调用方式如下：

% 在u=0.4位置插入一个新节点 curve_refined = nrbkntins(curve, 0.4); % 比较插入前后的控制点数量 fprintf('原始控制点数: %d\n', size(curve.coefs, 2)); fprintf('插入后控制点数: %d\n', size(curve_refined.coefs, 2)); % 验证形状一致性 p_orig = nrbeval(curve, 0.4); p_ref = nrbeval(curve_refined, 0.4); fprintf('u=0.4处坐标误差: %e\n', norm(p_orig - p_ref));

插入节点后控制点数量从 3 变成 4，但在 u=0.4 处求出的坐标误差应在 1e-15 量级，这是浮点误差而非算法误差。nrbkntins的第二个参数可以是标量（插入单个节点）或向量（批量插入多个节点），批量插入比逐个插入效率更高，因为函数内部一次性更新所有受影响的控制点。后续做曲线逼近时经常先用节点插入加密控制点网格，再调整控制点位置去逼近目标形状。

3.3 可视化常见陷阱

nrbplot的第二个参数是采样密度，取值太小会得到折线感很强的曲线，一般 50 到 100 够用。如果发现曲线首尾不经过控制点，先检查节点向量两端是否各重复了 p+1 次——这是最容易被忽视的问题。还有一个坑是控制点坐标和权重混在一起，MATLAB 里把 coefs(4,:) 忘了赋值时默认为 0，导致分母为零，nrbeval 会给出 NaN。每次构造完数据先用isnan(sum(coefs))做一次检查，能省不少排查时间。

4. NURBS 曲面实战：张量积构造与形状控制

4.1 从曲线到曲面：张量积原理

NURBS 曲面是两条曲线的张量积，控制点变成二维网格 P(i,j)，节点向量变成两组 U 和 V，分别对应曲面两个参数方向。曲面方程：

S(u,v) = Σi Σj Ni,p(u) · Nj,q(v) · wi,j · Pi,j / Σi Σj Ni,p(u) · Nj,q(v) · wi,j

这意味着曲面的每一行控制点可以看作一条 u 方向的 NURBS 曲线，每一列则是 v 方向的曲线。工具箱的 nrbmak 支持传 4×n×m 的三维数组作为控制点，第一维为 4 表示齐次坐标，n 和 m 分别是两个方向的控制点数量。

4.1.1 构造一个半球面

% 用NURBS精确构造单位半球面（半径r=1） r = 1; % 控制点网格 3x3，中心控制点内缩以模拟球面曲率 % 这里用一个简化的5x5网格演示张量积结构 ctrl_pts = zeros(4, 5, 5); for j = 1:5 for i = 1:5 theta = (i-1) * pi / 4; % u方向角度范围 0~pi phi = (j-1) * pi / 4; % v方向角度范围 0~pi ctrl_pts(1,i,j) = r * sin(phi) * cos(theta); ctrl_pts(2,i,j) = r * sin(phi) * sin(theta); ctrl_pts(3,i,j) = r * cos(phi); ctrl_pts(4,i,j) = 1; % 权重初始为1 end end % 两个方向的节点向量，二次曲面，两端重复3次 knots_u = [0 0 0 0.25 0.5 0.75 1 1 1]; knots_v = [0 0 0 0.25 0.5 0.75 1 1 1]; srf = nrbmak(ctrl_pts, {knots_u, knots_v}); % 求值并绘制曲面 [u, v] = meshgrid(linspace(0, 1, 30)); pts = nrbeval(srf, {u, v}); surf(squeeze(pts(1,:,:)), squeeze(pts(2,:,:)), squeeze(pts(3,:,:))); axis equal;

这里的权重全部取 1，因此实际上得到的是一个双三次 B 样条曲面。要让控制点网格精确逼近真实球面，需要调整权重为 cos(θ)·cos(φ) 的形式，工具箱的文档里对圆锥和球面给出了示例权重，EAN_rapport.pdf 中也对这一部分有数学推导。如果你拿到的曲面形状明显偏离预期，先检查节点向量内部节点值是否单调递增，这是 MATLAB 不会自动校验的硬性条件。

4.2 曲面拟合：从散点到 NURBS 曲面

逆向工程里最常见的需求是把测量点云拟合成 NURBS 曲面。整体思路分两步：先固定节点向量，再按最小二乘反算控制点。下面给出一个最简实现：

% 假设有 m×n 个散点数据，对应参数坐标 uq、vq % 构造两个方向的B样条基函数矩阵，然后求解线性方程组 function [ctrl_pts, U, V] = fit_nurbs_surface(data_pts, uq, vq, p, q, n, m) % data_pts: 3×m×n 数据点, uq,vq: 采样参数坐标 % 生成均匀节点向量 U = linspace(0, 1, n+p+2); V = linspace(0, 1, m+q+2); % 组装基函数矩阵（u方向和v方向） Nu = zeros(length(uq(:)), n); for k = 1:length(uq(:)) for i = 1:n Nu(k, i) = bspline_basis(i-1, p, uq(k), U); end end % Nv 同理... % 用 Kronecker 积构造张量积矩阵，求最小二乘解 A = kron(Nv, Nu); P = A \ data_pts(:); % 解线性方程组 ctrl_pts = reshape(P, [3, n, m]); end

这段代码的核心是把张量积曲面的双线性系统拆成 Kronecker 积形式，一次求解得到所有控制点。控制点数量 n、m 选多少直接决定拟合效果：太少误差大，太多会过拟合。我一般从数据量的 1/4 开始尝试，逐步加密节点，观察拟合误差变化曲线。如果误差在某个节点数之后不再明显下降，说明已经到达合理范围。

4.2.1 参数化对拟合效果的影响

散点数据的参数坐标 uq、vq 分配方式直接影响拟合质量。均匀参数化最简单，但数据点分布不均时容易出现震荡；弦长参数化（按相邻点欧氏距离累计分配）更稳定，是这个工具箱场景下推荐的做法：

% 弦长参数化示例：将一条曲线的点坐标转换为参数值 function t = chordal_param(pts) n = size(pts, 2); d = sqrt(sum(diff(pts, 1, 2).^2, 1)); t = [0 cumsum(d)]; t = t / t(end); % 归一化到[0,1] end

弦长参数化的好处是让相邻点之间的参数间隔与几何距离成正比，避免曲率大的区域参数间隔过密导致过拟合。替换掉均匀参数化之后，同样节点数下拟合误差通常能降低 30% 到 50%。

4.3 常用排错手段

用这个工具箱最容易踩的坑集中在三类：第一是节点向量不合法，内部节点值重复次数超过 p+1 会导致基函数奇异，排查方法是打印节点向量肉眼检查；第二是控制点数组维度不对，二维曲面必须传 {U, V} 元胞数组而非两个独立矩阵，传错时nrbmak会报维度不匹配；第三是求值范围越界，nrbeval在参数值等于节点向量最后一个元素时有时返回 NaN，因为基函数递推在末端开区间和闭区间的处理逻辑不同，处理办法是把求值参数稍微内缩一点，比如 linspace(0, 1-eps, n)。

提示：几何建模领域对 NURBS 的求值稳定性要求很高，日常开发时如果只是做数据可视化，可以用nrbdeval替代nrbeval，它在密集采样点场景下做了缓存优化，速度提升明显。

5. 曲面质量分析与 STL 导出

前面拿到了拟合曲面，怎么验证它适合做后续的有限元分析或 3D 打印？这里给一个实用的曲面光滑性检查技巧：计算曲面在等参线网格上的法向量并绘制高斯曲率分布图。高斯曲率突变的地方就是几何异常区域，往往对应控制点权重突变或节点向量异常。

% 计算NURBS曲面的一阶和二阶导，进而求高斯曲率 nrb_deriv = nrbderiv(srf); % 返回曲面的一阶导和二阶导信息 [u, v] = meshgrid(linspace(0.01, 0.99, 50)); [~, dudv] = nrbdeval(srf, nrb_deriv, {u, v}); % dudv 包含 Su, Sv, Suu, Suv, Svv Su = squeeze(dudv{1}(1,:,:)); Sv = squeeze(dudv{1}(2,:,:)); E = sum(Su.^2, 3); F = sum(Su.*Sv, 3); G = sum(Sv.^2, 3); % 高斯曲率 K = (L*N - M^2) / (E*G - F^2)，需要进一步求 L,M,N surf(u, v, K); % 绘制曲率热力图

运行这段代码前先确认nrbderiv在当前工具箱版本里返回的是完整的导函数张量还是仅返回采样值，两种实现对应的索引方式差异很大。我拿到一个不熟悉的 NURBS 工具箱时，第一步永远是disp(nrbderiv(srf))看输出结构。

STL 导出方面，常见做法是先在 u、v 方向生成足够密的采样网格，然后用surf2stl函数写出二进制 STL 文件。网格越密面片越多，建议先粗采样确认法向正确再加密。一个经验值：曲率变化平缓的区域 30×30 网格足够，急转弯处单独加密到 60×60。导出后的 STL 文件在 MeshLab 或 Cura 里打开检查有无破面，如果出现法向翻转，回 MATLAB 检查surf2stl输出的三角面片顶点顺序是否为逆时针——这是 STL 法向一致性的唯一判据。

本文还有配套的精品资源，点击获取