曲线插值：大直若曲

声明：以下内容出现了一些 TeX 数学公式，它们显示不出来，这不是我的错 :)

更正：段错误网站有个 Bug，它会自动检测文本里是否出现了 TeX 数学公式，然后才能进行公式渲染，但检测规则是必须出现行间公式，即

$$ E = mc^2 $$

这样的公式。

在一些与几何模型设计相关的工作中，设计人员通常会先构造出一些型值点，而这些点往往来自于现实中的一些经验或测量结果。为了美观，或者为了便于加工，或者有些情况下是为了让设计转化为现实之后能够满足现实中的性能需求，例如飞机、汽车、船舶等交通工具都要考虑如何在形体设计上降低流体的阻力，总之，在多数情况下会期望曲线/曲面能够尽量光滑地穿过这些型值点。

这种问题在 CAGD 中，就是曲线插值问题，即空间中存在一组点，$P_1,\cdots,P_n$，如何用一条多项式曲线光滑的穿过它们。这是计算机辅助几何设计（CAGD）领域的最基本的问题之一。

控制点

由于曲线必须通过 $P_0,\cdots,P_n$，即任何一个 $P_i$，$i\in\{0,1,\cdots,n\}$ 都对曲线的总体形状具有决定作用，因此将 $P_i$ 称为控制点，取其可控制曲线形状之意。

控制点的数量称为待插值曲线的阶数，即曲线的次数加 1。

在一维空间里生活是什么感觉

在思考曲线插值这个问题之前，不妨先思考或回顾一下线性插值。或者不妨假设，我们现在也只会解决线性插值问题。

线性插值，就是求穿过点 $P_0$ 与 $P_1$ 的一条直线。不要去回想我们曾经学过的在直角坐标系上定义的直线方程了，原因很简单，直线之所以是直线，并非直角坐标系决定的，而是它的内蕴决定的。

过 $P_0$ 与 $P_1$ 的一条直线，它的内蕴是 $P(t) = \frac{t_1 - t}{t_1 - t_0}P_0 + \frac{t - t_0}{t_1 - t_0}P_1$。为什么这样说呢，因为无论是直线还是曲线，它们构成的是一维空间，这里的 $t$ 就是这个空间中任意一个点的「坐标」。

如果我生活在这条直线之内，那么我只知道我在这个空间中的任何走动只会导致 $t$ 发生变化，除此之外，我一无所知。至于我在这条直线所嵌入的那个空间内，我是 $P_0$ 还是 $P_1$ 或者其他什么东西，那都与我无关。当然，如果你生活与 $P_0$ 与 $P_1$ 同一维度的空间里，当我自认为经过 $t = t_0$ 时，在你看来，我经过的是 $P_0$，而当我自认为经过 $t = t_1$ 时，在你看来则是经过 $P_1$。

我甚至不知道这个空间是不是弯曲的，因为在我看来，它始终都是直的，而且只有两个方向，除非我发现我向一个方向走啊走，有一天发现我又走回了我出发的地方，这时我可以断言我所在的空间是弯曲的。

我不知道我的空间弯曲了

现在，你在你生存的世界里要对不共线的三个点 $P_0$、$P_1$、$P_2$ 进行曲线插值，插值的结果就是让我生存的空间变得弯曲，然而我却对此一无所知。当我在这个空间中走动的时候，当我经过 $t_0$、$t_1$、$t_2$ 时，我一直觉得是在一个平直的空间里行走，而在你看来，我经过的则是 $P_0$、$P_1$、$P_2$ 。

那么，你应当如何给我制造这样的空间呢？很简单，就是做三次线性插值。第一次是对 $P_0$ 与 $P_1$ 进行线性插值，结果为 $P_{01}(t)$。第二次是对 $P_1$ 与 $P_2$ 进行线性插值，结果为 $P_{12}(t)$。第三次是对 $P_{01}(t)$ 与 $P_{12}(t)$ 插值，结果为 $P_{012}(t)$。 结果就是，你在你的世界里构造了一条弯曲的 2 次光滑曲线 $P_{012}(t)$，并让它无比精确地穿过了 $P_0$、$P_1$、$P_2$，然而在我看来，我的空间依然平直。

大直若曲

$P_{01}(t)$ 与 $P_{12}(t)$ 很容易理解，在你的空间里，它是两条相交的直线，$P_1$ 是交点。那么，如何对 $P_{01}(t)$ 与 $P_{12}(t)$ 进行线性插值呢？显然，在 $t_0 \le t \le t_1$ 时，你希望我位于 $P_{01}(t)$ 中，而 $t_1 \le t \le t_2$ 时，你希望我位于 $P_{12}(t)$ 中，能同时满足这些要求的对 $P_{01}(t)$ 与 $P_{12}(t)$ 的线性插值形式只有一种，即 $P_{012}(t) = \frac{t_2 - t}{t_2 - t_0}P_{01}(t) + \frac{t - t_0}{t_2 - t_0}P_{12}(t)$。

这似乎没什么难理解的，在你们的世界里，曾经有个叫欧几里得的老先生说过，过相异两点，能作且只能作一条直线。对于 $P_{012}(t)$ 而言，按照这位老先生说的去理解就没错。$P_{01}(t)$ 与 $P_{12}(t)$ 虽然是直线，但是只要固定 $t$，它们就得老老实实地「坍缩」为两个点，而此时 $P_{012}(t)$ 是过这两个点的直线上的一点。若不固定 $t$，而是让它递增或递减的发生变化，那么 $P_{012}(t)$ 就可以形成一条曲线。因此，我们也可以像欧老先生那样说话，过相交的两条直线，通过线性插值能作且只能作一条二次曲线。

还可以说得更广泛，过两条相异且存在部分区域重叠的 n 次曲线，通过线性插值，能作且只能作一条 n + 1 次曲线。例如，若沿袭上文对 $P(t)$ 的标记方式，那么对 $P_{012}(t)$ 与 $P_{123}(t)$ 这两条相异且存在部分区域重叠的二次曲线进行线性插值，可以得到 $P_{0123}(t)$，它是一条三次曲线，过 $P_0$、$P_1$、$P_2$、$P_3$ 四个点。

像 $P_{01\cdots n}(t)$ 这样的函数，在中学数学课本里称为多项式。现在，它的系数包含着 n + 1 个点，它就被称为多项式曲线了。数学家们已经证明了，过 n + 1 个点的 n 次多项式曲线不仅存在，而且惟一。正是基于这一结论，在上文我才敢胡说八道，开欧几里得老先生的玩笑。

你看，就一直在玩线性插值的戏法，结果能变出来在你看来挺弯曲的线，而对于活在一维空里的我看来，它还是一条很直的线……我们这个空间，也有个老先生，他实在是太老了，所以叫老子，他说，大直若曲，大巧若拙。此言真是不虚啊。后来，还有个叫姜太公的老先生，用直的鱼钩钓了条大鱼。

Neville 算法

有位老先生，名叫 Eric Harold Neville，是一个数学家。他提出一个算法，很简洁地描述出上述的曲线插值过程。这个算法以他的姓氏命名，叫 Neville 算法。这位老先生，我国的人知道的似乎并不太多。在百度上搜索「Neville 算法」，结果连前 3 页都未能占满。

这个算法很简单，用图的形式描述起来会很直观，不过它却是计算机算法类教科书说不清、道不白的大名鼎鼎的动态规划算法。以 2 次曲线插值三个点为例，下面便是 Neville 算法的动态规划图：

我在 Hamal 项目中提供了这个算法的 Python 实现，见 https://github.com/liyanrui/h...

这个脚本的用法如下：

$ python3 neville --curve="foo-curve" --sample=100 --nodes=foo-nodes.asc foo.asc

foo.asc 为控制点集文件，其中每一行是控制点的坐标，例如：

1 1 0
3 2 0
4 3 0
7 2 0

foo-nodes.asc 为控制点对应的 t 值，例如：

0
1
2
3

neville 脚本会根据 --sample 选项的值作为采样次数对插值曲线进行均匀采样。例如，当 --sample 指定的值为 10 时，neville 脚本会以 0, 1/10, 2/10, ..., 10/10 对插值曲线进行求值，从而得到 10 个点，这些点被保存在 foo-curve.asc 文件，而它们的连接关系被保存在 foo-curve-graph.asc 文件。这两份文件的名字，前缀部分皆为 foo-curve，由 --curve 选项给定。

可视化

可以使用 hamal 脚本将插值曲线的采样点集、采样点之间的关系（形成一组首尾相连的线段集合的关系）以及控制点等数据转换为 POV-Ray 场景文件，交由 povray 进行图形渲染。

对于上一节 neville 脚本示例的输入数据及输出结果，使用 hamal 脚本转换为 POV-Ray 场景文件及图形渲染的命令如下：

$ python3 hamal --object="neville" \
    --curve=foo-curve-graph.asc \
    --control-points=foo.asc --line-width=0.015 \
    --control-point-size=0.04 foo-curve.asc
$ povray +P neville.pov

结果如下图所示：

关于 hamal 脚本的其他用法，见「hamal 指南」。

下一篇：曲线里的时间、位移与质量