一、局部曲线插值预备知识

给定 { Q k } ( k = 0 , 1 , ⋯ , n ) \{\pmb Q_k\}(k=0,1,\cdots,n) {Qk​}(k=0,1,⋯,n)，局部曲线插值是指创建 $n$ 条多项式或有理曲线段：$C_i(u)(i=0,1,\cdots ,n-1)$，要求相邻曲线段在连接点处满足事先指定的连续阶，采用多项式构造每个曲线段，然后选择一个合适的节点矢量得到 B 样条曲线。

为了生成 B 样条曲线段 $C_i(u)$，需要计算内控制点：二次的有一个；三次的有两个。这些控制点都在与曲线相切于 $Q_k$ 点的直线上，因此我们需要知道每一个 $Q_k$ 处的切向矢量 $T_k$。
$$\begin{array}{rl}& q_k=Q_k-Q_{k-1}\\ & V_k=(1-{\alpha }_k)q_k+{\alpha }_k{q}_{k+1}\\ & T_k=\frac{V_k}{\mid V_k\mid }\end{array}\text{\hspace{1em}(1)}$$

其中，$T$ 表示单位长度切矢量。插值参数 ${\alpha }_k$ 的计算方法如下（五点法）：
$${\alpha }_k=\frac{\mid q_{k-1}\times q_k\mid }{\mid q_{k-1}\times q_k\mid +\mid q_{k+1}\times q_{k+2}\mid },k=2,\cdots ,n-2\text{\hspace{1em}(2)}$$

需要注意的是：当（2）式分母为 0 意味着要么 $Q_k$ 是一个角点，要么 $Q_{k-2}$ 到 $Q_{k+2}$ 是一条直线段。在这种情况下，${\alpha }_k$ 的值有多种选择方式，可以按如下方式选取：

• ${\alpha }_k=1$，这时有 $V_k=q_{k+1}$；它会在 $Q_k$ 处产生一个角点（当希望在 $Q_k$ 处保留角点时）；
• ${\alpha }_k=\frac{1}{2}$，有 $V_k=\frac{1}{2}(q_k+q_{k+1})$；这种选择将会消除掉角点（使其变光滑）（当不希望保留角点时）。

需要对端点处进行特殊的处理：
$$q_0=2q_1-q_2,q_{-1}=2q_0-q_1,q_{n+1}=2q_n-q_{n-1},q_{n+2}=2q_{n+1}-q_n\text{\hspace{1em}(3)}$$

将（3）式代入（1）式和（2）式，即可得到端点处的切矢 $T_0,T_1$ 和 $T_{n-1},T_n$。

def TangentialVector(Q, flag):"""计算数据点 Q 每一点处的切矢:param Q: 数据点:param flag: 0:保留角点；1：消除角点:return: 切矢 T""""q_k = Q_k - Q_{k-1}  q 的下标为 -1 ~ n+2"n = Q.shape[1] - 1  # 数据点 Q 的最大下标q = np.zeros((Q.shape[0], n + 3))  # q 应该有 n + 4 个列元素，缺少一个 q[-1] 为了统一下标，后续会额外增加一个 q_neg"计算 q[1]~q[n]"for i in range(1, n + 1):  # 计算q[:, i] = Q[:, i] - Q[:, i - 1]q[:, 0] = 2 * q[:, 1] - q[:, 2]  # q[0]q[:, n + 1] = 2 * q[:, n] - q[:, n - 1]  # q[n+1]q[:, n + 2] = 2 * q[:, n + 1] - q[:, n]  # q[n+2]q_neg = 2 * q[:, 0] - q[:, 1]  # q[-1]"五点法计算插值参数 alpha[k]"V = np.zeros((Q.shape[0], n + 1))T = np.zeros((Q.shape[0], n + 1))alpha = np.zeros(n + 1)for k in range(1, n + 1):down = np.linalg.norm(np.cross(q[:, k - 1], q[:, k])) + np.linalg.norm(np.cross(q[:, k + 1], q[:, k + 2]))up = np.linalg.norm(np.cross(q[:, k - 1], q[:, k]))if down == 0:  # Q[k] 是一个角点，或者，Q[k-2] 到 Q[k+2] 是一条直线段if flag == 0:  # 保留角点V[:, k] = q[:, k + 1]T[:, k] = V[:, k] / np.linalg.norm(V[:, k])elif flag == 1:  # 消除角点V[:, k] = (q[:, k] + q[:, k + 1]) / 2T[:, k] = V[:, k] / np.linalg.norm(V[:, k])else:print("参数 flag 输入错误")else:alpha[k] = up / downV[:, k] = (1 - alpha[k]) * q[:, k] + alpha[k] * q[:, k + 1]T[:, k] = V[:, k] / np.linalg.norm(V[:, k])parm_up = np.linalg.norm(np.cross(q_neg, q[:, 0]))parm_down = (np.linalg.norm(np.cross(q_neg, q[:, 0])) + np.linalg.norm(np.cross(q[:, 1], q[:, 2])))if parm_down == 0:  # Q[k] 是一个角点，或者，Q[k-2] 到 Q[k+2] 是一条直线段if flag == 0:  # 保留角点V[:, 0] = q[:, 1]T[:, 0] = V[:, 0] / np.linalg.norm(V[:, 0])elif flag == 1:  # 消除角点V[:, 0] = (q[:, 0] + q[:, 1]) / 2T[:, 0] = V[:, 0] / np.linalg.norm(V[:, 0])else:print("参数 flag 输入错误")else:alpha[0] = parm_up / parm_downV[:, 0] = (1 - alpha[0]) * q[:, 0] + alpha[0] * q[:, 1]T[:, 0] = V[:, 0] / np.linalg.norm(V[:, 0])return T

二、局部抛物线插值

给定 $xy$ 平面上的点 { Q k } , k = 0 , 1 , ⋯ , n \{\pmb Q_k\},k=0,1,\cdots,n {Qk​},k=0,1,⋯,n，用上一节的方法来计算相应的 { T k } \{\pmb T_k\} {Tk​}。用 $L_k$ 表示由 $(Q_k,T_k)$ 定义的有向直线，用 $R_k$ 表示 $L_{k-1}$ 和 $L_k$ 的交点
$$R_k=Q_{k-1}+{\gamma }_{k-1}{T}_{k-1},R_k=Q_k+{\gamma }_k{T}_k\text{\hspace{1em}(4)}$$

这里需要保证
$${\gamma }_{k-1}>0,{\gamma }_k<0\text{\hspace{1em}(5)}$$

这样才能保证沿着切矢方向上存在交点。（稍后会补充不满足该条件时如何求 $R_k$）

若 $p$ 次贝塞尔曲线是定义在 U = { 0 , ⋯ , 0 , 1 , ⋯ , 1 } U=\{0,\cdots,0,1,\cdots,1\} U={0,⋯,0,1,⋯,1}（没有内节点）上的 B 样条曲线，则端点处的一阶导矢计算公式如下所示：
$$\begin{array}{rl}& C^{\prime }(0)=Q_0=\frac{p}{u_{p+1}}(P_1-P_0)\\ & C^{\prime }(1)=Q_{n-1}=\frac{p}{1-u_{m-p-1}}(P_n-P_{n-1})\end{array}\text{\hspace{1em}(6)}$$