样条曲线(下)之插值问题(贝塞尔曲线、B样条和一般样条曲线插值)

贝塞尔曲线插值与B样条插值

前言： 这篇是“样条曲线”的接续，前面主要集中在了理论部分，这篇文章主要内容是贝塞尔曲线与B样条是如何应用到插值中的。

前篇：样条曲线

0. 插值问题

问题是一切发展的动力。还是从问题出发，考虑怎么通过数学手段解决问题。

对于插值问题，就是已知一些离散点，插值出一些新的点，使得所有点形成一条光滑曲线。

如已知六个点 $P_0,P_1,\cdots ,P_5$ ，如图红色点 $(1,1),(3,6),(6,3),(8,0),(11,6),(12,12)$ ，如何通过插值生成类似于图中蓝色光滑曲线，而非僵硬的绿色多线段。

1. 贝塞尔曲线插值

1.1 曲线的数学描述

依据前篇的数学描述，假设一共 $n+1$ 个控制点 $P_0,P_1,\cdots ,P_n$ ，这 $n+1$ 个点可以确定一条 $n$ 次贝塞尔曲线。可以用如下式计算曲线上任意点坐标
$$B(t)=\sum _{i=0}^n{C}_n^i(1-t{)}^{n-i}{t}^i{P}_i$$
特别的，有几个常用贝塞尔曲线

二次贝塞尔曲线
$$B(t)=(1-t{)}^2{P}_0+2t(1-t)P_1+t^2{P}_2,t\in [0,1]$$
三次贝塞尔曲线
$$B(t)=(1-t{)}^3{P}_0+3t(1-t{)}^2{P}_1+3t^2(1-t)P_2+t^3{P}_3$$

1.2 曲线插值

这里以最常用的三次贝塞尔曲线为例，探究其插值应用。

对于贝塞尔曲线，有三个要点：

• 其是通过控制点来生成的，控制点不全在最终曲线上。
• 控制点首末的两点是最终曲线的端点（意味着首末控制点会在最终曲线上），且各自与相邻点的连线同最终曲线相切。
• 两个贝塞尔曲线如果平滑连接，则需要连接点与其左右相邻两端点共线。

基于要点2，可以在把两个已知点作为三次贝塞尔曲线的两个端点 $P_0,P_3$，然后想办法再指定两个控制点 $P_2,P_3$ 即可。

且控制点设定要满足贝塞尔曲线相连时平滑。

总结一下就是：

• 在两点中生成贝塞尔曲线，且这两个点作为三次贝塞尔曲线的两个端点。
• 再确定额外两个控制点，控制点设定满足两相邻的贝塞尔曲线连接点和其连接点左右控制点共线。

如上图是由两个贝塞尔曲线组成的一条平滑曲线，经过点 $P_1,P_2,P_3$ ，且 $P_1,P_2$ 是第一条三次贝塞尔曲线的两端点，$P_2,P_3$ 是第二条三次贝塞尔曲线的两端点。我们需要的是每个贝塞尔曲线再设定两个控制点即可。

最简单直接的思路就是列方程，解方程。如上，两条贝塞尔曲线，一共四个未知点(控制点)，需要四个方程。

• 令其连接点满足 $c^1$ 连续，也即 $P_1^2,P_2,P_2^1$ 共线。此处有一个方程。

• 对连接点处 $c^2$ 连续也加以考虑。此处也有一个方程。

  还差两个方程。限制两端点即可。比如二阶导为某一定值，或三阶导为某一定值。诸如此类，这一点在后面会有详细描述。B样条中，简单的三次样条中都有详细描述。

这种思路原理直接且简单，不再用代码实现。无非就是列方程、解方程的问题。这一点在B样条中有详细类似的过程，但比这复杂的多。

下面详细描述一种有趣的解决方案：

最重要的一点是满足 $P_1^2,P_2,P_2^1$ 三点共线，所以只要关注$P_1^2,P_2^1$ 的设定即可，因为$P_1^1,P_2^2$ 可以依据同样的规则设定。解决方案如下，$P_1^2,P_2^1$ 肯定不是凭空产生，需要依赖于已有的一些点。而其中可依赖的点有 $P_1,P_2,P_3$ 。

第一种思路可以在$P_1,P_2$ 连线上，控制一个比例因子生成一点，同样 $P_2,P_3$ 连线上亦如是。如上图中的蓝色两点，然后连接两点形成一条线段，平移线段，使其通过 $P_2$ ，平移后的两点就是我们需要的 $P_1^2,P_2^1$ 。

第二种思路可以直接在$P_1,P_3$ 连线上直接根据某特定比例因此生成两点，按照同样的平移步骤，平移后的该两点当作控制点。当然比例因此可以根据先验信息，比如参考线 $P_1{P}_2$ 和线 $P_2,P_3$ 长度比的数值。等等。

1.3 代码实现

import numpy as np
import matplotlib.pyplot as plt

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


def getControlPointList(pointsArray, k1=-1, k2=1):
	points = np.array(pointsArray)
	index = points.shape[0] - 2

	res = []
	for i in range(index):
		tmp = points[i:i+3]
		p1 = tmp[0]
		p2 = tmp[1]
		p3 = tmp[2]

		if k1 == -1:
			l1 = np.sqrt(np.sum((p1-p2)**2))
			l2 = np.sqrt(np.sum((p2-p3)**2))
			k1 = l1/(l1+l2)
			k2 = l2/(l1+l2)

		p01 = k1*p1 + (1-k1)*p2
		p02 = (1-k2)*p2 + k2*p3
		p00 = k2*p01 + (1-k2)*p02
		
		sub = p2 - p00
		p12 = p01 + sub
		p21 = p02 + sub

		res.append(p12)
		res.append(p21)
	pFirst = points[0] + 0.1*(res[0] - points[0])
	pEnd = points[-1] + 0.1*(res[-1] - points[-1])
	res.insert(0,pFirst)
	res.append(pEnd)

	return np.array(res)

if __name__ == '__main__':
	points = [[1,1],[3,6],[6,3],[8,0],[11,6],[12,12]]
	controlP = getControlPointList(points)
	l = len(points) - 1
	figure = plt.figure()
	t = np.linspace(0,1,50)
	for i in range(l):
		p = np.array([points[i], controlP[2*i], controlP[2*i+1], points[i+1]])
		interPoints = getInterpolationPoints(p, t)
		x = np.array(interPoints)[:,0]
		y = np.array(interPoints)[:,1]
		plt.plot(x,y)
	plt.scatter(np.array(points)[:,0],np.array(points)[:,1],color='gray')
	plt.show()

结果

分析

结果如预期一样，不会特别好，但可以确定该种方法是可行的。结果受比例因子 $k1,k2$ 影响很大，受原始点（型值点）的密度、趋势变化影响。

观察也知，该想法对于没有明显趋势反转（凸曲线变为凹曲线）的点效果还行，但对于转折点，效果不太行。分析其原因，是在起初分析时，使用的图是便没有考虑这种情况。如下，

而该问题的本质原因，是因为该思路只限制了连接点处 $c^1$ 连续，而对$c^2$ 连续没有考虑，所以就出现了结果图中连接点确实看起来很平滑，但后面的过渡不太好。这也侧面说明了插值问题中 $c^2$ 连续对于平滑是很重要的。

2. B样条插值

2.1 数学表达和一些补充

具体推导见前篇：样条曲线

对于 n+1 个控制点 $P_0,P_1,\cdots ,P_n$ ，有一个包含 m+1 个节点的列表（或节点向量）$t_0,t_1,\cdots ,t_m$ ，其 k 次B样条曲线表达式为（且m=n+k+1）
$$P(t)=\sum _{i=0}^n{B}_{i,k}(t)P_i$$
其中 $B_{i,k}(t)$ 称为 k 次 B 样条基函数，也叫调和函数。且 $B_{i,k}(t)$ 满足如下递推式（de Boor递推式）
$$\begin{gathered} k=0,B_{i,0}(t)=\{\begin{array}{c}0,t\in [t_i,t_i+1]\\ 1,Otherwise\end{array} \\ 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) \end{gathered}$$

前篇中在尝试各种节点列表时，得到了很多结果，但对结果总结分析较少，以下是一些非常重要的小总结和补充

• 前篇实例中许多曲线看起来很奇怪，是因为没有考虑定义域，也即由基函数完全支持的部分。该定义域在前篇中有提到，是 $[t_k,t_{n+1}]or[t_k,t_{m-k}]$ ，这点十分重要。在考虑定义域的情况下，曲线平滑且不会显示地很奇怪。

  定义域是依赖于节点列表的，值得一提的是，如果想要较广定义域$[t_k,t_{n+1}]$，从式子上来看，至少有两种思路

  • 使得定义域内 t 区间数量尽量多，这种情况在节点均匀分布时可用，别的情况可能没用。即 $k$ 尽量小，$n+1$ 尽量大，也即 $m-(n+1)$ 尽量小，$n+1$ 尽量大。也即 $m$ 与 n+1 尽量接近，n+1 尽量大。换句话说，节点个数尽量多，控制点数与节点数接近。

    使用该方法重新设定节点列表绘制前篇中第一个实例曲线
    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 20 , 2 20 , 8 20 , 14 20 , 19 20 , 1 , 1 } n = 4 , m = 7 , k = 2 d o m a i n = [ t 3 , t 5 ] = [ 2 20 , 19 20 ] controlPonts = [[50,50], [100,300], [300,100], [380,200], [400,600]]\\ knots = \{0,\frac{1}{20},\frac{2}{20},\frac{8}{20},\frac{14}{20},\frac{19}{20},1,1\} \\ n=4,m=7,k=2 \\ domain = [t_3,t_5] = [\frac{2}{20},\frac{19}{20}] controlPonts=[[50,50],[100,300],[300,100],[380,200],[400,600]]knots={0,201​,202​,208​,2014​,2019​,1,1}n=4,m=7,k=2domain=[t3​,t5​]=[202​,2019​]
    定义域已经几乎包括了所有区间，最终定义域内曲线如下

    

    曲线没有包含全部，但 $[\frac{2}{20},\frac{19}{20}]$ 也几乎包含全部了。Excellent！

  • 使得定义域内 t 区间跨度大，其实严格来说就是我们要解决的问题本质，第一种思路只是实现第二种的一种解决方案。跨度大，也即 $t_{n+1}-t_k$ 尽量大。如果 $m,n,k$ 确定的情况下，如何确保呢。可以改变节点列表中的 t 的分布。我们需要的是区间$[t_k,t_{m-k}]$ ，如果前 k+1 个 t 数据全部设为 0，而第 m-k 个 t 之后的数据全部设为 1。那么区间 $[t_k,t_{m-k}]$ 跨度不就是 $[0,1]$ 嘛，也就是整个曲线嘛，很棒的想法。从曲线的角度来看，此时就是我们绘制出的整个曲线。就没有之前绘图中出现的奇奇怪怪的曲线变化了。

    使用该方法重新设定节点列表绘制前篇中第一个实例曲线
    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 , 0 , 1 2 , 1 , 1 , 1 , 1 } n = 4 , m = 8 , k = 3 d o m a i n = [ 0 , 1 ] controlPonts = [[50,50], [100,300], [300,100], [380,200], [400,600]]\\ knots = \{0,0,0,0,\frac{1}{2},1,1,1,1\} \\ n=4,m=8,k=3 \\ domain = [0,1] controlPonts=[[50,50],[100,300],[300,100],[380,200],[400,600]]knots={0,0,0,0,21​,1,1,1,1}n=4,m=8,k=3domain=[0,1]
    

    不再有奇怪的点，生成的曲线也十分"完美"。

  第二种思路极其重要，如果 $[t_k,t_{m-k}]$ 区间内节点又均匀分布，该样条曲线也被称作为准均匀样条。在首末两端设定一定数量的 0,1 的方法在样条曲线中极其常见，因为它确实很实用啊。

• 前篇中也提到了一个有趣且及其重要的现象，即 当节点列表中首末两端 0,1 重复 k+1次时，那么生成的样条曲线首末两端同贝塞尔曲线一样，都是实际原始点（型值点）的首末两点。该样条曲线也有个特定的名字，叫 clamped B-spline curves 。结合上面第二种思路提到的内容，我们得到的样条曲线也就是整条曲线。

• 其中还有一个 Wrapping 问题，即，使得生成曲线形成闭环。待用到时再探究和总结（挖坑）。

2.2 曲线插值

贝塞尔曲线中应用插值时提到了生成曲线端点为数据点（型值点）的重要性，该特点是实现插值的基础，或者说是实现两点间插值的基础。幸运的是，B样条中有一种样条满足这样的条件，即上述的 clamped B-spline。需要节点首末的 $0,1$ 重复$k+1$ 次。

我们常用三次的样条插值，下面就以三次样条插值为例，对其进行详细分析。

2.2.1 三次 clamped B-样条

对于三次 camped B-样条，有以下要点

• 次数 k
  $$k=3$$

• 型值点个数 $n+1$ 与 节点个数 $m+1$ 满足

$$m=n+k+1\Rightarrow m-n=4$$

• camped 情况下，节点 $t_0,t_1,\cdots ,t_m$ 满足
  $$\{\begin{array}{c}\begin{array}{rl}& t_0=t_1=t_2=t_3=0\\ & t_4,t_5,\cdots ,t_n\in [0,1]\\ & t_{n+1}=t_{n+2}=t_{n+3}=t_{n+4}=1\end{array}\end{array}$$
  特别地，如果使用准均匀样条，即 $t_4,t_5,\cdots ,t_n$ 均匀分布（当然也可以使用别的分布方法），即有
  $$\begin{array}{rl}& t_4=\frac{4-3}{n+1-k}=\frac{1}{n-2}\\ & t_5=\frac{5-3}{n-2}=\frac{2}{n-2}\\ & \cdots \\ & t_j=\frac{j-3}{n-2}\\ & \cdots \\ & t_n=\frac{n-3}{n-2}\end{array}$$

为了形象展示，使用之前类似图表来表示

注意，其中橙色标注的区间为 $[0,0]or[1,1]$ 长度为 0，另，又因为其不在定义域内，我们不作计算，也即，其 0 次 基函数值始终记为 0。假设我们计算定义域内某区间 $[t_j,t_{j+1}]$ 的曲线，按照前篇介绍的递推计算公式，0 次基函数 $b_{j,0}=1$ 其余为 0 ，正如上表所示。

我们要得到 n+1 个点的权重，即 $B_{0,3},B_{1,3},\cdots ,B_{n,3}$ ，由上我们知道，对于一段区间（如 $[t_j,t_{j+1}]$），三次基函数只会有四个非零值，分别为 $B_{j,3},B_{j-1,3},B_{j-2,3},B_{j-3,3}$

如前篇所述，递推公式为
$$b_{j,k}(t)=\frac{t-t_j}{t_{j+k}-t_j}{b}_{j,k-1}(t)+\frac{t_{j+k+1}-t}{t_{j+k+1}-t_{j+1}}{b}_{j+1,k-1}(t)$$

四个权重求解具体如下，当 $t\in [j,j+1]$

• 0 次 基函数有一个非零 $b_{j,0}$

  $b_{j,0}=1$

• 1 次 基函数有两个非零 $b_{j,1},b_{j-1,1}$

  $\begin{array}{rl}{b}_{j,1}(t)& =\frac{t-t_j}{t_{j+1}-t_j}{b}_{j,0}(t)+\frac{t_{j+2}-t}{t_{j+2}-t_{j+1}}{b}_{j+1,0}(t)\\ & =\frac{t-t_j}{t_{j+1}-t_j}\end{array}$

  $\begin{array}{rl}{b}_{j-1,1}(t)& =\frac{t-t_{j-1}}{t_j-t_{j-1}}{b}_{j-1,0}(t)+\frac{t_{j+1}-t}{t_{j+1}-t_j}{b}_{j,0}(t)\\ & =\frac{t_{j+1}-t}{t_{j+1}-t_j}\end{array}$

• 2 次 基函数有三个非零 $b_{j,2},b_{j-1,2},b_{j-2,2}$
  $\begin{array}{rl}{b}_{j,2}(t)& =\frac{t-t_j}{t_{j+2}-t_j}{b}_{j,1}(t)+\frac{t_{j+3}-t}{t_{j+3}-t_{j+1}}{b}_{j+1,1}(t)\\ & =\frac{t-t_j}{t_{j+2}-t_j}\cdot \frac{t-t_j}{t_{j+1}-t_j}\\ & =\frac{(t-t_j{)}^2}{(t_{j+2}-t_j)(t_{j+1}-t_j)}\end{array}$

  $\begin{array}{rl}{b}_{j-1,2}(t)& =\frac{t-t_{j-1}}{t_{j+1}-t_{j-1}}{b}_{j-1,1}(t)+\frac{t_{j+2}-t}{t_{j+2}-t_j}{b}_{j,1}(t)\\ & =\frac{t-t_{j-1}}{t_{j+1}-t_{j-1}}\cdot \frac{t_{j+1}-t}{t_{j+1}-t_j}+\frac{t_{j+2}-t}{t_{j+2}-t_j}\cdot \frac{t-t_j}{t_{j+1}-t_j}\\ & =\frac{(t-t_{j-1})(t_{j+1}-t)}{(t_{j+1}-t_{j-1})(t_{j+1}-t_j)}+\frac{(t_{j+2}-t)(t-t_j)}{(t_{j+2}-t_j)(t_{j+1}-t_j)}\end{array}$

  $\begin{array}{rl}{b}_{j-2,2}(t)& =\frac{t-t_{j-2}}{t_j-t_{j-2}}{b}_{j-2,1}(t)+\frac{t_{j+1}-t}{t_{j+1}-t_{j-1}}{b}_{j-1,1}(t)\\ & =\frac{t_{j+1}-t}{t_{j+1}-t_{j-1}}\cdot \frac{t_{j+1}-t}{t_{j+1}-t_j}\\ & =\frac{(t_{j+1}-t{)}^2}{(t_{j+1}-t_{j-1})(t_{j+1}-t_j)}\end{array}$

• 3 次 基函数有四个非零 $b_{j,3},b_{j-1,3},b_{j-2,3},b_{j-3,3}$

  $\begin{array}{rl}{b}_{j,3}(t)& =\frac{t-t_j}{t_{j+3}-t_j}{b}_{j,2}(t)+\frac{t_{j+4}-t}{t_{j+4}-t_{j+1}}{b}_{j+1,2}(t)\\ & =\frac{t-t_j}{t_{j+3}-t_j}\cdot \frac{(t-t_j{)}^2}{(t_{j+2}-t_j)(t_{j+1}-t_j)}\\ & =\frac{(t-t_j{)}^3}{(t_{j+3}-t_j)(t_{j+2}-t_j)(t_{j+1}-t_j)}\end{array}$

  $\begin{array}{rl}{b}_{j-1,3}(t)& =\frac{t-t_{j-1}}{t_{j+2}-t_{j-1}}{b}_{j-1,2}(t)+\frac{t_{j+3}-t}{t_{j+3}-t_j}{b}_{j,2}(t)\\ & =\frac{t-t_{j-1}}{t_{j+2}-t_{j-1}}\cdot [\frac{(t-t_{j-1})(t_{j+1}-t)}{(t_{j+1}-t_{j-1})(t_{j+1}-t_j)}+\frac{(t_{j+2}-t)(t-t_j)}{(t_{j+2}-t_j)(t_{j+1}-t_j)}]+\frac{t_{j+3}-t}{t_{j+3}-t_j}\cdot \frac{(t-t_j{)}^2}{(t_{j+2}-t_j)(t_{j+1}-t_j)}\\ & =\frac{(t-t_{j-1}{)}^2(t_{j+1}-t)}{(t_{j+2}-t_{j-1})(t_{j+1}-t_{j-1})(t_{j+1}-t_j)}+\frac{(t-t_{j-1})(t_{j+2}-t)(t-t_j)}{(t_{j+2}-t_{j-1})(t_{j+2}-t_j)(t_{j+1}-t_j)}+\frac{(t_{j+3}-t)(t-t_j{)}^2}{(t_{j+3}-t_j)(t_{j+2}-t_j)(t_{j+1}-t_j)}\end{array}$

  $\begin{array}{rl}{b}_{j-2,3}(t)& =\frac{t-t_{j-2}}{t_{j+1}-t_{j-2}}{b}_{j-2,2}(t)+\frac{t_{j+2}-t}{t_{j+2}-t_{j-1}}{b}_{j-1,2}(t)\\ & =\frac{t-t_{j-2}}{t_{j+1}-t_{j-2}}\cdot \frac{(t_{j+1}-t)(t_{j+1}-t)}{(t_{j+1}-t_{j-1})(t_{j+1}-t_j)}+\frac{t_{j+2}-t}{t_{j+2}-t_{j-1}}\cdot [\frac{(t-t_{j-1})(t_{j+1}-t)}{(t_{j+1}-t_{j-1})(t_{j+1}-t_j)}+\frac{(t_{j+2}-t)(t-t_j)}{(t_{j+2}-t_j)(t_{j+1}-t_j)}]\\ & =\frac{(t-t_{j-2})(t_{j+1}-t{)}^2}{(t_{j+1}-t_{j-2})(t_{j+1}-t_{j-1})(t_{j+1}-t_j)}+\frac{(t_{j+2}-t)(t-t_{j-1})(t_{j+1}-t)}{(t_{j+2}-t_{j-1})(t_{j+1}-t_{j-1})(t_{j+1}-t_j)}+\frac{(t_{j+2}-t{)}^2(t-t_j)}{(t_{j+2}-t_{j-1})(t_{j+2}-t_j)(t_{j+1}-t_j)}\end{array}$

  $\begin{array}{rl}{b}_{j-3,3}(t)& =\frac{t-t_{j-3}}{t_j-t_{j-3}}{b}_{j-3,2}(t)+\frac{t_{j+1}-t}{t_{j+1}-t_{j-2}}{b}_{j-2,2}(t)\\ & =\frac{t_{j+1}-t}{t_{j+1}-t_{j-2}}\cdot \frac{(t_{j+1}-t{)}^2}{(t_{j+1}-t_{j-1})(t_{j+1}-t_j)}\\ & =\frac{(t_{j+1}-t{)}^3}{(t_{j+1}-t_{j-2})(t_{j+1}-t_{j-1})(t_{j+1}-t_j)}\end{array}$

也即
$$\begin{array}{rl}& B_{j,3}(t)=\frac{(t-t_j{)}^3}{(t_{j+3}-t_j)(t_{j+2}-t_j)(t_{j+1}-t_j)}\\ & B_{j-1,3}(t)=\frac{(t-t_{j-1}{)}^2(t_{j+1}-t)}{(t_{j+2}-t_{j-1})(t_{j+1}-t_{j-1})(t_{j+1}-t_j)}+\frac{(t-t_{j-1})(t_{j+2}-t)(t-t_j)}{(t_{j+2}-t_{j-1})(t_{j+2}-t_j)(t_{j+1}-t_j)}+\frac{(t_{j+3}-t)(t-t_j{)}^2}{(t_{j+3}-t_j)(t_{j+2}-t_j)(t_{j+1}-t_j)}\\ & B_{j-2,3}(t)=\frac{(t-t_{j-2})(t_{j+1}-t{)}^2}{(t_{j+1}-t_{j-2})(t_{j+1}-t_{j-1})(t_{j+1}-t_j)}+\frac{(t_{j+2}-t)(t-t_{j-1})(t_{j+1}-t)}{(t_{j+2}-t_{j-1})(t_{j+1}-t_{j-1})(t_{j+1}-t_j)}+\frac{(t_{j+2}-t{)}^2(t-t_j)}{(t_{j+2}-t_{j-1})(t_{j+2}-t_j)(t_{j+1}-t_j)}\\ & B_{j-3,3}(t)=\frac{(t_{j+1}-t{)}^3}{(t_{j+1}-t_{j-2})(t_{j+1}-t_{j-1})(t_{j+1}-t_j)}\end{array}\text{\hspace{1em}(*)}$$
由此，得到该区间的曲线表达式为
$$\begin{array}{rl}B(t)& =B_{0,3}{P}_0+B_{1,3}{P}_1+\cdots +B_{n,3}{P}_n\\ & =B_{j-3,3}(t)P_{j-3}+B_{j-2,3}(t)P_{j-2}+B_{j-1,3}(t)P_{j-1}+B_{j,3}(t)P_j\end{array}$$
根据公式可知，当节点 $t_0,t_1,\cdots ,t_n$ 确定后，$B_{j,3},B_{j-1,3},B_{j-2,3},B_{j-3,3}$ 都是关于 t 的三次函数。$B(t)$ 也随着 t $（\in [0,1]）$ 的取值而确定，也即确定了曲线上的一系列点 B(t)，正是这些点组成了样条曲线。

2.2.2 B样条插值方程的获取

clamped -B样条本身跟贝塞尔曲线很接近，当然可以考虑使用类似于上节贝塞尔曲线插值那样，在两点内依次插入 B 样条，通过某些设定确保曲线光滑。但这种方法会同上面一样很粗糙。

对于 B 样条，不同与贝塞尔曲线的一点就是，B 样条曲线有节点，而且节点可以映射到最终曲线上，节点对应坐标为 $B(t_j)$，依据以上公式计算。而且B样条本身就是一条光滑曲线，取代于上面贝塞尔曲线的两点间逐个插值，考虑使用B样条一次性生成所需曲线，且曲线通过原始点（型值点）。

而要满足生成的 B 样条直接穿过型值点，恰好可以利用上述节点性质（即节点在曲线上的有对应点 $B(t_j)$）。因此只需要考虑让 型值点等值于 节点对应点即可。即便这样，我们仍然有很大的参数设定自由度。

如，假设有一系列型值点中考虑其中一点 P’，我们使用三次B样条对所有型值点进行插值，让 P’ 等值于第 i 个节点对应点 $B(t_j)$ ，即
$$P‘=B(t_i)=B_{j-3,3}(t_j)P_{j-3}+B_{j-2,3}(t_j)P_{j-2}+B_{j-1,3}(t_j)P_{j-1}+B_{j,3}(t_j)P_j$$
该式中，P‘ 已知，控制点 $P_{j-3},P_{j-2},Pj-1,Pj$ 都未知，第 j 个节点值 $t_j$ 也未知，而且这些值都可以一定程度地"随意设置。

也就是说，我们有很大的自由度来通过手动设定控制点、节点分布等来满足上式。而只要满足上式，就会有生成的B样条曲线经过型值点 P’。

当然这是对一个型值点而言，考虑所有型值点的情况下，每个公式中的控制点会有重复，手动设定的自由度可能会下降，但大致感觉最终应该还是会有很多自由度。当然，可以用数学公式具体来计算说明，在此不赘述，仅为了定性说明问题。下面有详细分析求解过程。

同样以 三次B样条为例，假设已知有 l+1 个型值点为 $z_0,z_1,\cdots ,z_l$ ；我们要插值生成的三次B样条参数有：$n+1$ 个控制点 $P_0,P_1,\cdots ,P_n$ 和 $m+1$ 个节点 $t_0,t_1,\cdots ,t_m$ 。

首先需要明确的一点是，就目前来说，目标三次样条曲线中控制点个数没限制，节点数目没限制，仅有的一个限制是
$$m=n+k+1=n+4\text{\hspace{1em}(1)}$$

对于每个没有限制的点（控制点、节点）都意味着一个个的自由度。下面开始施加限制

先考虑最重要的一点，即，使型值点与节点（确切说是，其在目标曲线上的点）对应，也即 $z_i=B(t_i)$。可以用如下图表示

上图中表示了，$z_0=B(t_0),z_1=B(t_1),\cdots ,z_n=B(t_m)$，此时是最简单的对应关系，此时有 $m=l$。但此时会有个问题，因为有一个不可忽视的一点，样条曲线是有定义域的，对于三次样条，定义域为 $[t_3,t_{n+1}]$，因此上述对应关系是有问题的，至少 $z_0$ 应该与定义域内第一个值 $B(t_3)$ 对应。同样，最后一个 $z_l$ 与 定义域内最后一个值 $B(t_{n+1})$ 对应。

上面提到了定义域的问题，一个极好的解决方案是使用 clamped B-spline

对于三次的 clamped B-样条，有
$$\begin{gathered} t_0=t_1=t_2=t_3=0 \\ {t}_4,t_5,\cdots ,t_n\in [0,1] \\ {t}_{n+1}=t_{n+2}=t_{n+3}=t_{n+4}=1 \end{gathered}$$
而定义域为 $[t_3,t_{n+1}]$ ，因此我们采用的对应关系为
$$\begin{gathered} t_3\iff z_0 \\ {t}_5\iff z_1 \\ ⋮ \\ {t}_{n+1}\iff z_l \end{gathered}$$
也即有
$$\begin{gathered} z_0=B(t_3)=B(0) \\ ⋮ \\ {z}_j=B(t_{j+3}) \\ ⋮ \\ {z}_l=B(t_{n+1})=B(1)\text{\hspace{1em}(2)} \end{gathered}$$
此时，$n+1=l+3$ 即满足 $n=l+2$。

综上所述，根据限制，我们需要满足的关系有：

• 需要设定样条的控制点的数目比型值点多2个，即
  $$n=l+2$$

• 需要设定的节点数目比型值点多6个，即
  $$m=n+k+1=l+2+3+1=l+6$$

• 满足一系列关系式（共 $l+1$ 个）
  $$z_{j-3}=B(t_j)=B_{j-3,3}(t_j)P_{j-3}+B_{j-2,3}(t_j)P_{j-2}+B_{j-1,3}(t_j)P_{j-1}+B_{j,3}(t_j)P_j,j=3,4,\cdots ,l+3$$

  且对于 $B_{j,3},B_{j-1,3},B_{j-2,3},B_{j-3,3}$ ，带入（*）公式得
  $$\begin{array}{rl}& B_{j-3,3}(t_j)=\frac{(t_{j+1}-t_j{)}^2}{(t_{j+1}-t_{j-2})(t_{j+1}-t_{j-1})}\\ & B_{j-2,3}(t_j)=\frac{(t_j-t_{j-2})(t_{j+1}-t_j{)}^2}{(t_{j+1}-t_{j-2})(t_{j+1}-t_{j-1})(t_{j+1}-t_j)}+\frac{(t_{j+2}-t_j)(t_j-t_{j-1})}{(t_{j+2}-t_{j-1})(t_{j+1}-t_{j-1})}\\ & B_{j-1,3}(t_j)=\frac{(t_j-t_{j-1}{)}^2(t_{j+1}-t_j)}{(t_{j+2}-t_{j-1})(t_{j+1}-t_{j-1})(t_{j+1}-t_j)}\\ & B_{j,3}(t_j)=\frac{(t-t_j{)}^3}{(t_{j+3}-t_j)(t_{j+2}-t_j)(t_{j+1}-t_j)}=0\end{array}$$

  式子 $B_{j,3}=0$ ，当 t 列表确定时，$B_{j-1,3},B_{j-2,3},B_{j-3,3}$ 都为定值。上述关系式变为
  $$z_{j-3}=B(t_j)=B_{j-3,3}(t_j)P_{j-3}+B_{j-2,3}(t_j)P_{j-2}+B_{j-1,3}(t_j)P_{j-1},j=3,4,\cdots ,l+3$$

  展开
  $$\begin{gathered} z_0=B_{0,3}(t_3)P_0+B_{1,3}(t_3)P_1+B_{2,3}(t_3)P_2 \\ {z}_1=B_{1,3}(t_4)P_1+B_{2,3}(t_4)P_2+B_{3,3}(t_4)P_3 \\ ⋮ \\ {z}_l=B_{l,3}(t_{l+3})P_l+B_{l+1,3}(t_{l+3}{P}_{l+1}+B_{l+2,3}(t_{l+3})P_{l+2} \end{gathered}$$

  用矩阵更简洁地表示为
  $$\left[\begin{array}{c}{z}_0\\ z_1\\ ⋮\\ z_l\end{array}\right]=\left[\begin{array}{cccccc}{B}_{0,3}(t_3)& B_{1,3}(t_3)& B_{2,3}(t_3)& 0& \cdots & 0\\ 0& B_{1,3}(t_4)& B_{2,3}(t_4)& B_{3,3}(t_4)& \cdots & 0\\ ⋮& 0& \ddots & & & 0\\ 0& \cdots & 0& B_{l,3}(t_{l+3})& B_{l+1,3}(t_{l+3})& B_{l+2,3}(t_{l+3})\end{array}\right]\left[\begin{array}{c}{P}_0\\ P_1\\ ⋮\\ P_{l+2}\end{array}\right]$$

  且已知
  $$\begin{gathered} t_0=t_1=t_2=t_3=0 \\ {t}_{l+3}=t_{l+4}=t_{l+5}=t_{l+6}=1 \end{gathered}$$
  则可得
  $$\begin{array}{rl}& B_{0,3}(t_3)=\frac{(t_4-t_3{)}^2}{(t_4-t_1)(t_4-t_2)}=1\\ & B_{1,3}(t_3)=\frac{(t_3-t_1)(t_4-t_3{)}^2}{(t_4-t_1)(t_4-t_2)(t_4-t_3)}+\frac{(t_5-t_3)(t_3-t_2)}{(t_5-t_2)(t_4-t_2)}=0\\ & B_{2,3}(t_3)=\frac{(t_3-t_2{)}^2(t_4-t_3)}{(t_5-t_2)(t_4-t_2)(t_4-t_3)}=0\end{array}$$
  那么，对于方程组第一个式子，就变为
  $$z_0=B_{0,3}(t_3)P_0+B_{1,3}(t_3)P_1+B_{2,3}(t_3)P_2=P_0$$
  反过来说，要使得生成曲线的第一个点就是控制点的第一个点，就是使得 $B_{0,3}(t_3)=1,B_{1,3}(t)=0,B_{2,3}(t_3)=0$ ，根据以上公式，也即令
  $$t_1=t_2=t_3$$
  且已经有 $t_0=0$

  理论上，只要满足
  $$\begin{gathered} t_0=0 \\ {t}_1=t_2=t_3 \\ \text{\hspace{1em}(\#1)} \end{gathered}$$
  即可确保三次B样条曲线，首端点经过第一个控制点。

  推广到 k 次应该是，只要确保
  $$\begin{gathered} t_0=0 \\ {t}_1=t_2=\cdots =t_k \end{gathered}$$
  即可有首端点为第一个控制点。

  这也是为什么首末节点重复 $k+1$ 次（clamped B-spline curves）后，样条曲线首末端点为控制点的首末两点的原因。

  末端点可类推，但会有所差别，差别之处在于 $B_{j-1,3}(t_j)$
  $$\begin{array}{rl}& B_{l,3}(t_{l+3})=0\\ & B_{l+1,3}(t_{l+3})=0\\ & B_{l+2,3}(t_{l+3})=\frac{(t_{l+3}-t_{l+2}{)}^2(t_{l+4}-t_{l+3})}{(t_{l+5}-t_{l+2})(t_{l+4}-t_{l+2})(t_{l+4}-t_{l+3})}=\frac{(t_{l+3}-t_{l+2}{)}^2\ast 0}{(t_{l+3}-t_{l+2}{)}^2\ast 0}=\frac{0}{0}\end{array}$$
  所以如果使得
  $$\begin{gathered} t_{l+3}=t_{l+4}=t_{l+5} \\ {t}_{l+6}=1\text{\hspace{1em}(\#2)} \end{gathered}$$
  此时式子中分母会有 0 出现。编程的时候应该特别注意。但实际上，上面内容节点域都表达的不太准确，确切说，不是$[t_j,t_{j+1}]$，而是 $[t_j,t_{j+1})$ 。上面不做更改，仅再此说明。

  另，有趣的一点是，完整的最后一个等式中
  $$z_l=B_{l,3}(t_{l+3})P_l+B_{l+1,3}(t_{l+3}{P}_{l+1}+B_{l+2,3}(t_{l+3})P_{l+2}+B_{l+3,3}(t_{l+3})P_{l+3}$$
  实际上 $P_{l+3}$ 是不存在的，因为 $n=l+2$ ，所以控制点仅有 $P_0,P_1,\cdots ,P_{l+2}$

  因此严格来说，最后一个等式不应写入方程，末端点亦不在节点域内（左闭右开）。但利用极限思想，应该可以得到，如果满足式子（#2），曲线末端点会无限接近末控制点。

  结合以上分析，且对于最后一个式子，无论分母为零或是控制点溢出都说明该式子不存在，但要映射末端点等价，因此可以将最后一个式子写为
  $$z_l=P_{l+2}$$
  矩阵方程变为
  $$\left[\begin{array}{c}{z}_0\\ z_1\\ ⋮\\ z_l\end{array}\right]=\left[\begin{array}{cccccc}1& 0& 0& 0& \cdots & 0\\ 0& B_{1,3}(t_4)& B_{2,3}(t_4)& B_{3,3}(t_4)& \cdots & 0\\ ⋮& 0& \ddots & & & 0\\ 0& \cdots & 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}{P}_0\\ P_1\\ ⋮\\ P_{l+2}\end{array}\right]$$

  总共有 $l+1$ 个方程；已知 $t_3=0,t_{l+3}=1$，还有 $l-1$ 个 t 未确定，分别为 $t_4,t_5,\cdots ,t_{l+2}$ ；$l+3$ 个控制点，即 $l+3$ 个未知数；

  所以，现在的需求即是，找到 $l-1$ 个未确定的 t 和 $l+$3 个未确定的 P，只要满足这 $l+1$ 个方程组即可。

  就直观来讲，答案似乎不是唯一的，甚至有很多解。而这一点也符合现实，因为，假设有三个型值点，过这三个型值点的B样条曲线应该不会只有一条。

  问题是如何求解该方程? 获得 t 和 P，从而确定样条曲线。

2.2.3 方程求解，确定样条曲线

已知型值点 $z_0,z_1,\cdots ,t_l$ ，求满足下方程组的 $t_4,t_5,\cdots ,t_{l+2}$ 和 $P_0,P_1,\cdots ,P_{l+2}$ ,进而确定最终的经过型值点的插值曲线（三次B样条）
$$\left[\begin{array}{c}{z}_0\\ z_1\\ ⋮\\ z_l\end{array}\right]=\left[\begin{array}{cccccc}1& 0& 0& 0& \cdots & 0\\ 0& B_{1,3}(t_4)& B_{2,3}(t_4)& B_{3,3}(t_4)& \cdots & 0\\ ⋮& 0& \ddots & & & 0\\ 0& \cdots & 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}{P}_0\\ P_1\\ ⋮\\ P_{l+2}\end{array}\right]$$
$l+1$ 个方程，很多未知数。不妨先确定节点分布，就使用均匀节点分布，在 2.2.1 节有描述，此时
$$\begin{array}{rl}& t_0=t_1=t_2=t_3=0\\ & t_4=\frac{1}{l}\\ & t_5=\frac{2}{l}\\ & \cdots \\ & t_j=\frac{j-3}{l}\\ & \cdots \\ & t_{l+2}=\frac{l-1}{l}\\ & t_{l+3}=t_{l+4}=t_{l+5}=t_{l+6}=1\end{array}$$
将此节点带入到 $B_{j,k}(t_j)$ 求解，上面方程组变为
$$\left[\begin{array}{c}{z}_0\\ z_1\\ ⋮\\ z_l\end{array}\right]=\frac{1}{6}\left[\begin{array}{ccccccc}6& 0& 0& 0& 0& \cdots & 0\\ 0& 1& 4& 1& 0& \cdots & 0\\ 0& 0& 1& 4& 1& \cdots & 0\\ ⋮& 0& & \ddots & & & 0\\ 0& 0& \cdots & 1& 4& 1& 0\\ 0& 0& \cdots & 0& 0& 0& 6\end{array}\right]\left[\begin{array}{c}{P}_0\\ P_1\\ ⋮\\ P_{l+2}\end{array}\right]$$
知系数矩阵的秩为 $l+1$ ，增广矩阵秩为 $l+1$，未知数个数为 $l+3$ 。

则知，即便在已经手动设定节点分布情况下，对于该方程组，未知数个数仍旧大于增广矩阵的秩，即仍旧有无穷多个解。确切的说，还缺少两个方程限制，使得方程有且仅有一解。

这两个方程如何获得，或者如何再对其施加两个限制呢？既然是两个，何不对两个端点施加限制呢。

现在考虑首末区间的两段曲线，曲线表达式分别为
$$\begin{array}{rl}& when,t_i\in [3,4)\\ B(t)& =B_{0,3}(t)P_0+B_{1,3}(t)P_1+B_{2,3}(t)P_2+B_{3,3}(t)P_3\\ & =\frac{(t_4-t{)}^3}{t_4^3}{P}_0+[\frac{t(t_4-4{)}^2}{t_4^3}+\frac{t(t_4-t)(t_5-t)}{t_4^2{t}_5}+\frac{t(t-t_5{)}^2}{t_4{t}_5^2}]P_1+[\frac{t^2(t_4-t)}{t_4^2{t}_5}+\frac{t^2(t_5-t)}{t_4{t}_5^2}+\frac{t^2(t_6-t)}{t_4{t}_5{t}_6}]P_2+\frac{t_3^3}{t_4{t}_5{t}_6}{P}_3\\ & when,t_i\in [l+2,l+3)\\ B(t)& =B_{l-1,3}(t)P_{l-1}+B_{l,3}(t)P_l+B_{l+1,3}(t)P_{l+1}+B_{l+2,3}(t)P_{l+2}\\ & =a{P}_{l-1}+b{P}_l+c{P}_{l+1}+d{P}_{l+2}\end{array}$$
说明：$t\in [l+2,l+3)$ 式子没有详细写出，键盘敲公式太麻烦了。总之有如下特点，第一段曲线，是由 $P_0,P_1,P_2,P_3$ 决定，且权重是关于 t 的三次方程；第二段曲线是由 $P_{l-1},P_l,P_{l+1},P_{l+2}$ 决定，同样权重是关于 t 的三次方程。

有一种常用的边界条件是非扭结边界条件，具体为：使第一个点的三次导数和第二点的三次导数一样，最后一个点的三次导数和倒数第一个点一样。

很明显，当对曲线求三次导后，表达式应该是关于 P 的常量线性组合，即权重系数变为常数。如果使用非扭结边界条件，就会有另外两个方程
$$\begin{gathered} a_3{P}_0+b_3{P}_1+c_3{P}_2+d_3{P}_3=a_4{P}_1+b_4{P}_2+c_4{P}_3+d_4{P}_4 \\ {a}_{l+1}{P}_{l-2}+b_{l+1}{P}_{l-1}+c_{l+1}{P}_l+d_{l+1}{P}_{l+1}=a_{l+2}{P}_{l-1}+b_{l+2}{P}_l+c_{l+2}{P}_{l+1}+d_{l+2}{P}_{l+2} \end{gathered}$$
其中 P 前系数都为常数，可具体计算出来，同样不给出具体值了。

所以，所有方程都具备了，样条曲线也就能求出来了。且是唯一解。那么如何解?以下仅描述自己思路

先将上面两个方程写入到本节开头的矩阵中去，毕竟矩阵解方程是线代基础，而且用起来熟悉且方便。整理非扭结边界的两方程
$$\begin{gathered} a_3{P}_0+(b_3-a_4)P_1+(c_3-b_4)P_2+(d_3-c_4)P_3-d_4{P}_4=0 \\ {a}_{l+1}{P}_{l-2}+(b_{l+1}-a_{l+2})P_{l-1}+(c_{l+1}-b_{l+2})P_l+(d_{l+1}-c_{l+2})P_{l+1}-d_{l+2}{P}_{l+2}=0 \end{gathered}$$
加入到矩阵方程：
$$\left[\begin{array}{c}0\\ z_0\\ z_1\\ ⋮\\ z_l\\ 0\end{array}\right]=\frac{1}{6}\left[\begin{array}{cccccccc}{a}_3& b_3-a_4& c_3-b_4& d_3-c_4& -d_4& \cdots & 0& 0\\ 6& 0& 0& 0& 0& \cdots & 0& 0\\ 0& 1& 4& 1& 0& \cdots & 0& 0\\ 0& 0& 1& 4& 1& \cdots & 0& 0\\ ⋮& 0& & \ddots & & & 0& 0\\ 0& 0& \cdots & 0& 1& 4& 1& 0\\ 0& 0& \cdots & 0& 0& 0& 0& 6\\ 0& 0& \cdots & a_{l+1}& b_{l+1}-a_{l+2}& c_{l+1}-b_{l+2}& d_{l+1}-c_{l+2}& -d_{l+2}\end{array}\right]\left[\begin{array}{c}{P}_0\\ P_1\\ ⋮\\ P_{l+2}\end{array}\right]$$
$l+3$ 个方程，$l+3$ 个未知控制点。解方程即可的到所需插值曲线。

以上便是 B 样条用于插值的原理与解决方案。

关于方程的具体解法，可能有巧妙的算法，如果没有也能够使用线代中最朴实方法，通过变换增广矩阵来求出最终每个控制点的表达式。当然是关于 z 坐标的求解表达式。既然表达式算出来了，也就可以通过编程实现，在已知型值点 z 坐标的情况下，得到其3次B样条插值曲线。

已经迷糊于密集的数学字符、简单却繁琐的数学计算、麻烦的表达式输入，不再写代码实现。原理如此，实现不难。

3. 代码实现（挖坑）

仪式且礼貌性的一节

3. 三次样条插值

上面说明了三次B样条插值，而三次样条插值则思路直接且简单的多。也在此作个简要概述。

这里的三次样条，就是满足三次函数的曲线，即
$$y=a{t}^3+b{t}^2+ct+d$$
这简直比 B 样条直接多了。问题解决思路也直接有效。

以实例说明，如已知四（$orl$）个型值点，需要插值出以下的曲线。每两个点间插入一条三次样条，因此四个点需要三$(thatis,l-1)$个三次样条。每个三次样条有四个未知数$(a,b,c,d)$，所以一共有12$(thatis,4(l-1))$个未知数，因此需要12$(thatis,4(l-1))$个方程求解。

下面就是找方程，找限制条件：

• 每个样条曲线都经过其曲线所在的两端点。这就有 6 个方程了。而这一条件，也确保了 $c^0$ 连续。$(thatis,2(l-1))$
• 在每个连接点处，应该够平滑。那就让其在连接点的左右一阶导数相等。这就又有了 2 个方程。这一条件，确保了 $c^1$ 连续。$(thatis,(l-2))$
• 贝塞尔曲线中展示了 $c^2$ 连续的重要性，就在确保下 $c^2$ 连续。即连接点左右二阶导数想等。这就又有了 2 个方程。$(thatis,(l-2))$

上面确保了三个连续平滑条件，一共构建了 10 个方程$(thatis,4l-6)$ ，还差两个方程。还真跟上节B样条情况类似，缺少两个方程。其实条件已经限制的差不多了。如何解决呢？

• 可以考虑使用B样条中的非扭结边界，即 $a_1=a_2,a_2=a_3$

• 既然连接点的一、二阶导数都限制了，那就对两个端点处也限制一下嘛。一阶导数不好限制，因为直接设定一定值的话，也就确定了斜率，可能会极大地影响形状，那就限制二阶导数。二阶有点麻烦，就另三阶为 0 ，即 $a_1=0,a_3=0$

• 实在没有好的办法，就随意指定两个参数的值。可能也不失为一种方法。

方程确定了也就可以求解出各个样条了。

4. 总结

对贝塞尔曲线插值、B样条曲线插值、一般样条曲线插值进行了详细地理论挖掘。三种插值是共通的，毕竟都同属样条曲线，印象最深的几点在于

• B样条插值一步到位，而贝塞尔曲线、一般样条插值是逐段生成插值曲线
• 插值本质都是列方程、解方程，从而确定插值曲线
• 方程限制一般都是围绕 $c^0,c^1,c^2$ 连续展开
• 在三种连续设定完毕后，一般还会缺少两个限制条件，此时可以从端点入手。至于限制方法，可以灵活设定。

此致，终于完结了。。。