[gtk/curve-ops: 1/5] Add gsk_curve_intersect




commit bd95d0566e5d6df4356fb5b73254bcb0bdc03a7a
Author: Matthias Clasen <mclasen redhat com>
Date:   Sun Dec 6 22:09:48 2020 -0500

    Add gsk_curve_intersect
    
    Add a way to find the intersections of two curves.
    This will be used in stroking.

 gsk/gskcurve.c        | 420 +++++++++++++++++++++++++++++++++++++++++++++++++-
 gsk/gskcurveprivate.h |   6 +
 2 files changed, 424 insertions(+), 2 deletions(-)
---
diff --git a/gsk/gskcurve.c b/gsk/gskcurve.c
index c5feef335b..c0c36c0c7f 100644
--- a/gsk/gskcurve.c
+++ b/gsk/gskcurve.c
@@ -525,14 +525,104 @@ gsk_conic_curve_eval (const GskCurve   *curve,
     }
 }
 
+static void
+split_bezier3d_recurse (const graphene_point3d_t *p,
+                        int                       l,
+                        float                     t,
+                        graphene_point3d_t       *left,
+                        graphene_point3d_t       *right,
+                        int                      *lpos,
+                        int                      *rpos)
+{
+  if (l == 1)
+    {
+      left[*lpos] = p[0];
+      right[*rpos] = p[0];
+    }
+  else
+    {
+      graphene_point3d_t *np;
+      int i;
+
+      np = g_alloca (sizeof (graphene_point3d_t) * (l - 1));
+      for (i = 0; i < l - 1; i++)
+        {
+          if (i == 0)
+            {
+              left[*lpos] = p[i];
+              (*lpos)++;
+            }
+          if (i + 1 == l - 1)
+            {
+              right[*rpos] = p[i + 1];
+              (*rpos)--;
+            }
+          graphene_point3d_interpolate (&p[i], &p[i + 1], t, &np[i]);
+        }
+      split_bezier3d_recurse (np, l - 1, t, left, right, lpos, rpos);
+    }
+}
+
+static void
+split_bezier3d (const graphene_point3d_t *p,
+                int                       l,
+                float                     t,
+                graphene_point3d_t       *left,
+                graphene_point3d_t       *right)
+{
+  int lpos = 0;
+  int rpos = l - 1;
+  split_bezier3d_recurse (p, l, t, left, right, &lpos, &rpos);
+}
+
 static void
 gsk_conic_curve_split (const GskCurve   *curve,
                        float             progress,
                        GskCurve         *start,
                        GskCurve         *end)
 {
-  //const GskConicCurve *self = &curve->conic;
-  g_warning ("FIXME: Stop treating conics as lines");
+  const GskConicCurve *self = &curve->conic;
+  graphene_point3d_t p[3];
+  graphene_point3d_t l[3], r[3];
+  graphene_point_t left[4], right[4];
+  float w;
+
+  /* do de Casteljau in homogeneous coordinates... */
+  w = self->points[2].x;
+  p[0] = GRAPHENE_POINT3D_INIT (self->points[0].x, self->points[0].y, 1);
+  p[1] = GRAPHENE_POINT3D_INIT (self->points[1].x * w, self->points[1].y * w, w);
+  p[2] = GRAPHENE_POINT3D_INIT (self->points[3].x, self->points[3].y, 1);
+
+  split_bezier3d (p, 3, progress, l, r);
+
+  /* then project the control points down */
+  left[0] = GRAPHENE_POINT_INIT (l[0].x / l[0].z, l[0].y / l[0].z);
+  left[1] = GRAPHENE_POINT_INIT (l[1].x / l[1].z, l[1].y / l[1].z);
+  left[3] = GRAPHENE_POINT_INIT (l[2].x / l[2].z, l[2].y / l[2].z);
+
+  right[0] = GRAPHENE_POINT_INIT (r[0].x / r[0].z, r[0].y / r[0].z);
+  right[1] = GRAPHENE_POINT_INIT (r[1].x / r[1].z, r[1].y / r[1].z);
+  right[3] = GRAPHENE_POINT_INIT (r[2].x / r[2].z, r[2].y / r[2].z);
+
+  /* normalize the outer weights to be 1 by using
+   * the fact that weights w_i and c*w_i are equivalent
+   * for any nonzero constant c
+   */
+  for (int i = 0; i < 3; i++)
+    {
+      l[i].z /= l[0].z;
+      r[i].z /= r[2].z;
+    }
+
+  /* normalize the inner weight to be 1 by using
+   * the fact that w_0*w_2/w_1^2 is a constant for
+   * all equivalent weights.
+   */
+  left[2] = GRAPHENE_POINT_INIT (l[1].z / sqrt (l[2].z), 0);
+  right[2] = GRAPHENE_POINT_INIT (r[1].z / sqrt (r[0].z), 0);
+
+  gsk_curve_init (start, gsk_pathop_encode (GSK_PATH_CONIC, left));
+  gsk_curve_init (end, gsk_pathop_encode (GSK_PATH_CONIC, right));
 }
 
 /* taken from Skia, including the very descriptive name */
@@ -772,3 +862,329 @@ gsk_curve_get_bounds (const GskCurve  *curve,
 {
   get_class (curve->op)->get_bounds (curve, bounds);
 }
+
+/** Intersections **/
+
+static int
+line_intersect (const GskCurve   *curve1,
+                const GskCurve   *curve2,
+                float            *t1,
+                float            *t2,
+                graphene_point_t *p)
+{
+  const graphene_point_t *pts1 = curve1->line.points;
+  const graphene_point_t *pts2 = curve2->line.points;
+  float a1 = pts1[0].x - pts1[1].x;
+  float b1 = pts1[0].y - pts1[1].y;
+  float a2 = pts2[0].x - pts2[1].x;
+  float b2 = pts2[0].y - pts2[1].y;
+  float det = a1 * b2 - b1 * a2;
+
+  if (det != 0)
+    {
+      float tt =   ((pts1[0].x - pts2[0].x) * b2 - (pts1[0].y - pts2[0].y) * a2) / det;
+      float ss = - ((pts1[0].y - pts2[0].y) * a1 - (pts1[0].x - pts2[0].x) * b1) / det;
+
+      if (acceptable (tt) && acceptable (ss))
+        {
+          p->x = pts1[0].x + tt * (pts1[1].x - pts1[0].x);
+          p->y = pts1[0].y + tt * (pts1[1].y - pts1[0].y);
+
+          *t1 = tt;
+          *t2 = ss;
+
+          return 1;
+        }
+    }
+
+  return 0;
+}
+
+static void
+direction_vector (const graphene_point_t *p0,
+                  const graphene_point_t *p1,
+                  graphene_vec2_t        *t)
+{
+  graphene_vec2_init (t, p1->x - p0->x, p1->y - p0->y);
+  graphene_vec2_normalize (t, t);
+}
+
+static void
+align_points (const graphene_point_t *p,
+              const graphene_point_t *a,
+              const graphene_point_t *b,
+              graphene_point_t       *q,
+              int                     n)
+{
+  graphene_vec2_t n1;
+  float angle;
+
+  direction_vector (a, b, &n1);
+  angle = -atan2 (graphene_vec2_get_y (&n1), graphene_vec2_get_x (&n1));
+
+  for (int i = 0; i < n; i++)
+    {
+      q[i].x = (p[i].x - a->x) * cos (angle) - (p[i].y - a->y) * sin (angle);
+      q[i].y = (p[i].x - a->x) * sin (angle) + (p[i].y - a->y) * cos (angle);
+    }
+}
+
+static void
+find_point_on_line (const graphene_point_t *p1,
+                    const graphene_point_t *p2,
+                    const graphene_point_t *q,
+                    float                  *t)
+{
+  float tx = p2->x - p1->x;
+  float ty = p2->y - p1->y;
+  float sx = q->x - p1->x;
+  float sy = q->y - p1->y;
+
+  *t = (tx*sx + ty*sy) / (tx*tx + ty*ty);
+}
+
+static float
+cuberoot (float v)
+{
+  if (v < 0)
+    return -pow (-v, 1.f / 3);
+  return pow (v, 1.f / 3);
+}
+
+/* Solve P = 0 where P is
+ * P = (1-t)^3*pa + 3*t*(1-t)^2*pb + 3*t^2*(1-t)*pc + t^3*pd
+ */
+static int
+get_cubic_roots (float pa, float pb, float pc, float pd, float roots[3])
+{
+  float a, b, c, d;
+  float q, q2;
+  float p, p3;
+  float discriminant;
+  float u1, v1, sd;
+  int n_roots = 0;
+
+  d = -pa + 3*pb - 3*pc + pd;
+  a = 3*pa - 6*pb + 3*pc;
+  b = -3*pa + 3*pb;
+  c = pa;
+
+  if (fabs (d) < 0.0001)
+    {
+      if (fabs (a) < 0.0001)
+        {
+          if (fabs (b) < 0.0001)
+            return 0;
+
+          if (acceptable (-c / b))
+            {
+              roots[0] = -c / b;
+
+              return 1;
+            }
+
+          return 0;
+        }
+      q = sqrt (b*b - 4*a*c);
+      roots[n_roots] = (-b + q) / (2 * a);
+      if (acceptable (roots[n_roots]))
+        n_roots++;
+
+      roots[n_roots] = (-b - q) / (2 * a);
+      if (acceptable (roots[n_roots]))
+        n_roots++;
+
+      return n_roots;
+    }
+
+  a /= d;
+  b /= d;
+  c /= d;
+
+  p = (3*b - a*a)/3;
+  p3 = p/3;
+  q = (2*a*a*a - 9*a*b + 27*c)/27;
+  q2 = q/2;
+  discriminant = q2*q2 + p3*p3*p3;
+
+  if (discriminant < 0)
+    {
+      float mp3 = -p/3;
+      float mp33 = mp3*mp3*mp3;
+      float r = sqrt (mp33);
+      float t = -q / (2*r);
+      float cosphi = t < -1 ? -1 : (t > 1 ? 1 : t);
+      float phi = acos (cosphi);
+      float crtr = cuberoot (r);
+      float t1 = 2*crtr;
+
+      roots[n_roots] = t1 * cos (phi/3) - a/3;
+      if (acceptable (roots[n_roots]))
+        n_roots++;
+      roots[n_roots] = t1 * cos ((phi + 2*M_PI) / 3) - a/3;
+      if (acceptable (roots[n_roots]))
+        n_roots++;
+      roots[n_roots] = t1 * cos ((phi + 4*M_PI) / 3) - a/3;
+      if (acceptable (roots[n_roots]))
+        n_roots++;
+
+    return n_roots;
+  }
+
+  if (discriminant == 0)
+    {
+      u1 = q2 < 0 ? cuberoot (-q2) : -cuberoot (q2);
+      roots[n_roots] = 2*u1 - a/3;
+      if (acceptable (roots[n_roots]))
+        n_roots++;
+      roots[n_roots] = -u1 - a/3;
+      if (acceptable (roots[n_roots]))
+        n_roots++;
+
+      return n_roots;
+    }
+
+  sd = sqrt (discriminant);
+  u1 = cuberoot (sd - q2);
+  v1 = cuberoot (sd + q2);
+  roots[n_roots] = u1 - v1 - a/3;
+  if (acceptable (roots[n_roots]))
+    n_roots++;
+
+  return n_roots;
+}
+
+static int
+line_curve_intersect (const GskCurve   *curve1,
+                      const GskCurve   *curve2,
+                      float            *t1,
+                      float            *t2,
+                      graphene_point_t *p,
+                      int               n)
+{
+  const graphene_point_t *a = &curve1->line.points[0];
+  const graphene_point_t *b = &curve1->line.points[1];
+  graphene_point_t pts[4];
+  float t[3];
+  int m, i;
+
+  /* Rotate things to place curve1 on the x axis,
+   * then solve curve2 for y == 0.
+   */
+  align_points (curve2->curve.points, a, b, pts, 4);
+
+  m = get_cubic_roots (pts[0].y, pts[1].y, pts[2].y, pts[3].y, t);
+
+  m = MIN (m, n);
+  for (i = 0; i < m; i++)
+    {
+      t2[i] = t[i];
+      gsk_curve_eval (curve2, t[i], &p[i], NULL);
+      find_point_on_line (a, b, &p[i], &t1[i]);
+    }
+
+  return m;
+}
+
+static void
+curve_intersect_recurse (const GskCurve   *curve1,
+                         const GskCurve   *curve2,
+                         float             t1l,
+                         float             t1r,
+                         float             t2l,
+                         float             t2r,
+                         float            *t1,
+                         float            *t2,
+                         graphene_point_t *p,
+                         int               n,
+                         int              *pos)
+{
+  graphene_rect_t b1, b2;
+  float d1, d2;
+  GskCurve p11, p12, p21, p22;
+
+  if (*pos == n)
+    return;
+
+  gsk_curve_get_bounds (curve1, &b1);
+  gsk_curve_get_bounds (curve2, &b2);
+
+  if (!graphene_rect_intersection (&b1, &b2, NULL))
+    return;
+
+  d1 = (t1r - t1l) / 2;
+  d2 = (t2r - t2l) / 2;
+
+  if (b1.size.width < 0.1 && b1.size.height < 0.1 &&
+      b2.size.width < 0.1 && b2.size.height < 0.1)
+    {
+      graphene_point_t c;
+      t1[*pos] = t1l + d1;
+      t2[*pos] = t2l + d2;
+      gsk_curve_eval (curve1, 0.5, &c, NULL);
+
+      for (int i = 0; i < *pos; i++)
+        {
+          if (graphene_point_near (&c, &p[i], 0.1))
+            return;
+        }
+
+      p[*pos] = c;
+      (*pos)++;
+
+      return;
+    }
+
+  gsk_curve_split (curve1, 0.5, &p11, &p12);
+  gsk_curve_split (curve2, 0.5, &p21, &p22);
+
+  curve_intersect_recurse (&p11, &p21, t1l,      t1l + d1, t2l,      t2l + d2, t1, t2, p, n, pos);
+  curve_intersect_recurse (&p11, &p22, t1l,      t1l + d1, t2l + d2, t2r,      t1, t2, p, n, pos);
+  curve_intersect_recurse (&p12, &p21, t1l + d1, t1r,      t2l,      t2l + d2, t1, t2, p, n, pos);
+  curve_intersect_recurse (&p12, &p22, t1l + d1, t1r,      t2l + d2, t2r,      t1, t2, p, n, pos);
+}
+
+static int
+curve_intersect (const GskCurve   *curve1,
+                 const GskCurve   *curve2,
+                 float            *t1,
+                 float            *t2,
+                 graphene_point_t *p,
+                 int               n)
+{
+  int pos = 0;
+
+  curve_intersect_recurse (curve1, curve2, 0, 1, 0, 1, t1, t2, p, n, &pos);
+
+  return pos;
+}
+
+/* Place intersections between the curves in p, and their Bezier positions
+ * in t1 and t2, up to n. Return the number of intersections found.
+ *
+ * Note that two cubic Beziers can have up to 9 intersections.
+ */
+int
+gsk_curve_intersect (const GskCurve   *curve1,
+                     const GskCurve   *curve2,
+                     float            *t1,
+                     float            *t2,
+                     graphene_point_t *p,
+                     int               n)
+{
+  const GskCurveClass *c1 = get_class (curve1->op);
+  const GskCurveClass *c2 = get_class (curve2->op);
+
+  /* We special-case line-line and line-curve intersections,
+   * since we can solve them directly.
+   * Everything else is done via bisection.
+   */
+  if (c1 == &GSK_LINE_CURVE_CLASS && c2 == &GSK_LINE_CURVE_CLASS)
+    return line_intersect (curve1, curve2, t1, t2, p);
+  else if (c1 == &GSK_LINE_CURVE_CLASS && c2 == &GSK_CURVE_CURVE_CLASS)
+    return line_curve_intersect (curve1, curve2, t1, t2, p, n);
+  else if (c1 == &GSK_CURVE_CURVE_CLASS && c2 == &GSK_LINE_CURVE_CLASS)
+    return line_curve_intersect (curve2, curve1, t2, t1, p, n);
+  else
+    return curve_intersect (curve1, curve2, t1, t2, p, n);
+}
diff --git a/gsk/gskcurveprivate.h b/gsk/gskcurveprivate.h
index 746e735ecf..dfdd45b034 100644
--- a/gsk/gskcurveprivate.h
+++ b/gsk/gskcurveprivate.h
@@ -102,6 +102,12 @@ const graphene_point_t *gsk_curve_get_end_point                 (const GskCurve
 void                    gsk_curve_get_bounds                    (const GskCurve         *curve,
                                                                  graphene_rect_t        *bounds);
 
+int                     gsk_curve_intersect                     (const GskCurve         *curve1,
+                                                                 const GskCurve         *curve2,
+                                                                 float                  *t1,
+                                                                 float                  *t2,
+                                                                 graphene_point_t       *p,
+                                                                 int                     n);
 
 G_END_DECLS
 


[Date Prev][Date Next]   [Thread Prev][Thread Next]   [Thread Index] [Date Index] [Author Index]