贝塞尔曲线

Bezier曲线

贝塞尔曲线于 1962 年，由法国工程师皮埃尔·贝济埃（Pierre Bézier）所广泛发表，他运用贝塞尔曲线来为汽车的主体进行设计,贝塞尔曲线最初由保尔·德·卡斯特里奥于1959年运用德卡斯特里奥算法开发，以稳定数值的方法求出贝塞尔曲线.

先从几个简单的例子开始：

一阶贝塞尔曲线

对于一阶贝塞尔曲线为我们可以看到是一条直线，通过几何知识，很容易根据t的值得出线段上那个点的坐标。

给定点 $P_0$、$P_1$，线性贝塞尔曲线只是一条两点之间的直线。这条线由下式给出：

$$\mathbf{B}(t)=\mathbf{P}_0+(\mathbf{P}_1-\mathbf{P}_0)t=(1-t)\mathbf{P}_0+t\mathbf{P}_1,t\in[0,1]$$

一阶曲线就是很好理解, 就是根据t来的线性插值。$P_0$表示的是一个向量 $[x ,y]$, 其中$x$和$y$是分别按照这个公式来计算的。

二阶贝塞尔曲线

在平面内任选 3 个不共线的点，依次用线段连接。在第一条线段上任选一个点 $D$。计算该点到线段起点的距离$AD$，与该线段总长$AB$的比例。

根据上一步得到的比例，从第二条线段上找出对应的点 $E$，使得 $AD:AB = BE:BC$

这时候$DE$又是一条直线了, 就可以按照一阶的贝塞尔方程来进行线性插值了：

$$t= AD:AE$$

这时候就可以推出公式了：

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

$$P_1=(1-t)P_1+tP_2$$

整理一下公式, 得到二阶贝塞尔公式：

$$\mathbf{B}(t) =(1-t)^{2}\mathbf{P}_{0}+2t(1-t)\mathbf{P}_{1}+t^{2}\mathbf{P}_{2},t\in[0,1]$$

三阶贝塞尔曲线

四个点对应是三次的贝塞尔曲线. 分别在 $AB$、$BC$、$CD$ 之间采$EFG$点, $EFG$三个点对应着二阶贝塞尔, 在$EF$、$FG$之间采集$H$、$I$点来降阶为一阶贝塞尔曲线.

公式：

$$\mathbf{B}(t)=\mathbf{P}_0(1-t)^3+3\mathbf{P}_1t(1-t)^2+3\mathbf{P}_2t^2(1-t)+\mathbf{P}_3t^3,t\in[0,1]_0$$

动画效果

贝塞尔曲线公式

贝塞尔曲线是通过空间中的$n+1$个点$P_0, P_1, P_2, \ldots, P_n$来定义的，这些点称为控制点。控制点决定了曲线的形状，贝塞尔曲线定义如下：

$$\mathbf{C}(u)=\sum_{i=0}^nB_{n,i}(u)\mathbf{P}_i$$

其中系数定义如下：

$$B_{n,i}(u)=\frac{n!}{i!(n-i)!}u^{i}(1-u)^{n-i}$$

贝塞尔曲线上对应于参数位置为$u$的点，是所有控制点的"加权"平均值，其中权重是系数$B_{n,i}(u)$，其中$u$的定义域为$[0,1]$。函数$B_{n,i}(u), 0 \leq i \leq n$被称为贝塞尔基函数或伯恩斯坦多项式。

将控制点点$P_0P_1, P_2, \ldots, P_n$按此顺序连接，我们得到了一个控制多边形，它帮助我们描述了曲线的大致形状。

下图展示了由11个控制点定义的贝塞尔曲线，其中蓝色点是对应于$u=0.4$的曲线上的点。如图所示，曲线或多或少地跟随着多边形线。

---

问题：当$u$的定义域非[0,1]时如何处理？

贝塞尔曲线不总是定义在$[0,1]$这一标准区间上，有时其定义域可能为$[a,b]$。在这种情况下，就需要进行变量的转换。我们需要做的是，将$u$从区间$[a,b]$映射到$[0,1]$，然后在基函数中使用转换后的$u$值。具体的转换方法如下：

$$u'=\frac{u-a}{b-a}$$

将转换后的$u'$代入到基函数$B_{n,i}(u)$中，可以得到新的公式：

$$B_{n,i}(u')=\frac{n!}{i!(n-i)!}\left(\frac{u-a}{b-a}\right)^i\left(1-\frac{u-a}{b-a}\right)^{n-i}$$

这一变换后的基函数定义了在区间$[a,b]$上的贝塞尔曲线

贝塞尔曲线的特点和性质

• 贝塞尔曲线由$n+1$个控制点定义，曲线的度数是$n$：
  在每个基函数中，$u$的指数是$i + (n - i) = n$。因此，曲线的度数是$n$

• 曲线$C(u)$通过第一个控制点$P_0$和最后一个控制点$P_n$

• 非负性：所有基函数都是非负的。

• 单位分割：在任一特定的$u$值下，这些基函数值的总和为1。这保证了曲线上的任一点都是控制点加权平均的结果，且权重总和为1

  

  在上面的左侧图片中，展示了一个基于五个控制点构成的贝塞尔曲线。而右侧的图则展示了作为$u$值函数的五个基函数。

  这些图表揭示了当$u=0.5$时，这五个基函数的分布情况。图中最右侧的竖条图清晰地展现了如何将数值1等分为五部分，这就是所谓的单位分割。

  所有基函数的取值范围从0到1，它们的总和恰好为1，因此，这些基函数在计算加权平均时充当了权重的角色。更简单地说，计算C($u$)时，实际上是将控制点$P_i$按照$B_{n,i}(u)$的权重进行加权汇总。

• 凸包性质

  贝塞尔曲线的一个显著特点是，它完全被其$n + 1$个控制点所形成的凸包所包围。所谓凸包，指的是能够包含一组点集所有点的最小凸形状。

  在下面的示例中，11个控制点形成的凸包以灰色表示。值得注意的是，并不是所有控制点都位于凸包的边缘。比如，控制点3、4、5、6、8和9位于凸包内部。除了起点和终点，曲线始终位于这个凸包之内。

  

• 变差减少性

  在平面上，这个性质意味着没有任何直线能与贝塞尔曲线相交的次数超过它与曲线的控制多边形相交的次数。

  

  如上图所示，黄色直线与曲线相交3次，但与控制多边形相交了7次；
  红色直线与曲线相交5次，与控制多边形也相交7次；
  而青色直线与曲线及其控制多边形都只相交两次。
  你可以尝试绘制其他直线以验证这一性质。

  若曲线是空间中的曲线，则只需将“直线”改为“平面”。

  这一性质的意义何在？它告诉我们，与控制多边形相比，贝塞尔曲线的复杂性（即转弯和扭曲的程度）会更低。换言之，控制多边形相较于它所定义的贝塞尔曲线，会有更频繁的转弯和扭曲，因为任意直线与控制多边形相交的次数总是多于与曲线相交的次数。

• 仿射不变性

  对贝塞尔曲线进行仿射变换时，其变换结果能够通过控制点经过相同变换后的新位置来确定。这一特性极为有利。当我们需要对贝塞尔曲线执行几何或仿射变换操作时，直接对控制点施加变换即可。完成控制点的变换后，由这些调整后的控制点所定义的新贝塞尔曲线即为所求。因此，无需直接对曲线本身进行变换处理。

• 递归性

  递归性指其系数满足下式：

  $$B_{i,n}(t)=(1-t)B_{i,n-1}(t)+tB_{i-1,n-1}(t)\quad(i=0,1,\cdots,n)$$

• 一阶导数性质

  

  假设上图中贝塞尔的 $t$ 是由左到右从 $0$ 到 $1$ 增加的，那么贝塞尔曲线在 $t=0$ 时的导数是和 $P_0P_1$的斜率（导数）是相同，$t=1$ 时的导数是和$P_3P_4$ 的斜率（导数）是相同的

  

  这一点的性质可以用在贝塞尔曲线的拼接，只要保证三点一线中的中间点是两段贝塞尔曲线的连接点，就可以保证两端贝塞尔曲线的导数连续

控制点的移动对曲线的影响

当我们调整一个贝塞尔曲线的控制点位置时，曲线形状将发生变化。这引出一个问题：
控制点移动到新的位置时，贝塞尔曲线将如何变形？

假设控制点 $P_k$ 移动到 $P_k+\mathbf{v}$ 的新位置，向量 $\mathbf{v}$ 表示移动的方向和距离。如图所示：

原始贝塞尔曲线表示为：

$$\mathbf{C}(u)=\sum_{i=0}^nB_{n,i}(u)\mathbf{P}_i$$

由于新的贝塞尔曲线是由 $P_0, P_1, ..., P_k+\mathbf{v}, ..., P_n$定义的，新的贝塞尔曲线方程 $D(u)$ 可以写成：

$$\begin{aligned} \mathbf{D}(u)&=\sum_{i=0}^{k-1}B_{n,i}(u)\mathbf{P}_{i}+B_{n,k}(u)(\mathbf{P}_{k}+\mathbf{v})+\sum_{i=k+1}^{n}B_{n,i}(u)\mathbf{P}_{i} \\ &=\sum_{i=0}^{n}B_{n,i}(u)\mathbf{P}_{i}+B_{n,k}(u)\mathbf{v} \\ &=\mathbf{C}(u)+B_{n,k}(u)\mathbf{v} \end{aligned}$$

在上式中，因为只有第 $k$ 项控制点使用了 $P_k+\mathbf{v}$，所以我们得出新曲线是原始曲线与额外的项 $B_{n,k}(u)\mathbf{v}$ 的组合。

因此，新曲线上任一给定 $u$ 对应的点，可以通过将原始曲线上相同位置的点沿向量 $\mathbf{v}$ 的方向平移 $|B_{n,k}(u)\mathbf{v}|$ 的距离来获得。

更具体地来说，对于任意给定的 u 值，原曲线上对应的点为 $C(u)$，新曲线上对应的点为 $D(u)$，且 $D(u) = C(u) + B_{n,k}(u)\mathbf{v}$。

以下图说明这一变化：

图中，黑色与红色曲线分别代表由9个控制点定义的同一8次贝塞尔曲线的原始和变形状态。黑色曲线为原始状态。若其中一个控制点沿蓝色向量移动到新位置，黑色曲线便转变为红色曲线。

图中同时标出了$u=0.5$时两条曲线上的对应点。可以明显看出，$C(0.5)$ 沿相同方向移动至 $D(0.5)$。从 $C(0.5)$ 到 $D(0.5)$ 的距离等同于向量

$$B_{8,3}(0.5)\mathbf{v} = \frac{8!}{(3!(8-3)!)}×0.5^3(1-0.5)^{8-3}\mathbf{v} = 0.22\mathbf{v}$$

因此，该距离约为原始控制点3和下图所示的新控制点3之间距离的22%

从这一分析中，我们可以得出另一个重要的结论：
因为 $B_{n,k}(u)$ 在开区间 $(0,1)$ 内非零，因此 $B_{n,k}(u)\mathbf{v}$ 在此区间内不是零向量。这意味着除了端点 $C(0)$ 和 $C(1)$ 外，所有原始曲线上的点都会被移动到新的位置上。由此可见：

​ 移动控制点将会全面改变贝塞尔曲线的形状。

计算贝塞尔曲线上一点：De Casteljau 算法

在构造 Bézier 曲线后，如何找到曲线上参数为 $u$ 处的点 $C(u)$？

一种简单的方法是将 $u$ 代入贝塞尔曲线的公式：

$$\mathbf{C}(u)=\sum_{i=0}^nB_{n,i}(u)\mathbf{P}_i$$

其中系数定义如下：

$$B_{n,i}(u)=\frac{n!}{i!(n-i)!}u^{i}(1-u)^{n-i}$$

计算每个基函数的值，并计算每个基函数及其对应控制点的乘积，最后将它们相加。虽然这样可以，但不是数值稳定的（即在计算伯恩斯坦多项式的过程中可能引入数值误差）。

还记得本文开始几个小节关于一阶、二阶、三阶贝塞尔曲线的推导过程吗，这就是De Casteljau 算法：

首先，我们以编号的形式记录下控制点，如 $P_0$ 记作 $00$，$P_1$ 记作 $01$，依此类推，直至 $P_n$ 的编号为 $0n$。第一个数字 $0$ 表示第0阶或第一次迭代。在接下来的迭代中，这些数字会变成$1、2、3$ 等。

De Casteljau 算法的核心在于找到线段 AB 上某个点 $C$，该点按 $u:1-u$ 的比例将线段 $AB$ 划分。从而确立点 $C$ 的具体位置。

线段$AB$ 的向量可以表示为 $B - A$。由于 $u$ 是介于 $0$ 和$1$ 之间的比例，所以点 $C$ 的位置可以用 $u(B - A)$ 来描述。再加上 $A$ 点的坐标，我们可以得到 $C$ 的准确位置为

$$A + u(B - A) = (1 - u)A + uB$$

因此，对于任何给定的 $u$ 值， 我们都能准确地定位在 $A$ 与 $B$ 之间的一点 $C$。

如何找到贝塞尔曲线上任意一处参数 $u$ 的点 $C(u)$：

首先将控制点依次连接，得到一个折线：即 $00-01-02-03...-0n$ 多段线，然后根据上述公式在每对相邻控制点之间找到新点 $1i$，这个新点将原线段 $0i-0(i+1)$ 按 $u:1-u$ 的比例分割。通过这种方式，我们可以得到 $n$ 个新点：$10, 11, 12, ..., 1(n-1)$，它们共同构成了一条新的折线。

如上图所示，例如，当 $u = 0.4$ 时，通过插值可以得到点$10$位于起始控制点$00$和$01$之间， 其中$\frac{00-10}{10-01}=\frac{u}{1-u}$ 。($00$即点$0$，$00$中的第一个$0$表示$0$阶，第二个$0$表示第$0$个点。$01$即点$1$，$01$中的第一个$0$表示$0$阶，第二个$1$表示第$1$个点）

第二次迭代：点$11$位于点$01$和$02$之间，$\frac{01-11}{11-02}=\frac{u}{1-u}$以此类推，可以计算出点$12、13、14$。这些新产生的点都以蓝色标记。

每一次迭代，我们获得的"新控制点"的数量都比之前迭代的结果少1，直到最减少到唯一的点$n0$。这一点就是在给定 $u$ 值下，贝塞尔曲线上的对应点 $C(u)$。

举例：在上图中，
点$20$可以通过以 $u:1-u$ 的比例划分线段$10-11$得到，。
点$21$通过边$11-12$划分得到
点$22$通过边$12-13$划分得到
点$23$通过边$13-14$划分得到。
从而得到了由点$20,21,22,23$定义的三阶折线。
在折线段$20-21-22-23$的两两点之间进行插值，将得到个新的折线段$30-31-32$。继续迭代，将得到点$40$和$41$， 组成新的折线段$40-41$。
继续迭代，将得到点$50$，这就是曲线上的点 $C(0.4)$。

用一张图表示De Casteljau 算法的迭代过程：

从初始列，即第 $0$ 列，我们通过两两点之间的插值可以计算出第$1$列的"新点"，从第$1$列我们得到第$2$列，依此类推。最终，在应用了 $n$ 次之后，我们将得到唯一的点 $n0$，这就是贝塞尔曲线上的点$C(u)$。

递归关系：自顶向下的计算方法

上述计算可以递归地表达。初始时，令 $P_{0,j}$ 为 $P_j$，其中 $j = 0, 1, ..., n$。$P_{0,j}$ 表示第 0 列的第 $j$ 项。那么第 $i$ 列的第 $j$ 项的计算如下：

$$\left.\mathbf{P}_{i,j}=(1-u)\mathbf{P}_{i-1,j}+u\mathbf{P}_{i-1,j+1}\quad\left\{\begin{array}{l}i=1,2,\ldots,n\\j=0,1,\ldots,n-i\end{array}\right.\right.$$

基于这个公式，可能会立即想到以下递归算法：

function deCasteljau(i,j)
if i = 0 then
	return P_{0,j}
else
	return (1-u) * deCasteljau(i-1,j) + u * deCasteljau(i-1,j+1)

这个算法看起来简单，但是它极其低效。原因如下。我们从调用 $deCasteljau(n,0)$ 开始，计算 $P_{n,0}$。$else$ 分支将这个调用会递归产生两个新的调用，即$deCasteljau(n-1,0)$ 用于计算 $P_{n-1,0}$ 和 $deCasteljau(n-1,1)$ 用于计算 $P_{n-1,1}$。

对 $deCasteljau(n-1,0)$ 的调用。又会产生两个新的调用，$deCasteljau(n-2,0)$ 用于计算 $P_{n-2,0}$ 和 $deCasteljau(n-2,1)$ 用于计算 $P_{n-2,1}$。对 $deCasteljau(n-1,1)$ 的调用也会产生两个新的调用，$deCasteljau(n-2,1)$ 用于计算 $P_{n-2,1}$ 和 $deCasteljau(n-2,2)$ 用于计算 $P_{n-2,2}$。

因此，$deCasteljau(n-2,1)$ 被调用了两次。依次类推，会发现计算 $P_{i,j}$ 的函数被重复调用了，这是对性能非常浪费的。

三角形计算方案：自底向上的计算方法

以下面在由 8 个控制点 $00, 01, ..., 07$定义的 7 次贝塞尔曲线为例。给定一个在 $[0,1]$ 中的 $u$，我们如何计算这个贝塞尔曲线上对应的点呢？曲线上的点是由选定点形成的等边三角形的底边对面的顶点！

例如，如果选定的控制点点是 02, 03, 04 和 05，那么对应于 $u$ 的这四个控制点定义的曲线上的点是 32。即蓝色的三角形的顶点。如果选定的控制点是 11, 12 和 13，则曲线上的点是 31。即黄色的三角形的顶点。如果选定的点是 30, 31, 32, 33 和 34，曲线上的点是 70。

反之同理，70 是由控制点 60 和 61 定义的贝塞尔曲线上的点。它也可以是由 50, 51 和 52 控制点定义的曲线上的点，也可以是由 40, 41, 42 和 43 控制点定义的曲线上的点。

如此我们就可以自底向上的计算每一项$P_{i,j}$，以蓝色三角形为例，先计算最左列的$P_{0,2},P_{0,3},P_{0,4},P_{0,5}$，再计算第二列的$P_{1,2},P_{1,3},P_{1,4},P_{0,5}$，再计算第三列的$P_{2,2},P_{2,3}$，再计算最终点$P_{3,2}$，这样可以避免自定向下计算带来的重复项计算的问题。

贝塞尔曲线的导数

贝塞尔曲线的导数仍然是贝塞尔曲线

要在贝塞尔曲线上的某点计算切向量和法向量，首先要计算在该点的一阶和二阶导数。

由 $n + 1$ 个控制点 $P_0, P_1, ..., P_n$ 定义的贝塞尔曲线如下：

$$\mathbf{C}(u)=\sum_{i=0}^nB_{n,i}(u)\mathbf{P}_i$$

其中 $B_{n,i}(u)$ 定义如下：

$$B_{n,i}(u)=\frac{n!}{i!(n-i)!}u^{i}(1-u)^{n-i}$$

其中控制点$P_{i}$是常数且与变量 $u$ 相互独立，因此，计算贝塞尔曲线的导数 $C'(u)$可以简化为计算 $B_{n,i}(u)$ 的导数。

通过一些代数操作，我们可以得到：

$$\frac{\mathrm{d}}{\mathrm{d}u}B_{n,i}(u)=B_{n,i}^{\prime}(u)=n(B_{n-1,i-1}(u)-B_{n-1,i}(u))$$

因此曲线的导数 $C'(u)$为：

$$\frac{\mathrm{d}}{\mathrm{d}u}\mathbf{C}(u)=\mathbf{C}'(u)=\sum_{i=0}^{n-1}B_{n-1,i}(u)\{n(\mathbf{P}_{i+1}-\mathbf{P}_i)\}$$

设 :

$$Q_0 = n(P_1 - P_0)$$

$$Q_1 = n(P_2 - P_1)$$

$$Q_2 = n(P_3-P_2)$$

$$Q_{n-1} = n(P_n - P_{n-1})$$

上述方程简化为以下形式：

$$\mathbf{C}^{\prime}(u)=\sum_{i=0}^{n-1}B_{n-1,i}(u)\mathbf{Q}_{i}$$

结论：可以看出贝塞尔曲线 的导数$C'(u)$是由 $n$ 个控制点 $n(P_1 - P_0)$, $n(P_2 - P_1)$, $n(P_3 - P_2)$, …, $n(P_n - P_{n-1})$ 定义的 $n - 1$ 次贝塞尔曲线。

这个导数曲线通常被称为原贝塞尔曲线的轨迹图。

注意，$P_{i+1} - P_i$ 是从 $P_i$ 到 $P_{i+1}$ 的方向向量，而 $n (P_{i+1} - P_i)$ 是该方向向量的 $n$ 倍长。下面左图显示了一个7次贝塞尔曲线，右图显示了其导数，是一个6次贝塞尔曲线。

贝塞尔曲线与其第一条和最后一条边相切

令$u = 0$和$u = 1$，可以得到：

$$C'(0) = n(P_1 - P_0)$$

$$C'(1) = n(P_n - P_{n-1})$$

第一个公式意味着在$u = 0$时的切向量方向是$P_1 - P_0$乘以$n$。因此，指示方向的第一条边与贝塞尔曲线相切。
第二个公式意味着在$u = 1$时的切向量方向是$P_n - P_{n-1}$乘以$n$。因此，指示方向的最后一条边与贝塞尔曲线相切。下图很好地展示了这个属性。

用$C^1$连续性连接两个贝塞尔曲线

利用贝塞尔曲线与其第一条和最后一条边相切的这个特性，可以将两条或更多的贝塞尔曲线连接起来，组合成我们想要的形状。

设第一条曲线$C(u)$由$m + 1$个控制点$P_0, P_1, P_2, ..., P_m$定义，
设第二条曲线$D(u)$由$n + 1$个控制点$Q_0, Q_1, Q_2, ..., Q_n$定义，
如果我们想要将这两条贝塞尔曲线连接起来，那么$P_m$必须等于$Q_0$，这保证了$C^0$级连续

第一条曲线与其最后一条边相切，第二条曲线与其第一条边相切。因此，为了实现平滑过渡，$P_{m-1}$，$P_m = Q_0$和$Q_1$四个点必须在同一直线上，使得从$P_{m-1}$到$P_m$的方向与从$Q_0$到$Q_1$的方向相同，如下图所示：

以这种方式连接两个贝塞尔曲线看起来很平滑，但它仍然是$C^0$连续，而不是$C^1$连续。但是，它满足$G^1$连续（具有相同的切向量方向，但切向量大小不相等）。
为了保证$C^1$连续性，我们必须确保第一条曲线在$u = 1$处的切向量$C'(1)$，与第二条曲线在$u = 0$处的切向量$D'(0)$，是等的。也就是说，必须满足以下关系：

$$\mathbf{C}^{\prime}(1)=m(\mathbf{P}_m-\mathbf{P}_{m-1})=\mathbf{D}^{\prime}(0)=n(\mathbf{Q}_1-\mathbf{Q}_0)$$

这个关系说明要在连接点保证$C^1$连续性，第一条曲线的最后一条边的长度（即$|p_m - p_{m-1}|$）与第二条曲线的第一条边的长度（即$|q_1 - q_0|$）的比例必须是$n/m$。由于度数$m$和$n$是固定的，我们可以调整$p_{m-1}$或$q_1$在同一直线上的位置，使得上述关系得到满足。

上面左图中为两条贝塞尔曲线连接而成，其中浅色的折线定义了一条4次贝塞尔曲线，而深色的折线定义了一条5次贝塞尔曲线。
由于第一条曲线的最后一段和第二条曲线的第一段不在同一直线上，所以这两条曲线并未平滑相接。
上面右图展示了在连接点处与一条直线相切的两条贝塞尔曲线。然而，它们并不是$C^1$连续的。左边的曲线是4次的，而右边的曲线是7次的。
但是，左曲线的最后一段与第二条曲线的第一段的比例看起来接近1，它们是$G^1$连续的而不是$C^1$连续的。为了实现$C^1$连续性，我们应该增加左曲线的最后段的长度（或者减少右曲线的第一段的长度），使其比例接近$7/4=1.75$，从而满足$C^1$连续

构造闭合的贝塞尔曲线

利用上面的这种切线性质，如果让第一个和最后一个控制点相同（即$P_0 = P_n$）并且$P_1, P_0$和$P_{n-1}$共线，那么生成的贝塞尔曲线将是一个封闭曲线，并且在连接点处是$G^1$连续的。

如果要在$P_0$处达到$C^1$连续性，必须使 $P_1-P_0=P_n-P_{n-1}$（即曲线的第一段和最后一段有相同的长度并且$P_1, P_0 = P_n, P_{n-1}$共线）

导数与de Casteljau算法之间的关系

重新改写一下贝塞尔曲线的导数形式：

$$\begin{aligned} \mathrm{C}^{\prime}(u)&=\sum_{i=0}^{n-1}B_{n-1,i}(u)\{n(\mathbf{P}_{i+1}-\mathbf{P}_i)\} \\ &=n\left[\left(\sum_{i=0}^{n-1}B_{n,i}(u)\mathbf{P}_{i+1}\right)-\left(\sum_{i=0}^{n-1}B_{n,i}(u)\mathbf{P}_{i}\right)\right] \end{aligned}$$

贝塞尔曲线的导数可以看作是两条$n-1$次贝塞尔曲线的差
现在让这两条曲线分别为$C_1(u)$和$C_2(u)$：

$$\quad{C}_1(u)=\sum\limits_{i=0}^{n-1}B_{n,i}(u)\mathbf{P}_{i+1}\newline{C}_2(u)=\sum\limits_{i=0}^{n-1}B_{n,i}(u)\mathbf{P}_i$$

根据上面定义，可以知道第一条曲线$C_1(u)$由控制点$P_1, P_2, ..., P_n$定义，第二条曲线$C_2(u)$由控制点$P_0, P_1, ..., P_{n-1}$定义，因此，原曲线的导数定义为：

$$\mathbf{C}^{\prime}(u)=n(\mathbf{C}_1(u)-\mathbf{C}_2(u))$$

为了计算贝塞尔曲线的导数$C'(u)$，可以先使用de Casteljau算法计算$C_1(u)$和$C_2(u)$的值，然后，再计算它们的差，乘以n就得到了$C'(u)$。

计算$C_1(u)$和$C_2(u)$：

如下图，最左侧的列有所有给定的控制点，点$n0$代表原始贝塞尔曲线上的点$C(u)$，点$(n-1)0$和$(n-1)1$连成的线段是用来插值获取点$n0$的控制折线，即点$n0$位于$(n-1)0$和$(n-1)1$的线段上，并以$u:1-u$的比例将其分割。

由于曲线$C_1(u)$由控制点$01, 02, ..., 0n$定义，因此点$(n-1)1$是曲线$C_1(u)$上的点。同样地，由于$C_2(u)$由控制点$00, 01, ..., 0(n-1)$定义，因此点$(n-1)0$是曲线$C_2(u)$上的点，向量$C_1(u) - C_2(u)$表示从点 $(n-1)0$ 到点 $(n-1)1$ 的向量，而原贝塞尔曲线的导数$C'(u)$的是$n$倍的$C_1(u) - C_2(u)$。

这说明，从 $(n-1)0$ 到 $(n-1)1$ 的线段与$C(u)$处的切线向量方向相同，它在$C(u)$处与曲线相切！

结论：de Casteljau算法迭代中的最后一条折线，实际上是一条线段，在$C(u)$处与贝塞尔曲线相切。

下图是一个例子，所示的点为$u = 0.5$处。显然，de Casteljau迭代过程中的最后一条折线段在$C(0.5)$处与曲线相切。

计算贝塞尔的高阶导数

$C(u)$的一阶导数如下：

$$\mathbf{C}^{\prime}(u)=\sum_{i=0}^{n-1}B_{n-1,i}(u)\mathbf{Q}_i$$

将导数公式应用于上述贝塞尔曲线，得到以下结果，给出了原始贝塞尔曲线的二阶导数：

$$\begin{aligned} \mathbf{C}^{\prime\prime}&=\sum_{i=0}^{n-2}B_{n-2,i}(u)\{(n-1)(\mathbf{Q}_{i+1}-\mathbf{Q}_i)\}\\ &=\sum_{i=0}^{n-2}B_{n-2,i}(u)\{n(n-1)(\mathbf{P}_{i+2}-2\mathbf{P}_{i+1}+\mathbf{P}_i)\} \end{aligned}$$

为了简洁地表达高阶导数，可以使用有限差分的方式。定义$D_{i}^0$为控制点$P_i$本身，其中$0 \leq i \leq n$。定义第一层差分$D_{i}^1$为前一层的差分，如下：

$$\begin{aligned} \mathbf{D}_{0}^{1}&=\text{ D}_1^0-\mathbf{D}_0^0=\mathbf{P}_1-\mathbf{P}_0 \\ \mathrm{D}_{1}^{1}&=\text{ D}_2^0-\mathbf{D}_1^0=\mathbf{P}_2-\mathbf{P}_1 \\ &\vdots \\ \mathbf{D}_{n-2}^{1}&=\text{ D}_{n-1}^0-\mathbf{D}_{n-2}^0=\mathbf{P}_{n-1}-\mathbf{P}_{n-2} \\ \mathbf{D}_{n-1}^{1}&=\text{ D}_n^0-\mathbf{D}_{n-1}^0=\mathbf{P}_n-\mathbf{P}_{n-1} \end{aligned}$$

第一层差分只有$n$个点，比给定的控制点少一个，简写一下，第一层的差分如下所示：

$$\mathbf{D}_i^1=\mathbf{D}_{i+1}^0-\mathbf{D}_i^0\quad0\leq i\leq n-1$$

第二层差分被定义为第一层点的差分：

$$\mathbf{D}_i^2=\mathbf{D}_{i+1}^1-\mathbf{D}_i^1\quad0\le i\le n-2$$

这个第二层差分有$n-1$个点。重复这个过程，我们可以定义第$k$层差分如下。第$k$层差分有$n-k+1$个点。

$$\mathrm D_i^k=\mathrm D_{i+1}^{k-1}-\mathrm D_i^{k-1}\quad0\leq i\leq n-k$$

有了这些定义，可以简洁地表达更高阶的导数。首先，让我们使用$D_{i}^1$重新写下$C'(u)$：

$$\begin{aligned} \mathbf{C}^{\prime}(u)&=n\sum_{i=0}^{n-1}B_{n-1,i}(u)(\mathbf{P}_{i+1}-\mathbf{P}_{i}) \\ &=n\sum_{i=0}^{n-1}B_{n-1,i}(u)(\mathbf{D}_{i+1}^0-\mathbf{D}_i^0) \\ &=n\sum_{i=0}^{n-1}B_{n-1,i}(u)\mathbf{D}_i^1 \end{aligned}$$

将导数应用于$C'(u)$来计算$C''(u)$，得到了一个使用$D_{i}^2$表示$C''(u)$的公式：

$$\begin{aligned} \mathbf{C}''(u)&=n(n-1)\sum_{i=0}^{n-2}B_{n-2,i}(u)(\mathbf{D}_{i+1}^1-\mathbf{D}_i^1)\\ &=n(n-1)\sum_{i=0}^{n-2}B_{n-2,i}(u)\mathbf{D}_i^2 \end{aligned}$$

对$C''(u)$做同样的处理，我们得到了一个使用$D_{i}^3$表示$C'''(u)$的公式：

$$\begin{aligned} \mathbf{C}^{\prime\prime\prime}(u)&=n(n-1)(n-2)\sum_{i=0}^{n-3}B_{n-3,i}(u)(\mathbf{D}_{i+1}^2-\mathbf{D}_i^2)\\ &=n(n-1)(n-2)\sum_{i=0}^{n-3}B_{n-3,i}(u)\mathbf{D}_i^3 \end{aligned}$$

继续这个过程，我们可以计算出第$k$阶导数$C^{(k)}(u)$使用$D_{i}^k$的表示形式：

$$\begin{aligned} \mathbf{C}^{[k]}(u)&=n(n-1)(n-2)\cdots(n-k+1)\sum_{i=0}^{n-k}B_{n-k,i}(u)(\mathbf{D}_{i+1}^{k-1}-\mathbf{D}_{i}^{k-1})\\ &=n(n-1)(n-2)\cdots(n-k+1)\sum_{i=0}^{n-k}B_{n-k,i}(u)\mathbf{D}_{i}^{k} \end{aligned}$$

因此，要在特定的$u$处计算$C(u)$的第$k$阶导数，我们首先计算$D_{i}^k$的值，然后应用de Casteljau算法计算由$D_{i}^k$定义的贝塞尔曲线上对应于$u$的点。

如下图所示，我们已有第0列所有$D_{i}^0$的值（即控制点本身），然后使用两两差分来计算第1列的$D_{i}^1$的值。重复这个过程直到我们达到如下所示的第$k$列：

最后，使用de Casteljau算法计算出第$k$列点的$B_{n-k,i}(u)$的值，与$D_{i}^k$相乘后累加，再将结果乘以$n(n-1)(n-2)...(n-k+1)$，就得到了$C^{(k)}(u)$的值！

贝塞尔曲线的拆分

曲线拆分的过程是将一个贝塞尔曲线在特定点$C(u)$处一分为二，每一部分仍然是贝塞尔曲线。分割后的每段曲线都需要一套新的控制点来描述其形状，原有的控制点将不再适用。此外，由于原始贝塞尔曲线的阶次是$n$，被切割成两部分后，每个部分都是原贝塞尔曲线的子集，因此得到的新贝塞尔曲线的阶次也必须是$n$。

假设我们$n + 1$个控制点$P_0, P_1, P_2, ..., P_n$，参数 $0\leq u\leq1$，目标是找到两组新的$n+1$个控制点$Q_0, Q_1, Q_2, ..., Q_n$和$R_0, R_1, R_2, ..., R_n$。

这两组控制点各自定义的贝塞尔曲线分别对应原始曲线在$[0,u]$和$[u,1]$区间的部分。

虽然正确性证明有点繁琐，但实际上算法非常简单。事实上，de Casteljau 的算法用于计算曲线上点 $C(u)$ 时已经提供了所有必要的信息。在下面的左图中，展示了用 de Casteljau 算法计算 $C(u)$ 的所有中间步骤，右图显示了在点 $C(u)$ 处曲线拆分后的形状及相应的控制折线。

仔细比较这两幅图，左图中，左半段贝塞尔曲线的控制点由$P_{00}=P_0, P_{10},P_{20},P_{30},P_{40},P_{50}$和$P_{60} = C(u)$组成
右半段贝塞尔曲线的控制点由点$P_{60}=C(u),P_{51},P_{42},P_{33},P_{24},P_{15}$和$P_{06} = P_6$组成。
下图展示了这些点：

回想一下de Casteljau算法的三角形计算方案：对于给定的$u$，计算$C(u)$需要$n$次迭代，在计算过程中，可以收集每列的第一个点和最后一个点，并且最终，第一个点（或最后一个点）的集合对应于原始曲线在$[0,u]$（或$[u,1]$）上定义的片段。

因此，在下面的三角形方案中，箭头指向的顶边(红色：$00,10,20,30,40,50,60$)和箭头相反方向的底边(蓝色：$60,51,42,33,24,15,06$)分别对应位第一段曲线和第二段曲线段的控制点。有了控制点，也就可以得出拆分后的两段贝塞尔曲线（可以重新计算获得）

注意：点$50$和$51$定义的线段在点$60$处与曲线相切（见上图），即左曲线的最后一段（即线段$50-60$）与左曲线相切，右曲线的第一段（即线段$60-51$）与右曲线相切。

曲线拆分为何至关重要？
曲线拆分的应用范围广泛。例如，它可以被用来找出两个贝塞尔曲线的交汇点、对贝塞尔曲线进行渲染，以及简化曲线设计过程。设想我们设计了一条曲线，却发现它并不完全符合我们的预期。在这种情形下，我们可能会考虑在某些特定的点上将曲线分割为两个部分：一个部分达到了我们的满意度，而另一部分则没有。接下来，我们可以忽略掉那部分我们已经满意的，将注意力集中在仍需改进的部分上。

贝塞尔曲线的升阶

不改变贝塞尔曲线的形状，同时增加贝塞尔曲线的度数（即阶数/次数），称为贝塞尔曲线的升阶。

假设我们有一个由 $n+1$个控制点$P_0, P_1, P_2, ..., P_n$定义的$n$阶贝塞尔曲线，我们希望将这个曲线的度数提升到$n + 1$而不改变其形状。由于$n + 1$阶贝塞尔曲线需要由$n + 2$个控制点定义，因此，我们需要找到一个新的控制点集，显然，$P_0$ 和 $P_n$ 必须在新控制点集中，因为贝塞尔曲线会通过首尾控制点，曲线升阶不改变曲线的形状。因此，我们需要$n$个新的控制点，假设新的控制点集为 $Q_0, Q_1, Q_2, ..., Q_{n+1}$，根据上述要求有，$Q_0 = P_0$ 和 $Q_{n+1} = P_n$。其他控制点的计算方式如下：

$$\mathbf{Q}_i=\frac{i}{n+1}\mathbf{P}_{i-1}+\left(1-\frac{i}{n+1}\right)\mathbf{P}_i\quad1\le i\le n$$

以下是从$Q_1$到$Q_n$每个控制点的公式：

$$\begin{aligned} \mathbf{Q}_1&=\frac{1}{n+1}\mathbf{P}_0+\left(1-\frac{1}{n+1}\right)\mathbf{P}_1\\ \mathbf{Q}_2&=\frac{2}{n+1}\mathbf{P}_1+\left(1-\frac{2}{n+1}\right)\mathbf{P}_2\\ \mathbf{Q}_3&=\frac{3}{n+1}\mathbf{P}_2+\left(1-\frac{3}{n+1}\right)\mathbf{P}_3\\ &\vdots\\ \mathbf{Q}_{n-1}&=\frac{n-1}{n+1}\mathbf{P}_{n-2}+\left(1-\frac{n-1}{n+1}\right)\mathbf{P}_{n-1}\\ \mathbf{Q}_n&=\frac{n}{n+1}\mathbf{P}_{n-1}+\left(1-\frac{n}{n+1}\right)\mathbf{P}_n \end{aligned}$$

原始多段线的每一段上正好会有一个新的控制点。即：线段$P_{i-1}P_i$ 包含新控制点$Q_i$。

在de Casteljau算法中，线段$AB$上一点$C$，将$AB$按比例 $u:1-u$ 分割，可以写作 $C = (1 - u)A + uB$。在上式中，$Q_i$将线段$P_{i-1}P_i$ 按比例 $1-\frac{i}{n+1}:\frac{i}{n+1}$ 分割。这个比例不是一个常数，而是随着 i 的值变化的，但与de Casteljau算法很相似

获得了新的控制点集后，原控制点集就可以弃用。由于原始控制多段线的每一段都包含一个新控制点，用新控制点集合替换旧控制点集合的过程可以看作是在原始控制点处切掉角。

见下图，图中显示了一个4度贝塞尔曲线，控制点以红色矩形显示，控制多段线以虚线显示。将其度数增加到5后，新的控制多段线以实线段显示。显然，所有角都被切掉了。下表给出了原始控制多段线上每一段的比例。

i	1 - i/(n+1)	i/(n+1)
1	0.8	0.2
2	0.6	0.4
3	0.4	0.6
4	0.2	0.8

随着曲线度数的增加，控制点的数量也会增加。并且，新的控制多边形会逐渐向曲线靠拢。在下面的图形中，我们从一个度数为6、控制点为7个的贝塞尔曲线开始，然后，其度数增加到7、8、10、15和29。从图中可以看出，随着度数的增加，曲线的形状没有改变，控制多边形越来越接近曲线。最终，随着度数不断增加到无限，控制多边形将趋向于曲线，曲线则相当于控制多边形的一个极限位置。

总结

总的来说，要定义一个$n$阶的贝塞尔曲线，我们需要在空间中选择$n+1$个控制点，这些点大致描述了所需曲线的形状。

然后，如果它不符合我们的期望，我们可以移动控制点。随着一个或多个控制点的移动，贝塞尔曲线的形状相应地发生变化。但是，曲线始终位于由控制点定义的凸包内（凸包性质），并且生成的曲线形状比控制多边形简单（变差减少性质）。

除此之外，贝塞尔曲线还具有非负性、单位分割、仿射不变性、递归性等性质。贝塞尔曲线不具备局部修改性，移动控制点将会全面改变贝塞尔曲线的形状，这也是其逐渐被样条曲线和Nurbs曲线取代的原因。

详细介绍了De Casteljau算法，用于计算任意参数位置时，对应的贝塞尔曲线上一点，这是一种自底向上的三角形计算方案，能够避免自顶向下计算方案递归调用时带来的不必要的性能开销。

还介绍了贝塞尔曲线导数的性质，如何将两个贝塞尔曲线连接在一起并保证连接点处$C^1$连续，如何构造闭合的贝塞尔曲线以及计算高阶贝塞尔导数*。

最后还介绍了贝塞尔曲线的拆分和升阶方法，这个早期在工程中也是具有很大的应用。

文献引用

• Gerald Farin, Curves and Surfaces for CAGD: A Practical Guide, fourth edition, Academic Press, 1997.
• Josef Hoschek and Dieter Lasser, Fundamentals of Computer Aided Geometric Design, translated from the German 1989 edition by Larry L. Schumaker, A K Peters, 1993.
• Les Piegl and Wayne Tiller, The NURBS Book, second edition, Springer-Verlag, 1997.
• David F. Rogers and J. Alan Adams, Mathematical Elements for Computer Graphics, second edition, McGraw-Hill, 1990.

---

后续，将会继续介绍B样条曲线和Nurbs曲线，这两种类型的曲线会逐渐取代贝塞尔成为工程中的主流曲线定义方式。To be continued…