Spline

Bezier Curves

Lerp

由两个点组成的 Bezier Curve，定义表达式如下

$$P(t)=(1-t)P_0+tP_1$$

Quadratic Bezier Curve

由三个点组成的 Bezier Curve。该曲线中，每两个相邻点之间通过 Lerp 连接，得到中间点，而中间点经过 Lerp 处理之后可以得到最终的位置点

$$P(t)=(1-t)((1-t)P_0+tP_1)+t((1-t)P_1+tP_2)\\P(t)=(1-t)^2P_0+2t(1-t)P_1+t^2P_2$$

Cubix Bezier Curve

由四个点组成的 Bezier Curve。同上述 Quadratic Bezier Curve 类似，定义表达式如下

$$P(t)=(1-t)^3P_0+3t(1-t)^2P_1+3t^2(1-t)P_2+t^3P_3$$

上述的公式可以简化如下

$$P(t)=P_0+t(-3P_0+3P_1)+t^2(3P_0-6P_1+3P_2)+t^3(-P_0+3P_1-3P_2+P_3)\\P(t)=\begin{bmatrix}1&t&t^2&t^3\end{bmatrix}\begin{bmatrix}1&0&0&0\\-3&3&0&0\\3&-6&3&0\\-1&3&-3&1\end{bmatrix}\begin{bmatrix}P_0\\P_1\\P_2\\P_3\end{bmatrix}$$

其中中间的系数矩阵为特征矩阵，这个形式就是 Bernstein Polynomial Form

n-Degree Bezier Curve

$n$ 阶 Bezier Curve，由 $n+1$ 个点构成，表达式可以写作

$$P(t)=\sum_{i=0}^{n}(1-t)^{n-i}t^iP_i$$

De Casteljau’s Algorithm

Bezier Curve 通过这种方法来计算得到落点。对于一个 Cubix Bezier Curve 来说，计算公式如下

$$A=lerp(P_0,P_1,t)\\B=lerp(P_1,P_2,t)\\C=lerp(P_2,P_3,t)\\D=lerp(A,B,t)\\E=lerp(B,C,t)\\P=lerp(D,E,t)$$

其中 $lerp$ 表达式如下

$$lerp(A,B,t)=(1-t)A+tB$$

性质

1. 通过 Bezier Curve 的表达式可知，当有一个点改变时，将会影响整体的曲线状态
2. 曲线并不会通过设定的所有点，只能通过设定的起始点和终止点
3. 一旦点的数量多起来，计算量将会非常大，计算代价太过昂贵

Bezier Spline

由多段 Bezier Curves 首尾相接组成的样条曲线，一般会由 Cubix Bezier Curves 曲线连接而成

定义

1. 对于单段 Bezier Curves 来说，参数 $t$ 的范围是 $0\sim1$ ，在 Bezier Spline 中，参数范围将扩展到所有的曲线
2. 两段曲线之间的连接点被称为结点
3. 控制这些曲线的点被称之为控制点
4. 结点处的参数 $u$ 的值被称为结值
5. 节点之间的间距被称之为结间距
6. 对于均匀样条曲线来说，每一段结间隔都为 1，由 $n$ 段曲线组成的样条曲线的参数 $u$ 范围为 $0\sim n$ ，参数中整数部分表示曲线的索引，而小数部分表示每段曲线的局部 $t$ 值
7. 对于非均匀样条曲线来说，每一段的间隔都不是固定的，所以曲线的变化更加复杂

性质

1. 移动控制点只能在单段 Bezier Curves 中作用，移动结点也只能影响两段 Bezier Curves
2. 由于每一段曲线都是一段低阶数的 Bezier Curves，每次计算只需要计算该段 Bezier Curves 即可，所以计算量会比较小
3. 曲线通过所有的结点，但是并没有通过控制点。

结点

定义 $P_{i,j}$ 为第 $i$ 段曲线的第 $j$ 个控制点，并且假设每一段曲线由 $n+1$ 个控制点来控制。从而可以得到以下几种结点类型

1. broken： $\vec{P_{i-1,n-1}P_{i-1,n}}$ 与 $\vec{P_{i,0}P_{i,1}}$ 向量方向不一致
2. aligned： $\vec{P_{i-1,n-1}P_{i-1,n}}$ 与 $\vec{P_{i,0}P_{i,1}}$ 向量方向一致，但是其长度不同
3. mirrored： $\vec{P_{i-1,n-1}P_{i-1,n}}$ 与 $\vec{P_{i,0}P_{i,1}}$ 向量方向一致且长度一致

参数连续性

定义

对于上述的 Bezier Spline，在结点处可以发现，当参数 $u$ 从 0 开始逐渐增大时，当路过 broken 和 aligned 结点时，速度都会产生突变，也就是连续性问题。连续性是曲线连接程度的衡量标准

因此可以定义如下连续性的标准

1. 对于曲线中各部分曲线是不连接的，则没有连续性
2. 对于曲线是连接的，可以定义为 $C^0$ 连续
3. 如果曲线的一阶导数是连续的并且满足 $C^0$ 连续，则可定义为 $C^1$ 连续
4. 如果一条曲线的二阶导数是连续的并且满足 $C^1$ 连续，则可定义为 $C^2$ 连续
5. …

C1 连续性

定义第 $i$ 段曲线上的点为 $L_i(t)$ ，其中 $t$ 为每段曲线的参数。对于 $C^1$ 连续来说，所需要满足的条件为 $\dot{L_{i-1}}(1)=\dot{L_i}(0)$ ，也就是前一段曲线的终点的导数等于下一段曲线的起始点的导数。

对于一个由 $n$ 阶 Bezier Curves 组成的样条曲线来说，根据上述的式子可以得到 $P_{i,1}=2P_{i-1,n}-P_{i-1,n-1}$ ，根据此条件可以得到该结点的连接情况是属于 mirrored 类型的，所以也就是 mirrored 类型的结点处为 $C^1$ 连续的。但是如果要保证曲线的速度连续性，也就限制了 $P_{i,1}$ 点的位置，那么对曲线的控制也就受到了限制

C2 连续性

对于 $C^2$ 连续性来说，所需要满足的条件为 $\ddot{L_{i-1}}(1)=\ddot{L_i}(0)$

对于一个由 $n$ 阶 Bezier Curves 组成的样条曲线来说，根据上述的式子可以得到 $P_{i,2}=2P_{i-1,n-2}+4(P_{i-1,n}-P_{i-1,n-1})$ ，那么可以得到此时 $P_{i,2}$ 是完全受限制于 $P_{i-1,n-2},P_{i-1,n-1},P_{i-1,n}$ 点的，从而对曲线的控制造成限制

对于一个满足 $C^2$ 连续性的 Bezier Spline 来说，移动一个控制点的代价会非常大，这会导致整个曲线都收到一定程度的影响。而且想要满足 $C^2$ 连续性，就必须使用 4 阶及以上阶数的 Bezier Curves

几何连续性

曲率

给定曲线上的一点，有一个圆来描述该点的曲率。曲率用字母 $\kappa$ 表示，这个圆的半径为 $\frac{1}{\kappa}$

曲率的计算可以由该点的速度加速度行列式除以速度的三次方，如下

$$\kappa=\frac{\left|\begin{aligned}\dot{P}_x&&\ddot{P}_x\\\dot{P}_y&&\ddot{P}_y\end{aligned}\right|}{\Vert\dot{P}\Vert^3}$$

定义

由于参数连续性在几何设计中不太直观，并且依赖于参数的选择。而且对于上述中所说的 broken 结点满足 $C^0$ 连续性，而 mirrored 结点满足 $C^1$ 连续性。那么对于 aligned 点来说，这个点左右速度方向保持一致，只是速度的大小不同而已，此时定义该点的连续性为几何连续性

设 $\varphi(t)(a\leq t\leq b)$ 为一段给定的曲线，如果存在一个参数变换 $t=\rho(s)(a_1\leq s\leq b_1)$ ，使得 $\varphi(\rho(s))\in C^n[a_1,b_1])$ ，并且 $\frac{d\varphi(\rho(s))}{ds}\not=0$ ，则称曲线 $\varphi(t)(a\leq t\leq b)$ 为 $n$ 阶几何连续的曲线，记为 $\varphi(t)\in GC^n[a,b]$ 或者 $\varphi(t)\in G^n[a,b]$

性质

1. 条件 $\frac{d\varphi(\rho(s))}{ds}\not=0$ 保证了曲线上无奇点
2. 几何连续性于参数的选取无关，是曲线本身固有的几何性质
3. $G^n$ 的条件比 $C^n$ 的条件宽泛，曲线类型也更多

G1 连续性

对于 $G^1$ 连续性，需要满足的条件为 $\frac{\dot{L_{i-1}(1)}}{\Vert \dot{L_{i-1}(1)} \Vert}=\frac{\dot{L_{i}(0)}}{\Vert \dot{L_{i}(0)} \Vert}$ ，也就是连结点处的导数方向一致，并不考虑速度的大小。在几何上的表示为两曲线在连接点处有着公共的切线，即切线方向连续

与上述的 $C^1$ 连续性一样，用数学表达式可以得到控制点的关系式子，如下

$$P_{i,1}=P_{i-1,n}+(P_{i-1,n}-P_{i-1,n-1})\beta_1$$

其中 $\beta_1$ 可以是任意正值

G2 连续性

对于 $G^2$ 连续性，需要满足的条件自然是连结点处的加速度方向一致，不考虑加速度的大小。在几何上的表示为两曲线在连结点处有着公共的曲率圆，即曲率连续，曲率相等

根据曲率的计算公式，定义结点两端的曲率大小相等即可得到如下关系式

$$P_{i,2}=P_{i-1,n}+(P_{i-1,n}-P_{i-1,n-1})(2\beta_1+\beta_1^2+\frac{\beta_2}{2})+(P_{i-1,n-2}-P_{i-1,n-1})\beta_1^2$$

其中 $\beta_1,\beta_2$ 是任意的正值，此时曲线的轨迹也会受到一定限制，但是无论 $\beta_1,\beta_2$ 取任意值，都可以满足 $G^2$ 连续性的条件

与参数连续性的一些关联

1. $C^0$ 和 $G^0$ 连续性是一样的，都是曲线的连接
2. $C^n$ 连续性是 $G^n$ 连续性的充分条件，也就是当曲线满足 $C^n$ 连续性的时候，一定满足 $G^n$ 连续性
3. 对于 $G^n$ 连续性来说，曲线的速度不能为 0，否则曲线的曲率将会非常大，但是对于相对的 $C^n$ 连续性，曲线的速度是可以为 0 的

Hermite Spline

Hermite Curve

上述的 Cubix Bezier Curve 是给出了四个控制点来控制曲线的，但是对于 Hermite Curve 来说，给出的是端点坐标 $P_0,P_1$ 和端点处的速度向量 $V_0,V_1$

所以对于一段 Hermite Curve 就可以列出如下方程组

$$P(0)=P_0\\P(1)=P_1\\\dot{P}(0)=V_0\\\dot{P}(1)=V_1$$

有 4 个方程，所以对应的曲线的函数有 4 个未知数，如下

$$P(t)=at^3+bt^2+ct+d$$

方程求解可以得到如下的关系式子

$$P(t)=\begin{bmatrix}1&t&t^2&t^3\end{bmatrix}\begin{bmatrix}1&0&0&0\\0&1&0&0\\-3&-2&3&-1\\2&1&-2&1\end{bmatrix}\begin{bmatrix}P_0\\V_0\\P_1\\V_1\end{bmatrix}$$

性质

1. hermite spline 有着显式的导数
2. 通过分析可以得知，该曲线满足 $C^1$ 连续条件，但是不满足 $C^2$ 连续
3. 曲线穿过每个控制点，并且满足每个控制点上的速度

4. Hermite Curve 可以轻易的转换为 Bezier Curve。对于上述的点，转化之后的 Bezier Curve 的控制点为

   $$B_0=P_0\\B_1=P_0+\frac{V_0}{3}\\B_2=P_1-\frac{V_1}{3}\\B_3=P_1$$

5. 上述也说明 Bezier Curve 的连结点恰好对应于连结点处的速度的三分之一

Line Spline

定义

利用直线连接每一个控制点所形成的线

性质

1. 通过所有控制点
2. 弧长参数化很容易。弧长参数化是指以恒定速度的插值曲线
3. 计算代价很低，只是一个 Lerp Curve
4. 连续性只是 $C^0$ 连续的

Cardinal Spline

定义

对于 $n$ 个点来说，可以依次计算每个点处的速度，根据两个相邻点之间的连线来确定该点处的速度，如下

$$V_{i}=(P_{i+1}-P_{i-1})\beta_i$$

其中 $\beta$ 为比例因子，对于起点和终点来说，并没有两个相邻点，所以可以定义一个虚拟点 $P_{-1},P_{n+1}$

$$P_0-P_{-1}=P_1-P_0\\P_{n+1}-P_n=P_n-P_{n-1}$$

所以可以得到相应的起始点和终止点的速度

$$V_0=(P_1-P_{-1})\beta_0\\V_n=(P_{n+1}-P_{n-1})\beta_n$$

知道了每个点处的速度之后，可以使用 Hermite Curve 连接，从而可以得到一条曲线

如果对所有的速度使用同一个比例因子时，当比例因子 $\beta$ 不断减小时，曲线的连接处会不断变得尖锐从而不断趋近于 Line Spline，所以可以通过控制这个比例因子来控制曲线的锐度

对于一个四阶的 Cardinal Spline，其公式如下

$$P(t)=\begin{bmatrix}1&t&t^2&t^3\end{bmatrix}\begin{bmatrix}0&1&0&0\\-\beta&0&\beta&0\\2\beta&\beta-3&3-2\beta&-\beta\\-\beta&2-\beta&\beta-2&\beta\end{bmatrix}\begin{bmatrix}P_0\\P_1\\P_2\\P_3\end{bmatrix}$$

其中 $\beta$ 为比例因子，一般会将其称之为张力

性质

1. 穿过每一个控制点
2. 通过分析可以得知，该曲线满足 $C^1$ 连续条件，但是不满足 $C^2$ 连续
3. 曲线有着显式的导数

Catmull-Rom Spline

定义

对于上述所说的 Cardinal Spline，当将所有点的速度的比例因子调整为 $0.5$ 时，得到的曲线既没有特别平坦的地方，也没有不合理的急转，对于每个点也能平滑的通过

可以得到对应的矩阵表达式

$$P(t)=\frac{1}{2}\begin{bmatrix}1&t&t^2&t^3\end{bmatrix}\begin{bmatrix}0&2&0&0\\-1&0&1&0\\2&-5&4&-1\\-1&3&-3&1\end{bmatrix}\begin{bmatrix}P_0\\P_1\\P_2\\P_3\end{bmatrix}$$

性质

1. 曲线穿过所有控制点
2. 满足 $C^1$ 连续性
3. 不需要特别指定切线

C2 连续曲线

四阶 C2 连续曲线

首先定义曲线的矩阵表达式如下

$$P(t)=\begin{bmatrix}1&t&t^2&t^3\end{bmatrix}\begin{bmatrix}c_1&c_2&c_3&c_4\\c_5&c_6&c_7&c_8\\c_9&c_{10}&c_{11}&c_{12}\\c_{13}&c_{14}&c_{15}&c_{16}\end{bmatrix}\begin{bmatrix}P_0\\P_1\\P_2\\P_3\end{bmatrix}$$

由于点的位置为定值，所以可以判断前两个矩阵的情况即可，这两个矩阵的乘积称之为基础函数

为了满足 $C^2$ 连续性，定义如下方程，其中 $B_i(t)$ 表示第 $i$ 个点的基础函数的值

$$\left\{\begin{aligned}&B_0(1)=0\\&\dot{B_0}(1)=0\\&\ddot{B_0}(1)=0\\&B_1(1)=B_0(0)\\&\dot{B_1}(1)=\dot{B_0}(0)\\&\ddot{B_1}(1)=\ddot{B_0}(0)\\&B_2(1)=B_1(0)\\&\dot{B_2}(1)=\dot{B_1}(0)\\&\ddot{B_2}(1)=\ddot{B_1}(0)\\&B_3(1)=B_2(0)\\&\dot{B_3}(1)=\dot{B_2}(0)\\&\ddot{B_3}(1)=\ddot{B_2}(0)\\&B_3(1)=0\\&\dot{B_3}(1)=0\\&\ddot{B_3}(1)=0\\&B_0(t)+B_1(t)+B_2(t)+B_3(t)=1\end{aligned}\right.$$

解上述方程可以得到最终的解，如下

$$P(t)=\frac{1}{6}\begin{bmatrix}1&t&t^2&t^3\end{bmatrix}\begin{bmatrix}1&4&1&0\\-3&0&3&0\\3&-6&3&0\\-1&3&-3&1\end{bmatrix}\begin{bmatrix}P_0\\P_1\\P_2\\P_3\end{bmatrix}$$

性质

1. 曲线几乎没有通过任何控制点
2. 曲线达到了 $C^2$ 连续性

B-Spline

定义

B-Spline 指的是以 B 样条基函数为加权系数，对控制点进行线性组合所构造的曲线，实际上是分段连续的多项式曲线，具有 $n+1$ 个控制点

相对于 Bezier Spline 有几个优点

1. Bezier 一旦确定特征多边形就确定了曲线的阶数
2. Bezier 拼接复杂
3. Bezier 不能做局部修改，B-Spline 是整条曲线用一段一段的曲线拼接而成，采用分段连续多段式生成

解析定义

具有 $n+1$ 个控制顶点 $\{P_i|i=0,1,…,n\}$ ，定义在结点向量的 $u=\{u_i|i=0,…,n+k+1,u_i\leq u_{i+1}\}$ 的 $k$ 次 B-Spline 的解析定义如下

$$P(u)=\frac{\sum_{i=0}^{n}P_iB_i^k(u)}{\sum_{i=0}^nB_i^k(u)}\quad u\in[u_k,u_{n+1})$$

其中 $B_i^k$ 是控制点的加权系数

$$\left\{\begin{aligned}&B_i^0=\left\{\begin{aligned}&1&&u_i\leq u<u_{i+1}\\&0&&otherwise\end{aligned}\right.\\&B_i^k(u)=\frac{u-u_i}{u_{i+k}-u_i}B_i^{k-1}(u)+\frac{u_{i+k+1}-u}{u_{i+k+1}-u_{i+1}}B_{i+1}^{k-1}(u)\\&\frac{0}{0}=0\end{aligned}\right.$$

需要注意的是

1. 曲线参数取值 $u$ 一定不会小于前 $k$ 个节点值，不大于最后 1 个节点值
2. 多项式 $B_i^k(u)$ 是单位化的 B-Spline 基函数
3. B-Spline 的曲线控制顶点数量必须不小于曲线的阶数所需要的顶点数量，即 $n+1\geq k+1$
   1. 当 $n>k$ 时， $k$ 次 B-Spline 由 $k+1$ 段 $k$ 次多项式曲线构成
   2. 当 $n=k$ 时，B-Spline 退化为 Bezier Curve
4. $k$ 次 B-Spline 具有 $C^{k-1}$ 连续性
5. 当移动曲线的某一个控制顶点时，只会改变曲线的局部形状，而不会改变整条曲线的形状

struct Point {
  double x;
  double y;
  Point operator+ (const Point &p) { return (Point) { x + p.x, y + p.y }; }
};

const int k_ = 3;
const int n = 7;

double u[k_ + n + 1] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 };
Point  p[n + 1] = {
  {  0, 0 },
  { -1, 2 },
  {  0, 3 },
  {  2, 2 },
  {  4, 4 },
  {  5, 3 },
  {  5, 1 },
  {  8, 3 }
};

double base_func(const int &i, const int &k, const double &t, const double* u) {
  if (k == 0) {
    if (u[i] <= t && u[i + 1] > t) return 1.0;
    else return 0.0;
  } else {
    double len1 = (u[i + k] - u[i] == 0 ? 1 : u[i + k] - u[i]);
    double len2 = (u[i + k + 1] - u[i + 1] == 0 ? 1 : u[i + k + 1] - u[i + 1]);
    return (t - u[i]) / len1 * base_func(i, k - 1, t, u) +
    (u[i + k + 1] - t) / len2 * base_func(i + 1, k - 1, t, u);
  }
}

int main() {
  FILE* fp = fopen("./data.xlsx", "w");
  int   i;
  for (double t = 0; t < 10; t += 0.01) {
    double b[n + 1] = { 0 };
    double all_b = 0;
    for (i = 0; i < n + 1; ++i) {
      b[i] = base_func(i, k_ - 1, t, u);
      all_b += b[i];
    }
    all_b = (all_b == 0 ? 1 : all_b);
    Point p_now = { 0, 0 };
    for (i = 0; i < n + 1; ++i)
      p_now = p_now + (Point) { p[i].x * b[i] / all_b, p[i].y * b[i] / all_b };
    fprintf(fp, "%f\t%f\n", p_now.x, p_now.y);
  }
}

几何定义

B-Spline 曲线是根据控制顶点反复插值得到的，从低次到高次插值出 B-Spline 上的目标顶点的方法称之为 de Boor 算法

具有 $n+1$ 个控制顶点 $\{P_i|i=0,1,…,n\}$ ，定义在结点向量的 $u=\{u_i|i=0,…,n+k+1,u_i\leq u_{i+1}\}$ 的 $k$ 次 B-Spline 的几何定义如下

$$P(u)=P_l^k\quad u\in[u_l,u_{l+1}),l\in[k,n+1)$$

其中

$$P_i^j=\left\{\begin{aligned}&(1-\tau_i^j)P_{i-1}^{j-1}(u)+\tau_i^jP_i^{j-1}(u)&j>0\\&P_i&otherwise\\&\tau_i^j=\frac{u-u_i}{u_{i+k+1-j}-u_i}\end{aligned}\right.$$

性质

1. 未通过所有的控制点
2. $k$ 次 B-Spline 具有 $C^{k-1}$ 连续性

Non-Uniform B-Spline

非均匀 B-Spline。相对于 B-Spline 来说，每一段曲线中所用的时间并不是均匀的，可以随意改变，也就是对应的 $u(i+1)-u(i)$ 不为常数

NURBS

non-uniform rational basis spline 非均匀有理 B-Spline 与 B-Spline 相比，每一个控制顶点增加了一个加权系数 $\omega_i$ 来控制曲线的形状，具有 $n+1$ 个控制顶点 $\{P_i|i=0,1,…,n\}$ ，定义在结点向量的 $u=\{u_i|i=0,…,n+k+1,u_i\leq u_{i+1}\}$ 的 $k$ 次 NURBS 的解析定义如下

$$P(u)=\frac{\sum_{i=0}^{n}\omega_iP_iB_i^k(u)}{\sum_{i=0}^n\omega_iB_i^k(u)}\quad u\in[u_k,u_{n+1})\\=\sum_{i=0}^{n}P_i\frac{\omega_iB_i^k(u)}{\sum_{i=0}^n\omega_iB_i^k(u)}$$

需要注意的是

1. 多项式 $B_i^k(u)$ 是单位化的 B-Spline 基函数
2. 多项式 $\frac{\omega_iB_i^k(u)}{\sum_{i=0}^n\omega_iB_i^k(u)}$ 被称为有理 B-Spline 基函数
3. 当所有的加权系数 $\omega_i=1$ 时，NURBS 退化为 B-Spline
4. 当 $\omega_i=0$ 时，相应的控制点对曲线的形状没有影响
5. 当 $\omega_i\rightarrow\infty$ 时，曲线上的点 $P(u_{i+k})\rightarrow P_i$