贝塞尔三次曲线:匀加速运动

13
假设我有一个贝塞尔曲线B(u),如果我以恒定速率增加u参数,我无法沿着曲线获得恒定速度的移动,因为u参数和通过评估曲线获得的点之间的关系不是线性的。
我已经阅读并实现了David Eberly的文章。它解释了如何沿参数曲线以恒定速度移动。
假设我有一个函数F(t),它将时间值t和速度函数sigma作为输入,后者在时间t处返回速度值,我可以通过以恒定速率变化t参数来沿曲线获得恒定速度的移动:B(F(t)) 我正在使用的文章的核心是以下功能:
float umin, umax; // The curve parameter interval [umin,umax].
Point Y (float u); // The position Y(u), umin <= u <= umax.
Point DY (float u); // The derivative dY(u)/du, umin <= u <= umax.
float LengthDY (float u) { return Length(DY(u)); }
float ArcLength (float t) { return Integral(umin,u,LengthDY()); }
float L = ArcLength(umax); // The total length of the curve.
float tmin, tmax; // The user-specified time interval [tmin,tmax]
float Sigma (float t); // The user-specified speed at time t.

float GetU (float t) // tmin <= t <= tmax
{
  float h = (t - tmin)/n; // step size, `n' is application-specified
  float u = umin; // initial condition
  t = tmin; // initial condition
  for (int i = 1; i <= n; i++)
  {
    // The divisions here might be a problem if the divisors are
    // nearly zero.
    float k1 = h*Sigma(t)/LengthDY(u);
    float k2 = h*Sigma(t + h/2)/LengthDY(u + k1/2);
    float k3 = h*Sigma(t + h/2)/LengthDY(u + k2/2);
    float k4 = h*Sigma(t + h)/LengthDY(u + k3);
    t += h;
    u += (k1 + 2*(k2 + k3) + k4)/6;
  }
  return u;
}

它允许我使用提供的时间t和sigma函数计算曲线参数u。如果sigma代表均匀加速度,那么该函数在速度sigma恒定时运行良好。但是,如果sigma表示均匀加速度,则从中获取错误值。
以下是一个例子,它是一条直线Bezier曲线,其中P0和P1是控制点,T0 T1是切线。该曲线定义如下:
[x,y,z]= B(u) =(1–u)3P0 + 3(1–u)2uT0 + 3(1–u)u2T1 + u3P2 

enter image description here

假设我想知道在时间t = 3时曲线上的位置。 如果我以恒定速度移动:

float sigma(float t)
{
  return 1f;
}

和以下数据:

V0 = 1;
V1 = 1;
t0 = 0;
L = 10;

我可以通过分析计算出位置:

px = v0 * t = 1 * 3 = 3

如果我使用我的贝塞尔样条和上述算法解决相同的方程,并且n=5,我会得到:

px = 3.002595;

考虑数值逼近,该值相当精确(我进行了很多测试。我省略了细节,但是我的贝塞尔曲线实现很好,并且使用高斯积分计算曲线的长度非常准确)。

现在,如果我尝试将sigma定义为均匀加速度函数,则结果不佳。 考虑以下数据:

V0 = 1;
V1 = 2;
t0 = 0;
L = 10;

我可以使用线性运动方程计算粒子到达P1的时间:

L = 0.5 * (V0 + V1) * t1 =>
t1 = 2 * L / (V1 + V0) = 2 * 10 / 3 = 6.6666666

有了t,我可以计算加速度:

a = (V1 - V0) / (t1 - t0) = (2 - 1) / 6.6666666  = 0.15

我有所有数据来定义我的sigma函数:
float sigma (float t)
{
  float speed = V0 + a * t;
}

如果我解析地求解,我期望在时间 t = 3 后粒子的速度如下:
Vx = V0 + a * t = 1 + 0.15 * 3 = 1.45

并且位置将会是:

px = 0.5 * (V0 + Vx) * t = 0.5 * (1 + 1.45) * 3 = 3.675

但如果我使用上述算法进行计算,位置结果如下:

px = 4.358587

这与我所期望的非常不同。

对于这篇长文,如果有人耐心阅读,我会很高兴。

你有什么建议吗?我错过了什么?有人能告诉我我做错了什么吗?


编辑: 我正在尝试使用三维贝塞尔曲线。定义如下:

public Vector3 Bezier(float t)
{
    float a = 1f - t;
    float a_2 = a * a;
    float a_3 = a_2 *a;

    float t_2 = t * t;

    Vector3 point = (P0 * a_3) + (3f * a_2 * t * T0) + (3f * a * t_2 * T1) + t_2 * t * P1 ;

    return point;
}

和导数:

public Vector3 Derivative(float t)
{
    float a = 1f - t;
    float a_2 = a * a;
    float t_2 = t * t;
    float t6 = 6f*t;

    Vector3 der = -3f * a_2 * P0 + 3f * a_2 * T0 - t6 * a * T0 - 3f* t_2 * T1 + t6 * a * T1 + 3f * t_2 * P1;

    return der;
}

算法在t=6.6666...时会给出什么结果呢?它是值10,即L,还是其他值? - lmsteffan
2个回答

1
我的猜测是,n=5 对于手头的问题来说不够精确。它适用于恒定速度情况,并不意味着它也适用于恒定加速度情况。不幸的是,你必须自己定义一个折衷方案,以提供适合你需求和资源的n值。
无论如何,如果你真的有一个足够精确的恒定速度参数化X(u(t)),那么你可以将这个“时间”参数t重命名为“空间”(距离)参数s,这样你实际上拥有的是一个X(s),你只需要插入你需要的s(t):X(s(t))。在你的情况下(恒定加速度),s(t) = s0 + u t + a t2 / 2,其中u和a可以很容易地从你的输入数据中确定。

1
我认为在你的实现中未展示的函数Y和DY中可能只是有一个打字错误。我尝试了一个一维曲线,其中P0=0,T0=1,T1=9,P1=10,并且在n=5时得到了3.6963165,在n=30时改善为3.675044,在n=100时改善为3.6750002。
如果您的实现是二维的,请尝试使用P0=(0,0),T0=(1,0),T1=(9,0)和P1=(10,0)。然后再尝试P0=(0,0),T0=(0,1),T1=(0,9)和P1=(0,10)。
如果您正在使用C,请记住^运算符不表示指数。您必须使用pow(u,3)或u*u*u来获得u的立方。
尝试在每次迭代中打印出尽可能多的值。这是我得到的:
i=1
    h=0.6
    t=0.0
    u=0.0
    LengthDY(u)=3.0
    sigma(t)=1.0
    k1=0.2
    sigma(t+h/2)=1.045
    LengthDY(u+k1/2)=6.78
    k2=0.09247787
    LengthDY(u+k2/2)=4.8522377
    k3=0.12921873
    sigma(t+h)=1.09
    LengthDY(u+k3)=7.7258916
    k4=0.08465043
    t_new=0.6
    u_new=0.12134061
i=2
    h=0.6
    t=0.6
    u=0.12134061
    LengthDY(u)=7.4779167
    sigma(t)=1.09
    k1=0.08745752
    sigma(t+h/2)=1.135
    LengthDY(u+k1/2)=8.788503
    k2=0.0774876
    LengthDY(u+k2/2)=8.64721
    k3=0.078753725
    sigma(t+h)=1.1800001
    LengthDY(u+k3)=9.722377
    k4=0.07282171
    t_new=1.2
    u_new=0.20013426
i=3
    h=0.6
    t=1.2
    u=0.20013426
    LengthDY(u)=9.723383
    sigma(t)=1.1800001
    k1=0.072814174
    sigma(t+h/2)=1.225
    LengthDY(u+k1/2)=10.584761
    k2=0.069439456
    LengthDY(u+k2/2)=10.547299
    k3=0.069686085
    sigma(t+h)=1.27
    LengthDY(u+k3)=11.274727
    k4=0.06758479
    t_new=1.8000001
    u_new=0.26990926
i=4
    h=0.6
    t=1.8000001
    u=0.26990926
    LengthDY(u)=11.276448
    sigma(t)=1.27
    k1=0.06757447
    sigma(t+h/2)=1.315
    LengthDY(u+k1/2)=11.881528
    k2=0.06640561
    LengthDY(u+k2/2)=11.871877
    k3=0.066459596
    sigma(t+h)=1.36
    LengthDY(u+k3)=12.375444
    k4=0.06593703
    t_new=2.4
    u_new=0.3364496
i=5
    h=0.6
    t=2.4
    u=0.3364496
    LengthDY(u)=12.376553
    sigma(t)=1.36
    k1=0.06593113
    sigma(t+h/2)=1.405
    LengthDY(u+k1/2)=12.7838
    k2=0.06594283
    LengthDY(u+k2/2)=12.783864
    k3=0.0659425
    sigma(t+h)=1.45
    LengthDY(u+k3)=13.0998535
    k4=0.06641296
    t_new=3.0
    u_new=0.4024687

我曾经通过打印大量的变量,手动计算每个值并确保它们相同的方式来调试许多类似这样的程序。


谢谢提供数据。我会尝试进行更多的实验。实际上,我正在使用3D贝塞尔曲线。我将编辑帖子并添加代码。 - Heisenbug

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接