Essential Math Weblog Thoughts on math for games

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:

[cpp]
// FindParameterByDistance() – first pass bisection-Newton version
float
IvBezier::FindParameterByDistance( float t1, float s )
{
// set endpoints
float a = t1;
float b = mTimes[mCount-1];

// ensure that we remain within valid parameter space
if ( s >= ArcLength(t1, b) )
return b;
if ( s < = 0.0f ) return a; // make first guess float p = t1 + s*(mTimes[mCount-1]-mTimes[0])/mTotalLength; // iterate and look for zeros for ( u_int i = 0; i < 32; ++i ) { // compute function value and test against zero float func = ArcLength(t1, p) - s; if ( ::IvAbs(func) < 1.0e-03f ) { return p; } // update endpoints if ( func < 0.0f ) { a = p; } else { b = p; } // get speed along curve float speed = Velocity(p).Length(); // if something tells us we need to use bisection if ( bisection ) { // do bisection p = 0.5f*(a+b); } else { // otherwise Newton-Raphson p -= func/speed; } } // done iterating, return failure case return FLT_MAX; } [/cpp] 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: [cpp] // if speed along curve is zero if ( ::IsZero( speed ) ) { // do bisection p = 0.5f*(a+b); } else { // otherwise Newton-Raphson p -= func/speed; } [/cpp] 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: [cpp] p - func/speed < a || p - func/speed > b
[/cpp]
We can get rid of the division by multiplying through by speed:
[cpp]
p*speed – func < a*speed || p*speed - func > b*speed
[/cpp]
Or
[cpp]
(p-a)*speed < func || (p-b)*speed > func
[/cpp]

This test as it stands isn’t very floating point safe, so we can either add a little fuzziness by checking for
[cpp]
(p-a)*speed – func < kIntervalEps || (p-b)*speed - func > -kIntervalEps
[/cpp]
or the suggestion from Numerical Recipes is
[cpp]
((p-a)*speed – func)*((p-b)*speed – func) > 0.0f
[/cpp]

So our final iterative step becomes
[cpp]
// if result will lie outside [a,b]
if ( ((p-a)*speed – func)*((p-b)*speed – func) > 0.0f )
{
// do bisection
p = 0.5f*(a+b);
}
else
{
// otherwise Newton-Raphson
p -= func/speed;
}
[/cpp]

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.