贝塞尔曲线运动n阶追踪方程的数学原理及其匀速化方法和应用

前言

  首先声明一下本文讨论的终极目标——通过贝塞尔曲线实现可调控的数值缓动，也就是贝塞尔插值。之所以以数值为目的是因为它相较于实现某一个具体的案例而言，意义更为广泛，例如可以实现可控性很强的缓动动效，在本文的末尾会附上具体的贝塞尔运动案例分析与一些应用层面的介绍。文章的前一部分即推导部分会以较多的数学形式来加以表述，而中间关于贝塞尔运动的匀速化会则尽量用通俗易懂的语言来说明，最后附上一些相关代码以及应用和总结。
  其实这方面已经有几个月没有再接触了，之所以写这篇文章，是为了回顾与加深之前在这方面学习研究的印象，文章内容仅是个人的理解与一些思考和感悟，并不一定是绝对正确的，因此要是文章内容有什么错误或歧义，还得请各位点醒指出，必当虚心学习，不胜感激。

法国数学家Pierre Bézier

一、关于贝塞尔曲线 & 贝塞尔运动

1.简介

  贝塞尔曲线的介绍在网上到处都能找到，这里就不长篇大论地引入了。简而言之，贝塞尔曲线就是这样的一条曲线，它是依据四个位置任意的点坐标绘制出的一条光滑曲线。而贝塞尔追踪方程则是以点运动的形式描述了贝塞尔曲线的形成过程，将曲线描述为了一条随连续时间而形成的点迹。

2.作用

  由于用计算机画图大部分时间是操作鼠标来掌握线条的路径，与手绘的感觉和效果有很大的差别。即使是一位精明的画师能轻松绘出各种图形，拿到鼠标想随心所欲的画图也不是一件容易的事。这一点是计算机万万不能代替手工的工作，所以到目前为止人们只能颇感无奈。使用贝塞尔工具画图很大程度上弥补了这一缺憾。贝塞尔曲线是计算机图形图像造型的基本工具，是图形造型运用得最多的基本线条之一。

  它通过控制曲线上的四个点（起始点、终止点以及两个相互分离的中间点）来创造、编辑图形。其中起重要作用的是位于曲线中央的控制线。这条线是虚拟的，中间与贝塞尔曲线交叉，两端是控制端点。移动两端的端点时贝塞尔曲线改变曲线的曲率（弯曲的程度）；移动中间点（也就是移动虚拟的控制线）时，贝塞尔曲线在起始点和终止点锁定的情况下做均匀移动。注意，贝塞尔曲线上的所有控制点、节点均可编辑。这种“智能化”的矢量线条为艺术家提供了一种理想的图形编辑与创造的工具。

2.贝塞尔运动

  简单地说，贝塞尔运动就是按照贝塞尔曲线所描述的规律来完成一段运动。后文有详细案例，此处就先提上这么一句。

二、贝塞尔追踪方程的数学推导

1.一阶推导

  虽然说是依据四个点而绘制出的曲线，但实际上最少2个点就可以构成一条贝塞尔曲线，这样的贝塞尔曲线又叫一阶的贝塞尔曲线，也是无控制点的贝塞尔曲线。
  前面说到追踪方程描述了点迹，一般地我们将迹点运动的总时间记为1，并把这两个端点分别记作 P₀ 和 P₁ ，运动的迹点记作 M₀ ，就有 t 时刻对应 M₀ 的位置为 P₀ + (P₁ - P₀)t，将其写成方程形式：
$$\begin{array}{rl}{B}_1(t)& =P_0+(P_1-P_0)t\\ & =(1-t)P_0+t{P}_1\end{array}$$

2.二阶推导

  二阶在一阶的基础上增加一个点 P₂ ，此时不能够直接得出迹点的位置。 P₀ 和 P₁ 做一次一阶得到 M₀ 的位置，P₁ 和 P₂ 做一次一阶得到了 M₁ 的位置，此时出现了两个临时迹点 M₀ 和 M₁ ，需要再对 M₀ 和 M₁ 求一次一阶来得到真实迹点 M^/₀ 。即：
$$\begin{array}{rl}{B}_2(t)& =M_0+(M_1-M_0)t\\ & =P_0+(P_1-P_0)+[P_1+(P_2-P_1)t-P_0-(P_1-P_0)t]t\\ & =P_0+2(P_1-P_0)t+(P_0-2P_1+P_2)t^2\\ & =(1-t{)}^2{P}_0+2t(1-t)P_1+t^2{P}_2\end{array}$$

3.三阶推导

  仿照一阶到二阶的升阶过程，不难得到三阶推导过程：
B 2 ( t ) = M 0 ′ + ( M 1 ′ − M 0 ′ ) t = M 0 + ( M 1 − M 0 ) + [ M 1 + ( M 2 − M 1 ) t − M 0 − ( M 1 − M 0 ) t ] t = M 0 + 2 ( M 1 − M 0 ) t + ( M 0 − 2 M 1 + M 2 ) t 2 = P 0 + ( P 1 − P 0 ) t + 2 [ P 1 + ( P 2 − P 1 ) t − P 0 − ( P 1 − P 0 ) t ] t + { P 0 + ( P 1 − P 0 ) t − 2 [ P 1 + ( P 2 − P 1 ) t ] + P 2 + ( P 3 − P 2 ) t } t 2 = ( 1 − t ) 3 P 0 + 3 t ( 1 − t ) 2 P 1 + 3 t 2 ( 1 − t ) P 2 + t 3 P 3 \begin {aligned} B_{2}(t)&=M^\prime_0+(M^\prime_{1}-M^\prime_{0})t \\ &=M_0+(M_1-M_0)+[M_1+(M_2-M_1)t-M_0-(M_1-M_0)t]t \\ &=M_0+2(M_1-M_0)t+(M_0-2M_1+M_2)t^2 \\ &=P_0+(P_1-P_0)t+2[P_1+(P_2-P_1)t-P_0-(P_1-P_0)t]t+\{P_0+(P_1-P_0)t-2[P_1+(P_2-P_1)t]+P_2+(P_3-P_2)t\}t^2 \\ &=(1-t)^3P_0+3t(1-t)^2P_1+3t^2(1-t)P_2+t^3P_3 \\ \end {aligned} B2​(t)​=M0′​+(M1′​−M0′​)t=M0​+(M1​−M0​)+[M1​+(M2​−M1​)t−M0​−(M1​−M0​)t]t=M0​+2(M1​−M0​)t+(M0​−2M1​+M2​)t2=P0​+(P1​−P0​)t+2[P1​+(P2​−P1​)t−P0​−(P1​−P0​)t]t+{P0​+(P1​−P0​)t−2[P1​+(P2​−P1​)t]+P2​+(P3​−P2​)t}t2=(1−t)3P0​+3t(1−t)2P1​+3t2(1−t)P2​+t3P3​​

4.n阶

  综合上式列表，不难总结出一般性的规律。

阶数(n)	Bn(t)
1	B1(t) = (1-t)P0 + tP1
2	B2(t) = (1-t)2P0 + 2t(1-t)P1 + t2P2
3	B3(t) = (1-t)3P0 + 3t(1-t)2P1 + 3t2(1-t)P2 + t3P3
…	…

$$\begin{array}{rl}{B}_n(t)& =(1-t{)}^n{P}_0+\sum _{i=1}^{n-1}{C}_n^i{t}^i(1-t{)}^{n-i}{P}_i+t^n{P}_n(n\ge 2,n\in Z)\end{array}$$

三、匀速化处理

  前面介绍了贝塞尔曲线的数学表述形式，利用之我们可以根据任意个点的坐标画出曲线的直观路径，也可以用一组点来很好地拟合某一给定的曲线。现在，除了画来看之外，我们要用它来做点别的事情。什么事情？请看下图。

  可以看到，上图小球的运动是非常死板的，其位置的变动与时间成线性关系。如果我们用图像描述上述运动，考虑到时间是均匀变化的，小球的位置也是均匀变化的，因此图像对应的是一条线段。而每一时刻 t 对应线段上的点处的导数值就是此时刻的瞬时速率，因为直线的斜率是固定的，所以小球起始和末了都是同一个速度，导致运动开始的很突兀，结束的也很突兀。

  由此不难想到，如果我们想要使小球具有更加平滑的运动效果的话，只要对曲线稍加变形，让其更加符合让我们的动态审美。例如将上图的直线变形为下图。

  这样小球的速度就能够平缓地增大到一定数值进行运动，最后也是平滑地停止下来了。

  但是请注意，上图平滑的运动效果仅仅是我们想象出来的，因为最开始 position-t 图像是根据运动而画出来的，之后我们人为地修改了曲线，但并没有让原来的小球按照我们改过的曲线来运动。
  那么问题来了，如何让我们的小球按照我们绘制的曲线来运动呢？请继续往下看。

1.什么是匀速化 & 为什么要匀速化

  有些朋友看到这里可能按捺不住了，这难道还不简单吗？你刚才不是推导出了 n 阶的迹点追踪方程吗？均匀地传入 t 参数然后算出对应地position值赋给小球的位置不就完事了吗？
  然而实时并没有想象的那么轻松，回过头来看看我们之前推导出来的方程：
$$\begin{array}{rl}{B}_n(t)& =(1-t{)}^n{P}_0+\sum _{i=1}^{n-1}{C}_n^i{t}^i(1-t{)}^{n-i}{P}_i+t^n{P}_n(n\ge 2,n\in Z)\end{array}$$
  我们注意到，追踪方程的右侧并不是单纯的横坐标或是纵坐标，而是直接代以点的形式来表述。这表明了在实际计算的时候，需要按照需要的轴来代入点的坐标进行计算，例如需要计算以 P₀、P₁、P₂ 三个点确立的二阶贝塞尔曲线在 t=0.5 时的迹点的横坐标时，使用 P₀.x、P₁.x、P₂.x 代入B₂(t)的表达式即可计算输出对应的横坐标，纵坐标同理计算。
  也就是说代入均匀的时间 t ，输出的并不是我们想要的纵坐标，因为贝塞尔曲线描述的图像是在 y-x 图像下的，而我们的之前修改的运动图是 position-t 图，也就是 y-t 图，两者虽然看上去长的一一样，但不在同一个坐标系下，所具有的意义也就相去甚远。而按照这种不完全的映射关得到的所谓的 y-t 运动图像，我们称之为非匀速贝塞尔运动。
  那么究竟是什么原因导致了非匀速，或者说哪一步我们理解上容易出问题进而导致了错误的效果？究其因还是轴的偷换。下面我们不说具体的轴的名称了，就说横轴和纵轴。我们想要的结果是，横轴上的数匀速前进，向上投射到曲线上的某点，再映射到纵轴得到 position 坐标。

  可是贝塞尔曲线所在坐标系的横坐标是 x，并非 t！因此将匀速的 t 代入 B_n(t) 是不能实现我们要的效果的，我们需要代入匀速的 x。由于当 t 匀速时， x 本身是由 B_n(t) 算得的非匀速值，而通过匀速 t 对应非匀速的 x 求出其在匀速状态下应当对应的 t，就是匀速化的过程。

2.两种匀速化的办法

  实际操作的时候，我们将匀速变化的时间也就是正常流动的时间记为 t，虽说是匀速的时间但我们心里明白这个 t 对我们的需求来说是“非匀速”的，将 t 当作 y-x 系下的 x 坐标代入 B_n(t)，反解出其对应的真实的 r_t (即 real-t) 是多少，进而把 r_t 代回 B_n(t) 最终得到了我们要的 y 值，即 position 值。但想要反解 t 也不是那么轻松，毕竟 B_n(t) 的表达式有那么复杂。因此在这里匀速化的核心问题就是如何根据 P.x 反解 t。

3.方法一：利用辛普森积分法匀速化

  此方法以画曲为直为思路，算出贝塞尔曲线的近似长度，将其放倒为一条线段。水平轴的数值匀速变大，即可看作是迹点在线上匀速运动，又时间取做单位1，自然不难求出点运动的速度v，而后与正常传入的 t 相乘得到匀速运动的距离L。最后用二分法或黄金分割法求出一个 t’ 时刻使得与 t’ 对应的 0~x 范围内的曲线长度L’与 L 近似相等，此时所求得的 t 即是我们要找的realTime。
  此处给出了思路，下面附上求三阶贝塞尔曲线近似长度的算法，剩下的算法请感兴趣的读者自行完成。

// 求 0~t 段的贝塞尔曲线长度
double beze_length(Coordinate p[], double t)
{
	int TOTAL_SIMPSON_STEP;
	int stepCounts;
	int halfCounts;
	int i;
	double sum1, sum2, dStep;
	
	TOTAL_SIMPSON_STEP = 100000;
	stepCounts = TOTAL_SIMPSON_STEP*t;

	if(stepCounts == 0)
	{
		return 0;
	}
	if(stepCounts%2 == 0)
	{
		stepCounts++;
	}

	halfCounts = stepCounts/2;
	dStep = t / stepCounts;

	while(i<halfCounts)
	{
		sum1 += beze_speed (p, (2 * i + 1) * dStep);
		i++;
	}
	i=1;
	while(i<halfCounts)
	{
		sum2 += beze_speed (p, 2 * i * dStep);
		i++;
	}
	return ((beze_speed (p, 0) + beze_speed (p, 1) + 2 * sum2 + 4 * sum1) * dStep / 3);
}

double beze_speed_x (Coordinate p[], double t)
{
	// 三阶
	return -3 * p [0].x * pow (1 - t, 2) + 3 * p [1].x * pow (1 - t, 2) - 6 * p [1].x * (1 － t) * t ＋ 6 * p [2].x * (1 - t) * t - 3 * p [2].x * pow (t, 2) + 3 * p [3].x * pow (t, 2);
}

double beze_speed_y (Coordinate p[], double t)
{
	// 三阶
	return -3 * p [0].y * pow (1 - t, 2) + 3 * p [1].y * pow (1 - t, 2) - 6 * p [1].y * (1 － t) * t ＋ 6 * p [2].y * (1 - t) * t - 3 * p [2].y * pow (t, 2) + 3 * p [3].y * pow (t, 2);
}

//	求合速度
double beze_speed (Coordinate p[], double t)
{
	double vx = beze_speed_x (p, t);
	double vy = beze_speed_y (p, t);
	return sqrt(pow (vx, 2) + pow (vy, 2));
}

4.方法二：编程逼值实现匀速化

  此方法就是简单地不断二分来求 B₃(t) 的反函数值，与方法一相比较而言更为简单也更易于理解，下面附上完整代码。

#include <stdio.h>
#include <math.h>
struct Coordinate
{
	int x;
	int y;
};

// 声明函数
int b3(Coordinate p[], double t);
double t2rt(Coordinate p[], double t);

int main(void)
{
	// 定义四个贝塞尔点
	Coordinate p[4]={{100, 611}, {300, 411} ,{400, 311} ,{500, 411}};
	
	// 输出 t=0.5 对应的匀速 rt
	printf("%lf\n", t2rt(p, 0.5));
	
	return 0;
}

int b3(Coordinate p[], double t)
{
	// 三阶贝塞尔运算
	int retn;
	retn = pow(1-t, 3) * p[0].x + 3*t*pow(1-t, 2) * p[1].x + 3*pow(t, 2) * (1-t)*p[2].x + pow(t, 3)*p[3].x;
	return retn;
}

double t2rt(Coordinate p[], double t)
{
	// 定义真实时间与差时变量
	double realTime, deltaTime;
	
	// 曲线上的 x 坐标
	int bezierX;
	
	// 计算 t 对应曲线上匀速的 x 坐标
	int x = 100 + (p[3].x - p[0].x)*t;
	
	realTime = 1;
	do
	{
		// 半分
		if(deltaTime>0)
		{
			realTime -= (double)realTime/2;
		}
		else
		{
			realTime += (double)realTime/2;
		}

		// 计算本此 "rt" 对应的曲线上的 x 坐标
		bezierX = b3(p, realTime);
		
		// 计算差时
		deltaTime = bezierX - x;
	}
	// 差时逼近为0时跳出循环
	while(deltaTime != 0);

	return realTime;
}

  以上就是两种将贝塞尔曲线运动匀速化的方法，下面我们来看两个简单的应用。

四、匀速化后的应用

1.高自由度的缓动动效
  在将贝塞尔曲线应用于运动之前，一些动画的过度效果往往比较生硬和死板（如上面提到的线性运动）。而贝塞尔曲线运动则可以得到很好的视觉效果。一些网站提供了缓动函数表，使用这些已经封装好公开使用的缓动函数可以方便快捷地将一些流畅的动画效果引入到自己的项目当中实现一些简单的缓动动效。
  但是其所提供的缓动函数毕竟有限，在此之外的效果想要实现可能就麻烦了，例如博主曾经在其他地方推导介绍过利用
$$\begin{array}{rl}k(t)& =\frac{1}{2}[sin(\pi t-\frac{\pi }{2})+1]\end{array}$$
来简单地正弦软化线性数据，而应用本文所讲述的贝塞尔曲线运动则可以实现高自由度的可调整缓动效果，借此可实现动画效果与程序本体的逻辑分离，提高交互动效的开发自由度，也更加便于后期更新和维护。
（具体实例待更新）
2.物理效果的模拟
  考虑到某一物体(图形)在任意轴方向上的运动都可以通过贝塞尔曲线来方便地调节与控制，自然可以将其应用与物理效果的模拟上。例如令一小圆在水平方向上随时间匀速运动，而竖直方向上其y坐标的值按照二次形状的匀速化贝塞尔曲线来增大，这样就模拟了平抛运动，而通过自由地调节贝塞尔曲线的形状，可以非常轻松地模拟其他更加复杂的物理效果。
（具体实例待更新）

五、总结

待更新