B 样条曲线

一、B 样条曲线的定义和性质

$p$ 次 B 样条曲线的定义为
$$C(u)=\sum _{i=0}^n{N}_{i,p}(u)P_i,a\le u\le b\text{\hspace{1em}(1)}$$

这里 $P_i$ 是曲线的控制点，{Ni,p(u)}\{N_{i,p}(u)\}{Ni,p​(u)} 是定义在非周期（并且非均匀）节点矢量
$$U=(\underset{p+1}{\underbrace{a,\cdots ,a}},u_{p+1},\cdots ,u_{m-p-1},\underset{p+1}{\underbrace{b,\cdots ,b}})$$

（包含 $m+1$ 个节点）上的 $p$ 次 B 样条基函数。除非特别声明，通常取 $a=0,b=1$。由 $P_i$ 构成的多边形称为控制多边形。

对于固定的 $u$ 值，计算 B 样条曲线上的对应点需要三步：

1. 找到 $u$ 所在的节点区间；
2. 计算非零的基函数；
3. 将非零基函数的值与相应的控制点相乘，再求和。

例1. U={0,0,0,1,2,3,4,4,5,5,5},  u=52,  p=2\quad U=\{0,0,0,1,2,3,4,4,5,5,5\},\ \ u=\dfrac{5}{2},\ \ p=2U={0,0,0,1,2,3,4,4,5,5,5},  u=25​,  p=2
此时，$u\in [u_4,u_5)$，则
$$N_{2,2}\left(\frac{5}{2}\right)=\frac{1}{8},N_{3,2}\left(\frac{5}{2}\right)=\frac{6}{8},N_{4,2}\left(\frac{5}{2}\right)=\frac{1}{8}$$

与相应的控制点相乘并求和得
$$C\left(\frac{5}{2}\right)=\frac{1}{8}{P}_2+\frac{6}{8}{P}_3+\frac{1}{8}{P}_4$$

具体算法如下：

def CurvePoint(p, U, P, u):
    """
    计算 B 样条曲线上的点
    :param p: 基函数次数
    :param U: 节点矢量
    :param P: 控制点
    :param u: 固定值
    :return: 固定值 u 在 B 样条曲线上的对应值
    """
    span = FindSpan(p, u, U)
    N = BasisFuns(span, u, p, U)
    C = 0.0
    for i in range(p+1):
        C = C + N[i] * P[span - p + i]
    return C

二、B 样条曲线的导矢

令 $C^{(k)}(u)$ 表示 $C(u)$ 的 $k$ 阶导矢。对于固定的 $u$，我们可以通过计算基函数的 $k$ 阶导数得到 $C^{(k)}(u)$。特别地，有
$$C^{(k)}(u)=\sum _{i=0}^n{N}_{i,p}^{(k)}(u)P_i$$

例2. p=2,U={0,0,0,1,2,3,4,4,5,5,5},u=52\quad p=2,U=\{0,0,0,1,2,3,4,4,5,5,5\},u=\frac{5}{2}p=2,U={0,0,0,1,2,3,4,4,5,5,5},u=25​
$$\begin{array}{rl}& N_{2,2}^{\prime }\left(\frac{5}{2}\right)=0-\frac{2}{3-1}\frac{1}{2}=-\frac{1}{2}\\ & N_{3,2}^{\prime }\left(\frac{5}{2}\right)=\frac{2}{3-1}\frac{1}{2}-\frac{2}{4-2}\frac{1}{2}=0\\ & N_{4,2}^{\prime }\left(\frac{5}{2}\right)=\frac{2}{4-2}\frac{1}{2}-0=\frac{1}{2}\end{array}$$

于是有
$$C^{\prime }\left(\frac{5}{2}\right)=-\frac{1}{2}{P}_2+\frac{1}{2}{P}_4$$

以下给出对于固定的 $u$，计算 B 样条曲线上对应的点及其直到 $d$ 阶（包括 $d$ 阶）导矢的算法。我们允许 $d>p$，虽然（对非有理曲线而言）这些导矢都为零；但是，这些导矢对有理曲线是必要的。算法的输入是 $u,d$ 以及下面的一些数据定义的 B 样条曲线

$n$：控制点的最大下标（控制点的个数是 $n+1$）;
$p$：曲线的次数；
$U$：节点矢量；
$P$：控制点数组。

算法的输出是数组 CK

def CurveDerivsAlg(n, p, U, P, u, d):
    """
    计算固定的 u 时 B 样条曲线上对应的点及其 d 阶导矢
    :param n: 控制点的最大下标(控制点的个数是 n+1)
    :param p: 曲线的次数
    :param U: 节点矢量
    :param P: 控制点数组
    :param u: 固定点
    :param d: 导矢阶数
    :return: 数组 CK
    """
    du = min(d, p)
    CK = np.zeros(d+1)
    for k in range(p+1, d+1):
        CK[k] = 0.0
    span = FindSpan(p, u, U)
    nders = DersBasisFuns(span, u, p, n, U)
    for k in range(du+1):
        CK[k] = 0.0
        for j in range(p+1):
            CK[k] = CK[k] + nders[k, j] * P[span-p+j]
    
    return CK