曲线这个东西在计算机图形学里属于那种你天天用但未必真懂的基础设施。做UI的调个圆角、做动画的拉个缓动、做建模的捏个曲面背后全是贝塞尔和B样条在撑着。但很多人对它们的理解停留在拖控制点的层面一旦遇到需要自己实现曲线求值、判断曲率、做曲线拟合的场景就抓瞎了。这篇东西就是把我这些年在这两种曲线上踩过的坑、算过的公式、写过的代码从头到尾捋一遍。不管你是刚学计算机图形学的学生还是工作中需要手撸曲线算法的工程师应该都能从中找到点有用的东西。1. 从一条拉不直的曲线说起1.1 为什么需要参数曲线先说个最朴素的问题给你一堆点让你画一条光滑的曲线穿过它们你怎么画最直接的想法是多项式插值——找一个n-1次多项式让它精确穿过n个点。数学上完全可行但实际用起来是灾难。龙格现象了解一下高次多项式在区间边缘会产生剧烈振荡点越多振荡越厉害。你给10个点它给你画出一条上下翻飞的蛇。另一个思路是分段。每两个相邻点之间用一段低次曲线连起来但这样在连接处会出现折角——切线不连续看起来就是一段一段硬拼的。要让它光滑就得保证连接点处导数连续甚至二阶导数连续。这就是样条曲线的核心动机。参数曲线的好处在于它把x和y三维就是x、y、z都表示成独立参数t的函数P(t) (x(t), y(t), z(t))这样做的好处是曲线可以自交、可以竖直、可以闭合不受一个x对应一个y的限制。贝塞尔和B样条都是参数曲线的典型代表。1.2 贝塞尔和B样条的关系一句话说清很多人搞不清这俩的关系。我打个比方贝塞尔曲线像是全局投票制每个控制点对整条曲线都有影响改一个点整条曲线都跟着动B样条像是局部代表制每个控制点只影响它附近的一小段改一个点远处纹丝不动。这个区别直接决定了两者的适用场景。做设计工具、做字体轮廓贝塞尔够用且直观做工业建模、做复杂曲面拼接B样条才是正解。NURBS非均匀有理B样条则是在B样条基础上引入了权重能精确表示圆、椭圆这些二次曲线是CAD领域的绝对主力。2. 贝塞尔曲线直觉最好但暗藏陷阱2.1 德卡斯特里奥算法递归就是一切贝塞尔曲线的定义有好几种写法但真正适合写代码的是德卡斯特里奥de Casteljau算法。它的思想极其优雅在控制多边形上反复做线性插值最后收敛到曲线上的一个点。给定n1个控制点P0到Pn要算参数t处的点对每一对相邻点做插值 P_i^(1) (1-t) * P_i t * P_(i1) 然后对新的点集重复这个过程直到只剩一个点。用Python写出来大概长这样def de_casteljau(points, t): pts [p[:] for p in points] # 复制一份别改原数据 n len(pts) for r in range(1, n): for i in range(n - r): pts[i] [ (1 - t) * pts[i][0] t * pts[i 1][0], (1 - t) * pts[i][1] t * pts[i 1][1] ] return pts[0]这个算法的数值稳定性非常好因为它只做凸组合系数非负且和为1不会出现除零或者溢出。相比之下直接用伯恩斯坦基函数展开求值在控制点很多的时候数值误差会明显增大。提示如果你只是要算曲线上的点永远优先用德卡斯特里奥别去展开多项式。这是数值计算的基本素养。2.2 伯恩斯坦基函数理解贝塞尔性质的钥匙德卡斯特里奥算法好写但要理解贝塞尔曲线的性质得看伯恩斯坦基函数。n次贝塞尔曲线可以写成B(t) Σ C(n,i) * t^i * (1-t)^(n-i) * P_i其中C(n,i)是组合数。这些基函数有几个关键性质非负性每项都≥0单位分解所有基函数加起来恒等于1对称性B_i,n(t) B_(n-i),n(1-t)单位分解这个性质直接导致了贝塞尔曲线的一个重要特性——凸包性。曲线永远落在控制点构成的凸包内部。这个性质在碰撞检测、裁剪算法里非常有用如果凸包不相交曲线肯定不相交可以快速排除。但伯恩斯坦基函数是全局的每个基函数在整个[0,1]区间上都非零。这就是为什么贝塞尔曲线牵一发动全身——移动一个控制点所有基函数的值都变了整条曲线都受影响。2.3 曲率怎么算从一阶导到二阶导热词里有人问贝塞尔曲线有没有判断曲率的方法答案是当然有而且公式很标准。参数曲线的曲率公式是κ |xy - yx| / (x^2 y^2)^(3/2)对于贝塞尔曲线导数可以通过对控制点做差分来算。n次贝塞尔曲线的一阶导是n次实际是n-1次贝塞尔曲线控制点是n*(P_(i1) - P_i)。二阶导同理是对一阶导的控制点再做一次差分。def bezier_derivatives(points): n len(points) - 1 # 一阶导控制点 d1 [[n * (points[i1][0] - points[i][0]), n * (points[i1][1] - points[i][1])] for i in range(n)] # 二阶导控制点 d2 [[(n-1) * (d1[i1][0] - d1[i][0]), (n-1) * (d1[i1][1] - d1[i][1])] for i in range(n-1)] return d1, d2然后在参数t处分别求值代入曲率公式即可。这里有个坑在曲线端点处如果一阶导为零比如控制点重合曲率公式会除零。实际代码里要加保护或者用极限的方式处理。曲率在工程上有什么用最典型的是道路设计——曲率不能突变否则车辆行驶会不舒服还有字体渲染——曲率大的地方需要更密的采样点否则看起来会有棱角。2.4 贝塞尔的致命短板局部修改的代价假设你有一条20个控制点的贝塞尔曲线现在只想把中间某一段稍微调整一下。对不起做不到。你动任何一个控制点整条曲线都会变。这在交互式设计里是灾难性的。更麻烦的是贝塞尔曲线的次数等于控制点数减一。控制点一多次数就高高次贝塞尔曲线不仅计算量大而且数值稳定性差还容易出现控制点明明很规整曲线却扭来扭去的情况。这两个问题B样条都解决了。3. B样条工业级曲线的主力3.1 节点向量B样条的时间轴B样条最让人困惑的概念就是节点向量knot vector。我换个说法它定义了每个控制点什么时候上场、什么时候下场。一个p次B样条节点向量是一个非递减的实数序列U {u_0, u_1, ..., u_m}其中m n p 1n是控制点数减一。节点向量的长度 控制点数 次数 1。节点向量分几种类型类型特点用途均匀节点等间距理论分析简单场景准均匀两端重复p1次曲线过端点常用非均匀任意非递减最灵活NURBS基础分段贝塞尔内部节点重复p次等价于多段贝塞尔拼接准均匀节点向量是最常用的因为它保证曲线经过首末控制点同时内部又是均匀的。比如3次B样条4个控制点准均匀节点向量是{0,0,0,0,1,1,1,1}——这其实就是一条贝塞尔曲线。3.2 考克斯-德布尔递推B样条的核心算法B样条的基函数由考克斯-德布尔Cox-de Boor递推公式定义N_i,0(u) 1 if u_i u u_(i1), else 0 N_i,p(u) (u - u_i)/(u_(ip) - u_i) * N_i,p-1(u) (u_(ip1) - u)/(u_(ip1) - u_(i1)) * N_i1,p-1(u)这个递推看起来吓人但逻辑很清晰p次基函数是由两个p-1次基函数加权组合而成的。权重是两个线性函数分别对应从左边界靠近和从右边界远离。写代码的时候要注意分母为零的情况——当节点重复时分母可能为0。标准做法是约定0/0 0。def cox_de_boor(i, p, u, knots): if p 0: return 1.0 if knots[i] u knots[i1] else 0.0 left 0.0 right 0.0 denom1 knots[ip] - knots[i] if denom1 1e-12: left (u - knots[i]) / denom1 * cox_de_boor(i, p-1, u, knots) denom2 knots[ip1] - knots[i1] if denom2 1e-12: right (knots[ip1] - u) / denom2 * cox_de_boor(i1, p-1, u, knots) return left right这个递归实现直观但效率低实际工程中会用迭代版本或者预计算基函数表。3.3 局部支撑性B样条最值钱的性质B样条基函数有个关键性质叫局部支撑N_i,p(u)只在区间[u_i, u_(ip1))上非零。这意味着移动一个控制点只影响p1个节点区间内的曲线段。这个性质的价值怎么强调都不过分。它意味着修改局部不影响全局交互式编辑友好曲线次数可以独立于控制点数想加多少控制点都行次数不变数值计算可以局部化效率高举个例子一条3次B样条100个控制点。移动第50个控制点只有第47到第53个控制点对应的曲线段会变其余部分完全不动。这在贝塞尔曲线里是不可想象的。3.4 德布尔算法B样条版的德卡斯特里奥和贝塞尔有德卡斯特里奥一样B样条有德布尔de Boor算法。它的结构和德卡斯特里奥几乎一样只是插值参数不是固定的t而是根据节点向量来确定的。def de_boor(u, points, knots, p): # 找到u所在的节点区间 k find_span(u, points, knots, p) # 复制p1个控制点 d [points[k-pi][:] for i in range(p1)] for r in range(1, p1): for i in range(p, r-1, -1): alpha (u - knots[k-pi]) / (knots[k1i-r] - knots[k-pi]) d[i] [(1-alpha)*d[i-1][j] alpha*d[i][j] for j in range(len(d[i]))] return d[p]德布尔算法同样数值稳定是B样条求值的标准方法。find_span函数负责定位参数u落在哪个节点区间可以用二分查找加速。4. 从理论到代码手撸一个曲线库4.1 整体架构设计光看公式容易飘我把自己写的一个小曲线库的结构分享一下。核心就三个类BezierCurve贝塞尔曲线存控制点列表BSplineCurveB样条曲线存控制点、次数、节点向量CurveUtils工具函数采样、求导、曲率、弧长设计上有个取舍要不要把贝塞尔做成B样条的特例理论上可以分段贝塞尔就是内部节点重复p次的B样条但实际用起来贝塞尔有更简单的专用算法硬套B样条反而绕。所以我选择分开实现共享工具函数。4.2 采样与自适应细分把曲线画出来最直接的方法是均匀采样t从0到1每隔0.01取一个点连成折线。但这样有个问题曲率大的地方采样太稀看起来有棱角曲率小的地方采样太密浪费。更好的做法是自适应细分递归地把曲线段一分为二如果这一段足够平就停止否则继续分。判断平的标准可以用控制点到弦的最大距离。def adaptive_sample(curve, tolerance0.5, max_depth10): result [] def subdivide(t0, t1, depth): p0 curve.evaluate(t0) p1 curve.evaluate(t1) pm curve.evaluate((t0 t1) / 2) # 中点到弦的距离 d point_line_distance(pm, p0, p1) if d tolerance or depth max_depth: result.append((t0, p0)) result.append((t1, p1)) else: subdivide(t0, (t0t1)/2, depth1) subdivide((t0t1)/2, t1, depth1) subdivide(0, 1, 0) return result这个方法的采样密度自动适应曲率画出来的曲线既光滑又不会浪费顶点。实测下来tolerance取0.5像素左右视觉效果就很好了。4.3 弧长参数化匀速运动的关键如果你想让一个物体沿着曲线匀速运动直接用t做参数是不行的——t均匀变化时曲线上的实际速度不均匀。需要做弧长参数化。基本思路是先用数值积分算出曲线的总弧长然后建立t和弧长s的映射表运动时按s均匀推进再反查t。def build_arc_length_table(curve, samples1000): table [(0.0, 0.0)] prev curve.evaluate(0) total 0.0 for i in range(1, samples 1): t i / samples curr curve.evaluate(t) total distance(prev, curr) table.append((t, total)) prev curr return table, total def t_at_arc_length(table, total, s): target s * total # 二分查找 lo, hi 0, len(table) - 1 while lo hi - 1: mid (lo hi) // 2 if table[mid][1] target: lo mid else: hi mid # 线性插值 t0, s0 table[lo] t1, s1 table[hi] if s1 - s0 1e-12: return t0 return t0 (target - s0) / (s1 - s0) * (t1 - t0)这个表建一次可以反复用对于动画路径来说非常划算。注意采样数要够否则反查精度不够运动会有轻微抖动。4.4 曲线拟合从离散点到控制点实际工作中经常遇到反问题给一堆离散点求一条逼近它们的B样条曲线。这就是曲线拟合。最常用的方法是最小二乘拟合。步骤大致是对数据点做参数化均匀、弦长、向心都行确定节点向量通常用平均法或均匀法建立最小二乘方程组求解控制点弦长参数化是最常用的因为它考虑了点的实际分布u_0 0 u_i u_(i-1) |P_i - P_(i-1)| 最后归一化到[0,1]节点向量的选择对拟合质量影响很大。平均法averaging是比较稳妥的选择u_(jp) (1/p) * Σ u_(i) for i j to jp-1拟合的精度取决于控制点数量。控制点太少拟合不上去太多会过拟合曲线出现不必要的波动。实践中一般从较少控制点开始逐步增加直到误差满足要求。5. 那些文档里不会写的坑5.1 节点向量的端点重复问题写B样条代码最容易翻车的地方就是节点向量。准均匀节点向量要求两端各重复p1次比如3次B样条节点向量开头必须是{0,0,0,0}结尾是{1,1,1,1}。少一个重复曲线就不过端点多一个分母就出零。我见过有人把节点向量写成{0,0,0,1,1,1}然后纳闷为什么曲线不经过第一个控制点。原因就是3次B样条需要4个重复他只写了3个。注意节点向量的长度必须等于控制点数 次数 1。这个关系式是硬约束写代码时最好加个断言检查。5.2 曲率计算的数值陷阱前面提过曲率公式在导数为零时会炸。但实际中还有更隐蔽的问题当曲线接近直线时二阶导很小曲率计算会有很大的相对误差。我的做法是加一个阈值判断如果一阶导的模长小于某个epsilon直接返回曲率为零如果二阶导也很小同样返回零。虽然理论上不严谨但工程上足够用而且避免了NaN传播。另一个坑是曲率符号。二维曲线的曲率有正负之分表示弯曲方向。公式里的绝对值去掉分子xy - yx的符号就是曲率符号。做道路设计或者字体渲染时这个符号很重要。5.3 高次B样条的数值稳定性虽然B样条理论上次数可以任意高但实践中很少超过5次。原因和高次多项式一样数值稳定性下降而且高次基函数的支撑区间很宽局部性变差。工业界最常用的是3次B样条C2连续偶尔用4次或5次。如果你发现自己需要10次B样条大概率是建模方式有问题应该考虑用多段低次曲线拼接。5.4 重节点与连续性节点重复次数直接影响曲线的连续性。一个p次B样条如果某个内部节点的重复度是r那么曲线在该节点处的连续性降为C^(p-r)。重复度1C^(p-1)连续标准情况重复度pC^0连续曲线在该点有折角重复度p1曲线断开这个性质在建模时非常有用。想要一条有尖角的曲线把某个节点重复p次就行。想要曲线断开成两段重复p1次。5.5 浮点数比较的坑节点区间的查找、参数是否在[0,1]范围内、两个点是否重合——这些判断在浮点数下都不能用。我一般用1e-10作为epsilon但具体值要根据数据尺度调整。数据范围是毫米级和公里级epsilon肯定不一样。一个实用的技巧是所有参数化都归一化到[0,1]这样epsilon可以固定用1e-10不用每次调整。6. NURBSB样条的完全体6.1 权重到底在干什么NURBS在B样条基础上给每个控制点加了一个权重w_i。基函数变成R_i,p(u) N_i,p(u) * w_i / Σ N_j,p(u) * w_j权重的物理意义是吸引力权重越大曲线越靠近该控制点。当所有权重相等时NURBS退化为普通B样条。权重最强大的地方在于它能精确表示圆锥曲线。一个标准的圆用普通B样条只能逼近用NURBS可以精确表示。这是NURBS成为CAD标准的核心原因。6.2 齐次坐标NURBS的实现技巧直接按上面的公式实现NURBS每次求值都要算分母效率低。标准做法是用齐次坐标把控制点升维(x, y, z) - (wx, wy, wz, w)在四维空间做普通B样条求值把结果投影回三维(X, Y, Z, W) - (X/W, Y/W, Z/W)这样所有B样条的算法德布尔、求导、细分都能直接复用代码量大大减少。这个技巧在OpenGL的NURBS求值器里也是这么用的。6.3 NURBS的求导NURBS的导数不能直接套B样条的导数公式因为分母里有权重。标准做法是用商法则C(u) A(u) / w(u) C(u) (A(u) - w(u) * C(u)) / w(u)其中A(u)是齐次坐标下的B样条曲线w(u)是权重基函数的和。高阶导数可以递推。实现的时候先算齐次空间里的导数再转换回来。7. 实际项目中的选型建议7.1 什么时候用贝塞尔什么时候用B样条这个问题我被问过很多次。我的经验是场景推荐理由UI圆角、图标贝塞尔控制点少直观浏览器原生支持字体轮廓贝塞尔二次/三次行业标准渲染器优化好动画缓动贝塞尔cubic-bezier就够用工业建模NURBS精度要求高需要精确圆锥曲线路径规划B样条需要局部修改C2连续保证平滑数据拟合B样条控制点数灵活局部支撑简单说控制点少、需要直观交互用贝塞尔控制点多、需要局部修改、需要高阶连续用B样条。7.2 次数选择的经验法则2次C1连续够用但不够光滑适合快速原型3次C2连续工业标准99%的场景选它4次及以上特殊需求比如需要C3连续的气动外形3次是甜点因为C2连续对人眼来说已经足够光滑而且计算量适中。再高次数视觉上提升不明显计算和数值稳定性都变差。7.3 控制点数量的权衡控制点越多曲线越灵活但计算量线性增长过拟合风险增加交互编辑变复杂我的经验是先用尽量少的控制点看拟合误差不够再加。对于拟合问题一般从数据点数的1/3到1/2开始试。8. 调试曲线的实用手段8.1 可视化控制多边形和基函数调试曲线问题第一件事是把控制多边形画出来。曲线一定在控制多边形的凸包内如果曲线跑到外面去了肯定是代码有bug。第二件事是把基函数画出来。B样条的基函数图能直观显示每个控制点的影响范围。如果基函数图看起来不对比如有负值、不连续、支撑区间不对那节点向量肯定有问题。8.2 用已知特例验证写曲线库的时候我习惯用几个已知特例做单元测试所有控制点重合曲线应该退化成一个点控制点共线曲线应该是一条直线段准均匀节点向量控制点数次数1应该等价于贝塞尔曲线对称的控制点曲线应该对称这些测试能抓住大部分实现错误。8.3 数值验证导数用差分校验解析求导容易写错用数值差分校验是个好办法f(x) ≈ (f(xh) - f(x-h)) / (2h)h取1e-6左右比较解析导数和差分导数的相对误差。如果误差大于1e-3基本可以确定解析求导有问题。这个方法我用来抓过好几次符号错误。8.4 曲率的可视化曲率不好直接看但可以画曲率梳curvature comb在每个采样点处沿法线方向画一条长度正比于曲率的线段。曲率梳光滑说明曲线质量好曲率梳有尖刺说明曲率不连续曲率梳交叉说明有拐点。这个工具在汽车外形设计里是标配做任何高质量曲线都应该用。9. 性能优化从能用 to 好用9.1 预计算基函数表如果曲线要反复求值比如动画每帧都要算预计算基函数表能大幅提速。把[0,1]分成N份每份预先算好所有非零基函数的值求值时直接查表加插值。代价是内存和初始化时间但对于实时应用来说完全值得。N取1000左右精度和速度的平衡就很好。9.2 利用局部支撑做裁剪B样条的局部支撑意味着对于给定的参数u只有p1个基函数非零。求值时只需要遍历这p1个控制点而不是全部。控制点越多这个优化的收益越大。实现上先用二分查找定位u所在的节点区间然后只取对应的p1个控制点参与计算。9.3 SIMD与并行曲线求值天然适合并行不同参数t之间没有依赖。用SIMD指令一次算4个或8个点或者用多线程分摊采样任务都能获得接近线性的加速。不过要注意德布尔算法内部有依赖每一轮依赖上一轮的结果所以并行应该在参数维度上做而不是在算法内部做。10. 写在最后的一些个人体会曲线这东西公式看着多但核心思想就那么几个用基函数加权控制点用节点向量控制基函数的形状用递推算法保证数值稳定。把这三个东西吃透贝塞尔、B样条、NURBS就是同一套框架的不同配置。我刚开始学的时候最大的误区是死记公式。后来发现真正有用的是理解每个性质背后的几何直觉——为什么凸包性成立为什么局部支撑重要为什么重节点会降连续性。这些想明白了公式忘了都能推出来。另一个体会是一定要自己动手实现一遍。看十遍德布尔算法不如自己写一遍然后调通。调通过程中遇到的那些坑——节点向量长度不对、分母为零、端点不过点——才是真正长在身上的知识。最后说个实用的如果你只是要用曲线不想自己实现Python里scipy的interpolate模块有BSpline类JavaScript里d3有曲线生成器C里OpenCASCADE有完整的NURBS实现。但如果你要理解曲线还是得自己撸一遍。
贝塞尔与B样条曲线：从原理到代码实现，避开工程中的坑
曲线这个东西在计算机图形学里属于那种你天天用但未必真懂的基础设施。做UI的调个圆角、做动画的拉个缓动、做建模的捏个曲面背后全是贝塞尔和B样条在撑着。但很多人对它们的理解停留在拖控制点的层面一旦遇到需要自己实现曲线求值、判断曲率、做曲线拟合的场景就抓瞎了。这篇东西就是把我这些年在这两种曲线上踩过的坑、算过的公式、写过的代码从头到尾捋一遍。不管你是刚学计算机图形学的学生还是工作中需要手撸曲线算法的工程师应该都能从中找到点有用的东西。1. 从一条拉不直的曲线说起1.1 为什么需要参数曲线先说个最朴素的问题给你一堆点让你画一条光滑的曲线穿过它们你怎么画最直接的想法是多项式插值——找一个n-1次多项式让它精确穿过n个点。数学上完全可行但实际用起来是灾难。龙格现象了解一下高次多项式在区间边缘会产生剧烈振荡点越多振荡越厉害。你给10个点它给你画出一条上下翻飞的蛇。另一个思路是分段。每两个相邻点之间用一段低次曲线连起来但这样在连接处会出现折角——切线不连续看起来就是一段一段硬拼的。要让它光滑就得保证连接点处导数连续甚至二阶导数连续。这就是样条曲线的核心动机。参数曲线的好处在于它把x和y三维就是x、y、z都表示成独立参数t的函数P(t) (x(t), y(t), z(t))这样做的好处是曲线可以自交、可以竖直、可以闭合不受一个x对应一个y的限制。贝塞尔和B样条都是参数曲线的典型代表。1.2 贝塞尔和B样条的关系一句话说清很多人搞不清这俩的关系。我打个比方贝塞尔曲线像是全局投票制每个控制点对整条曲线都有影响改一个点整条曲线都跟着动B样条像是局部代表制每个控制点只影响它附近的一小段改一个点远处纹丝不动。这个区别直接决定了两者的适用场景。做设计工具、做字体轮廓贝塞尔够用且直观做工业建模、做复杂曲面拼接B样条才是正解。NURBS非均匀有理B样条则是在B样条基础上引入了权重能精确表示圆、椭圆这些二次曲线是CAD领域的绝对主力。2. 贝塞尔曲线直觉最好但暗藏陷阱2.1 德卡斯特里奥算法递归就是一切贝塞尔曲线的定义有好几种写法但真正适合写代码的是德卡斯特里奥de Casteljau算法。它的思想极其优雅在控制多边形上反复做线性插值最后收敛到曲线上的一个点。给定n1个控制点P0到Pn要算参数t处的点对每一对相邻点做插值 P_i^(1) (1-t) * P_i t * P_(i1) 然后对新的点集重复这个过程直到只剩一个点。用Python写出来大概长这样def de_casteljau(points, t): pts [p[:] for p in points] # 复制一份别改原数据 n len(pts) for r in range(1, n): for i in range(n - r): pts[i] [ (1 - t) * pts[i][0] t * pts[i 1][0], (1 - t) * pts[i][1] t * pts[i 1][1] ] return pts[0]这个算法的数值稳定性非常好因为它只做凸组合系数非负且和为1不会出现除零或者溢出。相比之下直接用伯恩斯坦基函数展开求值在控制点很多的时候数值误差会明显增大。提示如果你只是要算曲线上的点永远优先用德卡斯特里奥别去展开多项式。这是数值计算的基本素养。2.2 伯恩斯坦基函数理解贝塞尔性质的钥匙德卡斯特里奥算法好写但要理解贝塞尔曲线的性质得看伯恩斯坦基函数。n次贝塞尔曲线可以写成B(t) Σ C(n,i) * t^i * (1-t)^(n-i) * P_i其中C(n,i)是组合数。这些基函数有几个关键性质非负性每项都≥0单位分解所有基函数加起来恒等于1对称性B_i,n(t) B_(n-i),n(1-t)单位分解这个性质直接导致了贝塞尔曲线的一个重要特性——凸包性。曲线永远落在控制点构成的凸包内部。这个性质在碰撞检测、裁剪算法里非常有用如果凸包不相交曲线肯定不相交可以快速排除。但伯恩斯坦基函数是全局的每个基函数在整个[0,1]区间上都非零。这就是为什么贝塞尔曲线牵一发动全身——移动一个控制点所有基函数的值都变了整条曲线都受影响。2.3 曲率怎么算从一阶导到二阶导热词里有人问贝塞尔曲线有没有判断曲率的方法答案是当然有而且公式很标准。参数曲线的曲率公式是κ |xy - yx| / (x^2 y^2)^(3/2)对于贝塞尔曲线导数可以通过对控制点做差分来算。n次贝塞尔曲线的一阶导是n次实际是n-1次贝塞尔曲线控制点是n*(P_(i1) - P_i)。二阶导同理是对一阶导的控制点再做一次差分。def bezier_derivatives(points): n len(points) - 1 # 一阶导控制点 d1 [[n * (points[i1][0] - points[i][0]), n * (points[i1][1] - points[i][1])] for i in range(n)] # 二阶导控制点 d2 [[(n-1) * (d1[i1][0] - d1[i][0]), (n-1) * (d1[i1][1] - d1[i][1])] for i in range(n-1)] return d1, d2然后在参数t处分别求值代入曲率公式即可。这里有个坑在曲线端点处如果一阶导为零比如控制点重合曲率公式会除零。实际代码里要加保护或者用极限的方式处理。曲率在工程上有什么用最典型的是道路设计——曲率不能突变否则车辆行驶会不舒服还有字体渲染——曲率大的地方需要更密的采样点否则看起来会有棱角。2.4 贝塞尔的致命短板局部修改的代价假设你有一条20个控制点的贝塞尔曲线现在只想把中间某一段稍微调整一下。对不起做不到。你动任何一个控制点整条曲线都会变。这在交互式设计里是灾难性的。更麻烦的是贝塞尔曲线的次数等于控制点数减一。控制点一多次数就高高次贝塞尔曲线不仅计算量大而且数值稳定性差还容易出现控制点明明很规整曲线却扭来扭去的情况。这两个问题B样条都解决了。3. B样条工业级曲线的主力3.1 节点向量B样条的时间轴B样条最让人困惑的概念就是节点向量knot vector。我换个说法它定义了每个控制点什么时候上场、什么时候下场。一个p次B样条节点向量是一个非递减的实数序列U {u_0, u_1, ..., u_m}其中m n p 1n是控制点数减一。节点向量的长度 控制点数 次数 1。节点向量分几种类型类型特点用途均匀节点等间距理论分析简单场景准均匀两端重复p1次曲线过端点常用非均匀任意非递减最灵活NURBS基础分段贝塞尔内部节点重复p次等价于多段贝塞尔拼接准均匀节点向量是最常用的因为它保证曲线经过首末控制点同时内部又是均匀的。比如3次B样条4个控制点准均匀节点向量是{0,0,0,0,1,1,1,1}——这其实就是一条贝塞尔曲线。3.2 考克斯-德布尔递推B样条的核心算法B样条的基函数由考克斯-德布尔Cox-de Boor递推公式定义N_i,0(u) 1 if u_i u u_(i1), else 0 N_i,p(u) (u - u_i)/(u_(ip) - u_i) * N_i,p-1(u) (u_(ip1) - u)/(u_(ip1) - u_(i1)) * N_i1,p-1(u)这个递推看起来吓人但逻辑很清晰p次基函数是由两个p-1次基函数加权组合而成的。权重是两个线性函数分别对应从左边界靠近和从右边界远离。写代码的时候要注意分母为零的情况——当节点重复时分母可能为0。标准做法是约定0/0 0。def cox_de_boor(i, p, u, knots): if p 0: return 1.0 if knots[i] u knots[i1] else 0.0 left 0.0 right 0.0 denom1 knots[ip] - knots[i] if denom1 1e-12: left (u - knots[i]) / denom1 * cox_de_boor(i, p-1, u, knots) denom2 knots[ip1] - knots[i1] if denom2 1e-12: right (knots[ip1] - u) / denom2 * cox_de_boor(i1, p-1, u, knots) return left right这个递归实现直观但效率低实际工程中会用迭代版本或者预计算基函数表。3.3 局部支撑性B样条最值钱的性质B样条基函数有个关键性质叫局部支撑N_i,p(u)只在区间[u_i, u_(ip1))上非零。这意味着移动一个控制点只影响p1个节点区间内的曲线段。这个性质的价值怎么强调都不过分。它意味着修改局部不影响全局交互式编辑友好曲线次数可以独立于控制点数想加多少控制点都行次数不变数值计算可以局部化效率高举个例子一条3次B样条100个控制点。移动第50个控制点只有第47到第53个控制点对应的曲线段会变其余部分完全不动。这在贝塞尔曲线里是不可想象的。3.4 德布尔算法B样条版的德卡斯特里奥和贝塞尔有德卡斯特里奥一样B样条有德布尔de Boor算法。它的结构和德卡斯特里奥几乎一样只是插值参数不是固定的t而是根据节点向量来确定的。def de_boor(u, points, knots, p): # 找到u所在的节点区间 k find_span(u, points, knots, p) # 复制p1个控制点 d [points[k-pi][:] for i in range(p1)] for r in range(1, p1): for i in range(p, r-1, -1): alpha (u - knots[k-pi]) / (knots[k1i-r] - knots[k-pi]) d[i] [(1-alpha)*d[i-1][j] alpha*d[i][j] for j in range(len(d[i]))] return d[p]德布尔算法同样数值稳定是B样条求值的标准方法。find_span函数负责定位参数u落在哪个节点区间可以用二分查找加速。4. 从理论到代码手撸一个曲线库4.1 整体架构设计光看公式容易飘我把自己写的一个小曲线库的结构分享一下。核心就三个类BezierCurve贝塞尔曲线存控制点列表BSplineCurveB样条曲线存控制点、次数、节点向量CurveUtils工具函数采样、求导、曲率、弧长设计上有个取舍要不要把贝塞尔做成B样条的特例理论上可以分段贝塞尔就是内部节点重复p次的B样条但实际用起来贝塞尔有更简单的专用算法硬套B样条反而绕。所以我选择分开实现共享工具函数。4.2 采样与自适应细分把曲线画出来最直接的方法是均匀采样t从0到1每隔0.01取一个点连成折线。但这样有个问题曲率大的地方采样太稀看起来有棱角曲率小的地方采样太密浪费。更好的做法是自适应细分递归地把曲线段一分为二如果这一段足够平就停止否则继续分。判断平的标准可以用控制点到弦的最大距离。def adaptive_sample(curve, tolerance0.5, max_depth10): result [] def subdivide(t0, t1, depth): p0 curve.evaluate(t0) p1 curve.evaluate(t1) pm curve.evaluate((t0 t1) / 2) # 中点到弦的距离 d point_line_distance(pm, p0, p1) if d tolerance or depth max_depth: result.append((t0, p0)) result.append((t1, p1)) else: subdivide(t0, (t0t1)/2, depth1) subdivide((t0t1)/2, t1, depth1) subdivide(0, 1, 0) return result这个方法的采样密度自动适应曲率画出来的曲线既光滑又不会浪费顶点。实测下来tolerance取0.5像素左右视觉效果就很好了。4.3 弧长参数化匀速运动的关键如果你想让一个物体沿着曲线匀速运动直接用t做参数是不行的——t均匀变化时曲线上的实际速度不均匀。需要做弧长参数化。基本思路是先用数值积分算出曲线的总弧长然后建立t和弧长s的映射表运动时按s均匀推进再反查t。def build_arc_length_table(curve, samples1000): table [(0.0, 0.0)] prev curve.evaluate(0) total 0.0 for i in range(1, samples 1): t i / samples curr curve.evaluate(t) total distance(prev, curr) table.append((t, total)) prev curr return table, total def t_at_arc_length(table, total, s): target s * total # 二分查找 lo, hi 0, len(table) - 1 while lo hi - 1: mid (lo hi) // 2 if table[mid][1] target: lo mid else: hi mid # 线性插值 t0, s0 table[lo] t1, s1 table[hi] if s1 - s0 1e-12: return t0 return t0 (target - s0) / (s1 - s0) * (t1 - t0)这个表建一次可以反复用对于动画路径来说非常划算。注意采样数要够否则反查精度不够运动会有轻微抖动。4.4 曲线拟合从离散点到控制点实际工作中经常遇到反问题给一堆离散点求一条逼近它们的B样条曲线。这就是曲线拟合。最常用的方法是最小二乘拟合。步骤大致是对数据点做参数化均匀、弦长、向心都行确定节点向量通常用平均法或均匀法建立最小二乘方程组求解控制点弦长参数化是最常用的因为它考虑了点的实际分布u_0 0 u_i u_(i-1) |P_i - P_(i-1)| 最后归一化到[0,1]节点向量的选择对拟合质量影响很大。平均法averaging是比较稳妥的选择u_(jp) (1/p) * Σ u_(i) for i j to jp-1拟合的精度取决于控制点数量。控制点太少拟合不上去太多会过拟合曲线出现不必要的波动。实践中一般从较少控制点开始逐步增加直到误差满足要求。5. 那些文档里不会写的坑5.1 节点向量的端点重复问题写B样条代码最容易翻车的地方就是节点向量。准均匀节点向量要求两端各重复p1次比如3次B样条节点向量开头必须是{0,0,0,0}结尾是{1,1,1,1}。少一个重复曲线就不过端点多一个分母就出零。我见过有人把节点向量写成{0,0,0,1,1,1}然后纳闷为什么曲线不经过第一个控制点。原因就是3次B样条需要4个重复他只写了3个。注意节点向量的长度必须等于控制点数 次数 1。这个关系式是硬约束写代码时最好加个断言检查。5.2 曲率计算的数值陷阱前面提过曲率公式在导数为零时会炸。但实际中还有更隐蔽的问题当曲线接近直线时二阶导很小曲率计算会有很大的相对误差。我的做法是加一个阈值判断如果一阶导的模长小于某个epsilon直接返回曲率为零如果二阶导也很小同样返回零。虽然理论上不严谨但工程上足够用而且避免了NaN传播。另一个坑是曲率符号。二维曲线的曲率有正负之分表示弯曲方向。公式里的绝对值去掉分子xy - yx的符号就是曲率符号。做道路设计或者字体渲染时这个符号很重要。5.3 高次B样条的数值稳定性虽然B样条理论上次数可以任意高但实践中很少超过5次。原因和高次多项式一样数值稳定性下降而且高次基函数的支撑区间很宽局部性变差。工业界最常用的是3次B样条C2连续偶尔用4次或5次。如果你发现自己需要10次B样条大概率是建模方式有问题应该考虑用多段低次曲线拼接。5.4 重节点与连续性节点重复次数直接影响曲线的连续性。一个p次B样条如果某个内部节点的重复度是r那么曲线在该节点处的连续性降为C^(p-r)。重复度1C^(p-1)连续标准情况重复度pC^0连续曲线在该点有折角重复度p1曲线断开这个性质在建模时非常有用。想要一条有尖角的曲线把某个节点重复p次就行。想要曲线断开成两段重复p1次。5.5 浮点数比较的坑节点区间的查找、参数是否在[0,1]范围内、两个点是否重合——这些判断在浮点数下都不能用。我一般用1e-10作为epsilon但具体值要根据数据尺度调整。数据范围是毫米级和公里级epsilon肯定不一样。一个实用的技巧是所有参数化都归一化到[0,1]这样epsilon可以固定用1e-10不用每次调整。6. NURBSB样条的完全体6.1 权重到底在干什么NURBS在B样条基础上给每个控制点加了一个权重w_i。基函数变成R_i,p(u) N_i,p(u) * w_i / Σ N_j,p(u) * w_j权重的物理意义是吸引力权重越大曲线越靠近该控制点。当所有权重相等时NURBS退化为普通B样条。权重最强大的地方在于它能精确表示圆锥曲线。一个标准的圆用普通B样条只能逼近用NURBS可以精确表示。这是NURBS成为CAD标准的核心原因。6.2 齐次坐标NURBS的实现技巧直接按上面的公式实现NURBS每次求值都要算分母效率低。标准做法是用齐次坐标把控制点升维(x, y, z) - (wx, wy, wz, w)在四维空间做普通B样条求值把结果投影回三维(X, Y, Z, W) - (X/W, Y/W, Z/W)这样所有B样条的算法德布尔、求导、细分都能直接复用代码量大大减少。这个技巧在OpenGL的NURBS求值器里也是这么用的。6.3 NURBS的求导NURBS的导数不能直接套B样条的导数公式因为分母里有权重。标准做法是用商法则C(u) A(u) / w(u) C(u) (A(u) - w(u) * C(u)) / w(u)其中A(u)是齐次坐标下的B样条曲线w(u)是权重基函数的和。高阶导数可以递推。实现的时候先算齐次空间里的导数再转换回来。7. 实际项目中的选型建议7.1 什么时候用贝塞尔什么时候用B样条这个问题我被问过很多次。我的经验是场景推荐理由UI圆角、图标贝塞尔控制点少直观浏览器原生支持字体轮廓贝塞尔二次/三次行业标准渲染器优化好动画缓动贝塞尔cubic-bezier就够用工业建模NURBS精度要求高需要精确圆锥曲线路径规划B样条需要局部修改C2连续保证平滑数据拟合B样条控制点数灵活局部支撑简单说控制点少、需要直观交互用贝塞尔控制点多、需要局部修改、需要高阶连续用B样条。7.2 次数选择的经验法则2次C1连续够用但不够光滑适合快速原型3次C2连续工业标准99%的场景选它4次及以上特殊需求比如需要C3连续的气动外形3次是甜点因为C2连续对人眼来说已经足够光滑而且计算量适中。再高次数视觉上提升不明显计算和数值稳定性都变差。7.3 控制点数量的权衡控制点越多曲线越灵活但计算量线性增长过拟合风险增加交互编辑变复杂我的经验是先用尽量少的控制点看拟合误差不够再加。对于拟合问题一般从数据点数的1/3到1/2开始试。8. 调试曲线的实用手段8.1 可视化控制多边形和基函数调试曲线问题第一件事是把控制多边形画出来。曲线一定在控制多边形的凸包内如果曲线跑到外面去了肯定是代码有bug。第二件事是把基函数画出来。B样条的基函数图能直观显示每个控制点的影响范围。如果基函数图看起来不对比如有负值、不连续、支撑区间不对那节点向量肯定有问题。8.2 用已知特例验证写曲线库的时候我习惯用几个已知特例做单元测试所有控制点重合曲线应该退化成一个点控制点共线曲线应该是一条直线段准均匀节点向量控制点数次数1应该等价于贝塞尔曲线对称的控制点曲线应该对称这些测试能抓住大部分实现错误。8.3 数值验证导数用差分校验解析求导容易写错用数值差分校验是个好办法f(x) ≈ (f(xh) - f(x-h)) / (2h)h取1e-6左右比较解析导数和差分导数的相对误差。如果误差大于1e-3基本可以确定解析求导有问题。这个方法我用来抓过好几次符号错误。8.4 曲率的可视化曲率不好直接看但可以画曲率梳curvature comb在每个采样点处沿法线方向画一条长度正比于曲率的线段。曲率梳光滑说明曲线质量好曲率梳有尖刺说明曲率不连续曲率梳交叉说明有拐点。这个工具在汽车外形设计里是标配做任何高质量曲线都应该用。9. 性能优化从能用 to 好用9.1 预计算基函数表如果曲线要反复求值比如动画每帧都要算预计算基函数表能大幅提速。把[0,1]分成N份每份预先算好所有非零基函数的值求值时直接查表加插值。代价是内存和初始化时间但对于实时应用来说完全值得。N取1000左右精度和速度的平衡就很好。9.2 利用局部支撑做裁剪B样条的局部支撑意味着对于给定的参数u只有p1个基函数非零。求值时只需要遍历这p1个控制点而不是全部。控制点越多这个优化的收益越大。实现上先用二分查找定位u所在的节点区间然后只取对应的p1个控制点参与计算。9.3 SIMD与并行曲线求值天然适合并行不同参数t之间没有依赖。用SIMD指令一次算4个或8个点或者用多线程分摊采样任务都能获得接近线性的加速。不过要注意德布尔算法内部有依赖每一轮依赖上一轮的结果所以并行应该在参数维度上做而不是在算法内部做。10. 写在最后的一些个人体会曲线这东西公式看着多但核心思想就那么几个用基函数加权控制点用节点向量控制基函数的形状用递推算法保证数值稳定。把这三个东西吃透贝塞尔、B样条、NURBS就是同一套框架的不同配置。我刚开始学的时候最大的误区是死记公式。后来发现真正有用的是理解每个性质背后的几何直觉——为什么凸包性成立为什么局部支撑重要为什么重节点会降连续性。这些想明白了公式忘了都能推出来。另一个体会是一定要自己动手实现一遍。看十遍德布尔算法不如自己写一遍然后调通。调通过程中遇到的那些坑——节点向量长度不对、分母为零、端点不过点——才是真正长在身上的知识。最后说个实用的如果你只是要用曲线不想自己实现Python里scipy的interpolate模块有BSpline类JavaScript里d3有曲线生成器C里OpenCASCADE有完整的NURBS实现。但如果你要理解曲线还是得自己撸一遍。