update to latest 2geom (rev1497)

[inkscape.git] / src / 2geom / basic-intersection.cpp

diff --git a/src/2geom/basic-intersection.cpp b/src/2geom/basic-intersection.cpp
index 5375e5b58296c80009e83701c25bf009c9a5bb53..48b990445ed3c0296378237811ed247750aade2b 100644 (file)
--- a/src/2geom/basic-intersection.cpp
+++ b/src/2geom/basic-intersection.cpp
@@ -1,5 +1,8 @@
 #include <2geom/basic-intersection.h>
 #include <2geom/exception.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_multiroots.h>
+     
 
 unsigned intersect_steps = 0;
 
@@ -12,6 +15,7 @@ public:
     OldBezier() {
     }
     void split(double t, OldBezier &a, OldBezier &b) const;
+    Point operator()(double t) const;
     
     ~OldBezier() {}
 
@@ -152,7 +156,28 @@ void OldBezier::split(double t, OldBezier &left, OldBezier &right) const {
         right.p[j] = Vtemp[sz-1-j][j];
 }
 
-    
+/*
+ * split the curve at the midpoint, returning an array with the two parts
+ * Temporary storage is minimized by using part of the storage for the result
+ * to hold an intermediate value until it is no longer needed.
+ */
+Point OldBezier::operator()(double t) const {
+    const unsigned sz = p.size();
+    Geom::Point Vtemp[sz][sz];
+
+    /* Copy control points     */
+    std::copy(p.begin(), p.end(), Vtemp[0]);
+
+    /* Triangle computation    */
+    for (unsigned i = 1; i < sz; i++) {        
+        for (unsigned j = 0; j < sz - i; j++) {
+            Vtemp[i][j] = lerp(t, Vtemp[i-1][j], Vtemp[i-1][j+1]);
+        }
+    }
+    return Vtemp[sz-1][0];
+}
+
+
 /*
  * Test the bounding boxes of two OldBezier curves for interference.
  * Several observations:
@@ -324,7 +349,7 @@ double Lmax(Point p) {
 }
 
 unsigned wangs_theorem(OldBezier a) {
-    return 12; // seems a good approximation!
+    return 6; // seems a good approximation!
     double la1 = Lmax( ( a.p[2] - a.p[1] ) - (a.p[1] - a.p[0]) );
     double la2 = Lmax( ( a.p[3] - a.p[2] ) - (a.p[2] - a.p[1]) );
     double l0 = std::max(la1, la2);
@@ -337,6 +362,71 @@ unsigned wangs_theorem(OldBezier a) {
     return ra;
 }
 
+struct rparams
+{
+    OldBezier &A;
+    OldBezier &B;
+};
+     
+static int
+intersect_polish_f (const gsl_vector * x, void *params,
+              gsl_vector * f)
+{
+    const double x0 = gsl_vector_get (x, 0);
+    const double x1 = gsl_vector_get (x, 1);
+     
+    Geom::Point dx = ((struct rparams *) params)->A(x0) - 
+        ((struct rparams *) params)->B(x1);
+     
+    gsl_vector_set (f, 0, dx[0]);
+    gsl_vector_set (f, 1, dx[1]);
+     
+    return GSL_SUCCESS;
+}
+
+static void intersect_polish_root (OldBezier &A, double &s,
+                            OldBezier &B, double &t) {
+    const gsl_multiroot_fsolver_type *T;
+    gsl_multiroot_fsolver *sol;
+     
+    int status;
+    size_t iter = 0;
+     
+    const size_t n = 2;
+    struct rparams p = {A, B};
+    gsl_multiroot_function f = {&intersect_polish_f, n, &p};
+     
+    double x_init[2] = {s, t};
+    gsl_vector *x = gsl_vector_alloc (n);
+     
+    gsl_vector_set (x, 0, x_init[0]);
+    gsl_vector_set (x, 1, x_init[1]);
+     
+    T = gsl_multiroot_fsolver_hybrids;
+    sol = gsl_multiroot_fsolver_alloc (T, 2);
+    gsl_multiroot_fsolver_set (sol, &f, x);
+     
+    do
+    {
+        iter++;
+        status = gsl_multiroot_fsolver_iterate (sol);
+     
+        if (status)   /* check if solver is stuck */
+            break;
+     
+        status =
+            gsl_multiroot_test_residual (sol->f, 1e-12);
+    }
+    while (status == GSL_CONTINUE && iter < 1000);
+    
+    s = gsl_vector_get (sol->x, 0);
+    t = gsl_vector_get (sol->x, 1);
+    
+    gsl_multiroot_fsolver_free (sol);
+    gsl_vector_free (x);
+}
+
+
 std::vector<std::pair<double, double> > find_intersections( OldBezier a, OldBezier b)
 {
     std::vector<std::pair<double, double> > parameters;
@@ -346,6 +436,9 @@ std::vector<std::pair<double, double> > find_intersections( OldBezier a, OldBezi
                                b, 0., 1., wangs_theorem(b), 
                                parameters);
     }
+    for(unsigned i = 0; i < parameters.size(); i++)
+        intersect_polish_root(a, parameters[i].first,
+                              b, parameters[i].second);
     std::sort(parameters.begin(), parameters.end());
     return parameters;
 }