匀速二阶贝塞尔曲线

关于二阶贝塞尔曲线匀速运动的实现，网上有很多文章介绍：

匀速贝塞尔曲线运动的实现(一) | 我的博客和笔记

匀速贝塞尔曲线运动算法-CSDN博客

How to achieve uniform speed of movement on a bezier curve?

但每次艰难理解后，下次再看时又得重头再来，所以这次准备把自己的理解过程记录下来留待后观。

二阶贝塞尔曲线长度

给定3个点(P_0) (P_1) (P_2)，用(B(t))表示二阶贝塞尔曲线：

[begin{align} B(t)&=(1-t)^2P_0+2t(1-t)P_1+t^2P_2, &tin[0,1] \ begin{bmatrix} x(t) \ y(t) \ end{bmatrix} &=(1-t)^2 begin{bmatrix} x_0 \ y_0 \ end{bmatrix} +2t(1-t) begin{bmatrix} x_1 \ y_1 \ end{bmatrix} +t^2 begin{bmatrix} x_2 \ y_2 \ end{bmatrix} , &tin[0,1] end{align} ]

求二阶贝塞尔曲线相对于(t)的速度：

[begin{align} V(t)&=B'(t)\ &= begin{bmatrix} x'(t) \ y'(t) \ end{bmatrix}\ &=-2(1-t) begin{bmatrix} x_0 \ y_0 \ end{bmatrix} +2(1-2t) begin{bmatrix} x_1 \ y_1 \ end{bmatrix} +2t begin{bmatrix} x_2 \ y_2 \ end{bmatrix} , &tin[0,1] end{align} ]

速度(V(t))是一个二维向量，计算曲线长度不需要速度的方向，因此取其标量：

[begin{align} s(t) &= Vert B'(t) Vert \ &= sqrt{{x'(t)}^2 + {y'(t)}^2} \ end{align} ]

其中(x'(t))可以进一步变换为：

[begin{align} x'(t)&=-2x_0+2x_0t+2x_1-4x_1t+2x_2t \ &=(2x_0-4x_1+2x_2)t-2x_0+2x_1 end{align} ]

同理(y'(t))可以进一步变换为：

[begin{align} y'(t)&=-2y_0+2y_0t+2y_1-4y_1t+2y_2t \ &=(2y_0-4y_1+2y_2)t-2y_0+2y_1 end{align} ]

分别用(a_x) (b_x) (a_y) (b_y) 替换(x'(t))和(y'(t))中复杂的部分：

[begin{align} x'(t)&=a_xt+b_x \ a_x&=2x_0-4x_1+2x_2 \ b_x&=-2x_0+2x_1 \ end{align} ]

[begin{align} y'(t)&=a_yt+b_y \ a_y&=2y_0-4y_1+2y_2 \ b_y&=-2y_0+2y_1 \ end{align} ]

那么(s(t))可以进一步简化为：

[begin{align} s(t) &= sqrt{{x'(t)}^2 + {y'(t)}^2} \ &= sqrt{a_x^2t^2+2a_xb_xt+b_x^2+a_y^2t^2+2a_yb_yt+b_y^2} \ &= sqrt{(a_x^2+a_y^2)t^2+(2a_xb_x+2a_yb_y)t+b_x^2+b_y^2} end{align} ]

再一次用(A) (B) (C)代替(s(t))中复杂的部分：

[begin{align} s(t) &= sqrt{At^2+Bt+C} \ A &= a_x^2+a_y^2 \ B &= 2a_xb_x+2a_yb_y \ C &= b_x^2+b_y^2 end{align} ]

这里的(A) (B) (C)和[1]和[2]中的不太一样，但最终结果是一样的。对速度(s(t))进行积分可求得距离，可以使用积分计算器进行计算：

[begin{align} int s(t) &= int sqrt{At^2+Bt+C} \ &= frac{1}{8A^frac{3}{2}} big(2sqrt{A}(B+2At)sqrt{C+Bt+At^2}-(B^2-4AC)ln(B+2At+2sqrt{A}sqrt{C+Bt+At^2})big) end{align} ]

由于(t)的取值范围是0到1，因此二阶贝塞尔曲线的长度公式应该写为：

[begin{align} L(t) =& int_0^t s(t) \ =& int s(t) - int s(t) |_{t=0} \ =& frac{1}{8A^frac{3}{2}} big(2sqrt{A}(B+2At)sqrt{C+Bt+At^2}-(B^2-4AC)ln(B+2At+2sqrt{A}sqrt{C+Bt+At^2})big) \ &- frac{1}{8A^frac{3}{2}} big(2Bsqrt{AC}-(B^2-4AC)ln(B+2sqrt{AC})big) \ =& frac{1}{8A^frac{3}{2}} big(2sqrt{A}(B+2At)sqrt{C+Bt+At^2}-(B^2-4AC)ln(B+2At+2sqrt{A}sqrt{C+Bt+At^2}) \ &- 2Bsqrt{AC}+(B^2-4AC)ln(B+2sqrt{AC}) big) \ =& frac{1}{8A^frac{3}{2}} big(T_1T_0-T_2ln(T_1+T_0) \ &- BT_3+T_2ln(B+T_3) big) \ end{align} ]

其中(T_0) (T_1) (T_2)表示为：

[begin{align} T_0 =& 2sqrt{A}sqrt{C+Bt+At^2} \ T_1 =& B+2At \ T_2 =& B^2-4AC \ T_3 =& 2sqrt{AC} end{align} ]

为了避免括号层级太多，这里没有进一步提取公因子，最终保持两层括号，同时将替换复杂重复的算式替换为变量。((31))和((32))进行公因子提取变换后和[1][2]的结果是一致的。

匀速二阶贝塞尔曲线

现在已经有了二阶贝塞尔曲线的长度公式，要沿贝塞尔曲线进行运动，也就是每次移动的距离相等，首先(L(1.0))为贝塞尔曲线的总长度，如果移动(N)次，那么每次移动的距离为(frac{1}{N}L(1.0))。如果已知第(n)次的等距移动的(t)表示为(t_n)，那么(t_{n+1})可以表示为：

[begin{align} t_{n}=t_{n-1}+frac{frac{1}{N}L(1.0)}{s(t_{n-1})} quad,quad nin[1, N] end{align} ]

其中(s(t_n))为(t_n)处的速度，那么移动(frac{1}{N}L(1.0))的距离所需的时间为(frac{frac{1}{N}L(1.0)}{s(t_n)})，当然这种计算方式仅能求得近似值，这也是[3]给出的方法。如果仔细看[1][2]的代码实现对上式进行了修正：

[begin{align} t_n=t_{n-1}+frac{frac{n}{N}L(1.0)-L(t_{n-1})}{s(t_{n-1})} quad,quad nin[1, N] end{align} ]

但最终求得的解仍然是近似解。

Python实现

最终匀速贝塞尔曲线的Python实现代码如下：

import numpy as np import matplotlib.pyplot as plt import sys   def bezier_length(A, B, C, t):  # 二阶贝塞尔曲线长度     T0 = 2 * np.sqrt(A) * np.sqrt(C + B * t + A * t * t)     T1 = B + 2 * A * t     T2 = B * B - 4 * A * C     T3 = 2 * np.sqrt(A * C)     return 1 / (8 * np.pow(A, 1.5)) * (T0 * T1 - T2 * np.log(T1 + T0) - B * T3 + T2 * np.log(B + T3))   def bezier_velocity(A, B, C, t):  # 二阶贝塞尔曲线速度     return np.sqrt(C + B * t + A * t * t)   def bezier_movement(p0, p1, p2, N):  # 二阶贝塞尔曲线等距移动     a = 2 * p0 - 4 * p1 + 2 * p2     b = -2 * p0 + 2 * p1     A = a[0] * a[0] + a[1] * a[1]     B = 2 * a[0] * b[0] + 2 * a[1] * b[1]     C = b[0] * b[0] + b[1] * b[1]     L1 = bezier_length(A, B, C, 1)     tn = 0     res = np.zeros((N + 1, 2))     for n in range(N + 1): # 注意这里n的取值范围是[0, N]         res[n] = (1 - tn) * (1 - tn) * p0 + 2 * tn * (1 - tn) * p1 + tn * tn * p2         tn = tn + ((n + 1) / N * L1 - bezier_length(A, B, C, tn)) / bezier_velocity(A, B, C, tn)     return res   P0 = np.array([0, 0])  # 起点 P1 = np.array([3, 8])  # 控制点 P2 = np.array([6, 0])  # 终点  N = 10 if len(sys.argv) > 1:     N = int(sys.argv[1])  curve = bezier_movement(P0, P1, P2, N)  plt.figure(figsize=(10, 6))  plt.plot([P0[0], P1[0], P2[0]], [P0[1], P1[1], P2[1]], 'ro--', label='Control-Line') plt.plot(curve[:, 0], curve[:, 1], 'o', linewidth=2, label='Bezier movement')  plt.text(P0[0], P0[1], 'P0', fontsize=12, ha='right', va='top') plt.text(P1[0], P1[1], 'P1', fontsize=12, ha='center', va='bottom') plt.text(P2[0], P2[1], 'P2', fontsize=12, ha='left', va='top')  plt.title('Bezier Curve(N=' + str(N) + ')', fontsize=14) plt.xlabel('X-axis', fontsize=12) plt.ylabel('Y-axis', fontsize=12) plt.grid(True, linestyle='--', alpha=0.7) plt.legend() plt.axis('equal') plt.show() 

当(N=10)，(N=20)，(N=30)的运行结果如下：

在(N=10)时，最后一个点明显与端点不重合（运行Python脚本后，放大后也可观察），因此贝塞尔曲线匀速运动只能求得近似解。在使用中，需要根据实际情况选择合适的(N)，以及如何处理最后一个点。

参考

1. 匀速贝塞尔曲线运动的实现(一) | 我的博客和笔记

2. 匀速贝塞尔曲线运动算法-CSDN博客

3. How to achieve uniform speed of movement on a bezier curve?