计算机辅助几何设计(2)：Bézier样条曲线

1. 简介

之前我们已经介绍过Bézier曲线的原理与实现方法，但由于Bézier曲线的一些问题，例如：

高多项式次数
- $n+1$个Bézier控制点将生成一条$n$次多项式表达的Bézier曲线
全局，伪局部性
- Bézier曲线在始末两控制点进行插值，在中间控制点进行拟合，当控制点增加时，曲线形态难以控制

因此我们引入的Bézier样条曲线来解决上述两个问题。在了解Bézier样条曲线的工作原理之前，需要弄清楚几何设计中连续性的基本概念。

1.1. 连续性的基本概念

1.1.1. 参数连续性

给定两条曲线：
$$
\begin{align}
&\pmb x_1(t)\ \mathrm{over}\ [t_0,t_1]\\\\
&\pmb x_2(t)\ \mathrm{over}\ [t_1,t_2]
\end{align}
$$
若$\pmb x_1$和$\pmb x_2$在$t_1$处的0阶到$r$阶导数向量均重合，则若$\pmb x_1$和$\pmb x_2$在$t_1$处$C^r$连续

常见的参数连续性有：

$C^0$：位置变化连续
$C^1$：一阶导数在交界处连续（速度向量相同）
$C^2$：二阶导数在交界处连续（加速度向量相同）

1.1.2. 几何连续性

给定两条曲线：
$$
\begin{align}
&\pmb x_1(t)\ \mathrm{over}\ [t_0,t_1]\\\\
&\pmb x_2(t)\ \mathrm{over}\ [t_1,t_2]
\end{align}
$$
若$\pmb x_1$和$\pmb x_2$能够以某种方式重新参数化使得在$t_1$处$C^r$连续，则$\pmb x_1$和$\pmb x_2$在$t_1$处$G^r$连续

常见的几何连续性有：

$G^0=C^0$：位置变化连续性（连接性）
$G^1$：切线方向变化连续性（相同切线）
- 正则化切线变化连续
- 等价于曲线能够重新参数化到$C^1$
- 等价于单位速度参数化为$C^1$
$G^2$：曲率变化连续性（相同切线与曲率）
- 等价于曲线能够重新参数化到$C^2$
- 等价于单位速度参数化为$C^2$

1.1.3. 参数连续性 vs. 几何连续性

参数连续性	几何连续性
在该曲线上运动的粒子是否有光滑的轨迹？（位置、速度、加速度）	曲线本身是否光滑
取决于参数化方式	独立于参数化方式
应用：动画（物体移动、摄像头轨迹）	与建模相关（曲线设计）

1.2. Bézier样条曲线的设计思想

既然Bézier曲线全局性太强、阶次太高，那么优化思路就是降低全局性、降低阶次，一个有效的方法就是进行分段，把一条长的曲线分段成若干条低阶的Bézier曲线，在衔接处满足某种连接性，即可使得曲线满足控制点插值、低次性和局部性等优点，这就是Bézier样条曲线的设计思路。

设计Bézier样条曲线需要考虑的两个问题是：

参数化
连续性

2. Bézier样条参数化

在Bézier曲线的构造中，我们选择的参数为$t\in[0,1]$，贯穿整条曲线，称之为全局参数。在Bézier样条曲线中，我们通过分段，将曲线分为$[t_0,t_1]$，$[t_1,t_2]$，……，$[t_{n-1},t_n]$的曲线段，每一段对应一条低阶Bézier曲线，这条低阶Bézier曲线使用的参数$u\in[0,1]$则称为局部参数。当局部参数$u$从0变化到1时，全局参数从$t_i$变化到$t_{i+1}$。

通常情况下，节序列$t_0,t_1,\cdots,t_n$可以自定义，选取不同的节序列会得到不同的曲线形态，常用的参数化方法有以下几种：

Equidistant (Uniform) Parameterization
- $t_{i+1}-t_i=\mathrm{const}$
- 不考虑数据点之间的几何关系
Chordal Parameterization
- $t_{i+1}-t_i=\Vert \pmb k_{i+1}-\pmb k_i\Vert$
- 参数间隔与相邻控制点的距离成正比
Centripetal Parameterization
- $t_{i+1}-t_i=\sqrt{\Vert\pmb k_{i+1}-\pmb k_i\Vert}$
- 参数间隔与相邻控制点的距离的开方成正比
Foley Parameterization
- 涉及控制多边形的角度
$$
t_{i+1}-t_i=\Vert\pmb k_{i+1}-\pmb k_i\Vert\cdot\Big(1+\dfrac{3}{2}\dfrac{\hat\alpha_i\Vert\pmb k_i-\pmb k_{i-1}\Vert}{\Vert\pmb k_i-\pmb k_{i-1}\Vert+\Vert\pmb k_{i+1}-\pmb k_i\Vert}+\dfrac{3}{2}\dfrac{\hat\alpha_{i+1}\Vert \pmb k_{i+1}-\pmb k_i\Vert}{\Vert\pmb k_{i+1}-\pmb k_i\Vert+\Vert\pmb k_{i+2}-\pmb k_{i+1}\Vert} \Big)
$$
- 其中，$\hat\alpha_i=\min\Big(\pi-\alpha_i,\dfrac{\pi}{2}\Big)$，且$\alpha_i=\mathrm{angle}(\pmb k_{i-1},\pmb k_i,\pmb k_{i+1})$

3. Bézier样条的连续性

给定

$\pmb y(u)$：$[0,1]$上的Bézier曲线
$\pmb x(u)$：$[t_i,t_{i+1}]$上的Bézier曲线

令$u(t)=\frac{t-t_i}{t_{i+1}-t_i}$，则$\pmb x(t)=\pmb y(u(t))$

对$\pmb x(t)$进行求导有：
$$
\begin{align}
\dot{\pmb x}(t)&=\dot{\pmb y}(u(t))\cdot\dot u(t)=\dfrac{\dot{\pmb y}(u(t))}{t_{i+1}-t_i}\\\\
\ddot{\pmb x}(t)&=\ddot{\pmb y}(u(t))\cdot(\dot u(t))^2+\dot{\pmb y}(u(t))\cdot\ddot u(t)=\dfrac{\ddot{\pmb y}(u(t))}{(t_{i+1}-t_i)^2}\\
\cdots\\\\
\pmb x^{[n]}(t)&=\dfrac{\pmb y^{[n]}(u(t))}{(t_{i+1}-t_i)^n}
\end{align}
$$
由Bézier曲线解析式：
$$
\pmb f(t)=\sum_{i=0}^nB_i^n(t)\pmb b_i
$$
其端点值：
$$
\begin{align}
\pmb f(0)&=\pmb b_0\\\\
\pmb f(1)&=\pmb b_1
\end{align}
$$
端点一阶导数：
$$
\begin{align}
\pmb f’(0)&=n[\pmb b_1-\pmb b_0]\\\\
\pmb f’(1)&=n[\pmb b_n-\pmb b_{n-1}]
\end{align}
$$
端点二阶导数：
$$
\begin{align}
\pmb f’’(0)&=n(n-1)[\pmb b_2-2\pmb b_1+\pmb b_0]\\\\
\pmb f’’(1)&=n(n-1)[\pmb b_n-2\pmb b_{n-1}+\pmb b_{n-2}]
\end{align}
$$
代入到$\pmb x(t)=\pmb y(u(t))$中，有：
$$
\begin{align}
\pmb x(t_i)&=\pmb b_0\\\\
\pmb x(t_{i+1})&=\pmb b_1\\\\
\dot{\pmb x}(t_i)&=\dfrac{n(\pmb b_1-\pmb b_0)}{t_{i+1}-t_i}\\\\
\dot{\pmb x}(t_{i+1})&=\dfrac{n(\pmb b_n-\pmb b_{n-1})}{t_{i+1}-t_i}\\\\
\ddot{\pmb x}(t_i)&=\dfrac{n(n-1)(\pmb b_2-2\pmb b_1+\pmb b_0)}{(t_{i+1}-t_i)^2}\\\\
\ddot{\pmb x}(t_{i+1})&=\dfrac{n(n-1)(\pmb b_n-2\pmb b_{n-1}+\pmb b_{n-2})}{(t_{i+1}-t_i)^2}\
\end{align}
$$

现有$m+1$段$n$阶Bézier样条$\pmb x^{(0)},\pmb x^{(1)},\cdots,\pmb x^{(m)}$构成一条完整曲线$\pmb x$，Bézier样条$\pmb x^{(i)}$的控制点为：$\pmb b_0^{(i)},\pmb b_1^{(i)},\cdots,\pmb b_n^{(i)}$，对应全局参数范围$[t_i,t_{i+1}]$

$C^0$连续性

每个样条线段插值两端控制点，因此相邻曲线段两端的控制点必须重合以获得$C^0$连续性，即：
$$
\pmb b^{(i-1)}_{n}=\pmb b^{(i)}_0
$$
$C^1$连续性

$C^1$连续即该参数化下一阶导数相同，即：
$$
\dot{\pmb x}^{(i)}(t_{i+1})=\dot{\pmb x}^{(i+1)}(t_{i+1})
$$
代入表达式有：
$$
\frac{\pmb b^{(i)}n-\pmb b{n-1}^{(i)}}{t_{i+1}-t_i}=
\frac{\pmb b^{(i)}{1}-\pmb b{0}^{(i)}}{t_{i+2}-t_{i+1}}
$$

$C^2$连续性

$C^1$连续即该参数化下二阶导数相同，即：
$$
\ddot{\pmb x}^{(i)}(t_{i+1})=\ddot{\pmb x}^{(i+1)}(t_{i+1})
$$
代入表达式有：
$$
\frac{\pmb b_n^{(i)}-2\pmb b_{n-1}^{(i)}+\pmb b_{n-2}^{(i)}}{(t_{i+1}-t_i)^2}=
\frac{\pmb b_2^{(i+1)}-2\pmb b_1^{(i+1)}+\pmb b_0^{(i+1)}}{(t_{i+2}-t_{i+1})^2}
$$
令
$$
\pmb d^-=\pmb b_{n-1}^{(i)}+\frac{\Delta_{i+1}}{\Delta_i}(\pmb b_{n-1}^{(i)}-\pmb b_{n-2}^{(i)})
$$
和
$$
\pmb d^+=\pmb b_{1}^{(i+1)}-\frac{\Delta_{i}}{\Delta_{i+1}}(\pmb b_{2}^{(i+1)}-\pmb b_{1}^{(i+1)})
$$
则有$C^2$连续性等价于$C^1$连续性且$\pmb d^-=\pmb d^+$

$G^1$连续性

满足在$t=t_i$处$G^1$连续的条件是：
$$
\begin{aligned}
&\pmb x^{(i)}(t_i)=\pmb x^{(i+1)}(t_i)\\\\
&\dot{\pmb x}^{(i)}(t_i)=\dot{\pmb x}^{(i+1)}(t_i)\\\\
&\ddot{\pmb x}^{(i+1)}(t_i)-\ddot{\pmb x}^{(i)}(t_i)\parallel \dot{\pmb x}(t_i)
\end{aligned}
$$

$G^2$连续性

满足在$t=t_i$处$G^2$连续的条件是：

$G^1$连续
$\pmb b_{n-2}^{(i)},\pmb b_{n-1}^{(i)},\pmb b_{n}^{(i)}=\pmb b_{0}^{(i+1)},\pmb b_{1}^{(i+1)},\pmb b_{2}^{(i+1)}$五个向量共面
且面积满足：
$$
\dfrac{\mathrm{area}(\pmb b_{n-2}^{(i)},\pmb b_{n-1}^{(i)},\pmb b_{n}^{(i)})}{\mathrm {area}(\pmb b_{0}^{(i+1)},\pmb b_{1}^{(i+1)},\pmb b_{2}^{(i+1)})}=\dfrac{a^3}{b^3}
$$

4. 三次Bézier样条曲线

理论上，只要有足够的边界条件约束，我们可以选用任意阶的Bézier样条进行构造我们的曲线，但一般常用的Bézier样条不会超过三次（某些CAD/CAM软件可能会使用）：

零次：分段常数，不光滑
一次：分段线性，不够光滑
二次：二阶导数为常数，不够灵活
三次：图形学应用中最为常用

在IlumEngine中，也选用了三次Bézier样条曲线进行实现。

对于给定的控制点$\pmb k_0,\pmb k_1,\cdots,\pmb k_n$，以及参数化节序列$t_0,t_1,\cdots,t_n$，我们需要生成用于构造连续分段三次Bezier样条曲线的Bézier点$\pmb b_0,\cdots,\pmb b_{3n}$

该问题中，需要求解：

$3n+1$个未知点
$C^0$连续：$\pmb b_{3i}=\pmb k_i$，$i=0,\cdots,n$，共$n+1$个方程
$C^1$连续：$\dfrac{\pmb b_{3i}-\pmb b_{3i-1}}{t_i-t_{i-1}}=\dfrac{\pmb b_{3i+1}-\pmb b_{3i}}{t_{i+1}-t_{i}}$，$i=0,\cdots,n$，共$n-1$个方程
$C^2$连续：$\dfrac{\pmb b_{3i}-2\pmb b_{3i-1}+\pmb b_{3i-2}}{(t_i-t_{i-1})^2}=\dfrac{\pmb b_{3i+2}-2\pmb b_{3i+1}+\pmb b_{3i}}{(t_{i+1}-t_{i})^2}$，$i=0,\cdots,n$，共$n-1$个方程
两个端点条件

端点条件的选取，通常有两种：

Bessel End Condition
- $\pmb k_0$处的切向量等于插值${\pmb k_0,\pmb k_1,\pmb k_2}$的抛物线在$\pmb k_0$处的切向量
- 插值${\pmb k_0,\pmb k_1,\pmb k_2}$的抛物线方程：
  $$
  \pmb p(t)=
  \dfrac{(t_2-t)(t_1-t)}{(t_2-t_0)(t_1-t_0)}\pmb k_0+
  \dfrac{(t_2-t)(t-t_0)}{(t_2-t_1)(t_1-t_0)}\pmb k_1+
  \dfrac{(t_0-t)(t_1-t)}{(t_2-t_1)(t_2-t_0)}\pmb k_2
  $$
- 插值抛物线导数：
  $$
  \dot{\pmb p}(t_0)=
  -\dfrac{(t_2-t_0)+(t_1-t_0)}{(t_2-t_0)(t_1-t_0)}\pmb k_0+
  \dfrac{(t_2-t_0)}{(t_2-t_1)(t_1-t_0)}\pmb k_1-
  \dfrac{(t_1-t_0)}{(t_2-t_1)(t_2-t_0)}\pmb k_2
  $$
- $\pmb b_1$的位置：
  $$
  \pmb b_1=\pmb b_0+\dfrac{t_1-t_0}{3}\dot{\pmb p}(t_0)
  $$
Natural End Condition
$$
\begin{align}
\ddot {\pmb x}(t_0)&=0\Leftrightarrow \pmb b_1=\dfrac{\pmb b_2+\pmb b_0}{2}\\\\
\ddot{\pmb x}(t_n)&=0\Leftrightarrow \pmb b_{3n-1}=\dfrac{\pmb b_{3n-2}+\pmb b_{3n}}{2}
\end{align}
$$