From 64e73668a836f3dee54029af47c676dc387148d2 Mon Sep 17 00:00:00 2001
From: Matthias Clasen <mclasen@redhat.com>
Date: Sat, 12 Aug 2023 07:29:39 -0400
Subject: [PATCH] curve: Allow specifying epsilon for length->t

Passing FLT_EPSILON gets vert precise results
at a high cost.
---
 gsk/gskcontour.c                    |  2 +-
 gsk/gskcurve.c                      | 25 ++++++++++++++++++-------
 gsk/gskcurveprivate.h               |  3 ++-
 testsuite/gsk/curve-special-cases.c |  2 +-
 4 files changed, 22 insertions(+), 10 deletions(-)

diff --git a/gsk/gskcontour.c b/gsk/gskcontour.c
index ca0c2377df..632c99997d 100644
--- a/gsk/gskcontour.c
+++ b/gsk/gskcontour.c
@@ -761,7 +761,7 @@ gsk_standard_contour_get_point (const GskContour *contour,
       GskCurve curve;
 
       gsk_curve_init (&curve, self->ops[i]);
-      result->t = gsk_curve_at_length (&curve, distance);
+      result->t = gsk_curve_at_length (&curve, distance, 0.001);
 #else
       GskCurveLutEntry *e0, *e1;
 
diff --git a/gsk/gskcurve.c b/gsk/gskcurve.c
index 60d06465ef..1ecaa4718c 100644
--- a/gsk/gskcurve.c
+++ b/gsk/gskcurve.c
@@ -1553,6 +1553,10 @@ const double C[] = {
     0.0123412297999871995468056670700372915759,
 };
 
+/* Compute arclength by using Gauss quadrature on
+ *
+ * \int_0^z \sqrt{ (dx/dt)^2 + (dy/dt)^2 } dt
+ */
 float
 gsk_curve_get_length (const GskCurve *curve)
 {
@@ -1581,33 +1585,40 @@ gsk_curve_get_length (const GskCurve *curve)
   return z * sum;
 }
 
+/* Compute the inverse of the arclength using bisection,
+ * to a given precision
+ */
 float
 gsk_curve_at_length (const GskCurve *curve,
-                     float           length)
+                     float           length,
+                     float           epsilon)
 {
   float t1, t2, t, l;
   GskCurve c1;
+  //int loop_count = 0;
 
   t1 = 0;
   t2 = 1;
 
-  while (1)
+  while (t1 < t2)
     {
+      //loop_count++;
       t = (t1 + t2) / 2;
 
-      if (t1 == t2)
-        return t;
-
       gsk_curve_split (curve, t, &c1, NULL);
 
       l = gsk_curve_get_length (&c1);
-      if (fabs (length - l) < FLT_EPSILON)
-        return t;
+      if (fabs (length - l) < epsilon)
+        break;
       else if (l < length)
         t1 = t;
       else
         t2 = t;
     }
+
+  //g_print ("loop count %d\n", loop_count);
+
+  return t;
 }
 
 /* }}} */
diff --git a/gsk/gskcurveprivate.h b/gsk/gskcurveprivate.h
index 884d469f10..578a90422e 100644
--- a/gsk/gskcurveprivate.h
+++ b/gsk/gskcurveprivate.h
@@ -155,7 +155,8 @@ gboolean                gsk_curve_get_closest_point             (const GskCurve
 
 float                  gsk_curve_get_length                     (const GskCurve         *curve);
 float                  gsk_curve_at_length                      (const GskCurve         *curve,
-                                                                 float                   distance);
+                                                                 float                   distance,
+                                                                 float                   epsilon);
 
 G_END_DECLS
 
diff --git a/testsuite/gsk/curve-special-cases.c b/testsuite/gsk/curve-special-cases.c
index caf35a5af1..0af9ec9c1a 100644
--- a/testsuite/gsk/curve-special-cases.c
+++ b/testsuite/gsk/curve-special-cases.c
@@ -156,7 +156,7 @@ test_bad_split (void)
   l = gsk_curve_get_length (&c);
 
   t[0] = 0.5;
-  t[1] = gsk_curve_at_length (&c, 2);
+  t[1] = gsk_curve_at_length (&c, 2, FLT_EPSILON);
 
   for (unsigned int i = 0; i < G_N_ELEMENTS (t); i++)
     {