简介这是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 的完整数学描述。NURBSNon-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(i1) u U(i2) N 1; else N 0; end return; end % 计算左项和右项的系数分母为零时该项视为0 left_num u - U(i1); left_den U(ip1) - U(i1); if left_den ~ 0 left left_num / left_den * bspline_basis(i, p-1, u, U); else left 0; end right_num U(ip2) - u; right_den U(ip2) - U(i11); if right_den ~ 0 right right_num / right_den * bspline_basis(i1, p-1, u, U); else right 0; end N left right; end这段代码是 Cox-de Boor 递推的直接翻译需要注意 MATLAB 数组下标从 1 开始而节点向量在数学文献中通常从 0 开始所以代码里做了 U(i1) 的偏移。递归写法清晰但效率一般工具箱里的生产版本用的是迭代写法本质上结果一致。当分母为零时意味着节点重复度达到了 p1此时基函数在该区间没有支撑直接取零即可。2.2 工具箱目录结构与工程分工解压后你会看到 src 文件夹下分成 Nurbs curves 和 Nurbs Surfaces 两个子目录另外根目录有 AddPath.m 和 EAN_rapport.pdf 文档。这种组织方式对 MATLAB 工具箱来说是比较标准的曲线相关的脚本处理 1D 参数域问题曲面脚本处理 2D 张量积问题两者共用一套节点向量工具函数。功能模块对应目录/文件典型函数曲线构造与求值Nurbs curvesnrbmak, nrbeval, nrbplot曲面构造与求值Nurbs Surfacesnrbmak (2D版), nrbdeval节点操作src 公用函数nrbkntins, nrbdegelev路径设置AddPath.maddpath(genpath(pwd))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坐标二维曲线z0 1 0.5 1]; % 权重w第四行 knots [0 0 0 1 1 1]; % 二次曲线需要6个节点值两端各重复3次 curve nrbmak(coefs, knots); % 在u0到u1之间均匀取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 是控制点数3p 是次数2所以长度为 7 的节点向量里内部只有一个非零节点值 1两端重复 p13 次保证曲线经过首末控制点。3.2 节点插入与曲线细分节点插入是 NURBS 最常用的操作之一用于在不改变曲线形状的前提下增加控制点密度。工具调用方式如下% 在u0.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(u0.4处坐标误差: %e\n, norm(p_orig - p_ref));插入节点后控制点数量从 3 变成 4但在 u0.4 处求出的坐标误差应在 1e-15 量级这是浮点误差而非算法误差。nrbkntins的第二个参数可以是标量插入单个节点或向量批量插入多个节点批量插入比逐个插入效率更高因为函数内部一次性更新所有受影响的控制点。后续做曲线逼近时经常先用节点插入加密控制点网格再调整控制点位置去逼近目标形状。3.3 可视化常见陷阱nrbplot的第二个参数是采样密度取值太小会得到折线感很强的曲线一般 50 到 100 够用。如果发现曲线首尾不经过控制点先检查节点向量两端是否各重复了 p1 次——这是最容易被忽视的问题。还有一个坑是控制点坐标和权重混在一起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精确构造单位半球面半径r1 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, np2); V linspace(0, 1, mq2); % 组装基函数矩阵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 常用排错手段用这个工具箱最容易踩的坑集中在三类第一是节点向量不合法内部节点值重复次数超过 p1 会导致基函数奇异排查方法是打印节点向量肉眼检查第二是控制点数组维度不对二维曲面必须传 {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 法向一致性的唯一判据。本文还有配套的精品资源点击获取
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 的完整数学描述。NURBSNon-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(i1) u U(i2) N 1; else N 0; end return; end % 计算左项和右项的系数分母为零时该项视为0 left_num u - U(i1); left_den U(ip1) - U(i1); if left_den ~ 0 left left_num / left_den * bspline_basis(i, p-1, u, U); else left 0; end right_num U(ip2) - u; right_den U(ip2) - U(i11); if right_den ~ 0 right right_num / right_den * bspline_basis(i1, p-1, u, U); else right 0; end N left right; end这段代码是 Cox-de Boor 递推的直接翻译需要注意 MATLAB 数组下标从 1 开始而节点向量在数学文献中通常从 0 开始所以代码里做了 U(i1) 的偏移。递归写法清晰但效率一般工具箱里的生产版本用的是迭代写法本质上结果一致。当分母为零时意味着节点重复度达到了 p1此时基函数在该区间没有支撑直接取零即可。2.2 工具箱目录结构与工程分工解压后你会看到 src 文件夹下分成 Nurbs curves 和 Nurbs Surfaces 两个子目录另外根目录有 AddPath.m 和 EAN_rapport.pdf 文档。这种组织方式对 MATLAB 工具箱来说是比较标准的曲线相关的脚本处理 1D 参数域问题曲面脚本处理 2D 张量积问题两者共用一套节点向量工具函数。功能模块对应目录/文件典型函数曲线构造与求值Nurbs curvesnrbmak, nrbeval, nrbplot曲面构造与求值Nurbs Surfacesnrbmak (2D版), nrbdeval节点操作src 公用函数nrbkntins, nrbdegelev路径设置AddPath.maddpath(genpath(pwd))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坐标二维曲线z0 1 0.5 1]; % 权重w第四行 knots [0 0 0 1 1 1]; % 二次曲线需要6个节点值两端各重复3次 curve nrbmak(coefs, knots); % 在u0到u1之间均匀取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 是控制点数3p 是次数2所以长度为 7 的节点向量里内部只有一个非零节点值 1两端重复 p13 次保证曲线经过首末控制点。3.2 节点插入与曲线细分节点插入是 NURBS 最常用的操作之一用于在不改变曲线形状的前提下增加控制点密度。工具调用方式如下% 在u0.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(u0.4处坐标误差: %e\n, norm(p_orig - p_ref));插入节点后控制点数量从 3 变成 4但在 u0.4 处求出的坐标误差应在 1e-15 量级这是浮点误差而非算法误差。nrbkntins的第二个参数可以是标量插入单个节点或向量批量插入多个节点批量插入比逐个插入效率更高因为函数内部一次性更新所有受影响的控制点。后续做曲线逼近时经常先用节点插入加密控制点网格再调整控制点位置去逼近目标形状。3.3 可视化常见陷阱nrbplot的第二个参数是采样密度取值太小会得到折线感很强的曲线一般 50 到 100 够用。如果发现曲线首尾不经过控制点先检查节点向量两端是否各重复了 p1 次——这是最容易被忽视的问题。还有一个坑是控制点坐标和权重混在一起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精确构造单位半球面半径r1 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, np2); V linspace(0, 1, mq2); % 组装基函数矩阵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 常用排错手段用这个工具箱最容易踩的坑集中在三类第一是节点向量不合法内部节点值重复次数超过 p1 会导致基函数奇异排查方法是打印节点向量肉眼检查第二是控制点数组维度不对二维曲面必须传 {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 法向一致性的唯一判据。本文还有配套的精品资源点击获取
