[sdf] Add shortest distance finding functions.

* src/sdf/ftsdf.c (get_min_distance_line, get_min_distance_conic,
get_min_distance_cubic): New functions.  Note that
`get_min_distance_conic` comes with two implementations (using an
analytical and an iterative method, to be controlled with the
`USE_NEWTON_FOR_CONIC` macro).

2 files changed, 1091 insertions, 1 deletions

diff --git a/ChangeLog b/ChangeLog
index dc0d45d1e..5d2a01b9f 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,15 @@
 2020-08-18  Anuj Verma  <anujv@iitbhilai.ac.in>
 
+	[sdf] Add shortest distance finding functions.
+
+	* src/sdf/ftsdf.c (get_min_distance_line, get_min_distance_conic,
+	get_min_distance_cubic): New functions.  Note that
+	`get_min_distance_conic` comes with two implementations (using an
+	analytical and an iterative method, to be controlled with the
+	`USE_NEWTON_FOR_CONIC` macro).
+
+2020-08-18  Anuj Verma  <anujv@iitbhilai.ac.in>
+
 	[sdf] Add function to resolve corner distances.
 
 	* src/sdf/ftsdf.c (resolve_corner): New function.
diff --git a/src/sdf/ftsdf.c b/src/sdf/ftsdf.c
index 251215530..9f7f57d82 100644
--- a/src/sdf/ftsdf.c
+++ b/src/sdf/ftsdf.c
@@ -1631,7 +1631,7 @@
    *
    *   The orthogonality is simply the sinus of the two vectors (i.e.,
    *   x - o) and the corresponding direction.  We efficiently pre-compute
-   *   the orthogonality with the corresponding `get_min_distance_`
+   *   the orthogonality with the corresponding `get_min_distance_*`
    *   functions.
    *
    * @Input:
@@ -1660,4 +1660,1084 @@
     return FT_ABS( sdf1.cross ) > FT_ABS( sdf2.cross ) ? sdf1 : sdf2;
   }
 
+
+  /**************************************************************************
+   *
+   * @Function:
+   *   get_min_distance_line
+   *
+   * @Description:
+   *   Find the shortest distance from the `line` segment to a given `point`
+   *   and assign it to `out`.  Use it for line segments only.
+   *
+   * @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.
+   *
+   * @Output:
+   *   out ::
+   *     Signed distance from `point` to `line`.
+   *
+   * @Return:
+   *   FreeType error, 0 means success.
+   *
+   * @Note:
+   *   The `line' parameter must have an 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, we do the following.  Let's assume that
+     *
+     * ```
+     * a = start point of the line segment
+     * b = end point of the line segment
+     * p = point from which shortest distance is to be calculated
+     * ```
+     *
+     * (1) Write the parametric equation of the line.
+     *
+     *     ```
+     *     point_on_line = a + (b - a) * t   (t is the factor)
+     *     ```
+     *
+     * (2) Find the projection of point `p` on the line.  The projection
+     *     will be perpendicular to the line, which allows us to get the
+     *     solution by making the dot product zero.
+     *
+     *     ```
+     *     (point_on_line - a) . (p - point_on_line) = 0
+     *
+     *                (point_on_line)
+     *      (a) x-------o----------------x (b)
+     *                |_|
+     *                  |
+     *                  |
+     *                 (p)
+     *     ```
+     *
+     * (3) Simplification of the above equation yields the factor of
+     *     `point_on_line`:
+     *
+     *     ```
+     *     t = ((p - a) . (b - a)) / |b - a|^2
+     *     ```
+     *
+     * (4) We clamp factor `t` between [0.0f, 1.0f] because `point_on_line`
+     *     can be outside of the line segment:
+     *
+     *     ```
+     *                                          (point_on_line)
+     *     (a) x------------------------x (b) -----o---
+     *                                           |_|
+     *                                             |
+     *                                             |
+     *                                            (p)
+     *     ```
+     *
+     * (5) Finally, the distance we are interested in is
+     *
+     *     ```
+     *     |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 `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;
+  }
+
+
+  /**************************************************************************
+   *
+   * @Function:
+   *   get_min_distance_conic
+   *
+   * @Description:
+   *   Find the shortest distance from the `conic` Bezier curve to a given
+   *   `point` and assign it to `out`.  Use it for conic/quadratic curves
+   *   only.
+   *
+   * @Input:
+   *   conic ::
+   *     The conic Bezier curve to which the shortest distance is to be
+   *     computed.
+   *
+   *   point ::
+   *     Point from which the shortest distance is to be computed.
+   *
+   * @Output:
+   *   out ::
+   *     Signed distance from `point` to `conic`.
+   *
+   * @Return:
+   *     FreeType error, 0 means success.
+   *
+   * @Note:
+   *   The `conic` parameter must have an edge type of `SDF_EDGE_CONIC`.
+   *
+   */
+
+#if !USE_NEWTON_FOR_CONIC
+
+  /*
+   * The function uses an analytical method to find the shortest distance
+   * which is faster than the Newton-Raphson method, but has underflows at
+   * the moment.  Use Newton's method if you can see artifacts in the SDF.
+   */
+  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 the line segment algorithm.  The
+     * shortest distance is 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 endpoints, because the
+     * perpendicular may not be the shortest.
+     *
+     * Let's assume that
+     * ```
+     * p0 = first endpoint
+     * p1 = control point
+     * p2 = second endpoint
+     * p  = point from which shortest distance is to be calculated
+     * ```
+     *
+     * (1) The equation of a quadratic Bezier curve can be written as
+     *
+     *     ```
+     *     B(t) = (1 - t)^2 * p0 + 2(1 - t)t * p1 + t^2 * p2
+     *     ```
+     *
+     *     with `t` a factor in the range [0.0f, 1.0f].  This equation can
+     *     be rewritten as
+     *
+     *     ```
+     *     B(t) = t^2 * (p0 - 2p1 + p2) + 2t * (p1 - p0) + p0
+     *     ```
+     *
+     *     With
+     *
+     *     ```
+     *     A = p0 - 2p1 + p2
+     *     B = p1 - p0
+     *     ```
+     *
+     *     we have
+     *
+     *     ```
+     *     B(t) = t^2 * A + 2t * B + p0
+     *     ```
+     *
+     * (2) The derivative of the last equation above is
+     *
+     *     ```
+     *     B'(t) = 2 *(tA + B)
+     *     ```
+     *
+     * (3) 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)`) make 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
+     *
+     *     ```
+     *     at^3 + bt^2 + ct + d = 0
+     *     ```
+     *
+     *     with
+     *
+     *     ```
+     *     a = A.A
+     *     b = 3A.B
+     *     c = 2B.B + A.p0 - A.p
+     *     d = p0.B - p.B
+     *     ```
+     *
+     * (4) Now the roots of the equation can be computed using 'Cardano's
+     *     Cubic formula'; we clamp the roots in the 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;
+    }
+
+    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 curve, but Behdad
+       * Esfahbod proved the following lemma.
+       *
+       * Lemma:
+       *
+       * (1) 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.
+       *
+       * (2) 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 more details see
+       *
+       *   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 );
+
+      /* compute `cross` if not perpendicular */
+      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 /* USE_NEWTON_FOR_CONIC */
+
+  /*
+   * 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.
+   */
+  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 a conic curve.  It does not involve solving
+     * any cubic equation, that is why there is no risk of underflow.
+     *
+     * Let's assume that
+     *
+     * ```
+     * p0 = first endpoint
+     * p1 = control point
+     * p3 = second endpoint
+     * p  = point from which shortest distance is to be calculated
+     * ```
+     *
+     * (1) The equation of a quadratic Bezier curve can be written as
+     *
+     *     ```
+     *     B(t) = (1 - t)^2 * p0 + 2(1 - t)t * p1 + t^2 * p2
+     *     ```
+     *
+     *     with `t` the factor in the range [0.0f, 1.0f].  The above
+     *     equation can be rewritten as
+     *
+     *     ```
+     *     B(t) = t^2 * (p0 - 2p1 + p2) + 2t * (p1 - p0) + p0
+     *     ```
+     *
+     *     With
+     *
+     *     ```
+     *     A = p0 - 2p1 + p2
+     *     B = 2 * (p1 - p0)
+     *     ```
+     *
+     *     we have
+     *
+     *     ```
+     *     B(t) = t^2 * A + t * B + p0
+     *     ```
+     *
+     * (2) The derivative of the above equation is
+     *
+     *     ```
+     *     B'(t) = 2t * A + B
+     *     ```
+     *
+     * (3) The second derivative of the above equation is
+     *
+     *     ```
+     *     B''(t) = 2A
+     *     ```
+     *
+     * (4) The equation `P(t)` of the distance from point `p` to the curve
+     *     can be written as
+     *
+     *     ```
+     *     P(t) = t^2 * A + t^2 * B + p0 - p
+     *     ```
+     *
+     *     With
+     *
+     *     ```
+     *     C = p0 - p
+     *     ```
+     *
+     *     we have
+     *
+     *     ```
+     *     P(t) = t^2 * A + t * B + C
+     *     ```
+     *
+     * (5) Finally, the equation of the angle between `B(t)` and `P(t)` can
+     *     be written as
+     *
+     *     ```
+     *     Q(t) = P(t) . B'(t)
+     *     ```
+     *
+     * (6) Our task is to find a value of `t` such that the above equation
+     *     `Q(t)` becomes zero, this is, the point-to-curve vector makes
+     *     90~degrees with the curve.  We solve this with the Newton-Raphson
+     *     method.
+     *
+     * (7) We first assume an arbitary value of factor `t`, which we then
+     *     improve.
+     *
+     *     ```
+     *     t := Q(t) / Q'(t)
+     *     ```
+     *
+     *     Putting the 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))
+     *     ```
+     *
+     *     Note that `P'(t)` is the same as `B'(t)` because the constant is
+     *     gone due to the derivative.
+     *
+     * (8) 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;
+    }
+
+    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 );
+
+        /* 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 is 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) = 2t * A + 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 );
+
+      /* compute `cross` if not perpendicular */
+      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 /* USE_NEWTON_FOR_CONIC */
+
+
+  /**************************************************************************
+   *
+   * @Function:
+   *   get_min_distance_cubic
+   *
+   * @Description:
+   *   Find the shortest distance from the `cubic` Bezier curve to a given
+   *   `point` and assigns it to `out`.  Use it for cubic curves only.
+   *
+   * @Input:
+   *   cubic ::
+   *     The cubic Bezier curve to which the shortest distance is to be
+   *     computed.
+   *
+   *   point ::
+   *     Point from which the shortest distance is to be computed.
+   *
+   * @Output:
+   *   out ::
+   *     Signed distance from `point` to `cubic`.
+   *
+   * @Return:
+   *   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 an 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 quadratic curve algorithm.  The only
+     * difference is that while calculating factor `t`, instead of a cubic
+     * polynomial equation we have to find the roots of a 5th degree
+     * polynomial equation.  Solving this would require a significant amount
+     * of time, and still the results may not be accurate.  We are thus
+     * going to directly approximate the value of `t` using the Newton-Raphson
+     * method.
+     *
+     * Let's assume that
+     *
+     * ```
+     * p0 = first endpoint
+     * p1 = first control point
+     * p2 = second control point
+     * p3 = second endpoint
+     * p  = point from which shortest distance is to be calculated
+     * ```
+     *
+     * (1) The equation of a cubic Bezier curve can be written as
+     *
+     *     ```
+     *     B(t) = (1 - t)^3 * p0 + 3(1 - t)^2 t * p1 +
+     *              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) +
+     *              3t^2 * (p0 - 2p1 + p2) + 3t * (-p0 + p1) + p0
+     *     ```
+     *
+     *     With
+     *
+     *     ```
+     *     A = -p0 + 3p1 - 3p2 + p3
+     *     B = 3(p0 - 2p1 + p2)
+     *     C = 3(-p0 + p1)
+     *     ```
+     *
+     *     we have
+     *
+     *     ```
+     *     B(t) = t^3 * A + t^2 * B + t * C + p0
+     *     ```
+     *
+     * (2) The derivative of the above equation is
+     *
+     *     ```
+     *     B'(t) = 3t^2 * A + 2t * B + C
+     *     ```
+     *
+     * (3) The second derivative of the above equation is
+     *
+     *     ```
+     *     B''(t) = 6t * A + 2B
+     *     ```
+     *
+     * (4) The equation `P(t)` of the distance from point `p` to the curve
+     *     can be written as
+     *
+     *     ```
+     *     P(t) = t^3 * A + t^2 * B + t * C + p0 - p
+     *     ```
+     *
+     *     With
+     *
+     *     ```
+     *     D = p0 - p
+     *     ```
+     *
+     *     we have
+     *
+     *     ```
+     *     P(t) = t^3 * A + t^2 * B + t * C + D
+     *     ```
+     *
+     * (5) Finally the equation of the angle between `B(t)` and `P(t)` can
+     *     be written as
+     *
+     *     ```
+     *     Q(t) = P(t) . B'(t)
+     *     ```
+     *
+     * (6) Our task is to find a value of `t` such that the above equation
+     *     `Q(t)` becomes zero, this is, the point-to-curve vector makes
+     *     90~degree with curve.  We solve this with the Newton-Raphson
+     *     method.
+     *
+     * (7) We first assume an arbitary value of factor `t`, which we then
+     *     improve.
+     *
+     *     ```
+     *     t := Q(t) / Q'(t)
+     *     ```
+     *
+     *     Putting the 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))
+     *     ```
+     *
+     *     Note that `P'(t)` is the same as `B'(t)` because the constant is
+     *     gone due to the derivative.
+     *
+     * (8) 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_factor    = 0;            /* factor at shortest distance */
+    FT_16D16  min_factor_sq = 0;            /* factor at shortest distance */
+    FT_16D16  cross;                        /* to determine the sign       */
+    FT_16D16  min           = FT_INT_MAX;   /* shortest distance           */
+
+    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;
+    }
+
+    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 + t * C + 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 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 );
+
+      /* compute `cross` if not perpendicular */
+      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 */