样条曲线(Spline Curves)
前言: 关于样条曲线的一些总结和使用,其中包含贝塞尔曲线的介绍。内容实际是为了解决在给定控制点的条件下,如何确定一个光滑曲线的问题。样条曲线被广泛应用于模型的几何重构平滑中。
下启:样条曲线(下)之插值问题
文章目录
0. 基本概念
样条曲线(Spline Curves): 是给定一系列控制点而得到的一条曲线,曲线形状由这些点控制。一般分为插值样条和拟合样条。
插值:在原有数据点上进行填充生成曲线,曲线必经过原有数据点。
拟合:依据原有数据点,通过参数调整设置,使得生成曲线与原有点差距最小(最小二乘),因此曲线未必会经过原有数据点。
1. 起源
样条曲线起源于一个常见问题,即已知若干点的条件下,如何得到通过这些点的一条光滑曲线?
一个简单且行之有效的方法是,把这些点作为限制点,然后在这些限制点中放置一条具有弹性的金属片,最后金属片绕过这些点后的最终状态即为所需曲线。而最终得到的形状曲线,就是样条曲线。这也是该名字的由来,其中金属片就是样条,形成的曲线就是样条曲线。如下图
该想法虽然巧妙,但显然不具有推广性。因此问题就出来了,如何将其抽象出一个数学模型,从而在已知控制点条件下,仅仅通过数学公式从而获得平滑的样条。
再简化一些,问题便是,如何在两点间进行插值,整体上获得一个平滑曲线。
贝塞尔曲线实际上就是一个样条曲线。可以通过多个控制点,插值获得一个平滑曲线。何为贝塞尔曲线?
2. 贝塞尔曲线(Bézier curve)
1.1 直线段表示(2 points)
从最简单的问题入手,假设有两点 P 0 , P 1 P_0,P_1 P0,P1(端点),如何根据两点确定一条(曲)线?
最简单的问题也有最简单的解答,确定的这条线就是连接 P 0 , P 1 P_0,P_1 P0,P1 两点的线段。从计算机的角度而言,本没有连续的线,因为计算机数据本身就是离散的,因此计算机中线的描述也只是由无数个点组成。
那么,如何用计算机绘制出这条线呢,确切说,如何表达这条线段上任何位置的坐标呢?(其实本质上就是在两点间插值)
解答同样简单,线段上任意位置点坐标可以用一个比例参数 t 来表示,表达式为
P
x
=
P
0
+
(
P
1
−
P
0
)
t
=
(
1
−
t
)
P
0
+
t
P
1
,
t
∈
[
0
,
1
]
h
e
r
e
,
P
0
P
x
P
0
P
1
=
t
P_x = P_0 + (P_1-P_0)t = (1-t)P_0 + tP_1,\ \ \ \ t\in [0,1] \\ here,\ \ \frac{P_0 P_x}{P_0 P_1} = t
Px=P0+(P1−P0)t=(1−t)P0+tP1, t∈[0,1]here, P0P1P0Px=t
实际上这就是一次的贝塞尔曲线的表达式。之所以为称为一次,是因为式中 t 的最高幂次为 1。对于式子的不同理解可以为,
(
1
−
t
)
(1-t)
(1−t) 为
P
0
P_0
P0 权重,t 为
P
1
P_1
P1 权重,权重越大,越靠近对应的点。
1.2 抛物线三切线定理(3 points)
问题来到三个点,即如果知道三个点 P 0 , P 1 , P 2 P_0,P_1,P_2 P0,P1,P2,且 P 0 , P 2 P_0,P_2 P0,P2 为端点,那么如何确定一条曲线呢?(理解贝塞尔曲线时,可以把 P 1 P_1 P1看作是控制点,我们要的是在控制点下,对 P 0 , P 2 P_0,P_2 P0,P2 两点间的插值)
最简单的想法是,直接连接 P 0 , P 1 P_0,P_1 P0,P1 和 P 1 , P 2 P_1,P_2 P1,P2,得到两个线段就可以确定一个曲线,但这明显不是我们想要的,它不够平滑,往往一阶导还不连续(即不满足c1连续,下有详细叙述)。
第二个想法,可以考虑借用二次曲线,但如何在知道三点的情况下确定一个二次曲线呢?
最简单做法是列三个方程组,解方程,然后得到一个二次函数,也即确定了二次曲线,且通过已知的三点。同样,递推到更多已知点,也可以用相似的方法。对于 n 个点,可以使用 n-1 阶次的函数来确定一个唯一的曲线。这种方法也即是常用的多项式插值。前面,提到插值的特点在于,所有已知点都在最终确定的曲线上。
如何使用前述把 P 1 P_1 P1当作控制点的思路来解决这个问题呢?即没必要所有点都处在最终确定的曲线上,是否可以参考上面两点情况下,使用权重或是比例参数 t 来确定一条曲线呢?
现在的目的就是,使用三点确定任意一点
P
x
P_x
Px ,且最终由无数个
P
x
P_x
Px 点形成的曲线连续可导(c1连续)。依据上面两点坐标下确定曲线任意点的思路,异想天开一下。先将三个点划分为两组,即
P
0
,
P
1
P_0,P_1
P0,P1和
P
1
,
P
2
P_1,P_2
P1,P2。对于两组点,分别使用上面结论,那么在参数 t 下,两组分别确定了两个任意点,记为
P
i
,
P
j
P_i,P_j
Pi,Pj ,也即
P
i
=
(
1
−
t
)
P
0
+
P
1
P
j
=
(
1
−
t
)
P
1
+
p
2
P_i = (1-t)P_0 + P_1 \\ P_j = (1-t)P_1 + p_2
Pi=(1−t)P0+P1Pj=(1−t)P1+p2
可参考下图中的
P
i
,
P
j
P_i,P_j
Pi,Pj 。
但我们想要的是由
P
0
,
P
1
,
P
2
P_0,P_1,P_2
P0,P1,P2 三个点确定的唯一点
P
x
P_x
Px,而现在我们得到了两个点,那再使用一次 上面公式不就可以了。也即
P
x
=
(
1
−
t
)
P
i
+
t
P
j
P_x = (1-t)P_i + tP_j
Px=(1−t)Pi+tPj
问题解决了。
幸运的是,对于
t
∈
[
0
,
1
]
t\in [0,1]
t∈[0,1] 最终确定的线是可导、平滑的。想法有点异想天开,但还是要总结一下,在基于此想法下,用到了三次上面两点下任意点的求法,且有
P
0
P
i
P
0
P
1
=
P
1
P
j
P
1
P
2
=
P
j
P
x
P
i
P
j
=
t
\frac{P_0 P_i}{P_0 P_1} = \frac{P_1 P_j}{P_1 P_2} = \frac{P_j P_x}{P_i P_j} = t
P0P1P0Pi=P1P2P1Pj=PiPjPjPx=t
最终结果
P
x
=
(
1
−
t
)
P
i
+
t
P
j
=
(
1
−
t
)
[
(
1
−
t
)
P
0
+
t
P
1
]
+
t
[
(
1
−
t
)
P
1
+
t
P
2
]
=
(
1
−
t
)
2
P
0
+
2
t
(
1
−
t
)
P
1
+
t
2
P
2
P_x = (1-t)P_i + tP_j = (1-t)[(1-t)P_0 + tP_1] + t[(1-t)P_1 + tP_2] = (1-t)^2P_0 + 2t(1-t)P_1 + t^2 P_2
Px=(1−t)Pi+tPj=(1−t)[(1−t)P0+tP1]+t[(1−t)P1+tP2]=(1−t)2P0+2t(1−t)P1+t2P2
实际上,上述描述的就是抛物线的三切线定理(性质)。即对于一个抛物线,经过线上两点
P
0
,
P
2
P_0,P_2
P0,P2 的两个切线相交于
P
1
P_1
P1 点,经过在
P
0
,
P
2
P_0,P_2
P0,P2 间抛物线一点
P
x
P_x
Px 的切线,与上述确定的两条切线交于
P
i
,
P
j
P_i,P_j
Pi,Pj 。则有
P
0
P
i
P
0
P
i
=
P
1
P
j
P
1
P
2
=
P
j
P
x
P
i
P
j
\frac{P_0 P_i}{P_0 P_i} = \frac{P_1 P_j}{P_1 P_2} = \frac{P_j P_x}{P_i P_j}
P0PiP0Pi=P1P2P1Pj=PiPjPjPx
如果令该比值为 t ,三个方程,三个未知点。最终可以同样反推到上述结果
P
x
=
(
1
−
t
)
2
P
0
+
2
t
(
1
−
t
)
P
1
+
t
2
P
2
P_x = (1-t)^2P_0 + 2t(1-t)P_1 + t^2 P_2
Px=(1−t)2P0+2t(1−t)P1+t2P2
该结果,就是二次的贝塞尔曲线。同样二次的含义在于,t 的最高幂次为 2
前面"异想天开"的想法,在于说明,贝塞尔曲线是递推的,或者说是可以递推的。二阶贝塞尔曲线是两个一阶贝塞尔曲线的线性组合。同样的思维,三阶贝塞尔曲线由两个二阶贝塞尔曲线线性组合。因此我们可以确定三阶、四阶,以至于更高阶。
值得说明的是,photoshop中的钢笔工具就是应用的三次贝塞尔曲线。
1.3 通用公式
假设一共有 n + 1 n+1 n+1 个点 P 0 , P 1 , ⋯ , P n P_0,P_1,\cdots,P_n P0,P1,⋯,Pn ,这 n + 1 n+1 n+1 个点确定了 n 次的贝塞尔曲线。
n 阶贝塞尔曲线
B
n
(
t
)
B^{n}(t)
Bn(t) 可以由前 n 个点决定的 n-1 次贝塞尔曲线
B
n
−
1
(
t
∣
P
0
,
⋯
,
P
n
−
1
)
B^{n-1}(t|P_0,\cdots,P_{n-1})
Bn−1(t∣P0,⋯,Pn−1) 与 后 n 个点决定的 n-1 次贝塞尔曲线
B
n
−
1
(
t
∣
P
1
,
⋯
,
P
n
)
B^{n-1}(t|P_1,\cdots,P_n)
Bn−1(t∣P1,⋯,Pn) 线性组合递推而来,即
B
n
(
t
∣
P
0
,
P
1
,
⋯
,
P
n
)
=
(
1
−
t
)
B
n
−
1
(
t
∣
P
0
,
P
1
,
⋯
,
P
n
−
1
)
+
t
B
n
−
1
(
t
∣
P
1
,
P
2
,
⋯
,
P
n
)
B^{n}(t|P_0,P_1,\cdots,P_n) = (1-t)B^{n-1}(t|P_0,P_1,\cdots,P_{n-1}) + t B^{n-1}(t|P_1,P_2,\cdots,P_n)
Bn(t∣P0,P1,⋯,Pn)=(1−t)Bn−1(t∣P0,P1,⋯,Pn−1)+tBn−1(t∣P1,P2,⋯,Pn)
且一次贝塞尔曲线,即为由两点决定的线段,也即
B
1
(
t
∣
P
0
,
P
1
)
=
(
1
−
t
)
P
0
+
t
P
1
B^1(t|P_0,P_1) = (1-t) P_0 + tP_1
B1(t∣P0,P1)=(1−t)P0+tP1
可以得到 n 次贝塞尔曲线的表达通式为
B
(
t
)
=
∑
i
=
0
n
C
n
i
(
1
−
t
)
n
−
i
t
i
P
i
B(t) = \sum\limits_{i=0}^n C_n^i(1-t)^{n-i}t^i P_i
B(t)=i=0∑nCni(1−t)n−itiPi
1.4 代码实现
python
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import comb
def inputPoints():
controlPoints = []
num = 1
while True:
print('\nenter %dst control point:'%num)
x = input('x:')
y = input('y:')
#z = input('z:')
print('Point:[%f,%f]'%(float(x),float(y)))
i = input('Are you sure?(y or n)')
if i=='y' or i=='Y':
controlPoints.append([float(x),float(y)])
inp = input('continue entering points or not?(y or n)')
num = num + 1
if inp == 'n':
break
else:
continue
#print(controlPoints)
return controlPoints
def getInterpolationPoints(controlPoints, tList):
n = len(controlPoints)-1
interPoints = []
for t in tList:
Bt = np.zeros(2, np.float64)
for i in range(len(controlPoints)):
Bt = Bt + comb(n,i) * np.power(1-t,n-i) * np.power(t,i) * np.array(controlPoints[i])
interPoints.append(list(Bt))
return interPoints
if __name__ == '__main__':
# points = inputPoints()
points = [[1,1],[3,4],[5,5],[7,2]]
tList = np.linspace(0,1,50)
interPointsList = getInterpolationPoints(points, tList)
x = np.array(interPointsList)[:,0]
y = np.array(interPointsList)[:,1]
plt.plot(x,y,color='b')
plt.scatter(np.array(points)[:,0],np.array(points)[:,1],color='r')
plt.show()
结果
说明
- 实现使用了二维点,应该可以更高维
- 结果展示了三次贝塞尔曲线,由四个控制点(红色)生成
3. 曲线连续性描述
上面提到了c1连续,在计算机图像学中,有详细描述
为理解说明,假设有两个曲线 f ( x ) , x ∈ [ a , b ] f(x),x\in[a,b] f(x),x∈[a,b] 和 g ( x ) , x ∈ [ b , c ] g(x),x\in[b,c] g(x),x∈[b,c],
c 0 c^0 c0连续:即函数连续,不存在间断点。实例中 f ( b ) = g ( b ) f(b) = g(b) f(b)=g(b)
c 1 c^1 c1连续:函数的一阶可导且导数连续。实例中 f ′ ( b ) = g ′ ( b ) f'(b) = g'(b) f′(b)=g′(b)
c 2 c^2 c2连续:函数的二阶可导且导数连续。实例中 f ′ ′ ( b ) = g ′ ′ ( b ) f''(b) = g''(b) f′′(b)=g′′(b)
c n c^n cn连续:函数的 n n n阶可导且导数连续。实例中 f ( n ) ( b ) = g ( n ) ( b ) f^{(n)}(b) = g^{(n)}(b) f(n)(b)=g(n)(b)
更详细说明和实例可以参考维基:光滑函数
4. B样条基函数(B-Spline Basis Function)
4.1 多段贝塞尔曲线
在剖物线三切线定理中,最初想法是三点直接相连,得到曲线虽然 c 0 c^0 c0 连续,但不满足 c 1 c^1 c1 连续。为了解决这个问题,进而使用了二次的贝塞尔曲线。
问题是,对于一个复杂弯曲的曲线,我们希望使用一个贝塞尔曲线来插值获得目标曲线,为此我们可以通过增加控制点来进行插值。但目标曲线越复杂,需要的控制点就越多,而 贝塞尔曲线幂次 = 控制点个数 - 1,即需要的计算也越复杂。该方法虽然可行,但是不明智的,低效率的。另外对于单一的贝塞尔曲线,改变其中一个控制点,那么整条曲线都会随之改变。
因此对于复杂曲线,一般做法是,使用三次贝塞尔曲线(常用次)一段一段地拼接成目标曲线,正如 Ps 或 Ai 中使用钢笔工具画出物体轮廓所做的那样。 如果使用这种方法,确保最终整体曲线 c 1 c^1 c1 连续的条件是 在连接点出两侧的斜率相等,即连接点和其两侧控制点共线。
下图是由两个三次贝塞尔曲线组成的曲线
假设从左到右点依次为 P 0 , P 1 , ⋅ , P 7 P_0,P_1,\cdot,P_7 P0,P1,⋅,P7 ,确保两段曲线拼接起来 c 1 c^1 c1 连续的条件就是 P 2 , P 3 , P 4 P_2,P_3,P_4 P2,P3,P4 三点共线。
4.2 从贝塞尔曲线到B样条基函数
4.2.1 贝塞尔曲线
回到最初问题上,通过一系列点,获取一条光滑的曲线。也即通过这些控制点,生成一系列点坐标,这些点坐标形成光滑曲线。
对于贝塞尔曲线上点的生成,是通过如下方程函数
B
(
t
)
=
∑
i
=
0
n
C
n
i
(
1
−
t
)
n
−
i
t
i
P
i
,
t
∈
[
0
,
1
]
B(t) = \sum\limits_{i=0}^n C_n^i(1-t)^{n-i}t^i P_i\ ,\ \ \ \ t\in[0,1]
B(t)=i=0∑nCni(1−t)n−itiPi , t∈[0,1]
或者展开写成这样
B
(
t
)
=
W
t
,
n
0
P
0
+
W
t
,
n
1
P
1
+
⋯
+
W
t
,
n
n
P
n
B(t) = W_{t,n}^0 P_0 + W_{t,n}^1 P_1 + \cdots + W_{t,n}^n P_n
B(t)=Wt,n0P0+Wt,n1P1+⋯+Wt,nnPn
其中
W
t
,
n
0
W_{t,n}^0
Wt,n0 为
P
0
P_0
P0 前系数,是最高幂次为 n 的关于 t 的多项式。当 t 确定后,该值就为定值。
因此整个式子可以理解为 B(t) 插值点是这 n+1 个点施加各自的权重 W 后累加得到的。这也是为什么改变其中一个控制点,整个贝塞尔曲线都会受到影响。
其实对于样条曲线的生成,思路就是这样的,最终曲线的生成,就对于各个控制点施加权重即可。在贝塞尔曲线中,该权重是关于 t 的函数即
W
i
=
C
n
i
(
1
−
t
)
n
−
i
t
i
W^i = C_n^i(1-t)^{n-i}t^i
Wi=Cni(1−t)n−iti
在B样条中,只不过这个权重设置更复杂点,下面一点点剖析其B样条曲线形成的原理。
4.2.2 B样条
上面提到,生成曲线,本质上就是在控制点前增加一个权重,然后累加即可。那么B样条中权重是如何计算的呢?
先介绍几个概念,容易混淆的概念,不太理解其来由也没关系。下面会一点点指出其用处。
- 控制点:也就是控制曲线的点,等价于贝塞尔函数的控制点,通过控制点可以控制曲线形状。假设有 n+1 个控制点 P 0 , P 1 , P 2 , ⋯ , P n P_0,P_1,P_2,\cdots,P_{n} P0,P1,P2,⋯,Pn
- 节点:这个跟控制点没有关系,而是人为地将目标曲线分为若干个部分,其目的就是尽量使得各个部分有所影响但也有一定独立性,这也是为什么B样条中,有时一个控制点的改变,不会很大影响到整条曲线,而只影响到局部的原因,这是区别于贝塞尔曲线的一点。节点划分影响到权重计算,可以预想到的是,实现局部性的影响的原理应该是在生成某区间内的点时,某些控制点前的权重值会为0,即对该点没有贡献,所以才有上述特点。事实上,就是如此的。先了解有这个概念即可。假设我们划分了 m+1 个节点 t 0 , t 1 , ⋯ , t m t_0,t_1,\cdots,t_m t0,t1,⋯,tm,将曲线分成了 m 段
- 次与阶 :次的概念就是贝塞尔中次的概念,即权重中 t 的最高幂次。而 阶=次+1。这里只需要知道次这个概念即可。假设我们用 k 表示次,即 k 次。
下面用一个实例来说明
如上图,对应为
控制点有 n+1 = 5 个,n=4,即 P 0 , P 1 , P 2 , P 3 , P 4 P_0,P_1,P_2,P_3,P_4 P0,P1,P2,P3,P4
节点规定为 m+1 = 10 个,m=9,即 t 0 , t 1 , ⋯ , t 9 t_0,t_1,\cdots,t_9 t0,t1,⋯,t9 ,该节点将要生成的目标曲线分为了 9 份,这里的 t 取值一般为 0-1的一系列非递减数。如 0 代表起始位置,1 代表末位置,0.5 代表一半的位置。 t 0 , t 1 , ⋯ , t 9 t_0,t_1,\cdots,t_9 t0,t1,⋯,t9 组成的序列,叫做节点表。如等分的节点表 { 0 , 1 9 , 2 9 , 3 9 , 4 9 , 5 9 , 6 9 , 7 9 , 8 9 , 1 } \{0,\frac{1}{9},\frac{2}{9},\frac{3}{9},\frac{4}{9},\frac{5}{9},\frac{6}{9},\frac{7}{9},\frac{8}{9},1 \} {0,91,92,93,94,95,96,97,98,1} 。
次为 k 。实际上有个必须要满足的关系式为 m = n + k + 1 m=n+k+1 m=n+k+1 ,即 k = m − n − 1 = 4 k = m-n-1 = 4 k=m−n−1=4 ,先对这个关系式有个印象即可。
下面开始用上面例子说明 B样条 是如何工作的。
我们的目的就是获取最终样条曲线,而根据前述,也即是获得任意点的坐标值,即
t
∈
[
0
,
1
]
t\in[0,1]
t∈[0,1] 时,获得一系列的点,这些点最终组成了目标曲线。这些点的坐标可以表示为 各控制点
P
i
P_i
Pi 与其权重
W
i
W_i
Wi 的乘积再累加,即
B
(
t
)
=
∑
i
=
0
n
W
i
P
i
B(t) = \sum\limits_{i=0}^n W_iP_i
B(t)=i=0∑nWiPi
也即我们获得
W
i
W_i
Wi 即可。
W
i
W_i
Wi 是关于 t 的函数,最高幂次为 k。B样条中通常记为
B
i
,
k
(
t
)
B_{i,k}(t)
Bi,k(t) ,即 表示第 i 点的权重,关于 t 的函数,且最高幂次为 k。而这个权重函数
B
i
,
k
(
t
)
B_{i,k}(t)
Bi,k(t) ,在B样条里叫做 k次B样条基函数。
下面就是关键,权重函数 B i , k ( t ) B_{i,k}(t) Bi,k(t) 是如何取值的?这里就要用到前面的节点表的值了。针对上面实例,我们现在的目标是获取 5 个控制点 P 0 , P 1 , P 2 , P 3 , P 4 P_0,P_1,P_2,P_3,P_4 P0,P1,P2,P3,P4 对应的 5 个权重值 B 0 , 4 , B 1 , 4 , B 2 , 4 , B 3 , 4 , B 4 , 4 B_{0,4},B_{1,4},B_{2,4},B_{3,4},B_{4,4} B0,4,B1,4,B2,4,B3,4,B4,4
B样条是如何做的?先看如下一个列表
我们的目标是获得右侧五个值,这里 k 的取值假设是不知道的(虽然是 4 ),大可当作不知道。
表中 b i , k b_{i,k} bi,k 代表含义跟 B i , k B_{i,k} Bi,k 是一样的,因为不是我们需要的最终值,而是需要求解的中间递推值,所以用小写表示。对于B样条权重,其是根据节点来计算,理论上,m+1个节点,可以得到如表所示的,每个k对应 9(m=9) 个b值,但我们最终只需要5个值,所以就需要对应规则限制。B样条获取5个权重值过程和规则如下(B样条基函数)
k = 0 时,即 0 次时,如果 t ∈ [ t j , t j + 1 ] t\in [t_j,t_{j+1}] t∈[tj,tj+1] ,那么规定 b j , 0 = 1 b_{j,0} = 1 bj,0=1, 其余都为 0。如 t ∈ [ t 2 , t 3 ] t\in[t_2,t_3] t∈[t2,t3] 就有如下结果
k = 1 时,
b
j
,
1
b_{j,1}
bj,1 取值与低一次的相邻两节点相关,即与 0 次同域两端的节点相关。即
b
j
,
1
b_{j,1}
bj,1 求解与
b
j
,
0
,
b
j
+
1
,
0
b_{j,0},b_{j+1,0}
bj,0,bj+1,0 相关,如
b
1
,
1
b_{1,1}
b1,1 就是依据
b
1
,
0
,
b
2
,
0
b_{1,0},b_{2,0}
b1,0,b2,0 求得。求解关系形式为
b
j
,
1
=
A
(
t
)
b
j
,
0
+
B
(
t
)
b
j
+
1
,
0
b_{j,1} = A(t)\ b_{j,0} + B(t)\ b_{j+1,0}
bj,1=A(t) bj,0+B(t) bj+1,0
这里 A(t), B(t) 是关于 t 的一次幂函数。具体为
A
(
t
)
=
t
−
t
j
t
j
+
k
−
t
j
B
(
t
)
=
t
j
+
k
+
1
−
t
t
j
+
k
+
1
−
t
j
+
1
A(t) = \frac{t-t_j}{t_{j+k}-t_j} \\ B(t) = \frac{t_{j+k+1}-t}{t_{j+k+1}-t_{j+1}}
A(t)=tj+k−tjt−tjB(t)=tj+k+1−tj+1tj+k+1−t
这里又涉及到“莫名其秒的下标,但涉及下标的那些值都为定值,此时只要记得是关于 t 的一次函数即可。这里有点类似与贝塞尔函数中递推过程中的前两值的线性组合。
我们用箭头表示上述计算过程,如下
值得说明的两点:
- 在 b j , 1 b_{j,1} bj,1 中,有些值会取 0,因为在 0 次的 b 值有 0 。观察我们可以很容易的值,此时值为 0 的有 b 0 , 1 , b 3 , 1 , b 4 , 1 , b 5 , 1 , b 6 , 1 , b 7 , 1 b_{0,1},b_{3,1},b_{4,1},b_{5,1},b_{6,1},b_{7,1} b0,1,b3,1,b4,1,b5,1,b6,1,b7,1 (注意:这是在 t 取值在 [ t 2 , t 3 ] [t_2,t_3] [t2,t3] 范围的条件下)。而取值为非0的 b 1 , 1 , b 2 , 1 b_{1,1},b_{2,1} b1,1,b2,1 都是关于 t 的一次函数(因为 A ( t ) , B ( t ) A(t),B(t) A(t),B(t) 为 t 的一次函数)。前面提到想要使用局部控制曲线(而非全局)需要某些值为 0,所以是符合预期的。
- b 8 , 1 b_{8,1} b8,1 是无法求解的,因为它需要 b 9 , 0 b_{9,0} b9,0 ,而该值不存在,所以,在 k=1 时,能得到的 b j , 1 b_{j,1} bj,1 减少了 1 (k=1)个,共有8个。同样是令人振奋的结果,因为我们的目标是需要 5 个值(因为有5个控制点),多减少几个不就达到目标了嘛。
k = 2 时,同样的过程,不再以表列出,同样有些 b j , 2 b_{j,2} bj,2 会取值为 0。有些值为非 0, 不同的一点在于,此时系数又乘了一次 A,B 函数,因此那些非 0 值,此时是关于 t 的二次幂函数。另外得到的 b b b 减少了 **2(k=2)**个,还剩 m − k = 9 − 2 = 7 m-k = 9-2 = 7 m−k=9−2=7 个。
k = x时,同样,非0值会是关于 t 的 x 次函数,这也是为什么 k 会是 t 的最高次幂。b 会减少 k 个,还剩 m − k m-k m−k 。
我们最终那个需要 5 个(即 n+1 个)值,那么就需要有
m
−
k
=
n
+
1
m-k = n + 1
m−k=n+1 ,也即
m
=
k
+
n
+
1
m = k + n + 1
m=k+n+1 ,这便是上面说到的
m
,
n
,
k
m,n,k
m,n,k 关系式的由来。所以实例中
k
=
m
−
n
−
1
=
9
−
4
−
1
=
4
k = m-n-1 = 9-4-1 = 4
k=m−n−1=9−4−1=4 ,即当 k=4 时,恰好还剩 5 个 b 值,就是我们需要的B样条基函数B。此时,最高幂次为 4。最终图表示意如下:
该图表表示的是在 t 取值在 [ t 2 , t 3 ] [t_2,t_3] [t2,t3] 范围时,求得的 5 个控制点对应的权重 B j , 4 B_{j,4} Bj,4(基函数),最终结果 B 0 , 4 , B 1 , 4 , B 2 , 4 B_{0,4},B_{1,4},B_{2,4} B0,4,B1,4,B2,4 都为非0值,且是关于 t 的四次函数, B 3 , 4 , B 4 , 4 B_{3,4},B_{4,4} B3,4,B4,4 为 0。 也就是说,本例中, [ t 2 , t 3 ] [t_2,t_3] [t2,t3] 区间的曲线上点不受第四和第五个控制点变化的影响。
这就是所有的 B样条曲线生成的工作原理。
值得注意的几点性质:
- 对于 t ∈ [ t j , t j + 1 ] t\in[t_j,t_{j+1}] t∈[tj,tj+1] 时,最终非零 B j , k B_{j,k} Bj,k为 B j − k , k , ⋯ , B j , k B_{j-k,k},\cdots,B_{j,k} Bj−k,k,⋯,Bj,k (如果 j − k > = 0 j-k>=0 j−k>=0的话),即第 j − k j-k j−k 至第 j j j 个控制点控制影响着 [ t j , t j + 1 ] [t_j,t_{j+1}] [tj,tj+1] 区间内曲线点。可进行局部控制。
- 最终所得曲线结果一般不会经过控制点,因为最终得到的非0权重始终不止一个
- t i t_i ti 为非递减的,可以取重值
- 对于图中,
4.2.3 一般数学描述
将上面过程作为严谨的数学描述
对于 n+1 个控制点
P
0
,
P
1
,
⋯
,
P
n
P_0,P_1,\cdots,P_n
P0,P1,⋯,Pn ,有一个包含 m+1 个节点的列表(或节点向量)
t
0
,
t
1
,
⋯
,
t
m
{t_0,t_1,\cdots,t_{m}}
t0,t1,⋯,tm ,其 k 次B样条曲线表达式为(且m=n+k+1)
P
(
t
)
=
∑
i
=
0
n
B
i
,
k
(
t
)
P
i
P(t) = \sum\limits_{i=0}^n B_{i,k}(t)\ P_i
P(t)=i=0∑nBi,k(t) Pi
其中
B
i
,
k
(
t
)
B_{i,k}(t)
Bi,k(t) 称为 k 次 B 样条基函数,也叫调和函数。且
B
i
,
k
(
t
)
B_{i,k}(t)
Bi,k(t) 满足如下递推式(de Boor递推式)
k
=
0
,
B
i
,
0
(
t
)
=
{
1
,
t
∈
[
t
i
,
t
i
+
1
]
0
,
O
t
h
e
r
w
i
s
e
k
>
0
,
B
i
,
k
(
t
)
=
t
−
t
i
t
i
+
k
−
t
i
B
i
,
k
−
1
(
t
)
+
t
i
+
k
+
1
−
t
t
i
+
k
+
1
−
t
i
+
1
B
i
+
1
,
k
−
1
(
t
)
k = 0,\ \ \ \ B_{i,0}(t) = \left\{ \begin{matrix} 1, \ \ \ \ t\in[t_i,t_i+1] \\ 0, \ \ \ \ \ \ \ \ Otherwise \end{matrix} \right.\\ k > 0,\ \ \ \ B_{i,k}(t) = \frac{t-t_i}{t_{i+k}-t_i} B_{i,k-1}(t) + \frac{t_{i+k+1}-t}{t_{i+k+1}-t_{i+1}} B_{i+1,k-1}(t)
k=0, Bi,0(t)={1, t∈[ti,ti+1]0, Otherwisek>0, Bi,k(t)=ti+k−tit−tiBi,k−1(t)+ti+k+1−ti+1ti+k+1−tBi+1,k−1(t)
5. 分类
- 均匀 B 样条:节点均匀分布
- 准均匀 B 样条:在开始和结束处的节点可重复,中间节点均匀分布
- 非均匀 B 样条:节点非均匀分布
- 分段贝塞尔曲线
- PLUS(挖坑):B样条无法描述圆锥曲线,为解决此问题,产生了非均匀有理B样条(non-uniform rational b-spline, NURBS)
6. 代码实现及验证的一些重要问题
6.1 完整python代码
import numpy as np
import matplotlib.pyplot as plt
# 计算在某一特定t下的 B_{i,k}
def getBt(controlPoints, knots, t):
# calculate m,n,k
m = knots.shape[0]-1
n = controlPoints.shape[0]-1
k = m - n - 1
# initialize B by zeros
B = np.zeros((k+1, m))
# get t region
tStart = 0
for x in range(m+1):
if t==1:
tStart = m-1
if knots[x] > t:
tStart = x-1
break
# calculate B(t)
for _k in range(k+1):
if _k == 0:
B[_k, tStart] = 1
else:
for i in range(m-_k):
if knots[i+_k]-knots[i]== 0:
w1 = 0
else:
w1 = (t-knots[i])/(knots[i+_k]-knots[i])
if knots[i+_k+1]-knots[i+1] == 0:
w2 = 0
else:
w2 = (knots[i+_k+1]-t)/(knots[i+_k+1]-knots[i+1])
B[_k,i] = w1*B[_k-1, i] + w2*B[_k-1, i+1]
return B
# 绘制 B_{i,k}(t)函数
def plotBt(Bt,num, i,k):
print(k,i)
Bt = np.array(Bt)
tt = np.linspace(0,1,num)
yy = [Bt[t,k,i] for t in range(num)]
plt.plot(tt, yy)
# 根据最后一列(最高阶次)的 B(t),即权重,乘以控制点坐标,从而求出曲线上点坐标
def getPt(Bt, controlPoints):
Bt = np.array(Bt)
ptArray = Bt.reshape(-1,1) * controlPoints
pt = ptArray.sum(axis = 0)
return pt
# 绘制出生成的样条曲线: useReg 表示是否使用曲线有效定义域[t_k, t_{m-k}]
def main1(useReg = False):
controlPoints = np.array([[50,50], [100,300], [300,100], [380,200], [400,600]])
knots = np.array([0,1/9,2/9,3/9,4/9,5/9,6/9,7/9,8/9,1])
m = knots.shape[0]-1
n = controlPoints.shape[0]-1
k = m - n - 1
print('n:',n)
print('m:',m)
print('k:',k)
for t in np.linspace(0,1,100):
if useReg and not(t >= knots[k] and t<= knots[n+1]):
continue
Bt = getBt(controlPoints, knots, t)
Pt = getPt(Bt[k, :n+1], controlPoints)
plt.scatter(Pt[0],Pt[1],color='b')
plt.scatter(controlPoints[:,0], controlPoints[:,1],color = 'r')
plt.show()
# 绘制 B_{i,k} 变化图:如果不给定{i,k}则显示所有B{i,k}(t)图像
def main2(i=-1,k=-1):
controlPoints = np.array([[50,50], [100,300], [300,100], [380,200], [400,600]])
knots = np.array([0,1/9,2/9,3/9,4/9,5/9,6/9,7/9,8/9,1])
m = knots.shape[0]-1
n = controlPoints.shape[0]-1
k = m - n - 1
print('n:',n)
print('m:',m)
print('k:',k)
B = []
num = 100 # 离散点数目
for t in np.linspace(0,1,num):
Bt = getBt(controlPoints, knots, t)
B.append(list(Bt))
figure1 = plt.figure('B_{i,k}')
if i==-1:
fig = []
for i in range(n+1):
for k in range(k+1):
plotBt(B,num, i,k)
fig.append('B_{%d,%d}'%(i,k))
else:
plotBt(B,num, i,k)
fig.append('B_{%d,%d}'%(i,k))
plt.legend(fig)
plt.show()
if __name__ == '__main__':
main1()
main2()
6.2 结果
- b i , k ( t ) b_{i,k}(t) bi,k(t) 随 t 的变化
我们需要的是 B i , 4 = b i , 4 B_{i,4} = b_{i,4} Bi,4=bi,4 ,如下
且最终结果
P
(
t
)
=
B
0
,
4
(
t
)
P
0
+
B
1
,
4
(
t
)
P
1
+
B
2
,
4
(
t
)
P
2
+
B
3
,
4
(
t
)
P
3
+
B
4
,
4
(
t
)
P
4
P(t) = B_{0,4}(t)\ P_0 + B_{1,4}(t)\ P_1 + B_{2,4}(t)\ P_2 +B_{3,4}(t)\ P_3 +B_{4,4}(t)\ P_4
P(t)=B0,4(t) P0+B1,4(t) P1+B2,4(t) P2+B3,4(t) P3+B4,4(t) P4
观察曲线可以得知(其实是分析数据),当
t
∈
[
t
4
,
t
5
]
t\in[t_4,t_5]
t∈[t4,t5] 即
t
∈
[
4
/
9
,
5
/
5
]
t\in[4/9,5/5]
t∈[4/9,5/5] 时,
B
i
,
4
B_{i,4}
Bi,4 都为非零值,也就是说该区间的点是由基函数完全定义的。
该区间可以一般性的写为 [ t k , t m − k ] [t_k,t_{m-k}] [tk,tm−k] ,该域也称作某些特定生成曲线的定义域。
- 另外一些实际例子的结果
- E1(均匀)
c
o
n
t
r
o
l
P
o
n
t
s
=
[
[
50
,
50
]
,
[
100
,
300
]
,
[
300
,
100
]
,
[
380
,
200
]
,
[
400
,
600
]
]
k
n
o
t
s
=
{
0
,
1
9
,
2
9
,
3
9
,
4
9
,
5
9
,
6
9
,
7
9
,
8
9
,
1
}
n
=
4
,
m
=
9
,
k
=
4
controlPonts = [[50,50], [100,300], [300,100], [380,200], [400,600]]\\ knots = \{0,\frac{1}{9},\frac{2}{9},\frac{3}{9},\frac{4}{9},\frac{5}{9},\frac{6}{9},\frac{7}{9},\frac{8}{9},1\} \\ n=4,m=9,k=4
controlPonts=[[50,50],[100,300],[300,100],[380,200],[400,600]]knots={0,91,92,93,94,95,96,97,98,1}n=4,m=9,k=4
红色点为控制点,蓝色为要生成的目标曲线。
- E2(非均匀)
c o n t r o l P o n t s = [ [ 50 , 50 ] , [ 100 , 300 ] , [ 300 , 100 ] , [ 380 , 200 ] , [ 400 , 600 ] ] k n o t s = { 0 , 1 9 , 1.5 9 , 5 9 , 5.5 9 , 8 9 , 1 } n = 4 , m = 6 , k = 1 controlPonts = [[50,50], [100,300], [300,100], [380,200], [400,600]]\\ knots = \{0,\frac{1}{9},\frac{1.5}{9},\frac{5}{9},\frac{5.5}{9},\frac{8}{9},1\} \\ n=4,m=6,k=1 controlPonts=[[50,50],[100,300],[300,100],[380,200],[400,600]]knots={0,91,91.5,95,95.5,98,1}n=4,m=6,k=1
- E3(重复点)
c o n t r o l P o n t s = [ [ 50 , 50 ] , [ 100 , 300 ] , [ 300 , 100 ] , [ 380 , 200 ] , [ 400 , 600 ] ] k n o t s = { 0 , 0 , 2 9 , 4 9 , 5 9 , 1 , 1 , 1 } n = 4 , m = 7 , k = 2 controlPonts = [[50,50], [100,300], [300,100], [380,200], [400,600]]\\ knots = \{0,0,\frac{2}{9},\frac{4}{9},\frac{5}{9},1,1,1\} \\ n=4,m=7,k=2 controlPonts=[[50,50],[100,300],[300,100],[380,200],[400,600]]knots={0,0,92,94,95,1,1,1}n=4,m=7,k=2
- E4(重复点)
c
o
n
t
r
o
l
P
o
n
t
s
=
[
[
50
,
50
]
,
[
100
,
300
]
,
[
300
,
100
]
,
[
380
,
200
]
,
[
400
,
600
]
]
k
n
o
t
s
=
{
0
,
0
,
0
,
4
9
,
5
9
,
1
,
1
,
1
}
n
=
4
,
m
=
7
,
k
=
2
controlPonts = [[50,50], [100,300], [300,100], [380,200], [400,600]]\\ knots = \{0,0,0,\frac{4}{9},\frac{5}{9},1,1,1\} \\ n=4,m=7,k=2
controlPonts=[[50,50],[100,300],[300,100],[380,200],[400,600]]knots={0,0,0,94,95,1,1,1}n=4,m=7,k=2
6.3 一些观察结果
- 节点列表的选择相当重要,会影响到最终曲线效果。一般情况下,全均匀分布的貌似效果不好?或者应该在此情况下考虑其定义域,定义域 [ t k , t m − k ] [t_k,t_{m-k}] [tk,tm−k] 内的曲线效果良好。
- 准均匀的样条曲线看起来很好。一个特别之处是,当在两个端点(即 t=0和1)时,如果节点列表分别重复有 k+1 次,那么该端点就类似于贝塞尔曲线端点。即,曲线会经过两端点,且与端点和相邻点连线相切。
- 应该有一些比较好用基于经验化的节点列表,或是好用的幂次样条曲线。
详细样条曲线插值原理见下篇: 样条曲线(下)之插值问题