TL;DR;

在计算机图形学中，插值（Interpolation）和拟合（Curve Fitting）在表面建模，路径规划，粒子系统等子领域中应用广泛， 且具有相似形式。这两种方法都是从一个点序列（Point Sequence）生成曲面。区别在于 Interpolation 一定会经过所有 Point， 而 Curve Fitting 的 Point 会分布在生成的曲线两侧。

为了计算方便，大多数工业中的插值和拟合方法都是分段多项式函数（样条）来表示的。使用样条可以避免出现著名的 Runge’s phenomenon ，同时，样条还具有局部性，也就是对局部的插值函数的修改不会影响到整体，对于建模和工业设计的局部调整意义重大。

一般来说，样条曲线方程不高于三次，因为较低次数的曲线能达到更好的局部拟合效果，而三次方程是实现曲线 \( G^2 \) 连续的最低幂次。也就是说，图形学/工业设计上遇到的曲线基本都是分段三次曲线，因此我们首先关注三次样条曲线（Cubic Spline）的一些共有性质。

在大学本科教材“数值分析”中，一个插值n个节点的三次样条曲线的形式是

$$ S(x) = S_1(x), x \in [x_1, x_{2}] \\\ … \\\ S(x) = S_k(x), x \in [x_k, x_{k+1}] \\\ S(x) = S_n(x), x \in [x_{n-1}, x_n] \\\ $$

$$ S_k(x) = y_k + b_k(x - x_k) + c_k(x - x_k)^2 +d_k(x - x_k)^3 $$

我们知道单个样条实际上就是一个截断的三次曲线，每个曲线段有4个参数 \( y_k, b_k, c_k, d_k \)，因此一共需要给出 \( 4n - 4 \) 个 参数，用以控制每个分段的曲线形状。如果要求 \( C^0 \) 连续即曲线通过所有数据点\( y_k \), 可以给出 \( 2n - 2 \) 个条件, 如果要求 \( C^1 \) 连续 和 \( C^2 \) 连续可以给出 \( 2n-4 \)个条件, 因此我们一共有 4n - 6 个方程组，却要求解 4n - 4 个条件。多出的两个自由的使得三次样条的解有无数多个

显然，这个一般形式并不能直接在建模中直接使用。首先其求解比较复杂，其次很多情况下我们并不需要 \( C^2 \) 的连续级别，因此在图形学中往往采用更直接的方式来给出参数。

另外，样条形式不仅仅可以用于插值，还可以用于拟合，这种情况下虽然曲线不经过数据点，但往往通过一些权重函数来组合数据点，在拟合中，这些点一般称为控制点

在图形学中广泛应用的插值和拟合方法分别是 Cubic Hermite Interpolation 和 Non Uniform Rational B Spline(NURBS)。本文将分别介绍

Cubic Hermite Interpolation

在计算机图形学中应用最广的是插值方法是 Cubic Hermite 插值，这种插值方法需要提供节点的坐标和节点的一阶导数。其好处是节点坐标和节点导数都可以通过图形界面来调整。且由于节点的导数和坐标已经确定，Cubic Hermite 插值是天然 \( G^1 \) 连续的

Cubic Hermite Interpolation 的基本形式是：

$$ \boldsymbol{p}(t) = (2t^3-3t^2+1)\boldsymbol{p_k} + (t^3-2t^2+t)\boldsymbol{m_k} + (-2t^3+3t^2)\boldsymbol{p_{k+1}} +(t^3-t^2)\boldsymbol{m_{k+1}} $$

对其简单求导即可获得切线参数方程

$$ \boldsymbol{p^{’}}(t) = (6t^2-6t)\boldsymbol{p_k} + (3t^2-4t+1)\boldsymbol{m_k} + (-6t^2+6t)\boldsymbol{p_{k+1}} +(3t^2-2t)\boldsymbol{m_{k+1}} $$

Catmull–Rom Spline

我们首先介绍最简单的 Cubic Hermite Interpolation, 称为 Catmull–Rom Spline。Catmull–Rom Spline 的特点是只需要给一组点序列的位置，就可以生成穿过所有点的 \( G^1 \) 连续曲线。

Catmull-Rom Spline 生成的管道

Catmull-Rom Spline 生成的管道

这种曲线是由刚刚获得图灵奖的 Catmull 教授发明的，他也是渲染器 Renderman 的作者。 在 Catmull–Rom Spine中，第 \( k \) 个点的切线 \( m_k \) 来自前后两个点的斜率：

$$ m_k = \frac{\boldsymbol{p_{k+1}}-\boldsymbol{p_{k-1}}}{x_{k+1}-x_{k-1}} $$

Catmull–Rom Spline 的计算非常简单，所以在图形学中应用最多，通常用来控制相机的平滑运动路径。

Catmull-Rom 曲线有两个改进形式分别是 Centripetal Catmull–Rom spline 和 Chordal Catmull–Rom Spline。

使用 Threejs 画出这三种形式，你可以看到 Centripetal 和 Chordal 以及原始（Uniform）CR 曲线的区别，前两者更加圆滑一些

Kochannek-Bartels Spline

Catmull–Rom Spline 虽然简单，但也有很明显的限制，比如不能调整样条的弯曲程度。 Kochannek-Bartels Spline 可看做是 Catmull-Rom Spline 的一般化。设计者除了需要给出节点位置，还需给出额外三个参数，节点的切线长度因子 t，切线方向偏移因子 b，以及节点连续性因子 c。通过这三个参数控制每段样条两端节点的切线 \( \mathbf{d_i} \) 和 \( \mathbf{d_{i+1}} \)

$$ \mathbf{d_i} = \frac{(1-t)(1+b)(1+c)}{2}(\mathbf{p_i} - \mathbf{p_{i-1}}) + \frac{(1-t)(1-b)(1-c)}{2}(\mathbf{p_{i+1}}-\mathbf{p_i}) $$

$$ \mathbf{d_{i+1}} = \frac{(1-t)(1+b)(1-c)}{2}(\mathbf{p_{i+1}} - \mathbf{p_i}) + \frac{(1-t)(1-b)(1+c)}{2}(\mathbf{p_{i+2}}-\mathbf{p_{i+1}}) $$

上直观反应了三个参数对曲线形状的影响。可以认为 Kochanek–Bartels Spline 是对 Catmull-Rom Spline 的整体修正。

Non-Uniform Rational B Spline (NURBS)

nurbs 曲面和曲线

NURBS 已经成为计算机建模领域的标准技术之一。在工业领域的很多复杂曲面都是用 NURBS 技术建立的. 例如汽车的外形，螺旋桨等。NURBS 的表达能力很强，能够表达二次曲面（一般的B-Spline就不行）。

在游戏领域里，NURBS 一般应用在生产阶段，用来帮助生成最终的三角形网格游戏资产。目前 NURBS 还比较少见应用于游戏的运行时，但基于程序化生成技术的发展，将来有可能我们能见到从 NURBS 实时生成的游戏模型。

B-Spline

Uniform B-Spline 生成的回形针

Uniform B-Spline 生成的回形针

NURBS 是 B-Spline 的推广。这个 B 指代 Basis，意思是这种样条曲线是根据样条基的线性组合而成的。一个具有 \( n \) 个控制点的 \( p \) 次 B-Spline 的形式为：

$$ P(t) = \sum_{i=0}^{n} P_iN_{i, p}(t) $$

其中 \( P_i \) 是控制点。\( N_{i, p}(t) \) 称为样条基函数，是一个 \( p \) 次（\( p+1 \)阶）多项式。

目前采用最广的样条基函数是由 Cox-De Boor 公式递归定义的。

$$ N_{i,0}(x):= \begin{cases} &1 &,t_{i}\leq x<t_{i+1} \\\ &0 &,{otherwise} \end{cases} $$

$$ N_{i,p}(x):={\frac {x-t_{i}}{t_{i+p}-t_{i}}}N_{i,p-1}(x)+{\frac {t_{i+p+1}-x}{t_{i+p+1}-t_{i+1}}}N_{i+1,p-1}(x). $$

其中，\( (t_0, t_1, … t_m) \) 是一个在 \( [0, 1] \) 上非递减的 m 维实数序列，且 \( m \) 满足 \( m = n + p - 1 \)。节点向量可以视作曲线的参数分割点，也就是说，我们在曲线上可以找到m个点，其位置分别是 \( P(t_0), P(t_1), … P(t_m) \)。

为什么非递减？因为节点可以重复，重复节点可以降低曲线上这一点的连续级别

观察上式可知，\( N_{i, p}(t) \) 的非零区间是 \( [t_i, t_{i+p+1}) \)。其他都是零，正是这个性质保证了 B 样条曲线的局部性。局部性的好处就是我们在实际计算的时候，可以先找到参数 t 所在区间 \( [u_k, u_{k+1}] \), 然后只考虑该区间前面 p+1 个控制点和相应基函数，即能得到等价结果

进一步可以证明，两个节点 \( P_i \), \( P_{i+1} \) 之间的曲线上一点受到 \(p+1\) 个样条基函数的影响，即 \( N_{i-p} … N{i} \)

所以我们可以重写 B-Spline 为：

$$ P(t) = \sum_{i=k-p}^{k} P_iN_{i, p}(t) $$

目前为止，我们可以直接使用上式计算曲线，但也可以选择更加数值稳定的De Boor算法。

De Boor Algorithm

De Boor 算法利用了 B 样条的一个操作：节点插入

节点插入是指在节点向量中插入一个新的分割点，而不改变曲线的形状，由于 m = n + p 的限制，每插入一个节点，控制点的个数也会相应增加一个。由局部性可知，基函数在 u 位置非零的点有 \( p + 1 \)个，所以节点插入影响了\( p + 1 \)个控制点的位置。

很有意思的是，新的控制节点会落在 \( P_{k-p} \) 和 \( P_k \) 之间线段组成的多边形的边上，类似切角的效果

假设我们想在位置 \( t \) 插入新节点 \( t \in [t_k, t_{k+1}) \)，那么受到影响的控制点为 \(P_{k-p + 1}, P_{k-p + 2} … P_k \)，设 \( i \in [k-p+1, k] \)， 新控制点 \( Q_i \) 可以由两个相邻的原控制点线性插值得到

$$ Q_i = (1 - a_i) P_{k-1} + a_k P_k $$

其中

$$ a_i = \frac{t-t_i}{t_{i+p}-ti} $$

现在我们来思考一个问题：如果在同一个位置重复插入节点 \( u, u \in [u_k, u_{k+1}) \)，，直到该节点的重复度达到 \( p + 1 \)。会发生什么？recall Cox-de Boor 公式，在u上非零的基函数最终只有一个，影响该点的控制点只有一个，而系数也只能是1。于是我们得出了一个求曲线上参数为u的一点的方法：**参数为 u 的点的位置就是在该处重复插入节点直到重复度为 \( p+1 \)时得到的新控制点 \( Q_i \) **

这就是 De Boor 算法的主要内容。同一位置连续插入在图形上的表现就是不断对控制点多边形切角，最后到达曲线上一点，这个过程生成的多边形网格叫做 De Boor net。De Boor 算法的具体计算实际也是不断结算切角之后的新控制点。

节点插入算法不仅仅是提供了一个数值稳定的求值方法，实际上也可以按照工程师的需求来精细化曲线和曲面的局部特征。

具体算法可以在这份讲义里查到，这里就不在赘述

B-Spline To Piecewise Bezier

节点插入算法还有什么好处呢？可以证明，如果一个 B 样条的每个节点重复度都是 p + 1，那么这个B样条就是一个分段线性贝塞尔函数，首尾相接的连续4个点就是每段贝塞尔曲线的控制点。如果我们往一个B样条的节点向量上不断插值，直到每个节点都具有 p + 1 个重复度，我们就得到了该 B 样条曲线的等价分段贝塞尔表示。于是针对贝塞尔曲线的各种算法就可以直接应用了，如这篇快速求点到Bezier曲线最近点的文章就可以直接复用过来。这样我们就能计算出 B 样条的 SDF，想象空间就很大了。

What is Non Uniform and Rational?

  • Non Uniform: 这个很好理解，只要节点向量之间不等距，就是 Non Uniform了。Uniform 节点向量可以自动生成，但表达能力差一些。Non Uniform 的话，需要用户输入，但控制能力更强
  • Rational：B 样条已经很强大，但似乎不能表示一些最简单的曲线：思考一下，能用多项式参数方程来表示二维圆弧吗？直觉上，一个二维圆弧的参数方程应该用三角函数定义：\( x = r\cos{t}, y = r\sin{t} \) ；如果把 \( u = \tan{t} \) 代入，则得到 \( x = r\frac{1+u^2}{2u}, y = r\frac{1-u^2}{2u} \) 。这是一个有理函数，也就是说，我们必须构造一个有理函数，才能表示一个圆弧。

NURBS具体是怎么构造有理函数的呢？我们尝试给每个控制点加一个权重 \( w \)，把控制点改写为一个齐次坐标\( P_i = (x_i, y_i, z_i, 1) \), 然后乘系数 \( w_i \) 就得到 \( w_iP_i = (w_ix_i, w_iy_i, w_iz_i, w_i) \)。由齐次坐标的性质可知，这个变换不改变控制点的位置，于是 NURBS 的定义式可以由 B 样条的定义式改写为

$$ P(t) = \frac{\sum_{i=0}^n N_{i, p}(t) w_i P_i}{\sum_{i=0}^n N_{i, p}(t)w_i} $$

这显然是一个有理函数，分子和分母的次数都是 \( p \)，这也是NURBS 中 Rational 的由来

观察分子，我们发现这个分子实际上是一个 4D B-Spline 的表达式，分母只不过把这个 4D B-Spline 投影到了 3D。因此我们可以自然地拓展 De Boor 算法到 NURBS，更进一步，我们可以把 NURBS 转换为 4D Piecewise Bezier Curve，也就是 Rational 3D Bezier Curve。这样我们就有可能求得一点到 NURBS 的最近点，从而定义 NURBS 的 SDF