简介基于Matlab实现自由曲面造型方法的源码包面向计算机、电子信息工程、数学等专业学生用于课程设计、期末大作业或毕业设计中的曲面造型部分可重点参考NURBS旋转曲面、Bezier曲面、Coons曲面、Ferguson桥接等经典算法。包内共43个文件含24个m源码文件、18个asv自动备份文件及1个说明txtm代码为主程序与函数实现asv可辅助追踪修改过程txt用于快速了解文件用途。压缩包仅25KB轻量精简便于下载与解压。目前已有810人学习下载。源码覆盖控制网格构建、旋转母线生成、样条曲面构造与绘制等环节并提供半圆球面等示例脚本方便对照运行需要一定Matlab基础适合作为算法实现与调试验证的参考资料。1. 自由曲面造型的Matlab实现先选对参数化路线拿到“基于Matlab实现自由曲面造型方法源码”这类压缩包先别急着找主脚本双击。自由曲面造型在 Matlab 里其实有两条完全不同的路线一条是参数曲面用 Bezier、B 样条或 NURBS把曲面写成“控制点 基函数 节点向量”的加权和另一条是网格变形直接编辑三角网格顶点或者用拉普拉斯坐标做局部形变。前者生成的曲面天然有参数域适合做精细 A 面、叶片型线和光学自由曲面后者更像雕刻灵活但连续性难控制导到 CAD 系统里还要重新曲面化。源码的价值不在某个炫酷 GUI而在把“控制点如何决定曲面”这个核心关系打开给你看。理解这一点之后不管是正向设计还是逆向重构你的工作量都会落在控制点上而不是反复改网格。这篇文章按正向设计、逆向拟合、交互调参、导出验证四条线展开代码尽量不依赖额外工具箱基于 Matlab 基础环境就能跑适合正在做曲面重构、需要对曲率连续性和数据交接负责的工程师。2. 用B样条把自由曲面造型的最小源码跑通自由曲面造型最容易被轻视的一步是选曲面数学形式。Bezier、B 样条、NURBS 在 Matlab 里都能写但工程上我一般默认用三次 B 样条起步。先把这一层的原理和参数设对后面换成 NURBS 只是多一组权因子的事骨架不用动。2.1 为什么在Matlab里用B样条而不是BezierBezier 曲面是 Bernstein 基函数的张量积。它最大的问题是全局支撑移动任意一个控制点整个曲面都会受影响。自由曲面造型进入交互调形阶段后这点非常难受你只想动翼型前缘结果后缘也跟着鼓了一块。B 样条用节点向量把定义域切成若干区间每个控制点只在一段参数区间上起作用这叫局部支撑。所以你拖一个控制点影响范围是有限的、可控的。另一个实际原因是次数可以固定为 3而控制点数量任意增加Bezier 想增加控制点就要升阶曲面阶次很快失控。至于 NURBS它是在 B 样条基础上给每个控制点加一个权因子可以精确表示圆弧、圆柱等解析曲面。如果你要做的自由曲面造型不涉及标准二次曲面权因子只会让调参更头疼我一般把 NURBS 留到需要与 CAD 特征对齐时再用。2.2 节点向量与次数源码里最常改的三个参数在写代码之前先把自由曲面造型里必调的参数说清楚。B 样条曲面的定义域是两个方向节点向量 U、V加上两个次数 p、q以及控制点网格 cp(i, j, :)。初学时最容易乱的是节点向量长度。控制点数为 n1、次数为 p 时节点向量长度必须是 np2而且首尾都重复 p1 次这叫 clamped也是参数化曲面的标准做法。参数含义初值建议p, qu、v 方向曲面的最高次数双向都取 3光顺度不够再取 5n1, m1两个方向控制点数量比次数至少大 1正向设计从 6×6 起步U, V 节点向量参数区间划分控制局部影响范围首尾重复 p1 次内部均匀即可cp控制点网格三维坐标初始按设计意图摆成网格状内部节点均布是正向设计的默认选择。逆向拟合时内部节点位置应该尽量和数据点的参数分布对齐否则会出现某一段控制点扎堆、另一段曲面太平的问题。这个区别在源码里通常体现为一行linspace和按弦长累积的两种写法。2.3 张量积曲面求值一个能跑的源码模板这里给一个不依赖任何工具箱的最小求值函数。Bbasis用 Cox-de Boor 递推计算单方向基函数bspSurfacePoint再把两个方向的基函数做张量积得到曲面上对应参数 (u, v) 的点。function P bspSurfacePoint(cp, U, V, p, q, u, v) % cp: (n1) x (m1) x 3 控制点网格 n size(cp, 1) - 1; m size(cp, 2) - 1; % u 方向所有非零基函数 Nu zeros(1, n 1); for i 0:n Nu(i 1) Bbasis(i, p, u, U); end % v 方向所有非零基函数 Nv zeros(1, m 1); for j 0:m Nv(j 1) Bbasis(j, q, v, V); end % 张量积加权和控制点 P zeros(1, 3); for i 0:n for j 0:m P P Nu(i 1) * Nv(j 1) * squeeze(cp(i 1, j 1, :)).; end end end function N Bbasis(i, p, u, U) % 0 次基函数落在区间内为 1端点和半开区间单独处理 if p 0 if abs(u - U(end)) 1e-12 i numel(U) - p - 2 N 1; elseif u U(i 1) u U(i 2) N 1; else N 0; end return; end d1 U(i p 1) - U(i 1); d2 U(i p 2) - U(i 2); w1 0; w2 0; if d1 0 w1 (u - U(i 1)) / d1 * Bbasis(i, p - 1, u, U); end if d2 0 w2 (U(i p 2) - u) / d2 * Bbasis(i 1, p - 1, u, U); end N w1 w2; end代码里的d1、d2是相邻节点区间长度出现 0 表示重复节点对应项直接取 0避免除以 0。递归写法容易读但控制点超过 30×30 或者要做实时拖拽时建议改成先findSpan再 de Boor 循环求值的版本基函数只算局部一段速度会快一个量级。下面用 6×7 的控制点网格生成一张曲面cp zeros(6, 7, 3); [xg, yg] meshgrid(1:7, 1:6); cp(:,:,1) xg; cp(:,:,2) yg; cp(:,:,3) 0.35 * sin(xg / 7 * pi) .* cos(yg / 6 * pi); p 3; q 3; n size(cp, 1) - 1; m size(cp, 2) - 1; U [zeros(1, p 1), (1:n-p)/(n-p1), ones(1, p 1)]; V [zeros(1, q 1), (1:m-q)/(m-q1), ones(1, q 1)]; [ug, vg] meshgrid(linspace(0, 1, 60)); X zeros(size(ug)); Y X; Z X; for k 1:numel(ug) pt bspSurfacePoint(cp, U, V, p, q, ug(k), vg(k)); X(k) pt(1); Y(k) pt(2); Z(k) pt(3); end surf(X, Y, Z, FaceAlpha, 0.75, EdgeColor, none); hold on; mesh(cp(:,:,1), cp(:,:,2), cp(:,:,3), FaceColor, none, EdgeColor, [0.3 0.3 0.3]); view(3);运行后会看到一张平滑曲面包在控制点网格里。到这里自由曲面造型的正向链路已经通了控制点决定曲面形状节点向量决定控制点的影响范围次数决定整张曲面的连续性。把cp中某个控制点的 z 值改掉只有附近一小块曲面会变化这就是局部支撑的实际表现。2.4 结果检查先看曲面再看控制网曲面算出来后别急着进下一步。先做两个检查一是用mesh把控制网叠在曲面上确认控制网没有交叉、没有明显不合理的尖点二是检查曲面是否在控制点凸包内。B 样条曲面的每个点都是控制点的凸组合如果发现曲面跑到控制网外面十有八九是节点向量写错了或者控制点网格的维度顺序和期望不一致。max(max(abs(Z))) % 看高度是否在预期范围 plot3(cp(:,:,1), cp(:,:,2), cp(:,:,3), k.-, MarkerSize, 12); xlabel(x); ylabel(y); zlabel(z);如果只关心形状不关心参数化是否正确很容易把维度顺序颠倒的 bug 带进后续代码。在自由曲面造型源码里行列顺序一错拟合出来的曲面就是麻花。3. 散点到自由曲面的拟合源码参数反算与平滑控制正向设计是从控制点得到曲面逆向重构则是把点云或扫描数据变成控制点。这是自由曲面造型里另一个高频需求。实现方法分为插值和逼近。插值要求控制点数量和数据点数量相等曲面严格穿过每个点逼近则让控制点数量远小于数据点数量用最小二乘把噪声磨掉。扫描数据通常带噪声我会优先选逼近。3.1 参数化散点为什么不能直接画曲面散点要拟合 B 样条曲面第一步不是算基函数而是给每个散点分配一个参数坐标 (u, v)。只有拿到参数坐标才能把数据点写入曲面方程的张量积结构。对于按行列组织的扫描网格一种做法是直接用行列索引归一化作为参数更稳的做法是按行和列分别做弦长累积参数化避免数据点疏密不均时曲面出现抖动。% 对每一行做弦长参数化 u zeros(size(xg)); for r 1:size(xg, 1) seg sqrt(sum(diff(squeeze([xg(r,:); yg(r,:); zg(r,:)]), 1, 2).^2, 1)); u(r, :) cumsum([0, seg]); end u u / max(u(:));这里的假设是数据点已经按行、列组织好了。如果拿到的是完全无序的点云那就必须先做平面投影、测地距离或近邻图参数化这一步偏差比后面的拟合算法更影响结果。3.2 最小二乘加正则化的Matlab拟合源码下面代码把散点拟合成一张三次 B 样条曲面。核心是构造基函数矩阵 A然后对 x、y、z 三个坐标分量分别解一个带平滑正则的最小二乘问题。function [cp, U, V, fitP] fitBsplineSurface(P, u, v, p, q, nc, mc, lambda) % P: 数据点N x 3 % u, v: 每个数据点的参数坐标N x 1 % nc, mc: 两个方向控制点数量 % lambda: 平滑正则系数 U [zeros(1, p 1), (1:nc-p)/(nc-p1), ones(1, p 1)]; V [zeros(1, q 1), (1:mc-q)/(mc-q1), ones(1, q 1)]; A zeros(size(P, 1), nc * mc); for k 1:size(P, 1) Nu zeros(1, nc); Nv zeros(1, mc); for i 0:nc-1 Nu(i 1) Bbasis(i, p, u(k), U); end for j 0:mc-1 Nv(j 1) Bbasis(j, q, v(k), V); end A(k, :) kron(Nv, Nu); % 张量积基函数按控制点列优先展开 end % 二阶差分平滑矩阵惩罚相邻控制点之间过大的二阶变化 Dx sparse(diff(eye(nc), 2)); Dy sparse(diff(eye(mc), 2)); S kron(speye(mc), Dx * Dx) kron(Dy * Dy, speye(nc)); cpMat zeros(nc * mc, 3); for d 1:3 cpMat(:, d) (A * A lambda * S) \ (A * P(:, d)); end cp reshape(cpMat, [nc, mc, 3]); fitP A * cpMat; endA的每一行对应一个数据点每一列对应一个控制点记录的是该数据点在当前参数位置对每个控制点的敏感程度。kron(Nv, Nu)是张量积曲面结构的线性化写法。S是二阶差分矩阵它让相邻控制点不能变化太剧烈lambda 越大曲面越光顺但残差也越大。调用方式N 24; [xg, yg] meshgrid(linspace(-2, 2, N)); zg peaks(N) 0.08 * randn(N); P [xg(:), yg(:), zg(:)]; u (xg(:) 2) / 4; v (yg(:) 2) / 4; [cp, U, V, fitP] fitBsplineSurface(P, u, v, 3, 3, 8, 9, 5e-3); rmse sqrt(mean(sum((P - fitP).^2, 2)));拟合结束后cp就是你要的自由曲面造型控制点fitP是数据点对应参数位置下的曲面采样值。RMSE 用来判断拟合精度但我一般不会只看 RMSE而是把P和fitP叠在一起看最大偏差出现在哪里。梯度大的区域如果偏差集中说明那个方向的控制点数量不够。3.3 控制点数量与lambda的取值套路控制点数量和 lambda 没有一个万能公式但可以按数据规模和噪声水平给一组有效的初值。数据干净时控制点数量可以接近数据点数量数据带噪声时控制点数量要明显小于数据点数量否则噪声会被当成形状特征拟合进去。数据噪声水平控制点规模lambda 初值说明干净无噪声min(数据点数/4, 20) 左右0 或 1e-6保持拟合精度轻微噪声数据点数/41e-3 到 1e-2以光顺优先明显噪声数据点数/61e-2 到 0.1防止曲面波浪需要过点控制点数 数据点数0插值模式噪声大时不适用lambda 过小曲面会在噪声点附近出现凹陷lambda 过大真实特征也会被磨平。一个实用的检查方法是把 lambda 从 0 开始按 10 倍递增同时看 RMSE 和 max 误差随 lambda 的变化曲线。拐点附近的 lambda 通常就是可用的折中值。如果还想压制最大误差把\换成正则参数回归或lsqlin就可以在解里加入更多线性约束这是 Matlab 优化工具箱在自由曲面造型里最常见的用法。4. 交互拖拽控制点自由曲面造型过程中的参数调优正向造型的最后一步是手感。每次改控制点都重跑整个脚本看不到曲率变化也不适合现场调形。我会在 figure 上绑三个事件ButtonDownFcn选点WindowButtonMotionFcn移动点WindowButtonUpFcn释放点。这是 Matlab 原生交互里最直接的一种实现方式。4.1 可交互拖拽控制点的源码骨架三维视图里用鼠标直接改 z 值不直观所以我通常先在俯视图或侧视图里拖 x、yz 保持不变。需要改高度时切换到侧视图再拖一次。下面的骨架实现了“选最近控制点 移动 x、y 重绘曲面”的完整回路。function setupDragger(cp, U, V, p, q) fig gcf; ax gca; setappdata(ax, cp, cp); setappdata(ax, U, U); setappdata(ax, V, V); setappdata(ax, p, p); setappdata(ax, q, q); hControl plot3(ax, cp(:,:,1), cp(:,:,2), cp(:,:,3), ko); set(hControl, ButtonDownFcn, startDrag); function startDrag(~, ~) pt get(ax, CurrentPoint); cpd getappdata(ax, cp); dist (cpd(:,:,1) - pt(1,1)).^2 (cpd(:,:,2) - pt(1,2)).^2; [~, idx] min(dist(:)); [r, c] ind2sub(size(cpd(:,:,1)), idx); setappdata(ax, dragIdx, [r, c]); set(fig, WindowButtonMotionFcn, movePoint); set(fig, WindowButtonUpFcn, endDrag); end function movePoint(~, ~) idx getappdata(ax, dragIdx); cpd getappdata(ax, cp); pt get(ax, CurrentPoint); cpd(idx(1), idx(2), 1:2) pt(1, 1:2); setappdata(ax, cp, cpd); redrawSurface(ax, cpd, getappdata(ax,U), getappdata(ax,V), ... getappdata(ax,p), getappdata(ax,q)); end function endDrag(~, ~) set(fig, WindowButtonMotionFcn, ); set(fig, WindowButtonUpFcn, ); end endredrawSurface就是把第 2 章的求值逻辑再跑一遍用新的控制点刷新surf对象的XData、YData、ZData。这里没有重新生成曲线所以鼠标拖动时响应很快。如果控制点很多每次重算 60×60 的曲面也够了再大就要缓存基函数只在被移动控制点影响的局部范围重新求值。4.2 边界连续与切矢约束自由曲面造型在单个曲面内部一般没问题麻烦出在边界。两个曲面要拼接时只把边界控制点对齐还不够还要保证跨界切矢方向一致。拖拽边界控制点时如果只动第一行第二行控制点不动G1 连续性会立刻被破坏。所以我在交互逻辑里会先记录下来拖拽前的边界间隔向量d移动边界控制点时下面一行跟着平移到新位置d cp(2, :, :) - cp(1, :, :); cp(1, :, :) newBoundary; cp(2, :, :) cp(1, :, :) d;这能保证边界法向和切平面方向不突然改变。若要更高阶的 C2 连续就要把第三行控制点也纳入约束自由度更少调形也受限通常曲率连续曲面才会用到。4.3 调参时最容易踩的三个坑现象原因处理曲面在端点处露馅、不贴近控制网节点向量没有端点多重重数首尾节点重复 p1 次拖动一个控制点曲面整体都变控制点数量等于次数1局部支撑失效增加控制点数量保持次数不变曲面出现明显波浪lambda 过小或控制点过多加大 lambda减少控制点交互调整的本质是优化问题只是优化变量变成了控制点坐标。如果你想让曲面自动避开一个障碍物或者满足某个面积约束把鼠标拖拽换成fmincon也是同一条路径目标函数是曲率和控制点移动量约束是障碍物距离。源码里一旦有了控制点这个中间层正向设计和优化就共用一套数据结构。5. 导出与验证让Matlab自由曲面造型源码能交接曲面造型做完只是第一步能导出、能重建、能被下一个环节消费才算真正收口。这里重点说两个导出目标STL 网格和 MAT 参数文件。5.1 把曲面转成STL并导出先用 2.3 的求值函数生成足够密的采样网格再把网格拆成三角面片。曲面越复杂采样网格要越密一般展示用 60×60做结构仿真用 150×150 以上。STL 只记录三角面片不保存参数化信息所以它只适合渲染和打印不适合继续做曲面编辑。function writeStlBinary(filename, X, Y, Z) [nr, nc] size(X); V [X(:), Y(:), Z(:)]; F zeros((nr - 1) * (nc - 1) * 2, 3); k 1; for i 1:nr-1 for j 1:nc-1 a i (j - 1) * nr; b a 1; c i j * nr; d c 1; F(k, :) [a b c]; k k 1; F(k, :) [b d c]; k k 1; end end fid fopen(filename, w); fwrite(fid, zeros(1, 80), uint8); fwrite(fid, size(F, 1), uint32); for k 1:size(F, 1) v1 V(F(k, 1), :); v2 V(F(k, 2), :); v3 V(F(k, 3), :); nrm cross(v2 - v1, v3 - v1); nrm nrm / norm(nrm); fwrite(fid, [nrm, v1, v2, v3], float32); fwrite(fid, 0, uint16); end fclose(fid); end调用writeStlBinary(freeform.stl, X, Y, Z)即可生成二进制 STL。注意三角形顶点顺序要按右手定则排列否则法向量指向会反查看器里会出现黑面或透光现象。5.2 用MAT文件保存可复现的曲面参数STL 和 MAT 是互补关系。STL 给下游消费MAT 给自己和团队留后路。保存时不要只存 X、Y、Z 采样网格那样等于把自由曲面降级成了一堆点。控制点cp、节点向量U、V、次数p、q才是曲面的完整定义任何时刻都可以用第 2 章的bspSurfacePoint重新采样。save(freeform_surface.mat, cp, U, V, p, q);下次打开时load这个文件再调一次bspSurfacePoint曲面不丢。逆向拟合的结果还可以把 lambda 和 RMSE 一起存进去这样别人拿到源码能知道当前曲面是逼近还是插值可信度到多少。5.3 一个避免失真的技巧按弧长重采样如果原始数据点间距不均匀拟合后曲面的控制点分布也会跟着失配密集区域曲面抖动稀疏区域曲面过渡得太快。与其硬调 lambda不如先把数据沿采样方向按弧长重采样arc cumsum([0; sqrt(sum(diff(P, 1, 1).^2, 2))]); arc arc / arc(end); Pnew interp1(arc, P, linspace(0, 1, size(P, 1)));重采样后再做fitBsplineSurface控制点分布会均匀得多。这个技巧对车身曲面、螺旋桨叶片这类长条数据尤其有效。5.4 用离屏渲染验证曲面完整性最后导出前用figure(Visible,off)做一次离屏渲染确认曲面没有 NaN、没有被控制点网络穿透再出图存档。fig figure(Visible, off); surf(X, Y, Z, EdgeColor, none, FaceColor, interp); light(Position, [1 1 1]); camlight; print(fig, surface_preview.png, -dpng, -r150); close(fig);如果这一步出现大面积黑斑优先检查法向量朝向和数据点参数化如果只是局部褶皱回到 lambda 和控制点密度上继续调。本文还有配套的精品资源点击获取
简介基于Matlab实现自由曲面造型方法的源码包面向计算机、电子信息工程、数学等专业学生用于课程设计、期末大作业或毕业设计中的曲面造型部分可重点参考NURBS旋转曲面、Bezier曲面、Coons曲面、Ferguson桥接等经典算法。包内共43个文件含24个m源码文件、18个asv自动备份文件及1个说明txtm代码为主程序与函数实现asv可辅助追踪修改过程txt用于快速了解文件用途。压缩包仅25KB轻量精简便于下载与解压。目前已有810人学习下载。源码覆盖控制网格构建、旋转母线生成、样条曲面构造与绘制等环节并提供半圆球面等示例脚本方便对照运行需要一定Matlab基础适合作为算法实现与调试验证的参考资料。1. 自由曲面造型的Matlab实现先选对参数化路线拿到“基于Matlab实现自由曲面造型方法源码”这类压缩包先别急着找主脚本双击。自由曲面造型在 Matlab 里其实有两条完全不同的路线一条是参数曲面用 Bezier、B 样条或 NURBS把曲面写成“控制点 基函数 节点向量”的加权和另一条是网格变形直接编辑三角网格顶点或者用拉普拉斯坐标做局部形变。前者生成的曲面天然有参数域适合做精细 A 面、叶片型线和光学自由曲面后者更像雕刻灵活但连续性难控制导到 CAD 系统里还要重新曲面化。源码的价值不在某个炫酷 GUI而在把“控制点如何决定曲面”这个核心关系打开给你看。理解这一点之后不管是正向设计还是逆向重构你的工作量都会落在控制点上而不是反复改网格。这篇文章按正向设计、逆向拟合、交互调参、导出验证四条线展开代码尽量不依赖额外工具箱基于 Matlab 基础环境就能跑适合正在做曲面重构、需要对曲率连续性和数据交接负责的工程师。2. 用B样条把自由曲面造型的最小源码跑通自由曲面造型最容易被轻视的一步是选曲面数学形式。Bezier、B 样条、NURBS 在 Matlab 里都能写但工程上我一般默认用三次 B 样条起步。先把这一层的原理和参数设对后面换成 NURBS 只是多一组权因子的事骨架不用动。2.1 为什么在Matlab里用B样条而不是BezierBezier 曲面是 Bernstein 基函数的张量积。它最大的问题是全局支撑移动任意一个控制点整个曲面都会受影响。自由曲面造型进入交互调形阶段后这点非常难受你只想动翼型前缘结果后缘也跟着鼓了一块。B 样条用节点向量把定义域切成若干区间每个控制点只在一段参数区间上起作用这叫局部支撑。所以你拖一个控制点影响范围是有限的、可控的。另一个实际原因是次数可以固定为 3而控制点数量任意增加Bezier 想增加控制点就要升阶曲面阶次很快失控。至于 NURBS它是在 B 样条基础上给每个控制点加一个权因子可以精确表示圆弧、圆柱等解析曲面。如果你要做的自由曲面造型不涉及标准二次曲面权因子只会让调参更头疼我一般把 NURBS 留到需要与 CAD 特征对齐时再用。2.2 节点向量与次数源码里最常改的三个参数在写代码之前先把自由曲面造型里必调的参数说清楚。B 样条曲面的定义域是两个方向节点向量 U、V加上两个次数 p、q以及控制点网格 cp(i, j, :)。初学时最容易乱的是节点向量长度。控制点数为 n1、次数为 p 时节点向量长度必须是 np2而且首尾都重复 p1 次这叫 clamped也是参数化曲面的标准做法。参数含义初值建议p, qu、v 方向曲面的最高次数双向都取 3光顺度不够再取 5n1, m1两个方向控制点数量比次数至少大 1正向设计从 6×6 起步U, V 节点向量参数区间划分控制局部影响范围首尾重复 p1 次内部均匀即可cp控制点网格三维坐标初始按设计意图摆成网格状内部节点均布是正向设计的默认选择。逆向拟合时内部节点位置应该尽量和数据点的参数分布对齐否则会出现某一段控制点扎堆、另一段曲面太平的问题。这个区别在源码里通常体现为一行linspace和按弦长累积的两种写法。2.3 张量积曲面求值一个能跑的源码模板这里给一个不依赖任何工具箱的最小求值函数。Bbasis用 Cox-de Boor 递推计算单方向基函数bspSurfacePoint再把两个方向的基函数做张量积得到曲面上对应参数 (u, v) 的点。function P bspSurfacePoint(cp, U, V, p, q, u, v) % cp: (n1) x (m1) x 3 控制点网格 n size(cp, 1) - 1; m size(cp, 2) - 1; % u 方向所有非零基函数 Nu zeros(1, n 1); for i 0:n Nu(i 1) Bbasis(i, p, u, U); end % v 方向所有非零基函数 Nv zeros(1, m 1); for j 0:m Nv(j 1) Bbasis(j, q, v, V); end % 张量积加权和控制点 P zeros(1, 3); for i 0:n for j 0:m P P Nu(i 1) * Nv(j 1) * squeeze(cp(i 1, j 1, :)).; end end end function N Bbasis(i, p, u, U) % 0 次基函数落在区间内为 1端点和半开区间单独处理 if p 0 if abs(u - U(end)) 1e-12 i numel(U) - p - 2 N 1; elseif u U(i 1) u U(i 2) N 1; else N 0; end return; end d1 U(i p 1) - U(i 1); d2 U(i p 2) - U(i 2); w1 0; w2 0; if d1 0 w1 (u - U(i 1)) / d1 * Bbasis(i, p - 1, u, U); end if d2 0 w2 (U(i p 2) - u) / d2 * Bbasis(i 1, p - 1, u, U); end N w1 w2; end代码里的d1、d2是相邻节点区间长度出现 0 表示重复节点对应项直接取 0避免除以 0。递归写法容易读但控制点超过 30×30 或者要做实时拖拽时建议改成先findSpan再 de Boor 循环求值的版本基函数只算局部一段速度会快一个量级。下面用 6×7 的控制点网格生成一张曲面cp zeros(6, 7, 3); [xg, yg] meshgrid(1:7, 1:6); cp(:,:,1) xg; cp(:,:,2) yg; cp(:,:,3) 0.35 * sin(xg / 7 * pi) .* cos(yg / 6 * pi); p 3; q 3; n size(cp, 1) - 1; m size(cp, 2) - 1; U [zeros(1, p 1), (1:n-p)/(n-p1), ones(1, p 1)]; V [zeros(1, q 1), (1:m-q)/(m-q1), ones(1, q 1)]; [ug, vg] meshgrid(linspace(0, 1, 60)); X zeros(size(ug)); Y X; Z X; for k 1:numel(ug) pt bspSurfacePoint(cp, U, V, p, q, ug(k), vg(k)); X(k) pt(1); Y(k) pt(2); Z(k) pt(3); end surf(X, Y, Z, FaceAlpha, 0.75, EdgeColor, none); hold on; mesh(cp(:,:,1), cp(:,:,2), cp(:,:,3), FaceColor, none, EdgeColor, [0.3 0.3 0.3]); view(3);运行后会看到一张平滑曲面包在控制点网格里。到这里自由曲面造型的正向链路已经通了控制点决定曲面形状节点向量决定控制点的影响范围次数决定整张曲面的连续性。把cp中某个控制点的 z 值改掉只有附近一小块曲面会变化这就是局部支撑的实际表现。2.4 结果检查先看曲面再看控制网曲面算出来后别急着进下一步。先做两个检查一是用mesh把控制网叠在曲面上确认控制网没有交叉、没有明显不合理的尖点二是检查曲面是否在控制点凸包内。B 样条曲面的每个点都是控制点的凸组合如果发现曲面跑到控制网外面十有八九是节点向量写错了或者控制点网格的维度顺序和期望不一致。max(max(abs(Z))) % 看高度是否在预期范围 plot3(cp(:,:,1), cp(:,:,2), cp(:,:,3), k.-, MarkerSize, 12); xlabel(x); ylabel(y); zlabel(z);如果只关心形状不关心参数化是否正确很容易把维度顺序颠倒的 bug 带进后续代码。在自由曲面造型源码里行列顺序一错拟合出来的曲面就是麻花。3. 散点到自由曲面的拟合源码参数反算与平滑控制正向设计是从控制点得到曲面逆向重构则是把点云或扫描数据变成控制点。这是自由曲面造型里另一个高频需求。实现方法分为插值和逼近。插值要求控制点数量和数据点数量相等曲面严格穿过每个点逼近则让控制点数量远小于数据点数量用最小二乘把噪声磨掉。扫描数据通常带噪声我会优先选逼近。3.1 参数化散点为什么不能直接画曲面散点要拟合 B 样条曲面第一步不是算基函数而是给每个散点分配一个参数坐标 (u, v)。只有拿到参数坐标才能把数据点写入曲面方程的张量积结构。对于按行列组织的扫描网格一种做法是直接用行列索引归一化作为参数更稳的做法是按行和列分别做弦长累积参数化避免数据点疏密不均时曲面出现抖动。% 对每一行做弦长参数化 u zeros(size(xg)); for r 1:size(xg, 1) seg sqrt(sum(diff(squeeze([xg(r,:); yg(r,:); zg(r,:)]), 1, 2).^2, 1)); u(r, :) cumsum([0, seg]); end u u / max(u(:));这里的假设是数据点已经按行、列组织好了。如果拿到的是完全无序的点云那就必须先做平面投影、测地距离或近邻图参数化这一步偏差比后面的拟合算法更影响结果。3.2 最小二乘加正则化的Matlab拟合源码下面代码把散点拟合成一张三次 B 样条曲面。核心是构造基函数矩阵 A然后对 x、y、z 三个坐标分量分别解一个带平滑正则的最小二乘问题。function [cp, U, V, fitP] fitBsplineSurface(P, u, v, p, q, nc, mc, lambda) % P: 数据点N x 3 % u, v: 每个数据点的参数坐标N x 1 % nc, mc: 两个方向控制点数量 % lambda: 平滑正则系数 U [zeros(1, p 1), (1:nc-p)/(nc-p1), ones(1, p 1)]; V [zeros(1, q 1), (1:mc-q)/(mc-q1), ones(1, q 1)]; A zeros(size(P, 1), nc * mc); for k 1:size(P, 1) Nu zeros(1, nc); Nv zeros(1, mc); for i 0:nc-1 Nu(i 1) Bbasis(i, p, u(k), U); end for j 0:mc-1 Nv(j 1) Bbasis(j, q, v(k), V); end A(k, :) kron(Nv, Nu); % 张量积基函数按控制点列优先展开 end % 二阶差分平滑矩阵惩罚相邻控制点之间过大的二阶变化 Dx sparse(diff(eye(nc), 2)); Dy sparse(diff(eye(mc), 2)); S kron(speye(mc), Dx * Dx) kron(Dy * Dy, speye(nc)); cpMat zeros(nc * mc, 3); for d 1:3 cpMat(:, d) (A * A lambda * S) \ (A * P(:, d)); end cp reshape(cpMat, [nc, mc, 3]); fitP A * cpMat; endA的每一行对应一个数据点每一列对应一个控制点记录的是该数据点在当前参数位置对每个控制点的敏感程度。kron(Nv, Nu)是张量积曲面结构的线性化写法。S是二阶差分矩阵它让相邻控制点不能变化太剧烈lambda 越大曲面越光顺但残差也越大。调用方式N 24; [xg, yg] meshgrid(linspace(-2, 2, N)); zg peaks(N) 0.08 * randn(N); P [xg(:), yg(:), zg(:)]; u (xg(:) 2) / 4; v (yg(:) 2) / 4; [cp, U, V, fitP] fitBsplineSurface(P, u, v, 3, 3, 8, 9, 5e-3); rmse sqrt(mean(sum((P - fitP).^2, 2)));拟合结束后cp就是你要的自由曲面造型控制点fitP是数据点对应参数位置下的曲面采样值。RMSE 用来判断拟合精度但我一般不会只看 RMSE而是把P和fitP叠在一起看最大偏差出现在哪里。梯度大的区域如果偏差集中说明那个方向的控制点数量不够。3.3 控制点数量与lambda的取值套路控制点数量和 lambda 没有一个万能公式但可以按数据规模和噪声水平给一组有效的初值。数据干净时控制点数量可以接近数据点数量数据带噪声时控制点数量要明显小于数据点数量否则噪声会被当成形状特征拟合进去。数据噪声水平控制点规模lambda 初值说明干净无噪声min(数据点数/4, 20) 左右0 或 1e-6保持拟合精度轻微噪声数据点数/41e-3 到 1e-2以光顺优先明显噪声数据点数/61e-2 到 0.1防止曲面波浪需要过点控制点数 数据点数0插值模式噪声大时不适用lambda 过小曲面会在噪声点附近出现凹陷lambda 过大真实特征也会被磨平。一个实用的检查方法是把 lambda 从 0 开始按 10 倍递增同时看 RMSE 和 max 误差随 lambda 的变化曲线。拐点附近的 lambda 通常就是可用的折中值。如果还想压制最大误差把\换成正则参数回归或lsqlin就可以在解里加入更多线性约束这是 Matlab 优化工具箱在自由曲面造型里最常见的用法。4. 交互拖拽控制点自由曲面造型过程中的参数调优正向造型的最后一步是手感。每次改控制点都重跑整个脚本看不到曲率变化也不适合现场调形。我会在 figure 上绑三个事件ButtonDownFcn选点WindowButtonMotionFcn移动点WindowButtonUpFcn释放点。这是 Matlab 原生交互里最直接的一种实现方式。4.1 可交互拖拽控制点的源码骨架三维视图里用鼠标直接改 z 值不直观所以我通常先在俯视图或侧视图里拖 x、yz 保持不变。需要改高度时切换到侧视图再拖一次。下面的骨架实现了“选最近控制点 移动 x、y 重绘曲面”的完整回路。function setupDragger(cp, U, V, p, q) fig gcf; ax gca; setappdata(ax, cp, cp); setappdata(ax, U, U); setappdata(ax, V, V); setappdata(ax, p, p); setappdata(ax, q, q); hControl plot3(ax, cp(:,:,1), cp(:,:,2), cp(:,:,3), ko); set(hControl, ButtonDownFcn, startDrag); function startDrag(~, ~) pt get(ax, CurrentPoint); cpd getappdata(ax, cp); dist (cpd(:,:,1) - pt(1,1)).^2 (cpd(:,:,2) - pt(1,2)).^2; [~, idx] min(dist(:)); [r, c] ind2sub(size(cpd(:,:,1)), idx); setappdata(ax, dragIdx, [r, c]); set(fig, WindowButtonMotionFcn, movePoint); set(fig, WindowButtonUpFcn, endDrag); end function movePoint(~, ~) idx getappdata(ax, dragIdx); cpd getappdata(ax, cp); pt get(ax, CurrentPoint); cpd(idx(1), idx(2), 1:2) pt(1, 1:2); setappdata(ax, cp, cpd); redrawSurface(ax, cpd, getappdata(ax,U), getappdata(ax,V), ... getappdata(ax,p), getappdata(ax,q)); end function endDrag(~, ~) set(fig, WindowButtonMotionFcn, ); set(fig, WindowButtonUpFcn, ); end endredrawSurface就是把第 2 章的求值逻辑再跑一遍用新的控制点刷新surf对象的XData、YData、ZData。这里没有重新生成曲线所以鼠标拖动时响应很快。如果控制点很多每次重算 60×60 的曲面也够了再大就要缓存基函数只在被移动控制点影响的局部范围重新求值。4.2 边界连续与切矢约束自由曲面造型在单个曲面内部一般没问题麻烦出在边界。两个曲面要拼接时只把边界控制点对齐还不够还要保证跨界切矢方向一致。拖拽边界控制点时如果只动第一行第二行控制点不动G1 连续性会立刻被破坏。所以我在交互逻辑里会先记录下来拖拽前的边界间隔向量d移动边界控制点时下面一行跟着平移到新位置d cp(2, :, :) - cp(1, :, :); cp(1, :, :) newBoundary; cp(2, :, :) cp(1, :, :) d;这能保证边界法向和切平面方向不突然改变。若要更高阶的 C2 连续就要把第三行控制点也纳入约束自由度更少调形也受限通常曲率连续曲面才会用到。4.3 调参时最容易踩的三个坑现象原因处理曲面在端点处露馅、不贴近控制网节点向量没有端点多重重数首尾节点重复 p1 次拖动一个控制点曲面整体都变控制点数量等于次数1局部支撑失效增加控制点数量保持次数不变曲面出现明显波浪lambda 过小或控制点过多加大 lambda减少控制点交互调整的本质是优化问题只是优化变量变成了控制点坐标。如果你想让曲面自动避开一个障碍物或者满足某个面积约束把鼠标拖拽换成fmincon也是同一条路径目标函数是曲率和控制点移动量约束是障碍物距离。源码里一旦有了控制点这个中间层正向设计和优化就共用一套数据结构。5. 导出与验证让Matlab自由曲面造型源码能交接曲面造型做完只是第一步能导出、能重建、能被下一个环节消费才算真正收口。这里重点说两个导出目标STL 网格和 MAT 参数文件。5.1 把曲面转成STL并导出先用 2.3 的求值函数生成足够密的采样网格再把网格拆成三角面片。曲面越复杂采样网格要越密一般展示用 60×60做结构仿真用 150×150 以上。STL 只记录三角面片不保存参数化信息所以它只适合渲染和打印不适合继续做曲面编辑。function writeStlBinary(filename, X, Y, Z) [nr, nc] size(X); V [X(:), Y(:), Z(:)]; F zeros((nr - 1) * (nc - 1) * 2, 3); k 1; for i 1:nr-1 for j 1:nc-1 a i (j - 1) * nr; b a 1; c i j * nr; d c 1; F(k, :) [a b c]; k k 1; F(k, :) [b d c]; k k 1; end end fid fopen(filename, w); fwrite(fid, zeros(1, 80), uint8); fwrite(fid, size(F, 1), uint32); for k 1:size(F, 1) v1 V(F(k, 1), :); v2 V(F(k, 2), :); v3 V(F(k, 3), :); nrm cross(v2 - v1, v3 - v1); nrm nrm / norm(nrm); fwrite(fid, [nrm, v1, v2, v3], float32); fwrite(fid, 0, uint16); end fclose(fid); end调用writeStlBinary(freeform.stl, X, Y, Z)即可生成二进制 STL。注意三角形顶点顺序要按右手定则排列否则法向量指向会反查看器里会出现黑面或透光现象。5.2 用MAT文件保存可复现的曲面参数STL 和 MAT 是互补关系。STL 给下游消费MAT 给自己和团队留后路。保存时不要只存 X、Y、Z 采样网格那样等于把自由曲面降级成了一堆点。控制点cp、节点向量U、V、次数p、q才是曲面的完整定义任何时刻都可以用第 2 章的bspSurfacePoint重新采样。save(freeform_surface.mat, cp, U, V, p, q);下次打开时load这个文件再调一次bspSurfacePoint曲面不丢。逆向拟合的结果还可以把 lambda 和 RMSE 一起存进去这样别人拿到源码能知道当前曲面是逼近还是插值可信度到多少。5.3 一个避免失真的技巧按弧长重采样如果原始数据点间距不均匀拟合后曲面的控制点分布也会跟着失配密集区域曲面抖动稀疏区域曲面过渡得太快。与其硬调 lambda不如先把数据沿采样方向按弧长重采样arc cumsum([0; sqrt(sum(diff(P, 1, 1).^2, 2))]); arc arc / arc(end); Pnew interp1(arc, P, linspace(0, 1, size(P, 1)));重采样后再做fitBsplineSurface控制点分布会均匀得多。这个技巧对车身曲面、螺旋桨叶片这类长条数据尤其有效。5.4 用离屏渲染验证曲面完整性最后导出前用figure(Visible,off)做一次离屏渲染确认曲面没有 NaN、没有被控制点网络穿透再出图存档。fig figure(Visible, off); surf(X, Y, Z, EdgeColor, none, FaceColor, interp); light(Position, [1 1 1]); camlight; print(fig, surface_preview.png, -dpng, -r150); close(fig);如果这一步出现大面积黑斑优先检查法向量朝向和数据点参数化如果只是局部褶皱回到 lambda 和控制点密度上继续调。本文还有配套的精品资源点击获取
