B样条曲线的deboor算法

 众所周知，B样条曲线的通用表达式为：
$$C(u)=\sum _{i=0}^n{N}_{i,p}(u)Q_i$$
其中$N_{i,p}(u)$为基函数，$Q_i$为控制点。且基函数由下式递推算出：

 在实际使用中，如果已知所有的控制点$Q_i$，我们希望求出参数u处的B样条值来。这种给定参数u,计算参数曲线上对应点的运算称为求值(Evaluate)操作。为了方便后续讨论，我们这里设B样条曲线有n+1个控制点，阶数为p。
 由B样条曲线的表达式可以看出来，最简单的Evaluate方法就是先由$N_{0,0},\cdots ,N_{n+p,0}$递推地求出基函数$N_{0,p},\cdots ,N_{n,p}$，然后再用公式$C(u_j)={\sum }_{i=0}^n{N}_{i,p}(u)Q_i$把$C(u)$的值算出来，这个计算过程可以用下图理解，计算复杂度为O(n*(n+n+1))=O(np)。程序可以参考B样条值的公式计算方法。

 由公式直接算值虽然简单直观，但是程序复杂度太高，所以有人又提出了一种更快速的Evaluate方法-deboor计算方法。推导过程如下((假设我们已经提前找到了要求的参数$u\in [u_j,u_{j+1}$)：
$$\begin{array}{c}\begin{array}{rl}C(u)& =\sum _{i=0}^n{N}_{i,p}(u)Q_i\\ & =\sum _{i=j-p}^j{N}_{i,p}(u)Q_i\\ & =\sum _{i=j-p}^j(\frac{u-u_i}{u_{i+p}-u_i}{N}_{i,p-1}(u)+\frac{u_{i+p+1}-u}{u_{i+p+1}-u_{i+1}}{N}_{i+1,p-1}(u))Q_i\\ & =\sum _{i=j-p}^j\frac{u-u_i}{u_{i+p}-u_i}{N}_{i,p-1}(u)Q_i+\sum _{i=j-p}^j\frac{u_{i+p+1}-u}{u_{i+p+1}-u_{i+1}}{N}_{i+1,p-1}(u)Q_i\\ & =\sum _{i=j-p+1}^j\frac{u-u_i}{u_{i+p}-u_i}{N}_{i,p-1}(u)Q_i+\frac{u-u_{j-p}}{u_{j-p+p}-u_{j-p}}{N}_{j-p,p-1}(u)Q_j\\ & +\sum _{i=j-p+1}^j\frac{u_{i+p}-u}{u_{i+p}-u_i}{N}_{i,p-1}(u)Q_{i-1}+\frac{u_{j+p+1}-u}{u_{j+p+1}-u_{j+1}}{N}_{j+1,p-1}(u)Q_j\\ & =\sum _{i=j-p+1}^j\frac{u-u_i}{u_{i+p}-u_i}{N}_{i,p-1}(u)Q_i+\sum _{i=j-p+1}^j\frac{u_{i+p}-u}{u_{i+p}-u_i}{N}_{i,p-1}(u)Q_{i-1}\\ & (因为N_{j-p,p-1}(u)只在[u_{j-p},u_{j-p+p-1+1})不为0，\\ & N_{j+1,p-1}(u)只在[u_{j+1},u_{j+1+p-1+1})不为0，而u_j\in [u_j,u_{j+1}))\\ & =\sum _{i=j-p+1}^j{N}_{i,p-1}(u)(\frac{u-u_i}{u_{i+p}-u_i}{Q}_i+\frac{u_{i+p}-u}{u_{i+p}-u_i}{Q}_{i-1})\end{array}\end{array}$$
令：
$$\begin{array}{c}{Q}_{i,k}(u)=\{\begin{array}{rlr}& Q_i& k=0\\ & \frac{u-u_i}{u_{i+p+1-k}-u_i}{Q}_{i,k-1}(u)+\frac{u_{i+p+1-k}-u}{u_{i+p+1-k}-u_i}{Q}_{i-1,k-1}(u)& k=1,\cdots ,p\end{array}\end{array}$$
则上式可以写为：
$$\begin{array}{c}\begin{array}{rl}C(u)& =\sum _{i=j-p+1}^j{N}_{i,p-1}(u)Q_{i,1}(u)\\ & =\sum _{i=j-p+1}^j(\frac{u-u_i}{u_{i+p-1}-u_i}{N}_{i,p-2}(u)+\frac{u_{i+p}-u}{u_{i+p}-u_{i+1}}{N}_{i+1,p-2}(u))Q_{i,1}(u)\\ & =\sum _{i=j-p+1}^j\frac{u-u_i}{u_{i+p-1}-u_i}{N}_{i,p-2}(u)Q_{i,1}(u)+\sum _{i=j-p+1}^j\frac{u_{i+p}-u}{u_{i+p}-u_{i+1}}{N}_{i+1,p-2}(u))Q_{i,1}(u)\\ & =\sum _{i=j-p+2}^j\frac{u-u_i}{u_{i+p-1}-u_i}{N}_{i,p-2}(u)Q_{i,1}(u)\\ & +\frac{u-u_{j-p+1}}{u_{j-p+1+p-1}-u_{j-p+1}}{N}_{j-p+1,p-2}(u)Q_{j-p+1,1}(u)\\ & +\sum _{i=j-p+2}^j\frac{u_{i+p-1}-u}{u_{i+p-1}-u_i}{N}_{i,p-2}(u)Q_{i-1,1}(u)\\ & +\frac{u_{j+p}-u}{u_{j+p}-u_{j+1}}{N}_{j+1,p-2}(u)Q_{j,1}(u)\\ & =\sum _{i=j-p+2}^j\frac{u-u_i}{u_{i+p-1}-u_i}{N}_{i,p-2}(u)Q_{i,1}(u)\\ & +\sum _{i=j-p+2}^j\frac{u_{i+p-1}-u}{u_{i+p-1}-u_i}{N}_{i,p-2}(u)Q_{i-1,1}(u)\\ & (因为N_{j-p+1,p-2}(u)只在[u_{j-p+1},u_{j-p+1+p-2+1})不为0，\\ & N_{j+1,p-2}(u)只在[u_{j+1},u_{j+1+p-2+1})不为0，而u_j\in [u_j,u_{j+1}))\\ & =\sum _{i=j-p+2}^j{N}_{i,p-2}(u)(\frac{u-u_i}{u_{i+p-1}-u_i}{Q}_{i,1}(u)+\frac{u_{i+p-1}-u}{u_{i+p-1}-u_i}{Q}_{i-1,1}(u))\\ & =\sum _{i=j-p+2}^j{N}_{i,p-2}(u)Q_{i,2}(u)\\ & =\cdots \\ & =\sum _{i=j-p+k}^j{N}_{i,p-k}(u)Q_{i,k}(u)\\ & =\cdots \\ & =\sum _{i=j-p+p}^j{N}_{i,p-p}(u)Q_{i,p}(u)\\ & =N_{j,0}(u)Q_{j,p}(u)\\ & =Q_{j,p}(u)\end{array}\end{array}$$
 综上所述，$C(u)=Q_{j,p}(u),u\in [u_j,u_{j+1})$，我们只需要算出$Q_{j,p}(u)$就可以了，而$Q_{j,p}(u)$的递推计算公式我们已经在上面给出了，这里再写一遍：
$$\begin{array}{c}{Q}_{i,k}(u)=\{\begin{array}{rlr}& Q_i& k=0\\ & \frac{u-u_i}{u_{i+p+1-k}-u_i}{Q}_{i,k-1}(u)+\frac{u_{i+p+1-k}-u}{u_{i+p+1-k}-u_i}{Q}_{i-1,k-1}(u)& k=1,\cdots ,p\end{array}\end{array}$$
递推计算方法如下图所示，复杂度为$O(p^2)$，代码参考EvaluateDeboor。

（由于网站要梯子，所以这里附EvaluateDeboor的代码：）

def deBoor(k: int, x: int, t, c, p: int):
    """Evaluates S(x).

    Arguments
    ---------
    k: Index of knot interval that contains x.
    x: Position.
    t: Array of knot positions, needs to be padded as described above.
    c: Array of control points.
    p: Degree of B-spline.
    """
    d = [c[j + k - p] for j in range(0, p + 1)] 

    for r in range(1, p + 1):
        for j in range(p, r - 1, -1):
            alpha = (x - t[j + k - p]) / (t[j + 1 + k - r] - t[j + k - p]) 
            d[j] = (1.0 - alpha) * d[j - 1] + alpha * d[j]

    return d[p]

参考文献：

1. https://en.wikipedia.org/wiki/De_Boor%27s_algorithm
2. http://www.whudj.cn/?p=535