[gtk/wip/matthiasc/lottie-stroke: 16/16] Add code to split conics

[Matthias Clasen <matthiasc src gnome org>]

commit 378f1a531339c41fae940ba086c4026f5072b858
Author: Matthias Clasen <mclasen redhat com>
Date:   Wed Dec 2 00:15:33 2020 -0500

    Add code to split conics
    
    Using a traditional recursive approach.
    This code isn't used yet.

 gsk/gskpathstroke.c | 151 ++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 151 insertions(+)
---
diff --git a/gsk/gskpathstroke.c b/gsk/gskpathstroke.c
index 7f3f28581d..8bc892e2f1 100644
--- a/gsk/gskpathstroke.c
+++ b/gsk/gskpathstroke.c
@@ -142,6 +142,157 @@ split_bezier (const graphene_point_t *points,
   split_bezier_recurse (points, length, t, left, right, &lpos, &rpos);
 }
 
+static void
+get_rational_bezier (const graphene_point_t *p,
+                     float                  *w,
+                     int                     l,
+                     float                   t,
+                     graphene_point_t       *q)
+{
+  if (l == 1)
+    {
+      q->x = p[0].x;
+      q->y = p[0].y;
+    }
+  else
+    {
+      graphene_point_t *np;
+      float *nw;
+      int i;
+
+      np = g_alloca (sizeof (graphene_point3d_t) * (l - 1));
+      nw = g_alloca (sizeof (float) * (l - 1));
+
+      for (i = 0; i < l - 1; i++)
+        {
+          nw[i] = (1 - t) * w[i] + t * w[i + 1];
+          np[i].x = (1 - t) * (w[i]/nw[i]) * p[i].x + t * (w[i + 1]/nw[i]) * p[i + 1].x;
+          np[i].y = (1 - t) * (w[i]/nw[i]) * p[i].y + t * (w[i + 1]/nw[i]) * p[i + 1].y;
+        }
+      get_rational_bezier (np, nw, l - 1, t, q);
+    }
+}
+
+/* Given control points and weight for a rational quadratic Bezier
+ * and a t in the range [0,1], compute the point on the curve at t
+ */
+static void
+get_conic (const graphene_point_t points[3],
+           float                  weight,
+           float                  t,
+           graphene_point_t      *p)
+{
+  get_rational_bezier (points, (float [3]) { 1, weight, 1}, 3, t, p);
+}
+
+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);
+}
+
+/* Given control points and weight for a rational quadratic
+ * Bezier and t, create two sets of the same that give the
+ * same curve as the original and split the curve at t.
+ */
+static void
+split_conic (const graphene_point_t points[3],
+             float                  weight,
+             float                  t,
+             graphene_point_t       lp[3],
+             graphene_point_t       rp[3],
+             float                 *lw,
+             float                 *rw)
+{
+  graphene_point3d_t p[3];
+  graphene_point3d_t left[3], right[3];
+  int i;
+
+  /* do de Casteljau in homogeneous coordinates... */
+  for (i = 0; i < 3; i++)
+    {
+      p[i].x = points[i].x;
+      p[i].y = points[i].y;
+      p[i].z = 1;
+    }
+
+  p[1].x *= weight;
+  p[1].y *= weight;
+  p[1].z *= weight;
+
+  split_bezier3d (p, 3, t, left, right);
+
+  /* then project the control points down */
+  for (i = 0; i < 3; i++)
+    {
+      lp[i].x = left[i].x / left[i].z;
+      lp[i].y = left[i].y / left[i].z;
+      rp[i].x = right[i].x / right[i].z;
+      rp[i].y = right[i].y / right[i].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 (i = 0; i < 3; i++)
+    {
+      left[i].z /= left[0].z;
+      right[i].z /= right[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.
+   */
+  *lw = left[1].z / sqrt (left[2].z);
+  *rw = right[1].z / sqrt (right[0].z);
+}
+
 /* compute the angle between a, b and c in the range of [0, 360] */
 static float
 three_point_angle (const graphene_point_t *a,