CAGD 学习笔记 | 样条

2023-01-30 约 3481 字预计阅读 7 分钟

在上一篇笔记中, 我们介绍了 Bézier 曲线, 它是一条由若干控制点给出的光滑的多项式曲线, 而这些控制点对应的基函数是 Bernstein 基函数. 不难注意到, Bernstein 基函数 $B^n_k(t)$ 在 $\frac{k}{n}$ 处达到极大值, 在远离 $\frac{k}{n}$ 时不断减小直到 $0$. 因此每一个控制点在它附近对 Bézier 曲线的形状影响最大, 而在远处影响较小但始终非 $0$. 用户移动某一控制点时, 总会对曲线整体造成影响. 此外, 当控制点增多时, 会产生更高阶的多项式, 影响计算效率.

在实践中, 我们有时会希望控制点只在它的附近对曲线的形状有影响, 并且避免高阶多项式, 此时 Bézier 曲线将无法达到要求. 为此我们需要使用另一种曲线: 样条曲线 (Spline curve).

样条曲线指的是平滑的分段多项式曲线, 而这个名字来自于二战期间英国制造飞机时的技术, 使用若干段薄木条穿过固定的点来形成飞机上的曲线. 因为样条曲线是分段的, 所以可以使每一个控制点的影响范围在它的附近. 此外, 样条曲线的另一个优势在于可以控制曲线的阶数, 避免高次多项式的出现.

本文将介绍 Bézier 样条曲线和 B 样条曲线两种样条曲线, 它们的区别在于构造方法的不同, 但本质都是分段的多项式曲线.

Bézier 样条曲线

Bézier 样条曲线 (Bézier spline curve) 的想法很简单, 就是给定一系列的控制点, 构造一系列相连的 Bézier 曲线将它们相连, 并要求它们在相连点处有很好的连续性.

问题背景

对于空间中的 $n+1$ 个控制点 $\{p_0, p_1, \dots, p_n\}$, 和若干参数结点 $\{0 = t_0 < t_1 < \dots < t_n = 1\}$, 我们希望有一段参数曲线 $p(t)$, 满足

$$ p(t_k) = p_k,\ k = 0, \dots, n \tag{1} $$

对于每一个 $0 \leq k \leq n-1$ 在 $t_k < t < t_{k+1}$ 上 $p(t)$ 是一段 Bézier 曲线, 并且曲线 $p(t)$ 有较好的连续性 (有限阶可导).

计算一条 Bézier 样条曲线, 要做的就是根据连续性条件确定每一段 Bézier 曲线中控制点的位置. 其中每一段曲线的首尾两个控制点就是 $p_k, p_{k+1}$.

3 阶样条曲线

3 阶的 Bézier 样条曲线指的是, 每一段 Bézier 曲线都是 3 阶的, 有 4 个控制点的曲线, 它使得曲线直到 2 阶的导数都是连续的. 这也是实践中常用的样条曲线, 被称为 3 阶样条曲线 (Cubic spline curve). 我们希望

对于 $t_k < t < t_{k+1}$ 上的一段 Bézier 曲线, 它的控制点记为 $p_k, p_k^+, p_{k+1}^-, p_{k+1}$. 通过 $[t_k, t_{k+1}]$ 与 $[0, 1]$ 之间的参数变换, 可以计算出

$$ \begin{aligned} \lim_{t \rightarrow t_k +0} p' (t) & = \frac{m}{t_{k+1}-t_k} \left(p_k^+ - p_k\right) \\ \lim_{t \rightarrow t_{k+1} -0} p' (t) & = \frac{m}{t_{k+1}-t_k} \left(p_{k+1} - p_{k+1}^-\right) \\ \lim_{t \rightarrow t_k +0} p'' (t) & = \frac{m(m-1)}{(t_{k+1}-t_k)^2} \left(p_k - 2 p_k^+ + p_{k+1}^-\right) \\ \lim_{t \rightarrow t_{k+1} -0} p' (t) & = \frac{m(m-1)}{(t_{k+1}-t_k)^2} \left(p_k^+ - 2 p_{k+1}^- + p_{k+1}\right) \end{aligned} $$

因此曲线在 $t_k$ 处的高阶连续性条件为

$$ \begin{aligned} \frac{p_k - p_k^-}{t_{k}-t_{k-1}} & = \frac{p_k^+ - p_k}{t_{k+1}-t_k} \\ \frac{p_{k-1}^+ - 2 p_{k}^- + p_{k}}{(t_{k}-t_{k-1})^2} & = \frac{p_k - 2 p_k^+ + p_{k+1}^-}{(t_{k+1}-t_k)^2} \end{aligned} \tag{2} $$

此时我们有 $2n$ 个未确定的控制点和 $2(n-1)$ 个条件, 需要额外再增加 2 个条件以确定唯一的一条曲线, 以下为两种附加条件:

自然端点条件: $p''(0) = p''(1) = \mathbf{0}$;
Bessel 端点条件: $p'(0) = q'(0),\ p'(1) = r'(1)$, 其中 $q$ 插值前三个点 $(t_0, p_0), (t_1, p_1), (t_2, p_2)$ 的 2 次多项式曲线, $r$ 是插值后三个点的 2 次多项式曲线.

如下图所示, 不同的附加条件将导致曲线有不同的形状.

如下图所示的是采用自然端点条件绘制的 3 阶自然样条曲线. 其中蓝色折线为插值点, 空心黑色折线为各段 Bézier 曲线的控制多边形, 后者的位置是根据插值点的位置计算出来的.

观察如上的连续性条件可知, 可以发现 $p_k$ 并非只会影响与它相邻区间中控制点的位置, 因此移动它也不止会改变与它相邻的两段 Bézier 曲线的形状. 事实上, 在实践中我们有时可以舍弃二阶连续性的条件 (例如 Photoshop 软件中的钢笔工具) 以达到使各个插值点对曲线的影响限制在局部的效果, 这时将把各段 Bézier 曲线的控制点也开放给用户修改.

参数化

通常情况下, 绘制 Bézier 样条曲线时只有给定若干控制点的条件, 而参数结点 $\{0 = t_0 < t_1 < \dots < t_n = 1\}$ 的值需要另外确定, Bézier 样条曲线的参数化指的就是确定这一组数的值. 值得注意的是, 通常情况下的参数化是给定某一个曲线或曲面, 然后构造一个它与参数域之间的一一映射, 但这里的参数化则不同, 因为如式 (2) 所示, 所绘制的 Bézier 样条曲线的形状是与参数结点的选取相关的. 以下为几种参数化方法:

等距参数化 (Equidistant (uniform) parameterization): $t_{k+1} - t_k = \text{const}$;
向心参数化 (Centripetal parameterization) ¹ : $t_{k+1} - t_k = \sqrt{\|p_{k+1} - p_k\|}$.

如下图所示, 不同的参数化方法将导致曲线有不同的形状.

B 样条曲线

B 样条曲线 (B-spline curve) 的 “B” 代表 basis (基), 它的形式与 Bézier 曲线类似, 也是用控制点对基函数作线性组合得到的. B 样条曲线与 Bézier 曲线的不同之处在于它的基函数具有局部支撑性, 并且曲线的阶数与控制点的个数无关.

B 样条基函数

设 $(t_0, t_1, \dots, t_m)$ 为结点向量, $t_0 \leq t_1 \leq \dots \leq t_m$. 如下的 Cox-de Boor 递归公式 给出了 B 样条基函数 (B-spline basis function) 的定义

$$ \begin{aligned} N_{k, 0}(t) & = \begin{cases} 1, & t_k \leq t1 \text{ and } k=0, \dots, n \end{aligned} $$

$N_{k, d}$ 称为 $d$ 阶的第 $k$ 个 B 样条基函数. 不难发现, $N_{k,d}$ 在使它非负的每一个区间 $[t_{l}, t_{l+1}]$ 上都是 $d$ 次的多项式.

例子

根据定义, 当结点向量中不出现重复时, 计算 B 样条基函数显式表达式的 Mathematica 代码如下:

BsBF[t_, k_, d_, knotVector_] :=
  If[d == 0, If[t >= knotVector[[k]] && t < knotVector[[k + 1]], 1, 0],
    (t - knotVector[[k]])/(knotVector[[k + d - 1]] -
      knotVector[[k]]) BsBF[t, k, d - 1, knotVector] +
    (knotVector[[k + d]] - t)/(knotVector[[k + d]] -
      knotVector[[k + 1]]) BsBF[t, k + 1, d - 1, knotVector]]

取结点向量为 $t_i = i$, 计算 2 阶的第 2 个 B 样条基函数得到显式表达式为

$$ N_{2,2}(t) = \begin{cases} \frac{1}{2}t^2- t+\frac{1}{2} & 1\leq t<2; \\ -t^2+5 t-\frac{11}{2} & 2\leq t<3; \\ \frac{1}{2} t^2-4 t+8 & 3\leq t<4; \\ 0 & \text{otherwise}. \end{cases} $$

其图像为

$N_{2,2}(t) 的图像$

这一函数具有局部支撑性. 并且容易验证, 它的 1 阶导数连续, 而 2 阶导数不连续.

B 样条基函数的性质

B 样条基函数具有如下性质

$N_{k,d}(t) \geq 0$;
$d>0$ 时, $N_{k,d}(t)$ 的支撑集为 $(t_k, t_{k+1+d})$;
在任意一个区间 $[t_k, t_{k+1})$ 上至多只有 $d+1$ 个 $d$ 阶基函数非负;
$n+1$ 个 $d$ 阶的 B 样条基函数需要 $n+d+2$ 个结点;
单位分解性 (Partition of Unity): 若某一区间上恰有 $d+1$ 个 $d$ 阶基函数非负, 那么这全部的 $d+1$ 个非负的基函数之和为 1;
若一个结点重复 $p$ 次, 则基函数在这一点处是 $C^{d-p}$.

下图为定义在结点向量 $(0, 1, \dots, 10)$ 上的 7 个 3 阶 B 样条基函数及它们的和. 拥有 $d+1 = 4$ 个非负基函数的区间共有 4 个, 即 $[3, 4),\ \dots,\ [6,7)$, 并且它们的和恰好在这些区间上为 1.

我们知道, Bézier 曲线的 Bernstein 基函数满足单位分解性, 这使得曲线上的每个点都是控制点的凸组合. 同样地, 我们也希望 B 样条曲线也具备同样的性质. 因此对于有 $n+1$ 个 $d$ 阶的 B 样条基函数的曲线, 它的结点向量为 $(t_0, t_1, \dots, t_{n+d+1})$, 并且在区间 $[t_d, t_{n+1})$ 上这些基函数满足单位分解性, 取这一区间为曲线的参数域.

B 样条曲线的性质

根据 B 样条基函数的性质, 我们可以得到一些 B 样条曲线的性质.

中间结点只出现 1 次的 $d$ 阶 B 样条曲线是 $C^{d-1}$ 的;
控制点 $p_k$ 只影响曲线上 $t \in [t_k, t_{k+1+d}]$ 的部分;
凸包性, 仿射无关性.

de Boor 算法

在本节中, 我们将介绍如何用 de Boor 算法绘制 B 样条曲线.

根据 B 样条基函数的递归定义, 我们将直接推导出 de Boor 算法.

$$ p(t) = \sum_{k=0}^{n} N_{k, d} (t) p_{k} $$

其中控制点 $p_k$ 也被称为 de Boor 点.

$$ \begin{aligned} p(t) & = \sum_{k=0}^{n} N_{k, d} (t) p_{k} = \sum_{k=l-d}^{l} N_{k, d} (t) p_{k} \\ & = \sum_{k=l-d}^{l} \left(\frac{t-t_k}{t_{k+d}-t_k} N_{k, d-1}(t)+\frac{t_{k+d+1}-t}{t_{k+d+1}-t_{k+1}} N_{k+1, d-1}(t)\right) p_k \\ & = \sum_{k=l-d}^{l} \frac{t-t_k}{t_{k+d}-t_k} N_{k, d-1}(t) p_k + \sum_{k=l-d+1}^{l+1} \frac{t_{k+d}-t}{t_{k+d}-t_{k}} N_{k, d-1}(t) p_{k-1} \\ & = \sum_{k=l-d+1}^{l} \frac{t-t_k}{t_{k+d}-t_k} N_{k, d-1}(t) p_k + \sum_{k=l-d+1}^{l} \frac{t_{k+d}-t}{t_{k+d}-t_{k}} N_{k, d-1}(t) p_{k-1} \\ & = \sum_{k=l-d+1}^{l} \left( \frac{t-t_k}{t_{k+d}-t_k} p_k + \frac{t_{k+d}-t}{t_{k+d}-t_{k}} p_{k-1} \right) N_{k, d-1}(t) \end{aligned} $$

现令

$$ p_k^{(r)}(t) = \begin{cases} (1 - \alpha_k^{(r-1)} (t)) p_{k-1}^{(r-1)} (t) + \alpha_k^{(r-1)}(t) p_{k}^{(r-1)}(t) & \text{if } 1 \leq r \leq d \\ p_k & \text{if } r = 0 \end{cases} $$$$ \alpha_k^r = \frac{t-t_k}{t_{k+d-r}-t_k} $$$$ \begin{aligned} p(t) & = \sum_{k=l-d+1}^{l} \left( \frac{t_{k+d-1}-t}{t_{k+d-1}-t_{k}} p_{k-1} + \frac{t-t_k}{t_{k+d-1}-t_k} p_k \right) N_{k, d-1}(t) \\ & = \sum_{k=l-d+1}^{l} p_k^{(1)}(t) N_{k, d-1}(t) = \sum_{k=l-d+r}^{l} p_k^{(r)}(t) N_{k, d-r}(t) \\ & = p_l^{(d)}(t) \end{aligned} $$

根据上述讨论, de Boor 算法计算 B 样条曲线的算法如下:

find l satisfying t_l <= t < t_{l+1}
for k = l - d, ..., l
	q[k - l + d] = p[k]
for r = 1, ..., d
	for k = l, ..., l - d + r
		a = (t - t[k]) / (t[k + d - r] - t[k])
		q[k - d + r] = (1 - a) * q[k - 1 - d + r] + a * q[k - d + r]
output: q[d]

结点向量

和 Bézier 样条曲线类似地, B 样条曲线的形状也与结点向量的选取有关. 若某一 $d$ 阶的 B 样条曲线有 $n+1$ 个控制点, 则它有 $n+d+2$ 个参数结点, 以下是几种结点向量的选取方法:

均匀 B 样条: $t_{k+1} - t_k = \text{const}$;
准均匀 B 样条: 设 $s_k = \frac{k}{m}$, 节点向量为 $(t_0, t_1, \dots, t_{n+d+1}) = (\overbrace{s_0, s_0, \dots, s_0}^{\text{重复 } d+1 \text{ 次}}, s_1, \dots, s_{m-1}, \overbrace{s_m, s_m, \dots, s_m}^{\text{重复 } d+1 \text{ 次}})$, 比较等式两端分量个数可知 $m = n-d+1$.
分段 Bézier 曲线: 首尾结点重复度²为 $d+1$, 其它结点重复度为 $d$, 则 B 样条基函数退化为若干组 Bernstein 基函数;
非均匀 B 样条基: 最一般的情况, 只要求首尾结点重复度 $\leq d+1$, 其它结点重复度 $\leq d$.

下图依次为均匀 B 样条基函数, 准均匀 B 样条基函数, 分段 Bernstein 基函数的图像.

$$ p(t_0) = p_0,\ p(t_{n+d+1}) = p_n $$

插值

$$ p(s_k) = p_k,\ 0 \leq k \leq m $$$$ \begin{aligned} s_0 & = t_0 = t_1 = t_2 = t_3 \\ s_k & = t_{k+3},\ 0 \leq k \leq m \\ s_m & = t_{m+3} = t_{m+4} = t_{m+5} = t_{m+6} \end{aligned} $$$$ p(t) = \sum_{l = 0}^{m+2} N^3_{l}(t) d_l $$$$ p(s_k) = p(t_{k+3}) = \sum_{l = 0}^{m+2} N^3_{l}(t_{k+3}) d_l = \sum_{l = k}^{k+2} N_l^3(t_{k+3}) d_l = p_k $$

此时共有 $m+3$ 个变量和 $m+1$ 个条件. 添加自然端点条件 $p''(s_0) = p'' (s_m) = 0$ 即可求解出所需要的 3 阶准均匀 B 样条曲线.

如下图所示的是采用自然端点条件绘制的 3 阶准均匀 B 样条曲线. 其中蓝色折线为插值点, 空心黑色折线为 de Boor 点组成的控制多边形, 后者的位置是根据插值点的位置计算出来的.

B 样条曲线与 Bézier 样条曲线

容易发现上一节中 3 阶准均匀自然 B 样条是二阶连续的、首尾二阶导为 0 的分段 3 阶多项式曲线. 那么它和本文中前半部分中介绍的 3 阶自然样条之间有什么关系?

TODO: 3 阶准均匀自然 B 样条和 3 阶自然样条是一个东西吗？

这一概念的中文是笔者自行翻译的, 以英文名称为准. ↩︎
结点向量中某一结点的重复度指出现的次数, 出现 1 次的结点重复度为 1. ↩︎

目录