2013-03-06 63 views
12

比方說,我有一個Bezier curveB(u),如果我增加以恆定的速度我沒有獲得沿着曲線costant高速運動u參數,因爲u參數和點之間的關係得到的評估曲線不是線性的。貝塞爾三次曲線:以勻加速度移動

我已閱讀並實施了David Eberly的article。它解釋瞭如何沿着參數曲線以恆定速度移動。

假設我有一個功能F(t)其作爲輸入,在時間t返回速度值的時間值t和速度函數sigma,我可以得到沿曲線a costant速度運動,以恆定的速率變化噸參數: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; 
} 

它可以讓我計算出的曲線參數u使用所提供的時間t和西格瑪功能。 速度西格馬是costant現在功能正常工作。如果西格瑪代表統一加速度,我會從中得到錯誤的值。

下面是一個直貝塞爾曲線的例子,其中P0和P1是控制點,T0是切線。曲線的定義:

[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; 

我可以analitically計算位置:

px = v0 * t = 1 * 3 = 3 

如果我解決使用我的貝塞爾樣條相同的方程和上面的算法n =5我得到:

px = 3.002595; 

考慮到數值上的近似值,這個值非常精確(我做了很多測試。我省略了細節,但貝塞爾曲線的實現很好,曲線本身的長度使用Gaussian Quadrature進行精確計算。

現在,如果我嘗試將西格馬定義爲統一的加速功能,我會得到不好的結果。 考慮以下數據:

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 

我有所有數據來定義我的西格馬函數:

float sigma (float t) 
{ 
    float speed = V0 + a * t; 
} 

如果我analitically解決這個問題,我期望一個粒子的速度以下時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 

那是相當不同的fr嗡,我期待。

很抱歉,如果有人有足夠的耐心閱讀它,我會很高興。

你有什麼建議?我錯過了什麼?任何人都可以告訴我我做錯了什麼?


編輯: 我想用3D Bezier曲線。這樣定義的:

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; 
} 
+0

算法給你什麼t = 6.6666 ...?它是10的值,即L還是另一個? – lmsteffan 2013-03-15 22:42:47

回答

1

我的猜測是,n=5根本不給你足夠的精度爲手頭的問題。恆速情況下可以的事實並不意味着恆速情況也是如此。不幸的是,您必須爲自己定義妥協方案,該方案爲符合您的需求和資源的n提供了價值。

無論如何,如果你真的有一個恆定的速度參數化X(U(T))以足夠的精度工作,那麼你可以重新命名這個「時間」參數t「太空」(距離)參數s,所以你真正擁有的是X(s),你只需要插入你需要的s(t):x(s(t))。在你的情況下(恆定加速度),s(t)= s + u t + a t /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 

我已經只是打印出一噸的變量,然後用手計算每個值,並確保它們是相同的調試很多這樣的節目。

+0

感謝您的數據。我會嘗試做更多的實驗。其實我正在使用3D Bezier曲線。我將使用代碼編輯帖子。 – Heisenbug 2013-03-17 12:12:51