一、B 样条曲线的定义和性质
p
p
p 次 B 样条曲线的定义为
C
(
u
)
=
∑
i
=
0
n
N
i
,
p
(
u
)
P
i
,
a
≤
u
≤
b
(1)
\pmb C(u)=\sum_{i=0}^nN_{i,p}(u)\pmb P_i,\quad a\leq u\leq b\tag{1}
C(u)=i=0∑nNi,p(u)Pi,a≤u≤b(1)
这里
P
i
\pmb P_i
Pi 是曲线的控制点,
{
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
)
U=\bigg(\underbrace{a,\cdots,a}_{p+1},u_{p+1},\cdots,u_{m-p-1},\underbrace{b,\cdots,b}_{p+1}\bigg)
U=(p+1
a,⋯,a,up+1,⋯,um−p−1,p+1
b,⋯,b)
(包含 m + 1 m+1 m+1 个节点)上的 p p p 次 B 样条基函数。除非特别声明,通常取 a = 0 , b = 1 a=0,b=1 a=0,b=1。由 P i \pmb P_i Pi 构成的多边形称为控制多边形。
对于固定的 u u u 值,计算 B 样条曲线上的对应点需要三步:
- 找到 u u u 所在的节点区间;
- 计算非零的基函数;
- 将非零基函数的值与相应的控制点相乘,再求和。
例1.
U
=
{
0
,
0
,
0
,
1
,
2
,
3
,
4
,
4
,
5
,
5
,
5
}
,
u
=
5
2
,
p
=
2
\quad 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, p=2
此时,
u
∈
[
u
4
,
u
5
)
u\in[u_4,u_5)
u∈[u4,u5),则
N
2
,
2
(
5
2
)
=
1
8
,
N
3
,
2
(
5
2
)
=
6
8
,
N
4
,
2
(
5
2
)
=
1
8
N_{2,2}\Big(\dfrac{5}{2}\Big)=\frac{1}{8},\quad N_{3,2}\Big(\dfrac{5}{2}\Big)=\frac{6}{8},\quad N_{4,2}\Big(\dfrac{5}{2}\Big)=\frac{1}{8}
N2,2(25)=81,N3,2(25)=86,N4,2(25)=81
与相应的控制点相乘并求和得
C
(
5
2
)
=
1
8
P
2
+
6
8
P
3
+
1
8
P
4
\pmb C\Big(\dfrac{5}{2}\Big)=\frac{1}{8}\pmb P_2+\frac{6}{8}\pmb P_3+\frac{1}{8}\pmb P_4
C(25)=81P2+86P3+81P4
具体算法如下:
def CurvePoint(p, U, P, u):
"""
计算 B 样条曲线上的点
:param p: 基函数次数
:param U: 节点矢量
:param P: 控制点
:param u: 固定值
:return: 固定值 u 在 B 样条曲线上的对应值
"""
span = FindSpan(p, u, U)
N = BasisFuns(span, u, p, U)
C = 0.0
for i in range(p+1):
C = C + N[i] * P[span - p + i]
return C
二、B 样条曲线的导矢
令
C
(
k
)
(
u
)
\pmb C^{(k)}(u)
C(k)(u) 表示
C
(
u
)
\pmb C(u)
C(u) 的
k
k
k 阶导矢。对于固定的
u
u
u,我们可以通过计算基函数的
k
k
k 阶导数得到
C
(
k
)
(
u
)
\pmb C^{(k)}(u)
C(k)(u)。特别地,有
C
(
k
)
(
u
)
=
∑
i
=
0
n
N
i
,
p
(
k
)
(
u
)
P
i
\pmb C^{(k)}(u)=\sum_{i=0}^nN_{i,p}^{(k)}(u)\pmb P_i
C(k)(u)=i=0∑nNi,p(k)(u)Pi
例2.
p
=
2
,
U
=
{
0
,
0
,
0
,
1
,
2
,
3
,
4
,
4
,
5
,
5
,
5
}
,
u
=
5
2
\quad p=2,U=\{0,0,0,1,2,3,4,4,5,5,5\},u=\frac{5}{2}
p=2,U={0,0,0,1,2,3,4,4,5,5,5},u=25
N
2
,
2
′
(
5
2
)
=
0
−
2
3
−
1
1
2
=
−
1
2
N
3
,
2
′
(
5
2
)
=
2
3
−
1
1
2
−
2
4
−
2
1
2
=
0
N
4
,
2
′
(
5
2
)
=
2
4
−
2
1
2
−
0
=
1
2
\begin{aligned} &N^\prime_{2,2}\Big(\frac{5}{2}\Big)=0-\frac{2}{3-1}\frac{1}{2}=-\frac{1}{2}\\[2ex] &N^\prime_{3,2}\Big(\frac{5}{2}\Big)=\frac{2}{3-1}\frac{1}{2}-\frac{2}{4-2}\frac{1}{2}=0\\[2ex] &N^\prime_{4,2}\Big(\frac{5}{2}\Big)=\frac{2}{4-2}\frac{1}{2}-0=\frac{1}{2} \end{aligned}
N2,2′(25)=0−3−1221=−21N3,2′(25)=3−1221−4−2221=0N4,2′(25)=4−2221−0=21
于是有
C
′
(
5
2
)
=
−
1
2
P
2
+
1
2
P
4
\pmb C^\prime\Big(\frac{5}{2}\Big)=-\frac{1}{2}\pmb P_2+\frac{1}{2}\pmb P_4
C′(25)=−21P2+21P4
以下给出对于固定的 u u u,计算 B 样条曲线上对应的点及其直到 d d d 阶(包括 d d d 阶)导矢的算法。我们允许 d > p d>p d>p,虽然(对非有理曲线而言)这些导矢都为零;但是,这些导矢对有理曲线是必要的。算法的输入是 u , d u,d u,d 以及下面的一些数据定义的 B 样条曲线
n
n
n:控制点的最大下标(控制点的个数是
n
+
1
n+1
n+1);
p
p
p:曲线的次数;
U
U
U:节点矢量;
P
P
P:控制点数组。
算法的输出是数组 CK
def CurveDerivsAlg(n, p, U, P, u, d):
"""
计算固定的 u 时 B 样条曲线上对应的点及其 d 阶导矢
:param n: 控制点的最大下标(控制点的个数是 n+1)
:param p: 曲线的次数
:param U: 节点矢量
:param P: 控制点数组
:param u: 固定点
:param d: 导矢阶数
:return: 数组 CK
"""
du = min(d, p)
CK = np.zeros(d+1)
for k in range(p+1, d+1):
CK[k] = 0.0
span = FindSpan(p, u, U)
nders = DersBasisFuns(span, u, p, n, U)
for k in range(du+1):
CK[k] = 0.0
for j in range(p+1):
CK[k] = CK[k] + nders[k, j] * P[span-p+j]
return CK