[ft-devel] Fast bisecting BBox_Cubic_Check

[Alexei Podtelezhnikov]

                                       In memory of Paul Alexi

In the early days of FreeType2 it was decided that the extrema of a
cubic Bézier segment should be found as roots of a second-degree
derivative polynomial. An alternative bisection approach was perceived
as too slow to converge. Never mind that the square root of the
discriminant in the former algorithm is also iterative and slow. Here,
I present a reincarnation of the bisection algorithm that beats the
hell out of the current procedure. Indeed, the bisection algorithm
requires at most 17 iterations when working with 32-bit integers by
definition, while FT_SqrtFixed takes 24 iterations plus the overhead.

Benchmarking using src/tools/test_bbox.c showed that the proposed
algorithm is at least 2 times faster than the current one. Additional
benchmarking will be posted shortly.

References:
http://lists.nongnu.org/archive/html/freetype-devel/2000-11/msg00138.html
http://lists.nongnu.org/archive/html/freetype-devel/2001-04/msg00114.html
http://lists.nongnu.org/archive/html/freetype-devel/2001-04/msg00133.html
http://lists.nongnu.org/archive/html/freetype-devel/2001-04/msg00158.html

Code:

  static void
  BBox_Cubic_Check( FT_Pos   p1,
                    FT_Pos   p2,
                    FT_Pos   p3,
                    FT_Pos   p4,
                    FT_Pos*  min,
                    FT_Pos*  max )
  {
    FT_Pos  q1, q2, q3, q4;


    q1 = p1;
    q2 = p2;
    q3 = p3;
    q4 = p4;

    /* one of the off-points must be below the current minimum */
    /* for the conic segment to possibly cross the threshold   */
    while ( q2 > *max || q3 > *max )
    {
      /* determine which half contains the maximum and split */
      if (q1 + q2 > q3 + q4 ) /* first half */
      {
        q4 = ( q4 + 3 * q3 + 3 * q2 + q1 ) / 8;
        q3 = ( q3 + 2 * q2 + q1 ) / 4;
        q2 = ( q2 + q1 ) / 2;
      }
      else                    /* second half */
      {
        q1 = ( q1 + 3 * q2 + 3 * q3 + q4 ) / 8;
        q2 = ( q2 + 2 * q3 + q4 ) / 4;
        q3 = ( q3 + q4 ) / 2;
      }

      /* check if either end reached the maximum */
      if ( q1 == q2 )
      {
        *max = q1;
        break;
      }
      if ( q3 == q4 )
      {
        *max = q4;
        break;
      }
    }

    q1 = p1;
    q2 = p2;
    q3 = p3;
    q4 = p4;

    /* one of the off-points must be below the current minimum */
    /* for the conic segment to possibly cross the threshold   */
    while ( q2 < *min || q3 < *min )
    {
      /* determine which half contains the minimum and split */
      if (q1 + q2 < q3 + q4 ) /* first half */
      {
        q4 = ( q4 + 3 * q3 + 3 * q2 + q1 ) / 8;
        q3 = ( q3 + 2 * q2 + q1 ) / 4;
        q2 = ( q2 + q1 ) / 2;
      }
      else                    /* second half */
      {
        q1 = ( q1 + 3 * q2 + 3 * q3 + q4 ) / 8;
        q2 = ( q2 + 2 * q3 + q4 ) / 4;
        q3 = ( q3 + q4 ) / 2;
      }

      /* check if either end reached the minimum */
      if ( q1 == q2 )
      {
        *min = q1;
        break;
      }
      if ( q3 == q4 )
      {
        *min = q4;
        break;
      }
    }
  }