本文仅简述B样条曲线的公式推导,并给出了一种代码实现。在阅读本文之前,请确保你已经对B样条曲线的背景知识有所了解。相关知识可以通过以下课程进行学习:MOOC-计算机图形学-中国农业大学-赵明或者观看B站搬运版。
公式定义
给定如下参数:
-
n
+
1
n+1
n+1个控制点
P
i
(
i
=
0
,
1
,
⋯
,
n
)
P_i(i=0,1,\cdots,n)
Pi(i=0,1,⋯,n)
-
n
+
k
+
1
(
2
⩽
k
⩽
n
+
1
)
n+k+1(2\leqslant k\leqslant n+1)
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}
构造多项式函数
P
(
u
)
=
∑
i
=
0
n
P
i
B
i
,
k
(
u
)
,
u
∈
[
u
k
−
1
,
u
n
+
1
]
(1)
\mathcal{P}(u)=\sum_{i=0}^{n}P_i B_{i,k}(u), \quad u\in[u_{k-1},u_{n+1}] \tag{1}
P(u)=i=0∑nPiBi,k(u),u∈[uk−1,un+1](1)来逼近曲线。其中
B
i
,
k
(
u
)
B_{i,k}(u)
Bi,k(u)称为
k
k
k阶(
k
−
1
k-1
k−1次)B样条基函数,它是由DeBoor-Cox递推公式来定义的:
B
i
,
1
(
u
)
=
{
1
,
u
i
⩽
u
<
u
i
+
1
0
,
o
t
h
e
r
w
i
s
e
,
B
i
,
k
(
u
)
=
u
−
u
i
u
i
+
k
−
1
−
u
i
B
i
,
k
−
1
(
u
)
+
u
i
+
k
−
u
u
i
+
k
−
u
i
+
1
B
i
+
1
,
k
−
1
(
u
)
。
(2)
\begin{alignedat}{3} & B_{i,1}(u) = \begin{cases} 1, & u_i\leqslant u < u_{i+1} \\ 0, & otherwise \end{cases} \thinspace, \\ & 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) \thinspace。 \end{alignedat} \tag{2}
Bi,1(u)={1,0,ui⩽u<ui+1otherwise,Bi,k(u)=ui+k−1−uiu−uiBi,k−1(u)+ui+k−ui+1ui+k−uBi+1,k−1(u)。(2)递推式中若遇到
0
0
0除以
0
0
0,则除法结果取值为
0
0
0。
B
i
,
k
(
u
)
B_{i,k}(u)
Bi,k(u)可以看作是控制点
P
i
P_i
Pi在
u
u
u下的权值,因为B样条基函数具有权性:
∑
i
=
0
n
B
i
,
k
(
u
)
≡
1
,
∀
u
。
(3)
\sum_{i=0}^{n} B_{i,k}(u) \equiv 1, \quad \forall u \thinspace。 \tag{3}
i=0∑nBi,k(u)≡1,∀u。(3)
对于某个
u
∈
[
u
j
,
u
j
+
1
)
u\in [u_j, u_{j+1})
u∈[uj,uj+1),运用公式(2)递推
B
i
,
k
(
u
)
B_{i,k}(u)
Bi,k(u),使得最终的结果中只含
B
x
,
1
(
u
)
B_{x,1}(u)
Bx,1(u),我们可以推出
x
=
i
,
i
+
1
,
⋯
,
i
+
k
−
1
x=i,i+1,\cdots,i+k-1
x=i,i+1,⋯,i+k−1。又根据
B
x
,
1
(
u
)
B_{x,1}(u)
Bx,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
)
B_{i,k}(u)
Bi,k(u)有可能不为0。
举个例子来解释一下。以
i
=
7
,
k
=
3
i=7,k=3
i=7,k=3为例,运用公式(2)递推
B
7
,
3
(
u
)
B_{7,3}(u)
B7,3(u):
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
)
\begin{aligned} B_{7,3}(u) &= c_0B_{7,2}(u)+c_1B_{8,2}(u) \\ &= c_2B_{7,1}(u)+c_3B_{8,1}(u)+c_4B_{9,1}(u) \end{aligned}
B7,3(u)=c0B7,2(u)+c1B8,2(u)=c2B7,1(u)+c3B8,1(u)+c4B9,1(u)(其中
c
0
,
c
1
,
⋯
,
c
4
c_0,c_1,\cdots,c_4
c0,c1,⋯,c4均为系数),我们发现最终的结果中只含
B
7
,
1
(
u
)
,
B
8
,
1
(
u
)
,
B
9
,
1
(
u
)
B_{7,1}(u), B_{8,1}(u),B_{9,1}(u)
B7,1(u),B8,1(u),B9,1(u)。对于某个
u
∈
[
u
j
,
u
j
+
1
)
u\in [u_j, u_{j+1})
u∈[uj,uj+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
)
B_{7,1}(u), B_{8,1}(u),B_{9,1}(u)
B7,1(u),B8,1(u),B9,1(u)不全为0,从而使得
B
7
,
3
(
u
)
B_{7,3}(u)
B7,3(u)有可能不为0。
算法原理
对于某个
u
∈
[
u
j
,
u
j
+
1
]
u\in [u_j, u_{j+1}]
u∈[uj,uj+1],我们定义
d
i
,
0
=
P
i
,
i
=
j
−
(
k
−
1
)
,
⋯
,
j
,
(4)
d_{i,0} = P_i,\quad i = j-(k-1),\cdots,j \thinspace,\tag{4}
di,0=Pi,i=j−(k−1),⋯,j,(4)
α
i
,
r
=
u
−
u
i
u
i
+
k
−
r
−
u
i
,
d
i
,
r
=
(
1
−
α
i
,
r
)
d
i
−
1
,
r
−
1
+
α
i
,
r
d
i
,
r
−
1
(
r
=
1
,
2
,
⋯
,
k
−
1
i
=
j
−
(
k
−
1
)
+
r
,
⋯
,
j
)
。
(5)
\begin{alignedat}{3} \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{alignedat} \left(\begin{array}{l} r = 1,2,\cdots,k-1 \\ i = j-(k-1)+r,\cdots,j \end{array}\right) \thinspace。 \tag{5}
αi,rdi,r=ui+k−r−uiu−ui,=(1−αi,r)di−1,r−1+αi,rdi,r−1(r=1,2,⋯,k−1i=j−(k−1)+r,⋯,j)。(5)
这里对公式(5)中的
i
=
j
−
(
k
−
1
)
+
r
,
⋯
,
j
i = j-(k-1)+r,\cdots,j
i=j−(k−1)+r,⋯,j做出解释:递推
d
i
,
r
d_{i,r}
di,r使最终结果中只含
d
x
,
0
d_{x,0}
dx,0,得到
x
=
i
−
r
,
⋯
,
i
x=i-r,\cdots,i
x=i−r,⋯,i。由
d
x
,
0
d_{x,0}
dx,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}。
可以证明
P
(
u
)
=
d
j
,
k
−
1
。
(6)
\mathcal{P}(u) = d_{j,k-1} \thinspace。 \tag{6}
P(u)=dj,k−1。(6)
于是得到算法逻辑如下(假定输入只提供了
n
+
1
n+1
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]
[0,1]。
- 将
[
0
,
1
]
[0,1]
[0,1]均匀分隔为
n
+
k
n+k
n+k段(控制点个数=
n
+
1
n+1
n+1,曲线阶数=
k
k
k)。
- 分隔之后产生了
n
+
k
+
1
n+k+1
n+k+1个分隔点,依次记为
u
0
,
u
1
,
⋯
,
u
n
+
k
u_0,u_1,\cdots,u_{n+k}
u0,u1,⋯,un+k。于是有
u
j
=
j
/
(
n
+
k
)
u_j=j/(n+k)
uj=j/(n+k)。
- 在区间
[
u
k
−
1
,
u
n
+
1
]
[u_{k-1},u_{n+1}]
[uk−1,un+1]上大量采样。对于每个样本
u
u
u,计算
P
(
u
)
\mathcal{P}(u)
P(u)。
- 根据公式(6),我们可以通过计算
d
j
,
k
−
1
d_{j,k-1}
dj,k−1来计算
P
(
u
)
\mathcal{P}(u)
P(u)。
-
d
j
,
k
−
1
d_{j,k-1}
dj,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
补充:
P
(
u
)
\mathcal{P}(u)
P(u)定义域的推导
在公式(1)中,我们直接定义变量
u
u
u的取值范围为
[
u
k
−
1
,
u
n
+
1
]
[u_{k-1},u_{n+1}]
[uk−1,un+1]。这一节中我们将详细说明其原理。
根据递推公式(2),我们可以推知:
-
B
i
,
1
B_{i,1}
Bi,1涉及1个区间(
[
u
i
,
u
i
+
1
]
[u_i,u_{i+1}]
[ui,ui+1])2个控制点(
u
i
u_i
ui和
u
i
+
1
u_{i+1}
ui+1)
-
B
i
,
2
B_{i,2}
Bi,2由
B
i
,
1
B_{i,1}
Bi,1和
B
i
+
1
,
1
B_{i+1,1}
Bi+1,1组成,涉及2个区间3个控制点
-
B
i
,
3
B_{i,3}
Bi,3由
B
i
,
2
B_{i,2}
Bi,2和
B
i
+
1
,
2
B_{i+1,2}
Bi+1,2组成,涉及3个区间4个控制点
-
⋯
\cdots
⋯
-
B
i
,
k
B_{i,k}
Bi,k涉及
k
k
k个区间
k
+
1
k+1
k+1个控制点
从而可以定义
B
i
,
k
B_{i,k}
Bi,k的支承区间为
[
u
i
,
u
i
+
k
]
[u_i,u_{i+k}]
[ui,ui+k]。
举个例子,如上图,其中
n
=
4
,
k
=
4
n=4,k=4
n=4,k=4。根据支承区间的定义,我们有:
-
P
0
B
0
,
4
P_0B_{0,4}
P0B0,4涉及节点
u
0
u_0
u0到
u
4
u_4
u4
-
P
1
B
1
,
4
P_1B_{1,4}
P1B1,4涉及节点
u
1
u_1
u1到
u
5
u_5
u5
-
⋯
\cdots
⋯
-
P
4
B
4
,
4
P_4B_{4,4}
P4B4,4涉及节点
u
4
u_4
u4到
u
8
u_8
u8
对于每个区间
[
u
i
,
u
i
+
1
]
[u_i,u_{i+1}]
[ui,ui+1],它至多包含于
k
k
k个不同的支承区间中,也就是说,它至多与
k
k
k个控制点相关。下图形象地展示了这种关系。
我们的目标是,对于
P
(
u
)
\mathcal{P}(u)
P(u)的定义域中的每个区间,要使它与尽可能多的控制点相关。于是我们最终选择了
[
u
k
−
1
,
u
n
+
1
]
[u_{k-1},u_{n+1}]
[uk−1,un+1]。可以证明,
∀
[
u
i
,
u
i
+
1
]
⊆
[
u
k
−
1
,
u
n
+
1
]
\forall [u_i,u_{i+1}]\subseteq [u_{k-1},u_{n+1}]
∀[ui,ui+1]⊆[uk−1,un+1],
[
u
i
,
u
i
+
1
]
[u_i,u_{i+1}]
[ui,ui+1]都与
k
k
k个控制点相关。
更多细节
由于篇幅原因,文中有些证明并未详细展开叙述。更多的证明细节请参阅相关论文。例如C. DE BOOR, On calculating with B-splines, 1972等。