Update to 2Geom rev. 1113

(bzr r3622)

45 files changed, 3382 insertions, 1168 deletions

diff --git a/src/2geom/Makefile_insert b/src/2geom/Makefile_insert
index a100c418d..4a87c1a93 100644
--- a/src/2geom/Makefile_insert
+++ b/src/2geom/Makefile_insert
@@ -9,6 +9,8 @@
 	2geom/bezier-to-sbasis.h	\

 	2geom/bezier-utils.cpp	\

 	2geom/bezier-utils.h	\

+	2geom/chebyshev.cpp	\

+	2geom/chebyshev.h	\

 	2geom/choose.h	\

 	2geom/circle-circle.cpp	\

 	2geom/circulator.h	\

@@ -18,6 +20,10 @@
 	2geom/convex-cover.cpp	\

 	2geom/convex-cover.h	\

 	2geom/coord.h	\

+	2geom/crossing.cpp	\

+	2geom/crossing.cpp	\

+	2geom/d2-sbasis.cpp	\

+	2geom/d2-sbasis.h	\

 	2geom/d2.cpp	\

 	2geom/d2.h	\

 	2geom/geom.cpp	\

@@ -27,6 +33,9 @@
 	2geom/linear.h	\

 	2geom/matrix.cpp	\

 	2geom/matrix.h	\

+	2geom/ord.h	\

+	2geom/path-intersection.cpp	\

+	2geom/path-intersection.h	\

 	2geom/path.cpp	\

 	2geom/path.h	\

 	2geom/piecewise.cpp	\

@@ -44,6 +53,8 @@
 	2geom/quadtree.cpp	\

 	2geom/quadtree.h	\

 	2geom/rect.h	\

+	2geom/region.cpp	\

+	2geom/region.h	\

 	2geom/sbasis-2d.cpp	\

 	2geom/sbasis-2d.h	\

 	2geom/sbasis.cpp	\

@@ -57,6 +68,8 @@
 	2geom/sbasis-roots.cpp	\

 	2geom/sbasis-to-bezier.cpp	\

 	2geom/sbasis-to-bezier.h	\

+	2geom/shape.cpp	\

+	2geom/shape.h	\

 	2geom/solve-bezier-one-d.cpp	\

 	2geom/solve-bezier-parametric.cpp	\

 	2geom/solver.h	\

@@ -65,6 +78,8 @@
 	2geom/svg-path.h	\

 	2geom/svg-path-parser.cpp	\

 	2geom/svg-path-parser.h	\

+	2geom/sweep.cpp	\

+	2geom/sweep.h	\

 	2geom/transforms.cpp	\

 	2geom/transforms.h	\

 	2geom/utils.h	\

diff --git a/src/2geom/basic-intersection.cpp b/src/2geom/basic-intersection.cpp
index a5b827023..694760d5a 100644
--- a/src/2geom/basic-intersection.cpp
+++ b/src/2geom/basic-intersection.cpp
@@ -159,7 +159,7 @@ void OldBezier::split(double t, OldBezier &left, OldBezier &right) const {
  *	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 OldBezier 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.
diff --git a/src/2geom/bezier.h b/src/2geom/bezier.h
index 7949e006c..2c1b49213 100644
--- a/src/2geom/bezier.h
+++ b/src/2geom/bezier.h
@@ -3,6 +3,7 @@
  *
  * Copyright 2007  MenTaLguY <mental@rydia.net>
  * Copyright 2007  Michael Sloan <mgsloan@gmail.com>
+ * Copyright 2007  Nathan Hurst <njh@njhurst.com>
  *
  * This library is free software; you can redistribute it and/or
  * modify it either under the terms of the GNU Lesser General Public
@@ -35,107 +36,246 @@
 #include "coord.h"
 #include "isnan.h"
 #include "bezier-to-sbasis.h"
+#include "d2.h"
+#include "solver.h"
+#include <boost/optional/optional.hpp>
 
 namespace Geom {
 
+template<unsigned order>
+Coord subdivideArr(Coord t, Coord const *v, Coord *left, Coord *right) {
+    Coord vtemp[order+1][order+1];
+
+    /* Copy control points	*/
+    std::copy(v, v+order+1, vtemp[0]);
+
+    /* Triangle computation	*/
+    for (unsigned i = 1; i <= order; i++) {	
+        for (unsigned j = 0; j <= order - i; j++) {
+            vtemp[i][j] = lerp(t, vtemp[i-1][j], vtemp[i-1][j+1]);
+        }
+    }
+    if(left != NULL)
+        for (unsigned j = 0; j <= order; j++)
+            left[j]  = vtemp[j][0];
+    if(right != NULL)
+        for (unsigned j = 0; j <= order; j++)
+            right[j] = vtemp[order-j][j];
+
+    return (vtemp[order][0]);
+}
+
+
 template <unsigned order>
 class Bezier {
 private:
-  Coord c_[order+1];
+    Coord c_[order+1];
+
+    template<unsigned ord>
+    friend Bezier<ord> portion(const Bezier<ord> & a, Coord from, Coord to);
+
+    template<unsigned ord>
+    friend Interval bounds_fast(Bezier<ord> const & b);
+
+    template<unsigned ord>
+    friend Bezier<ord-1> derivative(const Bezier<ord> & a);
 
 protected:
-  Bezier(Coord c[]) {
-    std::copy(c, c+order+1, c_);
-  }
+    Bezier(Coord const c[]) {
+        std::copy(c, c+order+1, c_);
+    }
 
 public:
-  template <unsigned required_order>
-  static void assert_order(Bezier<required_order> const *) {}
-
-  Bezier() {}
-
-  //Construct an order-0 bezier (constant Bézier)
-  explicit Bezier(Coord c0) {
-    assert_order<0>(this);
-    c_[0] = c0;
-  }
-
-  //Construct an order-1 bezier (linear Bézier)
-  Bezier(Coord c0, Coord c1) {
-    assert_order<1>(this);
-    c_[0] = c0; c_[1] = c1;
-  }
-
-  //Construct an order-2 bezier (quadratic Bézier)
-  Bezier(Coord c0, Coord c1, Coord c2) {
-    assert_order<2>(this);
-    c_[0] = c0; c_[1] = c1; c_[2] = c2;
-  }
-
-  //Construct an order-3 bezier (cubic Bézier)
-  Bezier(Coord c0, Coord c1, Coord c2, Coord c3) {
-    assert_order<3>(this);
-    c_[0] = c0; c_[1] = c1; c_[2] = c2; c_[3] = c3;
-  }
-
-  inline unsigned degree() const { return order; }
-
-  //IMPL: FragmentConcept
-  typedef Coord output_type;
-  inline bool isZero() const { 
-     for(int i = 0; i <= order; i++) {
-       if(c_[i] != 0) return false;
-     }
-     return true;
-  }
-  inline bool isFinite() const {
-    for(int i = 0; i <= order; i++) {
-      if(!is_finite(c_[i])) return false;
+
+    //TODO: fix this assert - get it compiling!
+    //template <unsigned required_order>
+    //static void assert_order(Bezier<required_order> const *) {}
+
+    Bezier() {}
+
+    //Construct an order-0 bezier (constant Bézier)
+    explicit Bezier(Coord c0) {
+        //assert_order<0>(this);
+        c_[0] = c0;
+    }
+
+    //Construct an order-1 bezier (linear Bézier)
+    Bezier(Coord c0, Coord c1) {
+        //assert_order<1>(this);
+        c_[0] = c0; c_[1] = c1;
+    }
+
+    //Construct an order-2 bezier (quadratic Bézier)
+    Bezier(Coord c0, Coord c1, Coord c2) {
+        //assert_order<2>(this);
+        c_[0] = c0; c_[1] = c1; c_[2] = c2;
+    }
+
+    //Construct an order-3 bezier (cubic Bézier)
+    Bezier(Coord c0, Coord c1, Coord c2, Coord c3) {
+        //assert_order<3>(this);
+        c_[0] = c0; c_[1] = c1; c_[2] = c2; c_[3] = c3;
     }
-    return true;
-  }
-  inline Coord at0() const { return c_[0]; }
-  inline Coord at1() const { return c_[order]; }
 
-  inline SBasis toSBasis() const { return bezier_to_sbasis<order>(c_); }
+    inline unsigned degree() const { return order; }
 
-  inline Interval bounds_fast() const { return Interval::fromArray(c_, order+1); }
-  //TODO: better bounds exact
-  inline Interval bounds_exact() const { return toSBasis().bounds_exact(); }
-  inline Interval bounds_local(double u, double v) const { return toSBasis().bounds_local(u, v); }
+    //IMPL: FragmentConcept
+    typedef Coord output_type;
+    inline bool isZero() const { 
+        for(unsigned i = 0; i <= order; i++) {
+            if(c_[i] != 0) return false;
+        }
+        return true;
+    }
+    inline bool isFinite() const {
+        for(unsigned i = 0; i <= order; i++) {
+            if(!is_finite(c_[i])) return false;
+        }
+        return true;
+    }
+    inline Coord at0() const { return c_[0]; }
+    inline Coord at1() const { return c_[order]; }
 
-  //Only mutator
-  inline Coord &operator[](int index) { return c_[index]; }
-  inline Coord const &operator[](int index) const { return c_[index]; }
+    inline Coord valueAt(double t) const { return subdivideArr<order>(t, c_, NULL, NULL); }
+    inline Coord operator()(double t) const { return valueAt(t); }
 
-  Maybe<int> winding(Point p) const {
-    return sbasis_winding(toSBasis(), p);
-  }
+    inline SBasis toSBasis() const { return bezier_to_sbasis<order>(c_); }
+
+    //Only mutator
+    inline Coord &operator[](unsigned ix) { return c_[ix]; }
+    inline Coord const &operator[](unsigned ix) const { return c_[ix]; }
+    inline void setPoint(unsigned ix, double val) { c_[ix] = val; }
+    
+    /* This is inelegant, as it uses several extra stores.  I think there might be a way to
+     * evaluate roughly in situ. */
+    
+    std::vector<Coord> valueAndDerivatives(Coord t, unsigned n_derivs) const {
+        throw NotImplemented();
+        // Can't be implemented without a dynamic version of subdivide.
+        /*std::vector<Coord> val_n_der;
+        Coord d_[order+1];
+        for(int di = n_derivs; di > 0; di--) {
+            val_n_der.push_back(subdivideArr<di>(t, d_, NULL, NULL));
+            for(unsigned i = 0; i < di; i++) {
+                d[i] = order*(a._c[i+1] - a._c[i]);
+            }
+            }*/
+    }
   
-  Point pointAndDerivativesAt(Coord t, unsigned n_derivs, Point *derivs) const {
-    //TODO
-    return Point(0,0);
-  }
+    std::pair<Bezier<order>, Bezier<order> > subdivide(Coord t) const {
+        Bezier<order> a, b;
+        subdivideArr(t, order, c_, a.c_, b.c_);
+        return std::pair<Bezier<order>, Bezier<order> >(a, b);
+    }
+
+    std::vector<double> roots() const {
+        std::vector<double> solutions;
+        find_bernstein_roots(c_, order, solutions, 0, 0.0, 1.0);
+        return solutions;
+    }
 };
 
+//TODO: implement others
+template<unsigned order>
+Bezier<order> operator+(const Bezier<order> & a, double v) {
+    Bezier<order> result;
+    for(unsigned i = 0; i <= order; i++)
+        result[i] = a[i] + v;
+    return result;
+}
+template<unsigned order>
+Bezier<order> operator-(const Bezier<order> & a, double v) {
+    Bezier<order> result;
+    for(unsigned i = 0; i <= order; i++)
+        result[i] = a[i] - v;
+    return result;
+}
+template<unsigned order>
+Bezier<order> operator*(const Bezier<order> & a, double v) {
+    Bezier<order> result;
+    for(unsigned i = 0; i <= order; i++)
+        result[i] = a[i] * v;
+    return result;
+}
+template<unsigned order>
+Bezier<order> operator/(const Bezier<order> & a, double v) {
+    Bezier<order> result;
+    for(unsigned i = 0; i <= order; i++)
+        result[i] = a[i] / v;
+    return result;
+}
+
 template<unsigned order>
 Bezier<order> reverse(const Bezier<order> & a) {
-  Bezier<order> result;
-  for(int i = 0; i <= order; i++)
-    result[i] = a[order - i];
-  return result;
+    Bezier<order> result;
+    for(unsigned i = 0; i <= order; i++)
+        result[i] = a[order - i];
+    return result;
+}
+
+template<unsigned order>
+Bezier<order> portion(const Bezier<order> & a, double from, double to) {
+    //TODO: implement better?
+    Coord res[order+1];
+    if(from == 0) {
+        if(to == 1) { return Bezier<order>(a); }
+        subdivideArr<order>(to, a.c_, res, NULL);
+        return Bezier<order>(res);
+    }
+    subdivideArr<order>(from, a.c_, NULL, res);
+    if(to == 1) return Bezier<order>(res);
+    Coord res2[order+1];
+    subdivideArr<order>((to - from)/(1 - from), res, res2, NULL);
+    return Bezier<order>(res2);
+}
+
+template<unsigned order>
+std::vector<Point> bezier_points(const D2<Bezier<order> > & a) {
+    std::vector<Point> result;
+    for(unsigned i = 0; i <= order; i++) {
+        Point p;
+        for(unsigned d = 0; d < 2; d++) p[d] = a[d][i];
+        result.push_back(p);
+    }
+    return result;
+}
+
+template<unsigned order>
+Bezier<order-1> derivative(const Bezier<order> & a) {
+    Bezier<order-1> der;
+    
+    for(unsigned i = 0; i < order; i++) {
+        der.c_[i] = order*(a.c_[i+1] - a.c_[i]);
+    }
+    return der;
 }
 
 template<unsigned order>
-vector<Point> bezier_points(const D2<Bezier<order> > & a) {
-  vector<Point> result;
-  for(int i = 0; i <= order; i++) {
-    Point p;
-    for(unsigned d = 0; d < 2; d++) p[d] = a[d][i];
-    result[i] = p;
-  }
-  return result;
+inline Interval bounds_fast(Bezier<order> const & b) {
+    return Interval::fromArray(b.c_, order+1);
+}
+
+//TODO: better bounds exact
+template<unsigned order>
+inline Interval bounds_exact(Bezier<order> const & b) {
+    return bounds_exact(b.toSBasis());
+}
+
+template<unsigned order>
+inline Interval bounds_local(Bezier<order> const & b, Interval i) {
+    return bounds_fast(portion(b, i.min(), i.max()));
+    //return bounds_local(b.toSBasis(), i);
 }
 
 }
 #endif //SEEN_BEZIER_H
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
diff --git a/src/2geom/chebyshev.cpp b/src/2geom/chebyshev.cpp
new file mode 100644
index 000000000..b56c57eac
--- /dev/null
+++ b/src/2geom/chebyshev.cpp
@@ -0,0 +1,126 @@
+#include "chebyshev.h"
+
+#include "sbasis.h"
+#include "sbasis-poly.h"
+
+#include <vector>
+using std::vector;
+
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_chebyshev.h>
+
+namespace Geom{
+
+SBasis cheb(unsigned n) {
+    static std::vector<SBasis> basis;
+    if(basis.empty()) {
+        basis.push_back(Linear(1,1));
+        basis.push_back(Linear(0,1));
+    }
+    for(unsigned i = basis.size(); i <= n; i++) {
+        basis.push_back(Linear(0,2)*basis[i-1] - basis[i-2]);
+    }
+    
+    return basis[n];
+}
+
+SBasis cheb_series(unsigned n, double* cheb_coeff) {
+    SBasis r;
+    for(unsigned i = 0; i < n; i++) {
+        double cof = cheb_coeff[i];
+        //if(i == 0)
+        //cof /= 2;
+        r += cheb(i)*cof;
+    }
+    
+    return r;
+}
+
+SBasis clenshaw_series(unsigned m, double* cheb_coeff) {
+    /** b_n = a_n
+        b_n-1 = 2*x*b_n + a_n-1
+        b_n-k = 2*x*b_{n-k+1} + a_{n-k} - b_{n - k + 2}
+        b_0 = x*b_1 + a_0 - b_2
+    */
+    
+    double a = -1, b = 1;
+    SBasis d, dd;
+    SBasis y = (Linear(0, 2) - (a+b)) / (b-a);
+    SBasis y2 = 2*y;
+    for(int j = m - 1; j >= 1; j--) {
+        SBasis sv = d;
+        d = y2*d - dd + cheb_coeff[j];
+        dd = sv;
+    }
+    
+    return y*d - dd + 0.5*cheb_coeff[0];
+}
+
+SBasis chebyshev_approximant (double (*f)(double,void*), int order, Interval in, void* p) {
+    gsl_cheb_series *cs = gsl_cheb_alloc (order+2);
+
+    gsl_function F;
+
+    F.function = f;
+    F.params = p;
+
+    gsl_cheb_init (cs, &F, in[0], in[1]);
+    
+    SBasis r = compose(clenshaw_series(order, cs->c), Linear(-1,1));
+        
+    gsl_cheb_free (cs);
+    return r;
+}
+
+struct wrap {
+    double (*f)(double,void*);
+    void* pp;
+    double fa, fb;
+    Interval in;
+};
+
+double f_interp(double x, void* p) {
+    struct wrap *wr = (struct wrap *)p;
+    double z = (x - wr->in[0]) / (wr->in[1] - wr->in[0]);
+    return (wr->f)(x, wr->pp) - ((1 - z)*wr->fa + z*wr->fb);
+}
+
+SBasis chebyshev_approximant_interpolating (double (*f)(double,void*), 
+                                            int order, Interval in, void* p) {
+    double fa = f(in[0], p);
+    double fb = f(in[1], p);
+    struct wrap wr;
+    wr.fa = fa;
+    wr.fb = fb;
+    wr.in = in;
+    printf("%f %f\n", fa, fb);
+    wr.f = f;
+    wr.pp = p;
+    return compose(Linear(in[0], in[1]), Linear(fa, fb)) + chebyshev_approximant(f_interp, order, in, &wr) + Linear(fa, fb);
+}
+
+SBasis chebyshev(unsigned n) {
+    static std::vector<SBasis> basis;
+    if(basis.empty()) {
+        basis.push_back(Linear(1,1));
+        basis.push_back(Linear(0,1));
+    }
+    for(unsigned i = basis.size(); i <= n; i++) {
+        basis.push_back(Linear(0,2)*basis[i-1] - basis[i-2]);
+    }
+    
+    return basis[n];
+}
+
+};
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
diff --git a/src/2geom/chebyshev.h b/src/2geom/chebyshev.h
new file mode 100644
index 000000000..f7e82a5df
--- /dev/null
+++ b/src/2geom/chebyshev.h
@@ -0,0 +1,30 @@
+#ifndef _CHEBYSHEV
+#define _CHEBYSHEV
+
+#include "sbasis.h"
+#include "interval.h"
+
+/*** Conversion between Chebyshev approximation and SBasis.
+ * 
+ */
+
+namespace Geom{
+
+SBasis chebyshev_approximant (double (*f)(double,void*), int order, Interval in, void* p=0);
+SBasis chebyshev_approximant_interpolating (double (*f)(double,void*), int order, Interval in, void* p=0);
+SBasis chebyshev(unsigned n);
+
+};
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
+
+#endif
diff --git a/src/2geom/circulator.h b/src/2geom/circulator.h
index 4c7f83bad..512aac39f 100644
--- a/src/2geom/circulator.h
+++ b/src/2geom/circulator.h
@@ -38,18 +38,18 @@ namespace Geom {
 template <typename Iterator>
 class Circulator {
 public:
-    typedef std::random_access_iterator_tag std::iterator_category;
-    typedef std::iterator_traits<Iterator>::value_type value_type;
-    typedef std::iterator_traits<Iterator>::difference_type difference_type;
-    typedef std::iterator_traits<Iterator>::pointer pointer;
-    typedef std::iterator_traits<Iterator>::reference reference;
+    typedef std::random_access_iterator_tag iterator_category;
+    typedef typename std::iterator_traits<Iterator>::value_type value_type;
+    typedef typename std::iterator_traits<Iterator>::difference_type difference_type;
+    typedef typename std::iterator_traits<Iterator>::pointer pointer;
+    typedef typename std::iterator_traits<Iterator>::reference reference;
 
     Circulator(Iterator const &first,
                Iterator const &last,
                Iterator const &pos)
     : _first(first), _last(last), _pos(pos)
     {
-        match_random_access(std::iterator_category(first));
+        match_random_access(iterator_category(first));
     }
 
     reference operator*() const {
@@ -101,12 +101,12 @@ public:
         return _pos - other._pos;
     }
 
-    reference operator[n] const {
+    reference operator[](int n) const {
         return *_offset(n);
     }
 
 private:
-    void match_random_access(random_access_iterator_tag) {}
+    void match_random_access(iterator_category) {}
 
     Iterator _offset(int n) {
         difference_type range=( _last - _first );
diff --git a/src/2geom/concepts.h b/src/2geom/concepts.h
index 42192c445..bd020ab66 100644
--- a/src/2geom/concepts.h
+++ b/src/2geom/concepts.h
@@ -34,7 +34,7 @@
 #include "sbasis.h"
 #include "interval.h"
 #include "point.h"
-
+#include <vector>
 #include <boost/concept_check.hpp>
 
 namespace Geom {
@@ -66,6 +66,9 @@ struct FragmentConcept {
     bool b;
     BoundsType i;
     Interval dom;
+    std::vector<OutputType> v;
+    unsigned u;
+    SbType sb;
     void constraints() {
         t = T(o);
         b = t.isZero();
@@ -74,7 +77,10 @@ struct FragmentConcept {
         o = t.at1();
         o = t.valueAt(d);
         o = t(d);
-        SbType sb = t.toSBasis();
+        v = t.valueAndDerivatives(d, u);
+		//Is a pure derivative (ignoring others) accessor ever much faster?
+		//u = number of values returned. first val is value.
+        sb = t.toSBasis();
         t = reverse(t);
         i = bounds_fast(t);
         i = bounds_exact(t);
diff --git a/src/2geom/crossing.cpp b/src/2geom/crossing.cpp
new file mode 100644
index 000000000..8f29635e8
--- /dev/null
+++ b/src/2geom/crossing.cpp
@@ -0,0 +1,81 @@
+#include "crossing.h"
+
+namespace Geom {
+
+void merge_crossings(Crossings &a, Crossings &b, unsigned i) {
+    Crossings n;
+    sort_crossings(b, i);
+    n.resize(a.size() + b.size());
+    std::merge(a.begin(), a.end(), b.begin(), b.end(), n.begin(), CrossingOrder(i));
+    a = n;
+}
+
+void offset_crossings(Crossings &cr, double a, double b) {
+    for(unsigned i = 0; i < cr.size(); i++) {
+        cr[i].ta += a;
+        cr[i].tb += b;
+    }
+}
+
+Crossings reverse_ta(Crossings const &cr, std::vector<double> max) {
+    Crossings ret;
+    for(Crossings::const_iterator i = cr.begin(); i != cr.end(); ++i) {
+        double mx = max[i->a];
+        ret.push_back(Crossing(i->ta > mx+0.01 ? (1 - (i->ta - mx) + mx) : mx - i->ta,
+                               i->tb, !i->dir));
+    }
+    return ret;
+}
+
+Crossings reverse_tb(Crossings const &cr, unsigned split, std::vector<double> max) {
+    Crossings ret;
+    for(Crossings::const_iterator i = cr.begin(); i != cr.end(); ++i) {
+        double mx = max[i->b - split];
+        std::cout << i->b << "\n";
+        ret.push_back(Crossing(i->ta, i->tb > mx+0.01 ? (1 - (i->tb - mx) + mx) : mx - i->tb,
+                               !i->dir));
+    }
+    return ret;
+}
+
+CrossingSet reverse_ta(CrossingSet const &cr, unsigned split, std::vector<double> max) {
+    CrossingSet ret;
+    for(unsigned i = 0; i < cr.size(); i++) {
+        Crossings res = reverse_ta(cr[i], max);
+        if(i < split) std::reverse(res.begin(), res.end());
+        ret.push_back(res);
+    }
+    return ret;
+}
+
+CrossingSet reverse_tb(CrossingSet const &cr, unsigned split, std::vector<double> max) {
+    CrossingSet ret;
+    for(unsigned i = 0; i < cr.size(); i++) {
+        Crossings res = reverse_tb(cr[i], split, max);
+        if(i >= split) std::reverse(res.begin(), res.end());
+        ret.push_back(res);
+    }
+    return ret;
+}
+
+void clean(Crossings &cr_a, Crossings &cr_b) {
+/*    if(cr_a.empty()) return;
+    
+    //Remove anything with dupes
+    
+    for(Eraser<Crossings> i(&cr_a); !i.ended(); i++) {
+        const Crossing cur = *i;
+        Eraser<Crossings> next(i);
+        next++;
+        if(near(cur, *next)) {
+            cr_b.erase(std::find(cr_b.begin(), cr_b.end(), cur));
+            for(i = next; near(*i, cur); i++) {
+                cr_b.erase(std::find(cr_b.begin(), cr_b.end(), *i));
+            }
+            continue;
+        }
+    }
+*/
+}
+
+}
diff --git a/src/2geom/crossing.h b/src/2geom/crossing.h
new file mode 100644
index 000000000..72b2eea4b
--- /dev/null
+++ b/src/2geom/crossing.h
@@ -0,0 +1,103 @@
+#ifndef __GEOM_CROSSING_H
+#define __GEOM_CROSSING_H
+
+#include <vector>
+#include <set>
+#include "rect.h"
+#include "sweep.h"
+namespace Geom {
+
+struct Crossing {
+    bool dir; //True: along a, a becomes outside.
+    double ta, tb;  //time on a and b of crossing
+    unsigned a, b;  //storage of indices
+    Crossing() : dir(false), ta(0), tb(1), a(0), b(1) {}
+    Crossing(double t_a, double t_b, bool direction) : dir(direction), ta(t_a), tb(t_b), a(0), b(1) {}
+    Crossing(double t_a, double t_b, unsigned ai, unsigned bi, bool direction) : dir(direction), ta(t_a), tb(t_b), a(ai), b(bi) {}
+    bool operator==(const Crossing & other) const { return a == other.a && b == other.b && dir == other.dir && ta == other.ta && tb == other.tb; }
+    bool operator!=(const Crossing & other) const { return !(*this == other); }
+    unsigned getOther(unsigned cur) { return a == cur ? b : a; }
+};
+
+
+/*inline bool near(Crossing a, Crossing b) {
+    return near(a.ta, b.ta) && near(a.tb, b.tb);
+}
+
+struct NearF { bool operator()(Crossing a, Crossing b) { return near(a, b); } };
+*/
+
+struct CrossingOrder {
+    unsigned ix;
+    CrossingOrder(unsigned i) : ix(i) {}
+    bool operator()(Crossing a, Crossing b) {
+        return (ix == a.a ? a.ta : a.tb) <
+               (ix == b.a ? b.ta : b.tb);
+    }
+};
+
+typedef std::vector<Crossing> Crossings;
+typedef std::vector<Crossings> CrossingSet;
+
+template<typename C>
+std::vector<Rect> bounds(C const &a) {
+    std::vector<Rect> rs;
+    for(unsigned i = 0; i < a.size(); i++) rs.push_back(a[i].boundsFast());
+    return rs;
+}
+
+inline void sort_crossings(Crossings &cr, unsigned ix) { std::sort(cr.begin(), cr.end(), CrossingOrder(ix)); }
+
+template<typename T>
+struct Crosser {
+    virtual Crossings crossings(T const &a, T const &b) { return crossings(std::vector<T>(1,a), std::vector<T>(1,b))[0]; }
+    virtual CrossingSet crossings(std::vector<T> const &a, std::vector<T> const &b) {
+        CrossingSet results(a.size() + b.size(), Crossings());
+    
+        std::vector<std::vector<unsigned> > cull = sweep_bounds(bounds(a), bounds(b));
+        for(unsigned i = 0; i < cull.size(); i++) {
+            for(unsigned jx = 0; jx < cull[i].size(); jx++) {
+                unsigned j = cull[i][jx];
+                unsigned jc = j + a.size();
+                Crossings cr = crossings(a[i], b[j]);
+                for(unsigned k = 0; k < cr.size(); k++) { cr[k].a = i; cr[k].b = jc; }
+                
+                //Sort & add A-sorted crossings
+                sort_crossings(cr, i);
+                Crossings n(results[i].size() + cr.size());
+                std::merge(results[i].begin(), results[i].end(), cr.begin(), cr.end(), n.begin(), CrossingOrder(i));
+                results[i] = n;
+                
+                //Sort & add B-sorted crossings
+                sort_crossings(cr, jc);
+                n.resize(results[jc].size() + cr.size());
+                std::merge(results[jc].begin(), results[jc].end(), cr.begin(), cr.end(), n.begin(), CrossingOrder(jc));
+                results[jc] = n;
+            }
+        }
+        return results;
+    }
+};
+void merge_crossings(Crossings &a, Crossings &b, unsigned i);
+void offset_crossings(Crossings &cr, double a, double b);
+
+Crossings reverse_ta(Crossings const &cr, std::vector<double> max);
+Crossings reverse_tb(Crossings const &cr, unsigned split, std::vector<double> max);
+CrossingSet reverse_ta(CrossingSet const &cr, unsigned split, std::vector<double> max);
+CrossingSet reverse_tb(CrossingSet const &cr, unsigned split, std::vector<double> max);
+
+void clean(Crossings &cr_a, Crossings &cr_b);
+
+}
+
+#endif
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
diff --git a/src/2geom/d2-sbasis.cpp b/src/2geom/d2-sbasis.cpp
new file mode 100644
index 000000000..b68e3e6c5
--- /dev/null
+++ b/src/2geom/d2-sbasis.cpp
@@ -0,0 +1,148 @@
+#include "d2.h"
+/* One would think that we would include d2-sbasis.h, however,
+ * you cannot actually include it in anything - only d2 may import it.
+ * This is due to the trickinesses of template submatching. */
+
+namespace Geom {
+

+SBasis L2(D2<SBasis> const & a, unsigned k) { return sqrt(dot(a, a), k); }

+

+D2<SBasis> multiply(Linear const & a, D2<SBasis> const & b) {

+    return D2<SBasis>(multiply(a, b[X]), multiply(a, b[Y]));

+}

+

+D2<SBasis> multiply(SBasis const & a, D2<SBasis> const & b) {

+    return D2<SBasis>(multiply(a, b[X]), multiply(a, b[Y]));

+}

+

+D2<SBasis> truncate(D2<SBasis> const & a, unsigned terms) {

+    return D2<SBasis>(truncate(a[X], terms), truncate(a[Y], terms));

+}

+

+unsigned sbasis_size(D2<SBasis> const & a) {

+    return std::max((unsigned) a[0].size(), (unsigned) a[1].size());

+}

+

+//TODO: Is this sensical? shouldn't it be like pythagorean or something?

+double tail_error(D2<SBasis> const & a, unsigned tail) {

+    return std::max(a[0].tailError(tail), a[1].tailError(tail));

+}

+

+Piecewise<D2<SBasis> > sectionize(D2<Piecewise<SBasis> > const &a) {

+    Piecewise<SBasis> x = partition(a[0], a[1].cuts), y = partition(a[1], a[0].cuts);

+    assert(x.size() == y.size());

+    Piecewise<D2<SBasis> > ret;

+    for(unsigned i = 0; i < x.size(); i++)

+        ret.push_seg(D2<SBasis>(x[i], y[i]));

+    ret.cuts.insert(ret.cuts.end(), x.cuts.begin(), x.cuts.end());

+    return ret;

+}

+

+D2<Piecewise<SBasis> > make_cuts_independant(Piecewise<D2<SBasis> > const &a) {

+    D2<Piecewise<SBasis> > ret;

+    for(unsigned d = 0; d < 2; d++) {

+        for(unsigned i = 0; i < a.size(); i++)

+            ret[d].push_seg(a[i][d]);

+        ret[d].cuts.insert(ret[d].cuts.end(), a.cuts.begin(), a.cuts.end());

+    }

+    return ret;

+}

+

+Piecewise<D2<SBasis> > rot90(Piecewise<D2<SBasis> > const &M){

+  Piecewise<D2<SBasis> > result;

+  if (M.empty()) return M;

+  result.push_cut(M.cuts[0]);

+  for (unsigned i=0; i<M.size(); i++){

+    result.push(rot90(M[i]),M.cuts[i+1]);

+  }

+  return result;

+}

+

+Piecewise<SBasis> dot(Piecewise<D2<SBasis> > const &a, 

+		      Piecewise<D2<SBasis> > const &b){

+  Piecewise<SBasis > result;

+  if (a.empty() || b.empty()) return result;

+  Piecewise<D2<SBasis> > aa = partition(a,b.cuts);

+  Piecewise<D2<SBasis> > bb = partition(b,a.cuts);

+

+  result.push_cut(aa.cuts.front());

+  for (unsigned i=0; i<aa.size(); i++){

+    result.push(dot(aa.segs[i],bb.segs[i]),aa.cuts[i+1]);

+  }

+  return result;

+}

+

+Piecewise<SBasis> cross(Piecewise<D2<SBasis> > const &a, 

+			Piecewise<D2<SBasis> > const &b){

+  Piecewise<SBasis > result;

+  if (a.empty() || b.empty()) return result;

+  Piecewise<D2<SBasis> > aa = partition(a,b.cuts);

+  Piecewise<D2<SBasis> > bb = partition(b,a.cuts);

+

+  result.push_cut(aa.cuts.front());

+  for (unsigned i=0; i<a.size(); i++){

+    result.push(cross(aa.segs[i],bb.segs[i]),aa.cuts[i+1]);

+  }

+  return result;

+}

+

+/* Replaced by remove_short_cuts in piecewise.h
+//this recursively removes the shortest cut interval until none is shorter than tol.

+//TODO: code this in a more efficient way!

+Piecewise<D2<SBasis> > remove_short_cuts(Piecewise<D2<SBasis> > const &f, double tol){

+    double min = tol;

+    unsigned idx = f.size();

+    for(unsigned i=0; i<f.size(); i++){

+        if (min > f.cuts[i+1]-f.cuts[i]){

+            min = f.cuts[i+1]-f.cuts[i];

+            idx = int(i);

+        }

+    }

+    if (idx==f.size()){

+        return f;

+    }

+    if (f.size()==1) {

+        //removing this seg would result in an empty pw<d2<sb>>...

+        return f;

+    }

+    Piecewise<D2<SBasis> > new_f=f;

+    for (int dim=0; dim<2; dim++){

+        double v = Hat(f.segs.at(idx)[dim][0]);

+        //TODO: what about closed curves?

+        if (idx>0 && f.segs.at(idx-1).at1()==f.segs.at(idx).at0()) 

+            new_f.segs.at(idx-1)[dim][0][1] = v;

+        if (idx<f.size() && f.segs.at(idx+1).at0()==f.segs.at(idx).at1()) 

+            new_f.segs.at(idx+1)[dim][0][0] = v;

+    }

+    double t = (f.cuts.at(idx)+f.cuts.at(idx+1))/2;

+    new_f.cuts.at(idx+1) = t;    

+    

+    new_f.segs.erase(new_f.segs.begin()+idx);

+    new_f.cuts.erase(new_f.cuts.begin()+idx);        

+    return remove_short_cuts(new_f, tol);

+}

+*/
+

+//if tol>0, only force continuity where the jump is smaller than tol.

+Piecewise<D2<SBasis> > force_continuity(Piecewise<D2<SBasis> > const &f, 

+                                        double tol,

+                                        bool closed){

+    if (f.size()==0) return f;

+    Piecewise<D2<SBasis> > result=f;

+    unsigned cur   = (closed)? 0:1;

+    unsigned prev  = (closed)? f.size()-1:0;

+    while(cur<f.size()){

+        Point pt0 = f.segs[prev].at1();

+        Point pt1 = f.segs[cur ].at0();

+        if (tol<=0 || L2sq(pt0-pt1)<tol*tol){

+            pt0 = (pt0+pt1)/2;

+            for (unsigned dim=0; dim<2; dim++){

+                result.segs[prev][dim][0][1]=pt0[dim];

+                result.segs[cur ][dim][0][0]=pt0[dim];

+            }

+        }

+        prev = cur++;

+    }

+    return result;

+}


+}
diff --git a/src/2geom/d2-sbasis.h b/src/2geom/d2-sbasis.h
new file mode 100644
index 000000000..dab2559ca
--- /dev/null
+++ b/src/2geom/d2-sbasis.h
@@ -0,0 +1,88 @@
+#ifdef _2GEOM_D2  /*This is intentional: we don't actually want anyone to
+                    include this, other than D2.h.  If somone else tries, D2
+                    won't be defined.  If it is, this will already be included. */
+#ifndef __2GEOM_SBASIS_CURVE_H
+#define __2GEOM_SBASIS_CURVE_H
+
+#include "sbasis.h"

+#include "sbasis-2d.h"

+#include "piecewise.h"
+
+//TODO: implement intersect
+
+namespace Geom {
+

+inline D2<SBasis> compose(D2<SBasis> const & a, SBasis const & b) {

+    return D2<SBasis>(compose(a[X], b), compose(a[Y], b));

+}

+

+SBasis L2(D2<SBasis> const & a, unsigned k);

+double L2(D2<double> const & a);

+

+D2<SBasis> multiply(Linear const & a, D2<SBasis> const & b);

+inline D2<SBasis> operator*(Linear const & a, D2<SBasis> const & b) { return multiply(a, b); }

+D2<SBasis> multiply(SBasis const & a, D2<SBasis> const & b);

+inline D2<SBasis> operator*(SBasis const & a, D2<SBasis> const & b) { return multiply(a, b); }

+D2<SBasis> truncate(D2<SBasis> const & a, unsigned terms);

+

+unsigned sbasis_size(D2<SBasis> const & a);

+double tail_error(D2<SBasis> const & a, unsigned tail);

+

+//Piecewise<D2<SBasis> > specific decls:

+

+Piecewise<D2<SBasis> > sectionize(D2<Piecewise<SBasis> > const &a);

+D2<Piecewise<SBasis> > make_cuts_independant(Piecewise<D2<SBasis> > const &a);

+Piecewise<D2<SBasis> > rot90(Piecewise<D2<SBasis> > const &a);

+Piecewise<SBasis> dot(Piecewise<D2<SBasis> > const &a, Piecewise<D2<SBasis> > const &b);

+Piecewise<SBasis> cross(Piecewise<D2<SBasis> > const &a, Piecewise<D2<SBasis> > const &b);

+

+Piecewise<D2<SBasis> > force_continuity(Piecewise<D2<SBasis> > const &f, 

+                                        double tol=0,

+                                        bool closed=false);

+

+class CoordIterator

+: public std::iterator<std::input_iterator_tag, SBasis const>

+{

+public:

+  CoordIterator(std::vector<D2<SBasis> >::const_iterator const &iter, unsigned d) : impl_(iter), ix_(d) {}

+

+  inline bool operator==(CoordIterator const &other) { return other.impl_ == impl_; }

+  inline bool operator!=(CoordIterator const &other) { return other.impl_ != impl_; }

+

+  inline SBasis operator*() const {

+        return (*impl_)[ix_];

+  }

+

+  inline CoordIterator &operator++() {

+    ++impl_;

+    return *this;

+  }

+  inline CoordIterator operator++(int) {

+    CoordIterator old=*this;

+    ++(*this);

+    return old;

+  }

+

+private:

+  std::vector<D2<SBasis> >::const_iterator impl_;

+  unsigned ix_;

+};

+

+inline CoordIterator iterateCoord(Piecewise<D2<SBasis> > const &a, unsigned d) {

+    return CoordIterator(a.segs.begin(), d);

+}

+

+//bounds specializations with order

+inline Rect bounds_fast(D2<SBasis> const & s, unsigned order=0) {

+    return Rect(bounds_fast(s[X], order),

+                bounds_fast(s[Y], order));

+}

+inline Rect bounds_local(D2<SBasis> const & s, Interval i, unsigned order=0) {

+    return Rect(bounds_local(s[X], i, order),

+                bounds_local(s[Y], i, order));

+}
+
+}
+
+#endif
+#endif
diff --git a/src/2geom/d2.h b/src/2geom/d2.h
index 82070519c..696ad0191 100644
--- a/src/2geom/d2.h
+++ b/src/2geom/d2.h
@@ -83,6 +83,15 @@ class D2{
     Point valueAt(double t) const {

         boost::function_requires<FragmentConcept<T> >();

         return (*this)(t);

+    }
+    std::vector<Point > valueAndDerivatives(double t, unsigned count) const {
+        std::vector<Coord> x = f[X].valueAndDerivatives(t, count),
+                           y = f[Y].valueAndDerivatives(t, count);
+        std::vector<Point> res;
+        for(unsigned i = 0; i < count; i++) {
+            res.push_back(Point(x[i], y[i]));
+        }
+        return res;
     }

     D2<SBasis> toSBasis() const {

         boost::function_requires<FragmentConcept<T> >();

@@ -91,13 +100,18 @@ class D2{
 

     Point operator()(double t) const;

     Point operator()(double x, double y) const;

-};

-

+};


 template <typename T>

-D2<T> reverse(const D2<T> &a) {

+inline D2<T> reverse(const D2<T> &a) {

     boost::function_requires<FragmentConcept<T> >();

     return D2<T>(reverse(a[X]), reverse(a[Y]));

 }

+
+template <typename T>
+inline D2<T> portion(const D2<T> &a, Coord f, Coord t) {
+    boost::function_requires<FragmentConcept<T> >();
+    return D2<T>(portion(a[X], f, t), portion(a[Y], f, t));
+}
 

 //IMPL: boost::EqualityComparableConcept

 template <typename T>

@@ -347,24 +361,16 @@ D2<T>::operator()(double x, double y) const {
 template<typename T>

 D2<T> derivative(D2<T> const & a) {

     return D2<T>(derivative(a[X]), derivative(a[Y]));

-}

-

+}


 template<typename T>

 D2<T> integral(D2<T> const & a) {

     return D2<T>(integral(a[X]), integral(a[Y]));

-}

+}
 

 } //end namespace Geom

 

-

-

-//TODO: implement intersect

-

-

 #include "rect.h"

-#include "sbasis.h"

-#include "sbasis-2d.h"

-#include "piecewise.h"

+#include "d2-sbasis.h"

 

 namespace Geom{

 

@@ -383,83 +389,7 @@ template <typename T>
 Rect bounds_local(const D2<T> &a, const Interval &t) {

     boost::function_requires<FragmentConcept<T> >();        

     return Rect(bounds_local(a[X], t), bounds_local(a[Y], t));

-}

-

-//D2<SBasis> specific decls:

-

-inline D2<SBasis> compose(D2<SBasis> const & a, SBasis const & b) {

-    return D2<SBasis>(compose(a[X], b), compose(a[Y], b));

-}

-

-SBasis L2(D2<SBasis> const & a, unsigned k);

-double L2(D2<double> const & a);

-

-inline D2<SBasis> portion(D2<SBasis> const &M, double t0, double t1){

-    return D2<SBasis>(portion(M[0],t0,t1),portion(M[1],t0,t1));

-}

-

-D2<SBasis> multiply(Linear const & a, D2<SBasis> const & b);

-inline D2<SBasis> operator*(Linear const & a, D2<SBasis> const & b) { return multiply(a, b); }

-D2<SBasis> multiply(SBasis const & a, D2<SBasis> const & b);

-inline D2<SBasis> operator*(SBasis const & a, D2<SBasis> const & b) { return multiply(a, b); }

-D2<SBasis> truncate(D2<SBasis> const & a, unsigned terms);

-

-unsigned sbasis_size(D2<SBasis> const & a);

-double tail_error(D2<SBasis> const & a, unsigned tail);

-

-//Piecewise<D2<SBasis> > specific decls:

-

-Piecewise<D2<SBasis> > sectionize(D2<Piecewise<SBasis> > const &a);

-D2<Piecewise<SBasis> > make_cuts_independant(Piecewise<D2<SBasis> > const &a);

-Piecewise<D2<SBasis> > rot90(Piecewise<D2<SBasis> > const &a);

-Piecewise<SBasis> dot(Piecewise<D2<SBasis> > const &a, Piecewise<D2<SBasis> > const &b);

-Piecewise<SBasis> cross(Piecewise<D2<SBasis> > const &a, Piecewise<D2<SBasis> > const &b);

-

-Piecewise<D2<SBasis> > force_continuity(Piecewise<D2<SBasis> > const &f, 

-                                        double tol=0,

-                                        bool closed=false);

-

-class CoordIterator

-: public std::iterator<std::input_iterator_tag, SBasis const>

-{

-public:

-  CoordIterator(std::vector<D2<SBasis> >::const_iterator const &iter, unsigned d) : impl_(iter), ix_(d) {}

-

-  inline bool operator==(CoordIterator const &other) { return other.impl_ == impl_; }

-  inline bool operator!=(CoordIterator const &other) { return other.impl_ != impl_; }

-

-  inline SBasis operator*() const {

-        return (*impl_)[ix_];

-  }

-

-  inline CoordIterator &operator++() {

-    ++impl_;

-    return *this;

-  }

-  inline CoordIterator operator++(int) {

-    CoordIterator old=*this;

-    ++(*this);

-    return old;

-  }

-

-private:

-  std::vector<D2<SBasis> >::const_iterator impl_;

-  unsigned ix_;

-};

-

-inline CoordIterator iterateCoord(Piecewise<D2<SBasis> > const &a, unsigned d) {

-    return CoordIterator(a.segs.begin(), d);

-}

-

-//bounds specializations with order

-inline Rect bounds_fast(D2<SBasis> const & s, unsigned order=0) {

-    return Rect(bounds_fast(s[X], order),

-                bounds_fast(s[Y], order));

-}

-inline Rect bounds_local(D2<SBasis> const & s, Interval i, unsigned order=0) {

-    return Rect(bounds_local(s[X], i, order),

-                bounds_local(s[Y], i, order));

-}

+}


 };

 

 /*

diff --git a/src/2geom/interval.h b/src/2geom/interval.h
index 459f2cd49..19e08978c 100644
--- a/src/2geom/interval.h
+++ b/src/2geom/interval.h
@@ -59,57 +59,50 @@ public:
             _b[0] = v; _b[1] = u;
         }
     }
-
+    
     double operator[](unsigned i) const {
         assert(i < 2);
         return _b[i];
     }
-    double& operator[](unsigned i) { return _b[i]; }  //Trust the user...
-
-    Coord min() const { return _b[0]; }
-    Coord max() const { return _b[1]; }
-    Coord extent() const { return _b[1] - _b[0]; }
-    Coord middle() const { return (_b[1] + _b[0]) * 0.5; }
-
-    bool isEmpty() const { return _b[0] == _b[1]; }
-    bool contains(Coord val) const { return _b[0] <= val && val <= _b[1]; }
+    inline double& operator[](unsigned i) { return _b[i]; }  //Trust the user...
+    
+    inline Coord min() const { return _b[0]; }
+    inline Coord max() const { return _b[1]; }
+    inline Coord extent() const { return _b[1] - _b[0]; }
+    inline Coord middle() const { return (_b[1] + _b[0]) * 0.5; }
+    
+    inline bool isEmpty() const { return _b[0] == _b[1]; }
+    inline bool contains(Coord val) const { return _b[0] <= val && val <= _b[1]; }
     bool contains(const Interval & val) const { return _b[0] <= val._b[0] && val._b[1] <= _b[1]; }
     bool intersects(const Interval & val) const {
         return contains(val._b[0]) || contains(val._b[1]) || val.contains(*this);
     }
-
-    static Interval fromArray(const Coord* c, int n) {
-        assert(n > 0);
-        Interval result(c[0]);
-        for(int i = 1; i < n; i++) result.extendTo(c[i]);
-        return result;
-    }
-
-    bool operator==(Interval other) { return _b[0] == other._b[0] && _b[1] == other._b[1]; }
-    bool operator!=(Interval other) { return _b[0] != other._b[0] || _b[1] != other._b[1]; }
-
+    
+    inline bool operator==(Interval other) { return _b[0] == other._b[0] && _b[1] == other._b[1]; }
+    inline bool operator!=(Interval other) { return _b[0] != other._b[0] || _b[1] != other._b[1]; }
+    
     //IMPL: OffsetableConcept
     //TODO: rename output_type to something else in the concept
     typedef Coord output_type;
-    Interval operator+(Coord amnt) {
+    inline Interval operator+(Coord amnt) {
         return Interval(_b[0] + amnt, _b[1] + amnt);
     }
-    Interval operator-(Coord amnt) {
+    inline Interval operator-(Coord amnt) {
         return Interval(_b[0] - amnt, _b[1] - amnt);
     }
-    Interval operator+=(Coord amnt) {
+    inline Interval operator+=(Coord amnt) {
         _b[0] += amnt; _b[1] += amnt;
         return *this;
     }
-    Interval operator-=(Coord amnt) {
+    inline Interval operator-=(Coord amnt) {
         _b[0] -= amnt; _b[1] -= amnt;
         return *this;
     }
-
+    
     //IMPL: ScalableConcept
-    Interval operator-() const { return Interval(*this); }
-    Interval operator*(Coord s) const { return Interval(_b[0]*s, _b[1]*s); }
-    Interval operator/(Coord s) const { return Interval(_b[0]/s, _b[1]/s); }
+    inline Interval operator-() const { return Interval(*this); }
+    inline Interval operator*(Coord s) const { return Interval(_b[0]*s, _b[1]*s); }
+    inline Interval operator/(Coord s) const { return Interval(_b[0]/s, _b[1]/s); }
     Interval operator*=(Coord s) {
         if(s < 0) {
             Coord temp = _b[0];
@@ -133,7 +126,7 @@ public:
         }
         return *this;
     }
-
+    
     //TODO: NaN handleage for the next two?
     //TODO: Evaluate if wrap behaviour is proper.
     //If val > max, then rather than becoming a min==max range, it 'wraps' over
@@ -154,18 +147,25 @@ public:
             _b[1] = val;
         }
     }
-
-    void extendTo(Coord val) {
+    
+    inline void extendTo(Coord val) {
        if(val < _b[0]) _b[0] = val;
        if(val > _b[1]) _b[1] = val;  //no else, as we want to handle NaN
     }
-
-    void expandBy(double amnt) {
+    
+    static Interval fromArray(const Coord* c, int n) {
+        assert(n > 0);
+        Interval result(c[0]);
+        for(int i = 1; i < n; i++) result.extendTo(c[i]);
+        return result;
+    }
+    
+    inline void expandBy(double amnt) {
         _b[0] -= amnt;
         _b[1] += amnt;
     }
-
-    void unionWith(const Interval & a) {
+    
+    inline void unionWith(const Interval & a) {
         if(a._b[0] < _b[0]) _b[0] = a._b[0];
         if(a._b[1] > _b[1]) _b[1] = a._b[1];
     }
diff --git a/src/2geom/matrix.cpp b/src/2geom/matrix.cpp
index f53795943..c0e64ad33 100644
--- a/src/2geom/matrix.cpp
+++ b/src/2geom/matrix.cpp
@@ -94,7 +94,7 @@ void Matrix::setExpansionY(double val) {
     double exp_y = expansionY();
     if(!near(exp_y, 0.0)) {  //TODO: best way to deal with it is to skip op?
         double coef = val / expansionY();
-        for(unsigned i=2;i<4;i++) _c[i] *= coef;
+        for(unsigned i=2; i<4; i++) _c[i] *= coef;
     }
 }
 
@@ -105,6 +105,8 @@ void Matrix::setIdentity() {
     _c[4] = 0.0; _c[5] = 0.0;
 }
 
+//TODO: use eps
+
 bool Matrix::isIdentity(Coord const eps) const {
     return near(_c[0], 1.0) && near(_c[1], 0.0) &&
            near(_c[2], 0.0) && near(_c[3], 1.0) &&
@@ -151,6 +153,14 @@ bool Matrix::isRotation(Coord const eps) const {
            near(_c[0]*_c[0] + _c[1]*_c[1], 1.0);
 }
 
+bool Matrix::onlyScaleAndTranslation(Coord const eps) const {
+    return near(_c[0], _c[3]) && near(_c[1], 0) && near(_c[2], 0);
+}
+
+bool Matrix::flips() const {
+    return cross(xAxis(), yAxis()) > 0;
+}
+
 /** Returns the Scale/Rotate/skew part of the matrix without the translation part. */
 Matrix Matrix::without_translation() const {
     return Matrix(_c[0], _c[1], _c[2], _c[3], 0, 0);
diff --git a/src/2geom/matrix.h b/src/2geom/matrix.h
index 6a45b50c8..c39c99716 100644
--- a/src/2geom/matrix.h
+++ b/src/2geom/matrix.h
@@ -90,6 +90,9 @@ class Matrix {
     bool isRotation(double eps = EPSILON) const;
     bool isScale(double eps = EPSILON) const;
     bool isUniformScale(double eps = EPSILON) const;
+    bool onlyScaleAndTranslation(double eps = EPSILON) const;
+
+    bool flips() const;
 
     Matrix without_translation() const;
 
diff --git a/src/2geom/ord.h b/src/2geom/ord.h
new file mode 100644
index 000000000..d0e348aec
--- /dev/null
+++ b/src/2geom/ord.h
@@ -0,0 +1,37 @@
+
+#ifndef __2GEOM_ORD__
+#define __2GEOM_ORD__
+
+namespace {
+
+enum Cmp {

+  LESS_THAN=-1,

+  GREATER_THAN=1,

+  EQUAL_TO=0

+};

+
+inline Cmp operator-(Cmp x) {
+  switch(x) {
+  case LESS_THAN:
+    return GREATER_THAN;
+  case GREATER_THAN:
+    return LESS_THAN;
+  case EQUAL_TO:
+    return EQUAL_TO;
+  }
+}
+

+template <typename T1, typename T2>

+inline Cmp cmp(T1 const &a, T2 const &b) {

+  if ( a < b ) {

+    return LESS_THAN;

+  } else if ( b < a ) {

+    return GREATER_THAN;

+  } else {

+    return EQUAL_TO;

+  }

+}
+
+}
+
+#endif
diff --git a/src/2geom/path-intersection.cpp b/src/2geom/path-intersection.cpp
new file mode 100644
index 000000000..d641fcc08
--- /dev/null
+++ b/src/2geom/path-intersection.cpp
@@ -0,0 +1,588 @@
+#include "path-intersection.h"
+
+#include "ord.h"
+
+//for path_direction:
+#include "sbasis-geometric.h"
+
+namespace Geom {
+
+/* This function computes the winding of the path, given a reference point.
+ * Positive values correspond to counter-clockwise in the mathematical coordinate system,
+ * and clockwise in screen coordinates.  This particular implementation casts a ray in
+ * the positive x direction.  It iterates the path, checking for intersection with the
+ * bounding boxes.  If an intersection is found, the initial/final Y value of the curve is
+ * used to derive a delta on the winding value.  If the point is within the bounding box,
+ * the curve specific winding function is called.
+ */
+int winding(Path const &path, Point p) {
+  //start on a segment which is not a horizontal line with y = p[y]
+  Path::const_iterator start;
+  for(Path::const_iterator iter = path.begin(); ; ++iter) {
+    if(iter == path.end_closed()) { return 0; }
+    if(iter->initialPoint()[Y]!=p[Y])  { start = iter; break; }
+    if(iter->finalPoint()[Y]!=p[Y])    { start = iter; break; }
+    if(iter->boundsFast().height()!=0.){ start = iter; break; }
+  }
+  int wind = 0;
+  int cnt = 0;
+  bool starting = true;
+  for (Path::const_iterator iter = start; iter != start || starting
+       ; ++iter, iter = (iter == path.end_closed()) ? path.begin() : iter )
+  {
+    cnt++;
+    if(cnt > path.size()) return wind;  //some bug makes this required
+    starting = false;
+    Rect bounds = iter->boundsFast();
+    Coord x = p[X], y = p[Y];
+    
+    if(x > bounds.right() || !bounds[Y].contains(y)) continue; //ray doesn't intersect box
+    
+    Point final = iter->finalPoint();
+    Point initial = iter->initialPoint();
+    Cmp final_to_ray = cmp(final[Y], y);
+    Cmp initial_to_ray = cmp(initial[Y], y);
+    
+    // if y is included, these will have opposite values, giving order.
+    Cmp c = cmp(final_to_ray, initial_to_ray); 
+    if(x < bounds.left()) {
+        // ray goes through bbox
+        // winding delta determined by position of endpoints
+        if(final_to_ray != EQUAL_TO) {
+            wind += int(c); // GT = counter-clockwise = 1; LT = clockwise = -1; EQ = not-included = 0
+            //std::cout << int(c) << " ";
+            goto cont;
+        }
+    } else {
+        //inside bbox, use custom per-curve winding thingie
+        int delt = iter->winding(p);
+        wind += delt;
+        //std::cout << "n" << delt << " ";
+    }
+    //Handling the special case of an endpoint on the ray:
+    if(final[Y] == y) {
+        //Traverse segments until it breaks away from y
+        //99.9% of the time this will happen the first go
+        Path::const_iterator next = iter;
+        next++;
+        for(; ; next++) {
+            if(next == path.end_closed()) next = path.begin();
+            Rect bnds = next->boundsFast();
+            //TODO: X considerations
+            if(bnds.height() > 0) {
+                //It has diverged
+                if(bnds.contains(p)) {
+                    const double fudge = 0.01;
+                    if(cmp(y, next->valueAt(fudge, Y)) == initial_to_ray) {
+                        wind += int(c);
+                        std::cout << "!!!!!" << int(c) << " ";
+                    }
+                    iter = next; // No increment, as the rest of the thing hasn't been counted.
+                } else {
+                    Coord ny = next->initialPoint()[Y];
+                    if(cmp(y, ny) == initial_to_ray) {
+                        //Is a continuation through the ray, so counts windingwise
+                        wind += int(c);
+                        std::cout << "!!!!!" << int(c) << " ";
+                    }
+                    iter = ++next;
+                }
+                goto cont;
+            }
+            if(next==start) return wind;
+        }
+        //Looks like it looped, which means everything's flat
+        return 0;
+    }
+    
+    cont:(void)0;
+  }
+  return wind;
+}
+
+/* This function should only be applied to simple paths (regions), as otherwise
+ * a boolean winding direction is undefined.  It returns true for fill, false for
+ * hole.  Defaults to using the sign of area when it reaches funny cases.
+ */
+bool path_direction(Path const &p) {
+    //could probably be more efficient, but this is a quick job
+    double y = p.initialPoint()[Y];
+    double x = p.initialPoint()[X];
+    Cmp res = cmp(p[0].finalPoint()[Y], y);
+    goto doh;
+    for(unsigned i = 1; i <= p.size(); i++) {
+        Cmp final_to_ray = cmp(p[i].finalPoint()[Y], y);
+        Cmp initial_to_ray = cmp(p[i].initialPoint()[Y], y);
+        // if y is included, these will have opposite values, giving order.
+        Cmp c = cmp(final_to_ray, initial_to_ray);
+        if(c != EQUAL_TO) {
+            std::vector<double> rs = p[i].roots(y, Y);
+            for(unsigned j = 0; j < rs.size(); j++) {
+                double nx = p[i].valueAt(rs[j], X);
+                if(nx > x) {
+                    x = nx;
+                    res = c;
+                }
+            }
+        } else if(final_to_ray == EQUAL_TO) goto doh;
+    }
+    return res < 0;
+    
+    doh:
+        //Otherwise fallback on area
+        
+        Piecewise<D2<SBasis> > pw = p.toPwSb();
+        double area;
+        Point centre;
+        Geom::centroid(pw, centre, area);
+        return area > 0;
+}
+
+//pair intersect code based on njh's pair-intersect
+
+// A little sugar for appending a list to another
+template<typename T>
+void append(T &a, T const &b) {
+    a.insert(a.end(), b.begin(), b.end());
+}
+
+/* Finds the intersection between the lines defined by A0 & A1, and B0 & B1.
+ * Returns through the last 3 parameters, returning the t-values on the lines
+ * and the cross-product of the deltas (a useful byproduct).  The return value
+ * indicates if the time values are within their proper range on the line segments.
+ */
+bool
+linear_intersect(Point A0, Point A1, Point B0, Point B1,
+                 double &tA, double &tB, double &det) {
+    // kramers rule as cross products
+    Point Ad = A1 - A0,
+          Bd = B1 - B0,
+           d = B0 - A0;
+    det = cross(Ad, Bd);
+    if( 1.0 + det == 1.0 )
+        return false;
+    else
+    {
+        double detinv = 1.0 / det;
+        tA = cross(d, Bd) * detinv;
+        tB = cross(d, Ad) * detinv;
+        return tA >= 0. && tA <= 1. && tB >= 0. && tB <= 1.;
+    }
+}
+
+/* This uses the local bounds functions of curves to generically intersect two.
+ * It passes in the curves, time intervals, and keeps track of depth, while
+ * returning the results through the Crossings parameter.
+ */
+void pair_intersect(Curve const & A, double Al, double Ah, 
+                    Curve const & B, double Bl, double Bh,
+                    Crossings &ret,  unsigned depth=0) {
+   // std::cout << depth << "(" << Al << ", " << Ah << ")\n";
+    Rect Ar = A.boundsLocal(Interval(Al, Ah));
+    if(Ar.isEmpty()) return;
+
+    Rect Br = B.boundsLocal(Interval(Bl, Bh));
+    if(Br.isEmpty()) return;
+    
+    if(!Ar.intersects(Br)) return;
+    
+    //Checks the general linearity of the function
+    if((depth > 12) || (A.boundsLocal(Interval(Al, Ah), 1).maxExtent() < 0.1 
+                    &&  B.boundsLocal(Interval(Bl, Bh), 1).maxExtent() < 0.1)) {
+        double tA, tB, c;
+        if(linear_intersect(A.pointAt(Al), A.pointAt(Ah), 
+                            B.pointAt(Bl), B.pointAt(Bh), 
+                            tA, tB, c)) {
+            tA = tA * (Ah - Al) + Al;
+            tB = tB * (Bh - Bl) + Bl;
+            if(depth % 2)
+                ret.push_back(Crossing(tB, tA, c < 0));
+            else
+                ret.push_back(Crossing(tA, tB, c > 0));
+            return;
+        }
+    }
+    if(depth > 12) return;
+    double mid = (Bl + Bh)/2;
+    pair_intersect(B, Bl, mid,
+                    A, Al, Ah,
+                    ret, depth+1);
+    pair_intersect(B, mid, Bh,
+                    A, Al, Ah,
+                    ret, depth+1);
+}
+
+// A simple wrapper around pair_intersect
+Crossings SimpleCrosser::crossings(Curve const &a, Curve const &b) {
+    Crossings ret;
+    pair_intersect(a, 0, 1, b, 0, 1, ret);
+    return ret;
+}
+
+/* Takes two paths and time ranges on them, with the invariant that the
+ * paths are monotonic on the range.  Splits A when the linear intersection
+ * doesn't exist or is inaccurate.  Uses the fact that it is monotonic to
+ * do very fast local bounds.
+ */
+void mono_pair(Path const &A, double Al, double Ah,
+               Path const &B, double Bl, double Bh,
+               Crossings &ret, double tol, unsigned depth = 0) {
+    if( Al >= Ah || Bl >= Bh) return;
+    std::cout << " " << depth << "[" << Al << ", " << Ah << "]" << "[" << Bl << ", " << Bh << "]";
+
+    Point A0 = A.pointAt(Al), A1 = A.pointAt(Ah),
+          B0 = B.pointAt(Bl), B1 = B.pointAt(Bh);
+    //inline code that this implies? (without rect/interval construction)
+    if(!Rect(A0, A1).intersects(Rect(B0, B1)) || A0 == A1 || B0 == B1) return;
+    
+    //Checks the general linearity of the function
+    //if((depth > 12) || (A.boundsLocal(Interval(Al, Ah), 1).maxExtent() < 0.1 
+    //                &&  B.boundsLocal(Interval(Bl, Bh), 1).maxExtent() < 0.1)) {
+        double tA, tB, c;
+        if(linear_intersect(A0, A1, B0, B1,
+                            tA, tB, c)) {
+            tA = tA * (Ah - Al) + Al;
+            tB = tB * (Bh - Bl) + Bl;
+            if(depth % 2)
+                ret.push_back(Crossing(tB, tA, c < 0));
+            else
+                ret.push_back(Crossing(tA, tB, c > 0));
+            return;
+        }
+    //}
+    if(depth > 12) return;
+    double mid = (Bl + Bh)/2;
+    mono_pair(B, Bl, mid,
+              A, Al, Ah,
+              ret, depth+1);
+    mono_pair(B, mid, Bh,
+              A, Al, Ah,
+              ret, depth+1);
+}
+
+// This returns the times when the x or y derivative is 0 in the curve.
+std::vector<double> curve_mono_splits(Curve const &d) {
+    std::vector<double> rs = d.roots(0, X);
+    append(rs, d.roots(0, Y));
+    std::sort(rs.begin(), rs.end());
+    return rs;
+}
+
+// Convenience function to add a value to each entry in a vector of doubles.
+std::vector<double> offset_doubles(std::vector<double> const &x, double offs) {
+    std::vector<double> ret;
+    for(unsigned i = 0; i < x.size(); i++) {
+        ret.push_back(x[i] + offs);
+    }
+    return ret;
+}
+
+/* Finds all the monotonic splits for a path.  Only includes the split between
+ * curves if they switch derivative directions at that point.
+ */
+std::vector<double> path_mono_splits(Path const &p) {
+    std::vector<double> ret;
+    if(p.empty()) return ret;
+    ret.push_back(0);
+    
+    Curve* deriv = p[0].derivative();
+    append(ret, curve_mono_splits(*deriv));
+    delete deriv;
+    
+    bool pdx=2, pdy=2;  //Previous derivative direction
+    for(unsigned i = 0; i <= p.size(); i++) {
+        deriv = p[i].derivative();
+        std::vector<double> spl = offset_doubles(curve_mono_splits(*deriv), i);
+        delete deriv;
+        bool dx = p[i].initialPoint()[X] > (spl.empty()? p[i].finalPoint()[X] :
+                                                         p.valueAt(spl.front(), X));
+        bool dy = p[i].initialPoint()[Y] > (spl.empty()? p[i].finalPoint()[Y] :
+                                                         p.valueAt(spl.front(), Y));
+        //The direction changed, include the split time
+        if(dx != pdx || dy != pdy) {
+            ret.push_back(i);
+            pdx = dx; pdy = dy;
+        }
+        append(ret, spl);
+    }
+    return ret;
+}
+
+/* Applies path_mono_splits to multiple paths, and returns the results such that 
+ * time-set i corresponds to Path i.
+ */
+std::vector<std::vector<double> > paths_mono_splits(std::vector<Path> const &ps) {
+    std::vector<std::vector<double> > ret;
+    for(unsigned i = 0; i < ps.size(); i++)
+        ret.push_back(path_mono_splits(ps[i]));
+    return ret;
+}
+
+/* Processes the bounds for a list of paths and a list of splits on them, yielding a list of rects for each.
+ * Each entry i corresponds to path i of the input.  The number of rects in each entry is guaranteed to be the
+ * number of splits for that path, subtracted by one.
+ */
+std::vector<std::vector<Rect> > split_bounds(std::vector<Path> const &p, std::vector<std::vector<double> > splits) {
+    std::vector<std::vector<Rect> > ret;
+    for(unsigned i = 0; i < p.size(); i++) {
+        std::vector<Rect> res;
+        for(unsigned j = 1; j < splits[i].size(); j++)
+            res.push_back(Rect(p[i].pointAt(splits[i][j-1]), p[i].pointAt(splits[i][j])));
+        ret.push_back(res);
+    }
+    return ret;
+}
+
+/* This is the main routine of "MonoCrosser", and implements a monotonic strategy on multiple curves.
+ * Finds crossings between two sets of paths, yielding a CrossingSet.  [0, a.size()) of the return correspond
+ * to the sorted crossings of a with paths of b.  The rest of the return, [a.size(), a.size() + b.size()],
+ * corresponds to the sorted crossings of b with paths of a.
+ *
+ * This function does two sweeps, one on the bounds of each path, and after that cull, one on the curves within.
+ * This leads to a certain amount of code complexity, however, most of that is factored into the above functions
+ */
+CrossingSet MonoCrosser::crossings(std::vector<Path> const &a, std::vector<Path> const &b) {
+    if(b.empty()) return CrossingSet(a.size(), Crossings());
+    CrossingSet results(a.size() + b.size(), Crossings());
+    if(a.empty()) return results;
+    
+    std::vector<std::vector<double> > splits_a = paths_mono_splits(a), splits_b = paths_mono_splits(b);
+    std::vector<std::vector<Rect> > bounds_a = split_bounds(a, splits_a), bounds_b = split_bounds(b, splits_b);
+    
+    std::vector<Rect> bounds_a_union, bounds_b_union; 
+    for(unsigned i = 0; i < bounds_a.size(); i++) bounds_a_union.push_back(union_list(bounds_a[i]));
+    for(unsigned i = 0; i < bounds_b.size(); i++) bounds_b_union.push_back(union_list(bounds_b[i]));
+    
+    std::vector<std::vector<unsigned> > cull = sweep_bounds(bounds_a_union, bounds_b_union);
+    Crossings n;
+    for(unsigned i = 0; i < cull.size(); i++) {
+        for(unsigned jx = 0; jx < cull[i].size(); jx++) {
+            unsigned j = cull[i][jx];
+            unsigned jc = j + a.size();
+            Crossings res;
+            
+            //Sweep of the monotonic portions
+            std::vector<std::vector<unsigned> > cull2 = sweep_bounds(bounds_a[i], bounds_b[j]);
+            for(unsigned k = 0; k < cull2.size(); k++) {
+                for(unsigned lx = 0; lx < cull2[k].size(); lx++) {
+                    unsigned l = cull2[k][lx];
+                    mono_pair(a[i], splits_a[i][k-1], splits_a[i][k],
+                              b[j], splits_b[j][l-1], splits_b[j][l],
+                              res, .1);
+                }
+            }
+            
+            for(unsigned k = 0; k < res.size(); k++) { res[k].a = i; res[k].b = jc; }
+            
+            merge_crossings(results[i], res, i);
+            merge_crossings(results[i], res, jc);
+        }
+    }
+
+    return results;
+}
+
+/* This function is similar codewise to the MonoCrosser, the main difference is that it deals with
+ * only one set of paths and includes self intersection
+CrossingSet crossings_among(std::vector<Path> const &p) {
+    CrossingSet results(p.size(), Crossings());
+    if(p.empty()) return results;
+    
+    std::vector<std::vector<double> > splits = paths_mono_splits(p);
+    std::vector<std::vector<Rect> > prs = split_bounds(p, splits);
+    std::vector<Rect> rs;
+    for(unsigned i = 0; i < prs.size(); i++) rs.push_back(union_list(prs[i]));
+    
+    std::vector<std::vector<unsigned> > cull = sweep_bounds(rs);
+    
+    //we actually want to do the self-intersections, so add em in:
+    for(unsigned i = 0; i < cull.size(); i++) cull[i].push_back(i);
+    
+    for(unsigned i = 0; i < cull.size(); i++) {
+        for(unsigned jx = 0; jx < cull[i].size(); jx++) {
+            unsigned j = cull[i][jx];
+            Crossings res;
+            
+            //Sweep of the monotonic portions
+            std::vector<std::vector<unsigned> > cull2 = sweep_bounds(prs[i], prs[j]);
+            for(unsigned k = 0; k < cull2.size(); k++) {
+                for(unsigned lx = 0; lx < cull2[k].size(); lx++) {
+                    unsigned l = cull2[k][lx];
+                    mono_pair(p[i], splits[i][k-1], splits[i][k],
+                              p[j], splits[j][l-1], splits[j][l],
+                              res, .1);
+                }
+            }
+            
+            for(unsigned k = 0; k < res.size(); k++) { res[k].a = i; res[k].b = j; }
+            
+            merge_crossings(results[i], res, i);
+            merge_crossings(results[j], res, j);
+        }
+    }
+    
+    return results;
+}
+*/
+
+Crossings curve_self_crossings(Curve const &a) {
+    Crossings res;
+    std::vector<double> spl;
+    spl.push_back(0);
+    append(spl, curve_mono_splits(a));
+    spl.push_back(1);
+    for(unsigned i = 1; i < spl.size(); i++)
+        for(unsigned j = i+1; j < spl.size(); j++)
+            pair_intersect(a, spl[i-1], spl[i], a, spl[j-1], spl[j], res);
+    return res;
+}
+
+/*
+void mono_curve_intersect(Curve const & A, double Al, double Ah, 
+                          Curve const & B, double Bl, double Bh,
+                          Crossings &ret,  unsigned depth=0) {
+   // std::cout << depth << "(" << Al << ", " << Ah << ")\n";
+    Point A0 = A.pointAt(Al), A1 = A.pointAt(Ah),
+          B0 = B.pointAt(Bl), B1 = B.pointAt(Bh);
+    //inline code that this implies? (without rect/interval construction)
+    if(!Rect(A0, A1).intersects(Rect(B0, B1)) || A0 == A1 || B0 == B1) return;
+     
+    //Checks the general linearity of the function
+    if((depth > 12) || (A.boundsLocal(Interval(Al, Ah), 1).maxExtent() < 0.1 
+                    &&  B.boundsLocal(Interval(Bl, Bh), 1).maxExtent() < 0.1)) {
+        double tA, tB, c;
+        if(linear_intersect(A0, A1, B0, B1, tA, tB, c)) {
+            tA = tA * (Ah - Al) + Al;
+            tB = tB * (Bh - Bl) + Bl;
+            if(depth % 2)
+                ret.push_back(Crossing(tB, tA, c < 0));
+            else
+                ret.push_back(Crossing(tA, tB, c > 0));
+            return;
+        }
+    }
+    if(depth > 12) return;
+    double mid = (Bl + Bh)/2;
+    mono_curve_intersect(B, Bl, mid,
+                         A, Al, Ah,
+                         ret, depth+1);
+    mono_curve_intersect(B, mid, Bh,
+                         A, Al, Ah,
+                         ret, depth+1);
+}
+
+std::vector<std::vector<double> > curves_mono_splits(Path const &p) {
+    std::vector<std::vector<double> > ret;
+    for(unsigned i = 0; i <= p.size(); i++) {
+        std::vector<double> spl;
+        spl.push_back(0);
+        append(spl, curve_mono_splits(p[i]));
+        spl.push_back(1);
+        ret.push_back(spl);
+    }
+}
+
+std::vector<std::vector<Rect> > curves_split_bounds(Path const &p, std::vector<std::vector<double> > splits) {
+    std::vector<std::vector<Rect> > ret;
+    for(unsigned i = 0; i < splits.size(); i++) {
+        std::vector<Rect> res;
+        for(unsigned j = 1; j < splits[i].size(); j++)
+            res.push_back(Rect(p.pointAt(splits[i][j-1]+i), p.pointAt(splits[i][j]+i)));
+        ret.push_back(res);
+    }
+    return ret;
+}
+
+
+Crossings path_self_crossings(Path const &p) {
+    Crossings ret;
+    std::vector<std::vector<unsigned> > cull = sweep_bounds(bounds(p));
+    std::vector<std::vector<double> > spl = curves_mono_splits(p);
+    std::vector<std::vector<Rect> > bnds = curves_split_bounds(p, spl);
+    for(unsigned i = 0; i < cull.size(); i++) {
+        Crossings res;
+        for(unsigned k = 1; k < spl[i].size(); k++)
+            for(unsigned l = k+1; l < spl[i].size(); l++)
+                mono_curve_intersect(p[i], spl[i][k-1], spl[i][k], p[i], spl[i][l-1], spl[i][l], res);
+        offset_crossings(res, i, i);
+        append(ret, res);
+        for(unsigned jx = 0; jx < cull[i].size(); jx++) {
+            unsigned j = cull[i][jx];
+            res.clear();
+            
+            std::vector<std::vector<unsigned> > cull2 = sweep_bounds(bnds[i], bnds[j]);
+            for(unsigned k = 0; k < cull2.size(); k++) {
+                for(unsigned lx = 0; lx < cull2[k].size(); lx++) {
+                    unsigned l = cull2[k][lx];
+                    mono_curve_intersect(p[i], spl[i][k-1], spl[i][k], p[j], spl[j][l-1], spl[j][l], res);
+                }
+            }
+            
+            //if(fabs(int(i)-j) == 1 || fabs(int(i)-j) == p.size()-1) {
+                Crossings res2;
+                for(unsigned k = 0; k < res.size(); k++) {
+                    if(res[k].ta != 0 && res[k].ta != 1 && res[k].tb != 0 && res[k].tb != 1) {
+                        res.push_back(res[k]);
+                    }
+                }
+                res = res2;
+            //}
+            offset_crossings(res, i, j);
+            append(ret, res);
+        }
+    }
+    return ret;
+}
+*/
+
+Crossings path_self_crossings(Path const &p) {
+    Crossings ret;
+    std::vector<std::vector<unsigned> > cull = sweep_bounds(bounds(p));
+    for(unsigned i = 0; i < cull.size(); i++) {
+        Crossings res = curve_self_crossings(p[i]);
+        offset_crossings(res, i, i);
+        append(ret, res);
+        for(unsigned jx = 0; jx < cull[i].size(); jx++) {
+            unsigned j = cull[i][jx];
+            res.clear();
+            pair_intersect(p[i], 0, 1, p[j], 0, 1, res);
+            
+            //if(fabs(int(i)-j) == 1 || fabs(int(i)-j) == p.size()-1) {
+                Crossings res2;
+                for(unsigned k = 0; k < res.size(); k++) {
+                    if(res[k].ta != 0 && res[k].ta != 1 && res[k].tb != 0 && res[k].tb != 1) {
+                        res.push_back(res[k]);
+                    }
+                }
+                res = res2;
+            //}
+            offset_crossings(res, i, j);
+            append(ret, res);
+        }
+    }
+    return ret;
+}
+
+CrossingSet crossings_among(std::vector<Path> const &p) {
+    CrossingSet results(p.size(), Crossings());
+    if(p.empty()) return results;
+    
+    SimpleCrosser cc;
+    
+    std::vector<std::vector<unsigned> > cull = sweep_bounds(bounds(p));
+    for(unsigned i = 0; i < cull.size(); i++) {
+        Crossings res = path_self_crossings(p[i]);
+        for(unsigned k = 0; k < res.size(); k++) { res[k].a = res[k].b = i; }
+        merge_crossings(results[i], res, i);
+        for(unsigned jx = 0; jx < cull[i].size(); jx++) {
+            unsigned j = cull[i][jx];
+            
+            Crossings res = cc.crossings(p[i], p[j]);
+            for(unsigned k = 0; k < res.size(); k++) { res[k].a = i; res[k].b = j; }
+            merge_crossings(results[i], res, i);
+            merge_crossings(results[j], res, j);
+        }
+    }
+}
+
+}
diff --git a/src/2geom/path-intersection.h b/src/2geom/path-intersection.h
new file mode 100644
index 000000000..3401386e0
--- /dev/null
+++ b/src/2geom/path-intersection.h
@@ -0,0 +1,65 @@
+#ifndef __GEOM_PATH_INTERSECTION_H
+#define __GEOM_PATH_INTERSECTION_H
+
+#include "path.h"
+
+#include "crossing.h"
+
+#include "sweep.h"
+
+namespace Geom {
+
+int winding(Path const &path, Point p);
+bool path_direction(Path const &p);
+
+inline bool contains(Path const & p, Point i, bool evenodd = true) {
+    return (evenodd ? winding(p, i) % 2 : winding(p, i)) != 0;
+}
+
+template<typename T>
+Crossings curve_sweep(Path const &a, Path const &b) {
+    T t;
+    Crossings ret;
+    std::vector<Rect> bounds_a = bounds(a), bounds_b = bounds(b);
+    std::vector<std::vector<unsigned> > ixs = sweep_bounds(bounds_a, bounds_b);
+    for(unsigned i = 0; i < a.size(); i++) {
+        for(std::vector<unsigned>::iterator jp = ixs[i].begin(); jp != ixs[i].end(); jp++) {
+            Crossings cc = t.crossings(a[i], b[*jp]);
+            offset_crossings(cc, i, *jp);
+            ret.insert(ret.end(), cc.begin(), cc.end());
+        }
+    }
+    return ret;
+}
+
+struct SimpleCrosser : public Crosser<Path> {
+    Crossings crossings(Curve const &a, Curve const &b);
+    Crossings crossings(Path const &a, Path const &b) { return curve_sweep<SimpleCrosser>(a, b); }
+    CrossingSet crossings(std::vector<Path> const &a, std::vector<Path> const &b) { return Crosser<Path>::crossings(a, b); }
+};
+
+struct MonoCrosser : public Crosser<Path> {
+    Crossings crossings(Path const &a, Path const &b) { return crossings(std::vector<Path>(1,a), std::vector<Path>(1,b))[0]; }
+    CrossingSet crossings(std::vector<Path> const &a, std::vector<Path> const &b);
+};
+
+typedef SimpleCrosser DefaultCrosser;
+
+std::vector<double> path_mono_splits(Path const &p);
+
+CrossingSet crossings_among(std::vector<Path> const & p);
+inline Crossings self_crossings(Path const & a) { return crossings_among(std::vector<Path>(1, a))[0]; }
+
+inline Crossings crossings(Path const & a, Path const & b) {
+    DefaultCrosser c = DefaultCrosser();
+    return c.crossings(a, b);
+}
+
+inline CrossingSet crossings(std::vector<Path> const & a, std::vector<Path> const & b) {
+    DefaultCrosser c = DefaultCrosser();
+    return c.crossings(a, b);
+}
+
+}
+
+#endif
diff --git a/src/2geom/path.cpp b/src/2geom/path.cpp
index 91868eb7e..f05d3b0cf 100644
--- a/src/2geom/path.cpp
+++ b/src/2geom/path.cpp
@@ -1,367 +1,258 @@
-/*

- * Path - Series of continuous curves

- *   

- * Copyright 2007  MenTaLguY <mental@rydia.net>

- *     

- * This library is free software; you can redistribute it and/or 

- * modify it either under the terms of the GNU Lesser General Public

- * License version 2.1 as published by the Free Software Foundation

- * (the "LGPL") or, at your option, under the terms of the Mozilla

- * Public License Version 1.1 (the "MPL"). If you do not alter this 

- * notice, a recipient may use your version of this file under either

- * the MPL or the LGPL.

- *            

- * You should have received a copy of the LGPL along with this library

- * in the file COPYING-LGPL-2.1; if not, write to the Free Software 

- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

- * You should have received a copy of the MPL along with this library 

- * in the file COPYING-MPL-1.1

- *                

- * The contents of this file are subject to the Mozilla Public License

- * Version 1.1 (the "License"); you may not use this file except in

- * compliance with the License. You may obtain a copy of the License at

- * http://www.mozilla.org/MPL/

- *                   

- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY

- * OF ANY KIND, either express or implied. See the LGPL or the MPL for

- * the specific language governing rights and limitations.

- */

-

-#include "path.h"

-

-namespace Geom {

-

-namespace {

-

-enum Cmp {

-  LESS_THAN=-1,

-  GREATER_THAN=1,

-  EQUAL_TO=0

-};

-

-template <typename T1, typename T2>

-inline Cmp cmp(T1 const &a, T2 const &b) {

-  if ( a < b ) {

-    return LESS_THAN;

-  } else if ( b < a ) {

-    return GREATER_THAN;

-  } else {

-    return EQUAL_TO;

-  }

-}

-

-}

-

-boost::optional<int> CurveHelpers::sbasis_winding(D2<SBasis> const &sb, Point p) {

-  Interval ix = bounds_fast(sb[X]);

-

-  if ( p[X] > ix.max() ) { /* ray does not intersect bbox */

-    return 0;

-  }

-

-  SBasis fy = sb[Y];

-  fy -= p[Y];

-

-  if (fy.empty()) { /* coincident horizontal segment */

-    return boost::optional<int>();

-  }

-

-  if ( p[X] < ix.min() ) { /* ray does not originate in bbox */

-    double y = p[Y];

-    /* winding determined by position of endpoints */

-    Cmp initial_to_ray = cmp(fy[0][0], y);

-    Cmp final_to_ray = cmp(fy[0][1], y);

-    switch (cmp(final_to_ray, initial_to_ray)) {

-    case GREATER_THAN:

-      /* exclude final endpoint */

-      return ( final_to_ray != EQUAL_TO );

-    case LESS_THAN:

-      /* exclude final endpoint */

-      return -( final_to_ray != EQUAL_TO );

-    default:

-      /* any intersections cancel out */

-      return 0;

-    }

-  } else { /* ray originates in bbox */

-    std::vector<double> ts = roots(fy);

-

-    static const unsigned MAX_DERIVATIVES=8;

-    boost::optional<SBasis> ds[MAX_DERIVATIVES];

-    ds[0] = derivative(fy);

-

-    /* winding determined by summing signs of derivatives at intersections */

-    int winding=0;

-    for ( std::vector<double>::iterator ti = ts.begin()

-        ; ti != ts.end()

-        ; ++ti )

-    { 

-      double t = *ti;

-      if ( sb[X](t) >= p[X] ) { /* root is ray intersection */

-        for ( boost::optional<SBasis> *di = ds

-            ; di != ( ds + MAX_DERIVATIVES )

-            ; ++di )

-        {

-          if (!*di) {

-            *di = derivative(**(di-1));

-          }

-          switch (cmp((**di)(t), 0)) {

-          case GREATER_THAN:

-            if ( t < 1 ) { /* exclude final endpoint */

-              winding += 1;

-            }

-            goto next_root;

-          case LESS_THAN:

-            if ( t < 1 ) { /* exclude final endpoint */

-              winding -= 1;

-            }

-            goto next_root;

-          default: (void)0;

-            /* give up */

-          };

-        }

-      } 

-next_root: (void)0;

-    }

-    

-    return winding;

-  }

-}

-

-Rect BezierHelpers::bounds(unsigned degree, Point const *points) {

-  Point min=points[0];

-  Point max=points[0];

-  for ( unsigned i = 1 ; i <= degree ; ++i ) {

-    for ( unsigned axis = 0 ; axis < 2 ; ++axis ) {

-      min[axis] = std::min(min[axis], points[i][axis]);

-      max[axis] = std::max(max[axis], points[i][axis]);

-    }

-  }

-  return Rect(min, max);

-}

-

-Point BezierHelpers::point_and_derivatives_at(Coord t,

-                                              unsigned degree,

-                                              Point const *points,

-                                              unsigned n_derivs,

-                                              Point *derivs)

-{

-  return Point(0,0); // TODO

-}

-

-Geom::Point

-BezierHelpers::subdivideArr(Coord t,              // Parameter value

-                            unsigned degree,      // Degree of bezier curve

-                            Geom::Point const *V, // Control pts

-                            Geom::Point *Left,    // RETURN left half ctl pts

-                            Geom::Point *Right)   // RETURN right half ctl pts

-{

-    Geom::Point Vtemp[degree+1][degree+1];

-

-    /* Copy control points	*/

-    std::copy(V, V+degree+1, Vtemp[0]);

-

-    /* Triangle computation	*/

-    for (unsigned i = 1; i <= degree; i++) {	

-        for (unsigned j = 0; j <= degree - i; j++) {

-            Vtemp[i][j] = lerp(t, Vtemp[i-1][j], Vtemp[i-1][j+1]);

-        }

-    }

-    

-    for (unsigned j = 0; j <= degree; j++)

-        Left[j]  = Vtemp[j][0];

-    for (unsigned j = 0; j <= degree; j++)

-        Right[j] = Vtemp[degree-j][j];

-

-    return (Vtemp[degree][0]);

-}

-

-void Path::swap(Path &other) {

-  std::swap(curves_, other.curves_);

-  std::swap(closed_, other.closed_);

-  std::swap(*final_, *other.final_);

-  curves_[curves_.size()-1] = final_;

-  other.curves_[other.curves_.size()-1] = other.final_;

-}

-

-Rect Path::bounds_fast() const {

-  Rect bounds=front().bounds_fast();

-  for ( const_iterator iter=++begin(); iter != end() ; ++iter ) {

-    bounds.unionWith(iter->bounds_fast());

-  }

-  return bounds;

-}

-

-Rect Path::bounds_exact() const {

-  Rect bounds=front().bounds_exact();

-  for ( const_iterator iter=++begin(); iter != end() ; ++iter ) {

-    bounds.unionWith(iter->bounds_exact());

-  }

-  return bounds;

-}

-

-int Path::winding(Point p) const {

-  int winding = 0;

-  boost::optional<Cmp> ignore = boost::optional<Cmp>();

-  for ( const_iterator iter = begin()

-      ; iter != end_closed()

-      ; ++iter )

-  {

-    boost::optional<int> w = iter->winding(p);

-    if (w) {

-      winding += *w;

-      ignore = boost::optional<Cmp>();

-    } else {

-      Point initial = iter->initialPoint();

-      Point final = iter->finalPoint();

-      switch (cmp(initial[X], final[X])) {

-      case GREATER_THAN:

-        if ( !ignore || *ignore != GREATER_THAN ) { /* ignore repeated */

-          winding += 1;

-          ignore = GREATER_THAN;

-        }

-        break;

-      case LESS_THAN:

-        if ( !ignore || *ignore != LESS_THAN ) { /* ignore repeated */

-          if ( p[X] < final[X] ) { /* ignore final point */

-            winding -= 1;

-            ignore = LESS_THAN;

-          }

-        }

-        break;

-      case EQUAL_TO:

-        /* always ignore null segments */

-        break;

-      }

-    }

-  }

-  return winding;

-}

-

-void Path::append(Curve const &curve) {

-  if ( curves_.front() != final_ && curve.initialPoint() != (*final_)[0] ) {

-    throw ContinuityError();

-  }

-  do_append(curve.duplicate());

-}

-

-void Path::append(D2<SBasis> const &curve) {

-  if ( curves_.front() != final_ ) {

-    for ( int i = 0 ; i < 2 ; ++i ) {

-      if ( curve[i][0][0] != (*final_)[0][i] ) {

-        throw ContinuityError();

-      }

-    }

-  }

-  do_append(new SBasisCurve(curve));

-}

-

-void Path::do_update(Sequence::iterator first_replaced,

-                     Sequence::iterator last_replaced,

-                     Sequence::iterator first,

-                    Sequence::iterator last)

-{

-  // note: modifies the contents of [first,last)

-

-  check_continuity(first_replaced, last_replaced, first, last);

-  delete_range(first_replaced, last_replaced);

-  if ( ( last - first ) == ( last_replaced - first_replaced ) ) {

-    std::copy(first, last, first_replaced);

-  } else {

-    // this approach depends on std::vector's behavior WRT iterator stability

-    curves_.erase(first_replaced, last_replaced);

-    curves_.insert(first_replaced, first, last);

-  }

-

-  if ( curves_.front() != final_ ) {

-    (*final_)[0] = back().finalPoint();

-    (*final_)[1] = front().initialPoint();

-  }

-}

-

-void Path::do_append(Curve *curve) {

-  if ( curves_.front() == final_ ) {

-    (*final_)[1] = curve->initialPoint();

-  }

-  curves_.insert(curves_.end()-1, curve);

-  (*final_)[0] = curve->finalPoint();

-}

-

-void Path::delete_range(Sequence::iterator first, Sequence::iterator last) {

-  for ( Sequence::iterator iter=first ; iter != last ; ++iter ) {

-    delete *iter;

-  }

-}

-

-void Path::check_continuity(Sequence::iterator first_replaced,

-                            Sequence::iterator last_replaced,

-                            Sequence::iterator first,

-                            Sequence::iterator last)

-{

-  if ( first != last ) {

-    if ( first_replaced != curves_.begin() ) {

-      if ( (*first_replaced)->initialPoint() != (*first)->initialPoint() ) {

-        throw ContinuityError();

-      }

-    }

-    if ( last_replaced != (curves_.end()-1) ) {

-      if ( (*(last_replaced-1))->finalPoint() != (*(last-1))->finalPoint() ) {

-        throw ContinuityError();

-      }

-    }

-  } else if ( first_replaced != last_replaced && first_replaced != curves_.begin() && last_replaced != curves_.end()-1) {

-    if ( (*first_replaced)->initialPoint() !=

-         (*(last_replaced-1))->finalPoint() )

-    {

-      throw ContinuityError();

-    }

-  }

-}

-

-Rect SBasisCurve::bounds_fast() const {

-  throw NotImplemented();

-  return Rect(Point(0,0), Point(0,0));

-}

-

-Rect SBasisCurve::bounds_exact() const {

-  throw NotImplemented();

-  return Rect(Point(0,0), Point(0,0));

-}

-

-Point SBasisCurve::pointAndDerivativesAt(Coord t, unsigned n_derivs, Point *derivs) const {

-  throw NotImplemented();

-  return Point(0,0);

-}

-

-Path const &SBasisCurve::subdivide(Coord t, Path &out) const {

-  throw NotImplemented();

-}

-

-Rect SVGEllipticalArc::bounds_fast() const {

-    throw NotImplemented();

-}

-Rect SVGEllipticalArc::bounds_exact() const {

-    throw NotImplemented();

-}

-

-Point SVGEllipticalArc::pointAndDerivativesAt(Coord t, unsigned n_derivs, Point *derivs) const {

-    throw NotImplemented();

-}

-

-D2<SBasis> SVGEllipticalArc::sbasis() const {

-    throw NotImplemented();

-}

-

-}

-

-/*

-  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=2:tabstop=8:softtabstop=2 :

-*/

-

+/*
+ * Path - Series of continuous curves
+ *   
+ * Copyright 2007  MenTaLguY <mental@rydia.net>
+ *     
+ * This library is free software; you can redistribute it and/or 
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this 
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *            
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software 
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library 
+ * in the file COPYING-MPL-1.1
+ *                
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *                   
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+#include "path.h"
+
+#include "ord.h"
+
+namespace Geom {
+
+int CurveHelpers::root_winding(Curve const &c, Point p) {
+    std::vector<double> ts = c.roots(p[Y], Y);
+
+    if(ts.empty()) return 0;
+
+    double const fudge = 0.01; //fudge factor used on first and last
+
+    std::sort(ts.begin(), ts.end());
+
+    // winding determined by crossings at roots
+    int wind=0;
+    // previous time
+    double pt = ts.front() - fudge;
+    for ( std::vector<double>::iterator ti = ts.begin()
+        ; ti != ts.end()
+        ; ++ti )
+    {
+        double t = *ti;
+        if ( t <= 0. || t >= 1. ) continue; //skip endpoint roots 
+        if ( c.valueAt(t, X) > p[X] ) { // root is ray intersection
+            // Get t of next:
+            std::vector<double>::iterator next = ti;
+            next++;
+            double nt;
+            if(next == ts.end()) nt = t + fudge; else nt = *next;
+            
+            // Check before in time and after in time for positions
+            // Currently we're using the average times between next and previous segs
+            Cmp after_to_ray =  cmp(c.valueAt((t + nt) / 2, Y), p[Y]);
+            Cmp before_to_ray = cmp(c.valueAt((t + pt) / 2, Y), p[Y]);
+            // if y is included, these will have opposite values, giving order.
+            Cmp dt = cmp(after_to_ray, before_to_ray);
+            if(dt != EQUAL_TO) //Should always be true, but yah never know..
+                wind += dt;
+            pt = t;
+        }
+    }
+    
+    return wind;
+}
+
+void Path::swap(Path &other) {
+  std::swap(curves_, other.curves_);
+  std::swap(closed_, other.closed_);
+  std::swap(*final_, *other.final_);
+  curves_[curves_.size()-1] = final_;
+  other.curves_[other.curves_.size()-1] = other.final_;
+}
+
+Rect Path::boundsFast() const {
+  Rect bounds=front().boundsFast();
+  for ( const_iterator iter=++begin(); iter != end() ; ++iter ) {
+    bounds.unionWith(iter->boundsFast());
+  }
+  return bounds;
+}
+
+Rect Path::boundsExact() const {
+  Rect bounds=front().boundsExact();
+  for ( const_iterator iter=++begin(); iter != end() ; ++iter ) {
+    bounds.unionWith(iter->boundsExact());
+  }
+  return bounds;
+}
+
+template<typename iter>
+iter inc(iter const &x, unsigned n) {
+  iter ret = x;
+  for(unsigned i = 0; i < n; i++)
+    ret++;
+  return ret;
+}
+
+//This assumes that you can't be perfect in your t-vals, and as such, tweaks the start
+void Path::appendPortionTo(Path &ret, double from, double to) const {
+  assert(from >= 0 && to >= 0);
+  if(to == 0) to = size()+0.999999;
+  if(from == to) { return; }
+  double fi, ti;
+  double ff = modf(from, &fi), tf = modf(to, &ti);
+  if(tf == 0) { ti--; tf = 1; }
+  const_iterator fromi = inc(begin(), (unsigned)fi);
+  if(fi == ti && from < to) {
+    Curve *v = fromi->portion(ff, tf);
+    ret.append(*v);
+    delete v;
+    return;
+  }
+  const_iterator toi   = inc(begin(), (unsigned)ti);
+  if(ff != 1.) {
+    Curve *fromv = fromi->portion(ff, 1.);
+    //fromv->setInitial(ret.finalPoint());
+    ret.append(*fromv);
+    delete fromv;
+  }
+  if(from > to) {
+    const_iterator ender = end();
+    if(ender->initialPoint() == ender->finalPoint()) ender++;
+    ret.insert(ret.end(), ++fromi, ender);
+    ret.insert(ret.end(), begin(), toi);
+  } else {
+    ret.insert(ret.end(), ++fromi, toi);
+  }
+  Curve *tov = toi->portion(0., tf);
+  ret.append(*tov);
+  delete tov;
+}
+
+const double eps = .1;
+
+void Path::append(Curve const &curve) {
+  if ( curves_.front() != final_ && !near(curve.initialPoint(), (*final_)[0], eps) ) {
+    throw ContinuityError();
+  }
+  do_append(curve.duplicate());
+}
+
+void Path::append(D2<SBasis> const &curve) {
+  if ( curves_.front() != final_ ) {
+    for ( int i = 0 ; i < 2 ; ++i ) {
+      if ( !near(curve[i][0][0], (*final_)[0][i], eps) ) {
+        throw ContinuityError();
+      }
+    }
+  }
+  do_append(new SBasisCurve(curve));
+}
+
+void Path::do_update(Sequence::iterator first_replaced,
+                     Sequence::iterator last_replaced,
+                     Sequence::iterator first,
+                    Sequence::iterator last)
+{
+  // note: modifies the contents of [first,last)
+
+  check_continuity(first_replaced, last_replaced, first, last);
+  delete_range(first_replaced, last_replaced);
+  if ( ( last - first ) == ( last_replaced - first_replaced ) ) {
+    std::copy(first, last, first_replaced);
+  } else {
+    // this approach depends on std::vector's behavior WRT iterator stability
+    curves_.erase(first_replaced, last_replaced);
+    curves_.insert(first_replaced, first, last);
+  }
+
+  if ( curves_.front() != final_ ) {
+    final_->setPoint(0, back().finalPoint());
+    final_->setPoint(1, front().initialPoint());
+  }
+}
+
+void Path::do_append(Curve *curve) {
+  if ( curves_.front() == final_ ) {
+    final_->setPoint(1, curve->initialPoint());
+  }
+  curves_.insert(curves_.end()-1, curve);
+  final_->setPoint(0, curve->finalPoint());
+}
+
+void Path::delete_range(Sequence::iterator first, Sequence::iterator last) {
+  for ( Sequence::iterator iter=first ; iter != last ; ++iter ) {
+    delete *iter;
+  }
+}
+
+void Path::check_continuity(Sequence::iterator first_replaced,
+                            Sequence::iterator last_replaced,
+                            Sequence::iterator first,
+                            Sequence::iterator last)
+{
+  if ( first != last ) {
+    if ( first_replaced != curves_.begin() ) {
+      if ( !near( (*first_replaced)->initialPoint(), (*first)->initialPoint(), eps ) ) {
+        throw ContinuityError();
+      }
+    }
+    if ( last_replaced != (curves_.end()-1) ) {
+      if ( !near( (*(last_replaced-1))->finalPoint(), (*(last-1))->finalPoint(), eps ) ) {
+        throw ContinuityError();
+      }
+    }
+  } else if ( first_replaced != last_replaced && first_replaced != curves_.begin() && last_replaced != curves_.end()-1) {
+    if ( !near((*first_replaced)->initialPoint(), (*(last_replaced-1))->finalPoint(), eps ) ) {
+      throw ContinuityError();
+    }
+  }
+}
+
+Rect SVGEllipticalArc::boundsFast() const {
+    throw NotImplemented();
+}
+Rect SVGEllipticalArc::boundsExact() const {
+    throw NotImplemented();
+}
+Rect SVGEllipticalArc::boundsLocal(Interval i, unsigned deg) const {
+    //throw NotImplemented();
+}
+
+std::vector<Point> SVGEllipticalArc::pointAndDerivatives(Coord t, unsigned n) const {
+    throw NotImplemented();
+}
+
+std::vector<double> SVGEllipticalArc::roots(double v, Dim2 d) const {
+    //throw NotImplemented();
+}
+
+D2<SBasis> SVGEllipticalArc::toSBasis() const {
+    return D2<SBasis>(Linear(initial_[X], final_[X]), Linear(initial_[Y], final_[Y]));
+}
+
+}
+
+/*
+  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=2:tabstop=8:softtabstop=2 :
+*/
diff --git a/src/2geom/path.h b/src/2geom/path.h
index 31a7173b7..6f39eb7ef 100644
--- a/src/2geom/path.h
+++ b/src/2geom/path.h
@@ -35,15 +35,21 @@
 #include <algorithm>
 #include <exception>
 #include <stdexcept>
-#include <boost/optional/optional.hpp>
 #include "d2.h"
-#include "bezier-to-sbasis.h"
+#include "matrix.h"
+#include "bezier.h"
+#include "crossing.h"
 
 namespace Geom {
 
-class Path;
+class Curve;
 
-class Curve {
+struct CurveHelpers {
+protected:
+  static int root_winding(Curve const &c, Point p);
+};
+
+class Curve : private CurveHelpers {
 public:
   virtual ~Curve() {}
 
@@ -52,112 +58,211 @@ public:
 
   virtual Curve *duplicate() const = 0;
 
-  virtual Rect bounds_fast() const = 0;
-  virtual Rect bounds_exact() const = 0;
-
-  virtual boost::optional<int> winding(Point p) const = 0;
+  virtual Rect boundsFast() const = 0;
+  virtual Rect boundsExact() const = 0;
+  virtual Rect boundsLocal(Interval i, unsigned deg) const = 0;
+  Rect boundsLocal(Interval i) const { return boundsLocal(i, 0); }
+  
+  virtual std::vector<double> roots(double v, Dim2 d) const = 0;
 
-  virtual Path const &subdivide(Coord t, Path &out) const = 0;
+  virtual int winding(Point p) const { return root_winding(*this, p); }
 
-  Point pointAt(Coord t) const { return pointAndDerivativesAt(t, 0, NULL); }
-  virtual Point pointAndDerivativesAt(Coord t, unsigned n, Point *ds) const = 0;
-  virtual D2<SBasis> sbasis() const = 0;
+  //mental: review these
+  virtual Curve *portion(double f, double t) const = 0;
+  virtual Curve *reverse() const { return portion(1, 0); }
+  virtual Curve *derivative() const = 0;
+  
+  virtual void setInitial(Point v) = 0;
+  virtual void setFinal(Point v) = 0;
+  
+  virtual Curve *transformed(Matrix const &m) const = 0;
+  
+  virtual Point pointAt(Coord t) const { return pointAndDerivatives(t, 1).front(); }
+  virtual Coord valueAt(Coord t, Dim2 d) const { return pointAt(t)[d]; }
+  virtual std::vector<Point> pointAndDerivatives(Coord t, unsigned n) const = 0;
+  virtual D2<SBasis> toSBasis() const = 0;
 };
 
-struct CurveHelpers {
-protected:
-  static boost::optional<int> sbasis_winding(D2<SBasis> const &sbasis, Point p);
-};
+class SBasisCurve : public Curve {
+private:
+  SBasisCurve();
+  D2<SBasis> inner;
+public:
+  explicit SBasisCurve(D2<SBasis> const &sb) : inner(sb) {}
+  explicit SBasisCurve(Curve const &other) : inner(other.toSBasis()) {}
+  Curve *duplicate() const { return new SBasisCurve(*this); }
 
-struct BezierHelpers {
-protected:
-  static Rect bounds(unsigned degree, Point const *points);
-  static Point point_and_derivatives_at(Coord t,
-                                        unsigned degree, Point const *points,
-                                        unsigned n_derivs, Point *derivs);
-  static Point subdivideArr(Coord t, unsigned degree, Point const *V,
-                         Point *Left, Point *Right);
+  Point initialPoint() const    { return inner.at0(); }
+  Point finalPoint() const      { return inner.at1(); }
+  Point pointAt(Coord t) const  { return inner.valueAt(t); }
+  std::vector<Point> pointAndDerivatives(Coord t, unsigned n) const {
+      return inner.valueAndDerivatives(t, n);
+  }
+  double valueAt(Coord t, Dim2 d) const { return inner[d].valueAt(t); }
+  
+  void setInitial(Point v) { for(unsigned d = 0; d < 2; d++) { inner[d][0][0] = v[d]; } }
+  void setFinal(Point v)   { for(unsigned d = 0; d < 2; d++) { inner[d][0][1] = v[d]; } }
+
+  Rect boundsFast() const  { return bounds_fast(inner); }
+  Rect boundsExact() const { return bounds_exact(inner); }
+  Rect boundsLocal(Interval i, unsigned deg) const { return bounds_local(inner, i, deg); }
+
+  std::vector<double> roots(double v, Dim2 d) const { return Geom::roots(inner[d] - v); }
+
+  Curve *portion(double f, double t) const {
+    return new SBasisCurve(Geom::portion(inner, f, t));
+  }
+  
+  Curve *transformed(Matrix const &m) const {
+    return new SBasisCurve(inner * m);
+  }
+  
+  Curve *derivative() const {
+    return new SBasisCurve(Geom::derivative(inner));
+  }
+  
+  D2<SBasis> toSBasis() const { return inner; }
 
 };
 
-template <unsigned bezier_degree>
-class Bezier : public Curve, private CurveHelpers, private BezierHelpers {
+template <unsigned order>
+class BezierCurve : public Curve {
+private:
+  D2<Bezier<order> > inner;
 public:
   template <unsigned required_degree>
-  static void assert_degree(Bezier<required_degree> const *) {}
+  static void assert_degree(BezierCurve<required_degree> const *) {}
+
+  BezierCurve() {}
 
-  Bezier() {}
+  explicit BezierCurve(D2<Bezier<order> > const &x) : inner(x) {}
+  
+  BezierCurve(Bezier<order> x, Bezier<order> y) : inner(x, y) {}
 
   // default copy
   // default assign
 
-  Bezier(Point c0, Point c1) {
+  BezierCurve(Point c0, Point c1) {
     assert_degree<1>(this);
-    c_[0] = c0;
-    c_[1] = c1;
+    for(unsigned d = 0; d < 2; d++)
+        inner[d] = Bezier<order>(c0[d], c1[d]);
   }
 
-  Bezier(Point c0, Point c1, Point c2) {
+  BezierCurve(Point c0, Point c1, Point c2) {
     assert_degree<2>(this);
-    c_[0] = c0;
-    c_[1] = c1;
-    c_[2] = c2;
+    for(unsigned d = 0; d < 2; d++)
+        inner[d] = Bezier<order>(c0[d], c1[d], c2[d]);
   }
 
-  Bezier(Point c0, Point c1, Point c2, Point c3) {
+  BezierCurve(Point c0, Point c1, Point c2, Point c3) {
     assert_degree<3>(this);
-    c_[0] = c0;
-    c_[1] = c1;
-    c_[2] = c2;
-    c_[3] = c3;
+    for(unsigned d = 0; d < 2; d++)
+        inner[d] = Bezier<order>(c0[d], c1[d], c2[d], c3[d]);
   }
 
-  unsigned degree() const { return bezier_degree; }
+  unsigned degree() const { return order; }
 
-  Curve *duplicate() const { return new Bezier(*this); }
+  Curve *duplicate() const { return new BezierCurve(*this); }
 
-  Point initialPoint() const { return c_[0]; }
-  Point finalPoint() const { return c_[bezier_degree]; }
+  Point initialPoint() const { return inner.at0(); }
+  Point finalPoint() const { return inner.at1(); }
 
-  Point &operator[](int index) { return c_[index]; }
-  Point const &operator[](int index) const { return c_[index]; }
+  void setInitial(Point v) { setPoint(0, v); }
+  void setFinal(Point v)   { setPoint(1, v); }
 
-  Rect bounds_fast() const { return bounds(bezier_degree, c_); }
-  Rect bounds_exact() const { return bounds(bezier_degree, c_); }
+  void setPoint(unsigned ix, Point v) { inner[X].setPoint(ix, v[X]); inner[Y].setPoint(ix, v[Y]); }
+  Point const operator[](unsigned ix) const { return Point(inner[X][ix], inner[Y][ix]); }
+  
+  Rect boundsFast() const { return bounds_fast(inner); }
+  Rect boundsExact() const { return bounds_exact(inner); }
+  Rect boundsLocal(Interval i, unsigned deg) const {
+      if(i.min() == 0 && i.max() == 1) return boundsFast();
+      if(deg == 0) return bounds_local(inner, i);
+      // TODO: UUUUUUGGGLLY
+      if(deg == 1 && order > 1) return Rect(bounds_local(Geom::derivative(inner[X]), i),
+                                            bounds_local(Geom::derivative(inner[Y]), i));
+      return Rect(Interval(0,0), Interval(0,0));
+  }
+//TODO: local
 
-  boost::optional<int> winding(Point p) const {
-    return sbasis_winding(sbasis(), p);
+//TODO: implement next 3 natively
+  int winding(Point p) const {
+    return SBasisCurve(toSBasis()).winding(p);
   }
 
-  Path const &subdivide(Coord t, Path &out) const;
+  std::vector<double>
+  roots(double v, Dim2 d) const {
+      return (inner[d] - v).roots();
+  }
   
-  Point pointAndDerivativesAt(Coord t, unsigned n_derivs, Point *derivs)
-  const
-  {
-    return point_and_derivatives_at(t, bezier_degree, c_, n_derivs, derivs);
+  void setPoints(std::vector<Point> ps) {
+    for(unsigned i = 0; i <= order; i++) {
+      setPoint(i, ps[i]);
+    }
+  }
+  std::vector<Point> points() const { return bezier_points(inner); }
+  
+  std::pair<BezierCurve<order>, BezierCurve<order> > subdivide(Coord t) const {
+    std::pair<Bezier<order>, Bezier<order> > sx = inner[X].subdivide(t), sy = inner[Y].subdivide(t);
+    return std::pair<BezierCurve<order>, BezierCurve<order> >(
+               BezierCurve<order>(sx.first, sy.first),
+               BezierCurve<order>(sx.second, sy.second));
+  }
+  
+  Curve *portion(double f, double t) const {
+    return new BezierCurve(Geom::portion(inner, f, t));
   }
 
-  D2<SBasis> sbasis() const {
-    return handles_to_sbasis<bezier_degree>(c_);
+  Curve *reverse() const {
+    return new BezierCurve(Geom::reverse(inner));
   }
 
-protected:
-  Bezier(Point c[]) {
-    std::copy(c, c+bezier_degree+1, c_);
+  Curve *transformed(Matrix const &m) const {
+    BezierCurve *ret = new BezierCurve();
+    std::vector<Point> ps = points();
+    for(unsigned i = 0;  i <= order; i++) ps[i] = ps[i] * m;
+    ret->setPoints(ps);
+    return ret;
   }
 
-private:
-  Point c_[bezier_degree+1];
+  Curve *derivative() const {
+     if(order > 1)
+        return new BezierCurve<order-1>(Geom::derivative(inner[X]), Geom::derivative(inner[Y]));
+     else if (order == 1) {
+        double dx = inner[X][1] - inner[X][0], dy = inner[Y][1] - inner[Y][0];
+        if(dx == 0) return new BezierCurve<1>(Point(0,0), Point(0,0));
+        double slope = dy / dx;
+        Geom::Point pnt;
+        if(slope == 0) pnt = Geom::Point(0, 0); else pnt = Geom::Point(slope, 1./slope);
+        return new BezierCurve<1>(pnt, pnt);
+     }
+  }
+
+  Point pointAt(double t) const { return inner.valueAt(t); } 
+  std::vector<Point> pointAndDerivatives(Coord t, unsigned n) const { return inner.valueAndDerivatives(t, n); }
+
+  double valueAt(double t, Dim2 d) const { return inner[d].valueAt(t); }
+  
+  D2<SBasis> toSBasis() const {return inner.toSBasis(); }
+
+protected:
+  BezierCurve(Point c[]) {
+    Coord x[order+1], y[order+1];
+    for(unsigned i = 0; i <= order; i++) {
+        x[i] = c[i][X]; y[i] = c[i][Y];
+    }
+    inner = Bezier<order>(x, y);
+  }
 };
 
-// Bezier<0> is meaningless; specialize it out
-template<> class Bezier<0> { Bezier(); };
+// BezierCurve<0> is meaningless; specialize it out
+template<> class BezierCurve<0> : public BezierCurve<1> { public: BezierCurve(); BezierCurve(Bezier<0> x, Bezier<0> y); };
 
-typedef Bezier<1> LineSegment;
-typedef Bezier<2> QuadraticBezier;
-typedef Bezier<3> CubicBezier;
+typedef BezierCurve<1> LineSegment;
+typedef BezierCurve<2> QuadraticBezier;
+typedef BezierCurve<3> CubicBezier;
 
-class SVGEllipticalArc : public Curve, private CurveHelpers {
+class SVGEllipticalArc : public Curve {
 public:
   SVGEllipticalArc() {}
 
@@ -173,18 +278,55 @@ public:
   Point initialPoint() const { return initial_; }
   Point finalPoint() const { return final_; }
 
-  Rect bounds_fast() const;
-  Rect bounds_exact() const;
+  void setInitial(Point v) { initial_ = v; }
+  void setFinal(Point v) { final_ = v; }
+
+  //TODO: implement funcs
 
-  boost::optional<int> winding(Point p) const {
-    return sbasis_winding(sbasis(), p);
+  Rect boundsFast() const;
+  Rect boundsExact() const;
+  Rect boundsLocal(Interval i, unsigned deg) const;
+
+  int winding(Point p) const {
+    return SBasisCurve(toSBasis()).winding(p);
   }
+  
+  std::vector<double> roots(double v, Dim2 d) const;
 
-  Path const &subdivide(Coord t, Path &out) const;
+  inline std::pair<SVGEllipticalArc, SVGEllipticalArc>
+  subdivide(Coord t) {
+    SVGEllipticalArc a(*this), b(*this);
+    a.final_ = b.initial_ = pointAt(t);
+    return std::pair<SVGEllipticalArc, SVGEllipticalArc>(a, b);
+  }
+  
+  Curve *portion(double f, double t) const {
+    SVGEllipticalArc *ret = new SVGEllipticalArc (*this);
+    ret->initial_ = pointAt(f);
+    ret->final_ = pointAt(t);
+    return ret;
+  }
+  
+  Curve *reverse(double f, double t) const {
+    SVGEllipticalArc *ret = new SVGEllipticalArc (*this);
+    ret->initial_ = final_;
+    ret->final_ = initial_;
+    return ret;
+  }
   
-  Point pointAndDerivativesAt(Coord t, unsigned n_derivs, Point *derivs) const;
+  //TODO: this next def isn't right
+  Curve *transformed(Matrix const & m) const {
+    SVGEllipticalArc *ret = new SVGEllipticalArc (*this);
+    ret->initial_ = initial_ * m;
+    ret->final_ = final_ * m;
+    return ret;
+  }
+  
+  Curve *derivative() const { throw NotImplemented(); }
+  
+  std::vector<Point> pointAndDerivatives(Coord t, unsigned n) const;
 
-  D2<SBasis> sbasis() const;
+  D2<SBasis> toSBasis() const;
 
 private:
   Point initial_;
@@ -196,39 +338,6 @@ private:
   Point final_;
 };
 
-class SBasisCurve : public Curve, private CurveHelpers {
-private:
-  SBasisCurve();
-public:
-  explicit SBasisCurve(D2<SBasis> const &coeffs)
-  : coeffs_(coeffs) {}
-
-  Point initialPoint() const {
-    return Point(coeffs_[X][0][0], coeffs_[Y][0][0]);
-  }
-  Point finalPoint() const {
-    return Point(coeffs_[X][0][1], coeffs_[Y][0][1]);
-  }
-
-  Curve *duplicate() const { return new SBasisCurve(*this); }
-
-  Rect bounds_fast() const;
-  Rect bounds_exact() const;
-
-  boost::optional<int> winding(Point p) const {
-    return sbasis_winding(coeffs_, p);
-  }
-
-  Path const &subdivide(Coord t, Path &out) const;
-
-  Point pointAndDerivativesAt(Coord t, unsigned n_derivs, Point *derivs) const;
-
-  D2<SBasis> sbasis() const { return coeffs_; }
-  
-private:
-  D2<SBasis> coeffs_;
-};
-
 template <typename IteratorImpl>
 class BaseIterator
 : public std::iterator<std::forward_iterator_tag, Curve const>
@@ -253,6 +362,7 @@ public:
     ++impl_;
     return *this;
   }
+  
   BaseIterator operator++(int) {
     BaseIterator old=*this;
     ++(*this);
@@ -282,7 +392,7 @@ public:
   }
 
   Curve *operator*() const { return (*impl_)->duplicate(); }
-
+  
   DuplicatingIterator &operator++() {
     ++impl_;
     return *this;
@@ -378,20 +488,77 @@ public:
   bool closed() const { return closed_; }
   void close(bool closed=true) { closed_ = closed; }
 
-  int winding(Point p) const;
-
-  Rect bounds_fast() const;
-  Rect bounds_exact() const;
+  Rect boundsFast() const;
+  Rect boundsExact() const;
 
   Piecewise<D2<SBasis> > toPwSb() const {
     Piecewise<D2<SBasis> > ret;
     ret.push_cut(0);
-    for(unsigned i = 0; i < size() + (closed_ ? 1 : 0); i++) {
-      ret.push(curves_[i]->sbasis(), i+1);
+    unsigned i = 1;
+    for(const_iterator it = begin(); it != end_default(); ++it, i++) { 
+      ret.push(it->toSBasis(), i);
+    }
+    return ret;
+  }
+  
+  Path operator*(Matrix const &m) const {
+    Path ret;
+    for(const_iterator it = begin(); it != end(); ++it) {
+      Curve *temp = it->transformed(m);
+      //Possible point of discontinuity?
+      ret.append(*temp);
+      delete temp;
     }
     return ret;
   }
 
+  Point pointAt(double t) const {
+    if(empty()) return Point(0,0);
+    double i, f = modf(t, &i);
+    if(i == size() && f == 0) { i--; }
+    assert(i >= 0 && i <= size()); 
+    return (*this)[unsigned(i)].pointAt(f);
+  }
+
+  double valueAt(double t, Dim2 d) const {
+    if(empty()) return 0;
+    double i, f = modf(t, &i);
+    if(i == size() && f == 0) { i--; }
+    assert(i >= 0 && i <= size());
+    return (*this)[unsigned(i)].valueAt(f, d);
+  }
+
+  std::vector<double> roots(double v, Dim2 d) const {
+    std::vector<double> res;
+    for(const_iterator it = begin(); it != end_closed(); ++it) {
+      std::vector<double> temp = it->roots(v, d);
+      res.insert(res.end(), temp.begin(), temp.end());
+    }
+    return res;
+  }
+
+  void appendPortionTo(Path &p, double f, double t) const;
+  
+  Path portion(double f, double t) const {
+    Path ret;
+    ret.close(false);
+    appendPortionTo(ret, f, t);
+    return ret;
+  }
+  Path portion(Interval i) const { return portion(i.min(), i.max()); }
+  
+  Path reverse() const {
+    Path ret;
+    ret.close(closed_);
+    for(int i = size() - (closed_ ? 0 : 1); i >= 0; i--) {
+      //TODO: do we really delete?
+      Curve *temp = (*this)[i].reverse();
+      ret.append(*temp);
+      delete temp;
+    }
+    return ret;
+  } 
+  
   void insert(iterator pos, Curve const &curve) {
     Sequence source(1, curve.duplicate());
     try {
@@ -481,14 +648,14 @@ public:
 
   void start(Point p) {
     clear();
-    (*final_)[0] = (*final_)[1] = p;
+    final_->setPoint(0, p);
+    final_->setPoint(1, p);
   }
 
   Point initialPoint() const { return (*final_)[1]; }
   Point finalPoint() const { return (*final_)[0]; }
 
   void append(Curve const &curve);
-
   void append(D2<SBasis> const &curve);
 
   template <typename CurveType, typename A>
@@ -573,37 +740,52 @@ inline static Piecewise<D2<SBasis> > paths_to_pw(std::vector<Path> paths) {
     return ret;
 }
 
-template <unsigned bezier_degree>
-inline Path const &Bezier<bezier_degree>::subdivide(Coord t, Path &out) const {
-  Bezier a, b;
-  subdivideArr(t, bezier_degree, c_, a.c_, b.c_);
-  out.clear();
-  out.close(false);
-  out.append(a);
-  out.append(b);
-  return out;
-}
+/*
+class PathPortion : public Curve {
+  Path *source;
+  double f, t;
+  boost::optional<Path> result;
+  
+  public:
+  double from() const { return f; }
+  double to() const { return t; }
+    
+  explicit PathPortion(Path *s, double fp, double tp) : source(s), f(fp), t(tp) {}
+  Curve *duplicate() const { return new PathPortion(*this); }
+  
+  Point initialPoint() const { return source->pointAt(f); }
+  Point finalPoint() const { return source->pointAt(t); }
+  
+  Path actualPath() {
+    if(!result) *result = source->portion(f, t);
+    return *result;
+  }
+  
+  Rect boundsFast() const { return actualPath().boundsFast; }
+  Rect boundsExact() const { return actualPath().boundsFast; }
+  Rect boundsLocal(Interval i) const { throw NotImplemented(); }
 
-inline Path const &SVGEllipticalArc::subdivide(Coord t, Path &out) const {
-  SVGEllipticalArc a;
-  SVGEllipticalArc b;
-  a.rx_ = b.rx_ = rx_;
-  a.ry_ = b.ry_ = ry_;
-  a.x_axis_rotation_ = b.x_axis_rotation_ = x_axis_rotation_;
-  a.initial_ = initial_;
-  a.final_ = b.initial_ = pointAt(t);
-  b.final_ = final_;
-  out.clear();
-  out.close(false);
-  out.append(a);
-  out.append(b);
-  return out;
-}
+  std::vector<double> roots(double v, Dim2 d) const = 0;
 
-class Set {
-public:
-private:
+  virtual int winding(Point p) const { return root_winding(*this, p); }
+
+  virtual Curve *portion(double f, double t) const = 0;
+  virtual Curve *reverse() const { return portion(1, 0); }
+
+  virtual Crossings crossingsWith(Curve const & other) const;
+  
+  virtual void setInitial(Point v) = 0;
+  virtual void setFinal(Point v) = 0;
+  
+  virtual Curve *transformed(Matrix const &m) const = 0;
+  
+  virtual Point pointAt(Coord t) const { return pointAndDerivatives(t, 1).front(); }
+  virtual Coord valueAt(Coord t, Dim2 d) const { return pointAt(t)[d]; }
+  virtual std::vector<Point> pointAndDerivatives(Coord t, unsigned n) const = 0;
+  virtual D2<SBasis> toSBasis() const = 0;
+    
 };
+*/
 
 }
 
diff --git a/src/2geom/poly-dk-solve.cpp b/src/2geom/poly-dk-solve.cpp
index dde230dc4..136687ae8 100644
--- a/src/2geom/poly-dk-solve.cpp
+++ b/src/2geom/poly-dk-solve.cpp
@@ -1,64 +1,64 @@
-#include "poly-dk-solve.h"
-#include <iterator>
-
-/*** implementation of the Durand-Kerner method.  seems buggy*/
-
-std::complex<double> evalu(Poly const & p, std::complex<double> x) {
-    std::complex<double> result = 0;
-    std::complex<double> xx = 1;
-    
-    for(unsigned i = 0; i < p.size(); i++) {
-        result += p[i]*xx;
-        xx *= x;
-    }
-    return result;
-}
-
-std::vector<std::complex<double> > DK(Poly const & ply, const double tol) {
-    std::vector<std::complex<double> > roots;
-    const int N = ply.degree();
-    
-    std::complex<double> b(0.4, 0.9);
-    std::complex<double> p = 1;
-    for(int i = 0; i < N; i++) {
-        roots.push_back(p);
-        p *= b;
-    }
-    assert(roots.size() == ply.degree());
-    
-    double error = 0;
-    int i;
-    for( i = 0; i < 30; i++) {
-        error = 0;
-        for(int r_i = 0; r_i < N; r_i++) {
-            std::complex<double> denom = 1;
-            std::complex<double> R = roots[r_i];
-            for(int d_i = 0; d_i < N; d_i++) {
-                if(r_i != d_i)
-                    denom *= R-roots[d_i];
-            }
-            assert(norm(denom) != 0);
-	    std::complex<double> dr = evalu(ply, R)/denom;
-            error += norm(dr);
-            roots[r_i] = R - dr;
-        }
-        /*std::copy(roots.begin(), roots.end(), std::ostream_iterator<std::complex<double> >(std::cout, ",\t"));
-          std::cout << std::endl;*/
-        if(error < tol)
-            break;
-    }
-    //std::cout << error << ", " << i<< std::endl;
-    return roots;
-}
-
-
-/*
-  Local Variables:
-  mode:c++
-  c-file-style:"stroustrup"
-  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
-  indent-tabs-mode:nil
-  fill-column:99
-  End:
-*/
-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
+#include "poly-dk-solve.h"

+#include <iterator>

+

+/*** implementation of the Durand-Kerner method.  seems buggy*/

+

+std::complex<double> evalu(Poly const & p, std::complex<double> x) {

+    std::complex<double> result = 0;

+    std::complex<double> xx = 1;

+    

+    for(unsigned i = 0; i < p.size(); i++) {

+        result += p[i]*xx;

+        xx *= x;

+    }

+    return result;

+}

+

+std::vector<std::complex<double> > DK(Poly const & ply, const double tol) {

+    std::vector<std::complex<double> > roots;

+    const int N = ply.degree();

+    

+    std::complex<double> b(0.4, 0.9);

+    std::complex<double> p = 1;

+    for(int i = 0; i < N; i++) {

+        roots.push_back(p);

+        p *= b;

+    }

+    assert(roots.size() == ply.degree());

+    

+    double error = 0;

+    int i;

+    for( i = 0; i < 30; i++) {

+        error = 0;

+        for(int r_i = 0; r_i < N; r_i++) {

+            std::complex<double> denom = 1;

+            std::complex<double> R = roots[r_i];

+            for(int d_i = 0; d_i < N; d_i++) {

+                if(r_i != d_i)

+                    denom *= R-roots[d_i];

+            }

+            assert(norm(denom) != 0);

+	    std::complex<double> dr = evalu(ply, R)/denom;

+            error += norm(dr);

+            roots[r_i] = R - dr;

+        }

+        /*std::copy(roots.begin(), roots.end(), std::ostream_iterator<std::complex<double> >(std::cout, ",\t"));

+          std::cout << std::endl;*/

+        if(error < tol)

+            break;

+    }

+    //std::cout << error << ", " << i<< std::endl;

+    return roots;

+}

+

+

+/*

+  Local Variables:

+  mode:c++

+  c-file-style:"stroustrup"

+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))

+  indent-tabs-mode:nil

+  fill-column:99

+  End:

+*/

+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :

diff --git a/src/2geom/poly-laguerre-solve.cpp b/src/2geom/poly-laguerre-solve.cpp
index 83f049286..fb3c2075b 100644
--- a/src/2geom/poly-laguerre-solve.cpp
+++ b/src/2geom/poly-laguerre-solve.cpp
@@ -1,147 +1,147 @@
-#include "poly-laguerre-solve.h"
-#include <iterator>
-
-typedef std::complex<double> cdouble;
-
-cdouble laguerre_internal_complex(Poly const & p,
-                                  double x0,
-                                  double tol,
-                                  bool & quad_root) {
-    cdouble a = 2*tol;
-    cdouble xk = x0;
-    double n = p.degree();
-    quad_root = false;
-    const unsigned shuffle_rate = 10;
-    static double shuffle[] = {0, 0.5, 0.25, 0.75, 0.125, 0.375, 0.625, 0.875, 1.0};
-    unsigned shuffle_counter = 0;
-    while(std::norm(a) > (tol*tol)) {
-        //std::cout << "xk = " << xk << std::endl;
-        cdouble b = p.back();
-        cdouble d = 0, f = 0;
-        double err = abs(b);
-        double abx = abs(xk);
-        for(int j = p.size()-2; j >= 0; j--) {
-            f = xk*f + d;
-            d = xk*d + b;
-            b = xk*b + p[j];
-            err = abs(b) + abx*err;
-        }
-        
-        err *= 1e-7; // magic epsilon for convergence, should be computed from tol
-        
-        cdouble px = b;
-        if(abs(b) < err)
-            return xk;
-        //if(std::norm(px) < tol*tol)
-        //    return xk;
-        cdouble G = d / px;
-        cdouble H = G*G - f / px;
-        
-        //std::cout << "G = " << G << "H = " << H;
-        cdouble radicand = (n - 1)*(n*H-G*G);
-        //assert(radicand.real() > 0);
-        if(radicand.real() < 0)
-            quad_root = true;
-        //std::cout << "radicand = " << radicand << std::endl;
-        if(G.real() < 0) // here we try to maximise the denominator avoiding cancellation
-            a = - sqrt(radicand);
-        else
-            a = sqrt(radicand);
-        //std::cout << "a = " << a << std::endl;
-        a = n / (a + G);
-        //std::cout << "a = " << a << std::endl;
-        if(shuffle_counter % shuffle_rate == 0)
-            ;//a *= shuffle[shuffle_counter / shuffle_rate];
-        xk -= a;
-        shuffle_counter++;
-        if(shuffle_counter >= 90)
-            break;
-    }
-    //std::cout << "xk = " << xk << std::endl;
-    return xk;
-}
-
-double laguerre_internal(Poly const & p,
-                         Poly const & pp,
-                         Poly const & ppp,
-                         double x0,
-                         double tol,
-                         bool & quad_root) {
-    double a = 2*tol;
-    double xk = x0;
-    double n = p.degree();
-    quad_root = false;
-    while(a*a > (tol*tol)) {
-        //std::cout << "xk = " << xk << std::endl;
-        double px = p(xk);
-        if(px*px < tol*tol)
-            return xk;
-        double G = pp(xk) / px;
-        double H = G*G - ppp(xk) / px;
-        
-        //std::cout << "G = " << G << "H = " << H;
-        double radicand = (n - 1)*(n*H-G*G);
-        assert(radicand > 0);
-        //std::cout << "radicand = " << radicand << std::endl;
-        if(G < 0) // here we try to maximise the denominator avoiding cancellation
-            a = - sqrt(radicand);
-        else
-            a = sqrt(radicand);
-        //std::cout << "a = " << a << std::endl;
-        a = n / (a + G);
-        //std::cout << "a = " << a << std::endl;
-        xk -= a;
-    }
-    //std::cout << "xk = " << xk << std::endl;
-    return xk;
-}
-
-
-std::vector<cdouble > 
-laguerre(Poly p, const double tol) {
-    std::vector<cdouble > solutions;
-    //std::cout << "p = " << p << " = ";
-    while(p.size() > 1)
-    {
-        double x0 = 0;
-        bool quad_root = false;
-        cdouble sol = laguerre_internal_complex(p, x0, tol, quad_root);
-        //if(abs(sol) > 1) break;
-        Poly dvs;
-        if(quad_root) {
-            dvs.push_back((sol*conj(sol)).real());
-            dvs.push_back(-(sol + conj(sol)).real());
-            dvs.push_back(1.0);
-            //std::cout << "(" <<  dvs << ")";
-	    //solutions.push_back(sol);
-	    //solutions.push_back(conj(sol));
-        } else {
-            //std::cout << sol << std::endl;
-            dvs.push_back(-sol.real());
-            dvs.push_back(1.0);
-	    solutions.push_back(sol);
-            //std::cout << "(" <<  dvs << ")";
-        }
-        Poly r;
-        p = divide(p, dvs, r);
-        //std::cout << r << std::endl;
-    }
-    return solutions;
-}
-
-std::vector<double> 
-laguerre_real_interval(Poly const & ply, 
-		       const double lo, const double hi,
-		       const double tol) {
-}
-
-/*
-  Local Variables:
-  mode:c++
-  c-file-style:"stroustrup"
-  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
-  indent-tabs-mode:nil
-  fill-column:99
-  End:
-*/
-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
+#include "poly-laguerre-solve.h"

+#include <iterator>

+

+typedef std::complex<double> cdouble;

+

+cdouble laguerre_internal_complex(Poly const & p,

+                                  double x0,

+                                  double tol,

+                                  bool & quad_root) {

+    cdouble a = 2*tol;

+    cdouble xk = x0;

+    double n = p.degree();

+    quad_root = false;

+    const unsigned shuffle_rate = 10;

+    static double shuffle[] = {0, 0.5, 0.25, 0.75, 0.125, 0.375, 0.625, 0.875, 1.0};

+    unsigned shuffle_counter = 0;

+    while(std::norm(a) > (tol*tol)) {

+        //std::cout << "xk = " << xk << std::endl;

+        cdouble b = p.back();

+        cdouble d = 0, f = 0;

+        double err = abs(b);

+        double abx = abs(xk);

+        for(int j = p.size()-2; j >= 0; j--) {

+            f = xk*f + d;

+            d = xk*d + b;

+            b = xk*b + p[j];

+            err = abs(b) + abx*err;

+        }

+        

+        err *= 1e-7; // magic epsilon for convergence, should be computed from tol

+        

+        cdouble px = b;

+        if(abs(b) < err)

+            return xk;

+        //if(std::norm(px) < tol*tol)

+        //    return xk;

+        cdouble G = d / px;

+        cdouble H = G*G - f / px;

+        

+        //std::cout << "G = " << G << "H = " << H;

+        cdouble radicand = (n - 1)*(n*H-G*G);

+        //assert(radicand.real() > 0);

+        if(radicand.real() < 0)

+            quad_root = true;

+        //std::cout << "radicand = " << radicand << std::endl;

+        if(G.real() < 0) // here we try to maximise the denominator avoiding cancellation

+            a = - sqrt(radicand);

+        else

+            a = sqrt(radicand);

+        //std::cout << "a = " << a << std::endl;

+        a = n / (a + G);

+        //std::cout << "a = " << a << std::endl;

+        if(shuffle_counter % shuffle_rate == 0)

+            ;//a *= shuffle[shuffle_counter / shuffle_rate];

+        xk -= a;

+        shuffle_counter++;

+        if(shuffle_counter >= 90)

+            break;

+    }

+    //std::cout << "xk = " << xk << std::endl;

+    return xk;

+}

+

+double laguerre_internal(Poly const & p,

+                         Poly const & pp,

+                         Poly const & ppp,

+                         double x0,

+                         double tol,

+                         bool & quad_root) {

+    double a = 2*tol;

+    double xk = x0;

+    double n = p.degree();

+    quad_root = false;

+    while(a*a > (tol*tol)) {

+        //std::cout << "xk = " << xk << std::endl;

+        double px = p(xk);

+        if(px*px < tol*tol)

+            return xk;

+        double G = pp(xk) / px;

+        double H = G*G - ppp(xk) / px;

+        

+        //std::cout << "G = " << G << "H = " << H;

+        double radicand = (n - 1)*(n*H-G*G);

+        assert(radicand > 0);

+        //std::cout << "radicand = " << radicand << std::endl;

+        if(G < 0) // here we try to maximise the denominator avoiding cancellation

+            a = - sqrt(radicand);

+        else

+            a = sqrt(radicand);

+        //std::cout << "a = " << a << std::endl;

+        a = n / (a + G);

+        //std::cout << "a = " << a << std::endl;

+        xk -= a;

+    }

+    //std::cout << "xk = " << xk << std::endl;

+    return xk;

+}

+

+

+std::vector<cdouble > 

+laguerre(Poly p, const double tol) {

+    std::vector<cdouble > solutions;

+    //std::cout << "p = " << p << " = ";

+    while(p.size() > 1)

+    {

+        double x0 = 0;

+        bool quad_root = false;

+        cdouble sol = laguerre_internal_complex(p, x0, tol, quad_root);

+        //if(abs(sol) > 1) break;

+        Poly dvs;

+        if(quad_root) {

+            dvs.push_back((sol*conj(sol)).real());

+            dvs.push_back(-(sol + conj(sol)).real());

+            dvs.push_back(1.0);

+            //std::cout << "(" <<  dvs << ")";

+	    //solutions.push_back(sol);

+	    //solutions.push_back(conj(sol));

+        } else {

+            //std::cout << sol << std::endl;

+            dvs.push_back(-sol.real());

+            dvs.push_back(1.0);

+	    solutions.push_back(sol);

+            //std::cout << "(" <<  dvs << ")";

+        }

+        Poly r;

+        p = divide(p, dvs, r);

+        //std::cout << r << std::endl;

+    }

+    return solutions;

+}

+

+std::vector<double> 

+laguerre_real_interval(Poly const & ply, 

+		       const double lo, const double hi,

+		       const double tol) {

+}

+

+/*

+  Local Variables:

+  mode:c++

+  c-file-style:"stroustrup"

+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))

+  indent-tabs-mode:nil

+  fill-column:99

+  End:

+*/

+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :

diff --git a/src/2geom/poly.cpp b/src/2geom/poly.cpp
index 5c3a30002..17f286494 100644
--- a/src/2geom/poly.cpp
+++ b/src/2geom/poly.cpp
@@ -112,8 +112,8 @@ Poly derivative(Poly const & p) {
 Poly compose(Poly const & a, Poly const & b) {
     Poly result;
     
-    for(unsigned i = a.size()-1; i >=0; i--) {
-        result = Poly(a[i]) + result * b;
+    for(unsigned i = a.size(); i > 0; i--) {
+        result = Poly(a[i-1]) + result * b;
     }
     return result;
     
diff --git a/src/2geom/rect.h b/src/2geom/rect.h
index d96802014..f0811baf6 100644
--- a/src/2geom/rect.h
+++ b/src/2geom/rect.h
@@ -1,91 +1,153 @@
-//D2<Interval> specialization to Rect:

-

- /* Authors of original rect class:

- *   Lauris Kaplinski <lauris@kaplinski.com>

- *   Nathan Hurst <njh@mail.csse.monash.edu.au>

- *   bulia byak <buliabyak@users.sf.net>

- *   MenTaLguY <mental@rydia.net>

- */

+/*
+ * rect.h - D2<Interval> specialization to Rect
+ *
+ * Copyright 2007 Michael Sloan <mgsloan@gmail.com>
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, output to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ *
+ */
+
+/* Authors of original rect class:
+ *   Lauris Kaplinski <lauris@kaplinski.com>
+ *   Nathan Hurst <njh@mail.csse.monash.edu.au>
+ *   bulia byak <buliabyak@users.sf.net>
+ *   MenTaLguY <mental@rydia.net>
+ */
 
 #ifdef _2GEOM_D2  /*This is intentional: we don't actually want anyone to
-                    include this, other than D2.h */
+                    include this, other than D2.h.  If somone else tries, D2
+                    won't be defined.  If it is, this will already be included. */
 #ifndef _2GEOM_RECT
 #define _2GEOM_RECT
-

-#include "matrix.h"

+
+#include "matrix.h"
 #include <boost/optional/optional.hpp>
-

-namespace Geom {

-

-typedef D2<Interval> Rect;

+
+namespace Geom {
+
+typedef D2<Interval> Rect;
 
 Rect unify(const Rect &, const Rect &);
-

-template<>

-class D2<Interval> {

-  private:

-    Interval f[2];

-    D2<Interval>();// { f[X] = f[Y] = Interval(0, 0); }

-

-  public:

-    D2<Interval>(Interval const &a, Interval const &b) {

-        f[X] = a;

-        f[Y] = b;

-    }

-

-    D2<Interval>(Point const & a, Point const & b) {

-        f[X] = Interval(a[X], b[X]);

-        f[Y] = Interval(a[Y], b[Y]);

-    }

-

-    inline Interval& operator[](unsigned i)              { return f[i]; }

-    inline Interval const & operator[](unsigned i) const { return f[i]; }

-

-    inline Point min() const { return Point(f[X].min(), f[Y].min()); }

-    inline Point max() const { return Point(f[X].max(), f[Y].max()); }

-

-    /** returns the four corners of the rectangle in order

-     *  (clockwise if +Y is up, anticlockwise if +Y is down) */

-    Point corner(unsigned i) const {

-      switch(i % 4) {

-        case 0: return Point(f[X].min(), f[Y].min());

-        case 1: return Point(f[X].max(), f[Y].min());

-        case 2: return Point(f[X].max(), f[Y].max());
        case 3: return Point(f[X].min(), f[Y].max());

-      }

-    }

 
+template<>
+class D2<Interval> {
+  private:
+    Interval f[2];
+    D2<Interval>();// { f[X] = f[Y] = Interval(0, 0); }
+
+  public:
+    D2<Interval>(Interval const &a, Interval const &b) {
+        f[X] = a;
+        f[Y] = b;
+    }
+
+    D2<Interval>(Point const & a, Point const & b) {
+        f[X] = Interval(a[X], b[X]);
+        f[Y] = Interval(a[Y], b[Y]);
+    }
+
+    inline Interval& operator[](unsigned i)              { return f[i]; }
+    inline Interval const & operator[](unsigned i) const { return f[i]; }
+
+    inline Point min() const { return Point(f[X].min(), f[Y].min()); }
+    inline Point max() const { return Point(f[X].max(), f[Y].max()); }
+
+    /** returns the four corners of the rectangle in positive order
+     *  (clockwise if +Y is up, anticlockwise if +Y is down) */
+    Point corner(unsigned i) const {
+        switch(i % 4) {
+            case 0: return Point(f[X].min(), f[Y].min());
+            case 1: return Point(f[X].max(), f[Y].min());
+            case 2: return Point(f[X].max(), f[Y].max());
+            case 3: return Point(f[X].min(), f[Y].max());
+        }
+    }
+        
+    //We should probably remove these - they're coord sys gnostic
+    inline double top() const { return f[Y].min(); }
+    inline double bottom() const { return f[Y].max(); }
+    inline double left() const { return f[X].min(); }
+    inline double right() const { return f[X].max(); }
+    
     inline double width() const { return f[X].extent(); }
     inline double height() const { return f[Y].extent(); }
-

-    /** returns a vector from min to max. */

-    inline Point dimensions() const { return Point(f[X].extent(), f[Y].extent()); }

-    inline Point midpoint() const { return Point(f[X].middle(), f[Y].middle()); }

-

-    inline double area() const { return f[X].extent() * f[Y].extent(); }

-    inline double maxExtent() const { return std::max(f[X].extent(), f[Y].extent()); }

-

-    inline bool isEmpty()                 const { return f[X].isEmpty()        && f[Y].isEmpty(); }

-    inline bool intersects(Rect const &r) const { return f[X].intersects(r[X]) && f[Y].intersects(r[Y]); }

-    inline bool contains(Rect const &r)   const { return f[X].contains(r[X])   && f[Y].contains(r[Y]); }

-    inline bool contains(Point const &p)  const { return f[X].contains(p[X])   && f[Y].contains(p[Y]); }

-

-    inline void expandTo(Point p)        { f[X].extendTo(p[X]);  f[Y].extendTo(p[Y]); }

-    inline void unionWith(Rect const &b) { f[X].unionWith(b[X]); f[Y].unionWith(b[Y]); }

-

-    inline void expandBy(double amnt)    { f[X].expandBy(amnt);  f[Y].expandBy(amnt); }

-    inline void expandBy(Point const p)  { f[X].expandBy(p[X]);  f[Y].expandBy(p[Y]); }

-

+
+    /** returns a vector from min to max. */
+    inline Point dimensions() const { return Point(f[X].extent(), f[Y].extent()); }
+    inline Point midpoint() const { return Point(f[X].middle(), f[Y].middle()); }
+
+    inline double area() const { return f[X].extent() * f[Y].extent(); }
+    inline double maxExtent() const { return std::max(f[X].extent(), f[Y].extent()); }
+
+    inline bool isEmpty()                 const { 
+	return f[X].isEmpty()        && f[Y].isEmpty(); 
+    }
+    inline bool intersects(Rect const &r) const { 
+	return f[X].intersects(r[X]) && f[Y].intersects(r[Y]); 
+    }
+    inline bool contains(Rect const &r)   const { 
+	return f[X].contains(r[X]) && f[Y].contains(r[Y]); 
+    }
+    inline bool contains(Point const &p)  const {
+	return f[X].contains(p[X]) && f[Y].contains(p[Y]);
+    }
+
+    inline void expandTo(Point p)        { 
+	f[X].extendTo(p[X]);  f[Y].extendTo(p[Y]); 
+    }
+    inline void unionWith(Rect const &b) { 
+	f[X].unionWith(b[X]); f[Y].unionWith(b[Y]); 
+    }
+
+    inline void expandBy(double amnt)    { 
+	f[X].expandBy(amnt);  f[Y].expandBy(amnt); 
+    }
+    inline void expandBy(Point const p)  { 
+	f[X].expandBy(p[X]);  f[Y].expandBy(p[Y]); 
+    }
+
     /** Transforms the rect by m. Note that it gives correct results only for scales and translates,
-        in the case of rotations, the area of the rect will grow as it cannot rotate. */

-    inline Rect operator*(Matrix const m) const { return unify(Rect(corner(0) * m, corner(2) * m),
-                                                               Rect(corner(1) * m, corner(3) * m)); }

+        in the case of rotations, the area of the rect will grow as it cannot rotate. */
+    inline Rect operator*(Matrix const m) const { 
+        return unify(Rect(corner(0) * m, corner(2) * m),
+                     Rect(corner(1) * m, corner(3) * m));
+    }
 };
 
-inline Rect unify(const Rect & a, const Rect & b) {
+inline Rect unify(Rect const & a, Rect const & b) {
     return Rect(unify(a[X], b[X]), unify(a[Y], b[Y]));
 }
 
-inline boost::optional<Rect> intersect(const Rect & a, const Rect & b) {
+inline Rect union_list(std::vector<Rect> const &r) {
+    if(r.empty()) return Rect(Interval(0,0), Interval(0,0));
+    Rect ret = r[0];
+    for(unsigned i = 1; i < r.size(); i++)
+        ret.unionWith(r[i]);
+    return ret;
+}
+
+inline boost::optional<Rect> intersect(Rect const & a, Rect const & b) {
     boost::optional<Interval> x = intersect(a[X], b[X]);
     boost::optional<Interval> y = intersect(a[Y], b[Y]);
     return x && y ? boost::optional<Rect>(Rect(*x, *y)) : boost::optional<Rect>();
@@ -95,3 +157,13 @@ inline boost::optional<Rect> intersect(const Rect & a, const Rect & b) {
 
 #endif //_2GEOM_RECT
 #endif //_2GEOM_D2
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
diff --git a/src/2geom/region.cpp b/src/2geom/region.cpp
new file mode 100644
index 000000000..cfab3a35a
--- /dev/null
+++ b/src/2geom/region.cpp
@@ -0,0 +1,42 @@
+#include "region.h"
+#include "utils.h"
+
+#include "shape.h"
+
+namespace Geom {
+
+Regions sanitize_path(Path const &p) {
+    Regions results;
+    Crossings crs = self_crossings(p);
+    for(unsigned i = 0; i < crs.size(); i++) {
+        
+    }
+}
+
+Region Region::operator*(Matrix const &m) const {
+    Region r((m.flips() ? boundary.reverse() : boundary) * m, fill);
+    if(box && m.onlyScaleAndTranslation()) r.box = (*box) * m;
+    return r;
+}
+
+bool Region::invariants() const {
+    return self_crossings(boundary).empty();
+}
+
+unsigned outer_index(Regions const &ps) {
+    if(ps.size() <= 1 || ps[0].contains(ps[1])) {
+        return 0;
+    } else {
+        /* Since we've already shown that chunks[0] is not outside
+           it can be used as an exemplar inner. */
+        Point exemplar = Path(ps[0]).initialPoint();
+        for(unsigned i = 1; i < ps.size(); i++) {
+            if(ps[i].contains(exemplar)) {
+                return i;
+            }
+        }
+    }
+    return ps.size();
+}
+
+}
diff --git a/src/2geom/region.h b/src/2geom/region.h
new file mode 100644
index 000000000..4b434f1e5
--- /dev/null
+++ b/src/2geom/region.h
@@ -0,0 +1,80 @@
+#ifndef __2GEOM_REGION_H
+#define __2GEOM_REGION_H
+
+#include "path.h"
+#include "path-intersection.h"
+
+namespace Geom {
+
+class Shape;
+
+class Region {
+    friend Crossings crossings(Region const &a, Region const &b);
+    friend class Shape;
+    friend Shape shape_boolean(bool rev, Shape const & a, Shape const & b, CrossingSet const & crs);
+
+    Path boundary;
+    mutable boost::optional<Rect> box;
+    bool fill;
+  public:
+    Region() : fill(true) {}
+    explicit Region(Path const &p) : boundary(p) { fill = path_direction(p); }
+    Region(Path const &p, bool dir) : boundary(p), fill(dir) {}
+    Region(Path const &p, boost::optional<Rect> const &b) : boundary(p), box(b) { fill = path_direction(p); }
+    Region(Path const &p, boost::optional<Rect> const &b, bool dir) : boundary(p), box(b), fill(dir) {}
+    
+    unsigned size() const { return boundary.size(); }
+    
+    bool isFill() const { return fill; }    
+    Region asFill() const { if(fill) return Region(*this); else return inverse(); } 
+    Region asHole() const { if(fill) return inverse(); else return Region(*this); }
+    
+    operator Path() const { return boundary; }
+    Rect boundsFast() const {
+        if(!box) box = boost::optional<Rect>(boundary.boundsFast());
+        return *box;
+    }
+    bool contains(Point const &p) const {
+        if(box && !box->contains(p)) return false;
+        return Geom::contains(boundary, p);
+    }
+    bool contains(Region const &other) const { return contains(other.boundary.initialPoint()); }
+    
+    Region inverse() const { return Region(boundary.reverse(), box, !fill); }
+    
+    Region operator*(Matrix const &m) const;
+    
+    bool invariants() const;
+};
+
+typedef std::vector<Region> Regions;
+
+unsigned outer_index(Regions const &ps);
+
+//assumes they're already sanitized somewhat
+inline Regions regions_from_paths(std::vector<Path> const &ps) {
+    Regions res;
+    for(unsigned i = 0; i < ps.size(); i++)
+        res.push_back(Region(ps[i]));
+    return res;
+}
+
+inline std::vector<Path> paths_from_regions(Regions const &rs) {
+    std::vector<Path> res;
+    for(unsigned i = 0; i < rs.size(); i++)
+        res.push_back(rs[i]);
+    return res;
+}
+
+Regions sanitize_path(Path const &p);
+
+Regions region_boolean(bool rev, Region const & a, Region const & b, Crossings const &cr);
+Regions region_boolean(bool rev, Region const & a, Region const & b, Crossings const & cr_a, Crossings const & cr_b);
+
+inline Regions region_boolean(bool rev, Region const & a, Region const & b) {
+    return region_boolean(rev, a, b, crossings(a, b));
+}
+
+}
+
+#endif
diff --git a/src/2geom/sbasis-geometric.cpp b/src/2geom/sbasis-geometric.cpp
index 170323c6c..fde907857 100644
--- a/src/2geom/sbasis-geometric.cpp
+++ b/src/2geom/sbasis-geometric.cpp
@@ -350,7 +350,7 @@ unsigned Geom::centroid(Piecewise<D2<SBasis> > const &p, Point& centroid, double
     }
 // join ends
     centroid_tmp *= 2;
-    Point final = p[p.size()].at1(), initial = p[0].at0();
+    Point final = p[p.size()-1].at1(), initial = p[0].at0();
     const double ai = cross(final, initial);
     atmp += ai;
     centroid_tmp += (final + initial)*ai; // first moment.
diff --git a/src/2geom/sbasis-math.cpp b/src/2geom/sbasis-math.cpp
index d88c26832..db4cf5e08 100644
--- a/src/2geom/sbasis-math.cpp
+++ b/src/2geom/sbasis-math.cpp
@@ -231,11 +231,11 @@ Piecewise<SBasis> reciprocalOnDomain(Interval range, double tol){
 
     if (a<=tol){
         reciprocal_fn.push_cut(0);
-        int i0=(int) floor(log(tol)/log(R));
+        int i0=(int) floor(std::log(tol)/std::log(R));
         a=pow(R,i0);
         reciprocal_fn.push(Linear(1/a),a);
     }else{
-        int i0=(int) floor(log(a)/log(R));
+        int i0=(int) floor(std::log(a)/std::log(R));
         a=pow(R,i0);
         reciprocal_fn.cuts.push_back(a);
     }  
diff --git a/src/2geom/sbasis-math.h b/src/2geom/sbasis-math.h
index 529641068..8f4e1612d 100644
--- a/src/2geom/sbasis-math.h
+++ b/src/2geom/sbasis-math.h
@@ -70,6 +70,9 @@ Piecewise<SBasis> cos(          SBasis  const &f, double tol=1e-3, int order=3);
 Piecewise<SBasis> cos(Piecewise<SBasis> const &f, double tol=1e-3, int order=3);
 Piecewise<SBasis> sin(          SBasis  const &f, double tol=1e-3, int order=3);
 Piecewise<SBasis> sin(Piecewise<SBasis> const &f, double tol=1e-3, int order=3);
+//-Log---------------------------------------------------------------
+Piecewise<SBasis> log(          SBasis const &f, double tol=1e-3, int order=3);
+Piecewise<SBasis> log(Piecewise<SBasis>const &f, double tol=1e-3, int order=3);
 
 //--1/x------------------------------------------------------------
 //TODO: change this...
diff --git a/src/2geom/sbasis-poly.cpp b/src/2geom/sbasis-poly.cpp
index e0fa828f9..ff797920f 100644
--- a/src/2geom/sbasis-poly.cpp
+++ b/src/2geom/sbasis-poly.cpp
@@ -15,6 +15,8 @@ SBasis poly_to_sbasis(Poly const & p) {
 }
 
 Poly sbasis_to_poly(SBasis const & sb) {
+    if(sb.isZero())
+        return Poly();
     Poly S; // (1-x)x = -1*x^2 + 1*x + 0
     Poly A, B;
     B.push_back(0);
diff --git a/src/2geom/sbasis-roots.cpp b/src/2geom/sbasis-roots.cpp
index ad2dfbda4..c4ef3c16d 100644
--- a/src/2geom/sbasis-roots.cpp
+++ b/src/2geom/sbasis-roots.cpp
@@ -68,7 +68,7 @@ Interval bounds_fast(const SBasis &sb, int order) {
         double a=sb[j][0];
         double b=sb[j][1];
 
-        double t, v;
+        double v, t = 0;
         v = res[0];
         if (v<0) t = ((b-a)/v+1)*0.5;
         if (v>=0 || t<0 || t>1) {
@@ -95,7 +95,7 @@ Interval bounds_local(const SBasis &sb, const Interval &i, int order) {
         double a=sb[j][0];
         double b=sb[j][1];
 
-        double t;
+        double t = 0;
         if (lo<0) t = ((b-a)/lo+1)*0.5;
         if (lo>=0 || t<t0 || t>t1) {
             lo = std::min(a*(1-t0)+b*t0+lo*t0*(1-t0),a*(1-t1)+b*t1+lo*t1*(1-t1));
diff --git a/src/2geom/sbasis-to-bezier.cpp b/src/2geom/sbasis-to-bezier.cpp
index 613c39288..206f18931 100644
--- a/src/2geom/sbasis-to-bezier.cpp
+++ b/src/2geom/sbasis-to-bezier.cpp
@@ -184,34 +184,6 @@ path_from_sbasis(D2<SBasis> const &B, double tol) {
 //TODO: some of this logic should be lifted into svg-path
 std::vector<Geom::Path>
 path_from_piecewise(Geom::Piecewise<Geom::D2<Geom::SBasis> > const &B, double tol) {
-/*
-    Geom::PathBuilder pb;
-    if(B.size() == 0) return pb.peek();
-    Geom::Point start = B[0].at0();
-    pb.moveTo(start);
-    for(unsigned i = 0; ; i++) {
-        if(i+1 == B.size() || !near(B[i+1].at0(), B[i].at1(), tol)) {
-            //start of a new path
-            if(near(start, B[i].at1())) {
-                //it's closed
-                pb.closePath();
-                if(sbasis_size(B[i]) <= 1) {
-                    //last line seg already there
-                    goto no_add;
-                }
-            }
-            build_from_sbasis(pb, B[i], tol);
-          no_add:
-            if(i+1 >= B.size()) break;
-            start = B[i+1].at0();
-            pb.moveTo(start);
-            i++;
-        }
-        build_from_sbasis(pb, B[i], tol);
-    }
-    pb.finish();
-    return pb.peek();
-*/
     Geom::PathBuilder pb;
     if(B.size() == 0) return pb.peek();
     Geom::Point start = B[0].at0();
diff --git a/src/2geom/sbasis.h b/src/2geom/sbasis.h
index a27d4bfce..b0df4435c 100644
--- a/src/2geom/sbasis.h
+++ b/src/2geom/sbasis.h
@@ -39,6 +39,7 @@
 
 #include "linear.h"
 #include "interval.h"
+#include "utils.h"
 
 namespace Geom {
 
@@ -88,6 +89,12 @@ public:
     double operator()(double t) const {
         return valueAt(t);
     }
+
+    std::vector<double> valueAndDerivatives(double t, unsigned n) const {
+        //TODO
+        throw NotImplemented();
+    }
+
     SBasis toSBasis() const { return SBasis(*this); }
 
     double tailError(unsigned tail) const;
diff --git a/src/2geom/shape.cpp b/src/2geom/shape.cpp
new file mode 100644
index 000000000..670792521
--- /dev/null
+++ b/src/2geom/shape.cpp
@@ -0,0 +1,442 @@
+#include "shape.h"
+#include "utils.h"
+#include "sweep.h"
+
+#include <iostream>
+#include <algorithm>
+
+namespace Geom {
+
+// Utility funcs
+
+// Yes, xor is !=, but I'm pretty sure this is safer in the event of strange bools
+bool logical_xor (bool a, bool b) { return (a || b) && !(a && b); }
+
+// A little sugar for appending a list to another
+template<typename T>
+void append(T &a, T const &b) {
+    a.insert(a.end(), b.begin(), b.end());
+}
+
+/* Used within shape_boolean and related functions, as the name describes, finds the
+ * first false within the list of lists of booleans.
+ */
+void first_false(std::vector<std::vector<bool> > visited, unsigned &i, unsigned &j) {
+    for(i = 0, j = 0; i < visited.size(); i++) {
+        std::vector<bool>::iterator unvisited = std::find(visited[i].begin(), visited[i].end(), false);
+        if(unvisited != visited[i].end()) {
+            j = unvisited - visited[i].begin();
+            break;
+        }
+    }
+}
+
+// Finds a crossing in a list of them, given the sorting index.
+unsigned find_crossing(Crossings const &cr, Crossing x, unsigned i) {
+    return std::lower_bound(cr.begin(), cr.end(), x, CrossingOrder(i)) - cr.begin();
+}
+
+/* This function handles boolean ops on shapes.  The first parameter is a bool
+ * which determines its behavior in each combination of cases.  For proper
+ * fill information and noncrossing behavior, the fill data of the regions
+ * must be correct.  The boolean parameter determines whether the operation
+ * is a union or a subtraction.  Reversed paths represent inverse regions,
+ * where everything is included in the fill except for the insides.
+ *
+ * Here is a chart of the behavior under various circumstances:
+ * 
+ * rev = false (union)
+ *            A
+ *        F         H
+ * F  A+B -> F  A-B -> H
+ *B
+ * H  B-A -> H  AxB -> H
+ *
+ * rev = true (intersect)
+ *            A
+ *        F         H
+ * F  AxB -> F  B-A -> F
+ *B
+ * H  A-B -> F  A+B -> H
+ *
+ * F/H = Fill outer / Hole outer
+ * A/B specify operands
+ * + = union, - = subtraction, x = intersection
+ * -> read as "produces"
+ *
+ * This is the main function of boolops, yet its operation isn't very complicated.
+ * It traverses the crossings, and uses the crossing direction to decide whether
+ * the next segment should be taken from A or from B.  The second half of the
+ * function deals with figuring out what to do with bits that have no intersection.
+ */
+Shape shape_boolean(bool rev, Shape const & a, Shape const & b, CrossingSet const & crs) {
+    const Regions ac = a.content, bc = b.content;
+
+    //Keep track of which crossings we've hit.
+    std::vector<std::vector<bool> > visited;
+    for(unsigned i = 0; i < crs.size(); i++)
+        visited.push_back(std::vector<bool>(crs[i].size(), false));
+    
+    //bool const exception = 
+    
+    //Traverse the crossings, creating chunks
+    Regions chunks;
+    while(true) {
+        unsigned i, j;
+        first_false(visited, i, j);
+        if(i == visited.size()) break;
+        
+        Path res;
+        do {
+            Crossing cur = crs[i][j];
+            visited[i][j] = true;
+            
+            //get indices of the dual:
+            unsigned io = cur.getOther(i), jo = find_crossing(crs[io], cur, io);
+            if(jo < visited[io].size()) visited[io][jo] = true;
+            
+            //main driving logic
+            if(logical_xor(cur.dir, rev)) {
+                if(i >= ac.size()) { i = io; j = jo; }
+                j++;
+                if(j >= crs[i].size()) j = 0;
+                Crossing next = crs[i][j];
+                ac[next.a].boundary.appendPortionTo(res, cur.ta, next.ta);
+            } else {
+                if(i < ac.size()) { i = io; j = jo; }
+                j++;
+                if(j >= crs[i].size()) j = 0;
+                Crossing next = crs[i][j];
+                bc[next.b - ac.size()].boundary.appendPortionTo(res, cur.tb, next.tb);
+            }
+        } while (!visited[i][j]);
+        if(res.size() > 0) chunks.push_back(Region(res));
+    }
+    
+    //If true, then we are on the 'subtraction diagonal'
+    bool const on_sub = logical_xor(a.fill, b.fill);
+    //If true, then the hole must be inside the other to be included
+    bool const a_mode = logical_xor(logical_xor(!rev, a.fill), on_sub),
+               b_mode = logical_xor(logical_xor(!rev, b.fill), on_sub);
+    
+    //Handle unintersecting portions
+    for(unsigned i = 0; i < crs.size(); i++) {
+        if(crs[i].size() == 0) {
+            Region    r(i < ac.size() ? ac[i] : bc[i - ac.size()]);
+            bool   mode(i < ac.size() ? a_mode : b_mode);
+            
+            if(logical_xor(r.fill, i < ac.size() ? a.fill : b.fill)) {
+                //is an inner (fill is opposite the outside fill)
+                Point exemplar = r.boundary.initialPoint();
+                Regions const & others = i < ac.size() ? bc : ac;
+                for(unsigned j = 0; j < others.size(); j++) {
+                    if(others[j].contains(exemplar)) {
+                        //contained in another
+                        if(mode) chunks.push_back(r);
+                        goto skip;
+                    }
+                }
+            }
+            //disjoint
+            if(!mode) chunks.push_back(r);
+            skip: (void)0;
+        }
+    }
+    
+    return Shape(chunks);
+}
+
+// Just a convenience wrapper for shape_boolean, which handles the crossings
+Shape shape_boolean(bool rev, Shape const & a, Shape const & b) {
+    CrossingSet crs = crossings_between(a, b);
+    
+    return shape_boolean(rev, a, b, crs);
+}
+
+
+// Some utility functions for boolop:
+
+std::vector<double> region_sizes(Shape const &a) {
+    std::vector<double> ret;
+    for(unsigned i = 0; i < a.size(); i++) {
+        ret.push_back(double(a[i].size()));
+    }
+    return ret;
+}
+
+Shape shape_boolean_ra(bool rev, Shape const &a, Shape const &b, CrossingSet const &crs) {
+    return shape_boolean(rev, a.inverse(), b, reverse_ta(crs, a.size(), region_sizes(a)));
+}
+
+Shape shape_boolean_rb(bool rev, Shape const &a, Shape const &b, CrossingSet const &crs) {
+    return shape_boolean(rev, a, b.inverse(), reverse_tb(crs, a.size(), region_sizes(b)));
+}
+
+/* This is a function based on shape_boolean which allows boolean operations
+ * to be specified as a logic table.  This logic table is 4 bit-flags, which
+ * correspond to the elements of the 'truth table' for a particular operation.
+ * These flags are defined with the enums starting with BOOLOP_ .
+ */
+Shape boolop(Shape const &a, Shape const &b, unsigned flags, CrossingSet const &crs) {
+    flags &= 15;
+    if(flags <= BOOLOP_UNION) {
+        switch(flags) {
+            case BOOLOP_INTERSECT:    return shape_boolean(true, a, b, crs);
+            case BOOLOP_SUBTRACT_A_B: return shape_boolean_rb(true, a, b, crs);
+            case BOOLOP_IDENTITY_A:   return a;
+            case BOOLOP_SUBTRACT_B_A: return shape_boolean_ra(true, a, b, crs);
+            case BOOLOP_IDENTITY_B:   return b;
+            case BOOLOP_EXCLUSION: {
+                Shape res = shape_boolean_rb(true, a, b, crs);
+                append(res.content, shape_boolean_ra(true, a, b, crs).content);
+                return res;
+            }
+            case BOOLOP_UNION:        return shape_boolean(false, a, b);
+        }
+    } else {
+        switch(flags - BOOLOP_NEITHER) {
+            case BOOLOP_SUBTRACT_A_B: return shape_boolean_ra(false, a, b, crs);
+            case BOOLOP_SUBTRACT_B_A: return shape_boolean_rb(false, a, b, crs);
+            case BOOLOP_EXCLUSION: {
+                Shape res = shape_boolean_ra(false, a, b, CrossingSet(crs));
+                append(res.content, shape_boolean_rb(false, a, b, CrossingSet(crs)).content);
+                return res;
+            }
+        }
+        return boolop(a, b, ~flags, crs).inverse();
+    }
+    return Shape();
+}
+
+/* This version of the boolop function doesn't require a set of crossings, as
+ * it computes them for you.  This is more efficient in some cases, as the
+ * shape can be inverted before finding crossings.  In the special case of
+ * exclusion it uses the other version of boolop.
+ */
+Shape boolop(Shape const &a, Shape const &b, unsigned flags) {
+    flags &= 15;
+    if(flags <= BOOLOP_UNION) {
+        switch(flags) {
+            case BOOLOP_INTERSECT:    return shape_boolean(true, a, b);
+            case BOOLOP_SUBTRACT_A_B: return shape_boolean(true, a, b.inverse());
+            case BOOLOP_IDENTITY_A:   return a;
+            case BOOLOP_SUBTRACT_B_A: return shape_boolean(true, b, a.inverse());
+            case BOOLOP_IDENTITY_B:   return b;
+            case BOOLOP_EXCLUSION: {
+                Shape res = shape_boolean(true, a, b.inverse());
+                append(res.content, shape_boolean(true, b, a.inverse()).content);
+                return res;
+            } //return boolop(a, b, flags,  crossings_between(a, b));
+            case BOOLOP_UNION:        return shape_boolean(false, a, b);
+        }
+    } else {
+        switch(flags - BOOLOP_NEITHER) {
+            case BOOLOP_SUBTRACT_A_B: return shape_boolean(false, b, a.inverse());
+            case BOOLOP_SUBTRACT_B_A: return shape_boolean(false, a, b.inverse());
+            case BOOLOP_EXCLUSION: {
+                Shape res = shape_boolean(false, a, b.inverse());
+                append(res.content, shape_boolean(false, b, a.inverse()).content);
+                return res;
+            } //return boolop(a, b, flags, crossings_between(a, b));
+        }
+        return boolop(a, b, ~flags).inverse();
+    }
+    return Shape();
+}
+
+
+int paths_winding(std::vector<Path> const &ps, Point p) {
+    int ret;
+    for(unsigned i = 0; i < ps.size(); i++)
+        ret += winding(ps[i], p);
+    return ret;
+}
+
+std::vector<double> y_of_roots(std::vector<Path> const & ps, double x) {
+    std::vector<double> res;
+    for(unsigned i = 0; i < ps.size(); i++) {
+        std::vector<double> temp = ps[i].roots(x, X);
+        for(unsigned i = 0; i < temp.size(); i++)
+            res.push_back(ps[i].valueAt(temp[i], Y));
+    }
+    std::sort(res.begin(), res.end());
+    return res;
+}
+
+struct Edge {
+    unsigned ix;
+    double from, to;
+    bool rev;
+    int wind;
+    Edge(unsigned i, double ft, double tt, bool r, unsigned w) : ix(i), from(ft), to(tt), rev(r), wind(w) {}
+    Edge(unsigned i, double ft, double tt, bool r, std::vector<Path> const &ps) : ix(i), from(ft), to(tt), rev(r) {
+        //TODO: get the edge wind data some other way
+        Point p = ps[i].pointAt(ft);
+        std::vector<double> rs = y_of_roots(ps, p[X]);
+        unsigned interv = std::lower_bound(rs.begin(), rs.end(), p[Y]) - rs.begin();
+        wind = interv % 2;
+    }
+    double initial() { return rev ? to : from; }
+    double final() { return rev ? from : to; }
+    void addTo(Path &res, std::vector<Path> const &ps) {
+        if(rev) {
+            Path p = ps[ix].portion(to, from).reverse();
+            for(unsigned i = 0; i < p.size(); i++)
+                res.append(p[i]);
+        } else {
+            ps[ix].appendPortionTo(res, from, to);
+        }
+    }
+};
+
+typedef std::vector<Edge> Edges;
+
+double point_cosine(Point a, Point b, Point c) {
+    Point db = b - a, dc = c - a;
+    return cross(db, dc) / (db.length() * dc.length());
+}
+
+//sanitize
+Regions regionize_paths(std::vector<Path> const &ps, bool evenodd) {
+    CrossingSet crs = crossings_among(ps);
+    
+    Edges es;
+    
+    for(unsigned i = 0; i < crs.size(); i++) {
+        for(unsigned j = 0; j < crs[i].size(); j++) {
+            Crossing cur = crs[i][j];
+            int io = i, jo = j;
+            
+            jo++;
+            if(jo >= crs[io].size()) jo = 0;
+            //std::cout << io << ", " << jo << "\n";
+            if(cur.a == i)
+                es.push_back(Edge(i, cur.ta, crs[io][jo].ta, false, ps));
+            else
+                es.push_back(Edge(i, cur.tb, crs[io][jo].tb, false, ps));
+            
+            jo-=2;
+            if(jo < 0) jo += crs[io].size();
+            // std::cout << io << ", " << jo << "\n";
+            if(cur.a == i)
+                es.push_back(Edge(i, crs[io][jo].ta, cur.ta, true, ps));
+            else
+                es.push_back(Edge(i, crs[io][jo].tb, cur.tb, true, ps));
+        }
+    }
+    for(unsigned i = 0; i<crs.size(); i++) {
+        if(crs[i].empty()) {
+            es.push_back(Edge(i, 0, ps[i].size(), false, ps));
+            es.push_back(Edge(i, ps[i].size(), 0, true, ps));
+        }
+    }
+    
+    Edges es2;
+    //filter
+    for(unsigned i = 0; i < es.size(); i++) {
+        if(true) //(evenodd && es[i].wind % 2 == 0) || (!evenodd && es[i].wind == 0))
+            es2.push_back(es[i]);
+    }
+    es = es2;
+    
+    std::cout << es.size() << " edges\n";
+    
+    Regions chunks;
+    for(unsigned i = 0; i < es.size(); i++) {
+        Edge cur = es[i];
+        if(cur.rev)
+            chunks.push_back(Region(ps[cur.ix].portion(cur.from, cur.to).reverse(), cur.wind % 2 != 0));
+        else
+            chunks.push_back(Region(ps[cur.ix].portion(cur.from, cur.to), cur.wind % 2 != 0));
+    }
+    return chunks;
+    
+    //Regions chunks;
+    std::vector<bool> used(es2.size(), false);
+    while(true) {
+        unsigned i = std::find(used.begin(), used.end(), false) - used.begin();
+        if(i == used.size()) break;
+        Path res;
+        do {
+            es2[i].addTo(res, ps);
+            Point pnt = res.finalPoint();
+            std::vector<unsigned> poss;
+            for(unsigned j = 0; j < es2.size(); j++)
+                if(near(pnt, ps[es2[j].ix].pointAt(es2[j].initial()))) poss.push_back(j);
+            if(poss.empty()) break;
+            unsigned best = 0;
+            if(poss.size() > 1) {
+                double crossval = 10;
+                Point along = ps[i].pointAt(es2[i].final()+0.1);
+                for(unsigned j = 0; j < poss.size(); j++) {
+                    unsigned ix = poss[j];
+                    double val = point_cosine(pnt, along, ps[ix].pointAt(es2[ix].initial()+.01));
+                    if(val < crossval) {
+                        crossval = val;
+                        best = j;
+                    }
+                }
+            }
+            i = poss[best];
+        } while(!used[i]);
+        chunks.push_back(Region(res));
+    }
+    return chunks;
+}
+
+/* This transforms a shape by a matrix.  In the case that the matrix flips
+ * the shape, it reverses the paths in order to preserve the fill.
+ */
+Shape Shape::operator*(Matrix const &m) const {
+    Shape ret;
+    for(unsigned i = 0; i < size(); i++)
+        ret.content.push_back(content[i] * m);
+    ret.fill = fill;
+    return ret;
+}
+
+// Inverse is a boolean not, and simply reverses all the paths & fill flags
+Shape Shape::inverse() const {
+    Shape ret;
+    for(unsigned i = 0; i < size(); i++)
+        ret.content.push_back(content[i].inverse());
+    ret.fill = !fill;
+    return ret;
+}
+
+struct ContainmentOrder {
+    std::vector<Region> const *rs;
+    explicit ContainmentOrder(std::vector<Region> const *r) : rs(r) {}
+    bool operator()(unsigned a, unsigned b) const { return (*rs)[b].contains((*rs)[a]); }
+};
+
+bool Shape::contains(Point const &p) const {
+    std::vector<Rect> pnt;
+    pnt.push_back(Rect(p, p));
+    std::vector<std::vector<unsigned> > cull = sweep_bounds(pnt, bounds(*this));
+    if(cull[0].size() == 0) return !fill;
+    return content[*min_element(cull[0].begin(), cull[0].end(), ContainmentOrder(&content))].isFill();
+}
+
+bool Shape::inside_invariants() const {  //semi-slow & easy to violate
+    for(unsigned i = 0; i < size(); i++)
+        if( logical_xor(content[i].isFill(), contains(content[i].boundary.initialPoint())) ) return false;
+    return true;
+}
+bool Shape::region_invariants() const { //semi-slow
+    for(unsigned i = 0; i < size(); i++)
+        if(!content[i].invariants()) return false;
+    return true;
+}
+bool Shape::cross_invariants() const { //slow
+    CrossingSet crs; // = crossings_among(paths_from_regions(content));
+    for(unsigned i = 0; i < crs.size(); i++)
+        if(!crs[i].empty()) return false;
+    return true;
+}
+
+bool Shape::invariants() const {
+    return inside_invariants() && region_invariants() && cross_invariants();
+}
+
+}
diff --git a/src/2geom/shape.h b/src/2geom/shape.h
new file mode 100644
index 000000000..3700e9e5a
--- /dev/null
+++ b/src/2geom/shape.h
@@ -0,0 +1,90 @@
+#ifndef __2GEOM_SHAPE_H
+#define __2GEOM_SHAPE_H
+
+#include <vector>
+#include <set>
+
+#include "region.h"
+
+//TODO: BBOX optimizations
+
+namespace Geom {
+
+enum {
+  BOOLOP_JUST_A  = 1,
+  BOOLOP_JUST_B  = 2,
+  BOOLOP_BOTH    = 4,
+  BOOLOP_NEITHER = 8
+};
+
+enum {
+  BOOLOP_NULL         = 0,
+  BOOLOP_INTERSECT    = BOOLOP_BOTH,
+  BOOLOP_SUBTRACT_A_B = BOOLOP_JUST_B,
+  BOOLOP_IDENTITY_A   = BOOLOP_JUST_A | BOOLOP_BOTH,
+  BOOLOP_SUBTRACT_B_A = BOOLOP_JUST_A,
+  BOOLOP_IDENTITY_B   = BOOLOP_JUST_B | BOOLOP_BOTH,
+  BOOLOP_EXCLUSION    = BOOLOP_JUST_A | BOOLOP_JUST_B,
+  BOOLOP_UNION        = BOOLOP_JUST_A | BOOLOP_JUST_B | BOOLOP_BOTH
+};
+
+class Shape {
+    Regions content;
+    mutable bool fill;
+    //friend Shape shape_region_boolean(bool rev, Shape const & a, Region const & b);
+    friend CrossingSet crossings_between(Shape const &a, Shape const &b);
+    friend Shape shape_boolean(bool rev, Shape const &, Shape const &, CrossingSet const &);
+    friend Shape boolop(Shape const &a, Shape const &b, unsigned);
+    friend Shape boolop(Shape const &a, Shape const &b, unsigned, CrossingSet const &);
+
+  public:
+    Shape() : fill(true) {}
+    explicit Shape(Region const & r) {
+        content = Regions(1, r);
+        fill = r.fill;
+    }
+    explicit Shape(Regions const & r) : content(r) { update_fill(); }
+    explicit Shape(bool f) : fill(f) {}
+    Shape(Regions const & r, bool f) : content(r), fill(f) {}
+    
+    Regions getContent() const { return content; }
+    bool isFill() const { return fill; }
+    
+    unsigned size() const { return content.size(); }
+    const Region &operator[](unsigned ix) const { return content[ix]; }
+    
+    Shape inverse() const;
+    Shape operator*(Matrix const &m) const;
+    
+    bool contains(Point const &p) const;
+    
+    bool inside_invariants() const;  //semi-slow & easy to violate : checks that the insides are inside, the outsides are outside
+    bool region_invariants() const; //semi-slow                    : checks for self crossing
+    bool cross_invariants() const; //slow                          : checks that everything is disjoint
+    bool invariants() const;      //vera slow (combo rombo, checks the above)
+
+  private:     
+    void update_fill() const {
+        unsigned ix = outer_index(content);
+        if(ix < size())
+            fill = content[ix].fill;
+        else if(size() > 0)
+            fill = content.front().fill;
+        else
+            fill = true;
+    }
+};
+
+inline CrossingSet crossings_between(Shape const &a, Shape const &b) { return crossings(paths_from_regions(a.content), paths_from_regions(b.content)); }
+
+Shape shape_boolean(bool rev, Shape const &, Shape const &, CrossingSet const &);
+Shape shape_boolean(bool rev, Shape const &, Shape const &);
+
+Shape boolop(Shape const &, Shape const &, unsigned flags);
+Shape boolop(Shape const &, Shape const &, unsigned flags, CrossingSet &);
+
+Regions regionize_paths(std::vector<Path> const &ps, bool evenodd=true);
+
+}
+
+#endif
diff --git a/src/2geom/solve-bezier-one-d.cpp b/src/2geom/solve-bezier-one-d.cpp
index 558c10c0f..4f37cd049 100644
--- a/src/2geom/solve-bezier-one-d.cpp
+++ b/src/2geom/solve-bezier-one-d.cpp
@@ -35,7 +35,7 @@ const double BEPSILON = ldexp(1.0,-MAXDEPTH-1); /*Flatness control value */
  *    of the roots in the open interval (0, 1).  Return the number of roots found.
  */
 void
-find_bernstein_roots(double *w, /* The control points  */
+find_bernstein_roots(double const *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 */
@@ -43,13 +43,11 @@ find_bernstein_roots(double *w, /* The control points  */
 {  
     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;
@@ -75,18 +73,16 @@ find_bernstein_roots(double *w, /* The control points  */
             const double Ax = right_t - left_t;
             const double Ay = w[degree] - w[0];
 
-            solutions.push_back(left_t + Ax*w[0] / Ay);
+            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  */
+    const double split = 0.5;
     Bernstein(w, degree, split, Left, Right);
     
     double mid_t = left_t*(1-split) + right_t*split;
@@ -101,73 +97,6 @@ find_bernstein_roots(double *w, /* The control points  */
 }
 
 /*
- *  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.
diff --git a/src/2geom/solver.h b/src/2geom/solver.h
index c2e290251..e0d66fc2a 100644
--- a/src/2geom/solver.h
+++ b/src/2geom/solver.h
@@ -5,6 +5,8 @@
 
 namespace Geom{
 
+	class Point;
+
 unsigned
 crossing_count(Geom::Point const *V,	/*  Control pts of Bezier curve	*/
 	       unsigned degree);	/*  Degree of Bezier curve */
@@ -21,14 +23,7 @@ crossing_count(double const *V,	/*  Control pts of Bezier curve	*/
 	       double left_t, double right_t);
 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=0, double right_t=1);
-void
-find_bernstein_roots_buggy(
-    double *w, /* The control points  */
+    double const *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 */
diff --git a/src/2geom/svg-path-parser.cpp b/src/2geom/svg-path-parser.cpp
index 0063fcfe9..0bd15e8b9 100644
--- a/src/2geom/svg-path-parser.cpp
+++ b/src/2geom/svg-path-parser.cpp
@@ -1,4 +1,4 @@
-#line 1 "/home/michael/2geom/src/svg-path-parser.rl"
+#line 1 "/home/njh/svn/lib2geom/src/svg-path-parser.rl"
 /*
  * parse SVG path specifications
  *
@@ -129,7 +129,7 @@ private:
 };
 
 
-#line 133 "/home/michael/2geom/src/svg-path-parser.cpp"
+#line 133 "/home/njh/svn/lib2geom/src/svg-path-parser.cpp"
 static const char _svg_path_actions[] = {
 	0, 1, 0, 1, 1, 1, 2, 1, 
 	3, 1, 4, 1, 5, 1, 15, 1, 
@@ -905,8 +905,8 @@ static const short _svg_path_indicies[] = {
 	185, 187, 188, 189, 190, 191, 192, 193, 
 	194, 195, 196, 197, 198, 199, 200, 201, 
 	202, 203, 194, 179, 180, 205, 0, 575, 
-	575, 576, 576, 577, 575, 578, 0, 761, 
-	761, 762, 762, 763, 761, 764, 0, 781, 
+	575, 576, 576, 577, 575, 578, 0, 769, 
+	769, 770, 770, 771, 769, 772, 0, 781, 
 	781, 782, 782, 783, 781, 784, 0, 750, 
 	741, 0, 714, 0, 699, 699, 701, 702, 
 	715, 715, 699, 700, 714, 0, 738, 738, 
@@ -939,8 +939,8 @@ static const short _svg_path_indicies[] = {
 	294, 295, 295, 297, 320, 300, 301, 303, 
 	304, 305, 306, 307, 308, 309, 310, 311, 
 	312, 313, 314, 315, 316, 317, 318, 319, 
-	310, 295, 296, 321, 0, 757, 757, 758, 
-	758, 759, 757, 760, 0, 777, 777, 778, 
+	310, 295, 296, 321, 0, 765, 765, 766, 
+	766, 767, 765, 768, 0, 777, 777, 778, 
 	778, 779, 777, 780, 0, 749, 740, 0, 
 	712, 0, 694, 694, 696, 697, 713, 713, 
 	694, 695, 712, 0, 737, 737, 728, 730, 
@@ -1054,9 +1054,9 @@ static const short _svg_path_indicies[] = {
 	123, 125, 148, 128, 129, 131, 132, 133, 
 	134, 135, 136, 137, 138, 139, 140, 141, 
 	142, 143, 144, 145, 146, 147, 138, 123, 
-	124, 149, 0, 769, 769, 770, 770, 771, 
-	769, 772, 0, 765, 765, 766, 766, 767, 
-	765, 768, 0, 599, 599, 600, 600, 601, 
+	124, 149, 0, 761, 761, 762, 762, 763, 
+	761, 764, 0, 757, 757, 758, 758, 759, 
+	757, 760, 0, 599, 599, 600, 600, 601, 
 	599, 602, 0, 607, 607, 608, 608, 609, 
 	607, 610, 0, 208, 209, 209, 211, 629, 
 	214, 215, 628, 217, 218, 219, 220, 221, 
@@ -1345,9 +1345,9 @@ static const unsigned char _svg_path_trans_actions_wi[] = {
 	0, 1, 1, 1, 0, 1, 1, 1, 
 	0, 0, 0, 0, 0, 0, 0, 0, 
 	0, 0, 0, 0, 0, 0, 0, 3, 
-	17, 3, 17, 0, 0, 9, 59, 59, 
-	59, 9, 59, 59, 59, 11, 62, 62, 
-	62, 11, 62, 62, 62, 0, 1, 1, 
+	17, 3, 17, 0, 0, 11, 62, 62, 
+	62, 11, 62, 62, 62, 9, 59, 59, 
+	59, 9, 59, 59, 59, 0, 1, 1, 
 	1, 0, 1, 1, 1, 0, 1, 1, 
 	1, 0, 0, 0, 0, 0, 0, 0
 };
@@ -1356,7 +1356,7 @@ static const int svg_path_start = 0;
 
 static const int svg_path_first_final = 326;
 
-#line 133 "/home/michael/2geom/src/svg-path-parser.rl"
+#line 133 "/home/njh/svn/lib2geom/src/svg-path-parser.rl"
 
 
 void Parser::parse(char const *str)
@@ -1369,7 +1369,7 @@ throw(SVGPathParseError)
     _reset();
 
     
-#line 1373 "/home/michael/2geom/src/svg-path-parser.cpp"
+#line 1373 "/home/njh/svn/lib2geom/src/svg-path-parser.cpp"
 	{
 	cs = svg_path_start;
 	}
@@ -1445,13 +1445,13 @@ _match:
 		switch ( *_acts++ )
 		{
 	case 0:
-#line 145 "/home/michael/2geom/src/svg-path-parser.rl"
+#line 145 "/home/njh/svn/lib2geom/src/svg-path-parser.rl"
 	{
             start = p;
         }
 	break;
 	case 1:
-#line 149 "/home/michael/2geom/src/svg-path-parser.rl"
+#line 149 "/home/njh/svn/lib2geom/src/svg-path-parser.rl"
 	{
             char const *end=p;
             std::string buf(start, end);
@@ -1460,55 +1460,55 @@ _match:
         }
 	break;
 	case 2:
-#line 156 "/home/michael/2geom/src/svg-path-parser.rl"
+#line 156 "/home/njh/svn/lib2geom/src/svg-path-parser.rl"
 	{
             _push(1.0);
         }
 	break;
 	case 3:
-#line 160 "/home/michael/2geom/src/svg-path-parser.rl"
+#line 160 "/home/njh/svn/lib2geom/src/svg-path-parser.rl"
 	{
             _push(0.0);
         }
 	break;
 	case 4:
-#line 164 "/home/michael/2geom/src/svg-path-parser.rl"
+#line 164 "/home/njh/svn/lib2geom/src/svg-path-parser.rl"
 	{
             _absolute = true;
         }
 	break;
 	case 5:
-#line 168 "/home/michael/2geom/src/svg-path-parser.rl"
+#line 168 "/home/njh/svn/lib2geom/src/svg-path-parser.rl"
 	{
             _absolute = false;
         }
 	break;
 	case 6:
-#line 172 "/home/michael/2geom/src/svg-path-parser.rl"
+#line 172 "/home/njh/svn/lib2geom/src/svg-path-parser.rl"
 	{
             _moveTo(_pop_point());
         }
 	break;
 	case 7:
-#line 176 "/home/michael/2geom/src/svg-path-parser.rl"
+#line 176 "/home/njh/svn/lib2geom/src/svg-path-parser.rl"
 	{
             _lineTo(_pop_point());
         }
 	break;
 	case 8:
-#line 180 "/home/michael/2geom/src/svg-path-parser.rl"
+#line 180 "/home/njh/svn/lib2geom/src/svg-path-parser.rl"
 	{
             _lineTo(Point(_pop_coord(X), _current[Y]));
         }
 	break;
 	case 9:
-#line 184 "/home/michael/2geom/src/svg-path-parser.rl"
+#line 184 "/home/njh/svn/lib2geom/src/svg-path-parser.rl"
 	{
             _lineTo(Point(_current[X], _pop_coord(Y)));
         }
 	break;
 	case 10:
-#line 188 "/home/michael/2geom/src/svg-path-parser.rl"
+#line 188 "/home/njh/svn/lib2geom/src/svg-path-parser.rl"
 	{
             Point p = _pop_point();
             Point c1 = _pop_point();
@@ -1517,7 +1517,7 @@ _match:
         }
 	break;
 	case 11:
-#line 195 "/home/michael/2geom/src/svg-path-parser.rl"
+#line 195 "/home/njh/svn/lib2geom/src/svg-path-parser.rl"
 	{
             Point p = _pop_point();
             Point c1 = _pop_point();
@@ -1525,7 +1525,7 @@ _match:
         }
 	break;
 	case 12:
-#line 201 "/home/michael/2geom/src/svg-path-parser.rl"
+#line 201 "/home/njh/svn/lib2geom/src/svg-path-parser.rl"
 	{
             Point p = _pop_point();
             Point c = _pop_point();
@@ -1533,14 +1533,14 @@ _match:
         }
 	break;
 	case 13:
-#line 207 "/home/michael/2geom/src/svg-path-parser.rl"
+#line 207 "/home/njh/svn/lib2geom/src/svg-path-parser.rl"
 	{
             Point p = _pop_point();
             _quadTo(_quad_tangent, p);
         }
 	break;
 	case 14:
-#line 212 "/home/michael/2geom/src/svg-path-parser.rl"
+#line 212 "/home/njh/svn/lib2geom/src/svg-path-parser.rl"
 	{
             Point point = _pop_point();
             bool sweep = _pop_flag();
@@ -1553,16 +1553,16 @@ _match:
         }
 	break;
 	case 15:
-#line 223 "/home/michael/2geom/src/svg-path-parser.rl"
+#line 223 "/home/njh/svn/lib2geom/src/svg-path-parser.rl"
 	{
             _closePath();
         }
 	break;
 	case 16:
-#line 360 "/home/michael/2geom/src/svg-path-parser.rl"
+#line 360 "/home/njh/svn/lib2geom/src/svg-path-parser.rl"
 	{goto _out;}
 	break;
-#line 1566 "/home/michael/2geom/src/svg-path-parser.cpp"
+#line 1566 "/home/njh/svn/lib2geom/src/svg-path-parser.cpp"
 		}
 	}
 
@@ -1571,7 +1571,7 @@ _again:
 	goto _resume;
 	_out: {}
 	}
-#line 370 "/home/michael/2geom/src/svg-path-parser.rl"
+#line 370 "/home/njh/svn/lib2geom/src/svg-path-parser.rl"
 
 
     if ( cs < svg_path_first_final ) {
diff --git a/src/2geom/svg-path.cpp b/src/2geom/svg-path.cpp
index 141ddbcf3..312db9d23 100644
--- a/src/2geom/svg-path.cpp
+++ b/src/2geom/svg-path.cpp
@@ -34,7 +34,7 @@
 namespace Geom {
 
 void output(Curve const &curve, SVGPathSink &sink) {
-    std::vector<Point> pts = sbasis_to_bezier(curve.sbasis(), 2); //TODO: use something better!
+    std::vector<Point> pts = sbasis_to_bezier(curve.toSBasis(), 2); //TODO: use something better!
     sink.curveTo(pts[0], pts[1], pts[2]);
 }
 
diff --git a/src/2geom/svg-path.h b/src/2geom/svg-path.h
index d09002220..ebc23a9b5 100644
--- a/src/2geom/svg-path.h
+++ b/src/2geom/svg-path.h
@@ -46,6 +46,7 @@ public:
                        bool large_arc, bool sweep, Point p) = 0;
     virtual void closePath() = 0;
     virtual void finish() = 0;
+    virtual ~SVGPathSink() {}
 };
 
 void output_svg_path(Path &path, SVGPathSink &sink);
@@ -66,14 +67,14 @@ public:
         _path.appendNew<LineSegment>(p);
     }
 
-    void curveTo(Point c0, Point c1, Point p) {
-        _path.appendNew<CubicBezier>(c0, c1, p);
-    }
-
     void quadTo(Point c, Point p) {
         _path.appendNew<QuadraticBezier>(c, p);
     }
 
+    void curveTo(Point c0, Point c1, Point p) {
+        _path.appendNew<CubicBezier>(c0, c1, p);
+    }
+
     void arcTo(double rx, double ry, double angle,
                bool large_arc, bool sweep, Point p)
     {
diff --git a/src/2geom/sweep.cpp b/src/2geom/sweep.cpp
new file mode 100644
index 000000000..4329bc559
--- /dev/null
+++ b/src/2geom/sweep.cpp
@@ -0,0 +1,104 @@
+#include "sweep.h"

+

+#include <algorithm>

+

+namespace Geom {

+

+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());

+    

+    for(unsigned i = 0; i < rs.size(); i++) {

+        events.push_back(Event(rs[i].left(), i, false));

+        events.push_back(Event(rs[i].right(), i, true));

+    }

+    std::sort(events.begin(), events.end());

+

+    std::vector<unsigned> open;

+    for(unsigned i = 0; i < events.size(); i++) {

+        unsigned ix = events[i].ix;

+        if(events[i].closing) {

+            std::vector<unsigned>::iterator iter = std::find(open.begin(), open.end(), ix);

+            //if(iter != open.end())

+            open.erase(iter);

+        } else {

+            for(unsigned j = 0; j < open.size(); j++) {

+                unsigned jx = open[j];

+                if(rs[jx][Y].intersects(rs[ix][Y])) {

+                    pairs[jx].push_back(ix);

+                }

+            }

+            open.push_back(ix);

+        }

+    }

+    return pairs;

+}

+

+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;

+    std::vector<Event> events[2];

+    events[0].reserve(a.size()*2);

+    events[1].reserve(b.size()*2);

+    

+    for(unsigned n = 0; n < 2; n++) {

+        unsigned sz = n ? b.size() : a.size();

+        events[n].reserve(sz*2);

+        for(unsigned i = 0; i < sz; i++) {

+            events[n].push_back(Event(n ? b[i].left() : a[i].left(), i, false));

+            events[n].push_back(Event(n ? b[i].right() : a[i].right(), i, true));

+        }

+        std::sort(events[n].begin(), events[n].end());

+    }

+

+    std::vector<unsigned> open[2];

+    bool n = events[1].front() < events[0].front();

+    for(unsigned i[] = {0,0}; i[n] < events[n].size();) {

+        unsigned ix = events[n][i[n]].ix;

+        bool closing = events[n][i[n]].closing;

+        //std::cout << n << "[" << ix << "] - " << (closing ? "closer" : "opener") << "\n";

+        if(closing) {

+            open[n].erase(std::find(open[n].begin(), open[n].end(), ix));

+        } else {

+            if(n) {

+                //n = 1

+                //opening a B, add to all open a

+                for(unsigned j = 0; j < open[0].size(); j++) {

+                    unsigned jx = open[0][j];

+                    if(a[jx][Y].intersects(b[ix][Y])) {

+                        pairs[jx].push_back(ix);

+                    }

+                }

+            } else {

+                //n = 0

+                //opening an A, add all open b

+                for(unsigned j = 0; j < open[1].size(); j++) {

+                    unsigned jx = open[1][j];

+                    if(b[jx][Y].intersects(a[ix][Y])) {

+                        pairs[ix].push_back(jx);

+                    }

+                }

+            }

+            open[n].push_back(ix);

+        }

+        i[n]++;

+        n = (events[!n][i[!n]] < events[n][i[n]]) ? !n : n;

+    }

+    return pairs;

+}

+

+//Fake cull, until the switch to the real sweep is made.

+std::vector<std::vector<unsigned> > fake_cull(unsigned a, unsigned b) {

+    std::vector<std::vector<unsigned> > ret;

+    

+    std::vector<unsigned> all;

+    for(unsigned j = 0; j < b; j++)

+        all.push_back(j);

+    

+    for(unsigned i = 0; i < a; i++)

+        ret.push_back(all);

+    

+    return ret;

+}

+

+}

diff --git a/src/2geom/sweep.h b/src/2geom/sweep.h
new file mode 100644
index 000000000..9587cec36
--- /dev/null
+++ b/src/2geom/sweep.h
@@ -0,0 +1,29 @@
+#ifndef __2GEOM_SWEEP_H__
+#define __2GEOM_SWEEP_H__
+
+#include <vector>
+#include "d2.h"
+
+namespace Geom {
+
+struct Event {
+    double x;
+    unsigned ix;
+    bool closing;
+    Event(double pos, unsigned i, bool c) : x(pos), ix(i), closing(c) {}
+// Lexicographic ordering by x then closing
+    bool operator<(Event const &other) const {
+        if(x < other.x) return true;
+        if(x > other.x) return false;
+        return closing < other.closing;
+    }
+
+};
+std::vector<std::vector<unsigned> > sweep_bounds(std::vector<Rect>);
+std::vector<std::vector<unsigned> > sweep_bounds(std::vector<Rect>, std::vector<Rect>);
+
+std::vector<std::vector<unsigned> > fake_cull(unsigned a, unsigned b);
+
+}
+
+#endif
diff --git a/src/2geom/transforms.h b/src/2geom/transforms.h
index 1cd1d3e5f..94058bff0 100644
--- a/src/2geom/transforms.h
+++ b/src/2geom/transforms.h
@@ -28,7 +28,7 @@ class Translate {
   public:
     explicit Translate(Point const &p) : vec(p) {}
     explicit Translate(Coord const x, Coord const y) : vec(x, y) {}
-    inline operator Matrix() const { return Matrix(0, 0, 0, 0, vec[X], vec[Y]); }
+    inline operator Matrix() const { return Matrix(1, 0, 0, 1, vec[X], vec[Y]); }
 
     inline Coord operator[](Dim2 const dim) const { return vec[dim]; }
     inline Coord operator[](unsigned const dim) const { return vec[dim]; }
diff --git a/src/2geom/utils.h b/src/2geom/utils.h
index a2f906ff4..ca880640d 100644
--- a/src/2geom/utils.h
+++ b/src/2geom/utils.h
@@ -38,6 +38,9 @@ public:
   NotImplemented() : std::logic_error("method not implemented") {}
 };
 
+// proper logical xor
+inline bool logical_xor (bool a, bool b) { return (a || b) && !(a && b); }
+
 /** Sign function - indicates the sign of a numeric type.  -1 indicates negative, 1 indicates
  *  positive, and 0 indicates, well, 0.  Mathsy people will know this is basically the derivative
  *  of abs, except for the fact that it is defined on 0.
diff --git a/src/live_effects/n-art-bpath-2geom.cpp b/src/live_effects/n-art-bpath-2geom.cpp
index f04c80f2b..9e5b966ea 100644
--- a/src/live_effects/n-art-bpath-2geom.cpp
+++ b/src/live_effects/n-art-bpath-2geom.cpp
@@ -45,7 +45,7 @@ static void curve_to_svgd(std::ostream & f, Geom::Curve const* c) {
 //    }

     else { 

         //this case handles sbasis as well as all other curve types

-        Geom::Path sbasis_path = path_from_sbasis(c->sbasis(), 0.1);

+        Geom::Path sbasis_path = Geom::path_from_sbasis(c->toSBasis(), 0.1);

 

         //recurse to convert the new path resulting from the sbasis to svgd

         for(Geom::Path::iterator iter = sbasis_path.begin(); iter != sbasis_path.end(); ++iter) {