[sdf] Added shortest distance finding functions.

* src/sdf/ftsdf.c (get_min_distance_): Added function to find the closest distance
  from a point to a curves of all three different types (i.e. line segment, conic
  bezier and cubic bezier).

1 files changed, 937 insertions, 0 deletions

diff --git a/src/sdf/ftsdf.c b/src/sdf/ftsdf.c
index 64048058f..e38aa7401 100644
--- a/src/sdf/ftsdf.c
+++ b/src/sdf/ftsdf.c
@@ -1569,4 +1569,941 @@
     return FT_ABS( sdf1.cross ) > FT_ABS( sdf2.cross ) ? sdf1 : sdf2;
   }
 
+  /**************************************************************************
+   *
+   * @Function:
+   *   get_min_distance_line
+   *
+   * @Description:
+   *   This function find the shortest distance from the `line' to
+   *   a given `point' and assigns it to `out'. Only use it for line
+   *   segments.
+   *
+   * @Input:
+   *   line ::
+   *     The line segment to which the shortest distance is to be
+   *     computed.
+   *
+   *   point ::
+   *     Point from which the shortest distance is to be computed.
+   *
+   * @Return:
+   *   out ::
+   *     Signed distance from the `point' to the `line'.
+   *
+   *   FT_Error ::
+   *     FreeType error, 0 means success.
+   *
+   * @Note:
+   *   The `line' parameter must have a `edge_type' of `SDF_EDGE_LINE'.
+   *
+   */
+  static FT_Error
+  get_min_distance_line( SDF_Edge*             line,
+                         FT_26D6_Vec           point, 
+                         SDF_Signed_Distance*  out )
+  {
+    /* in order to calculate the shortest distance from a point to */
+    /* a line segment.                                             */
+    /*                                                             */
+    /* a = start point of the line segment                         */
+    /* b = end point of the line segment                           */
+    /* p = point from which shortest distance is to be calculated  */
+    /* ----------------------------------------------------------- */
+    /* => we first write the parametric equation of the line       */
+    /*    point_on_line = a + ( b - a ) * t ( t is the factor )    */
+    /*                                                             */
+    /* => next we find the projection of point p on the line. the  */
+    /*    projection will be perpendicular to the line, that is    */
+    /*    why we can find it by making the dot product zero.       */
+    /*    ( point_on_line - a ) . ( p - point_on_line ) = 0        */
+    /*                                                             */
+    /*                 ( point_on_line )                           */
+    /*    ( a ) x-------o----------------x ( b )                   */
+    /*                |_|                                          */
+    /*                  |                                          */
+    /*                  |                                          */
+    /*                ( p )                                        */
+    /*                                                             */
+    /* => by simplifying the above equation we get the factor of   */
+    /*    point_on_line such that                                  */
+    /*    t = ( ( p - a ) . ( b - a ) ) / ( |b - a| ^ 2 )          */
+    /*                                                             */
+    /* => we clamp the factor t between [0.0f, 1.0f], because the  */
+    /*    point_on_line can be outside the line segment.           */
+    /*                                                             */
+    /*                                        ( point_on_line )    */
+    /*    ( a ) x------------------------x ( b ) -----o---         */
+    /*                                              |_|            */
+    /*                                                |            */
+    /*                                                |            */
+    /*                                              ( p )          */
+    /*                                                             */
+    /* => finally the distance becomes | point_on_line - p |       */
+
+    FT_Error     error = FT_Err_Ok;
+
+    FT_Vector    a;                 /* start position */
+    FT_Vector    b;                 /* end position   */
+    FT_Vector    p;                 /* current point  */
+
+    FT_26D6_Vec  line_segment;      /* `b' - `a'*/
+    FT_26D6_Vec  p_sub_a;           /* `p' - `a' */
+
+    FT_26D6      sq_line_length;    /* squared length of `line_segment' */
+    FT_16D16     factor;            /* factor of the nearest point      */
+    FT_26D6      cross;             /* used to determine sign           */
+
+    FT_16D16_Vec nearest_point;     /* `point_on_line'       */
+    FT_16D16_Vec nearest_vector;    /* `p' - `nearest_point' */
+
+
+    if ( !line || !out )
+    {
+      error = FT_THROW( Invalid_Argument );
+      goto Exit;
+    }
+
+    if ( line->edge_type != SDF_EDGE_LINE )
+    {
+      error = FT_THROW( Invalid_Argument );
+      goto Exit;
+    }
+
+    a = line->start_pos;
+    b = line->end_pos;
+    p = point;
+
+    line_segment.x = b.x - a.x;
+    line_segment.y = b.y - a.y;
+
+    p_sub_a.x = p.x - a.x;
+    p_sub_a.y = p.y - a.y;
+
+    sq_line_length = ( line_segment.x * line_segment.x ) / 64 +
+                     ( line_segment.y * line_segment.y ) / 64;
+
+    /* currently factor is 26.6 */
+    factor = ( p_sub_a.x * line_segment.x ) / 64 +
+             ( p_sub_a.y * line_segment.y ) / 64;
+
+    /* now factor is 16.16 */
+    factor = FT_DivFix( factor, sq_line_length );
+
+    /* clamp the factor between 0.0 and 1.0 in fixed point */
+    if ( factor > FT_INT_16D16( 1 ) )
+      factor = FT_INT_16D16( 1 );
+    if ( factor < 0 )
+      factor = 0;
+
+    nearest_point.x = FT_MulFix( FT_26D6_16D16(line_segment.x),
+                                 factor );
+    nearest_point.y = FT_MulFix( FT_26D6_16D16(line_segment.y),
+                                 factor );
+
+    nearest_point.x = FT_26D6_16D16( a.x ) + nearest_point.x; 
+    nearest_point.y = FT_26D6_16D16( a.y ) + nearest_point.y;
+
+    nearest_vector.x = nearest_point.x - FT_26D6_16D16( p.x );
+    nearest_vector.y = nearest_point.y - FT_26D6_16D16( p.y );
+
+    cross = FT_MulFix( nearest_vector.x, line_segment.y ) -
+            FT_MulFix( nearest_vector.y, line_segment.x );
+
+    /* assign the output */
+    out->sign = cross < 0 ? 1 : -1;
+    out->distance = VECTOR_LENGTH_16D16( nearest_vector );
+
+    /* Instead of finding cross for checking corner we */
+    /* directly set it here. This is more efficient    */
+    /* because if the distance is perpendicular we can */
+    /* directly set it to 1.                           */
+    if ( factor != 0 && factor != FT_INT_16D16( 1 ) )
+      out->cross = FT_INT_16D16( 1 );
+    else
+    {
+      /* [OPTIMIZATION]: Pre-compute this direction. */
+      /* if not perpendicular then compute the cross */
+      FT_Vector_NormLen( &line_segment );
+      FT_Vector_NormLen( &nearest_vector );
+
+      out->cross = FT_MulFix( line_segment.x, nearest_vector.y ) -
+                   FT_MulFix( line_segment.y, nearest_vector.x );
+    }
+
+  Exit:
+    return error;
+  }
+
+#if !USE_NEWTON_FOR_CONIC
+
+  /**************************************************************************
+   *
+   * @Function:
+   *   get_min_distance_conic
+   *
+   * @Description:
+   *   This function find the shortest distance from the `conic' bezier
+   *   curve to a given `point' and assigns it to `out'. Only use it for
+   *   conic/quadratic curves.
+   *
+   * @Input:
+   *   conic ::
+   *     The conic bezier to which the shortest distance is to be
+   *     computed.
+   *
+   *   point ::
+   *     Point from which the shortest distance is to be computed.
+   *
+   * @Return:
+   *   out ::
+   *     Signed distance from the `point' to the `conic'.
+   *
+   *   FT_Error ::
+   *     FreeType error, 0 means success.
+   *
+   * @Note:
+   *   The function uses analytical method to find shortest distance
+   *   which is faster than the Newton-Raphson's method, but has
+   *   underflows at the moment. Use Newton's method if you can
+   *   see artifacts in the SDF.
+   *
+   *   The `conic' parameter must have a `edge_type' of `SDF_EDGE_CONIC'.
+   *
+   */
+  static FT_Error
+  get_min_distance_conic( SDF_Edge*             conic,
+                          FT_26D6_Vec           point,
+                          SDF_Signed_Distance*  out )
+  {
+      /* The procedure to find the shortest distance from a point to */
+      /* a quadratic bezier curve is similar to a line segment. the  */
+      /* shortest distance will be perpendicular to the bezier curve */
+      /* The only difference from line is that there can be more     */
+      /* than one perpendicular and we also have to check the endpo- */
+      /* -ints, because the perpendicular may not be the shortest.   */
+      /*                                                             */
+      /* p0 = first endpoint                                         */
+      /* p1 = control point                                          */
+      /* p2 = second endpoint                                        */
+      /* p = point from which shortest distance is to be calculated  */
+      /* ----------------------------------------------------------- */
+      /* => the equation of a quadratic bezier curve can be written  */
+      /*    B( t ) = ( ( 1 - t )^2 )p0 + 2( 1 - t )tp1 + t^2p2       */
+      /*    here t is the factor with range [0.0f, 1.0f]             */
+      /*    the above equation can be rewritten as                   */
+      /*    B( t ) = t^2( p0 - 2p1 + p2 ) + 2t( p1 - p0 ) + p0       */
+      /*                                                             */
+      /*    now let A = ( p0 - 2p1 + p2), B = ( p1 - p0 )            */
+      /*    B( t ) = t^2( A ) + 2t( B ) + p0                         */
+      /*                                                             */
+      /* => the derivative of the above equation is written as       */
+      /*    B'( t ) = 2( tA + B )                                    */
+      /*                                                             */
+      /* => now to find the shortest distance from p to B( t ), we   */
+      /*    find the point on the curve at which the shortest        */
+      /*    distance vector ( i.e. B( t ) - p ) and the direction    */
+      /*    ( i.e. B'( t )) makes 90 degrees. in other words we make */
+      /*    the dot product zero.                                    */
+      /*    ( B( t ) - p ).( B'( t ) ) = 0                           */
+      /*    ( t^2( A ) + 2t( B ) + p0 - p ).( 2( tA + B ) ) = 0      */
+      /*                                                             */
+      /*    after simplifying we get a cubic equation as             */
+      /*    at^3 + bt^2 + ct + d = 0                                 */
+      /*    a = ( A.A ), b = ( 3A.B ), c = ( 2B.B + A.p0 - A.p )     */
+      /*    d = ( p0.B - p.B )                                       */
+      /*                                                             */
+      /* => now the roots of the equation can be computed using the  */
+      /*    `Cardano's Cubic formula' we clamp the roots in range    */
+      /*    [0.0f, 1.0f].                                            */
+      /*                                                             */
+      /* [note]: B and B( t ) are different in the above equations   */
+
+    FT_Error     error = FT_Err_Ok;
+
+    FT_26D6_Vec  aA, bB;         /* A, B in the above comment             */
+    FT_26D6_Vec  nearest_point;  /* point on curve nearest to `point'     */
+    FT_26D6_Vec  direction;      /* direction of curve at `nearest_point' */
+
+    FT_26D6_Vec  p0, p1, p2;     /* control points of a conic curve       */
+    FT_26D6_Vec  p;              /* `point' to which shortest distance    */
+
+    FT_26D6      a, b, c, d;     /* cubic coefficients                    */
+
+    FT_16D16     roots[3] = { 0, 0, 0 }; /* real roots of the cubic eq    */
+    FT_16D16     min_factor;             /* factor at `nearest_point'     */
+    FT_16D16     cross;                  /* to determine the sign         */
+    FT_16D16     min = FT_INT_MAX;       /* shortest squared distance     */
+
+    FT_UShort    num_roots;              /* number of real roots of cubic */
+    FT_UShort    i;
+
+
+    if ( !conic || !out )
+    {
+      error = FT_THROW( Invalid_Argument );
+      goto Exit;
+    }
+
+    if ( conic->edge_type != SDF_EDGE_CONIC )
+    {
+      error = FT_THROW( Invalid_Argument );
+      goto Exit;
+    }
+
+    /* assign the values after checking pointer */
+    p0 = conic->start_pos;
+    p1 = conic->control_a;
+    p2 = conic->end_pos;
+    p  = point;
+
+    /* compute substitution coefficients */
+    aA.x = p0.x - 2 * p1.x + p2.x;
+    aA.y = p0.y - 2 * p1.y + p2.y;
+
+    bB.x = p1.x - p0.x;
+    bB.y = p1.y - p0.y;
+
+    /* compute cubic coefficients */
+    a = VEC_26D6_DOT( aA, aA );
+
+    b = 3 * VEC_26D6_DOT( aA, bB );
+
+    c = 2 * VEC_26D6_DOT( bB, bB ) +
+            VEC_26D6_DOT( aA, p0 ) -
+            VEC_26D6_DOT( aA, p );
+
+    d = VEC_26D6_DOT( p0, bB ) -
+        VEC_26D6_DOT( p, bB );
+
+    /* find the roots */
+    num_roots = solve_cubic_equation( a, b, c, d, roots );
+
+    if ( num_roots == 0 )
+    {
+      roots[0] = 0;
+      roots[1] = FT_INT_16D16( 1 );
+      num_roots = 2;
+    }
+
+    /* [OPTIMIZATION]: Check the roots, clamp them and discard */
+    /*                 duplicate roots.                        */
+
+    /* convert these values to 16.16 for further computation */
+    aA.x = FT_26D6_16D16( aA.x );
+    aA.y = FT_26D6_16D16( aA.y );
+
+    bB.x = FT_26D6_16D16( bB.x );
+    bB.y = FT_26D6_16D16( bB.y );
+
+    p0.x = FT_26D6_16D16( p0.x );
+    p0.y = FT_26D6_16D16( p0.y );
+
+    p.x = FT_26D6_16D16( p.x );
+    p.y = FT_26D6_16D16( p.y );
+
+    for ( i = 0; i < num_roots; i++ )
+    {
+      FT_16D16      t    = roots[i];
+      FT_16D16      t2   = 0;
+      FT_16D16      dist = 0;
+
+      FT_16D16_Vec  curve_point;
+      FT_16D16_Vec  dist_vector;
+
+      /* Ideally we should discard the roots which are outside the     */
+      /* range [0.0, 1.0] and check the endpoints of the bezier, but   */
+      /* Behdad gave me a lemma:                                       */
+      /* Lemma:                                                        */
+      /* * If the closest point on the curve [0, 1] is to the endpoint */
+      /*   at t = 1 and the cubic has no real roots at t = 1 then, the */
+      /*   cubic must have a real root at some t > 1.                  */
+      /* * Similarly,                                                  */
+      /*   If the closest point on the curve [0, 1] is to the endpoint */
+      /*   at t = 0 and the cubic has no real roots at t = 0 then, the */
+      /*   cubic must have a real root at some t < 0.                  */
+      /*                                                               */
+      /* Now because of this lemma we only need to clamp the roots and */
+      /* that will take care of the endpoints.                         */
+      /*                                                               */
+      /* For proof contact: behdad@behdad.org                          */
+      /* For more details check message:                               */
+      /* https://lists.nongnu.org/archive/html/freetype-devel/2020-06/msg00147.html */
+      if ( t < 0 )
+        t = 0;
+      if ( t > FT_INT_16D16( 1 ) )
+        t = FT_INT_16D16( 1 );
+
+      t2 = FT_MulFix( t, t );
+
+      /* B( t ) = t^2( A ) + 2t( B ) + p0 - p */
+      curve_point.x = FT_MulFix( aA.x, t2 ) +
+                      2 * FT_MulFix( bB.x, t ) + p0.x;
+      curve_point.y = FT_MulFix( aA.y, t2 ) +
+                      2 * FT_MulFix( bB.y, t ) + p0.y;
+
+      /* `curve_point' - `p' */
+      dist_vector.x = curve_point.x - p.x;
+      dist_vector.y = curve_point.y - p.y;
+
+      dist = VECTOR_LENGTH_16D16( dist_vector );
+
+      if ( dist < min )
+      {
+        min = dist;
+        nearest_point = curve_point;
+        min_factor = t;
+      }
+    }
+
+    /* B'( t ) = 2( tA + B ) */
+    direction.x = 2 * FT_MulFix( aA.x, min_factor ) + 2 * bB.x;
+    direction.y = 2 * FT_MulFix( aA.y, min_factor ) + 2 * bB.y;
+
+    /* determine the sign */
+    cross = FT_MulFix( nearest_point.x - p.x, direction.y ) -
+            FT_MulFix( nearest_point.y - p.y, direction.x );
+
+    /* assign the values */
+    out->distance = min;
+    out->sign = cross < 0 ? 1 : -1;
+
+    if ( min_factor != 0 && min_factor != FT_INT_16D16( 1 ) )
+      out->cross = FT_INT_16D16( 1 );   /* the two are perpendicular */
+    else
+    {
+      /* convert to nearest vector */
+      nearest_point.x -= FT_26D6_16D16( p.x );
+      nearest_point.y -= FT_26D6_16D16( p.y );
+
+      /* if not perpendicular then compute the cross */
+      FT_Vector_NormLen( &direction );
+      FT_Vector_NormLen( &nearest_point );
+
+      out->cross = FT_MulFix( direction.x, nearest_point.y ) -
+                   FT_MulFix( direction.y, nearest_point.x );
+    }
+  Exit:
+    return error;
+  }
+
+#else
+
+  /**************************************************************************
+   *
+   * @Function:
+   *   get_min_distance_conic
+   *
+   * @Description:
+   *   This function find the shortest distance from the `conic' bezier
+   *   curve to a given `point' and assigns it to `out'. Only use it for
+   *   conic/quadratic curves.
+   *
+   * @Input:
+   *   conic ::
+   *     The conic bezier to which the shortest distance is to be
+   *     computed.
+   *
+   *   point ::
+   *     Point from which the shortest distance is to be computed.
+   *
+   * @Return:
+   *   out ::
+   *     Signed distance from the `point' to the `conic'.
+   *
+   *   FT_Error ::
+   *     FreeType error, 0 means success.
+   *
+   * @Note:
+   *   The function uses Newton's approximation to find the shortest
+   *   distance, which is a bit slower than the analytical method but
+   *   doesn't cause underflow. Use is upto your needs.
+   *
+   *   The `conic' parameter must have a `edge_type' of `SDF_EDGE_CONIC'.
+   *
+   */
+  static FT_Error
+  get_min_distance_conic( SDF_Edge*             conic,
+                          FT_26D6_Vec           point,
+                          SDF_Signed_Distance*  out )
+  {
+    /* This method uses Newton-Raphson's approximation to find the */
+    /* shortest distance from a point to conic curve which does    */
+    /* not involve solving any cubic equation, that is why there   */
+    /* is no risk of underflow. The method is as follows:          */
+    /*                                                             */
+    /* p0 = first endpoint                                         */
+    /* p1 = control point                                          */
+    /* p3 = second endpoint                                        */
+    /* p  = point from which shortest distance is to be calculated */
+    /* ----------------------------------------------------------- */
+    /* => the equation of a quadratic bezier curve can be written  */
+    /*    B( t ) = ( ( 1 - t )^2 )p0 + 2( 1 - t )tp1 + t^2p2       */
+    /*    here t is the factor with range [0.0f, 1.0f]             */
+    /*    the above equation can be rewritten as                   */
+    /*    B( t ) = t^2( p0 - 2p1 + p2 ) + 2t( p1 - p0 ) + p0       */
+    /*                                                             */
+    /*    now let A = ( p0 - 2p1 + p2), B = 2( p1 - p0 )           */
+    /*    B( t ) = t^2( A ) + t( B ) + p0                          */
+    /*                                                             */
+    /* => the derivative of the above equation is written as       */
+    /*    B'( t ) = 2t( A ) + B                                    */
+    /*                                                             */
+    /* => further derivative of the above equation is written as   */
+    /*    B''( t ) = 2A                                            */
+    /*                                                             */
+    /* => the equation of distance from point `p' to the curve     */
+    /*    P( t ) can be written as                                 */
+    /*    P( t ) = t^2( A ) + t^2( B ) + p0 - p                    */
+    /*    Now let C = ( p0 - p )                                   */
+    /*    P( t ) = t^2( A ) + t( B ) + C                           */
+    /*                                                             */
+    /* => finally the equation of angle between curve B( t ) and   */
+    /*    point to curve distance P( t ) can be written as         */
+    /*    Q( t ) = P( t ).B'( t )                                  */
+    /*                                                             */
+    /* => now our task is to find a value of t such that the above */
+    /*    equation Q( t ) becomes zero. in other words the point   */
+    /*    to curve vector makes 90 degree with curve. this is done */
+    /*    by Newton-Raphson's method.                              */
+    /*                                                             */
+    /* => we first assume a arbitary value of the factor `t' and   */
+    /*    then we improve it using Newton's equation such as       */
+    /*                                                             */
+    /*    t -= Q( t ) / Q'( t )                                    */
+    /*    putting value of Q( t ) from the above equation gives    */
+    /*                                                             */
+    /*    t -= P( t ).B'( t ) / derivative( P( t ).B'( t ) )       */
+    /*    t -= P( t ).B'( t ) /                                    */
+    /*         ( P'( t )B'( t ) + P( t ).B''( t ) )                */
+    /*                                                             */
+    /*    P'( t ) is noting but B'( t ) because the constant are   */
+    /*    gone due to derivative                                   */
+    /*                                                             */
+    /* => finally we get the equation to improve the factor as     */
+    /*    t -= P( t ).B'( t ) /                                    */
+    /*         ( B'( t ).B'( t ) + P( t ).B''( t ) )               */
+    /*                                                             */
+    /* [note]: B and B( t ) are different in the above equations   */
+
+    FT_Error     error = FT_Err_Ok;
+
+    FT_26D6_Vec  aA, bB, cC;     /* A, B, C in the above comment          */
+    FT_26D6_Vec  nearest_point;  /* point on curve nearest to `point'     */
+    FT_26D6_Vec  direction;      /* direction of curve at `nearest_point' */
+
+    FT_26D6_Vec  p0, p1, p2;     /* control points of a conic curve       */
+    FT_26D6_Vec  p;              /* `point' to which shortest distance    */
+
+    FT_16D16     min_factor = 0;         /* factor at `nearest_point'     */
+    FT_16D16     cross;                  /* to determine the sign         */
+    FT_16D16     min = FT_INT_MAX;       /* shortest squared distance     */
+
+    FT_UShort    iterations;
+    FT_UShort    steps;
+
+
+    if ( !conic || !out )
+    {
+      error = FT_THROW( Invalid_Argument );
+      goto Exit;
+    }
+
+    if ( conic->edge_type != SDF_EDGE_CONIC )
+    {
+      error = FT_THROW( Invalid_Argument );
+      goto Exit;
+    }
+
+    /* assign the values after checking pointer */
+    p0 = conic->start_pos;
+    p1 = conic->control_a;
+    p2 = conic->end_pos;
+    p  = point;
+
+    /* compute substitution coefficients */
+    aA.x = p0.x - 2 * p1.x + p2.x;
+    aA.y = p0.y - 2 * p1.y + p2.y;
+
+    bB.x = 2 * ( p1.x - p0.x );
+    bB.y = 2 * ( p1.y - p0.y );
+
+    cC.x = p0.x;
+    cC.y = p0.y;
+
+    /* do newton's iterations */
+    for ( iterations = 0; iterations <= MAX_NEWTON_DIVISIONS; iterations++ )
+    {
+      FT_16D16  factor = FT_INT_16D16( iterations ) / MAX_NEWTON_DIVISIONS;
+      FT_16D16  factor2;
+      FT_16D16  length;
+
+      FT_16D16_Vec  curve_point; /* point on the curve  */
+      FT_16D16_Vec  dist_vector; /* `curve_point' - `p' */
+
+      FT_26D6_Vec   d1;          /* first  derivative   */
+      FT_26D6_Vec   d2;          /* second derivative   */
+
+      FT_16D16      temp1;
+      FT_16D16      temp2;
+
+      for ( steps = 0; steps < MAX_NEWTON_STEPS; steps++ )
+      {
+        factor2 = FT_MulFix( factor, factor );
+
+        /* B( t ) = t^2( A ) + t( B ) + p0 */
+        curve_point.x = FT_MulFix( aA.x, factor2 ) +
+                        FT_MulFix( bB.x, factor ) + cC.x;
+        curve_point.y = FT_MulFix( aA.y, factor2 ) +
+                        FT_MulFix( bB.y, factor ) + cC.y;
+
+        /* convert to 16.16 */
+        curve_point.x = FT_26D6_16D16( curve_point.x );
+        curve_point.y = FT_26D6_16D16( curve_point.y );
+
+        /* B( t ) = t^2( A ) + t( B ) + p0 - p. P( t ) in the comment */
+        dist_vector.x = curve_point.x - FT_26D6_16D16( p.x );
+        dist_vector.y = curve_point.y - FT_26D6_16D16( p.y );
+
+        length = VECTOR_LENGTH_16D16( dist_vector );
+
+        if ( length < min )
+        {
+          min = length;
+          min_factor = factor;
+          nearest_point = curve_point;
+        }
+
+        /* This the actual Newton's approximation.       */
+        /*    t -= P( t ).B'( t ) /                      */
+        /*         ( B'( t ).B'( t ) + P( t ).B''( t ) ) */
+
+        /* B'( t ) = 2tA + B */
+        d1.x = FT_MulFix( aA.x, 2 * factor ) + bB.x;
+        d1.y = FT_MulFix( aA.y, 2 * factor ) + bB.y;
+
+        /* B''( t ) = 2A */
+        d2.x = 2 * aA.x;
+        d2.y = 2 * aA.y;
+
+        dist_vector.x /= 1024;
+        dist_vector.y /= 1024;
+
+        /* temp1 = P( t ).B'( t ) */
+        temp1 = VEC_26D6_DOT( dist_vector, d1 );
+
+        /* temp2 = ( B'( t ).B'( t ) + P( t ).B''( t ) ) */
+        temp2 = VEC_26D6_DOT( d1, d1 ) +
+                VEC_26D6_DOT( dist_vector, d2 );
+
+        factor -= FT_DivFix( temp1, temp2 );
+
+        if ( factor < 0 || factor > FT_INT_16D16( 1 ) )
+          break;
+      }
+    }
+
+    /* B'( t ) = 2tA + B */
+    direction.x = 2 * FT_MulFix( aA.x, min_factor ) + bB.x;
+    direction.y = 2 * FT_MulFix( aA.y, min_factor ) + bB.y;
+
+    /* determine the sign */
+    cross = FT_MulFix( nearest_point.x - FT_26D6_16D16( p.x ), direction.y ) -
+            FT_MulFix( nearest_point.y - FT_26D6_16D16( p.y ), direction.x );
+
+    /* assign the values */
+    out->distance = min;
+    out->sign = cross < 0 ? 1 : -1;
+
+    if ( min_factor != 0 && min_factor != FT_INT_16D16( 1 ) )
+      out->cross = FT_INT_16D16( 1 );   /* the two are perpendicular */
+    else
+    {
+      /* convert to nearest vector */
+      nearest_point.x -= FT_26D6_16D16( p.x );
+      nearest_point.y -= FT_26D6_16D16( p.y );
+
+      /* if not perpendicular then compute the cross */
+      FT_Vector_NormLen( &direction );
+      FT_Vector_NormLen( &nearest_point );
+
+      out->cross = FT_MulFix( direction.x, nearest_point.y ) -
+                   FT_MulFix( direction.y, nearest_point.x );
+    }
+
+  Exit:
+    return error;
+  }
+
+#endif
+
+  /**************************************************************************
+   *
+   * @Function:
+   *   get_min_distance_cubic
+   *
+   * @Description:
+   *   This function find the shortest distance from the `cubic' bezier
+   *   curve to a given `point' and assigns it to `out'. Only use it for
+   *   cubic curves.
+   *
+   * @Input:
+   *   cubic ::
+   *     The cubic bezier to which the shortest distance is to be
+   *     computed.
+   *
+   *   point ::
+   *     Point from which the shortest distance is to be computed.
+   *
+   * @Return:
+   *   out ::
+   *     Signed distance from the `point' to the `cubic'.
+   *
+   *   FT_Error ::
+   *     FreeType error, 0 means success.
+   *
+   * @Note:
+   *   The function uses Newton's approximation to find the shortest
+   *   distance. Another way would be to divide the cubic into conic
+   *   or subdivide the curve into lines, but that is not implemented.
+   *
+   *   The `cubic' parameter must have a `edge_type' of `SDF_EDGE_CUBIC'.
+   *
+   */
+  static FT_Error
+  get_min_distance_cubic( SDF_Edge*             cubic,
+                          FT_26D6_Vec           point,
+                          SDF_Signed_Distance*  out )
+  {
+    /* the procedure to find the shortest distance from a point to */
+    /* a cubic bezier curve is similar to a quadratic curve.       */
+    /* The only difference is that while calculating the factor    */
+    /* `t', instead of a cubic polynomial equation we have to find */
+    /* the roots of a 5th degree polynomial equation.              */
+    /* But since solving a 5th degree polynomial equation require  */
+    /* significant amount of time and still the results may not be */
+    /* accurate, we are going to directly approximate the value of */
+    /* `t' using Newton-Raphson method                             */
+    /*                                                             */
+    /* p0 = first endpoint                                         */
+    /* p1 = first control point                                    */
+    /* p2 = second control point                                   */
+    /* p3 = second endpoint                                        */
+    /* p = point from which shortest distance is to be calculated  */
+    /* ----------------------------------------------------------- */
+    /* => the equation of a cubic bezier curve can be written as:  */
+    /*    B( t ) = ( ( 1 - t )^3 )p0 + 3( ( 1 - t )^2 )tp1 +       */
+    /*             3( 1 - t )( t^2 )p2 + ( t^3 )p3                 */
+    /*    The equation can be expanded and written as:             */
+    /*    B( t ) = ( t^3 )( -p0 + 3p1 - 3p2 + p3 ) +               */
+    /*             3( t^2 )( p0 - 2p1 + p2 ) + 3t( -p0 + p1 ) + p0 */
+    /*                                                             */
+    /*    Now let A = ( -p0 + 3p1 - 3p2 + p3 ),                    */
+    /*            B = 3( p0 - 2p1 + p2 ), C = 3( -p0 + p1 )        */
+    /*    B( t ) = t^3( A ) + t^2( B ) + tC + p0                   */
+    /*                                                             */
+    /* => the derivative of the above equation is written as       */
+    /*    B'( t ) = 3t^2( A ) + 2t( B ) + C                        */
+    /*                                                             */
+    /* => further derivative of the above equation is written as   */
+    /*    B''( t ) = 6t( A ) + 2B                                  */
+    /*                                                             */
+    /* => the equation of distance from point `p' to the curve     */
+    /*    P( t ) can be written as                                 */
+    /*    P( t ) = t^3( A ) + t^2( B ) + tC + p0 - p               */
+    /*    Now let D = ( p0 - p )                                   */
+    /*    P( t ) = t^3( A ) + t^2( B ) + tC + D                    */
+    /*                                                             */
+    /* => finally the equation of angle between curve B( t ) and   */
+    /*    point to curve distance P( t ) can be written as         */
+    /*    Q( t ) = P( t ).B'( t )                                  */
+    /*                                                             */
+    /* => now our task is to find a value of t such that the above */
+    /*    equation Q( t ) becomes zero. in other words the point   */
+    /*    to curve vector makes 90 degree with curve. this is done */
+    /*    by Newton-Raphson's method.                              */
+    /*                                                             */
+    /* => we first assume a arbitary value of the factor `t' and   */
+    /*    then we improve it using Newton's equation such as       */
+    /*                                                             */
+    /*    t -= Q( t ) / Q'( t )                                    */
+    /*    putting value of Q( t ) from the above equation gives    */
+    /*                                                             */
+    /*    t -= P( t ).B'( t ) / derivative( P( t ).B'( t ) )       */
+    /*    t -= P( t ).B'( t ) /                                    */
+    /*         ( P'( t )B'( t ) + P( t ).B''( t ) )                */
+    /*                                                             */
+    /*    P'( t ) is noting but B'( t ) because the constant are   */
+    /*    gone due to derivative                                   */
+    /*                                                             */
+    /* => finally we get the equation to improve the factor as     */
+    /*    t -= P( t ).B'( t ) /                                    */
+    /*         ( B'( t ).B'( t ) + P( t ).B''( t ) )               */
+    /*                                                             */
+    /* [note]: B and B( t ) are different in the above equations   */
+
+    FT_Error      error = FT_Err_Ok;
+
+    FT_26D6_Vec   aA, bB, cC, dD; /* A, B, C in the above comment          */
+    FT_16D16_Vec  nearest_point;  /* point on curve nearest to `point'     */
+    FT_16D16_Vec  direction;      /* direction of curve at `nearest_point' */
+
+    FT_26D6_Vec   p0, p1, p2, p3; /* control points of a cubic curve       */
+    FT_26D6_Vec   p;              /* `point' to which shortest distance    */
+
+    FT_16D16      min = FT_INT_MAX;  /* shortest distance           */
+    FT_16D16      min_factor = 0;    /* factor at shortest distance */
+    FT_16D16      min_factor_sq = 0; /* factor at shortest distance */
+    FT_16D16      cross;             /* to determine the sign       */
+
+    FT_UShort     iterations;
+    FT_UShort     steps;
+
+
+    if ( !cubic || !out )
+    {
+      error = FT_THROW( Invalid_Argument );
+      goto Exit;
+    }
+
+    if ( cubic->edge_type != SDF_EDGE_CUBIC )
+    {
+      error = FT_THROW( Invalid_Argument );
+      goto Exit;
+    }
+
+    /* assign the values after checking pointer */
+    p0 = cubic->start_pos;
+    p1 = cubic->control_a;
+    p2 = cubic->control_b;
+    p3 = cubic->end_pos;
+    p  = point;
+
+    /* compute substitution coefficients */
+    aA.x = -p0.x + 3 * ( p1.x - p2.x ) + p3.x;
+    aA.y = -p0.y + 3 * ( p1.y - p2.y ) + p3.y;
+
+    bB.x = 3 * ( p0.x - 2 * p1.x + p2.x );
+    bB.y = 3 * ( p0.y - 2 * p1.y + p2.y );
+
+    cC.x = 3 * ( p1.x - p0.x );
+    cC.y = 3 * ( p1.y - p0.y );
+
+    dD.x = p0.x;
+    dD.y = p0.y;
+
+    for ( iterations = 0; iterations <= MAX_NEWTON_DIVISIONS; iterations++ )
+    {
+      FT_16D16  factor  = FT_INT_16D16( iterations ) / MAX_NEWTON_DIVISIONS;
+
+      FT_16D16  factor2;         /* factor^2            */
+      FT_16D16  factor3;         /* factor^3            */
+      FT_16D16  length;
+
+      FT_16D16_Vec  curve_point; /* point on the curve  */
+      FT_16D16_Vec  dist_vector; /* `curve_point' - `p' */
+
+      FT_26D6_Vec   d1;          /* first  derivative   */
+      FT_26D6_Vec   d2;          /* second derivative   */
+
+      FT_16D16      temp1;
+      FT_16D16      temp2;
+
+
+      for ( steps = 0; steps < MAX_NEWTON_STEPS; steps++ )
+      {
+        factor2 = FT_MulFix( factor, factor );
+        factor3 = FT_MulFix( factor2, factor );
+
+        /* B( t ) = t^3( A ) + t^2( B ) + tC + D */
+        curve_point.x = FT_MulFix( aA.x, factor3 ) +
+                        FT_MulFix( bB.x, factor2 ) +
+                        FT_MulFix( cC.x, factor ) + dD.x;
+        curve_point.y = FT_MulFix( aA.y, factor3 ) +
+                        FT_MulFix( bB.y, factor2 ) +
+                        FT_MulFix( cC.y, factor ) + dD.y;
+
+        /* convert to 16.16 */
+        curve_point.x = FT_26D6_16D16( curve_point.x );
+        curve_point.y = FT_26D6_16D16( curve_point.y );
+
+        /* P( t ) in the comment */
+        dist_vector.x = curve_point.x - FT_26D6_16D16( p.x );
+        dist_vector.y = curve_point.y - FT_26D6_16D16( p.y );
+
+        length = VECTOR_LENGTH_16D16( dist_vector );
+
+        if ( length < min )
+        {
+          min = length;
+          min_factor = factor;
+          min_factor_sq = factor2;
+          nearest_point = curve_point;
+        }
+
+        /* This the actual Newton's approximation.       */
+        /*    t -= P( t ).B'( t ) /                      */
+        /*         ( B'( t ).B'( t ) + P( t ).B''( t ) ) */
+
+        /* B'( t ) = 3t^2( A ) + 2t( B ) + C */
+        d1.x = FT_MulFix( aA.x, 3 * factor2 ) + 
+               FT_MulFix( bB.x, 2 * factor ) + cC.x;
+        d1.y = FT_MulFix( aA.y, 3 * factor2 ) + 
+               FT_MulFix( bB.y, 2 * factor ) + cC.y;
+
+        /* B''( t ) = 6t( A ) + 2B */
+        d2.x = FT_MulFix( aA.x, 6 * factor ) + 2 * bB.x;
+        d2.y = FT_MulFix( aA.y, 6 * factor ) + 2 * bB.y;
+
+        dist_vector.x /= 1024;
+        dist_vector.y /= 1024;
+
+        /* temp1 = P( t ).B'( t ) */
+        temp1 = VEC_26D6_DOT( dist_vector, d1 );
+
+        /* temp2 = ( B'( t ).B'( t ) + P( t ).B''( t ) ) */
+        temp2 = VEC_26D6_DOT( d1, d1 ) +
+                VEC_26D6_DOT( dist_vector, d2 );
+
+        factor -= FT_DivFix( temp1, temp2 );
+
+        if ( factor < 0 || factor > FT_INT_16D16( 1 ) )
+          break;
+      }
+    }
+
+    /* B'( t ) = 3t^2( A ) + 2t( B ) + C */
+    direction.x = FT_MulFix( aA.x, 3 * min_factor_sq ) + 
+                  FT_MulFix( bB.x, 2 * min_factor ) + cC.x;
+    direction.y = FT_MulFix( aA.y, 3 * min_factor_sq ) + 
+                  FT_MulFix( bB.y, 2 * min_factor ) + cC.y;
+
+    /* determine the sign */
+    cross = FT_MulFix( nearest_point.x - FT_26D6_16D16( p.x ), direction.y ) -
+            FT_MulFix( nearest_point.y - FT_26D6_16D16( p.y ), direction.x );
+
+    /* assign the values */
+    out->distance = min;
+    out->sign = cross < 0 ? 1 : -1;
+
+    if ( min_factor != 0 && min_factor != FT_INT_16D16( 1 ) )
+      out->cross = FT_INT_16D16( 1 );   /* the two are perpendicular */
+    else
+    {
+      /* convert to nearest vector */
+      nearest_point.x -= FT_26D6_16D16( p.x );
+      nearest_point.y -= FT_26D6_16D16( p.y );
+
+      /* if not perpendicular then compute the cross */
+      FT_Vector_NormLen( &direction );
+      FT_Vector_NormLen( &nearest_point );
+
+      out->cross = FT_MulFix( direction.x, nearest_point.y ) -
+                   FT_MulFix( direction.y, nearest_point.x );
+    }
+  Exit:
+    return error;
+  }
+
+
 /* END */