B 样条基函数

EasonZzzzzzz

已于 2024-02-18 19:31:32 修改

阅读量2k

点赞数 20

分类专栏：数学之美文章标签：机器学习算法人工智能

于 2024-02-18 18:30:59 首次发布

本文链接：https://blog.csdn.net/m0_72748751/article/details/136152823

版权

数学之美专栏收录该内容

13 篇文章 2 订阅

订阅专栏

一、B 样条基函数的定义和性质

令 $U=\{u_0,u_1,\cdots,u_m\}$ 是一个单调不减的实数序列，即 $u_i\leq u_{i+1},i=0,1,\cdots,m-1$ 。其中， $u_i$ 称为节点， $U$ 称为节点矢量，用 $N_{i,p}(u)$ 表示第 $i$ 个 $p$ 次（ $p + 1$ 阶）B 样条基函数，其定义为（也称为Cox-deBoor 公式）
$\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}$

由此可知：

$N_{i,0}(u)$ 是一个阶梯函数，它在半开区间 $u\in [u_i,u_{i+1})$ 外都为零；
当 $p > 0$ 时， $N_{i,p}(u)$ 是两个 $p - 1$ 次基函数的线性组合；
计算一组基函数时需要事先指定节点矢量 $U$ 和次数 $p$ ；
（1）式中可能出现 $0/0$ ，我们规定 $0/0 = 0$ ；
半开区间 $u\in [u_i,u_{i+1})$ 称为第 $i$ 个节点区间，它的长度可以为零，因为相邻节点可以是相同的。

例1. 令 $U=\{u_0=0,u_1=0,u_2=0,u_3=1,u_4=1,u_5=1\},p=2$ ，我们分别计算 0 次、1 次和 2 次的 B 样条基函数。
$\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}$

二、B 样条基函数的导数

基函数的求导公式为
$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}$

一般求导公式为
$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}$

（4）式是（2）式的另一个推广，它根据基函数 $N_{i,p-k},\cdots,N_{i+k,p-k}$ 来计算 $N_{i,p}(u)$ 的 $k$ 阶导数：
$N^{(k)}_{i,p}=\frac{p!}{(p-k)!}\sum_{j=0}^ka_{k,j}N_{i+j,p-k}\tag{4}$

其中，
$\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}$

这里给出另一个计算 B 样条基函数各阶导数的公式
$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}$

一旦次数已给定，则节点矢量完全决定了基函数 $N_{i,p}(u)$ 。节点矢量的类型有几种，本文中我们只考虑非周期节点矢量，它们具有以下形式：
$U=\{\underbrace{a,\cdots,a}_{p+1},u_{p+1},\cdots,u_{m-p-1},\underbrace{b,\cdots,b}_{p+1}\}\tag{6}$

即第一个和最后一个节点具有复杂度 $p + 1$ 。

三、B 样条基函数的计算

令 $U=\{u_0,u_1,\cdots,u_m\}$ 为形如（6）式的节点矢量，假设我们感兴趣的是 $p$ 次基函数，并假设 $u$ 是固定的， $u\in[u_i,u_{i+1})$ 。我们给出以下五个算法，分别用来计算：

节点区间的下标 $i$ ；
$N_{i-p,p},\cdots,N_{i,p}$ $\big($ 基于（1）式 $\big)$ ；
$N_{i-p,p}^{(k)}(u),\cdots,N_{i,p}^{(k)}(u),k=0,1,\cdots,p$ ；当 $k > p$ 时导数为 0 $\big($ 这个算法基于（4）式 $\big)$
单个基函数 $N_{j,p}(u),j=0,1,\cdots,m-p-1$ ；
单个基函数的导数 $N_{j,p}^{(k)}(u),j=0,1,\cdots,m-p-1,k=0,1,\cdots,p$ $\big($ 基于（3）式 $\big)$

基函数和其导数计算的第一步就是确定 $u$ 所属的节点区间，这可以通过对节点矢量使用线性搜索或二分法搜索得到。这里我们利用二分搜索法。由于我们使用的区间形式为 $u\in[u_i,u_{i+1})$ ，在计算基函数时需要考虑的一个情况是 $u=u_m$ 的情形，这时，我们直接将节点区间的下标设置为 $m - p - 1$ ，因此，在这种情况下 $u\in(u_{m-p-1},u_{m-p}]$ 。下面的 FindSpan 算法返回 $u$ 所在的节点区间的下标 $i$ 。

def FindSpan(p, u, U):
    """
    确定参数 u 所在的节点区间的下标
    :param p: 基函数次数
    :param u: 固定点
    :param U: 节点矢量
    :return: 节点区间的下标
    """
    start, end = 0, len(U) - 1
    if u == U[-1]:
        return len(U) - p - 2
    if u == U[0]:
        return p
    while end != start + 1 and end != start:
        mid_index = (start + end) // 2
        if U[mid_index] == u:
            return mid_index
        if u > U[mid_index]:
            start = mid_index
        else:
            end = mid_index
    return start

现在我们考虑第二个算法，假设 $u$ 在第 $i$ 个节点区间内，计算所有非零 B 样条基函数。
例2. 设 $p=2,U=\{0,0,0,1,2,3,4,4,5,5,5\},u=\dfrac{5}{2}$ 。因为 $u\in [u_4,u_5)$ ，因此 $i = 4$ ，所以我们需要计算
$\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)$

把 $u=\dfrac{5}{2}$ 代入（1）式得到
$\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}$

在实际进行上述计算过程中，很容易发现计算时存在大量冗余的计算，从而设计 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.0
    left = 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] - u
        saved = 0.0
        for r in range(j):
            temp = N[r] / (right[r+1] + left[j-r])
            N[r] = saved + right[r+1] * temp
            saved = left[j-r] * temp
        N[j] = saved
    return N

现在考虑第三个算法，用来计算所有的基函数及其导数 $N_{r,p}^{(k)}(u),r=i-p,\cdots,i,k=0,1,\cdots,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.0
    left = 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] - u
        saved = 0.0
        for 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] * temp
            saved = left[j-r] * temp
        ndu[j, j] = saved

    for j in range(p+1):
        "载入基函数的值"
        ders[0, j] = ndu[j, p]
    "下面计算导数"
    for r in range(p+1):
        s1 = 0
        s2 = 1
        a[0, 0] = 1.0
        "循环计算 k 阶导数，k=1,2,…,n"
        for k in range(1, n+1):
            d = 0.0
            rk = r - k
            pk = p - k
            if r >= k:
                a[s2, 0] = a[s1, 0] / ndu[pk+1, rk]
                d = a[s2, 0] * ndu[rk, pk]
            if rk >= -1:
                j1 = 1
            else:
                j1 = -rk
            if r-1 <= pk:
                j2 = k - 1
            else:
                j2 = p - r
            for 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] = d
            j = s1
            s1 = s2
            s2 = j
    "对结果乘以正确的因子"
    r = p
    for k in range(1, n+1):
        for j in range(p+1):
            ders[k, j] *= r
            r *= (p - k)
    return ders

EasonZzzzzzz

关注

20
点赞
踩
31

收藏

觉得还不错? 一键收藏
1
评论
B 样条基函数

令Uu0u1⋯um是一个单调不减的实数序列，即ui≤ui1i01⋯m−1。其中，ui称为节点，U称为节点矢量，用Nipu表示第i个p次（p1阶）B 样条基函数，其定义为（也称为Ni0uNipu10ui≤uui1elseuip−uiu−uiNip−1uuip1−ui1ui。
复制链接

扫一扫