title: Bezier Cubic Curves: moving with uniform acceleration
body: <p>Let's say I have a <a href="http://en.wikipedia.org/wiki/B%C3%A9zier_curve" rel="noreferrer">Bezier curve</a> <code>B(u)</code>, if I increment <code>u</code> parameter at a constant rate I don't obtain a costant speed movement along the curve, because the relation between <code>u</code> parameter and the point obtained evaluating the curve is not linear.</p>
  
  <p>I've read and implemented a David Eberly's <a href="http://www.geometrictools.com/Documentation/MovingAlongCurveSpecifiedSpeed.pdf" rel="noreferrer">article</a>. It explains how to move at constant speed along a parametric curve. </p>
  
  <p>Suppose I have a function <code>F(t)</code> that takes as input a time value <code>t</code> and a speed function <code>sigma</code> that returns the speed value at time <code>t</code>, I can obtain a costant speed movement along the curve, varying t parameter at a constant rate: <code>B(F(t))</code></p>
  
  <p>The core of the article I'm using is the following function: </p>
  
  <pre><code>float umin, umax; // The curve parameter interval [umin,umax].
  Point Y (float u); // The position Y(u), umin &lt;= u &lt;= umax.
  Point DY (float u); // The derivative dY(u)/du, umin &lt;= u &lt;= 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 &lt;= t &lt;= 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 &lt;= 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;
  }
  </code></pre>
  
  <p>It allows me to get the curve parameter <code>u</code> calculated using the supplied time <code>t</code> and sigma function. 
  Now the function works fine when the speed sigma is costant. If sigma is represent a uniform accelartion I'm getting wrong values from it.</p>
  
  <p>Here's an example, of a straight Bezier curve, where P0 and P1 are the control points, T0 T1 the tangent. The curve's defined:</p>
  
  <pre><code>[x,y,z]= B(u) =(1–u)3P0 + 3(1–u)2uT0 + 3(1–u)u2T1 + u3P2 
  </code></pre>
  
  <p><img src="https://i.stack.imgur.com/BzPgK.png" alt="enter image description here"></p>
  
  <p>Let's say I want to know the position along the curve at time <code>t = 3</code>.
  If I a constant velocity:</p>
  
  <pre><code>float sigma(float t)
  {
    return 1f;
  }
  </code></pre>
  
  <p>and the following data:</p>
  
  <pre><code>V0 = 1;
  V1 = 1;
  t0 = 0;
  L = 10;
  </code></pre>
  
  <p>I can analitically calculate the position:</p>
  
  <pre><code>px = v0 * t = 1 * 3 = 3
  </code></pre>
  
  <p>If I solve the same equation using my Bezier spline and the algorithm above with <code>n =5</code> I obtain:</p>
  
  <pre><code>px = 3.002595;
  </code></pre>
  
  <p>Considering numerically approximation the value is quite precise (I did a lot of test on that. I omit details but Bezier my curves implementation is fine and the length of the curve itself is calculated quite precisely using <a href="http://en.wikipedia.org/wiki/Gaussian_quadrature" rel="noreferrer">Gaussian Quadrature</a>).</p>
  
  <p>Now If I try to define sigma as a uniform acceleration function, I get bad results.
  Consider the following data:</p>
  
  <pre><code>V0 = 1;
  V1 = 2;
  t0 = 0;
  L = 10;
  </code></pre>
  
  <p>I can calculate the time a particle will reach the P1 using linear motion equations:</p>
  
  <pre><code>L = 0.5 * (V0 + V1) * t1 =&gt;
  t1 = 2 * L / (V1 + V0) = 2 * 10 / 3 = 6.6666666
  </code></pre>
  
  <p>Having <code>t</code> I can calculate acceleration:</p>
  
  <pre><code>a = (V1 - V0) / (t1 - t0) = (2 - 1) / 6.6666666  = 0.15
  </code></pre>
  
  <p>I have all data to define my sigma function:</p>
  
  <pre><code>float sigma (float t)
  {
    float speed = V0 + a * t;
  }
  </code></pre>
  
  <p>If I analitically solve this I'd expect the following velocity of a particle after time <code>t =3</code> :</p>
  
  <pre><code>Vx = V0 + a * t = 1 + 0.15 * 3 = 1.45
  </code></pre>
  
  <p>and the position will be:</p>
  
  <pre><code>px = 0.5 * (V0 + Vx) * t = 0.5 * (1 + 1.45) * 3 = 3.675
  </code></pre>
  
  <p>But if I calculate it with the alorithm above, the position results:</p>
  
  <pre><code>px = 4.358587
  </code></pre>
  
  <p>that is quite different from what I'm expecting.</p>
  
  <p>Sorry for the long post, if anyone has enough patience to read it, I'd be glad.</p>
  
  <p>Do you have any suggestion?What am I missing?Anyone can tell me what I'm doing wrong?</p>
  
  <hr>
  
  <p>EDIT:
  I'm trying with 3D Bezier curve. Defined this way:</p>
  
  <pre><code>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;
  }
  </code></pre>
  
  <p>and the derivative:</p>
  
  <pre><code>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;
  }
  </code></pre>