Essential Math Weblog

6/7/2005

Reparameterizing a Curve, Part Three

Filed under: Erratical, Mathematical — Jim @ 9:08 pm

Now we come to the final part: merging Newton-Raphson and bisection together, and choosing the correct one depending on the situation.

The merge is fairly simple. First, we unroll the loop. In both cases we are computing a guess, then the function value, and finally determining whether our guess is correct. Then we start all over again. There is some housekeeping involving a speed computation for Newton-Raphson and adjusting the bounds for bisection, but otherwise, the general structure is the same. The main difference is that for Newton-Raphson our first guess is computed outside the loop since it's different from the iterative step, but there's no reason we can't use that initial guess for bisection too. It'll be a lopsided binary search the first time through, but the test value still lies within the interval.

In either case, we'll have to compute the speed and update the interval. We'll need the speed to determine which method to use, and the interval has to be kept up-to-date with our current guess. So with that in mind, our general merged structure looks like:

PLAIN TEXT
C++:
  1.  
  2. // FindParameterByDistance() - first pass bisection-Newton version
  3. float
  4. IvBezier::FindParameterByDistance( float t1, float s )
  5. {
  6.     // set endpoints
  7.     float a = t1;
  8.     float b = mTimes[mCount-1];
  9.    
  10.     // ensure that we remain within valid parameter space
  11.     if ( s >= ArcLength(t1, b) )
  12.         return b;
  13.     if ( s < = 0.0f )
  14.         return a;
  15.    
  16.     // make first guess
  17.     float p = t1 + s*(mTimes[mCount-1]-mTimes[0])/mTotalLength;
  18.    
  19.     // iterate and look for zeros
  20.     for ( u_int i = 0; i < 32; ++i )
  21.     {
  22.         // compute function value and test against zero
  23.         float func = ArcLength(t1, p) - s;
  24.         if ( ::IvAbs(func) < 1.0e-03f )
  25.         {
  26.             return p;
  27.         }
  28.        
  29.          // update endpoints
  30.         if ( func < 0.0f )
  31.         {
  32.             a = p;
  33.         }
  34.         else
  35.         {
  36.             b = p;
  37.         }
  38.        
  39.         // get speed along curve
  40.         float speed = Velocity(p).Length();
  41.        
  42.         // if something tells us we need to use bisection
  43.         if ( bisection )
  44.         {
  45.             // do bisection
  46.             p = 0.5f*(a+b);
  47.         }
  48.         else
  49.         {
  50.             // otherwise Newton-Raphson
  51.             p -= func/speed;
  52.         }
  53.     }
  54.    
  55.     // done iterating, return failure case
  56.     return FLT_MAX;
  57. }
  58.  

So the question remains, what criteria do we use for choosing bisection over Newton-Raphson? Well, our initial problem was that the speed could end up being near-zero, so we could just test for that:

PLAIN TEXT
C++:
  1.  
  2.         // if speed along curve is zero
  3.         if ( ::IsZero( speed ) )
  4.         {
  5.             // do bisection
  6.             p = 0.5f*(a+b);
  7.         }
  8.         else
  9.         {
  10.             // otherwise Newton-Raphson
  11.             p -= func/speed;
  12.         }
  13.  

The problem is that checking for zero is not enough. Even if the speed isn't zero but is sufficiently small, the value func/speed could be quite large -- large enough to cause us to step outside the bisection interval or even the parameter space of the curve. As you might guess, this will lead to some bad results. At the very least, it's unlikely we'll converge to the solution. So this is the case where we want to use bisection instead of Newton's method.

To detect it, we want to find the instances where the new p value lies outside the bisection interval. We'll use that instead of the parameter space of the curve because a) we want to maintain a state so that bisection will work later, if we need it to and b) we know the zero lies in that interval somewhere, and if we narrow our search we'll converge faster. So we want to find the cases where:

PLAIN TEXT
C++:
  1.  
  2. p - func/speed <  a ||  p - func/speed > b
  3.  

We can get rid of the division by multiplying through by speed:

PLAIN TEXT
C++:
  1.  
  2. p*speed - func < a*speed ||  p*speed - func > b*speed
  3.  

Or

PLAIN TEXT
C++:
  1.  
  2. (p-a)*speed < func ||  (p-b)*speed > func
  3.  

This test as it stands isn't very floating point safe, so we can either add a little fuzziness by checking for

PLAIN TEXT
C++:
  1.  
  2. (p-a)*speed - func < kIntervalEps ||  (p-b)*speed - func > -kIntervalEps
  3.  

or the suggestion from Numerical Recipes is

PLAIN TEXT
C++:
  1.  
  2. ((p-a)*speed - func)*((p-b)*speed - func) > 0.0f
  3.  

So our final iterative step becomes

PLAIN TEXT
C++:
  1.  
  2.         // if result will lie outside [a,b]
  3.         if ( ((p-a)*speed - func)*((p-b)*speed - func) > 0.0f )
  4.         {
  5.             // do bisection
  6.             p = 0.5f*(a+b);
  7.         }   
  8.         else
  9.         {
  10.             // otherwise Newton-Raphson
  11.             p -= func/speed;
  12.         }
  13.  

And that's pretty much it. It's possible to add an additional check to see if bisection will converge faster than Newton-Raphson (see Numerical Recipes), and use bisection in that case too. However, this does add a tiny bit more overhead, and for this application the convergence of Newton-Raphson is probably fine.

FYI, I tested this out by taking three neighboring control points of a piecewise Bezier curve -- two non-interpolating and the center one interpolating -- and making them coincident. A degenerate case, but someone might want to do it, so best to be prepared.

And that's it for reparameterization. I'm not sure what's up next -- there's a new ray-box algorithm in the latest Journal of Graphics Tools, so maybe I'll talk about that.

3 Comments »

  1. what stands mTimes for? if mTimes are the t-values of a certain knot mTimes[mCount-1] – mTimes[0] must be zero so that p is equal to t1 at p = t1 + s*(mTimes[mCount-1]-mTimes[0])/mTotalLength;

    Comment by stephan — 3/23/2008 @ 9:32 am

  2. The array mTimes contains the monotonically increasing t-values at each knot, as you say. However, I don’t follow — why must (mTimes[mCount-1] – mTimes[0]) == 0? Are you questioning the initialization of p? The attempt here is to find an estimate to the time following t1 that corresponds to travelling distance s. It doesn’t have to be perfect, just a good guess — the iterations will refine it. Maybe rewriting it like this will help:

    float p = t1 + (mTimes[mCount-1]-mTimes[0])*s/mTotalLength;

    We start at t1 and add a percentage of the total time based on a ratio of the given length over the total length. p will equal t1 when s == 0 (or if mTimes[mCount-1] == mTimes[0], but then the curve is either a point or doesn’t make sense). There are issues if s is too big or too small, but that’s covered by the conditionals before we set p (lines 10-14 in the first code segment).

    Comment by Jim — 3/23/2008 @ 11:19 am

  3. ok, that was a error in reasoning. i thought mTimes[mCount-1] (the last knot of the curve) would be zero but it is 1 unlike mTimes[0] which is the first knot of the curve with 0. so

    p = t1 + s*(mTimes[mCount-1]-mTimes[0])/mTotalLength;

    is equal to

    p = t1 + s*(1-0)/mTotalLength;

    right?

    Comment by stephan — 3/23/2008 @ 6:34 pm

RSS feed for comments on this post.

Leave a comment

Powered by WordPress

Bad Behavior has blocked 51 access attempts in the last 7 days.