点线面的插值

何为法律？永恒不变的理性，普适万物的规则。

# 曲线、曲面与多项式

对空间中的一条曲线或曲面，我们往往用参数方程来表示：

$$\mathbf{p}(u)=\begin{pmatrix} x(u)\\ y(u)\\ z(u) \end{pmatrix} \quad \mathbf{p}(u,v)=\begin{pmatrix} x(u,v)\\ y(u,v)\\ z(u,v) \end{pmatrix}$$

我们往往认为$u,v$ 是一个从 0 变化到 1 的实数标量，那么前者$\mathbf{p}(u)$ 是一条曲线，后者$\mathbf{p}(u,v)$ 是一个曲面。

然而这是一般情况。在实际应用中，相比于一般函数$x(u),y(u),z(u)$，我们更关心它们是多项式时的情况。参数方程改写为：

$$\mathbf{p}(u)=\sum_{k=0}^n u^k\mathbf{c}_k\quad \mathbf{p}(u,v)=\sum_{k=0}^n\sum_{l=0}^m u^k v^l\mathbf{c}_{k,l}$$

其中$\mathbf{c}_{k},\mathbf{c}_{k,l}$ 都是三维的点或向量，严格来说$\mathbf{c}_0$ 应为点，其余是向量。

# 插值

插值问题就是给出一组点$\mathbf{p}_0,\mathbf{p}_1,...,\mathbf{p}_n$，要求你构造一条曲线去经过这些点。实际中往往采用三次函数来插值，因为二次往往有些欠拟合，而更高次又增加了计算量，是一个经验的折中。

譬如给出四个点$\mathbf{p}_0,\mathbf{p}_1,\mathbf{p}_2,\mathbf{p}_3$，那么考虑先前的参数方程，一个要求是我们希望这四个点均匀地分布在曲线上，即：

$$\begin{cases} \mathbf{p}(0)=\mathbf{p}_0\\ \mathbf{p}(1/3)=\mathbf{p}_1\\ \mathbf{p}(2/3)=\mathbf{p}_2\\ \mathbf{p}(1)=\mathbf{p}_3 \end{cases}$$

这已经有四个约束了，写成矩阵形式:

$$\begin{pmatrix} \mathbf{p}_0\\ \mathbf{p}_1\\ \mathbf{p}_2\\ \mathbf{p}_3\end{pmatrix}=\begin{pmatrix} 1&&&\\ 1&\frac{1}{3}&\frac{1}{3}^2&\frac{1}{3}^3\\ 1&\frac{2}{3}&\frac{2}{3}^2&\frac{2}{3}^3\\ 1&1&1&1 \end{pmatrix} \begin{pmatrix} \mathbf{c}_0\\ \mathbf{c}_1\\ \mathbf{c}_2\\ \mathbf{c}_3 \end{pmatrix}\\ \Rightarrow\mathbf{p}=A\mathbf{c}$$

然后我们就可以解出$\mathbf{c}=A^{-1}\mathbf{p}$，然后就得到了参数化的曲线方程$\mathbf{p}(u)=\mathbf{u}^T\mathbf{c}$，其中$\mathbf{u}^T=(1,u,u^2,u^3)$。

再简化写一下，我们称$\mathbf{b}(u)=(A^{-1})^T\mathbf{u}$ 为调和多项式 (blending polynomial)，那么就可以直接写出参数化曲线和四个点的关系：

$$\begin{aligned} \mathbf{p}(u)&=\mathbf{u}^T\mathbf{c}\\ &=\mathbf{u}^TA^{-1}\mathbf{p}\\ &=\mathbf{b}(u)^T\mathbf{p} \end{aligned}$$

---

讨论完曲线可以讨论三次插值曲面。譬如给出了 16 个点：

我们希望插值出一个三次曲面平均地通过它们，即：

$$\mathbf{p}(\frac{i}{3},\frac{j}{3})=\mathbf{p}_{ij}\\ i,j=0,1,2,3$$

记$\mathbf{u}=(1,u,u^2,u^3)^T,\mathbf{v}=(1,v,v^2,v^3)^T$，于是有：

$$\mathbf{p}(u,v)=\mathbf{u}^TC\mathbf{v}$$

其中，每个$C$ 的元素$C_{ij}$ 都是一个三维向量。然后我们有 16 个约束函数$\mathbf{p}(i/3,j/3)=\mathbf{p}_{ij}$，写成矩阵就是：

$$A_iCA_j^T=\mathbf{p}_{ij},i,j=0,1,2,3\\ \Rightarrow ACA^T=\mathbf{p}$$

于是解得$C=A^{-1}\mathbf{p}(A^T)^{-1}$。所以：

$$\begin{aligned} \mathbf{p}(u,v)&=(\mathbf{u}^TA^{-1})\mathbf{p}((A^T)^{-1}\mathbf{v})\\ &=\mathbf{b}(u)^T \mathbf{p} \mathbf{b}(v)\\ \end{aligned}$$

---

然后讨论一下 Hermite 插值，相较于给出四个点（也就是四个约束），Hermite 插值给出了两个端点，以及两个端点的导数，即给出了$\mathbf{p}_0,\mathbf{p}_1,\mathbf{p}_0',\mathbf{p}_1'$，同样也是四个约束。三次 Hermite 插值也是一个三次多项式，满足下列四个约束：

$$\begin{cases} \mathbf{p}(0)=\mathbf{c}_0&=\mathbf{p}_0\\ \mathbf{p}(1)=\mathbf{c}_0+\mathbf{c}_1+\mathbf{c}_2+\mathbf{c}_3 &=\mathbf{p}_1\\ \mathbf{p}'(0)=\mathbf{c}_1&=\mathbf{p}_0'\\ \mathbf{p}'(1)=\mathbf{c}_1+2\mathbf{c}_2+3\mathbf{c}_3&=\mathbf{p}_1' \end{cases}$$

这样也可以解方程，反解出$\mathbf{c}_0,...,\mathbf{c}_3$ 即可。

# 参数连续性和几何连续性

考虑两个参数化的曲线$\mathbf{p}(u),\mathbf{q}(v)$，如果我们要它们接个头即$\mathbf{p}(1)=\mathbf{q}(0)$，然后可以考虑接头处的连续性：

• 参数连续性，记为$C^n$:

$$C^0:\begin{cases}x_p(1)=x_q(0)\\ y_p(1)=y_q(0)\\ z_p(1)=z_q(0)\end{cases}\quad C^1: \begin{cases}x_p(1)=x_q(0)\\ y_p(1)=y_q(0)\\ z_p(1)=z_q(0)\\x_p'(1)=x_q'(0)\\ y_p'(1)=y_q'(0)\\ z_p'(1)=z_q'(0)\end{cases}\quad C^2:...$$

看一下规律就是到是要求$n$ 阶导数对应的参数完全相同。很显然要求的$n$ 越大，则约束条件越多，需要的多项式次数也就越高。

• 几何连续性，记为$G^n$:

$$G^0:\mathbf{p}(1)=\alpha \mathbf{q}(0)\quad G^1:\begin{cases} \mathbf{p}(1)=\alpha \mathbf{q}(0)\\ \mathbf{p}'(1)=\alpha \mathbf{q}'(0) \end{cases}\quad G^2:...$$

看一下实际上会发现几何连续性比参数连续性实际上少了一个约束，即只需要连接处的微矢是共线的就好。特别地，$G^0$ 和$C^0$ 是等价的。

# 三次 Bezier 曲线和 Bezier 曲面

Bezier 曲线是一种特殊的 Hermite 曲线，他给出了四个点$\mathbf{p}_0,\mathbf{p}_1,\mathbf{p}_2,\mathbf{p}_3$ 和如下约束：

$$\begin{cases} \mathbf{p}(0)=\mathbf{c}_0&=\mathbf{p}_0\\ \mathbf{p}(1)=\mathbf{c}_0+\mathbf{c}_1+\mathbf{c}_2+\mathbf{c}_3 &=\mathbf{p}_1\\ \mathbf{p}'(0)=\mathbf{c}_1&=3\mathbf{p}_1-3\mathbf{p}_0\\ \mathbf{p}'(1)=\mathbf{c}_1+2\mathbf{c}_2+3\mathbf{c}_3&=3\mathbf{p}_3-3\mathbf{p}_2 \end{cases}$$

我们可以解得调和多项式：

$$\mathbf{b}(u)=\begin{pmatrix} (1-u)^3\\ 3u(1-u)^2\\ 3u^2(1-u)\\ u^3 \end{pmatrix}$$

于是$\mathbf{p}(u)=\mathbf{b}(u)^T\mathbf{p}$。

同理推广到曲面的话，有$\mathbf{p}(u,v)=\mathbf{b}(u)^T P\mathbf{b}(v)$，其中$P_{ij}=\mathbf{p}_{ij}$。对于曲面，一次偏导$\frac{\partial\mathbf{p}}{\partial u}(0,0)=3(\mathbf{p}_{10}-\mathbf{p}_{00})$ 和 Bezier 曲线含义相同，但是二阶偏导$\frac{\partial^2\mathbf{p}}{\partial u\partial v}$ 的数值其实反映了曲面的曲率，即弯曲程度。

很显然，点$\mathbf{p}(u,v)$ 的法向量为:$\frac{\partial\mathbf{p}}{\partial u}\times \frac{\partial\mathbf{p}}{\partial v}$。

# 三次 B 样条曲线

三次 B 样条曲线给出了四个点$\mathbf{p}_0,\mathbf{p}_1,\mathbf{p}_2,\mathbf{p}_3$ 和如下约束：

$$\begin{cases} \mathbf{p}(0)=\mathbf{c}_0&=\frac{1}{6}(\mathbf{p}_0+4\mathbf{p}_1+\mathbf{p}_2)\\ \mathbf{p}(1)=\mathbf{c}_0+\mathbf{c}_1+\mathbf{c}_2+\mathbf{c}_3 &=\frac{1}{6}(\mathbf{p}_1+4\mathbf{p}_2+\mathbf{p}_3)\\ \mathbf{p}'(0)=\mathbf{c}_1&=\frac{1}{2}(\mathbf{p}_2-\mathbf{p}_0)\\ \mathbf{p}'(1)=\mathbf{c}_1+2\mathbf{c}_2+3\mathbf{c}_3&=\frac{1}{2}(\mathbf{p}_3-\mathbf{p}_1) \end{cases}$$

解出来可得：

$$\mathbf{b}(u)=\frac{1}{6}\begin{pmatrix} (1-u)^3\\ 4-6u^2+3u^3\\ 1+3u+3u^2-3u^3\\ u^3 \end{pmatrix},\quad \mathbf{p}(u)=\mathbf{b}(u)^T\mathbf{p}$$

这样做的话，使得曲线会被包裹在四个控制点的凸包内部。而且如果有一系列控制点$\mathbf{p}_0,\mathbf{p}_1,...,\mathbf{p}_n$，我们可以先用$(\mathbf{p}_0,\mathbf{p}_1,\mathbf{p}_2,\mathbf{p}_3)$ 生成曲线，然后再用$(\mathbf{p}_1,\mathbf{p}_2,\mathbf{p}_3,\mathbf{p}_4)$ 生成曲线，以此类推，这样的曲线具有较好的连续性。

对于一系列连续的控制点$\mathbf{p}_0,\mathbf{p}_1,...,\mathbf{p}_n$，我们单独考虑一个控制点$\mathbf{p}_i$，它为整条曲线做出的贡献只有四段：

$$(\mathbf{p}_{i-3}, \mathbf{p}_{i-2}, \mathbf{p}_{i-1}, \mathbf{p}_{i})\\ (\mathbf{p}_{i-2}, \mathbf{p}_{i-1}, \mathbf{p}_{i}, \mathbf{p}_{i+1})\\ (\mathbf{p}_{i-1}, \mathbf{p}_{i}, \mathbf{p}_{i+1}, \mathbf{p}_{i+2})\\ (\mathbf{p}_{i}, \mathbf{p}_{i+1}, \mathbf{p}_{i+2}, \mathbf{p}_{i+3})$$

如果每一段的曲线参数$u$ 都有一个单位区间，那么$n-1$ 段曲线$u$ 的范围就是$[0,n-1]$。$\mathbf{p}_i$ 的贡献可以用一个函数来表示：

$$B_i(u)=\begin{cases} 0 & u<i-2\\ b_0(u+2) & i-2\leq u<i-1\\ b_1(u+1) & i-1\leq u<i\\ b_2(u) & i\leq u<i+1\\ b_3(u-1) & i+1\leq u<i+2\\ 0 & u\geq i+2 \end{cases}$$

然后有：

$$\mathbf{p}(u)=\sum_{i=0}^nB_i(u)\mathbf{p}_i$$

这样就不用一段一段计算，然后可能一个控制点要在四段中都要出现的麻烦情况了。

定性分析的话，可以看出一个控制点$\mathbf{p}_i$ 对整条曲线的贡献函数：

显然就是曲线离控制点越近，受控制点影响越大。

当然，也可以推广到 B-Spline 曲面：

$$\mathbf{p}(u,v)=\sum_{i=0}^3\sum_{j=0}^3b_i(u)b_j(v)\mathbf{p}_{i,j}$$

# 一般的 B-spline 曲线

现在考虑这样一个问题，给了一系列控制点$\mathbf{p}_0,...,\mathbf{p}_m$，我们想寻找一条曲线$\mathbf{p}(u),u\in[u_{min},u_{max}]$。使得这条曲线 "平滑、并且接近这些控制点"。假设我们有一系列$u$ 的值，称为结点（knots）：

$$u_{min}=u_0\leq u_1\leq...\leq u_n=u_{max}$$

然后我们就把整条曲线分成了$n$ 段，每段为一个$d$ 次多项式：

$$\mathbf{p}(u_k<u<u_{k+1})=\sum_{j=0}^d\mathbf{c}_{jk}u^j\\$$

类似先前的想法，与其去考虑求解$\mathbf{c}_{jk}$，我们直接关心每个控制点对曲线的贡献。即写成：

$$\mathbf{p}(u)=\sum_{i=0}^m B_{id}(u)\mathbf{p}_i$$

其中，$B_{id}(u)$ 表示它是一个$d$ 次多项式。

实际上，B-spline 的含义就是 Basis Splines，而$\{B_{id}(u)\}$ 就是一组 basis。尽管这样的 basis 有很多，但是最著名的一个应该是 Cox-deBoor recursion 定义的：

$$\begin{aligned} & B_{k0}=\begin{cases} 1 & u_k\leq u<u_{k+1}\\ 0 & otherwise \end{cases}\\ & B_{kd}=\frac{u-u_k}{u_{k+d}-u_k}B_{k,d-1}(u)+\frac{u_{k+d}-u}{u_{k+d+1}-u_{k+1}}B_{k+1,d-1}(u) \end{aligned}$$

当$d=0$ 时，每个控制点$\mathbf{p}_i$ 只能影响一个区间$[u_k,u_{k+1}]$，实际上，每个控制点可以影响$d+1$ 个区间，取决于 basis 的次数。我们可以画出 0 次、1 次、2 次的 basis 函数图像：

控制点$\mathbf{p}_i$ 会影响区间：$[u_i,u_{i+1}],[u_{i+1},u_{i+2}],...,[u_{i+d},u_{i+d+1}]$。

这部分结点数$n+1$ 和控制点数$m+1$ 是没什么关系的，而且结点的选择也是自由的，这样的话就比较容易导致混乱。

所以实际中我们常常要求结点是等距选取的 (或重叠)，此时形成了 uniform B-spline。对于三次的 B-spline，我们常常取结点为：

$$\{0,0,0,0,1,2,...,n-1,n,n,n,n\}$$

而且由于一个控制点可以影响$d+1$ 个区间，所以实际上控制点数$m+1$、区间数$n$、basis 多项式次数$d$ 有以下关系：

$$n-(m+1)=d$$

至此，构造一个 B-spline 的方法比较完备了：

1. 给出控制点$\mathbf{p}_0,...,\mathbf{p}_m$，以及多项式次数$d$。
2. 选取结点$\{u_0,...,u_{d+m+1}\}$。
3. 计算出 basis 多项式$\{B_{id}(u)\}$。
4. 计算出曲线$\mathbf{p}(u)=\sum_{i=0}^m B_{id}(u)\mathbf{p}_i$。

特别地，若结点为$\{0,0,0,0,1,1,1,1\}$ 将退化成 Bezier 曲线。

此外，如果我们想特别地增加或减小某个控制点的影响力，也可以增加一个权重函数$w_i$，生成的曲线为 Nonuniform rational B-spline( NURBS )：

$$\mathbf{p}(u)=\frac{\sum_{i=0}^m B_{id}(u)w_i\mathbf{p}_i}{\sum_{i=0}^m B_{id}(u)w_i}$$

# Catmull-Rom Splines

如果我们放宽限制，不要求曲线必须在控制点的凸包内的话，Catmull-Rom spline 也是一种简单、受欢迎的方法。它给出的限制条件为：

$$\begin{cases} \mathbf{p}(0)=\mathbf{c}_0&=\mathbf{p}_1\\ \mathbf{p}(1)=\mathbf{c}_0+\mathbf{c}_1+\mathbf{c}_2+\mathbf{c}_3 &=\mathbf{p}_2\\ \mathbf{p}'(0)=\mathbf{c}_1&=\frac{1}{2}(\mathbf{p}_2-\mathbf{p}_0)\\ \mathbf{p}'(1)=\mathbf{c}_1+2\mathbf{c}_2+3\mathbf{c}_3&=\frac{1}{2}(\mathbf{p}_3-\mathbf{p}_1) \end{cases}$$

# 插值曲线的渲染

相较于直接去计算矩阵和多项式，有一种更为优雅的方式去渲染 Bezier 曲线。
对于四个控制点$\mathbf{p}_0, \mathbf{p}_1, \mathbf{p}_2,\mathbf{p}_3$ 和曲线$\mathbf{p}(u),u\in[0,1]$，我们可以考虑构造两条子曲线，即$\mathbf{l}(u),u\in[0,1]$ 和$\mathbf{r}(u),u\in[0,1]$，分别对应原曲线$[0,\frac{1}{2}],[\frac{1}{2},1]$ 两部分。

构造方法也很简单，一图表示为：

如上图取得都是中点，然后$\mathbf{l}_0,\mathbf{l}_1,\mathbf{l}_2,\mathbf{l}_3$ 和$\mathbf{r}_0,\mathbf{r}_1,\mathbf{r}_2,\mathbf{r}_3$ 分别对应了左半边曲线的控制点和右半边曲线的控制点。这样我们就可以递归地分治去构造子曲线了，然后把它们首尾相连起来。

其实这样做法得正确性也很好证明，只需要验证 Bezier 曲线的约束，两条子曲线拼成的曲线和原曲线等价即可。注意参数缩减了一般。所以$3(\mathbf{l}_1-\mathbf{l}_0)$ 才是$\mathbf{p}'(0)$。