4. 基函数计算：德布尔-考克斯递推公式的推导与编程实现

好，咱们今天聊点硬核的——NURBS 的基函数计算。

说实话，我刚入行那会儿，看到这个递推公式就头疼。一堆下标、一堆阶数，感觉像在看天书。但后来我发现，这东西其实没那么玄乎。你只要抓住一个核心思想：递推。

说白了，高阶的基函数，是由低阶的基函数一层一层“搭”出来的。就像搭积木，你不可能直接搭出顶层，得从底层开始。

4.1 为什么需要递推？

你想想看，NURBS 曲线要光滑，基函数就得有连续性。而 B 样条的基函数，天然就支持这种连续性。但问题是，你不能直接用一个公式算出任意阶数的基函数值。

我个人的习惯是，把这个问题拆成两步：

1. 找到支撑区间：这个基函数在哪些参数 u 的范围内是非零的？
2. 递推计算：在这个区间内，用低阶的值算出高阶的值。

嗯，这里要注意：支撑区间这个概念特别重要。它决定了你的计算范围，也决定了曲线的局部性。

4.2 德布尔-考克斯递推公式

公式长这样，别怕，咱们一行一行拆解：

// 0 阶基函数（最简单）
if (p == 0) {
    if (u_i <= u && u < u_{i+1})
        N[i][0] = 1.0;
    else
        N[i][0] = 0.0;
}
// 高阶基函数（递推）
else {
    double left  = (u - u_i)   / (u_{i+p}   - u_i);
    double right = (u_{i+p+1} - u) / (u_{i+p+1} - u_{i+1});
    N[i][p] = left * N[i][p-1] + right * N[i+1][p-1];
}

你看，0 阶的时候就是个“开关”——u 落在节点区间里就是 1，否则是 0。到了高阶，就是两个低阶基函数的加权和。

我在项目中遇到过一个问题：当分母为零时怎么办？比如 u_i+p == u_i，这时候 left 项的分母是 0。按照数学定义，这种情况我们约定 0/0 = 0。代码里要加个判断：

double left = 0.0;
if (u_{i+p} != u_i) {
    left = (u - u_i) / (u_{i+p} - u_i);
}

我曾经因为没处理这个边界情况，导致插补器在节点重复的地方直接崩溃。嗯，从那以后我再也不敢偷懒了。

4.3 理解支撑区间

支撑区间这个概念，说白了就是“这个基函数能影响到哪些地方”。

举个例子，假设节点向量是 [0, 0, 0, 1, 2, 3, 4, 4, 4]，p=2（二次样条）。那么：

• N_0,2(u) 的支撑区间是 [u₀, u₃) = [0, 1)
• N_1,2(u) 的支撑区间是 [u₁, u₄) = [0, 2)
• N_2,2(u) 的支撑区间是 [u₂, u₅) = [0, 3)

你发现规律了吗？每个基函数只覆盖一小段 u 范围。这就是 NURBS 的局部支撑性——移动一个控制点，只影响曲线的一部分，不会全局变形。

4.4 编程实现要点

写代码的时候，我建议你注意这几点：

1. 二维数组存储：用 N[i][p] 表示第 i 个 p 阶基函数的值。先算 p=0 的所有值，再算 p=1，以此类推。
2. 只计算非零区间：对于给定的 u，只有少数几个基函数是非零的。你可以通过二分查找找到 u 所在的节点区间，然后只计算相关的基函数。
3. 边界处理：在节点向量的两端，通常会有重复节点（比如 [0,0,0,...]）。这时候要特别注意分母为零的情况。

// 计算所有非零基函数
void computeBasisFunctions(double u, int p, const vector<double>& knots, vector<double>& N) {
    // 找到 u 所在的节点区间
    int span = findSpan(u, p, knots);
    
    // 初始化 0 阶基函数
    N[0] = 1.0;
    
    // 递推计算高阶
    for (int k = 1; k <= p; k++) {
        double saved = 0.0;
        for (int r = 0; r < k; r++) {
            double left  = (u - knots[span - r]) / (knots[span + k - r] - knots[span - r]);
            double right = (knots[span + k + 1] - u) / (knots[span + k + 1] - knots[span - r + 1]);
            
            double temp = N[r];
            N[r] = saved + left * temp;
            saved = right * temp;
        }
        N[k] = saved;
    }
}

这段代码是优化过的版本，只计算当前 u 相关的基函数，效率比全量计算高很多。我在做实时插补时用的就是这种实现。

4.5 知识体系图

下面这张图展示了基函数计算的核心逻辑：

4.6 避坑指南

好了，基函数这块就聊到这儿。你只要把递推公式写对，支撑区间搞清楚，后面算曲线点、算导数就水到渠成了。

---