B样条曲线的公式推导及代码实现

本文仅简述B样条曲线的公式推导，并给出了一种代码实现。在阅读本文之前，请确保你已经对B样条曲线的背景知识有所了解。相关知识可以通过以下课程进行学习：MOOC-计算机图形学-中国农业大学-赵明或者观看B站搬运版。

公式定义

给定如下参数：

• $n+1$个控制点$P_i(i=0,1,\cdots ,n)$
• $n+k+1(2⩽k⩽n+1)$个参数： { u i ∣ i = 0 , 1 , ⋯   , n + k } \{u_i|i=0,1,\cdots,n+k\} {ui​∣i=0,1,⋯,n+k}

构造多项式函数
$$\mathcal{P}(u)=\sum _{i=0}^n{P}_i{B}_{i,k}(u),u\in [u_{k-1},u_{n+1}]\text{\hspace{1em}(1)}$$来逼近曲线。其中$B_{i,k}(u)$称为$k$阶（$k-1$次）B样条基函数，它是由DeBoor-Cox递推公式来定义的：
$$\begin{array}{rl}& B_{i,1}(u)=\{\begin{array}{ll}1,& u_i⩽u<u_{i+1}\\ 0,& otherwise\end{array}，\\ & B_{i,k}(u)=\frac{u-u_i}{u_{i+k-1}-u_i}{B}_{i,k-1}(u)+\frac{u_{i+k}-u}{u_{i+k}-u_{i+1}}{B}_{i+1,k-1}(u)。\end{array}\text{\hspace{1em}(2)}$$递推式中若遇到$0$除以$0$，则除法结果取值为$0$。

$B_{i,k}(u)$可以看作是控制点$P_i$在$u$下的权值，因为B样条基函数具有权性：
$$\sum _{i=0}^n{B}_{i,k}(u)\equiv 1,\forall u。\text{\hspace{1em}(3)}$$

对于某个$u\in [u_j,u_{j+1})$，运用公式(2)递推$B_{i,k}(u)$，使得最终的结果中只含$B_{x,1}(u)$，我们可以推出$x=i,i+1,\cdots ,i+k-1$。又根据$B_{x,1}(u)$的定义式(2)，我们发现仅当 j ∈ { i , i + 1 , ⋯ i + k − 1 } j\in\{i,i+1,\cdots i+k-1\} j∈{i,i+1,⋯i+k−1}，即 i ∈ { j − ( k − 1 ) , ⋯   , j } i\in\{j-(k-1),\cdots,j\} i∈{j−(k−1),⋯,j}时，$B_{i,k}(u)$有可能不为0。

举个例子来解释一下。以$i=7,k=3$为例，运用公式(2)递推$B_{7,3}(u)$：
$$\begin{array}{rl}{B}_{7,3}(u)& =c_0{B}_{7,2}(u)+c_1{B}_{8,2}(u)\\ & =c_2{B}_{7,1}(u)+c_3{B}_{8,1}(u)+c_4{B}_{9,1}(u)\end{array}$$（其中$c_0,c_1,\cdots ,c_4$均为系数），我们发现最终的结果中只含$B_{7,1}(u),B_{8,1}(u),B_{9,1}(u)$。对于某个$u\in [u_j,u_{j+1})$，仅当 j ∈ { 7 , 8 , 9 } j\in\{7,8,9\} j∈{7,8,9}时$B_{7,1}(u),B_{8,1}(u),B_{9,1}(u)$不全为0，从而使得$B_{7,3}(u)$有可能不为0。

算法原理

对于某个$u\in [u_j,u_{j+1}]$，我们定义
$$d_{i,0}=P_i,i=j-(k-1),\cdots ,j，\text{\hspace{1em}(4)}$$ $$\begin{array}{rl}{\alpha }_{i,r}& =\frac{u-u_i}{u_{i+k-r}-u_i},\\ d_{i,r}& =(1-{\alpha }_{i,r})d_{i-1,r-1}+{\alpha }_{i,r}{d}_{i,r-1}\end{array}\left(\begin{array}{l}r=1,2,\cdots ,k-1\\ i=j-(k-1)+r,\cdots ,j\end{array}\right)。\text{\hspace{1em}(5)}$$

这里对公式(5)中的$i=j-(k-1)+r,\cdots ,j$做出解释：递推$d_{i,r}$使最终结果中只含$d_{x,0}$，得到$x=i-r,\cdots ,i$。由$d_{x,0}$的定义式(4)知， ∀ x , x ∈ { j − ( k − 1 ) , ⋯   , j } \forall x, x\in\{j-(k-1),\cdots,j\} ∀x,x∈{j−(k−1),⋯,j}。因此有 i − r ∈ { j − ( k − 1 ) , ⋯   , j } , ⋯   , i ∈ { j − ( k − 1 ) , ⋯   , j } i-r\in\{j-(k-1),\cdots,j\},\cdots,i\in\{j-(k-1),\cdots,j\} i−r∈{j−(k−1),⋯,j},⋯,i∈{j−(k−1),⋯,j}，从而得到 i ∈ { j − ( k − 1 ) + r , ⋯   , j } i\in\{j-(k-1)+r,\cdots,j\} i∈{j−(k−1)+r,⋯,j}。

可以证明
$$\mathcal{P}(u)=d_{j,k-1}。\text{\hspace{1em}(6)}$$

于是得到算法逻辑如下（假定输入只提供了$n+1$个控制点，参数 u i ( i ∈ { 0 , 1 , ⋯   , n + k } ) u_i(i\in\{0,1,\cdots,n+k\}) ui​(i∈{0,1,⋯,n+k})由我们自定义选取）：

• 在实数域上任取一段区间，不妨取$[0,1]$。
• 将$[0,1]$均匀分隔为$n+k$段（控制点个数=$n+1$，曲线阶数=$k$）。
• 分隔之后产生了$n+k+1$个分隔点，依次记为$u_0,u_1,\cdots ,u_{n+k}$。于是有$u_j=j/(n+k)$。
• 在区间$[u_{k-1},u_{n+1}]$上大量采样。对于每个样本$u$，计算$\mathcal{P}(u)$。
• 根据公式(6)，我们可以通过计算$d_{j,k-1}$来计算$\mathcal{P}(u)$。
• $d_{j,k-1}$可以通过公式(5)进行递归计算。

算法实现

def draw_curve(p_list, k: int = 4):
	"""
    :param p_list: (list of list of int: [[x0, y0], [x1, y1], [x2, y2], ...]) 曲线的控制点坐标列表
    :param k: B样条曲线的阶数，默认为4阶，即3次均匀B样条曲线
    :return: (list of list of int: [[x_0, y_0], [x_1, y_1], [x_2, y_2], ...]) 绘制结果的像素点坐标列表
    """
    n = len(p_list) - 1
    if n == 0:
        return p_list
    elif n == 1:
        return draw_line(p_list) # 调用绘制直线的算法
        
    assert 2 <= k <= n + 1
    freq = 500
    result = []
    
    for j in range(k - 1, n + 1):
        seg = n + k
        u_j = j / seg
        for u in range(freq + 1):
            u = u / (freq * seg) + u_j
            tmp = [
                [p_list[i + j - (k-1)][0],
                 p_list[i + j - (k-1)][1]
                ] for i in range(k)
            ]
            for r in range(1, k):
                for i in range(k - 1, r - 1, -1):
                    ii = i + j - (k - 1)
                    u_i = ii / seg
                    alpha = (u - u_i) / ((ii+k-r)/seg - u_i)
                    tmp[i] = [
                        (1-alpha)*tmp[i-1][0] + alpha*tmp[i][0],
                        (1-alpha)*tmp[i-1][1] + alpha*tmp[i][1]
                    ]
            result.append(
                [round(tmp[k-1][0]), round(tmp[k-1][1])]
            )
    return result

补充：$\mathcal{P}(u)$定义域的推导

在公式(1)中，我们直接定义变量$u$的取值范围为$[u_{k-1},u_{n+1}]$。这一节中我们将详细说明其原理。

根据递推公式(2)，我们可以推知：

• $B_{i,1}$涉及1个区间（$[u_i,u_{i+1}]$）2个控制点（$u_i$和$u_{i+1}$）
• $B_{i,2}$由$B_{i,1}$和$B_{i+1,1}$组成，涉及2个区间3个控制点
• $B_{i,3}$由$B_{i,2}$和$B_{i+1,2}$组成，涉及3个区间4个控制点
• $\cdots$
• $B_{i,k}$涉及$k$个区间$k+1$个控制点

从而可以定义$B_{i,k}$的支承区间为$[u_i,u_{i+k}]$。

举个例子，如上图，其中$n=4,k=4$。根据支承区间的定义，我们有：

• $P_0{B}_{0,4}$涉及节点$u_0$到$u_4$
• $P_1{B}_{1,4}$涉及节点$u_1$到$u_5$
• $\cdots$
• $P_4{B}_{4,4}$涉及节点$u_4$到$u_8$

对于每个区间$[u_i,u_{i+1}]$，它至多包含于$k$个不同的支承区间中，也就是说，它至多与$k$个控制点相关。下图形象地展示了这种关系。

我们的目标是，对于$\mathcal{P}(u)$的定义域中的每个区间，要使它与尽可能多的控制点相关。于是我们最终选择了$[u_{k-1},u_{n+1}]$。可以证明，$\forall [u_i,u_{i+1}]\subseteq [u_{k-1},u_{n+1}]$，$[u_i,u_{i+1}]$都与$k$个控制点相关。

更多细节

由于篇幅原因，文中有些证明并未详细展开叙述。更多的证明细节请参阅相关论文。例如C. DE BOOR, On calculating with B-splines, 1972等。