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 InterpolationNon 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