Commit LivePathEffect branch to trunk!

(disabled extension/internal/bitmap/*.* in build.xml to fix compilation)

(bzr r3472)

1 files changed, 268 insertions, 0 deletions

diff --git a/src/2geom/solve-bezier-one-d.cpp b/src/2geom/solve-bezier-one-d.cpp
new file mode 100644
index 000000000..558c10c0f
--- /dev/null
+++ b/src/2geom/solve-bezier-one-d.cpp
@@ -0,0 +1,268 @@
+#include "solver.h"
+#include "point.h"
+#include <algorithm>
+
+/*** Find the zeros of the bernstein function.  The code subdivides until it is happy with the
+ * linearity of the function.  This requires an O(degree^2) subdivision for each step, even when
+ * there is only one solution.
+ */
+
+namespace Geom{
+
+template<class t>
+static int SGN(t x) { return (x > 0 ? 1 : (x < 0 ? -1 : 0)); } 
+
+/*
+ *  Forward declarations
+ */
+static void 
+Bernstein(double const *V,
+          unsigned degree,
+          double t,
+          double *Left,
+          double *Right);
+
+static unsigned 
+control_poly_flat_enough(double const *V, unsigned degree,
+			 double left_t, double right_t);
+
+const unsigned MAXDEPTH = 64;	/*  Maximum depth for recursion */
+
+const double BEPSILON = ldexp(1.0,-MAXDEPTH-1); /*Flatness control value */
+
+/*
+ *  find_bernstein_roots : Given an equation in Bernstein-Bernstein form, find all 
+ *    of the roots in the open interval (0, 1).  Return the number of roots found.
+ */
+void
+find_bernstein_roots(double *w, /* The control points  */
+                     unsigned degree,	/* The degree of the polynomial */
+                     std::vector<double> &solutions, /* RETURN candidate t-values */
+                     unsigned depth,	/* The depth of the recursion */
+                     double left_t, double right_t)
+{  
+    unsigned 	n_crossings = 0;	/*  Number of zero-crossings */
+    
+    double split = 0.5;
+    int old_sign = SGN(w[0]);
+    for (unsigned i = 1; i <= degree; i++) {
+        int sign = SGN(w[i]);
+        if (sign) {
+            if (sign != old_sign && old_sign) {
+                split = (i-0.5)/degree;
+               n_crossings++;
+            }
+            old_sign = sign;
+        }
+    }
+    
+    switch (n_crossings) {
+    case 0: 	/* No solutions here	*/
+        return;
+	
+    case 1:
+ 	/* Unique solution	*/
+        /* Stop recursion when the tree is deep enough	*/
+        /* if deep enough, return 1 solution at midpoint  */
+        if (depth >= MAXDEPTH) {
+            solutions.push_back((left_t + right_t) / 2.0);
+            return;
+        }
+        
+        // I thought secant method would be faster here, but it'aint. -- njh
+
+        if (control_poly_flat_enough(w, degree, left_t, right_t)) {
+            const double Ax = right_t - left_t;
+            const double Ay = w[degree] - w[0];
+
+            solutions.push_back(left_t + Ax*w[0] / Ay);
+            return;
+        }
+        break;
+    }
+
+    /* Otherwise, solve recursively after subdividing control polygon  */
+    
+    
+    /* This bit is very clever (if I say so myself) - rather than arbitrarily subdividing at the t = 0.5 point, we subdivide at the most likely point of change of direction.  This buys us a factor of 10 speed up.  We also avoid lots of stack frames by avoiding tail recursion. */
+    double Left[degree+1],	/* New left and right  */
+           Right[degree+1];	/* control polygons  */
+    Bernstein(w, degree, split, Left, Right);
+    
+    double mid_t = left_t*(1-split) + right_t*split;
+    
+    find_bernstein_roots(Left,  degree, solutions, depth+1, left_t, mid_t);
+            
+    /* Solution is exactly on the subdivision point. */
+    if (Right[0] == 0)
+        solutions.push_back(mid_t);
+        
+    find_bernstein_roots(Right, degree, solutions, depth+1, mid_t, right_t);
+}
+
+/*
+ *  find_bernstein_roots : Given an equation in Bernstein-Bernstein form, find all 
+ *    of the roots in the interval [0, 1].  Return the number of roots found.
+ */
+void
+find_bernstein_roots_buggy(double *w, /* The control points  */
+                     unsigned degree,	/* The degree of the polynomial */
+                     std::vector<double> &solutions, /* RETURN candidate t-values */
+                     unsigned depth,	/* The depth of the recursion */
+                     double left_t, double right_t)
+{  
+    unsigned 	n_crossings = 0;	/*  Number of zero-crossings */
+    
+    int old_sign = SGN(w[0]);
+    for (unsigned i = 1; i <= degree; i++) {
+        int sign = SGN(w[i]);
+        if (sign) {
+            if (sign != old_sign && old_sign)
+               n_crossings++;
+            old_sign = sign;
+        }
+    }
+    
+    switch (n_crossings) {
+    case 0: 	/* No solutions here	*/
+        return;
+	
+    case 1:
+ 	/* Unique solution	*/
+        /* Stop recursion when the tree is deep enough	*/
+        /* if deep enough, return 1 solution at midpoint  */
+        if (depth >= MAXDEPTH) {
+            solutions.push_back((left_t + right_t) / 2.0);
+            return;
+        }
+        
+        // I thought secant method would be faster here, but it'aint. -- njh
+
+        if (control_poly_flat_enough(w, degree, left_t, right_t)) {
+            const double Ax = right_t - left_t;
+            const double Ay = w[degree] - w[0];
+
+            solutions.push_back(left_t + Ax*w[0] / Ay);
+            return;
+        }
+        break;
+    }
+
+    /* Otherwise, solve recursively after subdividing control polygon  */
+    
+    double split = 0.5;
+    
+    double Left[degree+1],	/* New left and right  */
+           Right[degree+1];	/* control polygons  */
+    Bernstein(w, degree, split, Left, Right);
+    
+    double mid_t = left_t*(1-split) + right_t*split;
+    
+    find_bernstein_roots_buggy(Left,  degree, solutions, depth+1, left_t, mid_t);
+            
+    /* Solution is exactly on the subdivision point. */
+    if (Right[0] == 0)
+        solutions.push_back(mid_t);
+        
+    find_bernstein_roots_buggy(Right, degree, solutions, depth+1, mid_t, right_t);
+}
+
+/*
+ *  control_poly_flat_enough :
+ *	Check if the control polygon of a Bernstein curve is flat enough
+ *	for recursive subdivision to bottom out.
+ *
+ */
+static unsigned 
+control_poly_flat_enough(double const *V, /* Control points	*/
+			 unsigned degree,
+			 double left_t, double right_t)	/* Degree of polynomial	*/
+{
+    /* Find the perpendicular distance from each interior control point to line connecting V[0] and
+     * V[degree] */
+
+    /* Derive the implicit equation for line connecting first */
+    /*  and last control points */
+    const double a = V[0] - V[degree];
+    const double b = right_t - left_t;
+    const double c = left_t * V[degree] - right_t * V[0] + a * left_t;
+
+    double max_distance_above = 0.0;
+    double max_distance_below = 0.0;
+    double ii = 0, dii = 1./degree;
+    for (unsigned i = 1; i < degree; i++) {
+        ii += dii;
+        /* Compute distance from each of the points to that line */
+        const double d = (a + V[i]) * ii*b  + c;
+        double dist = d*d;
+    // Find the largest distance
+        if (d < 0.0)
+            max_distance_below = std::min(max_distance_below, -dist);
+        else
+            max_distance_above = std::max(max_distance_above, dist);
+    }
+    
+    const double abSquared = (a * a) + (b * b);
+
+    const double intercept_1 = -(c + max_distance_above / abSquared);
+    const double intercept_2 = -(c + max_distance_below / abSquared);
+
+    /* Compute bounding interval*/
+    const double left_intercept = std::min(intercept_1, intercept_2);
+    const double right_intercept = std::max(intercept_1, intercept_2);
+
+    const double error = 0.5 * (right_intercept - left_intercept);
+    
+    if (error < BEPSILON * a)
+        return 1;
+    
+    return 0;
+}
+
+
+
+/*
+ *  Bernstein : 
+ *	Evaluate a Bernstein function at a particular parameter value
+ *      Fill in control points for resulting sub-curves.
+ * 
+ */
+static void 
+Bernstein(double const *V, /* Control pts	*/
+          unsigned degree,	/* Degree of bernstein curve */
+          double t,	/* Parameter value */
+          double *Left,	/* RETURN left half ctl pts */
+          double *Right)	/* RETURN right half ctl pts */
+{
+    double Vtemp[degree+1][degree+1];
+
+    /* Copy control points	*/
+    std::copy(V, V+degree+1, Vtemp[0]);
+
+    /* Triangle computation	*/
+    const double omt = (1-t);
+    Left[0] = Vtemp[0][0];
+    Right[degree] = Vtemp[0][degree];
+    for (unsigned i = 1; i <= degree; i++) {
+        for (unsigned j = 0; j <= degree - i; j++) {
+            Vtemp[i][j] = omt*Vtemp[i-1][j] + t*Vtemp[i-1][j+1];
+        }
+        Left[i] = Vtemp[i][0];
+        Right[degree-i] = Vtemp[i][degree-i];
+    }
+}
+
+};
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(substatement-open . 0))
+  indent-tabs-mode:nil
+  c-brace-offset:0
+  fill-column:99
+  End:
+  vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4 :
+*/
+