Bezier curve

贝塞尔曲线的数学基础是早在 1912 年就广为人知的伯恩斯坦多项式。但直到 1959 年，当时就职于雪铁龙的法国数学家 Paul de Casteljau 才开始对它进行图形化应用的尝试，并提出了一种数值稳定的 de Casteljau 算法。然而贝塞尔曲线的得名，却是由于 1962 年另一位就职于雷诺的法国工程师 Pierre Bézier 的广泛宣传。他使用这种只需要很少的控制点就能够生成复杂平滑曲线的方法，来辅助汽车车体的工业设计。

正是因为控制简便却具有极强的描述能力，贝塞尔曲线在工业设计领域迅速得到了广泛的应用。不仅如此，在计算机图形学领域，尤其是矢量图形学，贝塞尔曲线也占有重要的地位。今天我们最常见的一些矢量绘图软件，如 Flash、Illustrator、CorelDraw 等，无一例外都提供了绘制贝塞尔曲线的功能。甚至像 Photoshop 这样的位图编辑软件，也把贝塞尔曲线作为仅有的矢量绘制工具（钢笔工具）包含其中。

Bezier曲线由一组控制点\(P_0\)到\(P_n\)定义，其中n称为其阶(n = 1表示线性，2表示二次，等等)。第一个控制点和最后一个控制点始终是曲线的端点;然而，中间控制点(如果有的话)通常不在曲线上。

一阶贝塞尔曲线(Linear Bezier curves)

给定控制点\(P_0\)和\(P_1\), 线性贝塞尔曲线是两点中的线段。

\[ B_t = P_0 + t(P_1 - P_0) = (1 - t)P_0 + tP_1 \qquad t \in[0, 1] \]

二阶贝塞尔曲线(Quadratic Bezier curves)

给定三个控制点\(P_0 \ P_1\)和\(P_2\), 取\(S_0\)为直线\(P_0 P_1\)一阶贝塞尔曲线，\(S_1\)为\(P_1 P_2\)一阶贝塞尔曲线，得：

\[ \begin{aligned} B_t &= (1 - t)\left[(1 -t)P_0\right] + t\left[ (1 -t)P_1 + tP_2\right] \qquad t \in[0, 1] \\ &= P_1 + (1-t)^2(P_0 - P_1) + t^2(P_2 - P_1) \end{aligned} \]

可推导出一阶导数为：

\[ \dot{b_t} = 2(1-t)(P_1 - P_0) + 2t(P_2 - P_1) \]

由此可见，二阶贝塞尔曲线在\(P_0\)和\(P_2\)处的切线相交于\(P_1\)。

可推导出二阶导数为：

\[ \ddot{B_t} = 2(P_2 - 2P_1 + P_0) \]

三阶贝塞尔曲线(Cubic Bezier curves)

给定四个控制点\(P_0 \ P_1 \ P_2\)和\(P_3\)，由以上推导，三阶贝塞尔曲线可以由二阶导出：

\[ B_t = (1 -t)B_t^{P_0,P_1,P_2} + tB_t^{P_1, P_2, P_3} \qquad t \in [0, 1] \]

展开可得：

\[ B_t = (1-t)^3P_0 + 3(1-t)^2tP_1 + 3(1-t)t^2P_2 + t^3P_3 \]

可推导出一阶导数为：

\[ \dot{B_t} = 3(1-t)^2(P_1 - P_0) + 6(1 -t)t(P_2 - P_1) + 3t^2(P_3 - P_2) \]

可推导出二阶导数为：

\[ \ddot{B_t} = 6(1-t)(P_2 - 2P_1 + P_0) + 6t(P_3 - 2P_2 + P_1) \]

N阶通用形式

\[ \begin{aligned} &B_t = \sum_{i = 0}^n b_i,n(t) P_i \qquad t \in [0,1] \\ &b_i,n(t) = \left( \begin{matrix} n \\ i \end{matrix}\right) t^i(1-t)^{n-i} \qquad i = 0, \cdots ,n \\ &\left( \begin{matrix} n \\ i \end{matrix}\right) = \frac{n!}{i!(n -i)!} \end{aligned} \]

其导数为：

\[ \dot{B_t} = n\sum_{i = 0}^{n-1}b_{i,n-1}(t)(P_{i+1} - P_i) \]

其中\(b_{i,n}\)称为Bernstein polynomial，\(P_i\)称为贝塞尔曲线控制点，由控制点组成的多边形为贝塞尔多边形，并且此多边形为贝塞尔曲线的凸包(convex hull)。

贝塞尔曲线性质

  • 曲线从\(P_0\)开始，到\(P_n\)结束，称为端点插值性质。
  • 当且仅当所有控制点共线时，曲线是一条直线。
  • 曲线的起点和终点分别与贝塞尔多边形的第一段和最后一段相切。
  • 一条曲线可以在任意点上分成两条子曲线，也可以分成任意多条子曲线，每条子曲线都是贝塞尔曲线。
  • 有些曲线看起来很简单，比如圆，不能用Bezier曲线或分段Bezier曲线精确地描述。
  • 每条二次贝塞尔曲线也是一条三次贝塞尔曲线，更一般的，每一条n阶贝塞尔曲线也是一条m阶贝塞尔曲线，其中m > n。比如n阶贝塞尔曲线具有控制点\(P_0, \cdots, P_n\)也是n+1阶贝塞尔曲线，其中\(P_k^{\prime} = \frac{k}{n + 1}P_{k-1} + (1+\frac{k}{n+1})P_k\)
  • 贝塞尔曲线具有变异递减性质。从直观的角度来看，这意味着一条Bezier曲线的“波动”不会超过它的控制点的多边形，而实际上它的“波动”可能比控制点的多边形还要小.
  • 对于某些控制点的选择，大于2阶的贝塞尔曲线可能会与自身相交，也可能有一个尖顶。

Uniform Cubic B-Splines in SE3

SE3定义为：

\[ { }_{a}^{b} \mathbf{T}=\left[\begin{array}{cc} { }^{b} \mathbf{R} & {\mathbf{p}}_{a} \\ { }^{a} & 1 \end{array}\right], \quad{ }_{a}^{b} \mathbf{T} \in \mathbb{S} \mathbb{E}(3), \quad{ }_{a}^{b} \mathbf{R} \in \mathbb{S} \mathbb{O}(3) \]

假设角速度和线速度在每个变换过程中都是恒定的。利用矩阵对数和指数，我们可以定义如下函数:

\[ {}^w_s \mathbf{T}(t)=\exp \left(\tilde{B}_{0, k}(t) \log \left( {}^w_0 \mathbf{T} \right)\right) \prod_{i=1}^{n} \exp \left(\tilde{B}_{i, k}(t) \log \left( {^w_{i-1} \mathbf{T}}^{-1} {}^w_{i} \mathbf{T} \right)\right) \]

其中\({ }^w_s \mathbf{T}(t)\)是t时刻插值的位姿，\({}^w_{i} \mathbf{T}\)为控制点，\(\tilde{B}_{i, k}(t)\)是基函数。

我们可以通过指定k = 4（degree=k - 1）来转换成三次b样条曲线。表示在最大值处有四个控制点，它们对样条曲线上的任何点都有贡献。

因此得到：

其中对于k=4，基函数为：

化解得到：

其中\(u(t) = (t - t_i)/ \Delta t \quad t_i = i \Delta t \quad u \in [0,1)\)。

简化表达为：

对时间求导：

因此可以得到：

Demo

代码实现参考open_vins, 其中蓝色曲线为knots插值之后的曲线，绿色曲线为模拟的IMU数据积分得到的曲线。

Reference

[1] https://en.wikipedia.org/wiki/Bézier_curve [2] https://www.desmos.com/calculator/xlpbe9bgll [3] K. Qin, “General matrix representations for b-splines,” The Visual Computer, vol. 16, no. 3, pp. 177–186, 2000. [4] https://github.com/HKUST-Aerial-Robotics/Fast-Planner [5] Continuous-Time Visual-Inertial Trajectory Estimation with Event Cameras