[gtk/path-stroke: 12/13] Generalize handling of degenerate curves

[Matthias Clasen <matthiasc src gnome org>]

commit c58e9d2b8b1e169558afb021dbff3fd2079fa892
Author: Matthias Clasen <mclasen redhat com>
Date:   Sat Dec 19 22:27:15 2020 -0500

    Generalize handling of degenerate curves
    
    Handle all cusps in the same way, including
    the previously handled degenerate cases.
    
    This makes both
    M 100 100 C 300 300 100 300 300 100
    and
    M 100 100 C 200 100 300 100 200 100
    come up with the expected round joins at
    the cusps.

 gsk/gskpathstroke.c | 260 +++++++++++++++++++++-------------------------------
 1 file changed, 104 insertions(+), 156 deletions(-)
---
diff --git a/gsk/gskpathstroke.c b/gsk/gskpathstroke.c
index 02ef7ea863..6225c51bf6 100644
--- a/gsk/gskpathstroke.c
+++ b/gsk/gskpathstroke.c
@@ -205,54 +205,29 @@ cubic_is_simple (const GskCurve *curve)
   const graphene_point_t *pts = curve->curve.points;
   float a1, a2, s;
   graphene_vec2_t t1, t2, t3;
-  graphene_vec2_t n1, n2;
-
-  get_tangent (&pts[0], &pts[1], &t1);
-  get_tangent (&pts[1], &pts[2], &t2);
-  get_tangent (&pts[2], &pts[3], &t3);
-  a1 = angle_between (&t1, &t2);
-  a2 = angle_between (&t2, &t3);
-
-  if ((a1 < 0 && a2 > 0) || (a1 > 0 && a2 < 0))
-    return FALSE;
-
-  get_normal (&pts[0], &pts[1], &n1);
-  get_normal (&pts[2], &pts[3], &n2);
-
-  s = graphene_vec2_dot (&n1, &n2);
-
-  if (fabs (acos (s)) >= M_PI / 3.f)
-    return FALSE;
-
-  return TRUE;
-}
 
-static gboolean
-curve_is_degenerate (const GskCurve *curve)
-{
-  if (curve->op == GSK_PATH_CURVE)
+  if (!graphene_point_near (&pts[0], &pts[1], 0.001) &&
+      !graphene_point_near (&pts[1], &pts[2], 0.001) &&
+      !graphene_point_near (&pts[2], &pts[3], 0.001))
     {
-      const graphene_point_t *pts = curve->curve.points;
-      float a1, a2;
-      graphene_vec2_t t1, t2, t3;
-
       get_tangent (&pts[0], &pts[1], &t1);
       get_tangent (&pts[1], &pts[2], &t2);
       get_tangent (&pts[2], &pts[3], &t3);
       a1 = angle_between (&t1, &t2);
       a2 = angle_between (&t2, &t3);
 
-      if (a1 < 0)
-        a1 += M_PI;
-      if (a2 < 0)
-        a2 += M_PI;
-
-      if (MIN (fabs (M_PI - a1), fabs (a1)) < DEG_TO_RAD (3) &&
-          MIN (fabs (M_PI - a2), fabs (a2)) < DEG_TO_RAD (3))
-        return TRUE;
+      if ((a1 < 0 && a2 > 0) || (a1 > 0 && a2 < 0))
+        return FALSE;
     }
 
-  return FALSE;
+  gsk_curve_get_start_tangent (curve, &t1);
+  gsk_curve_get_end_tangent (curve, &t2);
+  s = graphene_vec2_dot (&t1, &t2);
+
+  if (fabs (acos (s)) >= M_PI / 3.f)
+    return FALSE;
+
+  return TRUE;
 }
 
 static void
@@ -277,17 +252,13 @@ align_points (const graphene_point_t *p,
     }
 }
 
+/* find solutions for at^2 + bt + c = 0 in the interval [0,1] */
 static int
-get_cubic_extrema (float pa, float pb, float pc, float pd, float t[2])
+solve_quadratic (float a, float b, float c, float t[2])
 {
-  float a, b, c;
   float d, tt;
   int n = 0;
 
-  a = 3 * (pd - 3*pc + 3*pb - pa);
-  b = 6 * (pc - 2*pb + pa);
-  c = 3 * (pb - pa);
-
   if (fabs (a) > 0.0001)
     {
       if (b*b > 4*a*c)
@@ -317,24 +288,6 @@ get_cubic_extrema (float pa, float pb, float pc, float pd, float t[2])
   return n;
 }
 
-/* For degenerate curves, find interior points
- * where the function has a cusp. We do that by
- * aligning the control points to be all on the
- * x axis, and then looking for the t where x
- * has the min/max value.
- */
-static int
-find_turning_points (const GskCurve *curve,
-                     float     t[2])
-{
-  const graphene_point_t *pts = curve->curve.points;
-  graphene_point_t p[4];
-
-  align_points (pts, &pts[0], &pts[1], p, 4);
-
-  return get_cubic_extrema (p[0].x, p[1].x, p[2].x, p[3].x, t);
-}
-
 /* Get the points where the curvature of curve is
  * zero, or a maximum or minimum, inside the open
  * interval from 0 to 1.
@@ -347,8 +300,6 @@ cubic_curvature_points (const GskCurve *curve,
   graphene_point_t p[4];
   float a, b, c, d;
   float x, y, z;
-  float u2, u, tt;
-  int n_roots = 0;
 
   align_points (pts, &pts[0], &pts[3], p, 4);
 
@@ -361,28 +312,53 @@ cubic_curvature_points (const GskCurve *curve,
   y = 3*a - b - 3*c;
   z = c - a;
 
-  if (fabs (x) >= 0.001)
-    {
-      tt = -y / (2*x);
-      if (0 < tt && tt < 1)
-        t[n_roots++] = tt;
+  return solve_quadratic (x, y, z, t);
+}
 
-      u2 = y*y - 4*x*z;
-      if (u2 > 0.001)
-        {
-          u = sqrt (u2);
+/* Find cusps in the interior of [0,1]. According to Stone &
+ * dRose, A Geometric Characterization of Parametric Cubic
+ * curves, a necessary and sufficient condition is that the
+ * first derivative vanishes.
+ */
+static int
+find_cusps (const GskCurve *curve,
+            float           t[2])
+{
+  const graphene_point_t *pts = curve->curve.points;
+  graphene_point_t p[3];
+  float ax, bx, cx;
+  float ay, by, cy;
+  float tx[2];
+  int nx;
+  int n = 0;
 
-          tt = (-y + u) / (2*x);
-          if (0 < tt && tt < 1)
-            t[n_roots++] = tt;
+  p[0].x = 3 * (pts[1].x - pts[0].x);
+  p[0].y = 3 * (pts[1].y - pts[0].y);
+  p[1].x = 3 * (pts[2].x - pts[1].x);
+  p[1].y = 3 * (pts[2].y - pts[1].y);
+  p[2].x = 3 * (pts[3].x - pts[2].x);
+  p[2].y = 3 * (pts[3].y - pts[2].y);
 
-          tt = (-y - u) / (2*x);
-          if (0 < tt && tt < 1)
-            t[n_roots++] = tt;
-        }
+  ax = p[0].x - 2 * p[1].x + p[2].x;
+  bx = - 2 * p[0].x + 2 * p[1].x;
+  cx = p[0].x;
+
+  nx = solve_quadratic (ax, bx, cx, tx);
+
+  ay = p[0].y - 2 * p[1].y + p[2].y;
+  by = - 2 * p[0].y + 2 * p[1].y;
+  cy = p[0].y;
+
+  for (int i = 0; i < nx; i++)
+    {
+      float ti = tx[i];
+
+      if (0 < ti && ti < 1 &&
+          fabs (ay * ti * ti + by * ti + cy) < 0.001)
+        t[n++] = ti;
     }
 
-  return n_roots;
+  return n;
 }
 
 /* }}} */
@@ -553,6 +529,7 @@ typedef struct
   gboolean has_current_point;  // r0, l0 have been set from a move
   gboolean has_current_curve;  // c, l, r are set from a curve
   gboolean is_first_curve; // if c, l, r are the first segments we've seen
+  gboolean force_round_join; // force next join to be round
 
   GskCurve c;  // previous segment of the path
   GskCurve l;  // candidate for left contour of c
@@ -737,8 +714,27 @@ add_debug (StrokeData     *stroke_data,
           break;
         case GSK_PATH_CURVE:
           p = curve->curve.points;
+          path_builder_move_to_point (stroke_data->debug, &p[0]);
+          path_builder_line_to_point (stroke_data->debug, &p[1]);
+          path_builder_line_to_point (stroke_data->debug, &p[2]);
+          path_builder_line_to_point (stroke_data->debug, &p[3]);
           gsk_path_builder_add_circle (stroke_data->debug, &p[0], 3);
+          gsk_path_builder_add_circle (stroke_data->debug, &p[1], 2);
+          gsk_path_builder_add_circle (stroke_data->debug, &p[2], 2);
           gsk_path_builder_add_circle (stroke_data->debug, &p[3], 3);
+
+          {
+            float t[2];
+            int n;
+            graphene_point_t q;
+
+            n = find_cusps (curve, t);
+            for (int i = 0; i < n; i++)
+              {
+                gsk_curve_get_point (curve, t[i], &q);
+                gsk_path_builder_add_circle (stroke_data->debug, &q, 5);
+              }
+          }
           break;
         case GSK_PATH_CONIC:
           p = curve->conic.points;
@@ -771,23 +767,8 @@ add_curve (StrokeData     *stroke_data,
   gsk_curve_offset (curve, - stroke_data->stroke->line_width / 2, &l);
   gsk_curve_offset (curve, stroke_data->stroke->line_width / 2, &r);
 
-  add_segments (stroke_data, curve, &r, &l, FALSE);
-}
-
-static void
-add_curve_round_join (StrokeData     *stroke_data,
-                      const GskCurve *curve)
-{
-  GskCurve l, r;
-
-#ifdef STROKE_DEBUG
-  add_debug (stroke_data, curve);
-#endif
-
-  gsk_curve_offset (curve, - stroke_data->stroke->line_width / 2, &l);
-  gsk_curve_offset (curve, stroke_data->stroke->line_width / 2, &r);
-
-  add_segments (stroke_data, curve, &r, &l, TRUE);
+  add_segments (stroke_data, curve, &r, &l, stroke_data->force_round_join);
+  stroke_data->force_round_join = FALSE;
 }
 
 static int
@@ -798,52 +779,6 @@ cmpfloat (const void *p1, const void *p2)
   return *f1 < *f2 ? -1 : (*f1 > *f2 ? 1 : 0);
 }
 
-static void
-add_degenerate_curve (StrokeData     *stroke_data,
-                      const GskCurve *curve)
-{
-  float t[3];
-  int n;
-  GskCurve c;
-  graphene_point_t p[2];
-
-  n = find_turning_points (curve, t);
-
-  if (n == 0)
-    add_curve (stroke_data, curve);
-  else if (n == 1)
-    {
-      p[0] = *gsk_curve_get_start_point (curve);
-      gsk_curve_get_point (curve, t[0], &p[1]);
-      gsk_curve_init_foreach (&c, GSK_PATH_LINE, p, 2, 0);
-      add_curve (stroke_data, &c);
-
-      p[0] = p[1];
-      p[1] = *gsk_curve_get_end_point (curve);
-      gsk_curve_init_foreach (&c, GSK_PATH_LINE, p, 2, 0);
-      add_curve_round_join (stroke_data, &c);
-    }
-  else
-    {
-      qsort (t, n, sizeof (float), cmpfloat);
-
-      p[0] = *gsk_curve_get_start_point (curve);
-      gsk_curve_get_point (curve, t[0], &p[1]);
-      gsk_curve_init_foreach (&c, GSK_PATH_LINE, p, 2, 0);
-      add_curve (stroke_data, &c);
-
-      p[0] = p[1];
-      gsk_curve_get_point (curve, t[1], &p[1]);
-      gsk_curve_init_foreach (&c, GSK_PATH_LINE, p, 2, 0);
-      add_curve_round_join (stroke_data, &c);
-
-      p[0] = p[1];
-      p[1] = *gsk_curve_get_end_point (curve);
-      gsk_curve_init_foreach (&c, GSK_PATH_LINE, p, 2, 0);
-      add_curve_round_join (stroke_data, &c);
-    }
-}
-
 static void
 add_degenerate_conic (StrokeData     *stroke_data,
                       const GskCurve *curve)
@@ -859,7 +794,8 @@ add_degenerate_conic (StrokeData     *stroke_data,
   p[0] = p[1];
   p[1] = *gsk_curve_get_end_point (curve);
   gsk_curve_init_foreach (&c, GSK_PATH_LINE, p, 2, 0);
-  add_curve_round_join (stroke_data, &c);
+  stroke_data->force_round_join = TRUE;
+  add_curve (stroke_data, &c);
 }
 
 #define MAX_SUBDIVISION 8
@@ -869,15 +805,30 @@ subdivide_and_add_curve (StrokeData     *stroke_data,
                          const GskCurve *curve,
                          int             level)
 {
-  if (level == MAX_SUBDIVISION && curve_is_degenerate (curve))
-    add_degenerate_curve (stroke_data, curve);
-  else if (level == 0 || (level < MAX_SUBDIVISION && cubic_is_simple (curve)))
+  GskCurve c1, c2;
+  float t[5];
+  int n;
+
+  if (level == 0)
+    add_curve (stroke_data, curve);
+  else if (level < MAX_SUBDIVISION && cubic_is_simple (curve))
     add_curve (stroke_data, curve);
+  else if (level == MAX_SUBDIVISION && (n = find_cusps (curve, t)) > 0)
+    {
+      t[n++] = 0;
+      t[n++] = 1;
+      qsort (t, n, sizeof (float), cmpfloat);
+      for (int i = 0; i + 1 < n; i++)
+        {
+          gsk_curve_segment (curve, t[i], t[i + 1], &c1);
+          if (i > 0)
+            stroke_data->force_round_join = TRUE;
+          subdivide_and_add_curve (stroke_data, &c1, level - 1);
+        }
+    }
   else
     {
-      float t[5];
-      int n = 0;
-
+      n = 0;
       t[n++] = 0;
       t[n++] = 1;
 
@@ -889,19 +840,16 @@ subdivide_and_add_curve (StrokeData     *stroke_data,
 
       if (n == 2)
         {
-          GskCurve c1, c2;
-
           gsk_curve_split (curve, 0.5, &c1, &c2);
           subdivide_and_add_curve (stroke_data, &c1, level - 1);
           subdivide_and_add_curve (stroke_data, &c2, level - 1);
         }
       else
         {
-          GskCurve c;
           for (int i = 0; i + 1 < n; i++)
             {
-              gsk_curve_segment (curve, t[i], t[i+1], &c);
-              subdivide_and_add_curve (stroke_data, &c, level - 1);
+              gsk_curve_segment (curve, t[i], t[i+1], &c1);
+              subdivide_and_add_curve (stroke_data, &c1, level - 1);
             }
         }
     }