[gtk/curve-ops: 18/22] Add gsk_curve_intersect
- From: Matthias Clasen <matthiasc src gnome org>
- To: commits-list gnome org
- Cc:
- Subject: [gtk/curve-ops: 18/22] Add gsk_curve_intersect
- Date: Wed, 9 Dec 2020 05:31:31 +0000 (UTC)
commit ba9145727ffe4f2439ffdebdbba1f1f2a1215f0d
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 | 464 ++++++++++++++++++++++++++++++++++++++++++++++++++
gsk/gskcurveprivate.h | 6 +
2 files changed, 470 insertions(+)
---
diff --git a/gsk/gskcurve.c b/gsk/gskcurve.c
index eb6317b060..2987a293bd 100644
--- a/gsk/gskcurve.c
+++ b/gsk/gskcurve.c
@@ -1092,3 +1092,467 @@ gsk_curve_get_tight_bounds (const GskCurve *curve,
{
get_class (curve->op)->get_tight_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
+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;
+
+ get_tangent (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_get_point (curve2, t[i], &p[i]);
+ 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)
+{
+ GskCurve p11, p12, p21, p22;
+ graphene_rect_t b1, b2;
+ float d1, d2;
+
+ if (*pos == n)
+ return;
+
+ gsk_curve_get_tight_bounds (curve1, &b1);
+ gsk_curve_get_tight_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_get_point (curve1, 0.5, &c);
+
+ 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;
+}
+
+static void
+get_bounds (const GskCurve *curve,
+ float tl,
+ float tr,
+ graphene_rect_t *bounds)
+{
+ if (curve->op == GSK_PATH_CONIC)
+ {
+ const graphene_point_t *pts = curve->conic.points;
+ graphene_point3d_t c[3], l[3], r[3], rest[3];
+ graphene_point3d_t *cc;
+ graphene_point_t p[4];
+ GskCurve curve1;
+ float w;
+
+ w = pts[2].x;
+ c[0] = GRAPHENE_POINT3D_INIT (pts[0].x, pts[0].y, 1);
+ c[1] = GRAPHENE_POINT3D_INIT (pts[1].x * w, pts[1].y * w, w);
+ c[2] = GRAPHENE_POINT3D_INIT (pts[3].x, pts[3].y, 1);
+
+ cc = c;
+
+ if (tl > 0)
+ {
+ split_bezier3d (cc, 3, tl, l, r);
+ cc = r;
+ }
+
+ if (tr < 1)
+ {
+ split_bezier3d (cc, 3, (tr - tl) / (1 - tl), l, rest);
+ cc = l;
+ }
+
+ p[0] = GRAPHENE_POINT_INIT (cc[0].x / cc[0].z, cc[0].y / cc[0].z);
+ p[1] = GRAPHENE_POINT_INIT (cc[1].x / cc[1].z, cc[1].y / cc[1].z);
+ p[3] = GRAPHENE_POINT_INIT (cc[2].x / cc[2].z, cc[2].y / cc[2].z);
+
+ for (int i = 0; i < 3; i++)
+ cc[i].z /= cc[0].z;
+
+ p[2] = GRAPHENE_POINT_INIT (cc[1].z / sqrt (cc[2].z), 0);
+
+ gsk_curve_init (&curve1, gsk_pathop_encode (GSK_PATH_CONIC, p));
+ gsk_curve_get_tight_bounds (&curve1, bounds);
+ }
+ else
+ {
+ GskCurve c1, c2, c3;
+ const GskCurve *c;
+
+ c = curve;
+ if (tl > 0)
+ {
+ gsk_curve_split (c, tl, &c1, &c2);
+ c = &c2;
+ }
+
+ if (tr < 1)
+ {
+ gsk_curve_split (c, (tr - tl) / (1 - tl), &c3, &c1);
+ c = &c3;
+ }
+ gsk_curve_get_tight_bounds (c, bounds);
+ }
+}
+
+static void
+general_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;
+
+ if (*pos == n)
+ return;
+
+ get_bounds (curve1, t1l, t1r, &b1);
+ get_bounds (curve2, t2l, t2r, &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_get_point (curve1, t1[*pos], &c);
+
+ for (int i = 0; i < *pos; i++)
+ {
+ if (graphene_point_near (&c, &p[i], 0.1))
+ return;
+ }
+
+ p[*pos] = c;
+ (*pos)++;
+
+ return;
+ }
+
+ /* Note that in the conic case, we cannot just split the curves and
+ * pass the two halves down, since splitting changes the parametrization,
+ * and we need the t's to be valid parameters wrt to the original curve.
+ *
+ * So, instead, we determine the bounding boxes above by always starting
+ * from the original curve. That is a bit less efficient, but also works
+ * for conics.
+ */
+ general_intersect_recurse (curve1, curve2, t1l, t1l + d1, t2l, t2l + d2, t1, t2, p, n, pos);
+ general_intersect_recurse (curve1, curve2, t1l, t1l + d1, t2l + d2, t2r, t1, t2, p, n, pos);
+ general_intersect_recurse (curve1, curve2, t1l + d1, t1r, t2l, t2l + d2, t1, t2, p, n, pos);
+ general_intersect_recurse (curve1, curve2, t1l + d1, t1r, t2l + d2, t2r, t1, t2, p, n, pos);
+}
+
+static int
+general_intersect (const GskCurve *curve1,
+ const GskCurve *curve2,
+ float *t1,
+ float *t2,
+ graphene_point_t *p,
+ int n)
+{
+ int pos = 0;
+
+ general_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 if (c1 == &GSK_CURVE_CURVE_CLASS && c2 == &GSK_CURVE_CURVE_CLASS)
+ return curve_intersect (curve1, curve2, t1, t2, p, n);
+ else
+ return general_intersect (curve1, curve2, t1, t2, p, n);
+
+}
diff --git a/gsk/gskcurveprivate.h b/gsk/gskcurveprivate.h
index 571c7a7805..f64f3f4fae 100644
--- a/gsk/gskcurveprivate.h
+++ b/gsk/gskcurveprivate.h
@@ -113,6 +113,12 @@ void gsk_curve_get_bounds (const GskCurve
void gsk_curve_get_tight_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]
