样条曲线（Spline Curves）

前言： 关于样条曲线的一些总结和使用，其中包含贝塞尔曲线的介绍。内容实际是为了解决在给定控制点的条件下，如何确定一个光滑曲线的问题。样条曲线被广泛应用于模型的几何重构平滑中。
下启：样条曲线(下)之插值问题

文章目录

样条曲线（Spline Curves）
- 0. 基本概念
- 1. 起源
- 2. 贝塞尔曲线（Bézier curve）
- - 1.1 直线段表示（2 points）
  - 1.2 抛物线三切线定理（3 points）
  - 1.3 通用公式
  - 1.4 代码实现
- 3. 曲线连续性描述
- 4. B样条基函数（B-Spline Basis Function）
- - 4.1 多段贝塞尔曲线
  - 4.2 从贝塞尔曲线到B样条基函数
  - - 4.2.1 贝塞尔曲线
    - 4.2.2 B样条
    - 4.2.3 一般数学描述
- 5. 分类
- 6. 代码实现及验证的一些重要问题
- - 6.1 完整python代码
  - 6.2 结果
  - 6.3 一些观察结果

0. 基本概念

样条曲线(Spline Curves)：是给定一系列控制点而得到的一条曲线，曲线形状由这些点控制。一般分为插值样条和拟合样条。

插值：在原有数据点上进行填充生成曲线，曲线必经过原有数据点。

拟合：依据原有数据点，通过参数调整设置，使得生成曲线与原有点差距最小（最小二乘），因此曲线未必会经过原有数据点。

1. 起源

样条曲线起源于一个常见问题，即已知若干点的条件下，如何得到通过这些点的一条光滑曲线？

一个简单且行之有效的方法是，把这些点作为限制点，然后在这些限制点中放置一条具有弹性的金属片，最后金属片绕过这些点后的最终状态即为所需曲线。而最终得到的形状曲线，就是样条曲线。这也是该名字的由来，其中金属片就是样条，形成的曲线就是样条曲线。如下图

该想法虽然巧妙，但显然不具有推广性。因此问题就出来了，如何将其抽象出一个数学模型，从而在已知控制点条件下，仅仅通过数学公式从而获得平滑的样条。

再简化一些，问题便是，如何在两点间进行插值，整体上获得一个平滑曲线。

贝塞尔曲线实际上就是一个样条曲线。可以通过多个控制点，插值获得一个平滑曲线。何为贝塞尔曲线？

2. 贝塞尔曲线（Bézier curve）

1.1 直线段表示（2 points）

从最简单的问题入手，假设有两点P0,P1P_0,P_1P0,P1（端点），如何根据两点确定一条(曲)线？

最简单的问题也有最简单的解答，确定的这条线就是连接P0,P1P_0,P_1P0,P1 两点的线段。从计算机的角度而言，本没有连续的线，因为计算机数据本身就是离散的，因此计算机中线的描述也只是由无数个点组成。

那么，如何用计算机绘制出这条线呢，确切说，如何表达这条线段上任何位置的坐标呢？（其实本质上就是在两点间插值）

解答同样简单，线段上任意位置点坐标可以用一个比例参数 t 来表示，表达式为
Px=P0+(P1−P0)t=(1−t)P0+tP1,t∈[0,1]here,P0PxP0P1=tP_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) 为 P0P_0P0 权重，t 为 P1P_1P1 权重，权重越大，越靠近对应的点。

1.2 抛物线三切线定理（3 points）

问题来到三个点，即如果知道三个点 P0,P1,P2P_0,P_1,P_2P0,P1,P2，且 P0,P2P_0,P_2P0,P2 为端点，那么如何确定一条曲线呢？（理解贝塞尔曲线时，可以把P1P_1P1看作是控制点，我们要的是在控制点下，对P0,P2P_0,P_2P0,P2 两点间的插值）

最简单的想法是，直接连接 P0,P1P_0,P_1P0,P1 和 P1,P2P_1,P_2P1,P2，得到两个线段就可以确定一个曲线，但这明显不是我们想要的，它不够平滑，往往一阶导还不连续（即不满足c1连续，下有详细叙述）。

第二个想法，可以考虑借用二次曲线，但如何在知道三点的情况下确定一个二次曲线呢？

最简单做法是列三个方程组，解方程，然后得到一个二次函数，也即确定了二次曲线，且通过已知的三点。同样，递推到更多已知点，也可以用相似的方法。对于 n 个点，可以使用 n-1 阶次的函数来确定一个唯一的曲线。这种方法也即是常用的多项式插值。前面，提到插值的特点在于，所有已知点都在最终确定的曲线上。

如何使用前述把P1P_1P1当作控制点的思路来解决这个问题呢？即没必要所有点都处在最终确定的曲线上，是否可以参考上面两点情况下，使用权重或是比例参数 t 来确定一条曲线呢？

现在的目的就是，使用三点确定任意一点PxP_xPx ，且最终由无数个PxP_xPx 点形成的曲线连续可导（c1连续）。依据上面两点坐标下确定曲线任意点的思路，异想天开一下。先将三个点划分为两组，即P0,P1P_0,P_1P0,P1和P1,P2P_1,P_2P1,P2。对于两组点，分别使用上面结论，那么在参数 t 下，两组分别确定了两个任意点，记为 Pi,PjP_i,P_jPi,Pj ，也即
Pi=(1−t)P0+P1Pj=(1−t)P1+p2P_i = (1-t)P_0 + P_1 \\ P_j = (1-t)P_1 + p_2 Pi=(1−t)P0+P1Pj=(1−t)P1+p2
可参考下图中的 Pi,PjP_i,P_jPi,Pj 。

但我们想要的是由 P0,P1,P2P_0,P_1,P_2P0,P1,P2 三个点确定的唯一点 PxP_xPx，而现在我们得到了两个点，那再使用一次上面公式不就可以了。也即
Px=(1−t)Pi+tPjP_x = (1-t)P_i + tP_j Px=(1−t)Pi+tPj
问题解决了。

幸运的是，对于t∈[0,1]t\in [0,1]t∈[0,1] 最终确定的线是可导、平滑的。想法有点异想天开，但还是要总结一下，在基于此想法下，用到了三次上面两点下任意点的求法，且有
P0PiP0P1=P1PjP1P2=PjPxPiPj=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
最终结果
Px=(1−t)Pi+tPj=(1−t)[(1−t)P0+tP1]+t[(1−t)P1+tP2]=(1−t)2P0+2t(1−t)P1+t2P2P_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

实际上，上述描述的就是抛物线的三切线定理（性质）。即对于一个抛物线，经过线上两点 P0,P2P_0,P_2P0,P2 的两个切线相交于 P1P_1P1 点，经过在 P0,P2P_0,P_2P0,P2 间抛物线一点 PxP_xPx 的切线，与上述确定的两条切线交于 Pi,PjP_i,P_jPi,Pj 。则有
P0PiP0Pi=P1PjP1P2=PjPxPiPj\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 ，三个方程，三个未知点。最终可以同样反推到上述结果
Px=(1−t)2P0+2t(1−t)P1+t2P2P_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+1n+1n+1 个点 P0,P1,⋯,PnP_0,P_1,\cdots,P_nP0,P1,⋯,Pn ，这 n+1n+1n+1 个点确定了 n 次的贝塞尔曲线。

n 阶贝塞尔曲线 Bn(t)B^{n}(t)Bn(t) 可以由前 n 个点决定的 n-1 次贝塞尔曲线 Bn−1(t∣P0,⋯,Pn−1)B^{n-1}(t|P_0,\cdots,P_{n-1})Bn−1(t∣P0,⋯,Pn−1) 与后 n 个点决定的 n-1 次贝塞尔曲线 Bn−1(t∣P1,⋯,Pn)B^{n-1}(t|P_1,\cdots,P_n)Bn−1(t∣P1,⋯,Pn) 线性组合递推而来，即
Bn(t∣P0,P1,⋯,Pn)=（1−t）Bn−1(t∣P0,P1,⋯,Pn−1)+tBn−1(t∣P1,P2,⋯,Pn)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)
且一次贝塞尔曲线，即为由两点决定的线段，也即
B1(t∣P0,P1)=(1−t)P0+tP1B^1(t|P_0,P_1) = (1-t) P_0 + tP_1 B1(t∣P0,P1)=(1−t)P0+tP1
可以得到 n 次贝塞尔曲线的表达通式为
B(t)=∑i=0nCni(1−t)n−itiPiB(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 combdef inputPoints():controlPoints = []num = 1while 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 + 1if inp == 'n':breakelse:continue#print(controlPoints)return controlPointsdef getInterpolationPoints(controlPoints, tList):n = len(controlPoints)-1interPoints = []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 interPointsif __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]，

c0c^0c0连续：即函数连续，不存在间断点。实例中f(b)=g(b)f(b) = g(b)f(b)=g(b)

c1c^1c1连续：函数的一阶可导且导数连续。实例中 f′(b)=g′(b)f'(b) = g'(b)f′(b)=g′(b)

c2c^2c2连续：函数的二阶可导且导数连续。实例中 f′′(b)=g′′(b)f''(b) = g''(b)f′′(b)=g′′(b)

cnc^ncn连续：函数的nnn阶可导且导数连续。实例中 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 多段贝塞尔曲线

在剖物线三切线定理中，最初想法是三点直接相连，得到曲线虽然 c0c^0c0 连续，但不满足 c1c^1c1 连续。为了解决这个问题，进而使用了二次的贝塞尔曲线。

问题是，对于一个复杂弯曲的曲线，我们希望使用一个贝塞尔曲线来插值获得目标曲线，为此我们可以通过增加控制点来进行插值。但目标曲线越复杂，需要的控制点就越多，而贝塞尔曲线幂次 = 控制点个数 - 1，即需要的计算也越复杂。该方法虽然可行，但是不明智的，低效率的。另外对于单一的贝塞尔曲线，改变其中一个控制点，那么整条曲线都会随之改变。

因此对于复杂曲线，一般做法是，使用三次贝塞尔曲线（常用次）一段一段地拼接成目标曲线，正如 Ps 或 Ai 中使用钢笔工具画出物体轮廓所做的那样。如果使用这种方法，确保最终整体曲线 c1c^1c1 连续的条件是在连接点出两侧的斜率相等，即连接点和其两侧控制点共线。

下图是由两个三次贝塞尔曲线组成的曲线

假设从左到右点依次为 P0,P1,⋅,P7P_0,P_1,\cdot,P_7P0,P1,⋅,P7 ，确保两段曲线拼接起来 c1c^1c1 连续的条件就是 P2,P3,P4P_2,P_3,P_4P2,P3,P4 三点共线。

4.2 从贝塞尔曲线到B样条基函数

4.2.1 贝塞尔曲线

回到最初问题上，通过一系列点，获取一条光滑的曲线。也即通过这些控制点，生成一系列点坐标，这些点坐标形成光滑曲线。

对于贝塞尔曲线上点的生成，是通过如下方程函数
B(t)=∑i=0nCni(1−t)n−itiPi,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)=Wt,n0P0+Wt,n1P1+⋯+Wt,nnPnB(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
其中 Wt,n0W_{t,n}^0Wt,n0 为P0P_0P0 前系数，是最高幂次为 n 的关于 t 的多项式。当 t 确定后，该值就为定值。

因此整个式子可以理解为 B(t) 插值点是这 n+1 个点施加各自的权重 W 后累加得到的。这也是为什么改变其中一个控制点，整个贝塞尔曲线都会受到影响。

其实对于样条曲线的生成，思路就是这样的，最终曲线的生成，就对于各个控制点施加权重即可。在贝塞尔曲线中，该权重是关于 t 的函数即
Wi=Cni(1−t)n−itiW^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 个控制点P0,P1,P2,⋯,PnP_0,P_1,P_2,\cdots,P_{n}P0,P1,P2,⋯,Pn
节点：这个跟控制点没有关系，而是人为地将目标曲线分为若干个部分，其目的就是尽量使得各个部分有所影响但也有一定独立性，这也是为什么B样条中，有时一个控制点的改变，不会很大影响到整条曲线，而只影响到局部的原因，这是区别于贝塞尔曲线的一点。节点划分影响到权重计算，可以预想到的是，实现局部性的影响的原理应该是在生成某区间内的点时，某些控制点前的权重值会为0，即对该点没有贡献，所以才有上述特点。事实上，就是如此的。先了解有这个概念即可。假设我们划分了 m+1 个节点 t0,t1,⋯,tmt_0,t_1,\cdots,t_mt0,t1,⋯,tm，将曲线分成了 m 段
次与阶 ：次的概念就是贝塞尔中次的概念，即权重中 t 的最高幂次。而阶=次+1。这里只需要知道次这个概念即可。假设我们用 k 表示次，即 k 次。

下面用一个实例来说明

如上图，对应为

控制点有 n+1 = 5 个，n=4，即 P0,P1,P2,P3,P4P_0,P_1,P_2,P_3,P_4P0,P1,P2,P3,P4

节点规定为 m+1 = 10 个，m=9，即 t0,t1,⋯,t9t_0,t_1,\cdots,t_9t0,t1,⋯,t9 ，该节点将要生成的目标曲线分为了 9 份，这里的 t 取值一般为 0-1的一系列非递减数。如 0 代表起始位置，1 代表末位置，0.5 代表一半的位置。t0,t1,⋯,t9t_0,t_1,\cdots,t_9t0,t1,⋯,t9 组成的序列，叫做节点表。如等分的节点表 {0,19,29,39,49,59,69,79,89,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+1m=n+k+1m=n+k+1 ，即 k=m−n−1=4k = m-n-1 = 4k=m−n−1=4 ，先对这个关系式有个印象即可。

下面开始用上面例子说明 B样条是如何工作的。

我们的目的就是获取最终样条曲线，而根据前述，也即是获得任意点的坐标值，即 t∈[0,1]t\in[0,1]t∈[0,1] 时，获得一系列的点，这些点最终组成了目标曲线。这些点的坐标可以表示为各控制点 PiP_iPi 与其权重 WiW_iWi 的乘积再累加，即
B(t)=∑i=0nWiPiB(t) = \sum\limits_{i=0}^n W_iP_i B(t)=i=0∑nWiPi
也即我们获得 WiW_iWi 即可。WiW_iWi 是关于 t 的函数，最高幂次为 k。B样条中通常记为 Bi,k(t)B_{i,k}(t)Bi,k(t) ，即表示第 i 点的权重，关于 t 的函数，且最高幂次为 k。而这个权重函数 Bi,k(t)B_{i,k}(t)Bi,k(t) ，在B样条里叫做 k次B样条基函数。

下面就是关键，权重函数 Bi,k(t)B_{i,k}(t)Bi,k(t) 是如何取值的？这里就要用到前面的节点表的值了。针对上面实例，我们现在的目标是获取 5 个控制点P0,P1,P2,P3,P4P_0,P_1,P_2,P_3,P_4P0,P1,P2,P3,P4 对应的 5 个权重值 B0,4,B1,4,B2,4,B3,4,B4,4B_{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 ），大可当作不知道。

表中 bi,kb_{i,k}bi,k 代表含义跟 Bi,kB_{i,k}Bi,k 是一样的，因为不是我们需要的最终值，而是需要求解的中间递推值，所以用小写表示。对于B样条权重，其是根据节点来计算，理论上，m+1个节点，可以得到如表所示的，每个k对应 9(m=9) 个b值，但我们最终只需要5个值，所以就需要对应规则限制。B样条获取5个权重值过程和规则如下（B样条基函数）

k = 0 时，即 0 次时，如果 t∈[tj,tj+1]t\in [t_j,t_{j+1}]t∈[tj,tj+1] ，那么规定 bj,0=1b_{j,0} = 1bj,0=1，其余都为 0。如 t∈[t2,t3]t\in[t_2,t_3]t∈[t2,t3] 就有如下结果

k = 1 时，bj,1b_{j,1}bj,1 取值与低一次的相邻两节点相关，即与 0 次同域两端的节点相关。即 bj,1b_{j,1}bj,1 求解与bj,0,bj+1,0b_{j,0},b_{j+1,0}bj,0,bj+1,0 相关，如 b1,1b_{1,1}b1,1 就是依据 b1,0,b2,0b_{1,0},b_{2,0}b1,0,b2,0 求得。求解关系形式为
bj,1=A(t)bj,0+B(t)bj+1,0b_{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−tjtj+k−tjB(t)=tj+k+1−ttj+k+1−tj+1A(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 的一次函数即可。这里有点类似与贝塞尔函数中递推过程中的前两值的线性组合。

我们用箭头表示上述计算过程，如下

值得说明的两点：

在bj,1b_{j,1}bj,1 中，有些值会取 0，因为在 0 次的 b 值有 0 。观察我们可以很容易的值，此时值为 0 的有 b0,1,b3,1,b4,1,b5,1,b6,1,b7,1b_{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 取值在 [t2,t3][t_2,t_3][t2,t3] 范围的条件下)。而取值为非0的 b1,1,b2,1b_{1,1},b_{2,1}b1,1,b2,1 都是关于 t 的一次函数（因为A(t),B(t)A(t),B(t)A(t),B(t) 为 t 的一次函数）。前面提到想要使用局部控制曲线（而非全局）需要某些值为 0，所以是符合预期的。
b8,1b_{8,1}b8,1 是无法求解的，因为它需要 b9,0b_{9,0}b9,0 ，而该值不存在，所以，在 k=1 时，能得到的 bj,1b_{j,1}bj,1 减少了 1 （k=1）个，共有8个。同样是令人振奋的结果，因为我们的目标是需要 5 个值（因为有5个控制点），多减少几个不就达到目标了嘛。

k = 2 时，同样的过程，不再以表列出，同样有些 bj,2b_{j,2}bj,2 会取值为 0。有些值为非 0, 不同的一点在于，此时系数又乘了一次 A,B 函数，因此那些非 0 值，此时是关于 t 的二次幂函数。另外得到的 bbb 减少了 **2（k=2）**个，还剩 m−k=9−2=7m-k = 9-2 = 7m−k=9−2=7 个。

k = x时，同样，非0值会是关于 t 的 x 次函数，这也是为什么 k 会是 t 的最高次幂。b 会减少 k 个，还剩 m−km-km−k 。

我们最终那个需要 5 个（即 n+1 个）值，那么就需要有 m−k=n+1m-k = n + 1m−k=n+1 ，也即 m=k+n+1m = k + n + 1m=k+n+1 ，这便是上面说到的 m,n,km,n,km,n,k 关系式的由来。所以实例中 k=m−n−1=9−4−1=4k = m-n-1 = 9-4-1 = 4k=m−n−1=9−4−1=4 ，即当 k=4 时，恰好还剩 5 个 b 值，就是我们需要的B样条基函数B。此时，最高幂次为 4。最终图表示意如下：

该图表表示的是在 t 取值在 [t2,t3][t_2,t_3][t2,t3] 范围时，求得的 5 个控制点对应的权重Bj,4B_{j,4}Bj,4（基函数），最终结果B0,4,B1,4,B2,4B_{0,4},B_{1,4},B_{2,4}B0,4,B1,4,B2,4 都为非0值，且是关于 t 的四次函数，B3,4,B4,4B_{3,4},B_{4,4}B3,4,B4,4 为 0。也就是说，本例中，[t2,t3][t_2,t_3][t2,t3] 区间的曲线上点不受第四和第五个控制点变化的影响。

这就是所有的 B样条曲线生成的工作原理。

值得注意的几点性质：

对于 t∈[tj,tj+1]t\in[t_j,t_{j+1}]t∈[tj,tj+1] 时，最终非零 Bj,kB_{j,k}Bj,k为 Bj−k,k,⋯,Bj,kB_{j-k,k},\cdots,B_{j,k}Bj−k,k,⋯,Bj,k （如果j−k>=0j-k>=0j−k>=0的话），即第 j−kj-kj−k 至第 jjj 个控制点控制影响着 [tj,tj+1][t_j,t_{j+1}][tj,tj+1] 区间内曲线点。可进行局部控制。
最终所得曲线结果一般不会经过控制点，因为最终得到的非0权重始终不止一个
tit_iti 为非递减的，可以取重值
对于图中，

4.2.3 一般数学描述

将上面过程作为严谨的数学描述

对于 n+1 个控制点 P0,P1,⋯,PnP_0,P_1,\cdots,P_nP0,P1,⋯,Pn ，有一个包含 m+1 个节点的列表（或节点向量）t0,t1,⋯,tm{t_0,t_1,\cdots,t_{m}}t0,t1,⋯,tm ，其 k 次B样条曲线表达式为（且m=n+k+1）
P(t)=∑i=0nBi,k(t)PiP(t) = \sum\limits_{i=0}^n B_{i,k}(t)\ P_i P(t)=i=0∑nBi,k(t) Pi
其中 Bi,k(t)B_{i,k}(t)Bi,k(t) 称为 k 次 B 样条基函数，也叫调和函数。且 Bi,k(t)B_{i,k}(t)Bi,k(t) 满足如下递推式（de Boor递推式）
k=0,Bi,0(t)={1,t∈[ti,ti+1]0,Otherwisek>0,Bi,k(t)=t−titi+k−tiBi,k−1(t)+ti+k+1−tti+k+1−ti+1Bi+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,km = knots.shape[0]-1n = controlPoints.shape[0]-1k = m - n - 1# initialize B by zeros B = np.zeros((k+1, m))# get t regiontStart = 0for x in range(m+1):if t==1:tStart = m-1if knots[x] > t:tStart = x-1break# calculate B(t)for _k in range(k+1):if _k == 0:B[_k, tStart] = 1else:for i in range(m-_k):if knots[i+_k]-knots[i]== 0:w1 = 0else:w1 = (t-knots[i])/(knots[i+_k]-knots[i]) if knots[i+_k+1]-knots[i+1] == 0:w2 = 0else: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) * controlPointspt = 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]-1n = controlPoints.shape[0]-1k = m - n - 1print('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]):continueBt = 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]-1n = controlPoints.shape[0]-1k = m - n - 1print('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 结果

bi,k(t)b_{i,k}(t)bi,k(t) 随 t 的变化

我们需要的是 Bi,4=bi,4B_{i,4} = b_{i,4}Bi,4=bi,4 ,如下

且最终结果
P(t)=B0,4(t)P0+B1,4(t)P1+B2,4(t)P2+B3,4(t)P3+B4,4(t)P4P(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∈[t4,t5]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] 时，Bi,4B_{i,4}Bi,4 都为非零值，也就是说该区间的点是由基函数完全定义的。

该区间可以一般性的写为 [tk,tm−k][t_k,t_{m-k}][tk,tm−k] ，该域也称作某些特定生成曲线的定义域。

另外一些实际例子的结果
- E1（均匀）

controlPonts=[[50,50],[100,300],[300,100],[380,200],[400,600]]knots={0,19,29,39,49,59,69,79,89,1}n=4,m=9,k=4controlPonts = [[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（非均匀）

controlPonts=[[50,50],[100,300],[300,100],[380,200],[400,600]]knots={0,19,1.59,59,5.59,89,1}n=4,m=6,k=1controlPonts = [[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（重复点）

controlPonts=[[50,50],[100,300],[300,100],[380,200],[400,600]]knots={0,0,29,49,59,1,1,1}n=4,m=7,k=2controlPonts = [[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（重复点）

controlPonts=[[50,50],[100,300],[300,100],[380,200],[400,600]]knots={0,0,0,49,59,1,1,1}n=4,m=7,k=2controlPonts = [[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 一些观察结果

节点列表的选择相当重要，会影响到最终曲线效果。一般情况下，全均匀分布的貌似效果不好？或者应该在此情况下考虑其定义域，定义域[tk,tm−k][t_k,t_{m-k}][tk,tm−k] 内的曲线效果良好。
准均匀的样条曲线看起来很好。一个特别之处是，当在两个端点(即 t=0和1)时，如果节点列表分别重复有 k+1 次，那么该端点就类似于贝塞尔曲线端点。即，曲线会经过两端点，且与端点和相邻点连线相切。
应该有一些比较好用基于经验化的节点列表，或是好用的幂次样条曲线。

详细样条曲线插值原理见下篇: 样条曲线(下)之插值问题

详解样条曲线（上）（包含贝塞尔曲线）相关推荐

一文弄懂元学习 (Meta Learing)（附代码实战）《繁凡的深度学习笔记》第 15 章元学习详解（上）万字中文综述
<繁凡的深度学习笔记>第 15 章元学习详解 (上)万字中文综述(DL笔记整理系列) 3043331995@qq.com https://fanfansann.blog.csdn.net ...
Drools 规则语言详解（上）
http://www.blogjava.net/guangnian0412/archive/2006/06/09/51574.html http://www.blogjava.net/guangnia ...
JavaScript数据结构与算法——链表详解（上）
注:与之前JavaScript数据结构与算法系列博客不同的是,从这篇开始,此系列博客采用es6语法编写,这样在学数据结构的同时还能对ECMAScript6有进一步的认识,如需先了解es6语法请浏览ht ...
JavaScript数据结构与算法——列表详解（上）
列表是一组有序的数据,每个数组中的数据项称为元素.数组相关知识不够了解的伙伴可以阅读本人上篇博客在JavaScript中,列表的元素可以是任意数据类型.列表中可以保存不定数量的元素,实际使用时元素的数 ...
中台详解（上）——什么是中台
中台详解(上)--什么是中台人人都是产品经理编辑导语:中台这一概念,这两年在国内大热,不少企业接连开始组织架构的调整,意图建设中台:中台到底是什么?哪类企业可以搭建中台?本文作者对中台进行了非常详 ...
php读取图片文件流,详解php文件包含原理（读取文件源码、图片马、各种协议、远程getshell等）...
详解php文件包含原理(读取文件源码.图片马.各种协议.远程getshell等) 作者是namezz (看完图相当于做了一轮实验系列) 现有文件代码如下 1.png (21.16 KB, 下载次数: ...
include详解 shell_详解php文件包含原理（读取文件源码、图片马、各种协议、远程getshell等） ......
详解php文件包含原理(读取文件源码.图片马.各种协议.远程getshell等) 作者是namezz (看完图相当于做了一轮实验系列) 现有文件代码如下 include和include_once.re ...
Python 装饰器详解（上）
Python 装饰器详解(上) 转自:https://blog.csdn.net/qq_27825451/article/details/84396970,博主仅对其中 demo 实现中不适合pyth ...
format函数_Python学习教程：Python3之字符串格式化format函数详解（上）
Python学习教程:Python3之字符串格式化format函数详解(上) 概述在Python3中,字符串格式化操作通过format()方法或者f'string'实现.而相比于老版的字符串格式化方 ...
我的世界刷猪人塔java版_我的世界速攻猪人塔详解史上最牛的经验塔
我的世界速攻猪人塔详解史上最牛的经验塔.那下面给大家分享的这个是一个可以让所有经验塔自叹不如的速攻猪人塔哦~那到底这个塔是什么呢?那下面就给大家详细的介绍一下吧!有感兴趣的玩家不妨进来看看哦~希望大 ...

详解样条曲线（上）（包含贝塞尔曲线）