merge in 2geom rev. 1154

(bzr r3756)

21 files changed, 998 insertions, 774 deletions

diff --git a/src/2geom/bezier-to-sbasis.h b/src/2geom/bezier-to-sbasis.h
index 03c99a9bd..f88913584 100644
--- a/src/2geom/bezier-to-sbasis.h
+++ b/src/2geom/bezier-to-sbasis.h
@@ -38,41 +38,25 @@
 
 namespace Geom{
 
-template <unsigned order>
-struct bezier_to_sbasis_impl {
-    static inline SBasis compute(Coord const *handles) {
-        return multiply(Linear(1, 0), bezier_to_sbasis_impl<order-1>::compute(handles)) +
-               multiply(Linear(0, 1), bezier_to_sbasis_impl<order-1>::compute(handles+1));
-    }
-};
-
-template <>
-struct bezier_to_sbasis_impl<1> {
-    static inline SBasis compute(Coord const *handles) {
+inline SBasis bezier_to_sbasis(Coord const *handles, unsigned order) {
+    if(order == 0)
+        return Linear(handles[0]);
+    else if(order == 1)
         return Linear(handles[0], handles[1]);
-    }
-};
-
-template <>
-struct bezier_to_sbasis_impl<0> {
-    static inline SBasis compute(Coord const *handles) {
-        return Linear(handles[0], handles[0]);
-    }
-};
-
-template <unsigned order>
-inline SBasis bezier_to_sbasis(Coord const *handles) {
-    return bezier_to_sbasis_impl<order>::compute(handles);
+    else
+        return multiply(Linear(1, 0), bezier_to_sbasis(handles, order-1)) +
+            multiply(Linear(0, 1), bezier_to_sbasis(handles+1, order-1));
 }
 
-template <unsigned order, typename T>
-inline D2<SBasis> handles_to_sbasis(T const &handles) {
+
+template <typename T>
+inline D2<SBasis> handles_to_sbasis(T const &handles, unsigned order) {
     double v[2][order+1];
     for(unsigned i = 0; i <= order; i++)
         for(unsigned j = 0; j < 2; j++)
              v[j][i] = handles[i][j];
-    return D2<SBasis>(bezier_to_sbasis<order>(v[0]),
-                      bezier_to_sbasis<order>(v[1]));
+    return D2<SBasis>(bezier_to_sbasis(v[0], order),
+                      bezier_to_sbasis(v[1], order));
 }
 
 };
diff --git a/src/2geom/bezier.h b/src/2geom/bezier.h
index 2c1b49213..e0d274e3b 100644
--- a/src/2geom/bezier.h
+++ b/src/2geom/bezier.h
@@ -1,4 +1,4 @@
-/*
+ /*
  * bezier.h
  *
  * Copyright 2007  MenTaLguY <mental@rydia.net>
@@ -34,6 +34,7 @@
 #define SEEN_BEZIER_H
 
 #include "coord.h"
+#include <valarray>
 #include "isnan.h"
 #include "bezier-to-sbasis.h"
 #include "d2.h"
@@ -42,8 +43,7 @@
 
 namespace Geom {
 
-template<unsigned order>
-Coord subdivideArr(Coord t, Coord const *v, Coord *left, Coord *right) {
+inline Coord subdivideArr(Coord t, Coord const *v, Coord *left, Coord *right, unsigned order) {
     Coord vtemp[order+1][order+1];
 
     /* Copy control points	*/
@@ -66,80 +66,93 @@ Coord subdivideArr(Coord t, Coord const *v, Coord *left, Coord *right) {
 }
 
 
-template <unsigned order>
 class Bezier {
 private:
-    Coord c_[order+1];
+    std::valarray<Coord> c_;
 
-    template<unsigned ord>
-    friend Bezier<ord> portion(const Bezier<ord> & a, Coord from, Coord to);
+    friend Bezier portion(const Bezier & a, Coord from, Coord to);
 
-    template<unsigned ord>
-    friend Interval bounds_fast(Bezier<ord> const & b);
+    friend Interval bounds_fast(Bezier const & b);
 
-    template<unsigned ord>
-    friend Bezier<ord-1> derivative(const Bezier<ord> & a);
+    friend Bezier derivative(const Bezier & a);
 
 protected:
-    Bezier(Coord const c[]) {
-        std::copy(c, c+order+1, c_);
+    Bezier(Coord const c[], unsigned ord) : c_(c, ord+1){
+        //std::copy(c, c+order()+1, &c_[0]);
     }
 
 public:
+    unsigned int order() const { return c_.size()-1;}
+    unsigned int size() const { return c_.size();}
+    
+    Bezier() :c_(0., 32) {}
+    Bezier(const Bezier& b) :c_(b.c_) {}
+    Bezier &operator=(Bezier const &other) {
+        if ( c_.size() != other.c_.size() ) {
+            c_.resize(other.c_.size());
+        }
+        c_ = other.c_;
+        return *this;
+    }
 
-    //TODO: fix this assert - get it compiling!
-    //template <unsigned required_order>
-    //static void assert_order(Bezier<required_order> const *) {}
+    struct Order {
+        unsigned order;
+        explicit Order(Bezier const &b) : order(b.order()) {}
+        explicit Order(unsigned o) : order(o) {}
+        operator unsigned() const { return order; }
+    };
 
-    Bezier() {}
+    //Construct an arbitrary order bezier
+    Bezier(Order ord) : c_(0., ord.order+1) {
+        assert(ord.order ==  order());
+    }
 
-    //Construct an order-0 bezier (constant Bézier)
-    explicit Bezier(Coord c0) {
-        //assert_order<0>(this);
+    explicit Bezier(Coord c0) : c_(0., 1) {
         c_[0] = c0;
     }
 
     //Construct an order-1 bezier (linear Bézier)
-    Bezier(Coord c0, Coord c1) {
-        //assert_order<1>(this);
+    Bezier(Coord c0, Coord c1) : c_(0., 2) {
         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);
+    Bezier(Coord c0, Coord c1, Coord c2) : c_(0., 3) {
         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);
+    Bezier(Coord c0, Coord c1, Coord c2, Coord c3) : c_(0., 4) {
         c_[0] = c0; c_[1] = c1; c_[2] = c2; c_[3] = c3;
     }
 
-    inline unsigned degree() const { return order; }
+    inline unsigned degree() const { return order(); }
 
     //IMPL: FragmentConcept
     typedef Coord output_type;
     inline bool isZero() const { 
-        for(unsigned i = 0; i <= order; i++) {
+        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++) {
+        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]; }
+    inline Coord at1() const { return c_[order()]; }
 
-    inline Coord valueAt(double t) const { return subdivideArr<order>(t, c_, NULL, NULL); }
+    inline Coord valueAt(double t) const { 
+        return subdivideArr(t, &c_[0], NULL, NULL, order()); 
+    }
     inline Coord operator()(double t) const { return valueAt(t); }
 
-    inline SBasis toSBasis() const { return bezier_to_sbasis<order>(c_); }
+    inline SBasis toSBasis() const { 
+        return bezier_to_sbasis(&c_[0], order());
+    }
 
     //Only mutator
     inline Coord &operator[](unsigned ix) { return c_[ix]; }
@@ -150,89 +163,91 @@ public:
      * 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]);
+        std::vector<Coord> val_n_der;
+        Coord d_[order()+1];
+        unsigned nn = n_derivs;
+        if(nn > order())
+            nn = order();
+        for(unsigned i = 0; i < size(); i++)
+            d_[i] = c_[i];
+        for(unsigned di = 0; di < nn; di++) {
+            val_n_der.push_back(subdivideArr(t, d_, NULL, NULL, order() - di));
+            for(unsigned i = 0; i < order() - di; i++) {
+                d_[i] = (order()-di)*(d_[i+1] - d_[i]);
             }
-            }*/
+        }
+        val_n_der.resize(n_derivs);
+        return val_n_der;
     }
   
-    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::pair<Bezier, Bezier > subdivide(Coord t) const {
+        Bezier a(Bezier::Order(*this)), b(Bezier::Order(*this));
+        subdivideArr(t, &c_[0], &a.c_[0], &b.c_[0], order());
+        return std::pair<Bezier, Bezier >(a, b);
     }
 
     std::vector<double> roots() const {
         std::vector<double> solutions;
-        find_bernstein_roots(c_, order, solutions, 0, 0.0, 1.0);
+        find_bernstein_roots(&c_[0], 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++)
+inline Bezier operator+(const Bezier & a, double v) {
+    Bezier result = Bezier(Bezier::Order(a));
+    for(unsigned i = 0; i <= a.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++)
+
+inline Bezier operator-(const Bezier & a, double v) {
+    Bezier result = Bezier(Bezier::Order(a));
+    for(unsigned i = 0; i <= a.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++)
+
+inline Bezier operator*(const Bezier & a, double v) {
+    Bezier result = Bezier(Bezier::Order(a));
+    for(unsigned i = 0; i <= a.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++)
+
+inline Bezier operator/(const Bezier & a, double v) {
+    Bezier result = Bezier(Bezier::Order(a));
+    for(unsigned i = 0; i <= a.order(); i++)
         result[i] = a[i] / v;
     return result;
 }
 
-template<unsigned order>
-Bezier<order> reverse(const Bezier<order> & a) {
-    Bezier<order> result;
-    for(unsigned i = 0; i <= order; i++)
-        result[i] = a[order - i];
+inline Bezier reverse(const Bezier & a) {
+    Bezier result = Bezier(Bezier::Order(a));
+    for(unsigned i = 0; i <= a.order(); i++)
+        result[i] = a[a.order() - i];
     return result;
 }
 
-template<unsigned order>
-Bezier<order> portion(const Bezier<order> & a, double from, double to) {
+inline Bezier portion(const Bezier & a, double from, double to) {
     //TODO: implement better?
-    Coord res[order+1];
+    Coord res[a.order()+1];
     if(from == 0) {
-        if(to == 1) { return Bezier<order>(a); }
-        subdivideArr<order>(to, a.c_, res, NULL);
-        return Bezier<order>(res);
+        if(to == 1) { return Bezier(a); }
+        subdivideArr(to, &a.c_[0], res, NULL, a.order());
+        return Bezier(res, a.order());
     }
-    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);
+    subdivideArr(from, &a.c_[0], NULL, res, a.order());
+    if(to == 1) return Bezier(res, a.order());
+    Coord res2[a.order()+1];
+    subdivideArr((to - from)/(1 - from), res, res2, NULL, a.order());
+    return Bezier(res2, a.order());
 }
 
-template<unsigned order>
-std::vector<Point> bezier_points(const D2<Bezier<order> > & a) {
+// XXX Todo: how to handle differing orders
+inline std::vector<Point> bezier_points(const D2<Bezier > & a) {
     std::vector<Point> result;
-    for(unsigned i = 0; i <= order; i++) {
+    for(unsigned i = 0; i <= a[0].order(); i++) {
         Point p;
         for(unsigned d = 0; d < 2; d++) p[d] = a[d][i];
         result.push_back(p);
@@ -240,33 +255,47 @@ std::vector<Point> bezier_points(const D2<Bezier<order> > & a) {
     return result;
 }
 
-template<unsigned order>
-Bezier<order-1> derivative(const Bezier<order> & a) {
-    Bezier<order-1> der;
+inline Bezier derivative(const Bezier & a) {
+    if(a.order() == 1) return Bezier(0.0);
+    Bezier der(Bezier::Order(a.order()-1));
     
-    for(unsigned i = 0; i < order; i++) {
-        der.c_[i] = order*(a.c_[i+1] - a.c_[i]);
+    for(unsigned i = 0; i < a.order(); i++) {
+        der.c_[i] = a.order()*(a.c_[i+1] - a.c_[i]);
     }
     return der;
 }
 
-template<unsigned order>
-inline Interval bounds_fast(Bezier<order> const & b) {
-    return Interval::fromArray(b.c_, order+1);
+inline Bezier integral(const Bezier & a) {
+    Bezier inte(Bezier::Order(a.order()+1));
+    
+    inte[0] = 0;
+    for(unsigned i = 0; i < inte.order(); i++) {
+        inte[i+1] = inte[i] + a[i]/(inte.order());
+    }
+    return inte;
+}
+
+inline Interval bounds_fast(Bezier const & b) {
+    return Interval::fromArray(&b.c_[0], b.size());
 }
 
 //TODO: better bounds exact
-template<unsigned order>
-inline Interval bounds_exact(Bezier<order> const & b) {
+inline Interval bounds_exact(Bezier const & b) {
     return bounds_exact(b.toSBasis());
 }
 
-template<unsigned order>
-inline Interval bounds_local(Bezier<order> const & b, Interval i) {
+inline Interval bounds_local(Bezier const & b, Interval i) {
     return bounds_fast(portion(b, i.min(), i.max()));
     //return bounds_local(b.toSBasis(), i);
 }
 
+inline std::ostream &operator<< (std::ostream &out_file, const Bezier & b) {
+    for(unsigned i = 0; i < b.size(); i++) {
+        out_file << b[i] << ", ";
+    }
+    return out_file;
+}
+
 }
 #endif //SEEN_BEZIER_H
 /*
diff --git a/src/2geom/crossing.h b/src/2geom/crossing.h
index 72b2eea4b..79912f024 100644
--- a/src/2geom/crossing.h
+++ b/src/2geom/crossing.h
@@ -16,7 +16,10 @@ struct Crossing {
     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; }
+    unsigned getOther(unsigned cur) const { return a == cur ? b : a; }
+    double getTime(unsigned cur) const { return a == cur ? ta : tb; }
+    double getOtherTime(unsigned cur) const { return a == cur ? tb : ta; }
+    bool onIx(unsigned ix) const { return a == ix || b == ix; }
 };
 
 
@@ -50,6 +53,7 @@ inline void sort_crossings(Crossings &cr, unsigned ix) { std::sort(cr.begin(), c
 
 template<typename T>
 struct Crosser {
+    virtual ~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());
diff --git a/src/2geom/d2.h b/src/2geom/d2.h
index 696ad0191..3a9e14bda 100644
--- a/src/2geom/d2.h
+++ b/src/2geom/d2.h
@@ -1,88 +1,88 @@
-/*

- * d2.h - Lifts one dimensional objects into 2d 

- *

- * 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.

- *

- */

-

-#ifndef _2GEOM_D2  //If this is change, change the guard in rect.h as well.

-#define _2GEOM_D2

-

-#include "point.h"

-#include "interval.h"

-#include "matrix.h"

-

-#include <boost/concept_check.hpp>

-#include "concepts.h"

-

-namespace Geom{

-

-template <class T>

-class D2{

-    //BOOST_CLASS_REQUIRE(T, boost, AssignableConcept);

-  private:

-    T f[2];

-

-  public:

-    D2() {f[X] = f[Y] = T();}

-    explicit D2(Point const &a) {

-        f[X] = T(a[X]); f[Y] = T(a[Y]);

-    }

-

-    D2(T const &a, T const &b) {

-        f[X] = a;

-        f[Y] = b;

-    }

-

-    //TODO: ask mental about operator= as seen in Point

-

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

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

-

-    //IMPL: FragmentConcept

-    typedef Point output_type;

-    bool isZero() const {

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

-        return f[X].isZero() && f[Y].isZero();

-    }

-    bool isFinite() const {

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

-        return f[X].isFinite() && f[Y].isFinite();

-    }

-    Point at0() const { 

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

-        return Point(f[X].at0(), f[Y].at0());

-    }

-    Point at1() const {

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

-        return Point(f[X].at1(), f[Y].at1());

-    }

-    Point valueAt(double t) const {

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

-        return (*this)(t);

+/*
+ * d2.h - Lifts one dimensional objects into 2d 
+ *
+ * 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.
+ *
+ */
+
+#ifndef _2GEOM_D2  //If this is change, change the guard in rect.h as well.
+#define _2GEOM_D2
+
+#include "point.h"
+#include "interval.h"
+#include "matrix.h"
+
+#include <boost/concept_check.hpp>
+#include "concepts.h"
+
+namespace Geom{
+
+template <class T>
+class D2{
+    //BOOST_CLASS_REQUIRE(T, boost, AssignableConcept);
+  private:
+    T f[2];
+
+  public:
+    D2() {f[X] = f[Y] = T();}
+    explicit D2(Point const &a) {
+        f[X] = T(a[X]); f[Y] = T(a[Y]);
+    }
+
+    D2(T const &a, T const &b) {
+        f[X] = a;
+        f[Y] = b;
+    }
+
+    //TODO: ask mental about operator= as seen in Point
+
+    T& operator[](unsigned i)              { return f[i]; }
+    T const & operator[](unsigned i) const { return f[i]; }
+
+    //IMPL: FragmentConcept
+    typedef Point output_type;
+    bool isZero() const {
+        boost::function_requires<FragmentConcept<T> >();
+        return f[X].isZero() && f[Y].isZero();
+    }
+    bool isFinite() const {
+        boost::function_requires<FragmentConcept<T> >();
+        return f[X].isFinite() && f[Y].isFinite();
+    }
+    Point at0() const { 
+        boost::function_requires<FragmentConcept<T> >();
+        return Point(f[X].at0(), f[Y].at0());
+    }
+    Point at1() const {
+        boost::function_requires<FragmentConcept<T> >();
+        return Point(f[X].at1(), f[Y].at1());
+    }
+    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),
@@ -92,314 +92,314 @@ class D2{
             res.push_back(Point(x[i], y[i]));
         }
         return res;
-    }

-    D2<SBasis> toSBasis() const {

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

-        return D2<SBasis>(f[X].toSBasis(), f[Y].toSBasis());

-    }

-

-    Point operator()(double t) const;

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

-};


-template <typename T>

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

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

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

-}

+    }
+    D2<SBasis> toSBasis() const {
+        boost::function_requires<FragmentConcept<T> >();
+        return D2<SBasis>(f[X].toSBasis(), f[Y].toSBasis());
+    }
+
+    Point operator()(double t) const;
+    Point operator()(double x, double y) const;
+};
+template <typename T>
+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>

-inline bool

-operator==(D2<T> const &a, D2<T> const &b) {

-    boost::function_requires<boost::EqualityComparableConcept<T> >();

-    return a[0]==b[0] && a[1]==b[1];

-}

-template <typename T>

-inline bool

-operator!=(D2<T> const &a, D2<T> const &b) {

-    boost::function_requires<boost::EqualityComparableConcept<T> >();

-    return a[0]!=b[0] || a[1]!=b[1];

-}

-

-//IMPL: NearConcept

-template <typename T>

-inline bool

-near(D2<T> const &a, D2<T> const &b, double tol) {

-    boost::function_requires<NearConcept<T> >();

-    return near(a[0], b[0]) && near(a[1], b[1]);

-}

-

-//IMPL: AddableConcept

-template <typename T>

-inline D2<T>

-operator+(D2<T> const &a, D2<T> const &b) {

-    boost::function_requires<AddableConcept<T> >();

-

-    D2<T> r;

-    for(unsigned i = 0; i < 2; i++)

-        r[i] = a[i] + b[i];

-    return r;

-}

-template <typename T>

-inline D2<T>

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

-    boost::function_requires<AddableConcept<T> >();

-

-    D2<T> r;

-    for(unsigned i = 0; i < 2; i++)

-        r[i] = a[i] - b[i];

-    return r;

-}

-template <typename T>

-inline D2<T>

-operator+=(D2<T> &a, D2<T> const &b) {

-    boost::function_requires<AddableConcept<T> >();

-

-    for(unsigned i = 0; i < 2; i++)

-        a[i] += b[i];

-    return a;

-}

-template <typename T>

-inline D2<T>

-operator-=(D2<T> &a, D2<T> const & b) {

-    boost::function_requires<AddableConcept<T> >();

-

-    for(unsigned i = 0; i < 2; i++)

-        a[i] -= b[i];

-    return a;

-}

-

-//IMPL: ScalableConcept

-template <typename T>

-inline D2<T>

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

-    boost::function_requires<ScalableConcept<T> >();

-    D2<T> r;

-    for(unsigned i = 0; i < 2; i++)

-        r[i] = -a[i];

-    return r;

-}

-template <typename T>

-inline D2<T>

-operator*(D2<T> const & a, Point const & b) {

-    boost::function_requires<ScalableConcept<T> >();

-

-    D2<T> r;

-    for(unsigned i = 0; i < 2; i++)

-        r[i] = a[i] * b[i];

-    return r;

-}

-template <typename T>

-inline D2<T>

-operator/(D2<T> const & a, Point const & b) {

-    boost::function_requires<ScalableConcept<T> >();

-    //TODO: b==0?

-    D2<T> r;

-    for(unsigned i = 0; i < 2; i++)

-        r[i] = a[i] / b[i];

-    return r;

-}

-template <typename T>

-inline D2<T>

-operator*=(D2<T> &a, Point const & b) {

-    boost::function_requires<ScalableConcept<T> >();

-

-    for(unsigned i = 0; i < 2; i++)

-        a[i] *= b[i];

-    return a;

-}

-template <typename T>

-inline D2<T>

-operator/=(D2<T> &a, Point const & b) {

-    boost::function_requires<ScalableConcept<T> >();

-    //TODO: b==0?

-    for(unsigned i = 0; i < 2; i++)

-        a[i] /= b[i];

-    return a;

-}

-

-template <typename T>

-inline D2<T> operator*(D2<T> const & a, double b) { return D2<T>(a[0]*b, a[1]*b); }

-template <typename T> 

-inline D2<T> operator*=(D2<T> & a, double b) { a[0] *= b; a[1] *= b; return a; }

-template <typename T>

-inline D2<T> operator/(D2<T> const & a, double b) { return D2<T>(a[0]/b, a[1]/b); }

-template <typename T> 

-inline D2<T> operator/=(D2<T> & a, double b) { a[0] /= b; a[1] /= b; return a; }

-

-template<typename T>

-D2<T> operator*(D2<T> const &v, Matrix const &m) {

-    boost::function_requires<AddableConcept<T> >();

-    boost::function_requires<ScalableConcept<T> >();

-    D2<T> ret;

-    for(unsigned i = 0; i < 2; i++)

-        ret[i] = v[X] * m[i] + v[Y] * m[i + 2] + m[i + 4];

-    return ret;

-}

-

-//IMPL: OffsetableConcept

-template <typename T>

-inline D2<T>

-operator+(D2<T> const & a, Point b) {

-    boost::function_requires<OffsetableConcept<T> >();

-    D2<T> r;

-    for(unsigned i = 0; i < 2; i++)

-        r[i] = a[i] + b[i];

-    return r;

-}

-template <typename T>

-inline D2<T>

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

-    boost::function_requires<OffsetableConcept<T> >();

-    D2<T> r;

-    for(unsigned i = 0; i < 2; i++)

-        r[i] = a[i] - b[i];

-    return r;

-}

-template <typename T>

-inline D2<T>

-operator+=(D2<T> & a, Point b) {

-    boost::function_requires<OffsetableConcept<T> >();

-    for(unsigned i = 0; i < 2; i++)

-        a[i] += b[i];

-    return a;

-}

-template <typename T>

-inline D2<T>

-operator-=(D2<T> & a, Point b) {

-    boost::function_requires<OffsetableConcept<T> >();

-    for(unsigned i = 0; i < 2; i++)

-        a[i] -= b[i];

-    return a;

-}

-

-template <typename T>

-inline T

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

-    boost::function_requires<AddableConcept<T> >();

-    boost::function_requires<MultiplicableConcept<T> >();

-

-    T r;

-    for(unsigned i = 0; i < 2; i++)

-        r += a[i] * b[i];

-    return r;

-}

-

-template <typename T>

-inline T

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

-    boost::function_requires<ScalableConcept<T> >();

-    boost::function_requires<MultiplicableConcept<T> >();

-

-    return a[1] * b[0] - a[0] * b[1];

-}

-

-

-//equivalent to cw/ccw, for use in situations where rotation direction doesn't matter.

-template <typename T>

-inline D2<T>

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

-    boost::function_requires<ScalableConcept<T> >();

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

-}

-

-//TODO: concepterize the following functions

-template <typename T>

-inline D2<T>

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

-    D2<T> r;

-    for(unsigned i = 0; i < 2; i++)

-        r[i] = compose(a[i],b);

-    return r;

-}

-

-template <typename T>

-inline D2<T>

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

-    D2<T> r;

-    for(unsigned i = 0; i < 2; i++)

-        r[i] = compose(a[i],b[i]);

-    return r;

-}

-

-template <typename T>

-inline D2<T>

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

-    D2<T> r;

-    for(unsigned i = 0; i < 2; i++)

-        r[i] = compose(a,b[i]);

-    return r;

-}

-

-

-template<typename T>

-inline Point

-D2<T>::operator()(double t) const {

-    Point p;

-    for(unsigned i = 0; i < 2; i++)

-       p[i] = (*this)[i](t);

-    return p;

-}

-

-//TODO: we might want to have this take a Point as the parameter.

-template<typename T>

-inline Point

-D2<T>::operator()(double x, double y) const {

-    Point p;

-    for(unsigned i = 0; i < 2; i++)

-       p[i] = (*this)[i](x, y);

-    return p;

-}

-

-

-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

-

-#include "rect.h"

-#include "d2-sbasis.h"

-

-namespace Geom{

-

-//Some D2 Fragment implementation which requires rect:

-template <typename T>

-Rect bounds_fast(const D2<T> &a) {

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

-    return Rect(bounds_fast(a[X]), bounds_fast(a[Y]));

-}

-template <typename T>

-Rect bounds_exact(const D2<T> &a) {

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

-    return Rect(bounds_exact(a[X]), bounds_exact(a[Y]));

-}

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

-}


-};

-

-/*

-  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

+
+//IMPL: boost::EqualityComparableConcept
+template <typename T>
+inline bool
+operator==(D2<T> const &a, D2<T> const &b) {
+    boost::function_requires<boost::EqualityComparableConcept<T> >();
+    return a[0]==b[0] && a[1]==b[1];
+}
+template <typename T>
+inline bool
+operator!=(D2<T> const &a, D2<T> const &b) {
+    boost::function_requires<boost::EqualityComparableConcept<T> >();
+    return a[0]!=b[0] || a[1]!=b[1];
+}
+
+//IMPL: NearConcept
+template <typename T>
+inline bool
+near(D2<T> const &a, D2<T> const &b, double tol) {
+    boost::function_requires<NearConcept<T> >();
+    return near(a[0], b[0]) && near(a[1], b[1]);
+}
+
+//IMPL: AddableConcept
+template <typename T>
+inline D2<T>
+operator+(D2<T> const &a, D2<T> const &b) {
+    boost::function_requires<AddableConcept<T> >();
+
+    D2<T> r;
+    for(unsigned i = 0; i < 2; i++)
+        r[i] = a[i] + b[i];
+    return r;
+}
+template <typename T>
+inline D2<T>
+operator-(D2<T> const &a, D2<T> const &b) {
+    boost::function_requires<AddableConcept<T> >();
+
+    D2<T> r;
+    for(unsigned i = 0; i < 2; i++)
+        r[i] = a[i] - b[i];
+    return r;
+}
+template <typename T>
+inline D2<T>
+operator+=(D2<T> &a, D2<T> const &b) {
+    boost::function_requires<AddableConcept<T> >();
+
+    for(unsigned i = 0; i < 2; i++)
+        a[i] += b[i];
+    return a;
+}
+template <typename T>
+inline D2<T>
+operator-=(D2<T> &a, D2<T> const & b) {
+    boost::function_requires<AddableConcept<T> >();
+
+    for(unsigned i = 0; i < 2; i++)
+        a[i] -= b[i];
+    return a;
+}
+
+//IMPL: ScalableConcept
+template <typename T>
+inline D2<T>
+operator-(D2<T> const & a) {
+    boost::function_requires<ScalableConcept<T> >();
+    D2<T> r;
+    for(unsigned i = 0; i < 2; i++)
+        r[i] = -a[i];
+    return r;
+}
+template <typename T>
+inline D2<T>
+operator*(D2<T> const & a, Point const & b) {
+    boost::function_requires<ScalableConcept<T> >();
+
+    D2<T> r;
+    for(unsigned i = 0; i < 2; i++)
+        r[i] = a[i] * b[i];
+    return r;
+}
+template <typename T>
+inline D2<T>
+operator/(D2<T> const & a, Point const & b) {
+    boost::function_requires<ScalableConcept<T> >();
+    //TODO: b==0?
+    D2<T> r;
+    for(unsigned i = 0; i < 2; i++)
+        r[i] = a[i] / b[i];
+    return r;
+}
+template <typename T>
+inline D2<T>
+operator*=(D2<T> &a, Point const & b) {
+    boost::function_requires<ScalableConcept<T> >();
+
+    for(unsigned i = 0; i < 2; i++)
+        a[i] *= b[i];
+    return a;
+}
+template <typename T>
+inline D2<T>
+operator/=(D2<T> &a, Point const & b) {
+    boost::function_requires<ScalableConcept<T> >();
+    //TODO: b==0?
+    for(unsigned i = 0; i < 2; i++)
+        a[i] /= b[i];
+    return a;
+}
+
+template <typename T>
+inline D2<T> operator*(D2<T> const & a, double b) { return D2<T>(a[0]*b, a[1]*b); }
+template <typename T> 
+inline D2<T> operator*=(D2<T> & a, double b) { a[0] *= b; a[1] *= b; return a; }
+template <typename T>
+inline D2<T> operator/(D2<T> const & a, double b) { return D2<T>(a[0]/b, a[1]/b); }
+template <typename T> 
+inline D2<T> operator/=(D2<T> & a, double b) { a[0] /= b; a[1] /= b; return a; }
+
+template<typename T>
+D2<T> operator*(D2<T> const &v, Matrix const &m) {
+    boost::function_requires<AddableConcept<T> >();
+    boost::function_requires<ScalableConcept<T> >();
+    D2<T> ret;
+    for(unsigned i = 0; i < 2; i++)
+        ret[i] = v[X] * m[i] + v[Y] * m[i + 2] + m[i + 4];
+    return ret;
+}
+
+//IMPL: OffsetableConcept
+template <typename T>
+inline D2<T>
+operator+(D2<T> const & a, Point b) {
+    boost::function_requires<OffsetableConcept<T> >();
+    D2<T> r;
+    for(unsigned i = 0; i < 2; i++)
+        r[i] = a[i] + b[i];
+    return r;
+}
+template <typename T>
+inline D2<T>
+operator-(D2<T> const & a, Point b) {
+    boost::function_requires<OffsetableConcept<T> >();
+    D2<T> r;
+    for(unsigned i = 0; i < 2; i++)
+        r[i] = a[i] - b[i];
+    return r;
+}
+template <typename T>
+inline D2<T>
+operator+=(D2<T> & a, Point b) {
+    boost::function_requires<OffsetableConcept<T> >();
+    for(unsigned i = 0; i < 2; i++)
+        a[i] += b[i];
+    return a;
+}
+template <typename T>
+inline D2<T>
+operator-=(D2<T> & a, Point b) {
+    boost::function_requires<OffsetableConcept<T> >();
+    for(unsigned i = 0; i < 2; i++)
+        a[i] -= b[i];
+    return a;
+}
+
+template <typename T>
+inline T
+dot(D2<T> const & a, D2<T> const & b) {
+    boost::function_requires<AddableConcept<T> >();
+    boost::function_requires<MultiplicableConcept<T> >();
+
+    T r;
+    for(unsigned i = 0; i < 2; i++)
+        r += a[i] * b[i];
+    return r;
+}
+
+template <typename T>
+inline T
+cross(D2<T> const & a, D2<T> const & b) {
+    boost::function_requires<ScalableConcept<T> >();
+    boost::function_requires<MultiplicableConcept<T> >();
+
+    return a[1] * b[0] - a[0] * b[1];
+}
+
+
+//equivalent to cw/ccw, for use in situations where rotation direction doesn't matter.
+template <typename T>
+inline D2<T>
+rot90(D2<T> const & a) {
+    boost::function_requires<ScalableConcept<T> >();
+    return D2<T>(-a[Y], a[X]);
+}
+
+//TODO: concepterize the following functions
+template <typename T>
+inline D2<T>
+compose(D2<T> const & a, T const & b) {
+    D2<T> r;
+    for(unsigned i = 0; i < 2; i++)
+        r[i] = compose(a[i],b);
+    return r;
+}
+
+template <typename T>
+inline D2<T>
+compose_each(D2<T> const & a, D2<T> const & b) {
+    D2<T> r;
+    for(unsigned i = 0; i < 2; i++)
+        r[i] = compose(a[i],b[i]);
+    return r;
+}
+
+template <typename T>
+inline D2<T>
+compose_each(T const & a, D2<T> const & b) {
+    D2<T> r;
+    for(unsigned i = 0; i < 2; i++)
+        r[i] = compose(a,b[i]);
+    return r;
+}
+
+
+template<typename T>
+inline Point
+D2<T>::operator()(double t) const {
+    Point p;
+    for(unsigned i = 0; i < 2; i++)
+       p[i] = (*this)[i](t);
+    return p;
+}
+
+//TODO: we might want to have this take a Point as the parameter.
+template<typename T>
+inline Point
+D2<T>::operator()(double x, double y) const {
+    Point p;
+    for(unsigned i = 0; i < 2; i++)
+       p[i] = (*this)[i](x, y);
+    return p;
+}
+
+
+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
+
+#include "rect.h"
+#include "d2-sbasis.h"
+
+namespace Geom{
+
+//Some D2 Fragment implementation which requires rect:
+template <typename T>
+Rect bounds_fast(const D2<T> &a) {
+    boost::function_requires<FragmentConcept<T> >();        
+    return Rect(bounds_fast(a[X]), bounds_fast(a[Y]));
+}
+template <typename T>
+Rect bounds_exact(const D2<T> &a) {
+    boost::function_requires<FragmentConcept<T> >();        
+    return Rect(bounds_exact(a[X]), bounds_exact(a[Y]));
+}
+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));
+}
+};
+
+/*
+  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/path-intersection.cpp b/src/2geom/path-intersection.cpp
index d641fcc08..2ad78e42f 100644
--- a/src/2geom/path-intersection.cpp
+++ b/src/2geom/path-intersection.cpp
@@ -105,6 +105,7 @@ int winding(Path const &path, Point p) {
  * hole.  Defaults to using the sign of area when it reaches funny cases.
  */
 bool path_direction(Path const &p) {
+    if(p.empty()) return false;
     //could probably be more efficient, but this is a quick job
     double y = p.initialPoint()[Y];
     double x = p.initialPoint()[X];
@@ -187,8 +188,8 @@ void pair_intersect(Curve const & A, double Al, double Ah,
     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)) {
+    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), 
@@ -425,6 +426,7 @@ CrossingSet crossings_among(std::vector<Path> const &p) {
 }
 */
 
+
 Crossings curve_self_crossings(Curve const &a) {
     Crossings res;
     std::vector<double> spl;
@@ -493,7 +495,6 @@ std::vector<std::vector<Rect> > curves_split_bounds(Path const &p, std::vector<s
     return ret;
 }
 
-
 Crossings path_self_crossings(Path const &p) {
     Crossings ret;
     std::vector<std::vector<unsigned> > cull = sweep_bounds(bounds(p));
@@ -535,7 +536,7 @@ Crossings path_self_crossings(Path const &p) {
 }
 */
 
-Crossings path_self_crossings(Path const &p) {
+Crossings 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++) {
@@ -551,7 +552,7 @@ Crossings path_self_crossings(Path const &p) {
                 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]);
+                        res2.push_back(res[k]);
                     }
                 }
                 res = res2;
@@ -563,6 +564,11 @@ Crossings path_self_crossings(Path const &p) {
     return ret;
 }
 
+void flip_crossings(Crossings &crs) {
+    for(unsigned i = 0; i < crs.size(); i++)
+        crs[i] = Crossing(crs[i].tb, crs[i].ta, crs[i].b, crs[i].a, !crs[i].dir);
+}
+
 CrossingSet crossings_among(std::vector<Path> const &p) {
     CrossingSet results(p.size(), Crossings());
     if(p.empty()) return results;
@@ -571,9 +577,11 @@ CrossingSet crossings_among(std::vector<Path> const &p) {
     
     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]);
+        Crossings res = 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);
+        flip_crossings(res);
+        merge_crossings(results[i], res, i);
         for(unsigned jx = 0; jx < cull[i].size(); jx++) {
             unsigned j = cull[i][jx];
             
@@ -583,6 +591,7 @@ CrossingSet crossings_among(std::vector<Path> const &p) {
             merge_crossings(results[j], res, j);
         }
     }
+    return results;
 }
 
 }
diff --git a/src/2geom/path-intersection.h b/src/2geom/path-intersection.h
index 3401386e0..6be04ad33 100644
--- a/src/2geom/path-intersection.h
+++ b/src/2geom/path-intersection.h
@@ -48,7 +48,7 @@ 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]; }
+Crossings self_crossings(Path const & a);
 
 inline Crossings crossings(Path const & a, Path const & b) {
     DefaultCrosser c = DefaultCrosser();
diff --git a/src/2geom/path.cpp b/src/2geom/path.cpp
index f05d3b0cf..79dc0a5f4 100644
--- a/src/2geom/path.cpp
+++ b/src/2geom/path.cpp
@@ -121,14 +121,14 @@ void Path::appendPortionTo(Path &ret, double from, double to) const {
     delete v;
     return;
   }
-  const_iterator toi   = inc(begin(), (unsigned)ti);
+  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) {
+  if(from >= to) {
     const_iterator ender = end();
     if(ender->initialPoint() == ender->finalPoint()) ender++;
     ret.insert(ret.end(), ++fromi, ender);
diff --git a/src/2geom/path.h b/src/2geom/path.h
index 6f39eb7ef..f4897fecc 100644
--- a/src/2geom/path.h
+++ b/src/2geom/path.h
@@ -128,16 +128,17 @@ public:
 template <unsigned order>
 class BezierCurve : public Curve {
 private:
-  D2<Bezier<order> > inner;
+  D2<Bezier > inner;
 public:
   template <unsigned required_degree>
   static void assert_degree(BezierCurve<required_degree> const *) {}
 
-  BezierCurve() {}
+  BezierCurve() : inner(Bezier::Order(order), Bezier::Order(order)) {
+  }
 
-  explicit BezierCurve(D2<Bezier<order> > const &x) : inner(x) {}
+  explicit BezierCurve(D2<Bezier > const &x) : inner(x) {}
   
-  BezierCurve(Bezier<order> x, Bezier<order> y) : inner(x, y) {}
+  BezierCurve(Bezier x, Bezier y) : inner(x, y) {}
 
   // default copy
   // default assign
@@ -145,19 +146,19 @@ public:
   BezierCurve(Point c0, Point c1) {
     assert_degree<1>(this);
     for(unsigned d = 0; d < 2; d++)
-        inner[d] = Bezier<order>(c0[d], c1[d]);
+        inner[d] = Bezier(c0[d], c1[d]);
   }
 
   BezierCurve(Point c0, Point c1, Point c2) {
     assert_degree<2>(this);
     for(unsigned d = 0; d < 2; d++)
-        inner[d] = Bezier<order>(c0[d], c1[d], c2[d]);
+        inner[d] = Bezier(c0[d], c1[d], c2[d]);
   }
 
   BezierCurve(Point c0, Point c1, Point c2, Point c3) {
     assert_degree<3>(this);
     for(unsigned d = 0; d < 2; d++)
-        inner[d] = Bezier<order>(c0[d], c1[d], c2[d], c3[d]);
+        inner[d] = Bezier(c0[d], c1[d], c2[d], c3[d]);
   }
 
   unsigned degree() const { return order; }
@@ -203,7 +204,7 @@ public:
   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);
+    std::pair<Bezier, Bezier > 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));
@@ -251,12 +252,12 @@ protected:
     for(unsigned i = 0; i <= order; i++) {
         x[i] = c[i][X]; y[i] = c[i][Y];
     }
-    inner = Bezier<order>(x, y);
+    inner = Bezier(x, y);
   }
 };
 
 // BezierCurve<0> is meaningless; specialize it out
-template<> class BezierCurve<0> : public BezierCurve<1> { public: BezierCurve(); BezierCurve(Bezier<0> x, Bezier<0> y); };
+template<> class BezierCurve<0> : public BezierCurve<1> { public: BezierCurve(); BezierCurve(Bezier x, Bezier y); };
 
 typedef BezierCurve<1> LineSegment;
 typedef BezierCurve<2> QuadraticBezier;
@@ -530,9 +531,10 @@ public:
 
   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());
+    for(unsigned i = 0; i <= size(); i++) {
+      std::vector<double> temp = (*this)[i].roots(v, d);
+      for(unsigned j = 0; j < temp.size(); j++)
+        res.push_back(temp[j] + i);
     }
     return res;
   }
diff --git a/src/2geom/rect.h b/src/2geom/rect.h
index f0811baf6..6a6d979a4 100644
--- a/src/2geom/rect.h
+++ b/src/2geom/rect.h
@@ -77,10 +77,10 @@ class D2<Interval> {
      *  (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());
+            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());
+            default: return Point(f[X].min(), f[Y].max());
         }
     }
         
diff --git a/src/2geom/region.cpp b/src/2geom/region.cpp
index cfab3a35a..116cc72fd 100644
--- a/src/2geom/region.cpp
+++ b/src/2geom/region.cpp
@@ -5,14 +5,6 @@
 
 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;
diff --git a/src/2geom/region.h b/src/2geom/region.h
index 4b434f1e5..e7eaa808b 100644
--- a/src/2geom/region.h
+++ b/src/2geom/region.h
@@ -34,12 +34,17 @@ class Region {
         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()); }
     
+    bool includes(Point const &p) const {
+        return logical_xor(!fill, contains(p));
+    }
+    
     Region inverse() const { return Region(boundary.reverse(), box, !fill); }
     
     Region operator*(Matrix const &m) const;
diff --git a/src/2geom/sbasis-geometric.cpp b/src/2geom/sbasis-geometric.cpp
index fde907857..b30c3a655 100644
--- a/src/2geom/sbasis-geometric.cpp
+++ b/src/2geom/sbasis-geometric.cpp
@@ -10,7 +10,7 @@
  *
  * The functions defined in this header related to 2d geometric operations such as arc length,
  * unit_vector, curvature, and centroid.  Most are built on top of unit_vector, which takes an
- * arbitrary D2 and returns an D2 with unit length with the same direction.
+ * arbitrary D2 and returns a D2 with unit length with the same direction.
  *
  * Todo/think about:
  *  arclength D2 -> sbasis (giving arclength function)
@@ -198,13 +198,13 @@ Geom::unitVector(D2<SBasis> const &V_in, double tol, unsigned order){
         r_eqn2 = Linear(1)-(a*a+b*b);
     }
     
-    //our candidat is:
+    //our candidate is:
     D2<SBasis> unitV;
     unitV[0] =  b;
     unitV[1] = -a;
 
     //is it good?
-    double rel_tol = max(1.,max(V_in[0].tailError(0),V_in[1].tailError(0)))*tol;
+    double rel_tol = std::max(1.,std::max(V_in[0].tailError(0),V_in[1].tailError(0)))*tol;
 
     if (r_eqn1.tailError(order)>rel_tol || r_eqn2.tailError(order)>tol){
         //if not: subdivide and concat results.
diff --git a/src/2geom/sbasis-math.cpp b/src/2geom/sbasis-math.cpp
index db4cf5e08..0140862b5 100644
--- a/src/2geom/sbasis-math.cpp
+++ b/src/2geom/sbasis-math.cpp
@@ -56,34 +56,34 @@ Piecewise<SBasis> abs(Piecewise<SBasis> const &f){
     return absf;
 }
 
-//-maxSb(x,y), minSb(x,y)--------------------------------------------------------
-Piecewise<SBasis> maxSb(          SBasis  const &f,           SBasis  const &g){
-    return maxSb(Piecewise<SBasis>(f),Piecewise<SBasis>(g));
+//-max(x,y), min(x,y)--------------------------------------------------------
+Piecewise<SBasis> max(          SBasis  const &f,           SBasis  const &g){
+    return max(Piecewise<SBasis>(f),Piecewise<SBasis>(g));
 }
-Piecewise<SBasis> maxSb(Piecewise<SBasis> const &f,           SBasis  const &g){
-    return maxSb(f,Piecewise<SBasis>(g));
+Piecewise<SBasis> max(Piecewise<SBasis> const &f,           SBasis  const &g){
+    return max(f,Piecewise<SBasis>(g));
 }
-Piecewise<SBasis> maxSb(          SBasis  const &f, Piecewise<SBasis> const &g){
-    return maxSb(Piecewise<SBasis>(f),g);
+Piecewise<SBasis> max(          SBasis  const &f, Piecewise<SBasis> const &g){
+    return max(Piecewise<SBasis>(f),g);
 }
-Piecewise<SBasis> maxSb(Piecewise<SBasis> const &f, Piecewise<SBasis> const &g){
-    Piecewise<SBasis> maxSb=partition(f,roots(f-g));
-    Piecewise<SBasis> gg =partition(g,maxSb.cuts);
-    maxSb = partition(maxSb,gg.cuts);
-    for (unsigned i=0; i<maxSb.size(); i++){
-        if (maxSb.segs[i](.5)<gg.segs[i](.5)) maxSb.segs[i]=gg.segs[i];
+Piecewise<SBasis> max(Piecewise<SBasis> const &f, Piecewise<SBasis> const &g){
+    Piecewise<SBasis> max=partition(f,roots(f-g));
+    Piecewise<SBasis> gg =partition(g,max.cuts);
+    max = partition(max,gg.cuts);
+    for (unsigned i=0; i<max.size(); i++){
+        if (max.segs[i](.5)<gg.segs[i](.5)) max.segs[i]=gg.segs[i];
     }
-    return maxSb;
+    return max;
 }
 
 Piecewise<SBasis> 
-minSb(          SBasis  const &f,           SBasis  const &g){ return -maxSb(-f,-g); }
+min(          SBasis  const &f,           SBasis  const &g){ return -max(-f,-g); }
 Piecewise<SBasis> 
-minSb(Piecewise<SBasis> const &f,           SBasis  const &g){ return -maxSb(-f,-g); }
+min(Piecewise<SBasis> const &f,           SBasis  const &g){ return -max(-f,-g); }
 Piecewise<SBasis> 
-minSb(          SBasis  const &f, Piecewise<SBasis> const &g){ return -maxSb(-f,-g); }
+min(          SBasis  const &f, Piecewise<SBasis> const &g){ return -max(-f,-g); }
 Piecewise<SBasis> 
-minSb(Piecewise<SBasis> const &f, Piecewise<SBasis> const &g){ return -maxSb(-f,-g); }
+min(Piecewise<SBasis> const &f, Piecewise<SBasis> const &g){ return -max(-f,-g); }
 
 
 //-sign(x)---------------------------------------------------------------
@@ -140,12 +140,12 @@ static Piecewise<SBasis> sqrt_internal(SBasis const &f,
 }
 
 Piecewise<SBasis> sqrt(SBasis const &f, double tol, int order){
-    return sqrt(maxSb(f,Linear(tol*tol)),tol,order);
+    return sqrt(max(f,Linear(tol*tol)),tol,order);
 }
 
 Piecewise<SBasis> sqrt(Piecewise<SBasis> const &f, double tol, int order){
     Piecewise<SBasis> result;
-    Piecewise<SBasis> ff=maxSb(f,Linear(tol*tol));
+    Piecewise<SBasis> ff=max(f,Linear(tol*tol));
 
     for (unsigned i=0; i<ff.size(); i++){
         Piecewise<SBasis> sqrtfi = sqrt_internal(ff.segs[i],tol,order);
diff --git a/src/2geom/sbasis-math.h b/src/2geom/sbasis-math.h
index 8f4e1612d..c20b09885 100644
--- a/src/2geom/sbasis-math.h
+++ b/src/2geom/sbasis-math.h
@@ -48,14 +48,14 @@ Piecewise<SBasis> abs(          SBasis const &f);
 Piecewise<SBasis> abs(Piecewise<SBasis>const &f);
 
 //- max(f,g), min(f,g) ----------------------------------------------
-Piecewise<SBasis> maxSb(          SBasis  const &f,           SBasis  const &g);
-Piecewise<SBasis> maxSb(Piecewise<SBasis> const &f,           SBasis  const &g);
-Piecewise<SBasis> maxSb(          SBasis  const &f, Piecewise<SBasis> const &g);
-Piecewise<SBasis> maxSb(Piecewise<SBasis> const &f, Piecewise<SBasis> const &g);
-Piecewise<SBasis> minSb(          SBasis  const &f,           SBasis  const &g);
-Piecewise<SBasis> minSb(Piecewise<SBasis> const &f,           SBasis  const &g);
-Piecewise<SBasis> minSb(          SBasis  const &f, Piecewise<SBasis> const &g);
-Piecewise<SBasis> minSb(Piecewise<SBasis> const &f, Piecewise<SBasis> const &g);
+Piecewise<SBasis> max(          SBasis  const &f,           SBasis  const &g);
+Piecewise<SBasis> max(Piecewise<SBasis> const &f,           SBasis  const &g);
+Piecewise<SBasis> max(          SBasis  const &f, Piecewise<SBasis> const &g);
+Piecewise<SBasis> max(Piecewise<SBasis> const &f, Piecewise<SBasis> const &g);
+Piecewise<SBasis> min(          SBasis  const &f,           SBasis  const &g);
+Piecewise<SBasis> min(Piecewise<SBasis> const &f,           SBasis  const &g);
+Piecewise<SBasis> min(          SBasis  const &f, Piecewise<SBasis> const &g);
+Piecewise<SBasis> min(Piecewise<SBasis> const &f, Piecewise<SBasis> const &g);
 
 //-sign(x)---------------------------------------------------------------
 Piecewise<SBasis> signSb(          SBasis const &f);
diff --git a/src/2geom/sbasis-roots.cpp b/src/2geom/sbasis-roots.cpp
index c4ef3c16d..52d3ef6a9 100644
--- a/src/2geom/sbasis-roots.cpp
+++ b/src/2geom/sbasis-roots.cpp
@@ -332,10 +332,8 @@ void subdiv_sbasis(SBasis const & s,
 
 std::vector<double> roots(SBasis const & s) {
     if(s.size() == 0) return std::vector<double>();
-    std::vector<double> b = sbasis_to_bezier(s), r;
     
-    find_bernstein_roots(&b[0], b.size()-1, r, 0, 0., 1.);
-    return r;
+    return sbasis_to_bezier(s).roots();
 }
 
 };
diff --git a/src/2geom/sbasis-to-bezier.cpp b/src/2geom/sbasis-to-bezier.cpp
index 206f18931..2484af18d 100644
--- a/src/2geom/sbasis-to-bezier.cpp
+++ b/src/2geom/sbasis-to-bezier.cpp
@@ -32,9 +32,8 @@ double W(unsigned n, unsigned j, unsigned k) {
 }
 
 // this produces a degree 2q bezier from a degree k sbasis
-std::vector<double>
+Bezier
 sbasis_to_bezier(SBasis const &B, unsigned q) {
-    std::vector<double> result;
     if(q == 0) {
         q = B.size();
         /*if(B.back()[0] == B.back()[1]) {
@@ -42,7 +41,7 @@ sbasis_to_bezier(SBasis const &B, unsigned q) {
             }*/
     }
     unsigned n = q*2;
-    result.resize(n, 0);
+    Bezier result = Bezier(Bezier::Order(n-1));
     if(q > B.size())
         q = B.size();
     n--;
diff --git a/src/2geom/sbasis-to-bezier.h b/src/2geom/sbasis-to-bezier.h
index d9eaabe7e..e566d7156 100644
--- a/src/2geom/sbasis-to-bezier.h
+++ b/src/2geom/sbasis-to-bezier.h
@@ -6,7 +6,7 @@
 
 namespace Geom{
 // this produces a degree k bezier from a degree k sbasis
-std::vector<double>
+Bezier
 sbasis_to_bezier(SBasis const &B, unsigned q = 0);
 
 std::vector<Geom::Point>
diff --git a/src/2geom/shape.cpp b/src/2geom/shape.cpp
index 670792521..54218d4d9 100644
--- a/src/2geom/shape.cpp
+++ b/src/2geom/shape.cpp
@@ -1,23 +1,39 @@
 #include "shape.h"
 #include "utils.h"
 #include "sweep.h"
+#include "ord.h"
 
 #include <iostream>
 #include <algorithm>
+#include <cstdlib>
 
 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());
 }
 
+//Orders a list of indices according to their containment within eachother.
+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]); }
+};
+
+//Returns the list of regions containing a particular point.  Useful in tandem with ContainmentOrder
+std::vector<unsigned> Shape::containment_list(Point p) const {
+    std::vector<Rect> pnt;
+    pnt.push_back(Rect(p, p));
+    std::vector<std::vector<unsigned> > cull = sweep_bounds(pnt, bounds(*this));
+    std::vector<unsigned> containers;
+    if(cull[0].size() == 0) return containers;
+    for(unsigned i = 0; i < cull[0].size(); i++)
+        if(content[cull[0][i]].contains(p)) containers.push_back(cull[0][i]);
+    return containers;
+}
+
 /* Used within shape_boolean and related functions, as the name describes, finds the
  * first false within the list of lists of booleans.
  */
@@ -77,8 +93,6 @@ Shape shape_boolean(bool rev, Shape const & a, Shape const & b, CrossingSet cons
     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) {
@@ -115,35 +129,32 @@ Shape shape_boolean(bool rev, Shape const & a, Shape const & b, CrossingSet cons
     
     //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);
+    //If true, outer paths are filled
+    bool const res_fill = rev ? (on_sub || (a.fill && b.fill)) : (a.fill && b.fill);
     
     //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);
+            bool env;
+            bool on_a = i < ac.size();
+            Region const & r(on_a ? ac[i] : bc[i - ac.size()]);
+            Shape const & other(on_a ? b : a);
             
-            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;
-                    }
-                }
+            std::vector<unsigned> containers = other.containment_list(r.boundary.initialPoint());
+            if(containers.empty()) {
+                //not included in any container, the environment fill is the opposite of the outer fill
+                env = !res_fill;
+                if(on_sub && logical_xor(other.fill, res_fill)) env = !env;  //If on the subtractor, invert the environment fill
+            } else {
+                //environment fill is the same as the inner-most container
+                std::vector<unsigned>::iterator cit = std::min_element(containers.begin(), containers.end(), ContainmentOrder(&other.content));
+                env = other[*cit].isFill();
             }
-            //disjoint
-            if(!mode) chunks.push_back(r);
-            skip: (void)0;
+            if(!logical_xor(rev, env)) chunks.push_back(r); //When unioning, environment must be hole for inclusion, when intersecting, it must be filled
         }
     }
     
-    return Shape(chunks);
+    return Shape(chunks, res_fill);
 }
 
 // Just a convenience wrapper for shape_boolean, which handles the crossings
@@ -176,8 +187,11 @@ Shape shape_boolean_rb(bool rev, Shape const &a, Shape const &b, CrossingSet con
  * 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_ .
+ *
+ * NOTE: currently doesn't work, as the CrossingSet reversal functions crash
  */
 Shape boolop(Shape const &a, Shape const &b, unsigned flags, CrossingSet const &crs) {
+    throw NotImplemented();
     flags &= 15;
     if(flags <= BOOLOP_UNION) {
         switch(flags) {
@@ -194,6 +208,7 @@ Shape boolop(Shape const &a, Shape const &b, unsigned flags, CrossingSet const &
             case BOOLOP_UNION:        return shape_boolean(false, a, b);
         }
     } else {
+        flags = ~flags & 15;
         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);
@@ -203,7 +218,7 @@ Shape boolop(Shape const &a, Shape const &b, unsigned flags, CrossingSet const &
                 return res;
             }
         }
-        return boolop(a, b, ~flags, crs).inverse();
+        return boolop(a, b, flags, crs).inverse();
     }
     return Shape();
 }
@@ -230,7 +245,8 @@ Shape boolop(Shape const &a, Shape const &b, unsigned flags) {
             case BOOLOP_UNION:        return shape_boolean(false, a, b);
         }
     } else {
-        switch(flags - BOOLOP_NEITHER) {
+        flags = ~flags & 15;
+        switch(flags) {
             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: {
@@ -239,150 +255,323 @@ Shape boolop(Shape const &a, Shape const &b, unsigned flags) {
                 return res;
             } //return boolop(a, b, flags, crossings_between(a, b));
         }
-        return boolop(a, b, ~flags).inverse();
+        return boolop(a, b, flags).inverse();
     }
     return Shape();
 }
 
-
 int paths_winding(std::vector<Path> const &ps, Point p) {
-    int ret;
+    int ret = 0;
     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;
+void add_to_shape(Shape &s, Path const &p, bool fill) {
+    if(fill)
+        s.content.push_back(Region(p).asFill());
+    else
+        s.content.push_back(Region(p).asHole());
+}
+
+int inner_winding(Path const & p, std::vector<Path> const &ps) {
+    Point pnt = p.initialPoint();
+    return paths_winding(ps, pnt) - winding(p, pnt) + 1;
+}
+
+double fudgerize(double d, bool rev) {
+    double ret = rev ? d - 0.01 : d + 0.01;
+    if(ret < 0) ret = 0;
+    return ret;
+}
+
+unsigned pick_coincident(unsigned ix, unsigned jx, bool &rev, std::vector<Path> const &ps, CrossingSet const &crs) {
+    unsigned ex_jx = jx;
+    unsigned oix = crs[ix][jx].getOther(ix);
+    double otime = crs[ix][jx].getTime(oix);
+    Point cross_point = ps[oix].pointAt(otime),
+          along = ps[oix].pointAt(fudgerize(otime, rev)) - cross_point,
+          prev = -along;
+    bool ex_dir = rev;
+    for(unsigned k = jx; k < crs[ix].size(); k++) {
+        unsigned koix = crs[ix][k].getOther(ix);
+        if(koix == oix) {
+            if(!near(otime, crs[ix][k].getTime(oix))) break;
+            for(unsigned dir = 0; dir < 2; dir++) {
+                Point val = ps[ix].pointAt(fudgerize(crs[ix][k].getTime(ix), dir)) - cross_point;
+                Cmp to_prev = cmp(cross(val, prev), 0);
+                Cmp from_along = cmp(cross(along, val), 0);
+                Cmp c = cmp(from_along, to_prev);
+                if(c == EQUAL_TO && from_along == LESS_THAN) {
+                    ex_jx = k;
+                    prev = val;
+                    ex_dir = dir;
+                }
+            }
+        }
     }
-    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);
+    rev = ex_dir;
+    return ex_jx;
+}
+
+unsigned crossing_along(double t, unsigned ix, unsigned jx, bool dir, Crossings const & crs) {
+    Crossing cur = Crossing(t, t, ix, ix, false);
+    if(jx < crs.size()) {
+        double ct = crs[jx].getTime(ix);
+        if(t == ct) {
+            cur = crs[jx];
+            if(cur.a == cur.b) {
+                if(jx+1 <= crs.size() && crs[jx+1].getOther(ix) == ix) return jx+1;
+                if(jx > 0 && crs[jx-1].getOther(ix) == ix) return jx-1;
+            }
         }
     }
-};
+    if(!dir) {
+        jx = std::upper_bound(crs.begin(), crs.end(), cur, CrossingOrder(ix)) - crs.begin();
+    } else {
+        jx = std::lower_bound(crs.begin(), crs.end(), cur, CrossingOrder(ix)) - crs.begin();
+        if(jx == 0) jx = crs.size() - 1; else jx--;
+        jx = std::lower_bound(crs.begin(), crs.end(), crs[jx], CrossingOrder(ix)) - crs.begin();
+    }
+    if(jx >= crs.size()) jx = 0;
+    return jx;
+}
 
-typedef std::vector<Edge> Edges;
+void crossing_dual(unsigned &i, unsigned &j, CrossingSet const & crs) {
+    Crossing cur = crs[i][j];
+    i = cur.getOther(i);
+    std::cout << i << "\n";
+    if(crs[i].empty())
+        j = 0;
+    else
+        j = std::lower_bound(crs[i].begin(), crs[i].end(), cur, CrossingOrder(i)) - crs[i].begin();
+}
 
-double point_cosine(Point a, Point b, Point c) {
-    Point db = b - a, dc = c - a;
-    return cross(db, dc) / (db.length() * dc.length());
+//locate a crossing on the outside, by casting a ray through the middle of the bbox
+void outer_crossing(unsigned &ix, unsigned &jx, bool & dir, std::vector<Path> const & ps, CrossingSet const & crs) {
+    Rect bounds = ps[ix].boundsFast();
+    double ry = bounds[Y].middle();
+    double max_val = bounds.left(), max_t = 0;
+    ix = ps.size();
+    for(unsigned i = 0; i < ps.size(); i++) {
+        if(!crs[i].empty()) {
+            std::vector<double> rts = ps[i].roots(ry, Y);
+            for(unsigned j = 0; j < rts.size(); j++) {
+                double val = ps[i].valueAt(rts[j], X);
+                if(val > max_val) {
+                    ix = i;
+                    max_val = val;
+                    max_t = rts[j];
+                }
+            }
+        }
+    }
+    if(ix != ps.size()) {
+        dir = ps[ix].valueAt(max_t + 0.01, Y) >
+              ps[ix].valueAt(max_t - 0.01, Y);
+        jx = crossing_along(max_t, ix, jx, dir, crs[ix]);
+    }
 }
 
-//sanitize
-Regions regionize_paths(std::vector<Path> const &ps, bool evenodd) {
-    CrossingSet crs = crossings_among(ps);
+std::vector<Path> inner_sanitize(std::vector<Path> const & ps) {
+    CrossingSet crs(crossings_among(ps));
     
-    Edges es;
+    Regions chunks;
     
-    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;
+    std::vector<bool> used_path(ps.size(), false);
+    std::vector<std::vector<bool> > visited;
+    for(unsigned i = 0; i < crs.size(); i++)
+        visited.push_back(std::vector<bool>(crs[i].size(), false));
+    
+    std::vector<Path> result_paths;
+    
+    while(true) {
+        unsigned ix = 0, jx = 0;
+        bool dir = false;
+
+        //find an outer crossing by trying various paths and checking if the crossings are used
+        for(; ix < crs.size(); ix++) {
+            //TODO: optimize so it doesn't unecessarily check stuff
+            bool cont = true;
+            for(unsigned j = 0; j < crs[ix].size(); j++) {
+                if(!visited[ix][j]) { cont = false; break; }
+            }
+            if(cont) continue;
+            unsigned rix = ix, rjx = jx;
+            outer_crossing(rix, rjx, dir, ps, crs);
+            if(rix >= crs.size() || visited[rix][rjx]) continue;
+            ix = rix; jx = rjx;
+            break;
+        }
+        if(ix == crs.size()) break;
+        crossing_dual(ix, jx, crs);
+
+        dir = !dir;
+
+        Path res;
+        do {
+            visited[ix][jx] = true;
+            //unsigned nix = ix, njx = jx;
+            //crossing_dual(nix, njx, crs);
+            //visited[nix][njx] = true;
+            unsigned fix = ix, fjx = jx;
             
-            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));
+            bool new_dir = dir;
             
-            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));
-        }
+            jx = crossing_along(crs[ix][jx].getTime(ix), ix, jx, dir, crs[ix]);
+            if(crs[ix][jx].a != crs[ix][jx].b) crossing_dual(ix, jx, crs); else new_dir = !new_dir;
+            jx = pick_coincident(ix, jx, new_dir, ps, crs);
+            
+            //unsigned nix = ix, njx = jx;
+            //crossing_dual(nix, njx, crs);
+            
+            Crossing from = crs[fix][fjx],
+                     to = crs[ix][jx];
+            if(dir) {
+                // backwards
+                std::cout << "r" << ix << "[" << from.getTime(ix)  << ", " << to.getTime(ix) << "]\n";
+                Path p = ps[ix].portion(from.getTime(ix), to.getTime(ix)).reverse();
+                for(unsigned i = 0; i < p.size(); i++)
+                    res.append(p[i]);
+            } else {
+                // forwards
+                std::cout << "f" << ix << "[" << from.getTime(ix) << ", " << to.getTime(ix) << "]\n";
+                ps[ix].appendPortionTo(res, from.getTime(ix), to.getTime(ix));
+            }
+            dir = new_dir;
+        } while(!visited[ix][jx]);
+        std::cout << "added " << res.size() << "\n";
+        result_paths.push_back(res);
     }
-    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));
-        }
+    for(unsigned i = 0; i < crs.size(); i++) {
+        if(crs[i].empty() && !used_path[i])
+            result_paths.push_back(ps[i]);
     }
-    
-    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]);
+    return result_paths;
+}
+
+Shape sanitize(std::vector<Path> const & ps) {
+    std::vector<Path> res;
+    for(unsigned i = 0; i < ps.size(); i++) {
+        append(res, inner_sanitize(std::vector<Path>(1, ps[i])));
+    }
+    return stopgap_cleaner(res);
+}
+
+/*  WIP sanitizer:
+unsigned pick_coincident(unsigned ix, unsigned jx, bool pref, bool &rev, std::vector<Path> const &ps, CrossingSet const &crs) {
+    unsigned ex_jx = jx;
+    unsigned oix = crs[ix][jx].getOther(ix);
+    double otime = crs[ix][jx].getTime(oix);
+    Point cross_point = ps[oix].pointAt(otime),
+          along = ps[oix].pointAt(otime + (rev ? -0.01 : 0.01)) - cross_point,
+          prev = -along;
+    bool ex_dir = rev;
+    for(unsigned k = jx; k < crs[ix].size(); k++) {
+        unsigned koix = crs[ix][k].getOther(ix);
+        if(koix == oix) {
+            if(!near(otime, crs[ix][k].getTime(oix))) break;
+            for(unsigned dir = 0; dir < 2; dir++) {
+                Point val = ps[ix].pointAt(crs[ix][k].getTime(ix) + (dir ? -0.01 : 0.01)) - cross_point;
+                Cmp to_prev = cmp(cross(val, prev), 0);
+                Cmp from_along = cmp(cross(along, val), 0);
+                Cmp c = cmp(from_along, to_prev);
+                if(c == EQUAL_TO && (from_along == LESS_THAN) == pref) {
+                    ex_jx = k;
+                    prev = val;
+                    ex_dir = dir;
+                }
+            }
+        }
     }
-    es = es2;
+    rev = ex_dir;
+    return ex_jx;
+}
+
+unsigned corner_index(unsigned &i) {
+    div_t div_res = div(i, 4);
+    i = div_res.quot;
+    return div_res.rem;
+}
+
+bool corner_direction(unsigned ix, unsigned jc, unsigned corner, CrossingSet const &crs) {
+    if(crs[ix][jc].a == ix) return corner > 1; else return corner %2 == 1;
+}
+
+Shape sanitize(std::vector<Path> const & ps) {
+    CrossingSet crs = crossings_among(ps);
     
-    std::cout << es.size() << " edges\n";
+    //Keep track of which CORNERS we've hit.
+    // FF FR RF RR, first is A dir, second B dir
+    std::vector<std::vector<bool> > visited;
+    for(unsigned i = 0; i < crs.size(); i++)
+        visited.push_back(std::vector<bool>(crs[i].size()*4, false));
     
     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;
+        unsigned i, j;
+        first_false(visited, i, j);
+        unsigned corner = corner_index(j);
+        
+        if(i == visited.size()) break;
+        
+        bool dir = corner_direction(i, j, corner, crs);
+        
+        //Figure out whether we hug the path cw or ccw, based on the orientation of the initial corner:        
+        unsigned oix = crs[i][j].getOther(i);
+        double otime = crs[i][j].getTime(oix);
+        bool odir = (oix == crs[i][j].a) ? corner > 1 : corner % 2 == 1;
+        Point cross_point = ps[oix].pointAt(otime),
+              along = ps[oix].pointAt(otime + (odir ? -0.01 : 0.01)) - cross_point,
+                val = ps[i].pointAt(crs[i][j].getTime(i) + (dir ? -0.01 : 0.01)) - cross_point;
+        
+        Cmp from_along = cmp(cross(along, val), 0);
+        bool cw = from_along == LESS_THAN;
+        std::cout << "cw = " << cw << "\n";
         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;
-                    }
-                }
+            Crossing cur = crs[i][j];
+            visited[i][j*4+corner] = true;
+            
+            unsigned fix = i, fjx = j;
+            crossing_dual(i, j, crs);
+            visited[i][j*4+corner] = true;
+            i = fix; j = fjx;
+            
+            j = crossing_along(crs[i][j].getTime(i), i, j, dir, crs[i]);
+            
+            crossing_dual(i, j, crs);
+            
+            bool new_dir = dir;
+            pick_coincident(i, j, cw, new_dir, ps, crs);
+            
+            Crossing from = crs[fix][fjx],
+                     to = crs[i][j];
+            if(dir) {
+                // backwards
+                std::cout << "r" << i << "[" << to.getTime(i)  << ", " << from.getTime(i) << "]\n";
+                Path p = ps[i].portion(to.getTime(i) + 0.001, from.getTime(i)).reverse();
+                for(unsigned k = 0; k < p.size(); k++)
+                    res.append(p[k]);
+            } else {
+                // forwards
+                std::cout << "f" << i << "[" << from.getTime(i) << ", " << to.getTime(i) << "]\n";
+                ps[i].appendPortionTo(res, from.getTime(i) + 0.001, to.getTime(i));
             }
-            i = poss[best];
-        } while(!used[i]);
+            if(i == to.a)
+                corner = (new_dir ? 2 : 0) + (dir ? 1 : 0);
+            else
+                corner = (new_dir ? 1 : 0) + (dir ? 2 : 0);
+            dir = new_dir;
+        } while(!visited[i][j*4+corner]);
         chunks.push_back(Region(res));
+//        if(use) {
+//            chunks.push_back(Region(res, true));
+//        }
     }
-    return chunks;
-}
+    return Shape(chunks);
+//    return ret;
+} */
 
 /* 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.
@@ -404,18 +593,19 @@ Shape Shape::inverse() const {
     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();
+    std::vector<unsigned> containers = containment_list(p);
+    if(containers.empty()) return !isFill();
+    unsigned ix = *min_element(containers.begin(), containers.end(), ContainmentOrder(&content));
+    return content[ix].isFill();
+}
+
+Shape stopgap_cleaner(std::vector<Path> const &ps) {
+    if(ps.empty()) return Shape(false);
+    Shape ret;
+    for(unsigned i = 0; i < ps.size(); i++)
+        add_to_shape(ret, ps[i], inner_winding(ps[i], ps) % 2 != 0);
+    return ret;
 }
 
 bool Shape::inside_invariants() const {  //semi-slow & easy to violate
diff --git a/src/2geom/shape.h b/src/2geom/shape.h
index 3700e9e5a..b9194537c 100644
--- a/src/2geom/shape.h
+++ b/src/2geom/shape.h
@@ -36,7 +36,7 @@ class Shape {
     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 &);
-
+    friend void add_to_shape(Shape &s, Path const &p, bool);
   public:
     Shape() : fill(true) {}
     explicit Shape(Region const & r) {
@@ -61,9 +61,10 @@ class Shape {
     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)
+    bool invariants() const;      //vera slow (combo, checks the above)
 
-  private:     
+  private:
+    std::vector<unsigned> containment_list(Point p) const;
     void update_fill() const {
         unsigned ix = outer_index(content);
         if(ix < size())
@@ -80,10 +81,21 @@ inline CrossingSet crossings_between(Shape const &a, Shape const &b) { return cr
 Shape shape_boolean(bool rev, Shape const &, Shape const &, CrossingSet const &);
 Shape shape_boolean(bool rev, Shape const &, Shape const &);
 
+//unsigned pick_coincident(unsigned ix, unsigned jx, bool &rev, std::vector<Path> const &ps, CrossingSet const &crs);
+//void outer_crossing(unsigned &ix, unsigned &jx, bool & dir, std::vector<Path> const & ps, CrossingSet const & crs);
+void crossing_dual(unsigned &i, unsigned &j, CrossingSet const & crs);
+unsigned crossing_along(double t, unsigned ix, unsigned jx, bool dir, Crossings const & crs);
+
 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);
+Shape sanitize(std::vector<Path> const &ps);
+
+Shape stopgap_cleaner(std::vector<Path> const &ps);
+
+inline std::vector<Path> desanitize(Shape const & s) {
+    return paths_from_regions(s.getContent());
+}
 
 }
 
diff --git a/src/2geom/svg-path-parser.cpp b/src/2geom/svg-path-parser.cpp
index 0bd15e8b9..ba0fbe7da 100644
--- a/src/2geom/svg-path-parser.cpp
+++ b/src/2geom/svg-path-parser.cpp
@@ -1,4 +1,4 @@
-#line 1 "/home/njh/svn/lib2geom/src/svg-path-parser.rl"
+#line 1 "/home/mental/trees/lib2geom/src/svg-path-parser.rl"
 /*
  * parse SVG path specifications
  *
@@ -129,7 +129,7 @@ private:
 };
 
 
-#line 133 "/home/njh/svn/lib2geom/src/svg-path-parser.cpp"
+#line 133 "/home/mental/trees/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, 769, 
-	769, 770, 770, 771, 769, 772, 0, 781, 
+	575, 576, 576, 577, 575, 578, 0, 761, 
+	761, 762, 762, 763, 761, 764, 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, 765, 765, 766, 
-	766, 767, 765, 768, 0, 777, 777, 778, 
+	310, 295, 296, 321, 0, 757, 757, 758, 
+	758, 759, 757, 760, 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, 761, 761, 762, 762, 763, 
-	761, 764, 0, 757, 757, 758, 758, 759, 
-	757, 760, 0, 599, 599, 600, 600, 601, 
+	124, 149, 0, 769, 769, 770, 770, 771, 
+	769, 772, 0, 765, 765, 766, 766, 767, 
+	765, 768, 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, 11, 62, 62, 
-	62, 11, 62, 62, 62, 9, 59, 59, 
-	59, 9, 59, 59, 59, 0, 1, 1, 
+	17, 3, 17, 0, 0, 9, 59, 59, 
+	59, 9, 59, 59, 59, 11, 62, 62, 
+	62, 11, 62, 62, 62, 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/njh/svn/lib2geom/src/svg-path-parser.rl"
+#line 133 "/home/mental/trees/lib2geom/src/svg-path-parser.rl"
 
 
 void Parser::parse(char const *str)
@@ -1369,7 +1369,7 @@ throw(SVGPathParseError)
     _reset();
 
     
-#line 1373 "/home/njh/svn/lib2geom/src/svg-path-parser.cpp"
+#line 1373 "/home/mental/trees/lib2geom/src/svg-path-parser.cpp"
 	{
 	cs = svg_path_start;
 	}
@@ -1445,13 +1445,13 @@ _match:
 		switch ( *_acts++ )
 		{
 	case 0:
-#line 145 "/home/njh/svn/lib2geom/src/svg-path-parser.rl"
+#line 145 "/home/mental/trees/lib2geom/src/svg-path-parser.rl"
 	{
             start = p;
         }
 	break;
 	case 1:
-#line 149 "/home/njh/svn/lib2geom/src/svg-path-parser.rl"
+#line 149 "/home/mental/trees/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/njh/svn/lib2geom/src/svg-path-parser.rl"
+#line 156 "/home/mental/trees/lib2geom/src/svg-path-parser.rl"
 	{
             _push(1.0);
         }
 	break;
 	case 3:
-#line 160 "/home/njh/svn/lib2geom/src/svg-path-parser.rl"
+#line 160 "/home/mental/trees/lib2geom/src/svg-path-parser.rl"
 	{
             _push(0.0);
         }
 	break;
 	case 4:
-#line 164 "/home/njh/svn/lib2geom/src/svg-path-parser.rl"
+#line 164 "/home/mental/trees/lib2geom/src/svg-path-parser.rl"
 	{
             _absolute = true;
         }
 	break;
 	case 5:
-#line 168 "/home/njh/svn/lib2geom/src/svg-path-parser.rl"
+#line 168 "/home/mental/trees/lib2geom/src/svg-path-parser.rl"
 	{
             _absolute = false;
         }
 	break;
 	case 6:
-#line 172 "/home/njh/svn/lib2geom/src/svg-path-parser.rl"
+#line 172 "/home/mental/trees/lib2geom/src/svg-path-parser.rl"
 	{
             _moveTo(_pop_point());
         }
 	break;
 	case 7:
-#line 176 "/home/njh/svn/lib2geom/src/svg-path-parser.rl"
+#line 176 "/home/mental/trees/lib2geom/src/svg-path-parser.rl"
 	{
             _lineTo(_pop_point());
         }
 	break;
 	case 8:
-#line 180 "/home/njh/svn/lib2geom/src/svg-path-parser.rl"
+#line 180 "/home/mental/trees/lib2geom/src/svg-path-parser.rl"
 	{
             _lineTo(Point(_pop_coord(X), _current[Y]));
         }
 	break;
 	case 9:
-#line 184 "/home/njh/svn/lib2geom/src/svg-path-parser.rl"
+#line 184 "/home/mental/trees/lib2geom/src/svg-path-parser.rl"
 	{
             _lineTo(Point(_current[X], _pop_coord(Y)));
         }
 	break;
 	case 10:
-#line 188 "/home/njh/svn/lib2geom/src/svg-path-parser.rl"
+#line 188 "/home/mental/trees/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/njh/svn/lib2geom/src/svg-path-parser.rl"
+#line 195 "/home/mental/trees/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/njh/svn/lib2geom/src/svg-path-parser.rl"
+#line 201 "/home/mental/trees/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/njh/svn/lib2geom/src/svg-path-parser.rl"
+#line 207 "/home/mental/trees/lib2geom/src/svg-path-parser.rl"
 	{
             Point p = _pop_point();
             _quadTo(_quad_tangent, p);
         }
 	break;
 	case 14:
-#line 212 "/home/njh/svn/lib2geom/src/svg-path-parser.rl"
+#line 212 "/home/mental/trees/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/njh/svn/lib2geom/src/svg-path-parser.rl"
+#line 223 "/home/mental/trees/lib2geom/src/svg-path-parser.rl"
 	{
             _closePath();
         }
 	break;
 	case 16:
-#line 360 "/home/njh/svn/lib2geom/src/svg-path-parser.rl"
+#line 360 "/home/mental/trees/lib2geom/src/svg-path-parser.rl"
 	{goto _out;}
 	break;
-#line 1566 "/home/njh/svn/lib2geom/src/svg-path-parser.cpp"
+#line 1566 "/home/mental/trees/lib2geom/src/svg-path-parser.cpp"
 		}
 	}
 
@@ -1571,7 +1571,7 @@ _again:
 	goto _resume;
 	_out: {}
 	}
-#line 370 "/home/njh/svn/lib2geom/src/svg-path-parser.rl"
+#line 370 "/home/mental/trees/lib2geom/src/svg-path-parser.rl"
 
 
     if ( cs < svg_path_first_final ) {
diff --git a/src/2geom/utils.h b/src/2geom/utils.h
index ca880640d..2be995248 100644
--- a/src/2geom/utils.h
+++ b/src/2geom/utils.h
@@ -1,5 +1,5 @@
-#ifndef MATH_UTILS_HEADER
-#define MATH_UTILS_HEADER
+#ifndef UTILS_HEADER
+#define UTILS_HEADER
 
 /** Various utility functions.
  *