简介面向 MATLAB 用户与图形学初学者的曲线绘制函数包整合了贝塞尔曲线和 B 样条曲线的核心算法。贝塞尔部分基于 Casteljau 递推剖分实现控制点求值B 样条部分采用 De Boor 算法支持非均匀分段构造用户只需提供控制点坐标并设置曲线阶数或权重即可完成曲线的计算、绘制与可视化。两类曲线在实际使用中各有特点代码也体现了贝塞尔简单灵活与 B 样条局部控制更强的差异适合用于路径规划、动画设计及 CAD 形状建模等场景。压缩包内共 4 个文件全部为 .m 源码整体体积仅 2KB结构精简、函数划分清晰便于阅读和二次修改。目前已吸引 1105 人学习下载是理解两类曲线数学原理并快速上手 MATLAB 实现的实用参考。通过研读函数细节使用者能掌握控制点、基函数、递推求值之间的联系并在此基础上扩展出满足自身需求的曲线生成工具。1. 贝塞尔曲线与B样条曲线的Matlab函数B样条有官方户口贝塞尔要靠自己写很多人翻遍Matlab工具箱想在命令里找到一个直接叫“bezier”的函数翻到最后会发现B样条那边有spmak、fnplt、spcol一整套官方函数贝塞尔曲线反而像没户口的那个要么找到的是第三方代码要么得自己动手写。这个反差正好是下面这篇笔记的起点贝塞尔曲线与B样条曲线的Matlab函数到底怎么组织控制点怎么给、节点向量怎么设、官方函数能省多少事、手写函数该避开哪些坑都是绕不过去的问题。文章适合正在做轨迹规划、CAD小工具或控制多边形可视化的从业者读完可以拿到一套从求值到绘图的完整Matlab实现以及一组调参边界和踩坑记录。2. 先立理论贝塞尔与B样条的本质区别决定你在Matlab里怎么选函数2.1 贝塞尔曲线的整体控制性每个控制点都管着整条曲线贝塞尔曲线的数学定义是n次Bernstein基函数的线性组合P(t) Σ P_i * C(n,i) * t^i * (1-t)^(n-i)参数t落在[0,1]。关键在第i个基函数在参数区间内部处处非零这意味着你把任何一个控制点往旁边挪一点整条曲线的形状都会发生变化只是变化幅度随着t呈多项式衰减。这个性质在交互式绘图里叫全局拖拽直观但不好用你想微调曲线末端走向前端也会跟着动。从实现角度看这个定义也解释了为什么Matlab没有专门封装它。整个计算只需要组合数、幂运算和向量加点乘十几行代码就能写完官方再包一层反而别扭。我一般把Bernstein基写成独立小函数后面所有贝塞尔求值都复用它想换数据类型或者加导数输出只改这一个函数就够了。在轨迹规划场景里贝塞尔大多用来做一小段平滑过渡控制点往往只有四五个全局影响不是问题反而是优势。2.2 B样条曲线的局部支撑性节点向量比控制点更关键B样条把基函数换成了由节点向量knot vector驱动的分段多项式。每个基函数只在少数几个相邻节点区间内非零比如三次B样条基函数最多覆盖4个节点区间这个性质叫局部支撑。所以你的控制点只管得住自己附近那几段曲线远处的曲线几乎不受影响。这是CAD里B样条比贝塞尔更受欢迎的核心原因局部调整全局不动。代价是B样条多了一层“节点向量”配置knots的写法、重复度、单调性都直接影响曲线形状。节点向量给错轻则报错重则曲线不经过端点、首尾莫名拐弯。Matlab的曲线拟合工具箱为此提供了一整套操作spmak负责从节点向量和控制点构造样条对象fnplt负责绘制fnval负责求值fnder负责求导spcol负责算基函数矩阵。官方函数已经把最容易写错的递推过程封装好了要学的是怎么正确传参。2.3 一个关键换算贝塞尔其实是B样条在特定节点下的特例理论上一条n次贝塞尔曲线就是节点向量两端各重复n1次、中间没有内部节点的n次B样条。比如三次贝塞尔等价于knots[0 0 0 0 1 1 1 1]的三次B样条。这个性质我经常用来做交叉验证手写贝塞尔函数算出的点应该和spmak构造的退化B样条完全一致。后面第6章的验证脚本就建立在这个等价关系上。工程上我不建议真的用spmak去画贝塞尔两个理由。一是可读性差节点向量写成一排0和一排1读代码的人得反应一会儿才意识到这是三次贝塞尔二是没必要贝塞尔自己几行代码就搞定套B样条纯属增加理解成本。反过来B样条也不要试图用Bernstein基硬算阶数一高数值稳定性就绷不住。选型规则很朴素控制点少、做展示动画或需要导数闭合解时用贝塞尔函数控制点多、需要局部调整或拼接多条曲线时用B样条官方函数。下面这张表是我平时给同事的建议。需求场景推荐做法理由3到5个控制点画平滑曲线手写Bernstein基或de Casteljau代码短参数只有控制点和采样数局部修改某一小段曲线spmak加fnplt局部支撑只影响邻近节点区间大量离散点做平滑拟合spap2最小二乘或spapi插值样条拟合比贝塞尔逼近更可控两条曲线拼接要求导数连续B样条节点重复度控制重复节点让端点处导数对齐验证手写函数是否正确退化节点B样条交叉核对数学上等价直接数值比对3. 用Matlab写贝塞尔曲线求值函数从Bernstein直接到de Casteljau3.1 直接法Bernstein基函数与参数方向先解决基函数。把数学式子翻译成Matlab关键点是让t以列向量传入这样返回的基函数也是列向量后面和控制点做矩阵乘法时维度正好闭合。function B bernstein_basis(i, n, t) % i: 基函数序号从0到n % n: 贝塞尔曲线次数控制点数 n1 % t: 列向量参数采样点范围[0,1] B nchoosek(n, i) .* t.^i .* (1 - t).^(n - i); end逻辑说明nchoosek(n,i)返回组合数t是整段参数网格幂运算和(1-t)的幂逐项相乘最后得到每个采样点上的基函数值。参数上最容易出错的是t的方向如果传入的是行向量B也是行向量后面叠加控制点时维度就对不上我第一次写这个函数就在这里翻车折腾了十几分钟才发现是向量方向反了。直接法适用范围是n不超过15左右的曲线。阶数再高nchoosek(n,i)先膨胀成大整数随后t^i和(1-t)^(n-i)又可能下溢成0中间过程的数值波动会让曲线尾部出现毛刺。监控办法很简单把曲线端点的坐标和首尾控制点对比差超过1e-6就要警惕数值稳定性了。3.2 直接法组装的曲线函数控制点矩阵就是一行的坐标组合有了基函数曲线求值就是控制点的加权和。我习惯把控制点组织成(n1)行乘d列的矩阵二维曲线d2三维轨迹d3这样同一个函数既能算平面曲线也能直接扩展到三维路径。function P bezier_curve_direct(ctrl, numPts) % ctrl: 控制点矩阵(n1)行每行一个控制点 % numPts: 输出曲线点的数量 n size(ctrl, 1) - 1; t linspace(0, 1, numPts); % 列向量[0,1]之间均匀采样 P zeros(numPts, size(ctrl, 2)); for i 0:n B bernstein_basis(i, n, t); P P B * ctrl(i 1, :); end end逻辑说明外层循环遍历每个控制点bernstein_basis返回该控制点对应的基函数在这些t上的取值列向量B乘以单个控制点行向量结果是numPts乘d的矩阵累加后就是整条曲线。这里只解释一个参数n由控制点个数决定控制点4个n就是3对应三次贝塞尔想提高曲线复杂度只能加控制点而不是改某个“次数”参数这是新手最容易误解的地方。这种写法优点是直观缺点是每求一个点都要完整算一遍n1个基函数。numPts取200时性能无所谓但要做几千个点的密集采样或者实时拖动控制点效率就有点吃紧。更稳的做法是下面这套de Casteljau递推。3.3 递推法de Casteljau算法顺便白拿切线de Casteljau的核心思想是反复线性插值相邻控制点连线按比例t切分得到新控制点再对新控制点重复操作直到只剩一个点。这个点就是曲线在t处的值。它比直接法数值稳定得多而且最后一轮递推时那两个中间点的连线方向乘以次数n就是曲线在t处的切向量。这个副产品在轨迹规划里非常有用后面验证脚本就会用到。function [P, dP] bezier_decasteljau(ctrl, numPts) % ctrl: 控制点矩阵(n1)行 % numPts: 采样点数 % P: 曲线上点的坐标 % dP: 曲线对参数t的导矢即切线向量 n size(ctrl, 1) - 1; t linspace(0, 1, numPts); P zeros(numPts, size(ctrl, 2)); dP zeros(numPts, size(ctrl, 2)); for k 1:numPts pts ctrl; % 当前层控制点 for r 1:n if r n a pts(1, :); % 最后一轮递推前的两个点 b pts(2, :); end pts (1 - t(k)) * pts(1:end-1, :) t(k) * pts(2:end, :); end P(k, :) pts(1, :); if nargout 1 dP(k, :) n * (b - a); end end end逻辑说明最内层循环每走一步点数组行数减1。当r等于n时先备份当前层的第1、2行因为这一步递推之后数组只剩1行切线信息只能从上一层的两个点之间恢复。dP是曲线对参数t的导矢是速度向量不是单位切向量后面做弧长参数化还要再除以模长。参数numPts决定输出密度起步值建议200拖动交互场景可以降到50左右换实时性。3.4 统一入口函数留一个开关方便换算法实际项目里我很少让业务代码直接调用上面两个函数而是封装一个入口用mode参数切换算法。这样控制点个数超过15时可以直接切到递推法不用改任何业务代码。function [P, dP] bezier_curve(ctrl, numPts, mode) % mode: direct 使用Bernstein直接法casteljau 使用递推法 if nargin 3 mode casteljau; % 默认递推数值更稳 end if strcmp(mode, direct) if nargout 1 error(direct模式不支持求导请改用casteljau); end P bezier_curve_direct(ctrl, numPts); dP []; else [P, dP] bezier_decasteljau(ctrl, numPts); end end我一般默认用递推法不是因为直接法不能用而是递推多一个免费切线输出。调用示范ctrl [0 0; 1 2; 3 -1; 5 1]; [P,dP] bezier_curve(ctrl, 200, casteljau); 之后plot(P(:,1), P(:,2), b-)再叠加plot(ctrl(:,1), ctrl(:,2), ro--)就能同时看到曲线和控制多边形。控制点矩阵是4行2列输出P是200行2列dP同尺寸这套尺寸约定在后面对比B样条时也会沿用。4. 用Curve Fitting Toolbox画B样条spmak与fnplt的最小可跑组合4.1 最小可跑命令三行代码画一条三次B样条B样条官方函数的用法比贝塞尔更强调参数配对正确。核心函数spmak接收两个参数节点向量knots和控制点coefs两者关系是节点个数等于控制点个数加次数再加1。以6个控制点、三次样条为例节点向量需要63110个元素。% 6个控制点三次B样条order4节点向量10个 knots [0 0 0 0 1 2 3 3 3 3]; % 两端各重复4次保证过首尾控制点 ctrl [0 0; 1 2; 2 -1; 3 1; 4 0; 5 2]; spl spmak(knots, ctrl); fnplt(spl); hold on; plot(ctrl(:,1), ctrl(:,2), ro--);这段代码背后有两个容易忽略的细节。第一knots必须单调不减[0 0 0 0 1 2 3 3 3 3]这种两端重4次的写法称为Clamped节点它保证曲线从第一个控制点出发、在最后一个控制点结束如果两端不重复曲线端点会悬空看起来像飘在控制多边形外面。第二ctrl的行数必须和knots长度满足上面的计数关系少了直接报错多了会静默产生额外样条段数值结果完全不是预期形状。fnplt默认按样条参数区间接细化绘图拐弯剧烈的地方偶尔显得粗糙。我习惯改用fnval手动求值tt linspace(0, 3, 300); pts fnval(spl, tt); % 返回2×300第1行是x第2行是y plot(pts(1,:), pts(2,:), b-);逻辑说明fnval返回的数据布局是维度在行、点在列和贝塞尔函数里“每行一个点”的约定正好相反。跨函数对比时最容易在这里踩坑需要先转置再比对或者干脆统一在自己写的封装里做转置。这段代码里tt的范围要和样条定义域匹配定义域可以从fnbrk(spl, interval)查到超出定义域的求值会得到NaN。4.2 节点向量构造用augknt避免手写重复数错手写knots最烦的是数重复次数三次样条两端各重复4次写错一个数整条曲线就变形。有个官方函数叫augknt它按指定的重复度把断点序列扩展成Clamped节点向量省掉手数重复的功夫。knots augknt([0 1 2 3], 4); % 内部节点1、2各重复1次端点0、3各重复4次 ctrl [0 0; 1 2; 2 -1; 3 1; 4 0; 5 2]; spl spmak(knots, ctrl);augknt的第一个参数是断点序列第二个参数是端点的重复次数不是曲线次数本身。次数是3重复度是4因为B样条阶数order等于次数加1。内部断点各重复1次所以[0 1 2 3]配合重复度4生成10个节点正好匹配6个控制点。如果只写augknt([0 3], 4)生成的节点只有8个只够配4个控制点拿6个控制点去spmak必然报错。注意augknt第二个参数是重复度要和order一致。三次样条order是4就写4不是写3。这个对应关系写反是高频错误。调试时我推荐一条经验先确定控制点个数和次数再用numel(knots) size(ctrl,1) 3 1反推节点数量对不上就逐个检查内部断点有没有写漏。曲线拟合工具箱里所有B样条函数最终都绕不开这个配对关系花5分钟把计数看懂省下来的是排错时间。4.3 用spcol把基函数画出来理解局部支撑最快的办法说了半天局部支撑不如直接画一眼。spcol(knots, order, t)返回B样条基函数在参数t处的取值矩阵每一列对应一个基函数非零区间就告诉你局部支撑的范围。knots [0 0 0 0 1 2 3 3 3 3]; order 4; % 三次样条阶数等于次数加1 tt linspace(0, 3, 500); B spcol(knots, order, tt); % 500×66个基函数 plot(tt, B(:, 3), b-, tt, B(:, 4), r-);逻辑说明第3列和第4列的基函数只在1和2附近的节点区间内非零远离这个区间就趋近于0。也就是说想调整曲线第3段附近形状时只需动第3、第4个控制点其他控制点怎么挪都不影响这一段。这正是B样条在CAD里成为工业标准的底气也是你决定局部修改该用哪个函数的判断依据。spcol在工程里还有更实际的用途构造最小二乘拟合矩阵。逆向工程时把测量点代入spcol得到基函数矩阵再用矩阵左除解控制点比直接用现成拟合函数更可控因为节点位置可以自己指定一旦需要手动加密节点这个做法就特别顺手。5. 贝塞尔与B样条Matlab实现的避坑5类高频问题定位与处理5.1 现象贝塞尔曲线整体偏离控制多边形在中间“鼓包”有人拿到曲线后发现它并不贴近控制点怀疑代码写错了。实际上贝塞尔是逼近不是插值除了首尾两个控制点中间控制点通常不在曲线上曲线只是被“吸引”过去。这个鼓包在控制点间距不均匀时尤其明显。原因就是Bernstein基函数的全局支撑性所有控制点共同作用谁也不能单独锁定某个过点。解决方案是先分清需求如果只是画平滑示意保持现状就好如果要求曲线精确穿过一系列点请改用样条插值比如spapi(4, x, y)不要再继续调贝塞尔参数那是白花时间。5.2 现象spmak报错提示节点向量与控制点不匹配这是B样条函数最常见的报错原因几乎都是同一个节点数量不等于控制点数加次数再加1。排查时不要凭感觉数直接在命令行打印numel(knots)和size(ctrl,1)再看numel(knots) - size(ctrl,1) - 4是否等于0三次样条按4算。不等于0就去检查内部节点有没有漏写。另一个常见变体是knots非单调递增比如手写时把1改成0.5又改回1导致后一个节点小于前一个节点spmak一样会报错。解决的办法是把节点序列重新排序并去重检查或者干脆统一用augknt生成别手写。5.3 现象fnval求值结果和plot画出的曲线对不上fnval返回的数据布局是维度在行、点在列和大多数人习惯的“每行一个点”完全相反。很多人直接pts fnval(spl, tt)然后取pts(1,:)当成第一个点结果画出来乱七八糟。这个错位在跨函数对比时很隐蔽只报维度不报数值。解决办法是统一规范自己写的函数一律每行一个点官方函数返回值先转置再进自己的流程写成pts fnval(spl, tt)后面所有代码按这个约定写就不会再犯。这个坑我在给同事review代码时见过至少三次每次都是转置问题。5.4 现象B样条端点没有落在首尾控制点上或者走向完全相反Clamped节点向量两端各重复order次是曲线过端点控制点的前提。先检查knots首尾重复次数三次样条开头要有4个连续相同值结尾也要有4个连续相同值少了就会悬空。多余4个没问题少了必出问题。还有一种情况是控制点顺序给反了曲线本身没毛病但起点终点互换了位置。解决方法是把ctrl用flipud(ctrl)倒过来再重新spmak。判断是不是顺序问题最简单的是打印ctrl的第一行和最后一行看哪个是你想要的起点。5.5 现象控制点超过15个后贝塞尔曲线出现锯齿和NaNBernstein直接法在阶数较高时组合数和幂运算一起出问题表现为曲线尾部在0和无穷之间震荡甚至出现NaN。原因和3.1节说的一致中间计算先膨胀后下溢浮点精度扛不住。解决方法是换de Casteljau递推它有数值上的天然优势如果控制点几十上百直接换spmak做B样条拟合不要再和贝塞尔死磕。我一般以15为分界线少于15个点直接法够用超过15个点闭眼用de Casteljau控制点上百直接用spmak这是避免数值玄学最省心的策略。实际排查顺序我总结为一句先看端点再看中间点最后看导数。端点错是参数错中间点错是算法错导数错是布局或方向错。按这个顺序走大部分贝塞尔和B样条的实现问题都能在十分钟内定位。6. 工程验证把贝塞尔函数和B样条函数交叉核对少走调参弯路最后分享一个我每次落地曲线代码时必做的验证技巧。原理就是第2章的数学关系当B样条节点向量退化为两端各重复n1次、内部无节点时它就是一条n次贝塞尔曲线。利用这个等价关系可以把手写函数和官方函数放在同一套控制点上互检谁错立刻暴露。% 用同一组控制点构造贝塞尔型退化B样条 ctrl [0 0; 1 2; 3 -1; 5 1]; knots augknt([0 1], 4); % 三次两端重复4次无内部节点 spl spmak(knots, ctrl); [P_bez, dP_bez] bezier_curve(ctrl, 200, casteljau); P_bspl fnval(spl, linspace(0, 1, 200)); maxDiff max(abs(P_bez(:) - P_bspl(:))); disp(maxDiff); % 期望输出接近1e-12若在1e-3量级说明控制点或端点有问题跑通坐标对比后再进一步把de Casteljau算出的切线向量和官方fnder的导矢对比验证方向一致性。vector direction写反是这类函数最容易出的隐蔽错误数值上完全看不出来但后续轨迹规划会直接“倒车”。dP_bspl fnval(fnder(spl), linspace(0, 1, 200)); maxDiffTangent max(abs(dP_bez(:) - dP_bspl(:)));两条都通过手写函数才算真正可信。我早年在某跨平台系统里调试轨迹规划模块时就因为没做这个交叉验证把de Casteljau的切线方向写反了导致后续速度规划全线倒车排查两天才发现是b减a的方向写成了a减b。从那以后我的习惯就是所有手写曲线函数落地前必须和官方B样条函数对一次数值对不上的那一方九成是把符号或转置搞反了。这个习惯帮我省掉了很多调参时间希望也能帮到你让贝塞尔曲线与B样条曲线的Matlab函数真正成为项目里可以放心依赖的底层工具。本文还有配套的精品资源点击获取