一、NURBS曲线求值算法：从de Boor到曲率计算
各位同学好，我是老张。今天咱们聊聊NURBS曲线求值这件事。
说实话，我刚入行那会儿，看到NURBS公式就头疼。一堆控制点、权因子、节点向量，光记符号就够呛。但后来我发现，只要抓住de Boor算法这根主线，一切都清晰了。
你想想看，一条NURBS曲线，本质上就是分段多项式。de Boor算法就是用来高效计算这些分段多项式上的点、导数、曲率的。我习惯把它叫做「递推求值法」，因为每一步都依赖上一步的结果。
1.1 de Boor算法原理
de Boor算法是B样条曲线求值的标准方法。它基于一个核心思想：通过线性插值，逐步降低曲线的次数。
举个例子。一条3次B样条曲线，每个点由4个控制点决定。de Boor算法会先做一次插值，得到3个中间点；再插值一次，得到2个点；最后一次插值，得到最终曲线上的点。
这个过程，说白了就是「剥洋葱」。一层层剥下去，直到露出核心。
核心公式：
对于k次B样条曲线，给定节点向量U = [u₀, u₁, ..., u_{n+k+1}]，控制点P₀, P₁, ..., Pₙ，计算参数u处的点：
// 初始化
d_j = P_{j} (j = i-k, i-k+1, ..., i)
// 递推
for r = 1 to k:
for j = i-k+r to i:
α = (u - u_j) / (u_{j+k-r+1} - u_j)
d_j = (1-α) * d_{j-1} + α * d_j
// 结果
C(u) = d_i这里有个关键点：i是u所在的节点区间索引。怎么找？二分查找就行。我当年第一次实现时，就在这卡了三天——节点区间找错了，后面全白搭。
1.2 曲线点计算实现
好了，理论讲完，咱们直接上代码。这是我个人习惯的写法，清晰、可读性强。
def deBoor_point(knots, ctrl_pts, degree, u):
"""
de Boor算法计算曲线上的点
参数:
knots: 节点向量 (list)
ctrl_pts: 控制点 (list of (x,y,z))
degree: 曲线次数
u: 参数值
返回:
曲线上的点 (x,y,z)
"""
# 1. 找到u所在的节点区间
n = len(ctrl_pts) - 1
m = len(knots) - 1
# 二分查找
i = degree
while i < m - degree - 1 and u >= knots[i+1]:
i += 1
# 2. 初始化控制点副本
d = [ctrl_pts[j] for j in range(i-degree, i+1)]
# 3. 递推计算
for r in range(1, degree+1):
for j in range(degree, r-1, -1):
alpha = (u - knots[i-degree+j]) / (knots[i+j-r+1] - knots[i-degree+j])
d[j] = [(1-alpha)*d[j-1][k] + alpha*d[j][k] for k in range(3)]
return d[degree]避坑指南：
我曾经在递推循环里把j的遍历方向搞反了。从大到小还是从小到大？答案是从大到小。因为每次计算要用到上一轮的d[j-1]，如果从小到大更新，d[j-1]已经被覆盖了，结果全错。
1.3 曲线导数计算
光算点不够，很多时候我们需要导数。比如刀具路径规划，就得知道曲线在每一点的方向。
de Boor算法求导的思路很巧妙：对控制点做差分，然后降低一次次数。
具体来说，k次B样条的导数是k-1次B样条。新的控制点由原控制点的差分得到：
def deBoor_derivative(knots, ctrl_pts, degree, u):
"""
计算曲线在u处的一阶导数
"""
n = len(ctrl_pts) - 1
m = len(knots) - 1
# 找到节点区间
i = degree
while i < m - degree - 1 and u >= knots[i+1]:
i += 1
# 计算差分控制点
d_deriv = []
for j in range(n):
denom = knots[j+degree+1] - knots[j+1]
if denom != 0:
diff = [(ctrl_pts[j+1][k] - ctrl_pts[j][k]) * degree / denom
for k in range(3)]
d_deriv.append(diff)
else:
d_deriv.append([0, 0, 0])
# 对降次后的曲线用de Boor求值
return deBoor_point(knots[1:-1], d_deriv, degree-1, u)嗯，这里要注意：节点向量也要跟着变。去掉首尾节点，因为次数降低了，节点向量长度也变了。我刚开始做的时候忘了这茬，导数算出来总是不对。
1.4 曲线曲率计算
曲率，说白了就是曲线弯曲的程度。数控加工里，曲率直接决定了进给速度——曲率大的地方，速度得降下来，不然加工精度跟不上。
曲率公式：
κ = |C'(u) × C''(u)| / |C'(u)|³
需要一阶导和二阶导。二阶导怎么算？对一阶导再求一次导就行。
def curvature(knots, ctrl_pts, degree, u):
"""
计算曲线在u处的曲率
"""
# 一阶导
d1 = deBoor_derivative(knots, ctrl_pts, degree, u)
# 二阶导：对一阶导的结果再求导
# 注意：一阶导的曲线次数是degree-1
d2 = deBoor_derivative(knots[1:-1],
[deBoor_point(knots, ctrl_pts, degree, v)
for v in [u-0.001, u, u+0.001]],
degree-1, u)
# 叉积模长
cross = [
d1[1]*d2[2] - d1[2]*d2[1],
d1[2]*d2[0] - d1[0]*d2[2],
d1[0]*d2[1] - d1[1]*d2[0]
]
cross_norm = (cross[0]**2 + cross[1]**2 + cross[2]**2)**0.5
# 一阶导模长
d1_norm = (d1[0]**2 + d1[1]**2 + d1[2]**2)**1.5
if d1_norm == 0:
return 0
return cross_norm / d1_norm重要提醒：
曲率计算对数值精度很敏感。我建议用解析法求二阶导，而不是数值差分。上面代码为了演示用了数值法，实际项目中请用解析法——也就是对de Boor算法再求一次导。
1.5 知识体系总览
说了这么多，咱们用一张图把整个知识体系串起来：
1.6 实战经验总结
最后，分享几个我这些年踩过的坑：
- 节点向量重复度：首尾节点重复度必须是k+1，否则曲线端点不通过控制点。我见过有人设成k，结果曲线端点飘了。
- 数值稳定性：当节点区间很小时，分母可能接近0。加个epsilon判断，或者用双精度浮点。
- 性能优化：如果连续计算大量点，可以缓存节点区间查找结果。我做过测试，能快30%以上。
- 边界处理：u在节点向量边界时，要特殊处理。比如u等于最后一个节点，直接返回最后一个控制点。
我的建议：
刚开始学de Boor算法，别急着看公式。拿张纸，画几个控制点，手动算一遍。我当年就是这么干的，算完一遍，所有细节都刻在脑子里了。
好了，这一章就到这里。de Boor算法是NURBS的基石，搞懂了它，后面的NURBS曲线插补、曲面求值都会轻松很多。
公众号:蓝海资料掘金营，微信deep3321