update to 2geom rev.1723

diff --git a/src/2geom/basic-intersection.cpp b/src/2geom/basic-intersection.cpp
index 2a40e7f453973d94c7e932893f58c843144cda55..16b5c0240f59019caa250848a1a6ea2df6816991 100644 (file)
--- a/src/2geom/basic-intersection.cpp
+++ b/src/2geom/basic-intersection.cpp
@@ -1,6 +1,3 @@
-
-
-
 #include <2geom/basic-intersection.h>
 #include <2geom/sbasis-to-bezier.h>
 #include <2geom/exception.h>
@@ -10,82 +7,57 @@
 #include <gsl/gsl_multiroots.h>
 
 
-unsigned intersect_steps = 0;
-
 using std::vector;
 namespace Geom {
 
-class OldBezier {
-public:
-    std::vector<Geom::Point> p;
-    OldBezier() {
-    }
-    void split(double t, OldBezier &a, OldBezier &b) const;
-    Point operator()(double t) const;
-
-    ~OldBezier() {}
-
-    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
-        maxax = p.back()[X];
-        if( minax > maxax )
-            std::swap(minax, maxax);
-        for(unsigned i = 1; i < p.size()-1; i++) {
-            if( p[i][X] < minax )
-                minax = p[i][X];
-            else if( p[i][X] > maxax )
-                maxax = p[i][X];
-        }
-
-        minay = p[0][Y];        // These are the most likely to be extremal
-        maxay = p.back()[Y];
-        if( minay > maxay )
-            std::swap(minay, maxay);
-        for(unsigned i = 1; i < p.size()-1; i++) {
-            if( p[i][Y] < minay )
-                minay = p[i][Y];
-            else if( p[i][Y] > maxay )
-                maxay = p[i][Y];
-        }
-
-    }
-
-};
-
-static std::vector<std::pair<double, double> >
-find_intersections( OldBezier a, OldBezier b);
+namespace detail{ namespace bezier_clipping {
+void portion (std::vector<Point> & B, Interval const& I);
+}; };
 
-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,
-                    vector<Geom::Point> const & B) {
-    OldBezier a, b;
-    a.p = A;
-    b.p = B;
-    return find_intersections(a,b);
+void find_intersections(std::vector<std::pair<double, double> > &xs,
+                        D2<SBasis> const & A,
+                        D2<SBasis> const & B) {
+    vector<Point> BezA, BezB;
+    sbasis_to_bezier(BezA, A);
+    sbasis_to_bezier(BezB, B);
+    
+    xs.clear();
+    
+    find_intersections(xs, BezA, BezB);
 }
 
-std::vector<std::pair<double, double> >
-find_self_intersections(OldBezier const &/*Sb*/) {
-    THROW_NOTIMPLEMENTED();
-}
+/*
+ * 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.
+ */
+void split(vector<Point> const &p, double t, 
+           vector<Point> &left, vector<Point> &right) {
+    const unsigned sz = p.size();
+    Geom::Point Vtemp[sz][sz];
 
-std::vector<std::pair<double, double> >
-find_self_intersections(D2<SBasis> const & A) {
-    OldBezier Sb;
-    sbasis_to_bezier(Sb.p, A);
-    return find_self_intersections(Sb, A);
-}
+    /* 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]);
+        }
+    }
 
-static std::vector<std::pair<double, double> >
-find_self_intersections(OldBezier const &Sb, D2<SBasis> const & A) {
+    left.resize(sz);
+    right.resize(sz);
+    for (unsigned j = 0; j < sz; j++)
+        left[j]  = Vtemp[j][0];
+    for (unsigned j = 0; j < sz; j++)
+        right[j] = Vtemp[sz-1-j][j];
+}
 
 
+void
+find_self_intersections(std::vector<std::pair<double, double> > &xs,
+                        D2<SBasis> const & A) {
     vector<double> dr = roots(derivative(A[X]));
     {
         vector<double> dyr = roots(derivative(A[Y]));
@@ -97,21 +69,21 @@ find_self_intersections(OldBezier const &Sb, D2<SBasis> const & A) {
     sort(dr.begin(), dr.end());
     unique(dr.begin(), dr.end());
 
-    std::vector<std::pair<double, double> > all_si;
-
-    vector<OldBezier> pieces;
+    vector<vector<Point> > pieces;
     {
-        OldBezier in = Sb, l, r;
+        vector<Point> in, l, r;
+        sbasis_to_bezier(in, A);
         for(unsigned i = 0; i < dr.size()-1; i++) {
-            in.split((dr[i+1]-dr[i]) / (1 - dr[i]), l, r);
+            split(in, (dr[i+1]-dr[i]) / (1 - dr[i]), l, r);
             pieces.push_back(l);
             in = r;
         }
     }
     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 =
-                find_intersections( pieces[i], pieces[j]);
+            std::vector<std::pair<double, double> > section;
+            
+            find_intersections( section, pieces[i], pieces[j]);
             for(unsigned k = 0; k < section.size(); k++) {
                 double l = section[k].first;
                 double r = section[k].second;
@@ -119,287 +91,31 @@ find_self_intersections(OldBezier const &Sb, D2<SBasis> const & A) {
                 if(j == i+1)
                     if((l == 1) && (r == 0))
                         continue;
-                all_si.push_back(std::make_pair((1-l)*dr[i] + l*dr[i+1],
+                xs.push_back(std::make_pair((1-l)*dr[i] + l*dr[i+1],
                                                 (1-r)*dr[j] + r*dr[j+1]));
             }
         }
     }
 
-    // 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
- * to hold an intermediate value until it is no longer needed.
- */
-void OldBezier::split(double t, OldBezier &left, OldBezier &right) 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]);
-        }
-    }
-
-    left.p.resize(sz);
-    right.p.resize(sz);
-    for (unsigned j = 0; j < sz; j++)
-        left.p[j]  = Vtemp[j][0];
-    for (unsigned j = 0; j < sz; j++)
-        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
- * 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];
-}
-#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;
+    // Because i is in order, xs should be roughly already in order?
+    //sort(xs.begin(), xs.end());
+    //unique(xs.begin(), xs.end());
 }
 
 
-/*
- * Test the bounding boxes of two OldBezier curves for interference.
- * 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
- *     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
- *     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.
- *     Third, in the final comparison (the interference test) the comparisons
- *     should both permit equality.  We want to find intersections even if they
- *     occur at the ends of segments.
- *     Finally, there are tighter bounding boxes that can be derived. It isn't
- *     clear whether the higher probability of rejection (and hence fewer
- *     subdivisions and tests) is worth the extra work.
- */
-
-bool intersect_BB( OldBezier a, OldBezier b ) {
-    double minax, maxax, minay, maxay;
-    a.bounds(minax, maxax, minay, maxay);
-    double minbx, maxbx, minby, maxby;
-    b.bounds(minbx, maxbx, minby, maxby);
-    // Test bounding box of b against bounding box of a
-    // Not >= : need boundary case
-    return not( ( minax > maxbx ) || ( minay > maxby )
-                || ( minbx > maxax ) || ( minby > maxay ) );
-}
-
-/*
- * 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
- * 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
- * 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
- * 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
- * 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
- * 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
- * is robust: a near-tangential intersection will yield zero or two
- * intersections.
- */
-void recursively_intersect( OldBezier a, double t0, double t1, int deptha,
-                          OldBezier b, double u0, double u1, int depthb,
-                          std::vector<std::pair<double, double> > &parameters)
-{
-    intersect_steps ++;
-    if( deptha > 0 )
-    {
-        OldBezier A[2];
-        a.split(0.5, A[0], A[1]);
-       double tmid = (t0+t1)*0.5;
-       deptha--;
-       if( depthb > 0 )
-        {
-           OldBezier B[2];
-            b.split(0.5, B[0], B[1]);
-           double umid = (u0+u1)*0.5;
-           depthb--;
-           if( intersect_BB( A[0], B[0] ) )
-               recursively_intersect( A[0], t0, tmid, deptha,
-                                     B[0], u0, umid, depthb,
-                                     parameters );
-           if( intersect_BB( A[1], B[0] ) )
-               recursively_intersect( A[1], tmid, t1, deptha,
-                                     B[0], u0, umid, depthb,
-                                     parameters );
-           if( intersect_BB( A[0], B[1] ) )
-               recursively_intersect( A[0], t0, tmid, deptha,
-                                     B[1], umid, u1, depthb,
-                                     parameters );
-           if( intersect_BB( A[1], B[1] ) )
-               recursively_intersect( A[1], tmid, t1, deptha,
-                                     B[1], umid, u1, depthb,
-                                     parameters );
-        }
-       else
-        {
-           if( intersect_BB( A[0], b ) )
-               recursively_intersect( A[0], t0, tmid, deptha,
-                                     b, u0, u1, depthb,
-                                     parameters );
-           if( intersect_BB( A[1], b ) )
-               recursively_intersect( A[1], tmid, t1, deptha,
-                                     b, u0, u1, depthb,
-                                     parameters );
-        }
-    }
-    else
-       if( depthb > 0 )
-        {
-           OldBezier B[2];
-            b.split(0.5, B[0], B[1]);
-           double umid = (u0 + u1)*0.5;
-           depthb--;
-           if( intersect_BB( a, B[0] ) )
-               recursively_intersect( a, t0, t1, deptha,
-                                     B[0], u0, umid, depthb,
-                                     parameters );
-           if( intersect_BB( a, B[1] ) )
-               recursively_intersect( a, t0, t1, deptha,
-                                     B[0], umid, u1, depthb,
-                                     parameters );
-        }
-       else // Both segments are fully subdivided; now do line segments
-        {
-           double xlk = a.p.back()[X] - a.p[0][X];
-           double ylk = a.p.back()[Y] - a.p[0][Y];
-           double xnm = b.p.back()[X] - b.p[0][X];
-           double ynm = b.p.back()[Y] - b.p[0][Y];
-           double xmk = b.p[0][X] - a.p[0][X];
-           double ymk = b.p[0][Y] - a.p[0][Y];
-           double det = xnm * ylk - ynm * xlk;
-           if( 1.0 + det == 1.0 )
-               return;
-           else
-            {
-               double detinv = 1.0 / det;
-               double s = ( xnm * ymk - ynm *xmk ) * detinv;
-               double t = ( xlk * ymk - ylk * xmk ) * detinv;
-               if( ( s < 0.0 ) || ( s > 1.0 ) || ( t < 0.0 ) || ( t > 1.0 ) )
-                   return;
-               parameters.push_back(std::pair<double, double>(t0 + s * ( t1 - t0 ),
-                                                         u0 + t * ( u1 - u0 )));
-            }
-        }
-}
-
-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.
- * Assuming there is a possible intersection, recursively_intersect is
- * used to find all the parameters corresponding to intersection points.
- * these are then sorted and returned in an array.
- */
-
-double Lmax(Point p) {
-    return std::max(fabs(p[X]), fabs(p[Y]));
-}
-
-unsigned wangs_theorem(OldBezier a) {
-    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);
-    unsigned ra;
-    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 ) );
-    std::cout << ra << std::endl;
-    return ra;
-}
+#include <gsl/gsl_multiroots.h>
+ 
+ 
 
 struct rparams
 {
-    OldBezier &A;
-    OldBezier &B;
+    D2<SBasis> const &A;
+    D2<SBasis> const &B;
 };
 
 static int
 intersect_polish_f (const gsl_vector * x, void *params,
-              gsl_vector * f)
+                    gsl_vector * f)
 {
     const double x0 = gsl_vector_get (x, 0);
     const double x1 = gsl_vector_get (x, 1);
@@ -413,7 +129,7 @@ intersect_polish_f (const gsl_vector * x, void *params,
     return GSL_SUCCESS;
 }
 
-typedef union dbl_64{
+union dbl_64{
     long long i64;
     double d64;
 };
@@ -427,8 +143,8 @@ static double EpsilonBy(double value, int eps)
 }
 
 
-static void intersect_polish_root (OldBezier &A, double &s,
-                            OldBezier &B, double &t) {
+static void intersect_polish_root (D2<SBasis> const &A, double &s,
+                                   D2<SBasis> const &B, double &t) {
     const gsl_multiroot_fsolver_type *T;
     gsl_multiroot_fsolver *sol;
 
@@ -502,23 +218,46 @@ static void intersect_polish_root (OldBezier &A, double &s,
 }
 
 
-std::vector<std::pair<double, double> > find_intersections( OldBezier a, OldBezier b)
+void polish_intersections(std::vector<std::pair<double, double> > &xs, 
+                        D2<SBasis> const  &A, D2<SBasis> const &B)
 {
-    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),
-                               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;
+    for(unsigned i = 0; i < xs.size(); i++)
+        intersect_polish_root(A, xs[i].first,
+                              B, xs[i].second);
 }
 
 
+ /**
+  * Compute the Hausdorf distance from A to B only.
+  */
+
+
+#if 0
+/** Compute the value of a bezier
+    Todo: find a good palce for this.
+ */
+// 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;
+}
+#endif
+
 /**
  * Compute the Hausdorf distance from A to B only.
  */

diff --git a/src/2geom/basic-intersection.h b/src/2geom/basic-intersection.h
index 77374b5ff934de57ba7955fce1fa2be02eadf22c..783df370ec006458f102aeabc61ce631f67cd3f5 100644 (file)
--- a/src/2geom/basic-intersection.h
+++ b/src/2geom/basic-intersection.h
@@ -42,24 +42,39 @@
 #include <vector>
 #include <utility>
 
+#define USE_RECURSIVE_INTERSECTOR 0
 
 namespace Geom {
 
-std::vector<std::pair<double, double> >
-find_intersections( D2<SBasis> const & A,
+#ifdef USE_RECURSIVE_INTERSECTOR
+void
+find_intersections( std::vector<std::pair<double, double> > &xs,
+                    D2<SBasis> const & A,
                     D2<SBasis> const & B);
 
-std::vector<std::pair<double, double> >
-find_self_intersections(D2<SBasis> const & A);
+void
+find_self_intersections(std::vector<std::pair<double, double> > &xs,
+                        D2<SBasis> const & A);
 
 // Bezier form
-std::vector<std::pair<double, double> >
-find_intersections( std::vector<Point> const & A,
-                    std::vector<Point> const & B);
+void
+find_intersections( std::vector<std::pair<double, double> > &xs,
+                    std::vector<Point> const & A,
+                    std::vector<Point> const & B,
+                    double precision = 1e-5);
 
-std::vector<std::pair<double, double> >
-find_self_intersections(std::vector<Point> const & A);
+// Bezier form
+void
+find_intersections_bezier_clipping( std::vector<std::pair<double, double> > &xs,
+                                    std::vector<Point> const & A,
+                                    std::vector<Point> const & B,
+                                    double precision = 1e-5);
+
+void
+find_self_intersections(std::vector<std::pair<double, double> > &xs,
+                        std::vector<Point> const & A);
 
+#else
 
 /*
  * find_intersection
@@ -72,11 +87,11 @@ find_self_intersections(std::vector<Point> const & A);
  *  This routine is based on the Bezier Clipping Algorithm,
  *  see: Sederberg, Nishita, 1990 - Curve intersection using Bezier clipping
  */
-void find_intersections (std::vector< std::pair<double, double> > & xs,
+void find_intersections_clipping (std::vector< std::pair<double, double> > & xs,
                          std::vector<Point> const& A,
                          std::vector<Point> const& B,
                          double precision = 1e-5);
-
+#endif
 /*
  * find_collinear_normal
  *
@@ -93,6 +108,23 @@ void find_collinear_normal (std::vector< std::pair<double, double> >& xs,
                             std::vector<Point> const& B,
                             double precision = 1e-5);
 
+void polish_intersections(std::vector<std::pair<double, double> > &xs, 
+                          D2<SBasis> const &A,
+                          D2<SBasis> const &B);
+
+void find_intersections(std::vector<std::pair<double, double> >& xs,
+                        D2<SBasis> const & A,
+                        D2<SBasis> const & B);
+
+void find_self_intersections(std::vector<std::pair<double, double> >& xs,
+                             D2<SBasis> const & A);
+
+
+void find_self_intersections(std::vector<std::pair<double, double> >& xs,
+                             std::vector<Point> const & A);
+
+
+
 /**
  * Compute the Hausdorf distance from A to B only.
  */

diff --git a/src/2geom/bezier-curve.h b/src/2geom/bezier-curve.h
index 589335d22629b352681d1ceded639d0a390b209c..acf878ee89fddb768b3ee45928ea5bcf0dbd6724 100644 (file)
--- a/src/2geom/bezier-curve.h
+++ b/src/2geom/bezier-curve.h
@@ -121,6 +121,10 @@ public:
   int winding(Point p) const {
     return SBasisCurve(toSBasis()).winding(p);
   }
+  
+  virtual int degreesOfFreedom() const {
+    return 2*order;
+  }
 
   std::vector<double>
   roots(double v, Dim2 d) const {

diff --git a/src/2geom/circle.h b/src/2geom/circle.h
index 1f9871276cfcab45235b2c4d66dc61f9fe40ec4a..27d4fcc3fe4942969e2d05c5f212d3c3bf5bb576 100644 (file)
--- a/src/2geom/circle.h
+++ b/src/2geom/circle.h
@@ -105,7 +105,6 @@ class Circle
   private:
     Point m_centre;
     Coord m_ray;
-    Coord m_angle;
 };
 
 

diff --git a/src/2geom/convex-cover.cpp b/src/2geom/convex-cover.cpp
index c4d8df3385bf77d6036550d6878fd081b91e33ec..e8ea2280d9d77ed808ef50e2e5f9a5efe8346d19 100644 (file)
--- a/src/2geom/convex-cover.cpp
+++ b/src/2geom/convex-cover.cpp
@@ -179,6 +179,13 @@ int mod(int i, int l) {
  * Tests if a point is left (outside) of a particular segment, n. */
 bool
 ConvexHull::is_left(Point p, int n) {
+    return SignedTriangleArea((*this)[n], (*this)[n+1], p) >= 0;
+}
+
+/*** ConvexHull::strict_left
+ * Tests if a point is left (outside) of a particular segment, n. */
+bool
+ConvexHull::is_strict_left(Point p, int n) {
     return SignedTriangleArea((*this)[n], (*this)[n+1], p) > 0;
 }
 
@@ -193,6 +200,20 @@ ConvexHull::find_left(Point p) {
     }
     return -1;
 }
+
+
+/*** ConvexHull::find_positive
+ * May return any number n where the segment n -> n + 1 (possibly looped around) in the hull such
+ * that the point is on the wrong side to be within the hull.  Returns -1 if it is within the hull.*/
+int
+ConvexHull::find_strict_left(Point p) {
+    int l = boundary.size(); //Who knows if C++ is smart enough to optimize this?
+    for(int i = 0; i < l; i++) {
+        if(is_strict_left(p, i)) return i;
+    }
+    return -1;
+}
+
 //OPT: do a spread iteration - quasi-random with no repeats and full coverage.
 
 /*** ConvexHull::contains_point
@@ -203,6 +224,14 @@ ConvexHull::contains_point(Point p) {
     return find_left(p) == -1;
 }
 
+/*** ConvexHull::strict_contains_point
+ * In order to test whether a point is strictly inside (not on the boundary) a convex hull we can travel once around the outside making
+ * sure that each triangle made from an edge and the point has positive area. */
+bool
+ConvexHull::strict_contains_point(Point p) {
+    return find_strict_left(p) == -1;
+}
+
 /*** ConvexHull::add_point
  * to add a point we need to find whether the new point extends the boundary, and if so, what it
  * obscures.  Tarjan?  Jarvis?*/
@@ -219,9 +248,9 @@ ConvexHull::merge(Point p) {
 
     bool pushed = false;
 
-    bool pre = is_left(p, -1);
+    bool pre = is_strict_left(p, -1);
     for(int i = 0; i < l; i++) {
-        bool cur = is_left(p, i);
+        bool cur = is_strict_left(p, i);
         if(pre) {
             if(cur) {
                 if(!pushed) {

diff --git a/src/2geom/convex-cover.h b/src/2geom/convex-cover.h
index f221777463b511a028c10c7d321fd4ca9c6fcfdb..52410896500b6df8de161b9dea8012206835678a 100644 (file)
--- a/src/2geom/convex-cover.h
+++ b/src/2geom/convex-cover.h
@@ -65,6 +65,7 @@ public:
 
     void merge(Point p);
     bool contains_point(Point p);
+    bool strict_contains_point(Point p);
 
     inline Point operator[](int i) const {
 
@@ -122,7 +123,9 @@ public:
     Point const * furthest(Point direction) const;
 
     bool is_left(Point p, int n);
+    bool is_strict_left(Point p, int n);
     int find_left(Point p);
+    int find_strict_left(Point p);
     double narrowest_diameter(Point &a, Point &b, Point &c);
 
 };

diff --git a/src/2geom/crossing.h b/src/2geom/crossing.h
index 3ad034ac9b435a846a405e1906462a55e0fc6d43..7737b24a9685a8b7d3d128586cbdedbf9c449bdf 100644 (file)
--- a/src/2geom/crossing.h
+++ b/src/2geom/crossing.h
@@ -3,9 +3,10 @@
  * \brief  \todo brief description
  *
  * Authors:
- *      ? <?@?.?>
+ *      Michael Sloane <?@?.?>
+ *      Marco <?@?.?>
  * 
- * Copyright ?-?  authors
+ * Copyright 2006-2008  authors
  *
  * This library is free software; you can redistribute it and/or
  * modify it either under the terms of the GNU Lesser General Public
@@ -38,6 +39,7 @@
 #include <vector>
 #include <2geom/rect.h>
 #include <2geom/sweep.h>
+#include <boost/optional/optional.hpp>
 
 namespace Geom {
 
@@ -58,6 +60,8 @@ struct Crossing {
     bool onIx(unsigned ix) const { return a == ix || b == ix; }
 };
 
+typedef boost::optional<Crossing> OptCrossing;
+
 
 /*
 struct Edge {

diff --git a/src/2geom/curve.h b/src/2geom/curve.h
index 08cc2380ceb552c3a8cdb343ebc97988bf2acd66..ce1fec3a65a20d1a409ef38345e144cccb6ae0a7 100644 (file)
--- a/src/2geom/curve.h
+++ b/src/2geom/curve.h
@@ -83,6 +83,8 @@ public:
 
   virtual int winding(Point p) const { return root_winding(*this, p); }
 
+  virtual int degreesOfFreedom() const { return 0;}
+
   //mental: review these
   virtual Curve *portion(double f, double t) const = 0;
   virtual Curve *reverse() const { return portion(1, 0); }
@@ -133,8 +135,8 @@ public:
    */
   virtual Point unitTangentAt(Coord t, unsigned n = 3) const
   {
-       std::vector<Point> derivs = pointAndDerivatives(t, n);
-       for (unsigned deriv_n = 1; deriv_n < derivs.size(); deriv_n++) {
+    std::vector<Point> derivs = pointAndDerivatives(t, n);
+    for (unsigned deriv_n = 1; deriv_n < derivs.size(); deriv_n++) {
       Coord length = derivs[deriv_n].length();
       if ( ! are_near(length, 0) ) {
          // length of derivative is non-zero, so return unit vector

diff --git a/src/2geom/d2-sbasis.cpp b/src/2geom/d2-sbasis.cpp
index 55c7ef55e0d0138ecd3acecc972d11bc1d76c610..2c52e4782908d282dc16585f1267ca92ba78c388 100644 (file)
--- a/src/2geom/d2-sbasis.cpp
+++ b/src/2geom/d2-sbasis.cpp
@@ -113,12 +113,12 @@ Piecewise<D2<SBasis> > force_continuity(Piecewise<D2<SBasis> > const &f, double
                 SBasis &cur_sb =result.segs[cur][dim];
                 Coord const c=pt0[dim];
                 if (prev_sb.empty()) {
-                  prev_sb.push_back(Linear(0.0, c));
+                  prev_sb = SBasis(Linear(0.0, c));
                 } else {
                   prev_sb[0][1] = c;
                 }
                 if (cur_sb.empty()) {
-                  cur_sb.push_back(Linear(c, 0.0));
+                  cur_sb = SBasis(Linear(c, 0.0));
                 } else {
                   cur_sb[0][0] = c;
                 }
@@ -157,7 +157,7 @@ static void set_first_point(Piecewise<D2<SBasis> > &f, Point a){
     }
     for (unsigned dim=0; dim<2; dim++){
         if (f.segs.front()[dim].size() == 0){
-            f.segs.front()[dim].push_back(Linear(a[dim],0));
+            f.segs.front()[dim] = SBasis(Linear(a[dim],0));
         }else{
             f.segs.front()[dim][0][0] = a[dim];
         }
@@ -170,7 +170,7 @@ static void set_last_point(Piecewise<D2<SBasis> > &f, Point a){
     }
     for (unsigned dim=0; dim<2; dim++){
         if (f.segs.back()[dim].size() == 0){
-            f.segs.back()[dim].push_back(Linear(0,a[dim]));
+            f.segs.back()[dim] = SBasis(Linear(0,a[dim]));
         }else{
             f.segs.back()[dim][0][1] = a[dim];
         }

diff --git a/src/2geom/d2.h b/src/2geom/d2.h
index a14e3b0ebe027b4aa28ba12720b3a41d8f7e7744..f55f6a5969fe770de62134c9f8e904b66296f0e6 100644 (file)
--- a/src/2geom/d2.h
+++ b/src/2geom/d2.h
@@ -124,6 +124,12 @@ inline D2<T> portion(const D2<T> &a, Coord f, Coord t) {
     return D2<T>(portion(a[X], f, t), portion(a[Y], f, t));
 }
 
+template <typename T>
+inline D2<T> portion(const D2<T> &a, Interval i) {
+    boost::function_requires<FragmentConcept<T> >();
+    return D2<T>(portion(a[X], i), portion(a[Y], i));
+}
+
 //IMPL: boost::EqualityComparableConcept
 template <typename T>
 inline bool

diff --git a/src/2geom/ellipse.cpp b/src/2geom/ellipse.cpp
index 2690a091bfafab64276f10328a0e12981c82604c..10071d09a15f35b8ac984e5a86518606dbbb9e15 100644 (file)
--- a/src/2geom/ellipse.cpp
+++ b/src/2geom/ellipse.cpp
@@ -264,6 +264,12 @@ Ellipse Ellipse::transformed(Matrix const& m) const
     return e;
 }
 
+Ellipse::Ellipse(Geom::Circle const &c)
+{
+    m_centre = c.center();
+    m_ray[X] = m_ray[Y] = c.ray();
+}
+
 }  // end namespace Geom
 
 

diff --git a/src/2geom/ellipse.h b/src/2geom/ellipse.h
index d5b882bc4e952fb4e9823386505ad69dbc4c43f9..7ed04e51b396ff4fcd0bae92cebd8b653019fb5a 100644 (file)
--- a/src/2geom/ellipse.h
+++ b/src/2geom/ellipse.h
@@ -44,6 +44,7 @@ namespace Geom
 {
 
 class SVGEllipticalArc;
+class Circle;
 
 class Ellipse
 {
@@ -67,6 +68,8 @@ class Ellipse
     {
         set(points);
     }
+    
+    Ellipse(Geom::Circle const &c);
 
     void set(double cx, double cy, double rx, double ry, double a)
     {

diff --git a/src/2geom/elliptical-arc.h b/src/2geom/elliptical-arc.h
index 25c79a363e3900ab877efd58b1fe430d931e12e2..b0c0bd9df2b906e36e1730805ff0ac42abae4d1d 100644 (file)
--- a/src/2geom/elliptical-arc.h
+++ b/src/2geom/elliptical-arc.h
@@ -200,6 +200,8 @@ class EllipticalArc : public Curve
     {
        return SBasisCurve(toSBasis()).winding(p);
     }
+
+    int degreesOfFreedom() const { return 7;}
     
     Curve *derivative() const;
     

diff --git a/src/2geom/geom.cpp b/src/2geom/geom.cpp
index f9b1a664bcc0c451d06e017c7665019c8fc16684..5eade57f2d4bb73467ea1934e452a2bd0445d81c 100644 (file)
--- a/src/2geom/geom.cpp
+++ b/src/2geom/geom.cpp
@@ -1,5 +1,4 @@
 /**
- *  \file geom.cpp
  *  \brief Various geometrical calculations.
  */
 
@@ -85,7 +84,6 @@ line_intersection(Geom::Point const &n0, double const d0,
 
 
 
-
 /* ccw exists as a building block */
 int
 intersector_ccw(const Geom::Point& p0, const Geom::Point& p1,
@@ -314,13 +312,35 @@ rect_line_intersect(Geom::Point const &c0, Geom::Point const &c1,
     return results;
 }
 
-std::vector<Geom::Point>
+/** Determine whether \& where an (infinite) line intersects a rectangle.
+ *
+ * \a c0, \a c1 are diagonal corners of the rectangle and
+ * \a p1, \a p1 are distinct points on the line
+ *
+ * \return A list (possibly empty) of points of intersection.  If two such points (say \a r0 and \a
+ * r1) then it is guaranteed that the order of \a r0, \a r1 along the line is the same as the that
+ * of \a c0, \a c1 (i.e., the vectors \a r1 - \a r0 and \a p1 - \a p0 point into the same
+ * direction).
+ */
+boost::optional<LineSegment>
 rect_line_intersect(Geom::Rect &r,
-                    Geom::Point const &p0, Geom::Point const &p1)
+                    Geom::LineSegment ls)
 {
-    return rect_line_intersect(r.min(), r.max(), p0, p1);
+    std::vector<Point> results;
+    
+    results =  rect_line_intersect(r.min(), r.max(), ls[0], ls[1]);
+    if(results.size() == 2) {
+        return LineSegment(results[0], results[1]);
+    } 
+    return boost::optional<LineSegment>();
 }
 
+boost::optional<LineSegment>
+rect_line_intersect(Geom::Rect &r,
+                    Geom::Line l)
+{
+    return rect_line_intersect(r, l.segment(0, 1));
+}
 
 /**
  * polyCentroid: Calculates the centroid (xCentroid, yCentroid) and area of a polygon, given its
@@ -336,7 +356,7 @@ rect_line_intersect(Geom::Rect &r,
     * for now we require the path to be a polyline and assume it is closed.
 **/
 
-int centroid(std::vector<Geom::Point> p, Geom::Point& centroid, double &area) {
+int centroid(std::vector<Geom::Point> const &p, Geom::Point& centroid, double &area) {
     const unsigned n = p.size();
     if (n < 3)
         return 1;

diff --git a/src/2geom/geom.h b/src/2geom/geom.h
index d0af7d7d24a5971d95c0bd5d3fc4dd635fe0e3f0..9233696d74fc631d07380c55c56eb18e1a2b2242 100644 (file)
--- a/src/2geom/geom.h
+++ b/src/2geom/geom.h
@@ -37,8 +37,10 @@
 //TODO: move somewhere else
 
 #include <vector>
-#include <2geom/point.h>
-#include <2geom/rect.h>
+#include <2geom/forward.h>
+#include <boost/optional/optional.hpp>
+#include <2geom/bezier-curve.h>
+#include <2geom/line.h>
 
 namespace Geom {
 
@@ -55,6 +57,9 @@ intersector_ccw(const Geom::Point& p0, const Geom::Point& p1,
 
 /* intersectors */
 
+#if 0
+// Use the new routines provided in line.h
+
 IntersectorKind
 line_intersection(Geom::Point const &n0, double const d0,
                  Geom::Point const &n1, double const d1,
@@ -69,17 +74,23 @@ IntersectorKind
 line_twopoint_intersect(Geom::Point const &p00, Geom::Point const &p01,
                        Geom::Point const &p10, Geom::Point const &p11,
                        Geom::Point &result);
+#endif
 
+#if 0
 std::vector<Geom::Point>
 rect_line_intersect(Geom::Point const &E, Geom::Point const &F,
                     Geom::Point const &p0, Geom::Point const &p1);
+#endif
 
+boost::optional<Geom::LineSegment>
+rect_line_intersect(Geom::Rect &r,
+                    Geom::LineSegment ls);
 
-std::vector<Geom::Point>
+boost::optional<Geom::LineSegment>
 rect_line_intersect(Geom::Rect &r,
-                    Geom::Point const &p0, Geom::Point const &p1);
+                    Geom::Line l);
 
-int centroid(std::vector<Geom::Point> p, Geom::Point& centroid, double &area);
+int centroid(std::vector<Geom::Point> const &p, Geom::Point& centroid, double &area);
 
 }
 

diff --git a/src/2geom/hvlinesegment.h b/src/2geom/hvlinesegment.h
index ccb95afdbc0b90f33dabf6c0e18ee4127a42c491..ac91ec80a73909b6511d9d8a8fd29870fb882e6e 100644 (file)
--- a/src/2geom/hvlinesegment.h
+++ b/src/2geom/hvlinesegment.h
@@ -136,6 +136,8 @@ class HLineSegment : public Curve
         return m_line_seg.winding(p);
     }
 
+    int degreesOfFreedom() const { return 3;}
+
     std::vector<double>
     roots(double v, Dim2 d) const
     {
@@ -375,6 +377,8 @@ class VLineSegment : public Curve
         return m_line_seg.winding(p);
     }
 
+    int degreesOfFreedom() const { return 3;}
+    
     std::vector<double>
     roots(double v, Dim2 d) const
     {

diff --git a/src/2geom/interval.h b/src/2geom/interval.h
index 00952858603b72f6f0d2f27ac2ca3fdbeb681fd0..d4dae41b429951010b19e49fab988d7943e7ed60 100644 (file)
--- a/src/2geom/interval.h
+++ b/src/2geom/interval.h
@@ -258,10 +258,18 @@ inline OptInterval intersect(const Interval & a, const Interval & b) {
     Coord u = std::max(a.min(), b.min()),
           v = std::min(a.max(), b.max());
     //technically >= might be incorrect, but singulars suck
-    return u >= v ? OptInterval()
+    return u > v ? OptInterval()
                   : OptInterval(Interval(u, v));
 }
 
+#ifdef _GLIBCXX_IOSTREAM
+inline std::ostream &operator<< (std::ostream &os, 
+                                 const Geom::Interval &I) {
+    os << "Interval("<<I[0] << ", "<<I[1] << ")";
+    return os;
+}
+#endif
+
 }
 #endif //SEEN_INTERVAL_H
 

diff --git a/src/2geom/linear.h b/src/2geom/linear.h
index e8828b283f3df5064d6f4bcaf4a474246ff483ef..a7f4c8f21950364b3bd00a8c7f73ab69f056da7d 100644 (file)
--- a/src/2geom/linear.h
+++ b/src/2geom/linear.h
@@ -52,36 +52,12 @@ inline double lerp(double t, double a, double b) { return a*(1-t) + b*t; }
 
 class SBasis;
 
-class Hat{
-public:
-    Hat () {}
-    Hat(double d) :d(d) {}
-    operator double() const { return d; }
-    double d;
-};
-
-class Tri{
-public:
-    Tri () {}
-    Tri(double d) :d(d) {}
-    operator double() const { return d; }
-    double d;
-};
-
 class Linear{
 public:
     double a[2];
     Linear() {}
     Linear(double aa, double b) {a[0] = aa; a[1] = b;}
-    Linear(Hat h, Tri t) {
-        a[0] = double(h) - double(t)/2; 
-        a[1] = double(h) + double(t)/2;
-    }
-
-    Linear(Hat h) {
-        a[0] = double(h); 
-        a[1] = double(h);
-    }
+    Linear(double aa) {a[0] = aa; a[1] = aa;}
 
     double operator[](const int i) const {
         assert(i >= 0);
@@ -113,10 +89,10 @@ public:
     inline OptInterval bounds_fast() const { return bounds_exact(); }
     inline OptInterval bounds_local(double u, double v) const { return Interval(valueAt(u), valueAt(v)); }
 
-    operator Tri() const {
+    double tri() const {
         return a[1] - a[0];
     }
-    operator Hat() const {
+    double hat() const {
         return (a[1] + a[0])/2;
     }
 };

diff --git a/src/2geom/path-intersection.cpp b/src/2geom/path-intersection.cpp
index 5a4e7602056e1234948c82690f558d1d1e17a874..bfb1fb96c508f1820fbc3742aa87f5608598a0bf 100644 (file)
--- a/src/2geom/path-intersection.cpp
+++ b/src/2geom/path-intersection.cpp
@@ -78,7 +78,7 @@ int winding(Path const &path, Point p) {
                     const double fudge = 0.01;
                     if(cmp(y, next->valueAt(fudge, Y)) == initial_to_ray) {
                         wind += int(c);
-                        std::cout << "!!!!!" << int(c) << " ";
+                        //std::cout << "!!!!!" << int(c) << " ";
                     }
                     iter = next; // No increment, as the rest of the thing hasn't been counted.
                 } else {
@@ -86,7 +86,7 @@ int winding(Path const &path, Point p) {
                     if(cmp(y, ny) == initial_to_ray) {
                         //Is a continuation through the ray, so counts windingwise
                         wind += int(c);
-                        std::cout << "!!!!!" << int(c) << " ";
+                        //std::cout << "!!!!!" << int(c) << " ";
                     }
                     iter = ++next;
                 }
@@ -177,6 +177,7 @@ linear_intersect(Point A0, Point A1, Point B0, Point B1,
 }
 
 
+#if 0
 typedef union dbl_64{
     long long i64;
     double d64;
@@ -196,6 +197,7 @@ static double EpsilonOf(double value)
     else
         return value - s.d64;
 }
+#endif
 
 struct rparams {
     Curve const &A;

diff --git a/src/2geom/pathvector.h b/src/2geom/pathvector.h
index 9efae7c7319ada1bcade2c2b6428ca75d27bfb12..d1d785a07a1897dd06e8226772acaf054b3d0864 100644 (file)
--- a/src/2geom/pathvector.h
+++ b/src/2geom/pathvector.h
@@ -38,7 +38,6 @@
 
 #include <2geom/forward.h>
 #include <2geom/path.h>
-#include <2geom/rect.h>
 #include <2geom/transforms.h>
 
 namespace Geom {

diff --git a/src/2geom/piecewise.h b/src/2geom/piecewise.h
index bbd1f054aeb423f56f9550b65ab1179c4a91c1fa..31bf6872af447aa4d321b2c7511142ea7da6e9ca 100644 (file)
--- a/src/2geom/piecewise.h
+++ b/src/2geom/piecewise.h
@@ -754,8 +754,21 @@ Piecewise<T> reverse(Piecewise<T> const &f) {
 }
 
 
+/**
+ *  Interpolates between a and b.
+ *  \return a if t = 0, b if t = 1, or an interpolation between a and b for t in [0,1]
+ */
+template<typename T>
+Piecewise<T> lerp(Piecewise<T> const &a, Piecewise<T> b, double t) {
+    // Make sure both paths have the same number of segments and cuts at the same locations
+    b.setDomain(a.domain());
+    Piecewise<T> pA = partition(a, b.cuts);
+    Piecewise<T> pB = partition(b, a.cuts);
+
+    return (pA*(1-t)  +  pB*t);
 }
 
+}
 #endif //SEEN_GEOM_PW_SB_H
 /*
   Local Variables:

diff --git a/src/2geom/sbasis-2d.cpp b/src/2geom/sbasis-2d.cpp
index e3a407094c769a55c0afa223f3fa8c649ca83838..399fb8595faf0888ccef69e9135f7a9047973409 100644 (file)
--- a/src/2geom/sbasis-2d.cpp
+++ b/src/2geom/sbasis-2d.cpp
@@ -4,7 +4,7 @@
 namespace Geom{
 
 SBasis extract_u(SBasis2d const &a, double u) {
-    SBasis sb;
+    SBasis sb(a.vs, Linear());
     double s = u*(1-u);
     
     for(unsigned vi = 0; vi < a.vs; vi++) {
@@ -14,14 +14,14 @@ SBasis extract_u(SBasis2d const &a, double u) {
             bo += (extract_u(a.index(ui, vi), u))*sk;
             sk *= s;
         }
-        sb.push_back(bo);
+        sb[vi] = bo;
     }
     
     return sb;
 }
 
 SBasis extract_v(SBasis2d const &a, double v) {
-    SBasis sb;
+    SBasis sb(a.us, Linear());
     double s = v*(1-v);
     
     for(unsigned ui = 0; ui < a.us; ui++) {
@@ -31,7 +31,7 @@ SBasis extract_v(SBasis2d const &a, double v) {
             bo += (extract_v(a.index(ui, vi), v))*sk;
             sk *= s;
         }
-        sb.push_back(bo);
+        sb[ui] = bo;
     }
     
     return sb;
@@ -105,7 +105,6 @@ SBasis2d partial_derivative(SBasis2d const &f, int dim) {
  */
 D2<SBasis>
 sb2dsolve(SBasis2d const &f, Geom::Point const &A, Geom::Point const &B, unsigned degmax){
-    D2<SBasis>result(Linear(A[X],B[X]),Linear(A[Y],B[Y]));
     //g_warning("check f(A)= %f = f(B) = %f =0!", f.apply(A[X],A[Y]), f.apply(B[X],B[Y]));
 
     SBasis2d dfdu = partial_derivative(f, 0);
@@ -115,8 +114,11 @@ sb2dsolve(SBasis2d const &f, Geom::Point const &A, Geom::Point const &B, unsigne
     Geom::Point nA = dfA/(dfA[X]*dfA[X]+dfA[Y]*dfA[Y]);
     Geom::Point nB = dfB/(dfB[X]*dfB[X]+dfB[Y]*dfB[Y]);
 
+    D2<SBasis>result(SBasis(degmax, Linear()), SBasis(degmax, Linear()));
     double fact_k=1;
     double sign = 1.;
+    for(int dim = 0; dim < 2; dim++)
+        result[dim][0] = Linear(A[dim],B[dim]);
     for(unsigned k=1; k<degmax; k++){
         // these two lines make the solutions worse!
         //fact_k *= k;
@@ -128,8 +130,8 @@ sb2dsolve(SBasis2d const &f, Geom::Point const &A, Geom::Point const &B, unsigne
         double bx = -sign*reste[1]/fact_k*nB[X];
         double by = -sign*reste[1]/fact_k*nB[Y];
 
-        result[X].push_back(Linear(ax,bx));
-        result[Y].push_back(Linear(ay,by));
+        result[X][k] = Linear(ax,bx);
+        result[Y][k] = Linear(ay,by);
         //sign *= 3;
     }    
     return result;

diff --git a/src/2geom/sbasis-curve.h b/src/2geom/sbasis-curve.h
index b45e63eb83af318200aa342a88fc50889af0cfd8..893cd23aff680bbd2ceadb8ccb92114d0a17a7e6 100644 (file)
--- a/src/2geom/sbasis-curve.h
+++ b/src/2geom/sbasis-curve.h
@@ -99,6 +99,8 @@ public:
 
   D2<SBasis> toSBasis() const { return inner; }
 
+  virtual int degreesOfFreedom() const { return inner[0].degreesOfFreedom() + inner[1].degreesOfFreedom();
+  }
 };
 
 

diff --git a/src/2geom/sbasis-geometric.cpp b/src/2geom/sbasis-geometric.cpp
index 96a5ed0ce0dddb8c061c3f68c081e810f9c36a28..d27255749946703fcf8c616b0958b673ca0d0b0a 100644 (file)
--- a/src/2geom/sbasis-geometric.cpp
+++ b/src/2geom/sbasis-geometric.cpp
@@ -27,6 +27,11 @@ using namespace std;
 
 //Some utils first.
 //TODO: remove this!! 
+/** 
+ * Return a list of doubles that appear in both a and b to within error tol
+ * a, b, vector of double
+ * tol tolerance
+ */
 static vector<double> 
 vect_intersect(vector<double> const &a, vector<double> const &b, double tol=0.){
     vector<double> inter;
@@ -72,21 +77,7 @@ static SBasis divide_by_t0k(SBasis const &a, int k) {
 }
 
 static SBasis divide_by_t1k(SBasis const &a, int k) {
-    if(k < 0) {
-        SBasis c = Linear(1,0);
-        for (int i=2; i<=-k; i++){
-            c*=c;
-        }
-        c*=a;
-        return(c);
-    }else{
-        SBasis c = Linear(0,1);
-        for (int i=2; i<=k; i++){
-            c*=c;
-        }
-        c*=a;
-        return(divide_by_sk(c,k));
-    }
+    return divide_by_t0k(a, -k);
 }
 
 static D2<SBasis> RescaleForNonVanishingEnds(D2<SBasis> const &MM, double ZERO=1.e-4){
@@ -110,6 +101,14 @@ static D2<SBasis> RescaleForNonVanishingEnds(D2<SBasis> const &MM, double ZERO=1
     return M;
 }
 
+/*static D2<SBasis> RescaleForNonVanishing(D2<SBasis> const &MM, double ZERO=1.e-4){
+    std::vector<double> levels;
+    levels.push_back(-ZERO);
+    levels.push_back(ZERO);
+    //std::vector<std::vector<double> > mr = multi_roots(MM, levels);
+    }*/
+
+
 //=================================================================
 //TODO: what's this for?!?!
 Piecewise<D2<SBasis> > 
@@ -204,13 +203,15 @@ Geom::unitVector(D2<SBasis> const &V_in, double tol, unsigned order){
     D2<SBasis> V = RescaleForNonVanishingEnds(V_in);
     if (V[0].empty() && V[1].empty())
         return Piecewise<D2<SBasis> >(D2<SBasis>(Linear(1),SBasis()));
-    SBasis x = V[0], y = V[1], a, b;
+    SBasis x = V[0], y = V[1];
     SBasis r_eqn1, r_eqn2;
 
     Point v0 = unit_vector(V.at0());
     Point v1 = unit_vector(V.at1());
-    a.push_back(Linear(-v0[1],-v1[1]));
-    b.push_back(Linear( v0[0], v1[0]));
+    SBasis a(order+1, Linear());
+    a[0] = Linear(-v0[1],-v1[1]);
+    SBasis b(order+1, Linear());
+    b[0] = Linear( v0[0], v1[0]);
 
     r_eqn1 = -(a*x+b*y);
     r_eqn2 = Linear(1.)-(a*a+b*b);
@@ -231,8 +232,8 @@ Geom::unitVector(D2<SBasis> const &V_in, double tol, unsigned order){
         a1 = r1/dot(v1,V(1))*v1[0]-rr1/2*v1[1];
         b1 = r1/dot(v1,V(1))*v1[1]+rr1/2*v1[0];
 
-        a.push_back(Linear(a0,a1));        
-        b.push_back(Linear(b0,b1));
+        a[k] = Linear(a0,a1);
+        b[k] = Linear(b0,b1);
         //TODO: use "incremental" rather than explicit formulas.
         r_eqn1 = -(a*x+b*y);
         r_eqn2 = Linear(1)-(a*a+b*b);
@@ -561,10 +562,10 @@ std::vector<double>
 solve_lambda0(double a0,double a1,double c0,double c1,
              int insist_on_speeds_signs){
 
-    SBasis p;
-    p.push_back(Linear( a1*c0*c0+c1, a1*a0*(a0+ 2*c0) +a1*c0*c0 +c1 -1  ));
-    p.push_back(Linear( -a1*a0*(a0+2*c0), -a1*a0*(3*a0+2*c0) ));
-    p.push_back(Linear( a1*a0*a0 ));
+    SBasis p(3, Linear());
+    p[0] = Linear( a1*c0*c0+c1, a1*a0*(a0+ 2*c0) +a1*c0*c0 +c1 -1  );
+    p[1] = Linear( -a1*a0*(a0+2*c0), -a1*a0*(3*a0+2*c0) );
+    p[2] = Linear( a1*a0*a0 );
 
     OptInterval domain = find_bounds_for_lambda0(a0,a1,c0,c1,insist_on_speeds_signs);
     if ( !domain ) 
@@ -671,9 +672,11 @@ Geom::cubics_fitting_curvature(Point const &M0,   Point const &M1,
         Point V1 = lambda1[i]*dM1;
         D2<SBasis> cubic;
         for(unsigned dim=0;dim<2;dim++){
-            cubic[dim] = Linear(M0[dim],M1[dim]);
-            cubic[dim].push_back(Linear( M0[dim]-M1[dim]+V0[dim],
-                                         -M0[dim]+M1[dim]-V1[dim]));
+            SBasis c(2, Linear());
+            c[0] = Linear(M0[dim],M1[dim]);
+            c[1] = Linear( M0[dim]-M1[dim]+V0[dim],
+                           -M0[dim]+M1[dim]-V1[dim]);
+            cubic[dim] = c;
         }
 #if 0
            Piecewise<SBasis> k = curvature(result);

diff --git a/src/2geom/sbasis.cpp b/src/2geom/sbasis.cpp
index 0bd672c15f7b365af9d26503081b8da85101afb1..bdc40c9361051b3bc81300c6da8424e9609224cc 100644 (file)
--- a/src/2geom/sbasis.cpp
+++ b/src/2geom/sbasis.cpp
@@ -65,10 +65,10 @@ bool SBasis::isFinite() const {
 There is an elegant way to compute the value and n derivatives for a polynomial using a variant of horner's rule.  Someone will someday work out how for sbasis.
 */
 std::vector<double> SBasis::valueAndDerivatives(double t, unsigned n) const {
-    std::vector<double> ret(n);
-    ret.push_back(valueAt(t));
+    std::vector<double> ret(n+1);
+    ret[0] = valueAt(t);
     SBasis tmp = *this;
-    for(unsigned i = 0; i < n; i++) {
+    for(unsigned i = 1; i < n+1; i++) {
         tmp.derive();
         ret[i+1] = tmp.valueAt(t);
     }
@@ -82,18 +82,17 @@ std::vector<double> SBasis::valueAndDerivatives(double t, unsigned n) const {
 
 */
 SBasis operator+(const SBasis& a, const SBasis& b) {
-    SBasis result;
     const unsigned out_size = std::max(a.size(), b.size());
     const unsigned min_size = std::min(a.size(), b.size());
-    result.reserve(out_size);
+    SBasis result(out_size, Linear());
 
     for(unsigned i = 0; i < min_size; i++) {
-        result.push_back(a[i] + b[i]);
+        result[i] = a[i] + b[i];
     }
     for(unsigned i = min_size; i < a.size(); i++)
-        result.push_back(a[i]);
+        result[i] = a[i];
     for(unsigned i = min_size; i < b.size(); i++)
-        result.push_back(b[i]);
+        result[i] = b[i];
 
     assert(result.size() == out_size);
     return result;
@@ -105,18 +104,17 @@ SBasis operator+(const SBasis& a, const SBasis& b) {
 
 */
 SBasis operator-(const SBasis& a, const SBasis& b) {
-    SBasis result;
     const unsigned out_size = std::max(a.size(), b.size());
     const unsigned min_size = std::min(a.size(), b.size());
-    result.reserve(out_size);
+    SBasis result(out_size, Linear());
 
     for(unsigned i = 0; i < min_size; i++) {
-        result.push_back(a[i] - b[i]);
+        result[i] = a[i] - b[i];
     }
     for(unsigned i = min_size; i < a.size(); i++)
-        result.push_back(a[i]);
+        result[i] = a[i];
     for(unsigned i = min_size; i < b.size(); i++)
-        result.push_back(-b[i]);
+        result[i] = -b[i];
 
     assert(result.size() == out_size);
     return result;
@@ -130,12 +128,12 @@ SBasis operator-(const SBasis& a, const SBasis& b) {
 SBasis& operator+=(SBasis& a, const SBasis& b) {
     const unsigned out_size = std::max(a.size(), b.size());
     const unsigned min_size = std::min(a.size(), b.size());
-    a.reserve(out_size);
+    a.resize(out_size);
 
     for(unsigned i = 0; i < min_size; i++)
         a[i] += b[i];
     for(unsigned i = min_size; i < b.size(); i++)
-        a.push_back(b[i]);
+        a[i] = b[i];
 
     assert(a.size() == out_size);
     return a;
@@ -149,12 +147,12 @@ SBasis& operator+=(SBasis& a, const SBasis& b) {
 SBasis& operator-=(SBasis& a, const SBasis& b) {
     const unsigned out_size = std::max(a.size(), b.size());
     const unsigned min_size = std::min(a.size(), b.size());
-    a.reserve(out_size);
+    a.resize(out_size);
 
     for(unsigned i = 0; i < min_size; i++)
         a[i] -= b[i];
     for(unsigned i = min_size; i < b.size(); i++)
-        a.push_back(-b[i]);
+        a[i] = -b[i];
 
     assert(a.size() == out_size);
     return a;
@@ -166,10 +164,9 @@ SBasis& operator-=(SBasis& a, const SBasis& b) {
 
 */
 SBasis operator*(SBasis const &a, double k) {
-    SBasis c;
-    c.reserve(a.size());
+    SBasis c(a.size(), Linear());
     for(unsigned i = 0; i < a.size(); i++)
-        c.push_back(a[i] * k);
+        c[i] = a[i] * k;
     return c;
 }
 
@@ -195,12 +192,14 @@ SBasis& operator*=(SBasis& a, double b) {
 
 */
 SBasis shift(SBasis const &a, int sh) {
-    SBasis c = a;
-    if(sh > 0) {
-        c.insert(c.begin(), sh, Linear(0,0));
-    } else {
-        //TODO: truncate
-    }
+    size_t n = a.size()+sh;
+    SBasis c(n, Linear());
+    size_t m = std::max(0, sh);
+    
+    for(int i = 0; i < sh; i++)
+        c[i] = Linear(0,0);
+    for(size_t i = m, j = 0; i < n; i++, j++)
+        c[i] = a[j];
     return c;
 }
 
@@ -211,11 +210,13 @@ SBasis shift(SBasis const &a, int sh) {
 
 */
 SBasis shift(Linear const &a, int sh) {
-    SBasis c;
-    if(sh >= 0) {
-        c.insert(c.begin(), sh, Linear(0,0));
-        c.push_back(a);
-    }
+    size_t n = 1+sh;
+    SBasis c(n, Linear());
+    
+    for(int i = 0; i < sh; i++)
+        c[i] = Linear(0,0);
+    if(sh >= 0)
+        c[sh] = a;
     return c;
 }
 
@@ -230,8 +231,8 @@ SBasis multiply(SBasis const &a, SBasis const &b) {
     c.resize(a.size() + b.size(), Linear(0,0));
     for(unsigned j = 0; j < b.size(); j++) {
         for(unsigned i = j; i < a.size()+j; i++) {
-            double tri = Tri(b[j])*Tri(a[i-j]);
-            c[i+1/*shift*/] += Linear(Hat(-tri));
+            double tri = b[j].tri()*a[i-j].tri();
+            c[i+1/*shift*/] += Linear(-tri);
         }
     }
     for(unsigned j = 0; j < b.size(); j++) {
@@ -258,8 +259,8 @@ SBasis multiply_add(SBasis const &a, SBasis const &b, SBasis c) {
     c.resize(a.size() + b.size(), Linear(0,0));
     for(unsigned j = 0; j < b.size(); j++) {
         for(unsigned i = j; i < a.size()+j; i++) {
-            double tri = Tri(b[j])*Tri(a[i-j]);
-            c[i+1/*shift*/] += Linear(Hat(-tri));
+            double tri = b[j].tri()*a[i-j].tri();
+            c[i+1/*shift*/] += Linear(-tri);
         }
     }
     for(unsigned j = 0; j < b.size(); j++) {
@@ -296,12 +297,12 @@ SBasis integral(SBasis const &c) {
     a[0] = Linear(0,0);
 
     for(unsigned k = 1; k < c.size() + 1; k++) {
-        double ahat = -Tri(c[k-1])/(2*k);
-        a[k] = Hat(ahat);
+        double ahat = -c[k-1].tri()/(2*k);
+        a[k][0] = a[k][1] = ahat;
     }
     double aTri = 0;
     for(int k = c.size()-1; k >= 0; k--) {
-        aTri = (Hat(c[k]).d + (k+1)*aTri/2)/(2*k+1);
+        aTri = (c[k].hat() + (k+1)*aTri/2)/(2*k+1);
         a[k][0] -= aTri/2;
         a[k][1] += aTri/2;
     }
@@ -397,7 +398,7 @@ SBasis reciprocal(Linear const &a, int k) {
     SBasis c;
     assert(!a.isZero());
     c.resize(k, Linear(0,0));
-    double r_s0 = (Tri(a)*Tri(a))/(-a[0]*a[1]);
+    double r_s0 = (a.tri()*a.tri())/(-a[0]*a[1]);
     double r_s0k = 1;
     for(unsigned i = 0; i < (unsigned)k; i++) {
         c[i] = Linear(r_s0k/a[0], r_s0k/a[1]);
@@ -444,7 +445,7 @@ SBasis compose(SBasis const &a, SBasis const &b) {
     SBasis r;
 
     for(int i = a.size()-1; i >= 0; i--) {
-        r = multiply_add(r, s, SBasis(Linear(Hat(a[i][0]))) - b*a[i][0] + b*a[i][1]);
+        r = multiply_add(r, s, SBasis(Linear(a[i][0])) - b*a[i][0] + b*a[i][1]);
     }
     return r;
 }
@@ -460,7 +461,7 @@ SBasis compose(SBasis const &a, SBasis const &b, unsigned k) {
     SBasis r;
 
     for(int i = a.size()-1; i >= 0; i--) {
-        r = multiply_add(r, s, SBasis(Linear(Hat(a[i][0]))) - b*a[i][0] + b*a[i][1]);
+        r = multiply_add(r, s, SBasis(Linear(a[i][0])) - b*a[i][0] + b*a[i][1]);
     }
     r.truncate(k);
     return r;
@@ -499,11 +500,11 @@ SBasis inverse(SBasis a, int k) {
     if(a1 != 1) {
         a /= a1;
     }
-    SBasis c;                           // c(v) := 0
+    SBasis c(k, Linear());                           // c(v) := 0
     if(a.size() >= 2 && k == 2) {
-        c.push_back(Linear(0,1));
+        c[0] = Linear(0,1);
         Linear t1(1+a[1][0], 1-a[1][1]);    // t_1
-        c.push_back(Linear(-a[1][0]/t1[0], -a[1][1]/t1[1]));
+        c[1] = Linear(-a[1][0]/t1[0], -a[1][1]/t1[1]);
     } else if(a.size() >= 2) {                      // non linear
         SBasis r = Linear(0,1);             // r(u) := r_0(u) := u
         Linear t1(1./(1+a[1][0]), 1./(1-a[1][1]));    // 1./t_1
@@ -519,7 +520,7 @@ SBasis inverse(SBasis a, int k) {
         //assert(t1 == t[1]);
 #endif
 
-        c.resize(k+1, Linear(0,0));
+        //c.resize(k+1, Linear(0,0));
         for(unsigned i = 0; i < (unsigned)k; i++) {   // for i:=0 to k do
 #ifdef DEBUG_INVERSION
             std::cout << "-------" << i << ": ---------" <<std::endl;
@@ -566,10 +567,11 @@ SBasis inverse(SBasis a, int k) {
 It is recommended to use the piecewise version unless you have good reason.
 */
 SBasis sin(Linear b, int k) {
-    SBasis s = Linear(std::sin(b[0]), std::sin(b[1]));
-    Tri tr(s[0]);
-    double t2 = Tri(b);
-    s.push_back(Linear(std::cos(b[0])*t2 - tr, -std::cos(b[1])*t2 + tr));
+    SBasis s(k+2, Linear());
+    s[0] = Linear(std::sin(b[0]), std::sin(b[1]));
+    double tr = s[0].tri();
+    double t2 = b.tri();
+    s[1] = Linear(std::cos(b[0])*t2 - tr, -std::cos(b[1])*t2 + tr);
 
     t2 *= t2;
     for(int i = 0; i < k; i++) {
@@ -578,7 +580,7 @@ SBasis sin(Linear b, int k) {
         bo -= s[i]*(t2/(i+1));
 
 
-        s.push_back(bo/double(i+2));
+        s[i+2] = bo/double(i+2);
     }
 
     return s;
@@ -605,7 +607,7 @@ TODO: compute order according to tol?
 TODO: requires g(0)=0 & g(1)=1 atm... adaptation to other cases should be obvious!
 */
 SBasis compose_inverse(SBasis const &f, SBasis const &g, unsigned order, double zero){
-    SBasis result; //result
+    SBasis result(order, Linear()); //result
     SBasis r=f; //remainder
     SBasis Pk=Linear(1)-g,Qk=g,sg=Pk*Qk;
     Pk.truncate(order);
@@ -634,7 +636,7 @@ SBasis compose_inverse(SBasis const &f, SBasis const &g, unsigned order, double
             a=( q01*r10-q10*r01)/det;
             b=(-p01*r10+p10*r01)/det;
         }
-        result.push_back(Linear(a,b));
+        result[k] = Linear(a,b);
         r=r-Pk*a-Qk*b;
 
         Pk=Pk*sg;

diff --git a/src/2geom/sbasis.h b/src/2geom/sbasis.h
index f2681a722f644680e7a49d0c625e5ee551e61e5b..3038c17762458994547b5f22f56467a79b081f7e 100644 (file)
--- a/src/2geom/sbasis.h
+++ b/src/2geom/sbasis.h
@@ -43,24 +43,62 @@
 #include <2geom/utils.h>
 #include <2geom/exception.h>
 
+//#define USE_SBASISN 1
 
-#ifdef USE_SBASIS_OF
+
+#if defined(USE_SBASIS_OF)
 
 #include "sbasis-of.h"
 
+#elif defined(USE_SBASISN)
+
+#include "sbasisN.h"
+namespace Geom{
+
+/*** An empty SBasis is identically 0. */
+class SBasis : public SBasisN<1>;
+
+};
 #else
 
 namespace Geom{
 
 /*** An empty SBasis is identically 0. */
-class SBasis : public std::vector<Linear>{
+class SBasis{
+    std::vector<Linear> d;
+    void push_back(Linear const&l) { d.push_back(l); }
+
 public:
+    size_t size() const {return d.size();}
+    Linear operator[](unsigned i) const {
+        assert(i < size());
+        return d[i];
+    }
+    Linear& operator[](unsigned i) { return d.at(i); }
+    Linear const* begin() const { return (Linear const*)&*d.begin();}
+    Linear const* end() const { return (Linear const*)&*d.end();}
+    Linear* begin() { return (Linear*)&*d.begin();}
+    Linear* end() { return (Linear*)&*d.end();}
+    bool empty() const {return d.empty();}
+    Linear &back() {return d.back();}
+    Linear const &back() const {return d.back();}
+    void pop_back() { d.pop_back();}
+    void resize(unsigned n) { d.resize(n);}
+    void resize(unsigned n, Linear const& l) { d.resize(n, l);}
+    void reserve(unsigned n) { d.reserve(n);}
+    void clear() {d.clear();}
+    void insert(Linear* before, const Linear* src_begin, const Linear* src_end) { d.insert(std::vector<Linear>::iterator(before), src_begin, src_end);}
+    //void insert(Linear* aa, Linear* bb, Linear* cc} { d.insert(aa, bb, cc);}
+    Linear& at(unsigned i) { return d.at(i);}
+    //void insert(Linear* before, int& n, Linear const &l) { d.insert(std::vector<Linear>::iterator(before), n, l);}
+
+    
     SBasis() {}
     explicit SBasis(double a) {
         push_back(Linear(a,a));
     }
     SBasis(SBasis const & a) :
-        std::vector<Linear>(a)
+        d(a.d)
     {}
     SBasis(Linear const & bo) {
         push_back(bo);
@@ -68,6 +106,7 @@ public:
     SBasis(Linear* bo) {
         push_back(*bo);
     }
+    explicit SBasis(size_t n, Linear const&l) : d(n, l) {}
 
     //IMPL: FragmentConcept
     typedef double output_type;
@@ -93,6 +132,8 @@ public:
     inline double at1() const{
         if(empty()) return 0; else return (*this)[0][1];
     }
+    
+    int degreesOfFreedom() const { return size()*2;}
 
     double valueAt(double t) const {
         double s = t*(1-t);
@@ -119,14 +160,7 @@ public:
 // compute f(g)
     SBasis operator()(SBasis const & g) const;
 
-    Linear operator[](unsigned i) const {
-        assert(i < size());
-        return std::vector<Linear>::operator[](i);
-    }
-
 //MUTATOR PRISON
-    Linear& operator[](unsigned i) { return this->at(i); }
-
     //remove extra zeros
     void normalize() {
         while(!empty() && 0 == back()[0] && 0 == back()[1])
@@ -152,21 +186,20 @@ OptInterval bounds_local(SBasis const &a, const OptInterval &t, int order = 0);
 useful for reversing a parameteric curve.
 */
 inline SBasis reverse(SBasis const &a) {
-    SBasis result;
-    result.reserve(a.size());
+    SBasis result(a.size(), Linear());
+    
     for(unsigned k = 0; k < a.size(); k++)
-       result.push_back(reverse(a[k]));
+        result[k] = reverse(a[k]);
     return result;
 }
 
 //IMPL: ScalableConcept
 inline SBasis operator-(const SBasis& p) {
     if(p.isZero()) return SBasis();
-    SBasis result;
-    result.reserve(p.size());
+    SBasis result(p.size(), Linear());
         
     for(unsigned i = 0; i < p.size(); i++) {
-        result.push_back(-p[i]);
+        result[i] = -p[i];
     }
     return result;
 }
@@ -183,7 +216,7 @@ SBasis& operator+=(SBasis& a, const SBasis& b);
 SBasis& operator-=(SBasis& a, const SBasis& b);
 
 //TODO: remove?
-inline SBasis operator+(const SBasis & a, Linear const & b) {
+/*inline SBasis operator+(const SBasis & a, Linear const & b) {
     if(b.isZero()) return a;
     if(a.isZero()) return b;
     SBasis result(a);
@@ -209,7 +242,7 @@ inline SBasis& operator-=(SBasis& a, const Linear& b) {
     else
         a[0] -= b;
     return a;
-}
+    }*/
 
 //IMPL: OffsetableConcept
 inline SBasis operator+(const SBasis & a, double b) {
@@ -226,14 +259,14 @@ inline SBasis operator-(const SBasis & a, double b) {
 }
 inline SBasis& operator+=(SBasis& a, double b) {
     if(a.isZero())
-        a.push_back(Linear(b,b));
+        a = SBasis(Linear(b,b));
     else
         a[0] += b;
     return a;
 }
 inline SBasis& operator-=(SBasis& a, double b) {
     if(a.isZero())
-        a.push_back(Linear(-b,-b));
+        a = SBasis(Linear(-b,-b));
     else
         a[0] -= b;
     return a;
@@ -300,6 +333,7 @@ SBasis compose_inverse(SBasis const &f, SBasis const &g, unsigned order=2, doubl
 
 */
 inline SBasis portion(const SBasis &t, double from, double to) { return compose(t, Linear(from, to)); }
+inline SBasis portion(const SBasis &t, Interval ivl) { return compose(t, Linear(ivl[0], ivl[1])); }
 
 // compute f(g)
 inline SBasis

diff --git a/src/2geom/svg-elliptical-arc.h b/src/2geom/svg-elliptical-arc.h
index 92ec51b49b800bc5de0ce24e21a8f8ee2b901281..dad9000c11034758512383b920f9c08b1022a29f 100644 (file)
--- a/src/2geom/svg-elliptical-arc.h
+++ b/src/2geom/svg-elliptical-arc.h
@@ -248,6 +248,8 @@ class SVGEllipticalArc : public Curve
         else
             return SBasisCurve(toSBasis()).winding(p);
     }
+    
+    int degreesOfFreedom() const { return 5;}
 
     Curve *derivative() const;
 

diff --git a/src/2geom/sweep.cpp b/src/2geom/sweep.cpp
index 53b9532558ef7e1af19c92f2e32d7894e86eb585..227c47e098f41da0a621cc78ed0ac175ba9773c4 100644 (file)
--- a/src/2geom/sweep.cpp
+++ b/src/2geom/sweep.cpp
@@ -4,6 +4,15 @@
 
 namespace Geom {
 
+/**
+ * \brief Make a list of pairs of self intersections in a list of Rects.
+ * 
+ * \param rs: vector of Rect.
+ *
+ * [(A = rs[i], B = rs[j]) for i,J in enumerate(pairs) for j in J]
+ * then A.left <= B.left
+ */
+
 std::vector<std::vector<unsigned> > sweep_bounds(std::vector<Rect> rs) {
     std::vector<Event> events; events.reserve(rs.size()*2);
     std::vector<std::vector<unsigned> > pairs(rs.size());
@@ -34,6 +43,15 @@ std::vector<std::vector<unsigned> > sweep_bounds(std::vector<Rect> rs) {
     return pairs;
 }
 
+/**
+ * \brief Make a list of pairs of red-blue intersections between two lists of Rects.
+ * 
+ * \param a: vector of Rect.
+ * \param b: vector of Rect.
+ *
+ * [(A = rs[i], B = rs[j]) for i,J in enumerate(pairs) for j in J]
+ * then A.left <= B.left, A in a, B in b
+ */
 std::vector<std::vector<unsigned> > sweep_bounds(std::vector<Rect> a, std::vector<Rect> b) {
     std::vector<std::vector<unsigned> > pairs(a.size());
     if(a.empty() || b.empty()) return pairs;