c++math

Numerically stable Angle bisector algorithm


Is there any numerically stable angle bisector algorithm?

The problem is the following:

Actually I'm computing it in the following way:

The problem with my approach is that when the angle is almost 180° (AB almost parallel to BC) the bisector is very inaccurate (of course because mid point is almost coincident with B). The current algorithm is so inaccurate that sometimes the resulting bisector is almost parallel to one of the other 2 segments.

And yes there are no "cast" problems, all computations are done in single precision floating point.


Solution

  • It’s not trivial. Let’s say the two edge vectors are a and b:

    float2 a = A - B;
    float2 b = C - B;
    
    1. Compute the dot product float dp = dot( a, b )

    2. Normalize both vectors:

    float2 a_norm = normalize( a );
    float2 b_norm = normalize( b );
    
    1. Check the sign bit of the dot product. When the dp is non-negative, return normalize( a_norm + b_norm ); and you’re done.

    2. When the dot product is negative, you have obtuse angle between input vectors. Applying a naïve formula in this case would screw up the numerical precision. Need another way.

    float2 c = normalize( a_norm - b_norm );
    float dir = dot( a, rotate90( b ) );
    return ( dir < 0 ) ? rotate90( c ) : rotate270( c );
    

    Note - instead of the +, this is what gives the precision win. When the angle between a and b is greater than 90°, the angle between a and -b is less than 90°, and the length of a_norm - b_norm is large enough to give accurate direction. We just need to rotate it by 90° afterwards, in the correct direction.

    P.S. Rotating 2D vectors by multiples of 90° is lossless operation. Here’s pseudocode for rotate90 and rotate270 functions:

    float2 rotate90( float2 vec )
    {
        return float2( vec.y, -vec.x );
    }
    float2 rotate270( float2 vec )
    {
        return float2( -vec.y, vec.x );
    }
    

    Update 2025: Here’s an efficient C++ implementation of that algorithm for AMD64 processors. The compiled codes look decent: less than 40 instructions, no branches or memory transactions.