一、B 样条基函数的定义和性质
令 U = { u 0 , u 1 , ⋯ , u m } U=\{u_0,u_1,\cdots,u_m\} U={u0,u1,⋯,um} 是一个单调不减的实数序列，即 u i ≤ u i + 1 , i = 0 , 1 , ⋯ , m − 1 u_i\leq u_{i+1},i=0,1,\cdots,m-1 ui≤ui+1,i=0,1,⋯,m−1。其中， u i u_i ui 称为节点， U U U 称为节点矢量，用 N i , p ( u ) N_{i,p}(u) Ni,p(u) 表示第 i i i 个 p p p 次（ p + 1 p+1 p+1 阶）B 样条基函数，其定义为（也称为Cox-deBoor 公式）
N i , 0 ( u ) = { 1 , u i ≤ u < u i + 1 , 0 , e l s e N i , p ( u ) = u − u i u i + p − u i N i , p − 1 ( u ) + u i + p + 1 − u u i + p + 1 − u i + 1 N i + 1 , p − 1 ( u ) (1) \begin{aligned} N_{i,0}(u)&=\begin{cases}1,\quad &u_i\leq u< u_{i+1},\\ 0,&else \end{cases}\\[3ex] N_{i,p}(u)&=\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)\tag{1} \end{aligned} Ni,0(u)Ni,p(u)={1,0,ui≤u<ui+1,else=ui+p−uiu−uiNi,p−1(u)+ui+p+1−ui+1ui+p+1−uNi+1,p−1(u)(1)
由此可知：
- N i , 0 ( u ) N_{i,0}(u) Ni,0(u) 是一个阶梯函数，它在半开区间 u ∈ [ u i , u i + 1 ) u\in [u_i,u_{i+1}) u∈[ui,ui+1) 外都为零；
- 当 p > 0 p>0 p>0 时， N i , p ( u ) N_{i,p}(u) Ni,p(u) 是两个 p − 1 p-1 p−1 次基函数的线性组合；
- 计算一组基函数时需要事先指定节点矢量 U U U 和次数 p p p；
- （1）式中可能出现 0 / 0 0/0 0/0，我们规定 0 / 0 = 0 0/0=0 0/0=0；
- 半开区间 u ∈ [ u i , u i + 1 ) u\in [u_i,u_{i+1}) u∈[ui,ui+1) 称为第 i i i 个节点区间，它的长度可以为零，因为相邻节点可以是相同的。
例1. 令 U = { u 0 = 0 , u 1 = 0 , u 2 = 0 , u 3 = 1 , u 4 = 1 , u 5 = 1 } , p = 2 U=\{u_0=0,u_1=0,u_2=0,u_3=1,u_4=1,u_5=1\},p=2 U={u0=0,u1=0,u2=0,u3=1,u4=1,u5=1},p=2，我们分别计算 0 次、1 次和 2 次的 B 样条基函数。
N 0 , 0 = N 1 , 0 = 0 , − ∞ < u < ∞ N 2 , 0 = { 1 , 0 ≤ u < 1 0 , 其他 N 3 , 0 = N 4 , 0 = 0 , − ∞ < u < ∞ N 0 , 1 = u − 0 0 − 0 N 0 , 0 + 0 − u 0 − 0 N 1 , 0 = 0 , − ∞ < u < ∞ N 1 , 1 = u − 0 0 − 0 N 1 , 0 + 1 − u 1 − 0 N 2 , 0 = { 1 − u , 0 ≤ u < 1 0 , 其他 N 2 , 1 = u − 0 1 − 0 N 2 , 0 + 1 − u 1 − 1 N 3 , 0 = { u , 0 ≤ u < 1 0 , 其他 N 3 , 1 = u − 1 1 − 1 N 3 , 0 + 1 − u 1 − 1 N 4 , 0 = 0 , − ∞ < u < ∞ N 0 , 2 = u − 0 0 − 0 N 0 , 1 + 1 − u 1 − 0 N 1 , 1 = { ( 1 − u ) 2 , 0 ≤ u < 1 0 , 其他 N 1 , 2 = u − 0 1 − 0 N 1 , 1 + 1 − u 1 − 0 N 2 , 1 = { 2 u ( 1 − u ) , 0 ≤ u < 1 0 , 其他 N 2 , 2 = u − 0 1 − 0 N 2 , 1 + 1 − u 1 − 1 N 3 , 1 = { u 2 , 0 ≤ u < 1 0 , 其他 \begin{aligned} &N_{0,0}=N_{1,0}=0,\quad -\infty<u<\infty\\[2ex] &N_{2,0}=\begin{cases}1,\quad &0\leq u<1\\[2ex] 0,&其他 \end{cases}\\[2ex] &N_{3,0}=N_{4,0}=0,\quad -\infty<u<\infty\\[2ex] &N_{0,1}=\frac{u-0}{0-0}N_{0,0}+\frac{0-u}{0-0}N_{1,0}=0,\quad -\infty<u<\infty\\[2ex] &N_{1,1}=\frac{u-0}{0-0}N_{1,0}+\frac{1-u}{1-0}N_{2,0}=\begin{cases}1-u,\quad &0\leq u<1\\[2ex]0,&其他\end{cases}\\[2ex] &N_{2,1}=\frac{u-0}{1-0}N_{2,0}+\frac{1-u}{1-1}N_{3,0}=\begin{cases}u,\quad &0\leq u<1\\[2ex]0,&其他\end{cases}\\[2ex] &N_{3,1}=\frac{u-1}{1-1}N_{3,0}+\frac{1-u}{1-1}N_{4,0}=0,\quad -\infty<u<\infty\\[2ex] &N_{0,2}=\frac{u-0}{0-0}N_{0,1}+\frac{1-u}{1-0}N_{1,1}=\begin{cases}(1-u)^2,\quad &0\leq u<1\\[2ex]0,&其他\end{cases}\\[2ex] &N_{1,2}=\frac{u-0}{1-0}N_{1,1}+\frac{1-u}{1-0}N_{2,1}=\begin{cases}2u(1-u),\quad &0\leq u<1\\[2ex]0,&其他\end{cases}\\[2ex] &N_{2,2}=\frac{u-0}{1-0}N_{2,1}+\frac{1-u}{1-1}N_{3,1}=\begin{cases}u^2,\quad &0\leq u<1\\[2ex]0,&其他\end{cases}\\[2ex] \end{aligned} N0,0=N1,0=0,−∞<u<∞N2,0=⎩ ⎨ ⎧1,0,0≤u<1其他N3,0=N4,0=0,−∞<u<∞N0,1=0−0u−0N0,0+0−00−uN1,0=0,−∞<u<∞N1,1=0−0u−0N1,0+1−01−uN2,0=⎩ ⎨ ⎧1−u,0,0≤u<1其他N2,1=1−0u−0N2,0+1−11−uN3,0=⎩ ⎨ ⎧u,0,0≤u<1其他N3,1=1−1u−1N3,0+1−11−uN4,0=0,−∞<u<∞N0,2=0−0u−0N0,1+1−01−uN1,1=⎩ ⎨ ⎧(1−u)2,0,0≤u<1其他N1,2=1−0u−0N1,1+1−01−uN2,1=⎩ ⎨ ⎧2u(1−u),0,0≤u<1其他N2,2=1−0u−0N2,1+1−11−uN3,1=⎩ ⎨ ⎧u2,0,0≤u<1其他
二、B 样条基函数的导数
基函数的求导公式为
N i , p ′ = p u i + p − u i N i , p − 1 ( u ) + p u i + p + 1 − u i + 1 N i + 1 , p − 1 ( u ) (2) N^\prime_{i,p}=\frac{p}{u_{i+p}-u_i}N_{i,p-1}(u)+\frac{p}{u_{i+p+1}-u_{i+1}}N_{i+1,p-1}(u)\tag{2} Ni,p′=ui+p−uipNi,p−1(u)+ui+p+1−ui+1pNi+1,p−1(u)(2)
一般求导公式为
N i , p ( k ) = p ( N i , p − 1 ( k − 1 ) u i + p − u i − N i + 1 , p − 1 ( k − 1 ) u i + p + 1 − u i + 1 ) (3) N^{(k)}_{i,p}=p\Big(\frac{N_{i,p-1}^{(k-1)}}{u_{i+p}-u_i}-\frac{N_{i+1,p-1}^{(k-1)}}{u_{i+p+1}-u_{i+1}}\Big)\tag{3} Ni,p(k)=p(ui+p−uiNi,p−1(k−1)−ui+p+1−ui+1Ni+1,p−1(k−1))(3)
（4）式是（2）式的另一个推广，它根据基函数 N i , p − k , ⋯ , N i + k , p − k N_{i,p-k},\cdots,N_{i+k,p-k} Ni,p−k,⋯,Ni+k,p−k 来计算 N i , p ( u ) N_{i,p}(u) Ni,p(u) 的 k k k 阶导数：
N i , p ( k ) = p ! ( p − k ) ! ∑ j = 0 k a k , j N i + j , p − k (4) N^{(k)}_{i,p}=\frac{p!}{(p-k)!}\sum_{j=0}^ka_{k,j}N_{i+j,p-k}\tag{4} Ni,p(k)=(p−k)!p!j=0∑kak,jNi+j,p−k(4)
其中，
a 0 , 0 = 1 a k , 0 = a k − 1 , 0 u i + p − k + 1 − u i a k , j = a k − 1 , j − a k − 1 , j − 1 u i + p + j − k + 1 − u i + j , j = 1 , 2 , ⋯ , k − 1 a k , k = − a k − 1 , k − 1 u i + p + 1 − u i + k \begin{aligned} &a_{0,0}=1\\[2ex] &a_{k,0}=\frac{a_{k-1,0}}{u_{i+p-k+1}-u_i}\\[2ex] &a_{k,j}=\frac{a_{k-1,j}-a_{k-1,j-1}}{u_{i+p+j-k+1}-u_{i+j}},\quad j=1,2,\cdots,k-1\\[2ex] &a_{k,k}=\frac{-a_{k-1,k-1}}{u_{i+p+1}-u_{i+k}} \end{aligned} a0,0=1ak,0=ui+p−k+1−uiak−1,0ak,j=ui+p+j−k+1−ui+jak−1,j−ak−1,j−1,j=1,2,⋯,k−1ak,k=ui+p+1−ui+k−ak−1,k−1
这里给出另一个计算 B 样条基函数各阶导数的公式
N i , p ( k ) = p p − k ( u − u i u i + p − u i N i , p − 1 ( k ) − u i + p + 1 − u u i + p + 1 − u i + 1 N i + 1 , p − 1 ( k ) ) , k = 0 , 1 , ⋯ , p − 1 (5) N_{i,p}^{(k)}=\frac{p}{p-k}\Big(\frac{u-u_i}{u_{i+p}-u_i}N_{i,p-1}^{(k)}-\frac{u_{i+p+1}-u}{u_{i+p+1}-u_{i+1}}N_{i+1,p-1}^{(k)}\Big),\quad k=0,1,\cdots,p-1\tag{5} Ni,p(k)=p−kp(ui+p−uiu−uiNi,p−1(k)−ui+p+1−ui+1ui+p+1−uNi+1,p−1(k)),k=0,1,⋯,p−1(5)
一旦次数已给定，则节点矢量完全决定了基函数 N i , p ( u ) N_{i,p}(u) Ni,p(u)。节点矢量的类型有几种，本文中我们只考虑非周期节点矢量，它们具有以下形式：
U = { a , ⋯ , a ⏟ p + 1 , u p + 1 , ⋯ , u m − p − 1 , b , ⋯ , b ⏟ p + 1 } (6) U=\{\underbrace{a,\cdots,a}_{p+1},u_{p+1},\cdots,u_{m-p-1},\underbrace{b,\cdots,b}_{p+1}\}\tag{6} U={p+1 a,⋯,a,up+1,⋯,um−p−1,p+1 b,⋯,b}(6)
即第一个和最后一个节点具有复杂度 p + 1 p+1 p+1。
三、B 样条基函数的计算
令 U = { u 0 , u 1 , ⋯ , u m } U=\{u_0,u_1,\cdots,u_m\} U={u0,u1,⋯,um} 为形如（6）式的节点矢量，假设我们感兴趣的是 p p p 次基函数，并假设 u u u 是固定的， u ∈ [ u i , u i + 1 ) u\in[u_i,u_{i+1}) u∈[ui,ui+1)。我们给出以下五个算法，分别用来计算：
- 节点区间的下标 i i i；
- N i − p , p , ⋯ , N i , p N_{i-p,p},\cdots,N_{i,p} Ni−p,p,⋯,Ni,p ( \big( ( 基于（1）式 ) \big) )；
- N i − p , p ( k ) ( u ) , ⋯ , N i , p ( k ) ( u ) , k = 0 , 1 , ⋯ , p N_{i-p,p}^{(k)}(u),\cdots,N_{i,p}^{(k)}(u),k=0,1,\cdots,p Ni−p,p(k)(u),⋯,Ni,p(k)(u),k=0,1,⋯,p；当 k > p k>p k>p 时导数为 0 ( \big( ( 这个算法基于（4）式 ) \big) )
- 单个基函数 N j , p ( u ) , j = 0 , 1 , ⋯ , m − p − 1 N_{j,p}(u),j=0,1,\cdots,m-p-1 Nj,p(u),j=0,1,⋯,m−p−1；
- 单个基函数的导数 N j , p ( k ) ( u ) , j = 0 , 1 , ⋯ , m − p − 1 , k = 0 , 1 , ⋯ , p N_{j,p}^{(k)}(u),j=0,1,\cdots,m-p-1,k=0,1,\cdots,p Nj,p(k)(u),j=0,1,⋯,m−p−1,k=0,1,⋯,p ( \big( ( 基于（3）式 ) \big) )
基函数和其导数计算的第一步就是确定 u u u 所属的节点区间，这可以通过对节点矢量使用线性搜索或二分法搜索得到。这里我们利用二分搜索法。由于我们使用的区间形式为 u ∈ [ u i , u i + 1 ) u\in[u_i,u_{i+1}) u∈[ui,ui+1)，在计算基函数时需要考虑的一个情况是 u = u m u=u_m u=um 的情形，这时，我们直接将节点区间的下标设置为 m − p − 1 m-p-1 m−p−1，因此，在这种情况下 u ∈ ( u m − p − 1 , u m − p ] u\in(u_{m-p-1},u_{m-p}] u∈(um−p−1,um−p]。下面的 FindSpan 算法返回 u u u 所在的节点区间的下标 i i i。
def FindSpan(p, u, U):"""确定参数 u 所在的节点区间的下标:param p: 基函数次数:param u: 固定点:param U: 节点矢量:return: 节点区间的下标"""start, end = 0, len(U) - 1if u == U[-1]:return len(U) - p - 2if u == U[0]:return pwhile end != start + 1 and end != start:mid_index = (start + end) // 2if U[mid_index] == u:return mid_indexif u > U[mid_index]:start = mid_indexelse:end = mid_indexreturn start
现在我们考虑第二个算法，假设 u u u 在第 i i i 个节点区间内，计算所有非零 B 样条基函数。
例2. 设 p = 2 , U = { 0 , 0 , 0 , 1 , 2 , 3 , 4 , 4 , 5 , 5 , 5 } , u = 5 2 p=2,U=\{0,0,0,1,2,3,4,4,5,5,5\},u=\dfrac{5}{2} p=2,U={0,0,0,1,2,3,4,4,5,5,5},u=25。因为 u ∈ [ u 4 , u 5 ) u\in [u_4,u_5) u∈[u4,u5)，因此 i = 4 i=4 i=4，所以我们需要计算
N 2 , 2 ( 5 2 ) N 3 , 1 ( 5 2 ) N 4 , 0 ( 5 2 ) N 3 , 2 ( 5 2 ) N 4 , 1 ( 5 2 ) N 4 , 2 ( 5 2 ) \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad N_{2,2}\Big(\dfrac{5}{2}\Big)\\[2ex] N_{3,1}\Big(\dfrac{5}{2}\Big)\\[2ex] N_{4,0}\Big(\dfrac{5}{2}\Big)\quad\quad\quad\quad\quad\quad\quad\quad\quad N_{3,2}\Big(\dfrac{5}{2}\Big)\\[2ex] N_{4,1}\Big(\dfrac{5}{2}\Big)\\[2ex] \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad N_{4,2}\Big(\dfrac{5}{2}\Big) N2,2(25)N3,1(25)N4,0(25)N3,2(25)N4,1(25)N4,2(25)
把 u = 5 2 u=\dfrac{5}{2} u=25 代入（1）式得到
N 4 , 0 ( 5 2 ) = 1 N 3 , 1 ( 5 2 ) = 1 2 N 4 , 1 ( 5 2 ) = 1 2 N 2 , 2 ( 5 2 ) = 1 8 N 3 , 2 ( 5 2 ) = 6 8 N 4 , 2 ( 5 2 ) = 1 8 \begin{aligned} &N_{4,0}\Big(\dfrac{5}{2}\Big)=1\\[3ex] &N_{3,1}\Big(\dfrac{5}{2}\Big)=\dfrac{1}{2}\quad N_{4,1}\Big(\dfrac{5}{2}\Big)=\dfrac{1}{2}\\[3ex] &N_{2,2}\Big(\dfrac{5}{2}\Big)=\dfrac{1}{8}\quad N_{3,2}\Big(\dfrac{5}{2}\Big)=\dfrac{6}{8}\quad N_{4,2}\Big(\dfrac{5}{2}\Big)=\dfrac{1}{8} \end{aligned} N4,0(25)=1N3,1(25)=21N4,1(25)=21N2,2(25)=81N3,2(25)=86N4,2(25)=81
在实际进行上述计算过程中，很容易发现计算时存在大量冗余的计算，从而设计 BasisFuns() 算法来计算所有非零 B 样条基函数并把它们存储在数组 N 中。
def BasisFuns(i, u, p, U):"""计算所有非零 B 样条基函数的值:param i: 所在的节点区间:param u: 固定值:param p: 基函数次数:param U: 节点矢量:return: 基函数值数组"""N = np.zeros(p+1)N[0] = 1.0left = np.zeros(p+1)right = np.zeros(p+1)for j in range(1, p+1):left[j] = u - U[i+1-j]right[j] = U[i+j] - usaved = 0.0for r in range(j):temp = N[r] / (right[r+1] + left[j-r])N[r] = saved + right[r+1] * tempsaved = left[j-r] * tempN[j] = savedreturn N
现在考虑第三个算法，用来计算所有的基函数及其导数 N r , p ( k ) ( u ) , r = i − p , ⋯ , i , k = 0 , 1 , ⋯ , n N_{r,p}^{(k)}(u),r=i-p,\cdots,i,k=0,1,\cdots,n Nr,p(k)(u),r=i−p,⋯,i,k=0,1,⋯,n
def DersBasisFuns(i, u, p, n, U):"""计算非零 B 样条基函数及其导数:param i: 所在的节点区间:param u: 固定值:param p: 基函数次数:param n: n 阶导数:param U: 节点矢量:return: ders""""第一部分是对算法 BasisFuns() 进行修改，将函数值和节点差存到数组中"ndu = np.zeros((p+1, p+1))ders = np.zeros((n+1, p+1))a = np.zeros((2, p+1))ndu[0, 0] = 1.0left = np.zeros(p + 1)right = np.zeros(p + 1)for j in range(1, p+1):left[j] = u - U[i+1-j]right[j] = U[i+j] - usaved = 0.0for r in range(j):"下三角"ndu[j, r] = right[r+1] + left[j-r]temp = ndu[r, j-1] / ndu[j, r]"上三角"ndu[r, j] = saved + right[r+1] * tempsaved = left[j-r] * tempndu[j, j] = savedfor j in range(p+1):"载入基函数的值"ders[0, j] = ndu[j, p]"下面计算导数"for r in range(p+1):s1 = 0s2 = 1a[0, 0] = 1.0"循环计算 k 阶导数，k=1,2,…,n"for k in range(1, n+1):d = 0.0rk = r - kpk = p - kif r >= k:a[s2, 0] = a[s1, 0] / ndu[pk+1, rk]d = a[s2, 0] * ndu[rk, pk]if rk >= -1:j1 = 1else:j1 = -rkif r-1 <= pk:j2 = k - 1else:j2 = p - rfor j in range(j1, j2 + 1):a[s2, j] = (a[s1, j] - a[s1, j-1]) / ndu[pk+1, rk+j]d += a[s2, j] * ndu[rk+j, pk]if r <= pk:a[s2, k] = - a[s1, k-1] / ndu[pk+1, r]d += a[s2, k] * ndu[r, pk]ders[k, r] = dj = s1s1 = s2s2 = j"对结果乘以正确的因子"r = pfor k in range(1, n+1):for j in range(p+1):ders[k, j] *= rr *= (p - k)return ders
)
)
