Spline Curve.

样条（Spline）是工人用来在已知若干控制点的情况下画光滑曲线的工具，它是一根柔软但有弹性的金属片；使用时将金属片绕过控制点便可以构造一条光滑的曲线。

给定一系列控制点，通过这些点构造的光滑曲线称为样条曲线(Spline Curve)。样条曲线一般分为插值样条和拟合样条。

插值：在原有数据点上进行填充生成曲线，曲线必经过原有数据点。
拟合：依据原有数据点，通过参数调整设置，使得生成曲线与原有点差距最小，曲线未必会经过原有数据点。

1. 贝塞尔曲线 (Bézier Curve)

对于空间中的任意两点$P_0,P_1$，光滑曲线构造为连接两点的直线：

\[P_x= P_0+(P_1-P_0)t=(1-t)P_0+tP_1, t\in[0,1]\]

对于空间中的任意三点$P_0,P_1,P_2$，先分别构造连接两组点$P_0,P_1$、$P_1,P_2$的直线，再对两个直线上的点$P_i,P_j$进行插值：

\[\begin{aligned} P_i&= (1-t)P_0+tP_1, P_j= (1-t)P_1+tP_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^2P_2 \end{aligned}\]

上式分别为一阶和二阶贝塞尔曲线(Bézier Curve)，其中阶数是指变量$t$的最高幂次。贝塞尔曲线是递推的：$n$阶贝塞尔曲线（对应$n+1$个控制点）是两个$n-1$阶贝塞尔曲线的线性组合。

假设一共有$n+1$个控制点$P_0,P_1,\cdots,P_n$，这$n+1$个点确定了$n$阶贝塞尔曲线$B^{n}(t)$，可以由前$n$个点决定的$n-1$阶贝塞尔曲线$B^{n-1}(t|P_0,\cdots,P_{n-1})$与后$n$个点决定的$n-1$阶贝塞尔曲线$B^{n-1}(t|P_1,\cdots,P_n)$线性组合递推而来，即:

\[\begin{aligned} &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) \\ &= (1-t)^2B^{n-2}(t|P_0,P_1,\cdots,P_{n-2}) + 2t(1-t)B^{n-2}(t|P_1,P_2,\cdots,P_{n-1}) + t^2B^{n-2}(t|P_2,P_3,\cdots,P_n) \\ &= \sum\limits_{i=0}^n C_n^i(1-t)^{n-i}t^i B^0(t| P_i)= \sum\limits_{i=0}^n C_n^i(1-t)^{n-i}t^i P_i \end{aligned}\]

其中$C_n^i$是组合数，$B^0$是零阶曲线（即控制点本身）。使用Python绘制贝塞尔曲线：

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import comb

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 = [[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()

2. B样条曲线 (B-Spline Curve)

贝塞尔曲线幂次等于控制点个数减一，目标曲线越复杂，需要的控制点就越多，计算也越复杂。并且贝塞尔曲线改变其中一个控制点，整条曲线都会随之改变。

观察贝塞尔曲线的表达式，本质是不同控制点的加权累加。对于$n+1$个控制点$P_0,P_1,\cdots,P_n$，B样条曲线 (B-Spline Curve)额外引入$m+1$个节点$t_0,t_1,\cdots,t_m$将曲线分成$m$段，使得各段有所影响但也有一定独立性。在生成某区间内的点时，某些控制点权重置零，从而实现曲线控制的局部性。B样条曲线的表达式：

\[B(t) = \sum\limits_{i=0}^n B_{i,k}P_i\]

其中权重函数$B_{i,k}$是关于$t$的$k$阶函数，表示第$i$点的权重，称作$k$次B样条基函数。$k$次B样条基函数的计算可以通过列出节点表实现：

$k = 0$时 ($0$阶)，如果$t\in [t_j,t_{j+1}]$，则规定$b_{j,0} = 1$，其余都为$0$。

\[B_{i,0}(t) = \begin{cases} 1, & t∈[t_i,t_{i+1}] \\ 0, & \text{Otherwise} \end{cases}\]

$k = 1$时，$b_{j,1}$取值与低一次的相邻两节点$b_{j,0},b_{j+1,0}$相关，求解关系形式为关于 $t$ 的一次幂函数：

\[\begin{aligned} b_{j,1} &= A(t)\ b_{j,0} + B(t)\ b_{j+1,0} \\ 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}} \end{aligned}\]

重复上述过程，构造$k$次B样条基函数的递推公式（de Boor递推式），形式为关于 $t$ 的$k$次幂函数：

\[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)\]

下图给出了$t$取值在$[t_2,t_3]$范围的条件下的计算过程。

对于$t\in[t_j,t_{j+1}]$时，最终非零$B_{j,k}$为$B_{j-k,k},\cdots,B_{j,k}$，即第$j−k$至第$j$个控制点控制影响着$[t_j,t_{j+1}]$区间内曲线点。故B样条曲线可进行局部控制。并且最终得到的非0权重始终不止一个，因此最终所得曲线结果一般不会经过控制点。下图对$10$个节点下的$1$至$4$阶B样条基函数进行可视化：

注意到阶次每提高$1$阶，计算得到的$B_{\cdot,k}$会减少$1$个。因此若$k$次B样条基函数对应$n+1$个控制点，则节点数量$m+1$应满足：

\[m-k = n+1 \quad \to \quad m=n+k+1\]

节点值$t_i$为非递减的，可以取重复值：

均匀 B 样条（uniform B-spline curve）：节点均匀分布
准均匀 B 样条：在开始和结束处的节点可重复，中间节点均匀分布
非均匀 B 样条：节点非均匀分布
受限 B 样条（clamped B-spline curve）：节点列表中首末两端$0,1$重复$k+1$次，曲线会经过首末控制点

使用Python绘制B样条曲线：

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(t)，即权重，乘以控制点坐标，从而求出曲线上点坐标
def getPt(Bt, controlPoints):
	Bt = np.array(Bt)
	ptArray = Bt.reshape(-1,1) * controlPoints
	pt = ptArray.sum(axis = 0)
	return pt

if __name__ == '__main__':
	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
    
	for t in np.linspace(0,1,100):
		if 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()

3. NURBS曲线 (NURBS Curve)

B样条无法描述圆锥曲线，为解决此问题，引入非均匀有理B样条（non-uniform rational b-spline, NURBS）。

非均匀：节点距离不相等
有理：控制点权值不相等

NURBS在B样条的基础上为每个控制点$p_i$加入了一个权重$w_i$：

\[C(t) = \sum\limits_{i=0}^n \left(\frac{B_{i,k}(t)w_i}{\sum_{j=0}^nB_{j,k}(t)w_j}\right)P_i\]