update 2geom (rev. 1569)

(bzr r6748)

1 files changed, 126 insertions, 51 deletions

diff --git a/src/2geom/basic-intersection.cpp b/src/2geom/basic-intersection.cpp
index 48b990445..334c23fea 100644
--- a/src/2geom/basic-intersection.cpp
+++ b/src/2geom/basic-intersection.cpp
@@ -1,8 +1,14 @@
+
+
+
 #include <2geom/basic-intersection.h>
+#include <2geom/sbasis-to-bezier.h>
 #include <2geom/exception.h>
+
+
 #include <gsl/gsl_vector.h>
 #include <gsl/gsl_multiroots.h>
-     
+
 
 unsigned intersect_steps = 0;
 
@@ -16,10 +22,10 @@ public:
     }
     void split(double t, OldBezier &a, OldBezier &b) const;
     Point operator()(double t) const;
-    
+
     ~OldBezier() {}
 
-    void bounds(double &minax, double &maxax, 
+    void bounds(double &minax, double &maxax,
                 double &minay, double &maxay) {
         // Compute bounding box for a
         minax = p[0][X];	 // These are the most likely to be extremal
@@ -45,17 +51,17 @@ public:
         }
 
     }
-    
+
 };
 
 static std::vector<std::pair<double, double> >
 find_intersections( OldBezier a, OldBezier b);
 
-static std::vector<std::pair<double, double> > 
+static std::vector<std::pair<double, double> >
 find_self_intersections(OldBezier const &Sb, D2<SBasis> const & A);
 
 std::vector<std::pair<double, double> >
-find_intersections( vector<Geom::Point> const & A, 
+find_intersections( vector<Geom::Point> const & A,
                     vector<Geom::Point> const & B) {
     OldBezier a, b;
     a.p = A;
@@ -63,23 +69,23 @@ find_intersections( vector<Geom::Point> const & A,
     return find_intersections(a,b);
 }
 
-std::vector<std::pair<double, double> > 
+std::vector<std::pair<double, double> >
 find_self_intersections(OldBezier const &/*Sb*/) {
     THROW_NOTIMPLEMENTED();
 }
 
-std::vector<std::pair<double, double> > 
+std::vector<std::pair<double, double> >
 find_self_intersections(D2<SBasis> const & A) {
     OldBezier Sb;
-    Sb.p = sbasis_to_bezier(A);
+    sbasis_to_bezier(Sb.p, A);
     return find_self_intersections(Sb, A);
 }
 
 
-static std::vector<std::pair<double, double> > 
+static std::vector<std::pair<double, double> >
 find_self_intersections(OldBezier const &Sb, D2<SBasis> const & A) {
 
-    
+
     vector<double> dr = roots(derivative(A[X]));
     {
         vector<double> dyr = roots(derivative(A[Y]));
@@ -90,9 +96,9 @@ find_self_intersections(OldBezier const &Sb, D2<SBasis> const & A) {
     // We want to be sure that we have no empty segments
     sort(dr.begin(), dr.end());
     unique(dr.begin(), dr.end());
-    
+
     std::vector<std::pair<double, double> > all_si;
-    
+
     vector<OldBezier> pieces;
     {
         OldBezier in = Sb, l, r;
@@ -104,7 +110,7 @@ find_self_intersections(OldBezier const &Sb, D2<SBasis> const & A) {
     }
     for(unsigned i = 0; i < dr.size()-1; i++) {
         for(unsigned j = i+1; j < dr.size()-1; j++) {
-            std::vector<std::pair<double, double> > section = 
+            std::vector<std::pair<double, double> > section =
                 find_intersections( pieces[i], pieces[j]);
             for(unsigned k = 0; k < section.size(); k++) {
                 double l = section[k].first;
@@ -118,17 +124,17 @@ find_self_intersections(OldBezier const &Sb, D2<SBasis> const & A) {
             }
         }
     }
-    
+
     // Because i is in order, all_si should be roughly already in order?
     //sort(all_si.begin(), all_si.end());
     //unique(all_si.begin(), all_si.end());
-    
+
     return all_si;
 }
 
 /* The value of 1.0 / (1L<<14) is enough for most applications */
 const double INV_EPS = (1L<<14);
-    
+
 /*
  * 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
@@ -142,12 +148,12 @@ void OldBezier::split(double t, OldBezier &left, OldBezier &right) const {
     std::copy(p.begin(), p.end(), Vtemp[0]);
 
     /* Triangle computation	*/
-    for (unsigned i = 1; i < sz; i++) {	
+    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]);
         }
     }
-    
+
     left.p.resize(sz);
     right.p.resize(sz);
     for (unsigned j = 0; j < sz; j++)
@@ -156,6 +162,7 @@ void OldBezier::split(double t, OldBezier &left, OldBezier &right) const {
         right.p[j] = Vtemp[sz-1-j][j];
 }
 
+#if 0
 /*
  * 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
@@ -169,13 +176,35 @@ Point OldBezier::operator()(double t) const {
     std::copy(p.begin(), p.end(), Vtemp[0]);
 
     /* Triangle computation	*/
-    for (unsigned i = 1; i < sz; i++) {	
+    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];
 }
+#endif
+
+// suggested by Sederberg.
+Point OldBezier::operator()(double t) const {
+    int n = p.size()-1;
+    double u, bc, tn, tmp;
+    int i;
+    Point r;
+    for(int dim = 0; dim < 2; dim++) {
+        u = 1.0 - t;
+        bc = 1;
+        tn = 1;
+        tmp = p[0][dim]*u;
+        for(i=1; i<n; i++){
+            tn = tn*t;
+            bc = bc*(n-i+1)/i;
+            tmp = (tmp + tn*bc*p[i][dim])*u;
+        }
+        r[dim] = (tmp + tn*t*p[n][dim]);
+    }
+    return r;
+}
 
 
 /*
@@ -183,11 +212,11 @@ Point OldBezier::operator()(double t) const {
  * Several observations:
  *	First, it is cheaper to compute the bounding box of the second curve
  *	and test its bounding box for interference than to use a more direct
- *	approach of comparing all control points of the second curve with 
+ *	approach of comparing all control points of the second curve with
  *	the various edges of the bounding box of the first curve to test
  * 	for interference.
  *	Second, after a few subdivisions it is highly probable that two corners
- *	of the bounding box of a given Bezier curve are the first and last 
+ *	of the bounding box of a given Bezier curve are the first and last
  *	control point.  Once this happens once, it happens for all subsequent
  *	subcurves.  It might be worth putting in a test and then short-circuit
  *	code for further subdivision levels.
@@ -209,39 +238,39 @@ bool intersect_BB( OldBezier a, OldBezier b ) {
     return not( ( minax > maxbx ) || ( minay > maxby )
                 || ( minbx > maxax ) || ( minby > maxay ) );
 }
-	
-/* 
- * Recursively intersect two curves keeping track of their real parameters 
+
+/*
+ * Recursively intersect two curves keeping track of their real parameters
  * and depths of intersection.
  * The results are returned in a 2-D array of doubles indicating the parameters
  * for which intersections are found.  The parameters are in the order the
  * intersections were found, which is probably not in sorted order.
- * When an intersection is found, the parameter value for each of the two 
+ * When an intersection is found, the parameter value for each of the two
  * is stored in the index elements array, and the index is incremented.
- * 
+ *
  * If either of the curves has subdivisions left before it is straight
  *	(depth > 0)
  * that curve (possibly both) is (are) subdivided at its (their) midpoint(s).
  * the depth(s) is (are) decremented, and the parameter value(s) corresponding
  * to the midpoints(s) is (are) computed.
- * Then each of the subcurves of one curve is intersected with each of the 
+ * Then each of the subcurves of one curve is intersected with each of the
  * subcurves of the other curve, first by testing the bounding boxes for
  * interference.  If there is any bounding box interference, the corresponding
  * subcurves are recursively intersected.
- * 
+ *
  * If neither curve has subdivisions left, the line segments from the first
- * to last control point of each segment are intersected.  (Actually the 
+ * to last control point of each segment are intersected.  (Actually the
  * only the parameter value corresponding to the intersection point is found).
  *
  * The apriori flatness test is probably more efficient than testing at each
  * level of recursion, although a test after three or four levels would
- * probably be worthwhile, since many curves become flat faster than their 
+ * probably be worthwhile, since many curves become flat faster than their
  * asymptotic rate for the first few levels of recursion.
  *
  * The bounding box test fails much more frequently than it succeeds, providing
  * substantial pruning of the search space.
  *
- * Each (sub)curve is subdivided only once, hence it is not possible that for 
+ * Each (sub)curve is subdivided only once, hence it is not possible that for
  * one final line intersection test the subdivision was at one level, while
  * for another final line intersection test the subdivision (of the same curve)
  * was at another.  Since the line segments share endpoints, the intersection
@@ -335,7 +364,7 @@ void recursively_intersect( OldBezier a, double t0, double t1, int deptha,
 }
 
 inline double log4( double x ) { return log(x)/log(4.); }
-    
+
 /*
  * Wang's theorem is used to estimate the level of subdivision required,
  * but only if the bounding boxes interfere at the top level.
@@ -354,7 +383,7 @@ unsigned wangs_theorem(OldBezier a) {
     double la2 = Lmax( ( a.p[3] - a.p[2] ) - (a.p[2] - a.p[1]) );
     double l0 = std::max(la1, la2);
     unsigned ra;
-    if( l0 * 0.75 * M_SQRT2 + 1.0 == 1.0 ) 
+    if( l0 * 0.75 * M_SQRT2 + 1.0 == 1.0 )
         ra = 0;
     else
         ra = (unsigned)ceil( log4( M_SQRT2 * 6.0 / 8.0 * INV_EPS * l0 ) );
@@ -367,63 +396,109 @@ 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) - 
+
+    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;
 }
 
+typedef union dbl_64{
+    long long i64;
+    double d64;
+};
+
+static double EpsilonBy(double value, int eps)
+{
+    dbl_64 s;
+    s.d64 = value;
+    s.i64 += eps;
+    return s.d64;
+}
+
+
 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);
+    
+    // This code does a neighbourhood search for minor improvements.
+    double best_v = L1(A(s) - B(t));
+    //std::cout  << "------\n" <<  best_v << std::endl;
+    Point best(s,t);
+    while (true) {
+        Point trial = best;
+        double trial_v = best_v;
+        for(int nsi = -1; nsi < 2; nsi++) {
+        for(int nti = -1; nti < 2; nti++) {
+            Point n(EpsilonBy(best[0], nsi),
+                    EpsilonBy(best[1], nti));
+            double c = L1(A(n[0]) - B(n[1]));
+            //std::cout << c << "; ";
+            if (c < trial_v) {
+                trial = n;
+                trial_v = c;
+            }
+        }
+        }
+        if(trial == best) {
+            //std::cout << "\n" << s << " -> " << s - best[0] << std::endl;
+            //std::cout << t << " -> " << t - best[1] << std::endl;
+            //std::cout << best_v << std::endl;
+            s = best[0];
+            t = best[1];
+            return;
+        } else {
+            best = trial;
+            best_v = trial_v;
+        }
+    }
 }
 
 
@@ -432,8 +507,8 @@ std::vector<std::pair<double, double> > find_intersections( OldBezier a, OldBezi
     std::vector<std::pair<double, double> > parameters;
     if( intersect_BB( a, b ) )
     {
-	recursively_intersect( a, 0., 1., wangs_theorem(a), 
-                               b, 0., 1., wangs_theorem(b), 
+	recursively_intersect( a, 0., 1., wangs_theorem(a),
+                               b, 0., 1., wangs_theorem(b),
                                parameters);
     }
     for(unsigned i = 0; i < parameters.size(); i++)