做网站需要用到的符号语言,自己做视频网站的流程,免费设立网站,黑马程序员教程一、B 样条基函数的定义和性质
令 U { u 0 , u 1 , ⋯ , u m } U\{u_0,u_1,\cdots,u_m\} U{u0,u1,⋯,um} 是一个单调不减的实数序列#xff0c;即 u i ≤ u i 1 , i 0 , 1 , ⋯ , m − 1 u_i\leq u_{i1},i0,1,\cdots,m-1 ui≤ui1,i0,1,⋯,m−1。其中#xff…一、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_{i1},i0,1,\cdots,m-1 ui≤ui1,i0,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 p1 p1 阶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_{i1},\\ 0,else \end{cases}\\[3ex] N_{i,p}(u)\frac{u-u_i}{u_{ip}-u_i}N_{i,p-1}(u)\frac{u_{ip1}-u}{u_{ip1}-u_{i1}}N_{i1,p-1}(u)\tag{1} \end{aligned} Ni,0(u)Ni,p(u){1,0,ui≤uui1,elseuip−uiu−uiNi,p−1(u)uip1−ui1uip1−uNi1,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_{i1}) u∈[ui,ui1) 外都为零当 p 0 p0 p0 时 N i , p ( u ) N_{i,p}(u) Ni,p(u) 是两个 p − 1 p-1 p−1 次基函数的线性组合计算一组基函数时需要事先指定节点矢量 U U U 和次数 p p p1式中可能出现 0 / 0 0/0 0/0我们规定 0 / 0 0 0/00 0/00半开区间 u ∈ [ u i , u i 1 ) u\in [u_i,u_{i1}) u∈[ui,ui1) 称为第 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_00,u_10,u_20,u_31,u_41,u_51\},p2 U{u00,u10,u20,u31,u41,u51},p2我们分别计算 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 -\inftyu\infty\\[2ex] N_{2,0}\begin{cases}1,\quad 0\leq u1\\[2ex] 0,其他 \end{cases}\\[2ex] N_{3,0}N_{4,0}0,\quad -\inftyu\infty\\[2ex] N_{0,1}\frac{u-0}{0-0}N_{0,0}\frac{0-u}{0-0}N_{1,0}0,\quad -\inftyu\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 u1\\[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 u1\\[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 -\inftyu\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 u1\\[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 u1\\[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 u1\\[2ex]0,其他\end{cases}\\[2ex] \end{aligned} N0,0N1,00,−∞u∞N2,0⎩ ⎨ ⎧1,0,0≤u1其他N3,0N4,00,−∞u∞N0,10−0u−0N0,00−00−uN1,00,−∞u∞N1,10−0u−0N1,01−01−uN2,0⎩ ⎨ ⎧1−u,0,0≤u1其他N2,11−0u−0N2,01−11−uN3,0⎩ ⎨ ⎧u,0,0≤u1其他N3,11−1u−1N3,01−11−uN4,00,−∞u∞N0,20−0u−0N0,11−01−uN1,1⎩ ⎨ ⎧(1−u)2,0,0≤u1其他N1,21−0u−0N1,11−01−uN2,1⎩ ⎨ ⎧2u(1−u),0,0≤u1其他N2,21−0u−0N2,11−11−uN3,1⎩ ⎨ ⎧u2,0,0≤u1其他
二、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_{ip}-u_i}N_{i,p-1}(u)\frac{p}{u_{ip1}-u_{i1}}N_{i1,p-1}(u)\tag{2} Ni,p′uip−uipNi,p−1(u)uip1−ui1pNi1,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_{ip}-u_i}-\frac{N_{i1,p-1}^{(k-1)}}{u_{ip1}-u_{i1}}\Big)\tag{3} Ni,p(k)p(uip−uiNi,p−1(k−1)−uip1−ui1Ni1,p−1(k−1))(3)
4式是2式的另一个推广它根据基函数 N i , p − k , ⋯ , N i k , p − k N_{i,p-k},\cdots,N_{ik,p-k} Ni,p−k,⋯,Nik,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_{j0}^ka_{k,j}N_{ij,p-k}\tag{4} Ni,p(k)(p−k)!p!j0∑kak,jNij,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_{ip-k1}-u_i}\\[2ex] a_{k,j}\frac{a_{k-1,j}-a_{k-1,j-1}}{u_{ipj-k1}-u_{ij}},\quad j1,2,\cdots,k-1\\[2ex] a_{k,k}\frac{-a_{k-1,k-1}}{u_{ip1}-u_{ik}} \end{aligned} a0,01ak,0uip−k1−uiak−1,0ak,juipj−k1−uijak−1,j−ak−1,j−1,j1,2,⋯,k−1ak,kuip1−uik−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_{ip}-u_i}N_{i,p-1}^{(k)}-\frac{u_{ip1}-u}{u_{ip1}-u_{i1}}N_{i1,p-1}^{(k)}\Big),\quad k0,1,\cdots,p-1\tag{5} Ni,p(k)p−kp(uip−uiu−uiNi,p−1(k)−uip1−ui1uip1−uNi1,p−1(k)),k0,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}_{p1},u_{p1},\cdots,u_{m-p-1},\underbrace{b,\cdots,b}_{p1}\}\tag{6} U{p1 a,⋯,a,up1,⋯,um−p−1,p1 b,⋯,b}(6)
即第一个和最后一个节点具有复杂度 p 1 p1 p1。
三、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_{i1}) u∈[ui,ui1)。我们给出以下五个算法分别用来计算
节点区间的下标 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),k0,1,\cdots,p Ni−p,p(k)(u),⋯,Ni,p(k)(u),k0,1,⋯,p当 k p kp kp 时导数为 0 ( \big( ( 这个算法基于4式 ) \big) )单个基函数 N j , p ( u ) , j 0 , 1 , ⋯ , m − p − 1 N_{j,p}(u),j0,1,\cdots,m-p-1 Nj,p(u),j0,1,⋯,m−p−1单个基函数的导数 N j , p ( k ) ( u ) , j 0 , 1 , ⋯ , m − p − 1 , k 0 , 1 , ⋯ , p N_{j,p}^{(k)}(u),j0,1,\cdots,m-p-1,k0,1,\cdots,p Nj,p(k)(u),j0,1,⋯,m−p−1,k0,1,⋯,p ( \big( ( 基于3式 ) \big) )
基函数和其导数计算的第一步就是确定 u u u 所属的节点区间这可以通过对节点矢量使用线性搜索或二分法搜索得到。这里我们利用二分搜索法。由于我们使用的区间形式为 u ∈ [ u i , u i 1 ) u\in[u_i,u_{i1}) u∈[ui,ui1)在计算基函数时需要考虑的一个情况是 u u m uu_m uum 的情形这时我们直接将节点区间的下标设置为 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 p2,U\{0,0,0,1,2,3,4,4,5,5,5\},u\dfrac{5}{2} p2,U{0,0,0,1,2,3,4,4,5,5,5},u25。因为 u ∈ [ u 4 , u 5 ) u\in [u_4,u_5) u∈[u4,u5)因此 i 4 i4 i4所以我们需要计算 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} u25 代入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(p1)N[0] 1.0left np.zeros(p1)right np.zeros(p1)for j in range(1, p1):left[j] u - U[i1-j]right[j] U[ij] - usaved 0.0for r in range(j):temp N[r] / (right[r1] left[j-r])N[r] saved right[r1] * 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),ri-p,\cdots,i,k0,1,\cdots,n Nr,p(k)(u),ri−p,⋯,i,k0,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((p1, p1))ders np.zeros((n1, p1))a np.zeros((2, p1))ndu[0, 0] 1.0left np.zeros(p 1)right np.zeros(p 1)for j in range(1, p1):left[j] u - U[i1-j]right[j] U[ij] - usaved 0.0for r in range(j):下三角ndu[j, r] right[r1] left[j-r]temp ndu[r, j-1] / ndu[j, r]上三角ndu[r, j] saved right[r1] * tempsaved left[j-r] * tempndu[j, j] savedfor j in range(p1):载入基函数的值ders[0, j] ndu[j, p]下面计算导数for r in range(p1):s1 0s2 1a[0, 0] 1.0循环计算 k 阶导数k1,2,…,nfor k in range(1, n1):d 0.0rk r - kpk p - kif r k:a[s2, 0] a[s1, 0] / ndu[pk1, 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[pk1, rkj]d a[s2, j] * ndu[rkj, pk]if r pk:a[s2, k] - a[s1, k-1] / ndu[pk1, r]d a[s2, k] * ndu[r, pk]ders[k, r] dj s1s1 s2s2 j对结果乘以正确的因子r pfor k in range(1, n1):for j in range(p1):ders[k, j] * rr * (p - k)return ders
