update to trunk

(bzr r13708.1.39)

51 files changed, 1507 insertions, 1306 deletions

diff --git a/src/2geom/CMakeLists.txt b/src/2geom/CMakeLists.txt
index eb7abb614..97b47b630 100644
--- a/src/2geom/CMakeLists.txt
+++ b/src/2geom/CMakeLists.txt
@@ -48,7 +48,6 @@ set(2geom_SRC
 	toposweep.cpp
 	transforms.cpp
 	utils.cpp
-	viewbox.cpp
 
 
 	# -------
@@ -72,7 +71,6 @@ set(2geom_SRC
 	crossing.h
 	curve.h
 	curves.h
-	d2-sbasis.h
 	d2.h
 	ellipse.h
 	elliptical-arc.h
diff --git a/src/2geom/Makefile_insert b/src/2geom/Makefile_insert
index e3c6836fd..b56942caa 100644
--- a/src/2geom/Makefile_insert
+++ b/src/2geom/Makefile_insert
@@ -120,8 +120,6 @@
 	2geom/transforms.h \
 	2geom/utils.cpp \
 	2geom/utils.h \
-	2geom/viewbox.cpp \
-	2geom/viewbox.h \
 	2geom/numeric/fitting-model.h \
 	2geom/numeric/fitting-tool.h \
 	2geom/numeric/linear_system.h \
diff --git a/src/2geom/angle.h b/src/2geom/angle.h
index af4442e88..f0caaba6b 100644
--- a/src/2geom/angle.h
+++ b/src/2geom/angle.h
@@ -383,10 +383,10 @@ private:
 
 /** @brief Given an angle in degrees, return radians
  * @relates Angle */
-inline Coord deg_to_rad(Coord deg) { return deg*M_PI/180.0;}
+inline Coord rad_from_deg(Coord deg) { return deg*M_PI/180.0;}
 /** @brief Given an angle in radians, return degrees
  * @relates Angle */
-inline Coord rad_to_deg(Coord rad) { return rad*180.0/M_PI;}
+inline Coord deg_from_rad(Coord rad) { return rad*180.0/M_PI;}
 
 } // end namespace Geom
 
diff --git a/src/2geom/bezier-clipping.cpp b/src/2geom/bezier-clipping.cpp
index be8dd5a5f..03a6dfecd 100644
--- a/src/2geom/bezier-clipping.cpp
+++ b/src/2geom/bezier-clipping.cpp
@@ -788,7 +788,7 @@ void iterate<intersection_point_tag> (std::vector<Interval>& domsA,
 #endif
         dom = clip<intersection_point_tag>(*C1, *C2, precision);
 
-        if (dom.isEmpty())
+        if (dom.empty())
         {
 #if VERBOSE
             std::cerr << "dom: empty" << std::endl;
@@ -937,7 +937,7 @@ void iterate<collinear_normal_tag> (std::vector<Interval>& domsA,
 #endif
         dom = clip<collinear_normal_tag>(*C1, *C2, precision);
 
-        if (dom.isEmpty()) {
+        if (dom.empty()) {
 #if VERBOSE
             std::cerr << "dom: empty" << std::endl;
 #endif
diff --git a/src/2geom/bezier-curve.cpp b/src/2geom/bezier-curve.cpp
index 17221264b..866263047 100644
--- a/src/2geom/bezier-curve.cpp
+++ b/src/2geom/bezier-curve.cpp
@@ -115,6 +115,17 @@ BezierCurve::BezierCurve(std::vector<Point> const &pts)
     }
 }
 
+bool BezierCurve::isDegenerate() const
+{
+    for (unsigned d = 0; d < 2; ++d) {
+        Coord ic = inner[d][0];
+        for (unsigned i = 1; i < size(); ++i) {
+            if (inner[d][i] != ic) return false;
+        }
+    }
+    return true;
+}
+
 Coord BezierCurve::length(Coord tolerance) const
 {
     switch (order())
@@ -169,6 +180,29 @@ BezierCurve::intersect(Curve const &other, Coord eps) const
     return result;
 }
 
+bool BezierCurve::isNear(Curve const &c, Coord precision) const
+{
+    if (this == &c) return true;
+
+    BezierCurve const *other = dynamic_cast<BezierCurve const *>(&c);
+    if (!other) return false;
+
+    if (!are_near(inner.at0(), other->inner.at0(), precision)) return false;
+    if (!are_near(inner.at1(), other->inner.at1(), precision)) return false;
+
+    if (size() == other->size()) {
+        for (unsigned i = 1; i < order(); ++i) {
+            if (!are_near(inner.point(i), other->inner.point(i), precision)) {
+                return false;
+            }
+        }
+        return true;
+    } else {
+        // TODO: comparison after degree elevation
+        return false;
+    }
+}
+
 bool BezierCurve::operator==(Curve const &c) const
 {
     if (this == &c) return true;
@@ -281,6 +315,23 @@ std::vector<CurveIntersection> BezierCurveN<1>::intersect(Curve const &other, Co
 }
 
 template <>
+int BezierCurveN<1>::winding(Point const &p) const
+{
+    Point ip = inner.at0(), fp = inner.at1();
+    if (p[Y] == std::max(ip[Y], fp[Y])) return 0;
+
+    Point v = fp - ip;
+    assert(v[Y] != 0);
+    Coord t = (p[Y] - ip[Y]) / v[Y];
+    assert(t >= 0 && t <= 1);
+    Coord xcross = lerp(t, ip[X], fp[X]);
+    if (xcross > p[X]) {
+        return v[Y] > 0 ? 1 : -1;
+    }
+    return 0;
+}
+
+template <>
 void BezierCurveN<1>::feed(PathSink &sink, bool moveto_initial) const
 {
     if (moveto_initial) {
@@ -308,7 +359,7 @@ void BezierCurveN<3>::feed(PathSink &sink, bool moveto_initial) const
 }
 
 
-static Coord bezier_length_internal(std::vector<Point> &v1, Coord tolerance)
+static Coord bezier_length_internal(std::vector<Point> &v1, Coord tolerance, int level)
 {
     /* The Bezier length algorithm used in 2Geom utilizes a simple fact:
      * the Bezier curve is longer than the distance between its endpoints
@@ -317,13 +368,15 @@ static Coord bezier_length_internal(std::vector<Point> &v1, Coord tolerance)
      * error tolerance, we can be sure that the true value is no further than
      * 0.5 * tolerance from their arithmetic mean. When it's larger, we recursively
      * subdivide the Bezier curve into two parts and add their lengths.
+     * 
+     * We cap the maximum number of subdivisions at 256, which corresponds to 8 levels.
      */
     Coord lower = distance(v1.front(), v1.back());
     Coord upper = 0.0;
     for (size_t i = 0; i < v1.size() - 1; ++i) {
         upper += distance(v1[i], v1[i+1]);
     }
-    if (upper - lower < 2*tolerance) {
+    if (upper - lower <= 2*tolerance || level >= 8) {
         return (lower + upper) / 2;
     }
         
@@ -381,7 +434,8 @@ static Coord bezier_length_internal(std::vector<Point> &v1, Coord tolerance)
         v2[i] = v1[0];
     }
 
-    return bezier_length_internal(v1, 0.5*tolerance) + bezier_length_internal(v2, 0.5*tolerance);
+    return bezier_length_internal(v1, 0.5 * tolerance, level + 1) +
+           bezier_length_internal(v2, 0.5 * tolerance, level + 1);
 }
 
 /** @brief Compute the length of a bezier curve given by a vector of its control points
@@ -390,17 +444,17 @@ Coord bezier_length(std::vector<Point> const &points, Coord tolerance)
 {
     if (points.size() < 2) return 0.0;
     std::vector<Point> v1 = points;
-    return bezier_length_internal(v1, tolerance);
+    return bezier_length_internal(v1, tolerance, 0);
 }
 
-/** @brief Compute the length of a quadratic bezier curve given by its control points
- * @relatesalso QuadraticBezier */
-Coord bezier_length(Point a0, Point a1, Point a2, Coord tolerance)
+static Coord bezier_length_internal(Point a0, Point a1, Point a2, Coord tolerance, int level)
 {
     Coord lower = distance(a0, a2);
     Coord upper = distance(a0, a1) + distance(a1, a2);
 
-    if (upper - lower < 2*tolerance) return (lower + upper)/2;
+    if (upper - lower <= 2*tolerance || level >= 8) {
+        return (lower + upper) / 2;
+    }
 
     Point // Casteljau subdivision
         // b0 = a0,
@@ -408,17 +462,25 @@ Coord bezier_length(Point a0, Point a1, Point a2, Coord tolerance)
         b1 = 0.5*(a0 + a1),
         c1 = 0.5*(a1 + a2),
         b2 = 0.5*(b1 + c1); // == c2
-    return bezier_length(a0, b1, b2, 0.5*tolerance) + bezier_length(b2, c1, a2, 0.5*tolerance);
+    return bezier_length_internal(a0, b1, b2, 0.5 * tolerance, level + 1) +
+           bezier_length_internal(b2, c1, a2, 0.5 * tolerance, level + 1);
 }
 
-/** @brief Compute the length of a cubic bezier curve given by its control points
- * @relatesalso CubicBezier */
-Coord bezier_length(Point a0, Point a1, Point a2, Point a3, Coord tolerance)
+/** @brief Compute the length of a quadratic bezier curve given by its control points
+ * @relatesalso QuadraticBezier */
+Coord bezier_length(Point a0, Point a1, Point a2, Coord tolerance)
+{
+    return bezier_length_internal(a0, a1, a2, tolerance, 0);
+}
+
+static Coord bezier_length_internal(Point a0, Point a1, Point a2, Point a3, Coord tolerance, int level)
 {
     Coord lower = distance(a0, a3);
     Coord upper = distance(a0, a1) + distance(a1, a2) + distance(a2, a3);
 
-    if (upper - lower < 2*tolerance) return (lower + upper)/2;
+    if (upper - lower <= 2*tolerance || level >= 8) {
+        return (lower + upper) / 2;
+    }
 
     Point // Casteljau subdivision
         // b0 = a0,
@@ -429,7 +491,15 @@ Coord bezier_length(Point a0, Point a1, Point a2, Point a3, Coord tolerance)
         b2 = 0.5*(b1 + t0),
         c2 = 0.5*(t0 + c1),
         b3 = 0.5*(b2 + c2); // == c3
-    return bezier_length(a0, b1, b2, b3, 0.5*tolerance) + bezier_length(b3, c2, c1, a3, 0.5*tolerance);
+    return bezier_length_internal(a0, b1, b2, b3, 0.5 * tolerance, level + 1) +
+           bezier_length_internal(b3, c2, c1, a3, 0.5 * tolerance, level + 1);
+}
+
+/** @brief Compute the length of a cubic bezier curve given by its control points
+ * @relatesalso CubicBezier */
+Coord bezier_length(Point a0, Point a1, Point a2, Point a3, Coord tolerance)
+{
+    return bezier_length_internal(a0, a1, a2, a3, tolerance, 0);
 }
 
 } // end namespace Geom
diff --git a/src/2geom/bezier-curve.h b/src/2geom/bezier-curve.h
index 9ac4d7b4d..9416ba78c 100644
--- a/src/2geom/bezier-curve.h
+++ b/src/2geom/bezier-curve.h
@@ -104,7 +104,7 @@ public:
     // implementation of virtual methods goes here
     virtual Point initialPoint() const { return inner.at0(); }
     virtual Point finalPoint() const { return inner.at1(); }
-    virtual bool isDegenerate() const { return inner.isConstant(0); }
+    virtual bool isDegenerate() const;
     virtual bool isLineSegment() const { return size() == 2; }
     virtual void setInitial(Point const &v) { setPoint(0, v); }
     virtual void setFinal(Point const &v) { setPoint(order(), v); }
@@ -166,6 +166,7 @@ public:
     }
     virtual Coord valueAt(Coord t, Dim2 d) const { return inner[d].valueAt(t); }
     virtual D2<SBasis> toSBasis() const {return inner.toSBasis(); }
+    virtual bool isNear(Curve const &c, Coord precision) const;
     virtual bool operator==(Curve const &c) const;
     virtual void feed(PathSink &sink, bool) const;
 };
@@ -244,6 +245,10 @@ public:
                    BezierCurveN(sx.second, sy.second));
     }
 
+    virtual bool isDegenerate() const {
+        return BezierCurve::isDegenerate();
+    }
+
     virtual bool isLineSegment() const {
         return size() == 2;
     }
@@ -274,6 +279,9 @@ public:
         // call super. this is implemented only to allow specializations
         return BezierCurve::intersect(other, eps);
     }
+    virtual int winding(Point const &p) const {
+        return Curve::winding(p);
+    }
     virtual void feed(PathSink &sink, bool moveto_initial) const {
         // call super. this is implemented only to allow specializations
         BezierCurve::feed(sink, moveto_initial);
@@ -304,10 +312,14 @@ Curve *BezierCurveN<degree>::derivative() const {
 }
 
 // optimized specializations
+template <> inline bool BezierCurveN<1>::isDegenerate() const {
+    return inner[X][0] == inner[X][1] && inner[Y][0] == inner[Y][1];
+}
 template <> inline bool BezierCurveN<1>::isLineSegment() const { return true; }
 template <> Curve *BezierCurveN<1>::derivative() const;
 template <> Coord BezierCurveN<1>::nearestTime(Point const &, Coord, Coord) const;
 template <> std::vector<CurveIntersection> BezierCurveN<1>::intersect(Curve const &, Coord) const;
+template <> int BezierCurveN<1>::winding(Point const &) const;
 template <> void BezierCurveN<1>::feed(PathSink &sink, bool moveto_initial) const;
 template <> void BezierCurveN<2>::feed(PathSink &sink, bool moveto_initial) const;
 template <> void BezierCurveN<3>::feed(PathSink &sink, bool moveto_initial) const;
diff --git a/src/2geom/circle.cpp b/src/2geom/circle.cpp
index 553981a72..934a8d3ab 100644
--- a/src/2geom/circle.cpp
+++ b/src/2geom/circle.cpp
@@ -144,7 +144,7 @@ bool Circle::contains(Circle const &other) const
 bool Circle::intersects(Line const &l) const
 {
     // http://mathworld.wolfram.com/Circle-LineIntersection.html
-    Coord dr = l.versor().length();
+    Coord dr = l.vector().length();
     Coord r = _radius;
     Coord D = cross(l.initialPoint(), l.finalPoint());
     Coord delta = r*r * dr*dr - D*D;
@@ -163,9 +163,9 @@ bool Circle::intersects(Circle const &other) const
 std::vector<ShapeIntersection> Circle::intersect(Line const &l) const
 {
     // http://mathworld.wolfram.com/Circle-LineIntersection.html
-    Coord dr = l.versor().length();
-    Coord dx = l.versor().x();
-    Coord dy = l.versor().y();
+    Coord dr = l.vector().length();
+    Coord dx = l.vector().x();
+    Coord dy = l.vector().y();
     Coord D = cross(l.initialPoint() - _center, l.finalPoint() - _center);
     Coord delta = _radius*_radius * dr*dr - D*D;
 
diff --git a/src/2geom/conicsec.cpp b/src/2geom/conicsec.cpp
index 089db71a4..2d396ba30 100644
--- a/src/2geom/conicsec.cpp
+++ b/src/2geom/conicsec.cpp
@@ -114,8 +114,8 @@ double RatQuad::lambda() const {
 RatQuad RatQuad::fromPointsTangents(Point P0, Point dP0,
                        Point P,
                        Point P2, Point dP2) {
-  Line Line0 = Line::from_origin_and_versor(P0, dP0);
-  Line Line2 = Line::from_origin_and_versor(P2, dP2);
+  Line Line0 = Line::from_origin_and_vector(P0, dP0);
+  Line Line2 = Line::from_origin_and_vector(P2, dP2);
   try {
     OptCrossing oc = intersection(Line0, Line2);
     if(!oc) // what to do?
@@ -276,8 +276,8 @@ std::vector<Point> decompose_degenerate(xAx const & C1, xAx const & C2, xAx cons
             n1 = Point(b-d, 1);
         }
 
-        Line L0 = Line::from_origin_and_versor(B0, rot90(n0));
-        Line L1 = Line::from_origin_and_versor(B0, rot90(n1));
+        Line L0 = Line::from_origin_and_vector(B0, rot90(n0));
+        Line L1 = Line::from_origin_and_vector(B0, rot90(n1));
 
         std::vector<double> rts = C1.roots(L0);
         for(unsigned i = 0; i < rts.size(); i++) {
@@ -327,12 +327,12 @@ std::vector<Point> decompose_degenerate(xAx const & C1, xAx const & C2, xAx cons
          */
         assert(L2sq(g) != 0);
 
-        Line Lx = Line::from_origin_and_versor(trial_pt, g); // a line along the gradient
+        Line Lx = Line::from_origin_and_vector(trial_pt, g); // a line along the gradient
         std::vector<double> rts = xC0.roots(Lx);
         for(unsigned i = 0; i < rts.size(); i++) {
             Point P0 = Lx.pointAt(rts[i]);
             //std::cout << P0 << "\n";
-            Line L = Line::from_origin_and_versor(P0, rot90(g));
+            Line L = Line::from_origin_and_vector(P0, rot90(g));
             std::vector<double> cnrts;
             // It's very likely that at least one of the conics is degenerate, this will hopefully pick the more generate of the two.
             if(fabs(C1.hessian().det()) > fabs(C2.hessian().det()))
@@ -515,7 +515,7 @@ xAx xAx::operator*(double const &b) const {
       if(L2sq(dA) <= 1e-10) { // perhaps a single point?
           return boost::optional<RatQuad> ();
       }
-      LineSegment ls = intersection(Line::from_origin_and_versor(A, dA), bnd);
+      LineSegment ls = intersection(Line::from_origin_and_vector(A, dA), bnd);
       return RatQuad::fromPointsTangents(A, dA, ls.pointAt(0.5), ls[1], dA);
   }
   else if(crs.size() >= 2 && crs.size() < 4) {
diff --git a/src/2geom/coord.cpp b/src/2geom/coord.cpp
index 8b5e28586..8c7485883 100644
--- a/src/2geom/coord.cpp
+++ b/src/2geom/coord.cpp
@@ -111,12 +111,6 @@ namespace Geom {
 
 namespace {
 
-inline int StrLength(const char* string) {
-  size_t length = strlen(string);
-  ASSERT(length == static_cast<size_t>(static_cast<int>(length)));
-  return static_cast<int>(length);
-}
-
 template <typename T>
 class Vector {
  public:
@@ -1733,8 +1727,6 @@ enum FastDtoaMode {
   FAST_DTOA_PRECISION
 };
 
-static const int kFastDtoaMaximalLength = 17;
-
 bool FastDtoa(double d,
               FastDtoaMode mode,
               int requested_digits,
@@ -2350,7 +2342,6 @@ bool FastFixedDtoa(double v,
   return true;
 }
 
-static const int kMaxExactDoubleIntegerDecimalDigits = 15;
 static const int kMaxUint64DecimalDigits = 19;
 
 static const int kMaxDecimalPower = 309;
diff --git a/src/2geom/crossing.h b/src/2geom/crossing.h
index 425fa58f5..0f007b192 100644
--- a/src/2geom/crossing.h
+++ b/src/2geom/crossing.h
@@ -119,7 +119,7 @@ struct CrossingOrder {
                    (ix == b.a ? b.ta : b.tb);
         else
             return (ix == a.a ? a.ta : a.tb) >
-                   (ix == a.a ? a.ta : a.tb);
+                   (ix == b.a ? b.ta : b.tb);
     }
 };
 
diff --git a/src/2geom/curve.cpp b/src/2geom/curve.cpp
index b45228514..387a9180b 100644
--- a/src/2geom/curve.cpp
+++ b/src/2geom/curve.cpp
@@ -69,14 +69,14 @@ int Curve::winding(Point const &p) const
         // skip endpoint roots when they are local maxima on the Y axis
         // this follows the convention used in other winding routines,
         // i.e. that the bottommost coordinate is not part of the shape
-        bool ingore_0 = unitTangentAt(0)[Y] <= 0;
+        bool ignore_0 = unitTangentAt(0)[Y] <= 0;
         bool ignore_1 = unitTangentAt(1)[Y] >= 0;
 
         int wind = 0;
         for (std::size_t i = 0; i < ts.size(); ++i) {
             Coord t = ts[i];
             //std::cout << t << std::endl;
-            if ((t == 0 && ingore_0) || (t == 1 && ignore_1)) continue;
+            if ((t == 0 && ignore_0) || (t == 1 && ignore_1)) continue;
             if (valueAt(t, X) > p[X]) { // root is ray intersection
                 Point tangent = unitTangentAt(t);
                 if (tangent[Y] > 0) {
diff --git a/src/2geom/curve.h b/src/2geom/curve.h
index abbdb1100..41db9ca8a 100644
--- a/src/2geom/curve.h
+++ b/src/2geom/curve.h
@@ -333,6 +333,9 @@ public:
      * @return True if the curves are identical, false otherwise */
     virtual bool operator==(Curve const &c) const = 0;
 
+    /** @brief Test whether two curves are approximately the same. */
+    virtual bool isNear(Curve const &c, Coord precision) const = 0;
+
     /** @brief Feed the curve to a PathSink */
     virtual void feed(PathSink &sink, bool moveto_initial) const;
     /// @}
diff --git a/src/2geom/d2-sbasis.cpp b/src/2geom/d2-sbasis.cpp
index ebec16fdd..07eccce76 100644
--- a/src/2geom/d2-sbasis.cpp
+++ b/src/2geom/d2-sbasis.cpp
@@ -1,7 +1,41 @@
+/**
+ * \file
+ * \brief Some two-dimensional SBasis operations
+ *//*
+ * Authors:
+ *   MenTaLguy <mental@rydia.net>
+ *   Jean-François Barraud <jf.barraud@gmail.com>
+ *   Johan Engelen <j.b.c.engelen@alumnus.utwente.nl>
+ *   
+ * Copyright 2007-2012 Authors
+ *
+ * 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.
+ *
+ */
+
 #include <2geom/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. */
+#include <2geom/piecewise.h>
 
 namespace Geom {
 
@@ -147,12 +181,12 @@ Piecewise<D2<SBasis> > force_continuity(Piecewise<D2<SBasis> > const &f, double
                 SBasis &prev_sb=result.segs[prev][dim];
                 SBasis &cur_sb =result.segs[cur][dim];
                 Coord const c=pt0[dim];
-                if (prev_sb.empty()) {
+                if (prev_sb.isZero(0)) {
                   prev_sb = SBasis(Linear(0.0, c));
                 } else {
                   prev_sb[0][1] = c;
                 }
-                if (cur_sb.empty()) {
+                if (cur_sb.isZero(0)) {
                   cur_sb = SBasis(Linear(c, 0.0));
                 } else {
                   cur_sb[0][0] = c;
@@ -198,30 +232,22 @@ Point unitTangentAt(D2<SBasis> const & a, Coord t, unsigned n)
     return Point (0,0);
 }
 
-static void set_first_point(Piecewise<D2<SBasis> > &f, Point a){
+static void set_first_point(Piecewise<D2<SBasis> > &f, Point const &a){
     if ( f.empty() ){
         f.concat(Piecewise<D2<SBasis> >(D2<SBasis>(SBasis(Linear(a[X])), SBasis(Linear(a[Y])))));
         return;
     }
     for (unsigned dim=0; dim<2; dim++){
-        if (f.segs.front()[dim].size() == 0){
-            f.segs.front()[dim] = SBasis(Linear(a[dim],0));
-        }else{
-            f.segs.front()[dim][0][0] = a[dim];
-        }
+        f.segs.front()[dim][0][0] = a[dim];
     }
 }
-static void set_last_point(Piecewise<D2<SBasis> > &f, Point a){
+static void set_last_point(Piecewise<D2<SBasis> > &f, Point const &a){
     if ( f.empty() ){
         f.concat(Piecewise<D2<SBasis> >(D2<SBasis>(SBasis(Linear(a[X])), SBasis(Linear(a[Y])))));
         return;
     }
     for (unsigned dim=0; dim<2; dim++){
-        if (f.segs.back()[dim].size() == 0){
-            f.segs.back()[dim] = SBasis(Linear(0,a[dim]));
-        }else{
-            f.segs.back()[dim][0][1] = a[dim];
-        }
+        f.segs.back()[dim][0][1] = a[dim];
     }
 }
 
diff --git a/src/2geom/d2-sbasis.h b/src/2geom/d2-sbasis.h
deleted file mode 100644
index a0769b314..000000000
--- a/src/2geom/d2-sbasis.h
+++ /dev/null
@@ -1,165 +0,0 @@
-/**
- * \file
- * \brief  Do not include this file
- *
- * We don't actually want anyone to
- * include this, other than D2.h.
- *//*
- * Authors:
- *      ? <?@?.?>
- * 
- * Copyright ?-?  authors
- *
- * 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.
- *
- */
-
-#ifdef LIB2GEOM_SEEN_D2_H  /*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 LIB2GEOM_SEEN_D2_SBASIS_H
-#define LIB2GEOM_SEEN_D2_SBASIS_H
-
-#include <2geom/sbasis.h>
-#include <2geom/sbasis-2d.h>
-#include <2geom/piecewise.h>
-#include <2geom/affine.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_independent(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> dot(Piecewise<D2<SBasis> > const &a, Point const &b);
-Piecewise<SBasis> cross(Piecewise<D2<SBasis> > const &a, Piecewise<D2<SBasis> > const &b);
-
-Piecewise<D2<SBasis> > operator*(Piecewise<D2<SBasis> > const &a, Affine const &m);
-
-Piecewise<D2<SBasis> > force_continuity(Piecewise<D2<SBasis> > const &f, double tol=0, bool closed=false);
-
-std::vector<Piecewise<D2<SBasis> > > fuse_nearby_ends(std::vector<Piecewise<D2<SBasis> > > const &f, double tol=0);
-
-std::vector<Geom::Piecewise<Geom::D2<Geom::SBasis> > > split_at_discontinuities (Geom::Piecewise<Geom::D2<Geom::SBasis> > const & pwsbin, double tol = .0001);
-
-/**
-* note that -unitTangentAt(reverse(a),0.) == unitTangentAt(a,1.) but the former is more reliable (sometimes the sign is wrong for the latter)
-*/
-Point unitTangentAt(D2<SBasis> const & a, Coord t, unsigned n = 3);
-
-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 OptRect bounds_fast(D2<SBasis> const & s, unsigned order=0) {
-    OptRect retval;
-    OptInterval xint = bounds_fast(s[X], order);
-    if (xint) {
-        OptInterval yint = bounds_fast(s[Y], order);
-        if (yint) {
-            retval = Rect(*xint, *yint);
-        }
-    }
-    return retval;
-}
-inline OptRect bounds_local(D2<SBasis> const & s, OptInterval i, unsigned order=0) {
-    OptRect retval;
-    OptInterval xint = bounds_local(s[X], i, order);
-    OptInterval yint = bounds_local(s[Y], i, order);
-    if (xint && yint) {
-        retval = Rect(*xint, *yint);
-    }
-    return retval;
-}
-
-std::vector<Interval> level_set( D2<SBasis> const &f, Rect region);
-std::vector<Interval> level_set( D2<SBasis> const &f, Point p, double tol);
-std::vector<std::vector<Interval> > level_sets( D2<SBasis> const &f, std::vector<Rect> regions);
-std::vector<std::vector<Interval> > level_sets( D2<SBasis> const &f, std::vector<Point> pts, double tol);
-
-}
-
-#endif // LIB2GEOM_SEEN_D2_SBASIS_H
-#endif // LIB2GEOM_SEEN_D2_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:fileencoding=utf-8:textwidth=99 :
diff --git a/src/2geom/d2.h b/src/2geom/d2.h
index dca6fa614..cd67ec65b 100644
--- a/src/2geom/d2.h
+++ b/src/2geom/d2.h
@@ -1,8 +1,12 @@
 /**
  * \file
- * \brief   Lifts one dimensional objects into 2d 
+ * \brief Lifts one dimensional objects into 2D
  *//*
- * Copyright 2007 Michael Sloan <mgsloan@gmail.com>
+ * Authors:
+ *   Michael Sloan <mgsloan@gmail.com>
+ *   Krzysztof Kosiński <tweenk.pl@gmail.com>
+ *
+ * Copyright 2007-2015 Authors
  *
  * This library is free software; you can redistribute it and/or
  * modify it either under the terms of the GNU Lesser General Public
@@ -461,12 +465,6 @@ inline std::ostream &operator<< (std::ostream &out_file, const Geom::D2<T> &in_d
     return out_file;
 }
 
-} //end namespace Geom
-
-#include <2geom/d2-sbasis.h>
-
-namespace Geom {
-
 //Some D2 Fragment implementation which requires rect:
 template <typename T>
 OptRect bounds_fast(const D2<T> &a) {
@@ -484,6 +482,73 @@ OptRect bounds_local(const D2<T> &a, const OptInterval &t) {
     return OptRect(bounds_local(a[X], t), bounds_local(a[Y], t));
 }
 
+
+
+// SBasis-specific declarations
+
+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 declarations
+
+Piecewise<D2<SBasis> > sectionize(D2<Piecewise<SBasis> > const &a);
+D2<Piecewise<SBasis> > make_cuts_independent(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> dot(Piecewise<D2<SBasis> > const &a, Point const &b);
+Piecewise<SBasis> cross(Piecewise<D2<SBasis> > const &a, Piecewise<D2<SBasis> > const &b);
+
+Piecewise<D2<SBasis> > operator*(Piecewise<D2<SBasis> > const &a, Affine const &m);
+
+Piecewise<D2<SBasis> > force_continuity(Piecewise<D2<SBasis> > const &f, double tol=0, bool closed=false);
+
+std::vector<Piecewise<D2<SBasis> > > fuse_nearby_ends(std::vector<Piecewise<D2<SBasis> > > const &f, double tol=0);
+
+std::vector<Geom::Piecewise<Geom::D2<Geom::SBasis> > > split_at_discontinuities (Geom::Piecewise<Geom::D2<Geom::SBasis> > const & pwsbin, double tol = .0001);
+
+Point unitTangentAt(D2<SBasis> const & a, Coord t, unsigned n = 3);
+
+//bounds specializations with order
+inline OptRect bounds_fast(D2<SBasis> const & s, unsigned order=0) {
+    OptRect retval;
+    OptInterval xint = bounds_fast(s[X], order);
+    if (xint) {
+        OptInterval yint = bounds_fast(s[Y], order);
+        if (yint) {
+            retval = Rect(*xint, *yint);
+        }
+    }
+    return retval;
+}
+inline OptRect bounds_local(D2<SBasis> const & s, OptInterval i, unsigned order=0) {
+    OptRect retval;
+    OptInterval xint = bounds_local(s[X], i, order);
+    OptInterval yint = bounds_local(s[Y], i, order);
+    if (xint && yint) {
+        retval = Rect(*xint, *yint);
+    }
+    return retval;
+}
+
+std::vector<Interval> level_set( D2<SBasis> const &f, Rect region);
+std::vector<Interval> level_set( D2<SBasis> const &f, Point p, double tol);
+std::vector<std::vector<Interval> > level_sets( D2<SBasis> const &f, std::vector<Rect> regions);
+std::vector<std::vector<Interval> > level_sets( D2<SBasis> const &f, std::vector<Point> pts, double tol);
+
+
 } // end namespace Geom
 
 #endif
diff --git a/src/2geom/ellipse.cpp b/src/2geom/ellipse.cpp
index 4b5ad9762..f4704460f 100644
--- a/src/2geom/ellipse.cpp
+++ b/src/2geom/ellipse.cpp
@@ -142,6 +142,24 @@ LineSegment Ellipse::semiaxis(Dim2 d, int sign) const
     return ls;
 }
 
+Rect Ellipse::boundsExact() const
+{
+    Angle extremes[2][2];
+    double sinrot, cosrot;
+    sincos(_angle, sinrot, cosrot);
+
+    extremes[X][0] = std::atan2( -ray(Y) * sinrot, ray(X) * cosrot );
+    extremes[X][1] = extremes[X][0] + M_PI;
+    extremes[Y][0] = std::atan2( ray(Y) * cosrot, ray(X) * sinrot );
+    extremes[Y][1] = extremes[Y][0] + M_PI;
+
+    Rect result;
+    for (unsigned d = 0; d < 2; ++d) {
+        result[d] = Interval(valueAt(extremes[d][0], d ? Y : X),
+                             valueAt(extremes[d][1], d ? Y : X));
+    }
+    return result;
+}
 
 std::vector<double> Ellipse::coefficients() const
 {
@@ -396,45 +414,39 @@ std::vector<ShapeIntersection> Ellipse::intersect(Line const &line) const
     coefficients(A, B, C, D, E, F);
     Affine iuct = inverseUnitCircleTransform();
 
-    if (line.isHorizontal()) {
-        // substitute y into the ellipse equation and solve
-        Coord y = line.initialPoint()[Y];
-        std::vector<Coord> xs = solve_quadratic(A, B*y + D, C*y*y + E*y + F);
-        for (unsigned i = 0; i < xs.size(); ++i) {
-            Point p(xs[i], y);
-            result.push_back(ShapeIntersection(atan2(p * iuct), line.timeAt(p), p));
-        }
-        return result;
-    }
-    if (line.isVertical()) {
-        // substitute y into the ellipse equation and solve
-        Coord x = line.initialPoint()[X];
-        std::vector<Coord> ys = solve_quadratic(C, B*x + E, A*x*x + D*x + F);
-        for (unsigned i = 0; i < ys.size(); ++i) {
-            Point p(x, ys[i]);
-            result.push_back(ShapeIntersection(atan2(p * iuct), line.timeAt(p), p));
-        }
-        return result;
-    }
-
     // generic case
     Coord a, b, c;
     line.coefficients(a, b, c);
+    Point lv = line.vector();
 
-    // y = -a/b x - C/B
-    // TODO: when is it better to substitute X?
-    Coord q = -a/b;
-    Coord r = -c/b;
+    if (fabs(lv[X]) > fabs(lv[Y])) {
+        // y = -a/b x - c/b
+        Coord q = -a/b;
+        Coord r = -c/b;
 
-    // substitute that into the ellipse equation, making it quadratic in x
-    Coord I = A + B*q + C*q*q;          // x^2 terms
-    Coord J = B*r + C*2*q*r + D + E*q;  // x^1 terms
-    Coord K = C*r*r + E*r + F;          // x^0 terms
-    std::vector<Coord> xs = solve_quadratic(I, J, K);
+        // substitute that into the ellipse equation, making it quadratic in x
+        Coord I = A + B*q + C*q*q;          // x^2 terms
+        Coord J = B*r + C*2*q*r + D + E*q;  // x^1 terms
+        Coord K = C*r*r + E*r + F;          // x^0 terms
+        std::vector<Coord> xs = solve_quadratic(I, J, K);
 
-    for (unsigned i = 0; i < xs.size(); ++i) {
-        Point p(xs[i], q*xs[i] + r);
-        result.push_back(ShapeIntersection(atan2(p * iuct), line.timeAt(p), p));
+        for (unsigned i = 0; i < xs.size(); ++i) {
+            Point p(xs[i], q*xs[i] + r);
+            result.push_back(ShapeIntersection(atan2(p * iuct), line.timeAt(p), p));
+        }
+    } else {
+        Coord q = -b/a;
+        Coord r = -c/a;
+
+        Coord I = A*q*q + B*q + C;
+        Coord J = A*2*q*r + B*r + D*q + E;
+        Coord K = A*r*r + D*r + F;
+        std::vector<Coord> xs = solve_quadratic(I, J, K);
+
+        for (unsigned i = 0; i < xs.size(); ++i) {
+            Point p(q*xs[i] + r, xs[i]);
+            result.push_back(ShapeIntersection(atan2(p * iuct), line.timeAt(p), p));
+        }
     }
     return result;
 }
diff --git a/src/2geom/ellipse.h b/src/2geom/ellipse.h
index 4b1701cce..f6089c944 100644
--- a/src/2geom/ellipse.h
+++ b/src/2geom/ellipse.h
@@ -171,6 +171,9 @@ public:
     LineSegment axis(Dim2 d) const;
     LineSegment semiaxis(Dim2 d, int sign = 1) const;
 
+    /// Get the tight-fitting bounding box of the ellipse.
+    Rect boundsExact() const;
+
     /// Get the coefficients of the ellipse's implicit equation.
     std::vector<double> coefficients() const;
     void coefficients(Coord &A, Coord &B, Coord &C, Coord &D, Coord &E, Coord &F) const;
diff --git a/src/2geom/elliptical-arc.cpp b/src/2geom/elliptical-arc.cpp
index 25a78b4ad..a0e379a21 100644
--- a/src/2geom/elliptical-arc.cpp
+++ b/src/2geom/elliptical-arc.cpp
@@ -92,47 +92,58 @@ namespace Geom
  * @ingroup Curves
  */
 
+
+/** @brief Compute bounds of an elliptical arc.
+ * The bounds computation works as follows. The extreme X and Y points
+ * are either the endpoints or local minima / maxima of the ellipse.
+ * We already have endpoints, and we can find the local extremes
+ * by computing a partial derivative with respect to the angle
+ * and equating that to zero:
+ * \f{align*}{
+     x &= r_x \cos \varphi \cos \theta - r_y \sin \varphi \sin \theta + c_x \\ 
+     \frac{\partial x}{\partial \theta} &= -r_x \cos \varphi \sin \theta - r_y \sin \varphi \cos \theta = 0 \\ 
+     \frac{\sin \theta}{\cos \theta} &= \tan\theta = -\frac{r_y \sin \varphi}{r_x \cos \varphi} \\ 
+     \theta &= \tan^{-1} \frac{-r_y \sin \varphi}{r_x \cos \varphi}
+   \f}
+ * The local extremes correspond to two angles separated by \f$\pi\f$.
+ * Once we compute these angles, we check whether they belong to the arc,
+ * and if they do, we evaluate the ellipse at these angles.
+ * The bounding box of the arc is equal to the bounding box of the endpoints
+ * and the local extrema that belong to the arc.
+ */
 Rect EllipticalArc::boundsExact() const
 {
     if (isChord()) {
         return chord().boundsExact();
     }
 
-    using std::swap;
+    Coord coord[2][4] = {
+        { _initial_point[X], _final_point[X], 0, 0 },
+        { _initial_point[Y], _final_point[Y], 0, 0 }
+    };
+    int ncoord[2] = { 2, 2 };
 
-    // TODO: simplify / document what is going on here.
-    double extremes[4];
+    Angle extremes[2][2];
     double sinrot, cosrot;
     sincos(rotationAngle(), sinrot, cosrot);
 
-    extremes[0] = std::atan2( -ray(Y) * sinrot, ray(X) * cosrot );
-    extremes[1] = extremes[0] + M_PI;
-    if ( extremes[0] < 0 ) extremes[0] += 2*M_PI;
-    extremes[2] = std::atan2( ray(Y) * cosrot, ray(X) * sinrot );
-    extremes[3] = extremes[2] + M_PI;
-    if ( extremes[2] < 0 ) extremes[2] += 2*M_PI;
-
-
-    double arc_extremes[4];
-    arc_extremes[0] = initialPoint()[X];
-    arc_extremes[1] = finalPoint()[X];
-    if ( arc_extremes[0] < arc_extremes[1] )
-        swap(arc_extremes[0], arc_extremes[1]);
-    arc_extremes[2] = initialPoint()[Y];
-    arc_extremes[3] = finalPoint()[Y];
-    if ( arc_extremes[2] < arc_extremes[3] )
-        swap(arc_extremes[2], arc_extremes[3]);
-
-    if ( !are_near(initialPoint(), finalPoint()) ) {
-        for (unsigned i = 0; i < 4; ++i) {
-            if (containsAngle(extremes[i])) {
-                arc_extremes[i] = valueAtAngle(extremes[i], (i >> 1) ? Y : X);
+    extremes[X][0] = std::atan2( -ray(Y) * sinrot, ray(X) * cosrot );
+    extremes[X][1] = extremes[X][0] + M_PI;
+    extremes[Y][0] = std::atan2( ray(Y) * cosrot, ray(X) * sinrot );
+    extremes[Y][1] = extremes[Y][0] + M_PI;
+
+    for (unsigned d = 0; d < 2; ++d) {
+        for (unsigned i = 0; i < 2; ++i) {
+            if (containsAngle(extremes[d][i])) {
+                coord[d][ncoord[d]++] = valueAtAngle(extremes[d][i], d ? Y : X);
             }
         }
     }
 
-    return Rect( Point(arc_extremes[1], arc_extremes[3]) ,
-                 Point(arc_extremes[0], arc_extremes[2]) );
+    Interval xival = Interval::from_range(coord[X], coord[X] + ncoord[X]);
+    Interval yival = Interval::from_range(coord[Y], coord[Y] + ncoord[Y]);
+    Rect result(xival, yival);
+    return result;
 }
 
 
@@ -149,116 +160,14 @@ Coord EllipticalArc::valueAtAngle(Coord t, Dim2 d) const
 
 std::vector<Coord> EllipticalArc::roots(Coord v, Dim2 d) const
 {
-    if (isChord()) {
-        return chord().roots(v, d);
-    }
-
     std::vector<Coord> sol;
-    Interval unit_interval(0, 1);
 
-    if ( are_near(ray(X), 0) && are_near(ray(Y), 0) ) {
-        if ( center(d) == v )
-            sol.push_back(0);
+    if (isChord()) {
+        sol = chord().roots(v, d);
         return sol;
     }
 
-    static const char* msg[2][2] =
-    {
-        { "d == X; ray(X) == 0; "
-          "s = (v - center(X)) / ( -ray(Y) * std::sin(_rot_angle) ); "
-          "s should be contained in [-1,1]",
-          "d == X; ray(Y) == 0; "
-          "s = (v - center(X)) / ( ray(X) * std::cos(_rot_angle) ); "
-          "s should be contained in [-1,1]"
-        },
-        { "d == Y; ray(X) == 0; "
-          "s = (v - center(X)) / ( ray(Y) * std::cos(_rot_angle) ); "
-          "s should be contained in [-1,1]",
-          "d == Y; ray(Y) == 0; "
-          "s = (v - center(X)) / ( ray(X) * std::sin(_rot_angle) ); "
-          "s should be contained in [-1,1]"
-        },
-    };
-
-    for ( unsigned int dim = 0; dim < 2; ++dim )
-    {
-        if (ray((Dim2) dim) == 0)
-        {
-            if ( initialPoint()[d] == v && finalPoint()[d] == v )
-            {
-                THROW_INFINITESOLUTIONS(0);
-            }
-            if ( (initialPoint()[d] < finalPoint()[d])
-                 && (initialPoint()[d] > v || finalPoint()[d] < v) )
-            {
-                return sol;
-            }
-            if ( (initialPoint()[d] > finalPoint()[d])
-                 && (finalPoint()[d] > v || initialPoint()[d] < v) )
-            {
-                return sol;
-            }
-            double ray_prj = 0.0;
-            switch(d)
-            {
-                case X:
-                    switch(dim)
-                    {
-                        case X: ray_prj = -ray(Y) * std::sin(rotationAngle());
-                                break;
-                        case Y: ray_prj = ray(X) * std::cos(rotationAngle());
-                                break;
-                    }
-                    break;
-                case Y:
-                    switch(dim)
-                    {
-                        case X: ray_prj = ray(Y) * std::cos(rotationAngle());
-                                break;
-                        case Y: ray_prj = ray(X) * std::sin(rotationAngle());
-                                break;
-                    }
-                    break;
-            }
-
-            double s = (v - center(d)) / ray_prj;
-            if ( s < -1 || s > 1 )
-            {
-                THROW_LOGICALERROR(msg[d][dim]);
-            }
-            switch(dim)
-            {
-                case X:
-                    s = std::asin(s); // return a value in [-PI/2,PI/2]
-                    if ( logical_xor( sweep(), are_near(initialAngle(), M_PI/2) )  )
-                    {
-                        if ( s < 0 ) s += 2*M_PI;
-                    }
-                    else
-                    {
-                        s = M_PI - s;
-                        if (!(s < 2*M_PI) ) s -= 2*M_PI;
-                    }
-                    break;
-                case Y:
-                    s = std::acos(s); // return a value in [0,PI]
-                    if ( logical_xor( sweep(), are_near(initialAngle(), 0) ) )
-                    {
-                        s = 2*M_PI - s;
-                        if ( !(s < 2*M_PI) ) s -= 2*M_PI;
-                    }
-                    break;
-            }
-
-            //std::cerr << "s = " << rad_to_deg(s);
-            s = timeAtAngle(s);
-            //std::cerr << " -> t: " << s << std::endl;
-            if (unit_interval.contains(s)) {
-                sol.push_back(s);
-            }
-            return sol;
-        }
-    }
+    Interval unit_interval(0, 1);
 
     double rotx, roty;
     if (d == X) {
@@ -278,10 +187,10 @@ std::vector<Coord> EllipticalArc::roots(Coord v, Dim2 d) const
     //std::cerr << "b = " << b << std::endl;
     //std::cerr << "c = " << c << std::endl;
 
-    if ( are_near(a,0) )
+    if (a == 0)
     {
         sol.push_back(M_PI);
-        if ( !are_near(b,0) )
+        if (b != 0)
         {
             double s = 2 * std::atan(-c/(2*b));
             if ( s < 0 ) s += 2*M_PI;
@@ -292,8 +201,7 @@ std::vector<Coord> EllipticalArc::roots(Coord v, Dim2 d) const
     {
         double delta = b * b - a * c;
         //std::cerr << "delta = " << delta << std::endl;
-        if ( are_near(delta, 0) )
-        {
+        if (delta == 0) {
             double s = 2 * std::atan(-b/a);
             if ( s < 0 ) s += 2*M_PI;
             sol.push_back(s);
@@ -312,7 +220,7 @@ std::vector<Coord> EllipticalArc::roots(Coord v, Dim2 d) const
 
     std::vector<double> arc_sol;
     for (unsigned int i = 0; i < sol.size(); ++i ) {
-        //std::cerr << "s = " << rad_to_deg(sol[i]);
+        //std::cerr << "s = " << deg_from_rad(sol[i]);
         sol[i] = timeAtAngle(sol[i]);
         //std::cerr << " -> t: " << sol[i] << std::endl;
         if (unit_interval.contains(sol[i])) {
@@ -886,6 +794,25 @@ bool EllipticalArc::operator==(Curve const &c) const
     return true;
 }
 
+bool EllipticalArc::isNear(Curve const &c, Coord precision) const
+{
+    EllipticalArc const *other = dynamic_cast<EllipticalArc const *>(&c);
+    if (!other) {
+        if (isChord()) {
+            return c.isNear(chord(), precision);
+        }
+        return false;
+    }
+
+    if (!are_near(_initial_point, other->_initial_point, precision)) return false;
+    if (!are_near(_final_point, other->_final_point, precision)) return false;
+    if (isChord() && other->isChord()) return true;
+
+    if (sweep() != other->sweep()) return false;
+    if (!are_near(_ellipse, other->_ellipse, precision)) return false;
+    return true;
+}
+
 void EllipticalArc::feed(PathSink &sink, bool moveto_initial) const
 {
     if (moveto_initial) {
@@ -894,6 +821,104 @@ void EllipticalArc::feed(PathSink &sink, bool moveto_initial) const
     sink.arcTo(ray(X), ray(Y), rotationAngle(), _large_arc, sweep(), _final_point);
 }
 
+int EllipticalArc::winding(Point const &p) const
+{
+    using std::swap;
+
+    double sinrot, cosrot;
+    sincos(rotationAngle(), sinrot, cosrot);
+
+    Angle ymin_a = std::atan2( ray(Y) * cosrot, ray(X) * sinrot );
+    Angle ymax_a = ymin_a + M_PI;
+
+    Point ymin = pointAtAngle(ymin_a);
+    Point ymax = pointAtAngle(ymax_a);
+    if (ymin[Y] > ymax[Y]) {
+        swap(ymin, ymax);
+        swap(ymin_a, ymax_a);
+    }
+
+    Interval yspan(ymin[Y], ymax[Y]);
+    if (!yspan.lowerContains(p[Y])) return 0;
+
+    bool left = cross(ymax - ymin, p - ymin) > 0;
+    bool inside = _ellipse.contains(p);
+    bool includes_ymin = _angles.contains(ymin_a);
+    bool includes_ymax = _angles.contains(ymax_a);
+
+    AngleInterval rarc(ymin_a, ymax_a, true),
+                  larc(ymax_a, ymin_a, true);
+
+    // we'll compute the result for an arc in the direction of increasing angles
+    // and then negate if necessary
+    Angle ia = initialAngle(), fa = finalAngle();
+    Point ip = _initial_point, fp = _final_point;
+    if (!sweep()) {
+        swap(ia, fa);
+        swap(ip, fp);
+    }
+
+    bool initial_left = larc.contains(ia);
+    bool initial_right = !initial_left; //rarc.contains(ia);
+    bool final_right = rarc.contains(fa);
+    bool final_left = !final_right;//larc.contains(fa);
+
+    int result = 0;
+    if (inside || left) {
+        if (includes_ymin && final_right) {
+            Interval ival(ymin[Y], fp[Y]);
+            if (ival.lowerContains(p[Y])) {
+                ++result;
+            }
+        }
+        if (initial_right && final_right && !largeArc()) {
+            Interval ival(ip[Y], fp[Y]);
+            if (ival.lowerContains(p[Y])) {
+                ++result;
+            }
+        }
+        if (initial_right && includes_ymax) {
+            Interval ival(ip[Y], ymax[Y]);
+            if (ival.lowerContains(p[Y])) {
+                ++result;
+            }
+        }
+        if (!initial_right && !final_right && includes_ymin && includes_ymax) {
+            Interval ival(ymax[Y], ymin[Y]);
+            if (ival.lowerContains(p[Y])) {
+                ++result;
+            }
+        }
+    }
+    if (left && !inside) {
+        if (includes_ymin && initial_left) {
+            Interval ival(ymin[Y], ip[Y]);
+            if (ival.lowerContains(p[Y])) {
+                --result;
+            }
+        }
+        if (initial_left && final_left && !largeArc()) {
+            Interval ival(ip[Y], fp[Y]);
+            if (ival.lowerContains(p[Y])) {
+                --result;
+            }
+        }
+        if (final_left && includes_ymax) {
+            Interval ival(fp[Y], ymax[Y]);
+            if (ival.lowerContains(p[Y])) {
+                --result;
+            }
+        }
+        if (!initial_left && !final_left && includes_ymin && includes_ymax) {
+            Interval ival(ymax[Y], ymin[Y]);
+            if (ival.lowerContains(p[Y])) {
+                --result;
+            }
+        }
+    }
+    return sweep() ? result : -result;
+}
+
 std::ostream &operator<<(std::ostream &out, EllipticalArc const &ea)
 {
     out << "EllipticalArc("
diff --git a/src/2geom/elliptical-arc.h b/src/2geom/elliptical-arc.h
index fbb290dca..96a92eb88 100644
--- a/src/2geom/elliptical-arc.h
+++ b/src/2geom/elliptical-arc.h
@@ -215,7 +215,7 @@ public:
     /// Get the angular interval of the arc.
     AngleInterval angularInterval() const { return _angles; }
 
-    /// Evaluate the arc in the curve domain, i.e. \f$[0, 1]\$.
+    /// Evaluate the arc in the curve domain, i.e. \f$[0, 1]\f$.
     virtual Point pointAt(Coord t) const;
 
     /// Evaluate a single coordinate on the arc in the curve domain.
@@ -300,7 +300,9 @@ public:
     virtual Curve *portion(double f, double t) const;
     virtual Curve *reverse() const;
     virtual bool operator==(Curve const &c) const;
+    virtual bool isNear(Curve const &other, Coord precision) const;
     virtual void feed(PathSink &sink, bool moveto_initial) const;
+    virtual int winding(Point const &p) const;
 
 private:
     void _updateCenterAndAngles();
diff --git a/src/2geom/generic-interval.h b/src/2geom/generic-interval.h
index 716390e57..e7892e2dd 100644
--- a/src/2geom/generic-interval.h
+++ b/src/2geom/generic-interval.h
@@ -32,6 +32,7 @@
 #define LIB2GEOM_SEEN_GENERIC_INTERVAL_H
 
 #include <cassert>
+#include <iostream>
 #include <boost/none.hpp>
 #include <boost/optional.hpp>
 #include <2geom/coord.h>
@@ -291,8 +292,6 @@ public:
     /// @}
 
     /** @brief Check whether this interval is empty. */
-    bool isEmpty() { return !*this; };
-    /// Alias of isEmpty() for STL similarity.
     bool empty() { return !*this; }
 
     /** @brief Union with another interval, gracefully handling empty ones. */
@@ -340,14 +339,12 @@ inline GenericOptInterval<C> operator&(GenericInterval<C> const &a, GenericInter
     return GenericOptInterval<C>(a) & GenericOptInterval<C>(b);
 }
 
-#ifdef _GLIBCXX_IOSTREAM
 template <typename C>
 inline std::ostream &operator<< (std::ostream &os, 
                                  Geom::GenericInterval<C> const &I) {
     os << "Interval("<<I.min() << ", "<<I.max() << ")";
     return os;
 }
-#endif
 
 } // namespace Geom
 #endif // !LIB2GEOM_SEEN_GENERIC_INTERVAL_H
diff --git a/src/2geom/generic-rect.h b/src/2geom/generic-rect.h
index f7d722757..fdc7f369b 100644
--- a/src/2geom/generic-rect.h
+++ b/src/2geom/generic-rect.h
@@ -41,6 +41,7 @@
 #define LIB2GEOM_SEEN_GENERIC_RECT_H
 
 #include <limits>
+#include <iostream>
 #include <boost/optional.hpp>
 #include <2geom/coord.h>
 
@@ -393,7 +394,7 @@ public:
     /// @name Check other rectangles and points for inclusion.
     /// @{
     /** @brief Check for emptiness. */
-    inline bool isEmpty() const { return !*this; };
+    inline bool empty() const { return !*this; };
     /** @brief Check whether the rectangles have any common points.
      * Empty rectangles will not intersect with any other rectangle. */
     bool intersects(CRect const &r) const { return r.intersects(*this); }
@@ -513,13 +514,11 @@ inline bool GenericRect<C>::contains(OptCRect const &r) const {
     return !r || contains(*r);
 }
 
-#ifdef _GLIBCXX_IOSTREAM
 template <typename C>
 inline std::ostream &operator<<(std::ostream &out, GenericRect<C> const &r) {
-    out << "X: " << r[X] << "  Y: " << r[Y];
+    out << "Rect " << r[X] << " x " << r[Y];
     return out;
 }
-#endif
 
 } // end namespace Geom
 
diff --git a/src/2geom/int-point.h b/src/2geom/int-point.h
index 7c737b1d7..50d3a6720 100644
--- a/src/2geom/int-point.h
+++ b/src/2geom/int-point.h
@@ -60,15 +60,6 @@ public:
         _pt[X] = x;
         _pt[Y] = y;
     }
-    IntPoint(IntPoint const &p) {
-        _pt[X] = p._pt[X];
-        _pt[Y] = p._pt[Y];
-    }
-    IntPoint &operator=(IntPoint const &p) {
-        _pt[X] = p._pt[X];
-        _pt[Y] = p._pt[Y];
-        return *this;
-    }
     /// @}
 
     /// @name Access the coordinates of a point
@@ -92,6 +83,10 @@ public:
 
     /// @name Vector-like arithmetic operations
     /// @{
+    IntPoint operator-() const {
+        IntPoint ret(-_pt[X], -_pt[Y]);
+        return ret;
+    }
     IntPoint &operator+=(IntPoint const &o) {
         _pt[X] += o._pt[X];
         _pt[Y] += o._pt[Y];
diff --git a/src/2geom/intersection-graph.cpp b/src/2geom/intersection-graph.cpp
index e18561f67..d469d3ffc 100644
--- a/src/2geom/intersection-graph.cpp
+++ b/src/2geom/intersection-graph.cpp
@@ -34,6 +34,7 @@
 #include <2geom/intersection-graph.h>
 #include <2geom/path.h>
 #include <2geom/pathvector.h>
+#include <2geom/utils.h>
 #include <iostream>
 #include <iterator>
 
@@ -56,101 +57,174 @@ struct PathIntersectionGraph::IntersectionVertexLess {
  */
 
 PathIntersectionGraph::PathIntersectionGraph(PathVector const &a, PathVector const &b, Coord precision)
-    : _a(a)
-    , _b(b)
+    : _graph_valid(true)
 {
     if (a.empty() || b.empty()) return;
 
+    _pv[0] = a;
+    _pv[1] = b;
+
+    _prepareArguments();
+    bool has_intersections = _prepareIntersectionLists(precision);
+    if (!has_intersections) return;
+
+    _assignEdgeWindingParities(precision);
+    _assignComponentStatusFromDegenerateIntersections();
+    _removeDegenerateIntersections();
+    if (_graph_valid) {
+        _verify();
+    }
+}
+
+void PathIntersectionGraph::_prepareArguments()
+{
     // all paths must be closed, otherwise we will miss some intersections
-    for (std::size_t i = 0; i < a.size(); ++i) {
-        _a[i].close();
+    for (int w = 0; w < 2; ++w) {
+        for (std::size_t i = 0; i < _pv[w].size(); ++i) {
+            _pv[w][i].close();
+        }
     }
-    for (std::size_t i = 0; i < b.size(); ++i) {
-        _b[i].close();
+    // remove degenerate segments
+    for (int w = 0; w < 2; ++w) {
+        for (std::size_t i = _pv[w].size(); i > 0; --i) {
+            if (_pv[w][i-1].empty()) {
+                _pv[w].erase(_pv[w].begin() + (i-1));
+            }
+            for (std::size_t j = _pv[w][i-1].size(); j > 0; --j) {
+                if (_pv[w][i-1][j-1].isDegenerate()) {
+                    _pv[w][i-1].erase(_pv[w][i-1].begin() + (j-1));
+                }
+            }
+        }
     }
+}
 
-    std::vector<PVIntersection> pxs = _a.intersect(_b, precision);
+bool PathIntersectionGraph::_prepareIntersectionLists(Coord precision)
+{
+    std::vector<PVIntersection> pxs = _pv[0].intersect(_pv[1], precision);
     // NOTE: this early return means that the path data structures will not be created
     // if there are no intersections at all!
-    if (pxs.empty()) return;
+    if (pxs.empty()) return false;
 
     // prepare intersection lists for each path component
-    for (std::size_t i = 0; i < _a.size(); ++i) {
-        _apaths.push_back(new PathData());
-    }
-    for (std::size_t i = 0; i < _b.size(); ++i) {
-        _bpaths.push_back(new PathData());
+    for (unsigned w = 0; w < 2; ++w) {
+        for (std::size_t i = 0; i < _pv[w].size(); ++i) {
+            _components[w].push_back(new PathData(w, i));
+        }
     }
 
+    // create intersection vertices
     for (std::size_t i = 0; i < pxs.size(); ++i) {
         IntersectionVertex *xa, *xb;
         xa = new IntersectionVertex();
         xb = new IntersectionVertex();
-        xa->processed = xb->processed = false;
+        //xa->processed = xb->processed = false;
+        xa->which = 0; xb->which = 1;
         xa->pos = pxs[i].first;
         xb->pos = pxs[i].second;
         xa->p = xb->p = pxs[i].point();
         xa->neighbor = xb;
         xb->neighbor = xa;
+        xa->next_edge = xb->next_edge = OUTSIDE;
+        xa->defective = xb->defective = false;
         _xs.push_back(xa);
         _xs.push_back(xb);
-        _apaths[xa->pos.path_index].xlist.push_back(*xa);
-        _bpaths[xb->pos.path_index].xlist.push_back(*xb);
+        _components[0][xa->pos.path_index].xlist.push_back(*xa);
+        _components[1][xb->pos.path_index].xlist.push_back(*xb);
     }
 
-    for (std::size_t i = 0; i < _apaths.size(); ++i) {
-        _apaths[i].xlist.sort(IntersectionVertexLess());
-    }
-    for (std::size_t i = 0; i < _bpaths.size(); ++i) {
-        _bpaths[i].xlist.sort(IntersectionVertexLess());
+    // sort components according to time value of intersections
+    for (unsigned w = 0; w < 2; ++w) {
+        for (std::size_t i = 0; i < _components[w].size(); ++i) {
+            _components[w][i].xlist.sort(IntersectionVertexLess());
+        }
     }
 
-    typedef IntersectionList::iterator Iter;
-
-    // determine in/out/on flags using winding
-    for (unsigned npv = 0; npv < 2; ++npv) {
-        boost::ptr_vector<PathData> &ls = npv ? _bpaths : _apaths;
-        boost::ptr_vector<PathData> &ols = npv ? _apaths : _bpaths;
-        PathVector const &pv = npv ? b : a;
-        PathVector const &other = npv ? a : b;
-
-        for (unsigned li = 0; li < ls.size(); ++li) {
-            IntersectionList &xl = ls[li].xlist;
-            for (Iter i = xl.begin(); i != xl.end(); ++i) {
-                Iter n = boost::next(i);
-                if (n == xl.end()) {
-                    n = xl.begin();
-                }
+    return true;
+}
+
+void PathIntersectionGraph::_assignEdgeWindingParities(Coord precision)
+{
+    // determine the winding numbers of path portions between intersections
+    for (unsigned w = 0; w < 2; ++w) {
+        unsigned ow = (w+1) % 2;
+
+        for (unsigned li = 0; li < _components[w].size(); ++li) {
+            IntersectionList &xl = _components[w][li].xlist;
+            for (ILIter i = xl.begin(); i != xl.end(); ++i) {
+                ILIter n = cyclic_next(i, xl);
                 std::size_t pi = i->pos.path_index;
 
-                PathInterval ival = forward_interval(i->pos, n->pos, pv[pi].size());
+                PathInterval ival = forward_interval(i->pos, n->pos, _pv[w][pi].size());
                 PathTime mid = ival.inside(precision);
 
-                // TODO check for degenerate cases
-                int w = other.winding(pv[pi].pointAt(mid));
-                if (w % 2) {
-                    i->next = POINT_INSIDE;
-                    n->previous = POINT_INSIDE;
+                Point wpoint = _pv[w][pi].pointAt(mid);
+                _winding_points.push_back(wpoint);
+                int wdg = _pv[ow].winding(wpoint);
+                if (wdg % 2) {
+                    i->next_edge = INSIDE;
                 } else {
-                    i->next = POINT_OUTSIDE;
-                    n->previous = POINT_OUTSIDE;
+                    i->next_edge = OUTSIDE;
                 }
             }
+        }
+    }
+}
+
+void PathIntersectionGraph::_assignComponentStatusFromDegenerateIntersections()
+{
+    // If a path has only degenerate intersections, assign its status now.
+    // This protects against later accidentaly picking a point for winding
+    // determination that is exactly at a removed intersection.
+    for (unsigned w = 0; w < 2; ++w) {
+        for (unsigned li = 0; li < _components[w].size(); ++li) {
+            IntersectionList &xl = _components[w][li].xlist;
+            bool has_in = false;
+            bool has_out = false;
+            for (ILIter i = xl.begin(); i != xl.end(); ++i) {
+                has_in |= (i->next_edge == INSIDE);
+                has_out |= (i->next_edge == OUTSIDE);
+            }
+            if (has_in && !has_out) {
+                _components[w][li].status = INSIDE;
+            }
+            if (!has_in && has_out) {
+                _components[w][li].status = OUTSIDE;
+            }
+        }
+    }
+}
 
-            // remove intersections that don't change between in/out
-            // and assign exit / entry flags
-            for (Iter i = xl.begin(); i != xl.end();) {
-                if (i->previous == i->next) {
-                    IntersectionList &oxl = ols[i->neighbor->pos.path_index].xlist;
-                    oxl.erase(oxl.iterator_to(*i->neighbor));
-                    xl.erase(i++);
-                    if (i->next == POINT_INSIDE) {
-                        ++ls[li].removed_in;
-                    } else {
-                        ++ls[li].removed_out;
+void PathIntersectionGraph::_removeDegenerateIntersections()
+{
+    // remove intersections that don't change between in/out
+    for (unsigned w = 0; w < 2; ++w) {
+        for (unsigned li = 0; li < _components[w].size(); ++li) {
+            IntersectionList &xl = _components[w][li].xlist;
+            for (ILIter i = xl.begin(); i != xl.end();) {
+                ILIter n = cyclic_next(i, xl);
+                if (i->next_edge == n->next_edge) {
+                    bool last_node = (i == n);
+                    ILIter nn = _getNeighbor(n);
+                    IntersectionList &oxl = _getPathData(nn).xlist;
+
+                    // When exactly 3 out of 4 edges adjacent to an intersection
+                    // have the same winding, we have a defective intersection,
+                    // which is neither degenerate nor normal. Those can occur in paths
+                    // that contain overlapping segments. We cannot handle that case
+                    // for now, so throw an exception.
+                    if (cyclic_prior(nn, oxl)->next_edge != nn->next_edge) {
+                        _graph_valid = false;
+                        n->defective = true;
+                        nn->defective = true;
+                        ++i;
+                        continue;
                     }
+
+                    oxl.erase(nn);
+                    xl.erase(n);
+                    if (last_node) break;
                 } else {
-                    i->entry = ((i->next == POINT_INSIDE) && (i->previous == POINT_OUTSIDE));
                     ++i;
                 }
             }
@@ -158,6 +232,20 @@ PathIntersectionGraph::PathIntersectionGraph(PathVector const &a, PathVector con
     }
 }
 
+void PathIntersectionGraph::_verify()
+{
+    for (unsigned w = 0; w < 2; ++w) {
+        for (unsigned li = 0; li < _components[w].size(); ++li) {
+            IntersectionList &xl = _components[w][li].xlist;
+            assert(xl.size() % 2 == 0);
+            for (ILIter i = xl.begin(); i != xl.end(); ++i) {
+                ILIter j = cyclic_next(i, xl);
+                assert(i->next_edge != j->next_edge);
+            }
+        }
+    }
+}
+
 PathVector PathIntersectionGraph::getUnion()
 {
     PathVector result = _getResult(false, false);
@@ -192,8 +280,9 @@ PathVector PathIntersectionGraph::getBminusA()
 
 PathVector PathIntersectionGraph::getXOR()
 {
-    PathVector r1 = getAminusB();
-    PathVector r2 = getBminusA();
+    PathVector r1, r2;
+    r1 = getAminusB();
+    r2 = getBminusA();
     std::copy(r2.begin(), r2.end(), std::back_inserter(r1));
     return r1;
 }
@@ -201,37 +290,61 @@ PathVector PathIntersectionGraph::getXOR()
 std::size_t PathIntersectionGraph::size() const
 {
     std::size_t result = 0;
-    for (std::size_t i = 0; i < _apaths.size(); ++i) {
-        result += _apaths[i].xlist.size();
+    for (std::size_t i = 0; i < _components[0].size(); ++i) {
+        result += _components[0][i].xlist.size();
     }
     return result;
 }
 
-std::vector<Point> PathIntersectionGraph::intersectionPoints() const
+std::vector<Point> PathIntersectionGraph::intersectionPoints(bool defective) const
 {
     std::vector<Point> result;
 
-    typedef IntersectionList::const_iterator Iter;
-    for (std::size_t i = 0; i < _apaths.size(); ++i) {
-        for (Iter j = _apaths[i].xlist.begin(); j != _apaths[i].xlist.end(); ++j) {
-            result.push_back(j->p);
+    typedef IntersectionList::const_iterator CILIter;
+    for (std::size_t i = 0; i < _components[0].size(); ++i) {
+        for (CILIter j = _components[0][i].xlist.begin(); j != _components[0][i].xlist.end(); ++j) {
+            if (j->defective == defective) {
+                result.push_back(j->p);
+            }
         }
     }
     return result;
 }
 
+void PathIntersectionGraph::fragments(PathVector &in, PathVector &out) const
+{
+    typedef boost::ptr_vector<PathData>::const_iterator PIter;
+    for (unsigned w = 0; w < 2; ++w) {
+        for (PIter li = _components[w].begin(); li != _components[w].end(); ++li) {
+            for (CILIter k = li->xlist.begin(); k != li->xlist.end(); ++k) {
+                CILIter n = cyclic_next(k, li->xlist);
+                // TODO: investigate why non-contiguous paths are sometimes generated here
+                Path frag(k->p);
+                frag.setStitching(true);
+                PathInterval ival = forward_interval(k->pos, n->pos, _pv[w][k->pos.path_index].size());
+                _pv[w][k->pos.path_index].appendPortionTo(frag, ival, k->p, n->p);
+                if (k->next_edge == INSIDE) {
+                    in.push_back(frag);
+                } else {
+                    out.push_back(frag);
+                }
+            }
+        }
+    }
+}
+
 PathVector PathIntersectionGraph::_getResult(bool enter_a, bool enter_b)
 {
-    typedef IntersectionList::iterator Iter;
+    typedef boost::ptr_vector<PathData>::iterator PIter;
     PathVector result;
     if (_xs.empty()) return result;
 
     // reset processed status
-    for (unsigned npv = 0; npv < 2; ++npv) {
-        boost::ptr_vector<PathData> &ls = npv ? _bpaths : _apaths;
-        for (std::size_t li = 0; li < ls.size(); ++li) {
-            for (Iter k = ls[li].xlist.begin(); k != ls[li].xlist.end(); ++k) {
-                k->processed = false;
+    _ulist.clear();
+    for (unsigned w = 0; w < 2; ++w) {
+        for (PIter li = _components[w].begin(); li != _components[w].end(); ++li) {
+            for (ILIter k = li->xlist.begin(); k != li->xlist.end(); ++k) {
+                _ulist.push_back(*k);
             }
         }
     }
@@ -239,18 +352,17 @@ PathVector PathIntersectionGraph::_getResult(bool enter_a, bool enter_b)
     unsigned n_processed = 0;
 
     while (true) {
-        PathVector const *cur = &_a, *other = &_b;
-        boost::ptr_vector<PathData> *lscur = &_apaths, *lsother = &_bpaths;
-
-        // find unprocessed intersection
-        Iter i;
-        if (!_findUnprocessed(i)) break;
+        // get unprocessed intersection
+        if (_ulist.empty()) break;
+        IntersectionVertex &iv = _ulist.front();
+        unsigned w = iv.which;
+        ILIter i = _components[w][iv.pos.path_index].xlist.iterator_to(iv);
 
         result.push_back(Path(i->p));
         result.back().setStitching(true);
 
-        while (!i->processed) {
-            Iter prev = i;
+        while (i->_proc_hook.is_linked()) {
+            ILIter prev = i;
             std::size_t pi = i->pos.path_index;
             // determine which direction to go
             // union: always go outside
@@ -258,53 +370,55 @@ PathVector PathIntersectionGraph::_getResult(bool enter_a, bool enter_b)
             // a minus b: go inside in b, outside in a
             // b minus a: go inside in a, outside in b
             bool reverse = false;
-            if (cur == &_a) {
-                reverse = i->entry ^ enter_a;
+            if (w == 0) {
+                reverse = (i->next_edge == INSIDE) ^ enter_a;
             } else {
-                reverse = i->entry ^ enter_b;
+                reverse = (i->next_edge == INSIDE) ^ enter_b;
             }
 
             // get next intersection
             if (reverse) {
-                if (i == (*lscur)[pi].xlist.begin()) {
-                    i = (*lscur)[pi].xlist.end();
-                }
-                --i;
+                i = cyclic_prior(i, _components[w][pi].xlist);
             } else {
-                ++i;
-                if (i == (*lscur)[pi].xlist.end()) {
-                    i = (*lscur)[pi].xlist.begin();
-                }
+                i = cyclic_next(i, _components[w][pi].xlist);
             }
 
             // append portion of path
             PathInterval ival = PathInterval::from_direction(
                 prev->pos.asPathTime(), i->pos.asPathTime(),
-                reverse, (*cur)[pi].size());
+                reverse, _pv[i->which][pi].size());
 
-            (*cur)[pi].appendPortionTo(result.back(), ival, prev->p, i->p);
+            _pv[i->which][pi].appendPortionTo(result.back(), ival, prev->p, i->p);
 
             // mark both vertices as processed
-            prev->processed = true;
-            i->processed = true;
+            //prev->processed = true;
+            //i->processed = true;
             n_processed += 2;
+            if (prev->_proc_hook.is_linked()) {
+                _ulist.erase(_ulist.iterator_to(*prev));
+            }
+            if (i->_proc_hook.is_linked()) {
+                _ulist.erase(_ulist.iterator_to(*i));
+            }
 
             // switch to the other path
-            i = (*lsother)[i->neighbor->pos.path_index].xlist.iterator_to(*i->neighbor);
-            std::swap(lscur, lsother);
-            std::swap(cur, other);
+            i = _getNeighbor(i);
+            w = i->which;
         }
         result.back().close(true);
 
         assert(!result.back().empty());
     }
 
+    /*if (n_processed != size() * 2) {
+        std::cerr << "Processed " << n_processed << " intersections, expected " << (size() * 2) << std::endl;
+    }*/
     assert(n_processed == size() * 2);
 
     return result;
 }
 
-void PathIntersectionGraph::_handleNonintersectingPaths(PathVector &result, int which, bool inside)
+void PathIntersectionGraph::_handleNonintersectingPaths(PathVector &result, unsigned which, bool inside)
 {
     /* Every component that has any intersections will be processed by _getResult.
      * Here we take care of paths that don't have any intersections. They are either
@@ -312,68 +426,53 @@ void PathIntersectionGraph::_handleNonintersectingPaths(PathVector &result, int
      * evaluating the winding rule at the initial point. If inside is true and
      * the path is inside, we add it to the result.
      */
-    boost::ptr_vector<PathData> const &ls = which ? _bpaths : _apaths;
-    PathVector const &cur = which ? _b : _a;
-    PathVector const &other = which ? _a : _b;
+    unsigned w = which;
+    unsigned ow = (w+1) % 2;
 
-    for (std::size_t i = 0; i < cur.size(); ++i) {
+    for (std::size_t i = 0; i < _pv[w].size(); ++i) {
         // the path data vector might have been left empty if there were no intersections at all
-        bool has_path_data = !ls.empty();
+        bool has_path_data = !_components[w].empty();
         // Skip if the path has intersections
-        if (has_path_data && !ls[i].xlist.empty()) continue;
+        if (has_path_data && !_components[w][i].xlist.empty()) continue;
         bool path_inside = false;
 
-        // If the path had any intersections removed, use the result of that,
-        // since one of those might have been at the initial point.
-        // Also, it saves time.
-        if (has_path_data && ls[i].removed_in != 0) {
+        // Use the in/out determination from constructor, if available
+        if (has_path_data && _components[w][i].status == INSIDE) {
             path_inside = true;
-        } else if (has_path_data && ls[i].removed_out != 0) {
+        } else if (has_path_data && _components[w][i].status == OUTSIDE) {
             path_inside = false;
         } else {
-            int w = other.winding(cur[i].initialPoint());
-            path_inside = w % 2 != 0;
+            int wdg = _pv[ow].winding(_pv[w][i].initialPoint());
+            path_inside = wdg % 2 != 0;
         }
 
         if (path_inside == inside) {
-            result.push_back(cur[i]);
+            result.push_back(_pv[w][i]);
         }
     }
 }
 
-bool PathIntersectionGraph::_findUnprocessed(IntersectionList::iterator &result)
+PathIntersectionGraph::ILIter PathIntersectionGraph::_getNeighbor(ILIter iter)
 {
-    typedef IntersectionList::iterator Iter;
-
-    Iter it, last;
-
-    for (std::size_t k = 0; k < _apaths.size(); ++k) {
-        it = _apaths[k].xlist.begin();
-        last = _apaths[k].xlist.end();
-        for (; it != last; ++it) {
-            if (!it->_hook.is_linked()) {
-                // this intersection was removed since it did not change inside/outside status
-                continue;
-            }
-            if (!it->processed) {
-                result = it;
-                return true;
-            }
-        }
-    }
-
-    return false;
+    unsigned ow = (iter->which + 1) % 2;
+    return _components[ow][iter->neighbor->pos.path_index].xlist.iterator_to(*iter->neighbor);
 }
 
+PathIntersectionGraph::PathData &
+PathIntersectionGraph::_getPathData(ILIter iter)
+{
+    return _components[iter->which][iter->pos.path_index];
+}
 
 std::ostream &operator<<(std::ostream &os, PathIntersectionGraph const &pig)
 {
+    typedef PathIntersectionGraph::IntersectionList::const_iterator CILIter;
     os << "Intersection graph:\n"
        << pig._xs.size()/2 << " total intersections\n"
        << pig.size() << " considered intersections\n";
-    for (std::size_t i = 0; i < pig._apaths.size(); ++i) {
-        typedef PathIntersectionGraph::IntersectionList::const_iterator Iter;
-        for (Iter j = pig._apaths[i].xlist.begin(); j != pig._apaths[i].xlist.end(); ++j) {
+    for (std::size_t i = 0; i < pig._components[0].size(); ++i) {
+        PathIntersectionGraph::IntersectionList const &xl = pig._components[0][i].xlist;
+        for (CILIter j = xl.begin(); j != xl.end(); ++j) {
             os << j->pos << " - " << j->neighbor->pos << " @ " << j->p << "\n";
         }
     }
diff --git a/src/2geom/intersection-graph.h b/src/2geom/intersection-graph.h
index bd9aaee81..332c83f18 100644
--- a/src/2geom/intersection-graph.h
+++ b/src/2geom/intersection-graph.h
@@ -58,23 +58,29 @@ public:
 
     /// Returns the number of intersections used when computing Boolean operations.
     std::size_t size() const;
-    std::vector<Point> intersectionPoints() const;
+    std::vector<Point> intersectionPoints(bool defective = false) const;
+    std::vector<Point> windingPoints() const {
+        return _winding_points;
+    }
+    void fragments(PathVector &in, PathVector &out) const;
+    bool valid() const { return _graph_valid; }
 
 private:
     enum InOutFlag {
-        POINT_INSIDE,
-        POINT_OUTSIDE
+        INSIDE,
+        OUTSIDE,
+        BOTH
     };
 
     struct IntersectionVertex {
         boost::intrusive::list_member_hook<> _hook;
+        boost::intrusive::list_member_hook<> _proc_hook;
         PathVectorTime pos;
         Point p; // guarantees that endpoints are exact
         IntersectionVertex *neighbor;
-        bool entry; // going in +t direction enters the other path
-        InOutFlag previous;
-        InOutFlag next;
-        bool processed; // TODO: use intrusive unprocessed list instead
+        InOutFlag next_edge;
+        unsigned which;
+        bool defective;
     };
 
     typedef boost::intrusive::list
@@ -86,27 +92,50 @@ private:
             >
         > IntersectionList;
 
+    typedef boost::intrusive::list
+        < IntersectionVertex
+        , boost::intrusive::member_hook
+            < IntersectionVertex
+            , boost::intrusive::list_member_hook<>
+            , &IntersectionVertex::_proc_hook
+            >
+        > UnprocessedList;
+
     struct PathData {
         IntersectionList xlist;
-        unsigned removed_in;
-        unsigned removed_out;
-
-        PathData()
-            : removed_in(0)
-            , removed_out(0)
+        std::size_t path_index;
+        int which;
+        InOutFlag status;
+
+        PathData(int w, std::size_t pi)
+            : path_index(pi)
+            , which(w)
+            , status(BOTH)
         {}
     };
 
     struct IntersectionVertexLess;
+    typedef IntersectionList::iterator ILIter;
+    typedef IntersectionList::const_iterator CILIter;
 
     PathVector _getResult(bool enter_a, bool enter_b);
-    void _handleNonintersectingPaths(PathVector &result, int which, bool inside);
-    bool _findUnprocessed(IntersectionList::iterator &result);
-
-    PathVector _a, _b;
+    void _handleNonintersectingPaths(PathVector &result, unsigned which, bool inside);
+    void _prepareArguments();
+    bool _prepareIntersectionLists(Coord precision);
+    void _assignEdgeWindingParities(Coord precision);
+    void _assignComponentStatusFromDegenerateIntersections();
+    void _removeDegenerateIntersections();
+    void _verify();
+
+    ILIter _getNeighbor(ILIter iter);
+    PathData &_getPathData(ILIter iter);
+
+    PathVector _pv[2];
     boost::ptr_vector<IntersectionVertex> _xs;
-    boost::ptr_vector<PathData> _apaths;
-    boost::ptr_vector<PathData> _bpaths;
+    boost::ptr_vector<PathData> _components[2];
+    UnprocessedList _ulist;
+    bool _graph_valid;
+    std::vector<Point> _winding_points;
 
     friend std::ostream &operator<<(std::ostream &, PathIntersectionGraph const &);
 };
diff --git a/src/2geom/line.cpp b/src/2geom/line.cpp
index bada8ef38..a9cc0e251 100644
--- a/src/2geom/line.cpp
+++ b/src/2geom/line.cpp
@@ -103,7 +103,7 @@ void Line::setCoefficients (Coord a, Coord b, Coord c)
 
 void Line::coefficients(Coord &a, Coord &b, Coord &c) const
 {
-    Point v = versor().cw();
+    Point v = vector().cw();
     a = v[X];
     b = v[Y];
     c = cross(_initial, _final);
@@ -139,7 +139,7 @@ std::vector<Coord> Line::roots(Coord v, Dim2 d) const {
 Coord Line::root(Coord v, Dim2 d) const
 {
     assert(d == X || d == Y);
-    Point vs = versor();
+    Point vs = vector();
     if (vs[d] != 0) {
         return (v - _initial[d]) / vs[d];
     } else {
@@ -149,7 +149,7 @@ Coord Line::root(Coord v, Dim2 d) const
 
 boost::optional<LineSegment> Line::clip(Rect const &r) const
 {
-    Point v = versor();
+    Point v = vector();
     // handle horizontal and vertical lines first,
     // since the root-based code below will break for them
     for (unsigned i = 0; i < 2; ++i) {
@@ -221,7 +221,7 @@ boost::optional<LineSegment> Line::clip(Rect const &r) const
  * @see timeAtProjection */
 Coord Line::timeAt(Point const &p) const
 {
-    Point v = versor();
+    Point v = vector();
     // degenerate case
     if (v[X] == 0 && v[Y] == 0) {
         return 0;
@@ -235,12 +235,26 @@ Coord Line::timeAt(Point const &p) const
     }
 }
 
+/** @brief Create a transformation that maps one line to another.
+ * This will return a transformation \f$A\f$ such that
+ * \f$L_1(t) \cdot A = L_2(t)\f$, where \f$L_1\f$ is this line
+ * and \f$L_2\f$ is the line passed as the parameter. The returned
+ * transformation will preserve angles. */
+Affine Line::transformTo(Line const &other) const
+{
+    Affine result = Translate(-_initial);
+    result *= Rotate(angle_between(vector(), other.vector()));
+    result *= Scale(other.vector().length() / vector().length());
+    result *= Translate(other._initial);
+    return result;
+}
+
 std::vector<ShapeIntersection> Line::intersect(Line const &other) const
 {
     std::vector<ShapeIntersection> result;
 
-    Point v1 = versor();
-    Point v2 = other.versor();
+    Point v1 = vector();
+    Point v2 = other.vector();
     Coord cp = cross(v1, v2);
     if (cp == 0) return result;
 
@@ -319,8 +333,8 @@ OptCrossing intersection_impl(Ray const& r1, Line const& l2, unsigned int i)
     using std::swap;
 
     OptCrossing crossing =
-        intersection_impl(r1.versor(), r1.origin(),
-                          l2.versor(), l2.origin() );
+        intersection_impl(r1.vector(), r1.origin(),
+                          l2.vector(), l2.origin() );
 
     if (crossing) {
         if (crossing->ta < 0) {
@@ -349,7 +363,7 @@ OptCrossing intersection_impl( LineSegment const& ls1,
     OptCrossing crossing =
         intersection_impl(ls1.finalPoint() - ls1.initialPoint(),
                           ls1.initialPoint(),
-                          l2.versor(),
+                          l2.vector(),
                           l2.origin() );
 
     if (crossing) {
@@ -382,7 +396,7 @@ OptCrossing intersection_impl( LineSegment const& ls1,
     OptCrossing crossing =
         intersection_impl( direction,
                            ls1.initialPoint(),
-                           r2.versor(),
+                           r2.vector(),
                            r2.origin() );
 
     if (crossing) {
@@ -400,7 +414,7 @@ OptCrossing intersection_impl( LineSegment const& ls1,
     }
 
     if ( are_near(r2.origin(), ls1) ) {
-        bool eqvs = (dot(direction, r2.versor()) > 0);
+        bool eqvs = (dot(direction, r2.vector()) > 0);
         if ( are_near(ls1.initialPoint(), r2.origin()) && !eqvs)  {
             crossing->ta = crossing->tb = 0;
             return crossing;
@@ -431,8 +445,8 @@ OptCrossing intersection_impl( LineSegment const& ls1,
 OptCrossing intersection(Line const& l1, Line const& l2)
 {
     OptCrossing c = detail::intersection_impl(
-        l1.versor(), l1.origin(),
-        l2.versor(), l2.origin());
+        l1.vector(), l1.origin(),
+        l2.vector(), l2.origin());
 
     if (!c && distance(l1.origin(), l2) == 0) {
         THROW_INFINITESOLUTIONS();
@@ -443,8 +457,8 @@ OptCrossing intersection(Line const& l1, Line const& l2)
 OptCrossing intersection(Ray const& r1, Ray const& r2)
 {
     OptCrossing crossing =
-    detail::intersection_impl( r1.versor(), r1.origin(),
-                               r2.versor(), r2.origin() );
+    detail::intersection_impl( r1.vector(), r1.origin(),
+                               r2.vector(), r2.origin() );
 
     if (crossing)
     {
@@ -463,7 +477,7 @@ OptCrossing intersection(Ray const& r1, Ray const& r2)
     if ( are_near(r1.origin(), r2) || are_near(r2.origin(), r1) )
     {
         if ( are_near(r1.origin(), r2.origin())
-             && !are_near(r1.versor(), r2.versor()) )
+             && !are_near(r1.vector(), r2.vector()) )
         {
             crossing->ta = crossing->tb = 0;
             return crossing;
@@ -568,7 +582,7 @@ Line make_angle_bisector_line(Line const& l1, Line const& l2)
     }
     Point O = l1.pointAt(crossing->ta);
     Point A = l1.pointAt(crossing->ta + 1);
-    double angle = angle_between(l1.versor(), l2.versor());
+    double angle = angle_between(l1.vector(), l2.vector());
     Point B = (angle > 0) ? l2.pointAt(crossing->tb + 1)
         : l2.pointAt(crossing->tb - 1);
 
diff --git a/src/2geom/line.h b/src/2geom/line.h
index 7d4766e12..62b3652f1 100644
--- a/src/2geom/line.h
+++ b/src/2geom/line.h
@@ -98,7 +98,7 @@ public:
     /// Create a line by extending a ray.
     explicit Line(Ray const &r)
         : _initial(r.origin())
-        , _final(r.origin() + r.versor())
+        , _final(r.origin() + r.vector())
     {}
 
     /// Create a line normal to a vector at a specified distance from origin.
@@ -111,7 +111,7 @@ public:
      * Note that each line direction has two possible unit vectors.
      * @param o Point through which the line will pass
      * @param v Unit vector of the line's direction */
-    static Line from_origin_and_versor(Point const &o, Point const &v) {
+    static Line from_origin_and_vector(Point const &o, Point const &v) {
         Line l(o, o + v);
         return l;
     }
@@ -126,9 +126,12 @@ public:
 
     /// Get the line's origin point.
     Point origin() const { return _initial; }
-    /** @brief Get the line's direction vector.
-     * Note that the retrieved vector is not normalized to unit length. */
-    Point versor() const { return _final - _initial; }
+    /** @brief Get the line's raw direction vector.
+     * The retrieved vector is normalized to unit length. */
+    Point vector() const { return _final - _initial; }
+    /** @brief Get the line's normalized direction vector.
+     * The retrieved vector is normalized to unit length. */
+    Point versor() const { return (_final - _initial).normalized(); }
     /// Angle the line makes with the X axis, in mathematical convention.
     Coord angle() const {
         Point d = _final - _initial;
@@ -147,7 +150,7 @@ public:
     }
     /** @brief Set the speed of the line.
      * Origin remains unchanged. */
-    void setVersor(Point const &v) {
+    void setVector(Point const &v) {
         _final = _initial + v;
     }
 
@@ -239,7 +242,7 @@ public:
      * @return Time value corresponding to a point closest to @c p. */
     Coord timeAtProjection(Point const& p) const {
         if ( isDegenerate() ) return 0;
-        Point v = versor();
+        Point v = vector();
         return dot(p - _initial, v) / dot(v, v);
     }
 
@@ -284,22 +287,22 @@ public:
     boost::optional<LineSegment> clip(Rect const &r) const;
 
     /** @brief Create a ray starting at the specified time value.
-     * The created ray will go in the direction of the line's versor (in the direction
+     * The created ray will go in the direction of the line's vector (in the direction
      * of increasing time values).
      * @param t Time value where the ray should start
-     * @return Ray starting at t and going in the direction of the versor */
+     * @return Ray starting at t and going in the direction of the vector */
     Ray ray(Coord t) {
         Ray result;
         result.setOrigin(pointAt(t));
-        result.setVersor(versor());
+        result.setVector(vector());
         return result;
     }
 
     /** @brief Create a derivative of the line.
      * The new line will always be degenerate. Its origin will be equal to this
-     * line's versor. */
+     * line's vector. */
     Line derivative() const {
-        Point v = versor();
+        Point v = vector();
         Line result(v, v);
         return result;
     }
@@ -314,7 +317,7 @@ public:
      * If Y grows upwards, then this is the left normal. If Y grows downwards,
      * then this is the right normal. */
     Point normal() const {
-        return rot90(versor()).normalized();
+        return rot90(vector()).normalized();
     }
 
     // what does this do?
@@ -326,7 +329,7 @@ public:
 
     /// Compute an affine matrix representing a reflection about the line.
     Affine reflection() const {
-        Point v = versor().normalized();
+        Point v = versor();
         Coord x2 = v[X]*v[X], y2 = v[Y]*v[Y], xy = v[X]*v[Y];
         Affine m(x2-y2, 2.*xy,
                  2.*xy, y2-x2,
@@ -344,7 +347,7 @@ public:
      * lower precision than e.g. a shear transform.
      * @param d Which coordinate of points on the line should be zero after the transformation */
     Affine rotationToZero(Dim2 d) const {
-        Point v = versor();
+        Point v = vector();
         if (d == X) {
             std::swap(v[X], v[Y]);
         } else {
@@ -359,6 +362,8 @@ public:
         Affine m = rotationToZero(other_dimension(d));
         return m;
     }
+
+    Affine transformTo(Line const &other) const;
     /// @}
 
     std::vector<ShapeIntersection> intersect(Line const &other) const;
@@ -419,7 +424,7 @@ bool are_near(Point const &p, Line const &line, double eps = EPSILON)
 inline
 bool are_parallel(Line const &l1, Line const &l2, double eps = EPSILON)
 {
-    return are_near(cross(l1.versor(), l2.versor()), 0, eps);
+    return are_near(cross(l1.vector(), l2.vector()), 0, eps);
 }
 
 /** @brief Test whether two lines are approximately the same.
@@ -440,7 +445,7 @@ bool are_same(Line const &l1, Line const &l2, double eps = EPSILON)
 inline
 bool are_orthogonal(Line const &l1, Line const &l2, double eps = EPSILON)
 {
-    return are_near(dot(l1.versor(), l2.versor()), 0, eps);
+    return are_near(dot(l1.vector(), l2.vector()), 0, eps);
 }
 
 // evaluate the angle between l1 and l2 rotating l1 in cw direction
@@ -449,7 +454,7 @@ bool are_orthogonal(Line const &l1, Line const &l2, double eps = EPSILON)
 inline
 double angle_between(Line const& l1, Line const& l2)
 {
-    double angle = angle_between(l1.versor(), l2.versor());
+    double angle = angle_between(l1.vector(), l2.vector());
     if (angle < 0) angle += M_PI;
     if (angle == M_PI) angle = 0;
     return angle;
@@ -472,7 +477,7 @@ bool are_near(Point const &p, LineSegment const &seg, double eps = EPSILON)
 inline
 Line make_orthogonal_line(Point const &p, Line const &line)
 {
-    Point d = line.versor().cw();
+    Point d = line.vector().cw();
     Line l(p, p + d);
     return l;
 }
diff --git a/src/2geom/linear.h b/src/2geom/linear.h
index b0306fb9f..316ed8639 100644
--- a/src/2geom/linear.h
+++ b/src/2geom/linear.h
@@ -5,8 +5,9 @@
  *  Authors:
  *   Nathan Hurst <njh@mail.csse.monash.edu.au>
  *   Michael Sloan <mgsloan@gmail.com>
+ *   Krzysztof Kosiński <tweenk.pl@gmail.com>
  *
- * Copyright (C) 2006-2007 authors
+ * Copyright (C) 2006-2015 authors
  *
  * This library is free software; you can redistribute it and/or
  * modify it either under the terms of the GNU Lesser General Public
@@ -38,14 +39,6 @@
 #include <2geom/interval.h>
 #include <2geom/math-utils.h>
 
-//#define USE_SBASIS_OF
-
-#ifdef USE_SBASIS_OF
-
-#include "linear-of.h"
-
-#else
-
 namespace Geom {
 
 class SBasis;
@@ -54,26 +47,31 @@ class SBasis;
  * @brief Function that interpolates linearly between two values.
  * @ingroup Fragments
  */
-class Linear {
+class Linear
+    : boost::additive< Linear
+    , boost::arithmetic< Linear, Coord
+    , boost::equality_comparable< Linear
+      > > >
+{
 public:
-    double a[2];
+    Coord a[2];
     Linear() {a[0]=0; a[1]=0;}
-    Linear(double aa, double b) {a[0] = aa; a[1] = b;}
-    Linear(double aa) {a[0] = aa; a[1] = aa;}
+    Linear(Coord aa, Coord b) {a[0] = aa; a[1] = b;}
+    Linear(Coord aa) {a[0] = aa; a[1] = aa;}
 
-    double operator[](unsigned i) const {
+    Coord operator[](unsigned i) const {
         assert(i < 2);
         return a[i];
     }
-    double &operator[](unsigned i) {
+    Coord &operator[](unsigned i) {
         assert(i < 2);
         return a[i];
     }
 
     //IMPL: FragmentConcept
-    typedef double output_type;
-    bool isZero(double eps=EPSILON) const { return are_near(a[0], 0., eps) && are_near(a[1], 0., eps); }
-    bool isConstant(double eps=EPSILON) const { return are_near(a[0], a[1], eps); }
+    typedef Coord output_type;
+    bool isZero(Coord eps=EPSILON) const { return are_near(a[0], 0., eps) && are_near(a[1], 0., eps); }
+    bool isConstant(Coord eps=EPSILON) const { return are_near(a[0], a[1], eps); }
     bool isFinite() const { return IS_FINITE(a[0]) && IS_FINITE(a[1]); }
 
     Coord at0() const { return a[0]; }
@@ -81,10 +79,10 @@ public:
     Coord at1() const { return a[1]; }
     Coord &at1() { return a[1]; }
 
-    double valueAt(double t) const { return lerp(t, a[0], a[1]); }
-    double operator()(double t) const { return valueAt(t); }
+    Coord valueAt(Coord t) const { return lerp(t, a[0], a[1]); }
+    Coord operator()(Coord t) const { return valueAt(t); }
 
-    // not very useful, but required for ShapeConcept
+    // not very useful, but required for FragmentConcept
     std::vector<Coord> valueAndDerivatives(Coord t, unsigned n) {
         std::vector<Coord> result(n+1, 0.0);
         result[0] = valueAt(t);
@@ -107,6 +105,44 @@ public:
     double hat() const {
         return (a[1] + a[0])/2;
     }
+
+    // addition of other Linears
+    Linear &operator+=(Linear const &other) {
+        a[0] += other.a[0];
+        a[1] += other.a[1];
+        return *this;
+    }
+    Linear &operator-=(Linear const &other) {
+        a[0] -= other.a[0];
+        a[1] -= other.a[1];
+        return *this;
+    }
+
+    // 
+    Linear &operator+=(Coord x) {
+        a[0] += x; a[1] += x;
+        return *this;
+    }
+    Linear &operator-=(Coord x) {
+        a[0] -= x; a[1] -= x;
+        return *this;
+    }
+    Linear &operator*=(Coord x) {
+        a[0] *= x; a[1] *= x;
+        return *this;
+    }
+    Linear &operator/=(Coord x) {
+        a[0] /= x; a[1] /= x;
+        return *this;
+    }
+    Linear operator-() const {
+        Linear ret(-a[0], -a[1]);
+        return ret;
+    }
+
+    bool operator==(Linear const &other) const {
+        return a[0] == other.a[0] && a[1] == other.a[1];
+    }
 };
 
 inline Linear reverse(Linear const &a) { return Linear(a[1], a[0]); }
@@ -115,64 +151,7 @@ inline Linear portion(Linear const &a, Coord from, Coord to) {
     return result;
 }
 
-//IMPL: AddableConcept
-inline Linear operator+(Linear const & a, Linear const & b) {
-    return Linear(a[0] + b[0], a[1] + b[1]);
-}
-inline Linear operator-(Linear const & a, Linear const & b) {
-    return Linear(a[0] - b[0], a[1] - b[1]);
-}
-inline Linear& operator+=(Linear & a, Linear const & b) {
-    a[0] += b[0]; a[1] += b[1];
-    return a;
-}
-inline Linear& operator-=(Linear & a, Linear const & b) {
-    a[0] -= b[0]; a[1] -= b[1];
-    return a;
-}
-//IMPL: OffsetableConcept
-inline Linear operator+(Linear const & a, double b) {
-    return Linear(a[0] + b, a[1] + b);
-}
-inline Linear operator-(Linear const & a, double b) {
-    return Linear(a[0] - b, a[1] - b);
-}
-inline Linear& operator+=(Linear & a, double b) {
-    a[0] += b; a[1] += b;
-    return a;
-}
-inline Linear& operator-=(Linear & a, double b) {
-    a[0] -= b; a[1] -= b;
-    return a;
-}
-//IMPL: boost::EqualityComparableConcept
-inline bool operator==(Linear const & a, Linear const & b) {
-    return a[0] == b[0] && a[1] == b[1];
-}
-inline bool operator!=(Linear const & a, Linear const & b) {
-    return a[0] != b[0] || a[1] != b[1];
-}
-//IMPL: ScalableConcept
-inline Linear operator-(Linear const &a) {
-    return Linear(-a[0], -a[1]);
-}
-inline Linear operator*(Linear const & a, double b) {
-    return Linear(a[0]*b, a[1]*b);
-}
-inline Linear operator/(Linear const & a, double b) {
-    return Linear(a[0]/b, a[1]/b);
-}
-inline Linear operator*=(Linear & a, double b) {
-    a[0] *= b; a[1] *= b;
-    return a;
-}
-inline Linear operator/=(Linear & a, double b) {
-    a[0] /= b; a[1] /= b;
-    return a;
-}
-
-}
-#endif
+} // end namespace Geom
 
 #endif //LIB2GEOM_SEEN_LINEAR_H
 
diff --git a/src/2geom/numeric/fitting-model.h b/src/2geom/numeric/fitting-model.h
index fb96d1d2a..b2ca92ad8 100644
--- a/src/2geom/numeric/fitting-model.h
+++ b/src/2geom/numeric/fitting-model.h
@@ -372,7 +372,6 @@ class LFMSBasis
 
     void instance(SBasis & sb, ConstVectorView const& raw_data) const
     {
-        sb.clear();
         sb.resize(m_order+1);
         for (unsigned int i = 0, k = 0; i < raw_data.size(); i+=2, ++k)
         {
diff --git a/src/2geom/path-sink.cpp b/src/2geom/path-sink.cpp
index 3b8d407f8..77301b716 100644
--- a/src/2geom/path-sink.cpp
+++ b/src/2geom/path-sink.cpp
@@ -73,8 +73,8 @@ void PathSink::feed(Rect const &r) {
 void PathSink::feed(Circle const &e) {
     Coord r = e.radius();
     Point c = e.center();
-    Point a = c + Point(0, c[Y] + r);
-    Point b = c + Point(0, c[Y] - r);
+    Point a = c + Point(0, +r);
+    Point b = c + Point(0, -r);
 
     moveTo(a);
     arcTo(r, r, 0, false, false, b);
diff --git a/src/2geom/path-sink.h b/src/2geom/path-sink.h
index 17ede18a4..3bdb00783 100644
--- a/src/2geom/path-sink.h
+++ b/src/2geom/path-sink.h
@@ -179,8 +179,10 @@ public:
     }
 
     void closePath() {
-        _path.close();
-        flush();
+        if (_in_path) {
+            _path.close();
+            flush();
+        }
     }
 
     void flush() {
diff --git a/src/2geom/path.cpp b/src/2geom/path.cpp
index 871441751..2e2bce9f1 100644
--- a/src/2geom/path.cpp
+++ b/src/2geom/path.cpp
@@ -98,6 +98,24 @@ bool PathInterval::contains(PathTime const &pos) const {
     }
 }
 
+PathInterval::size_type PathInterval::curveCount() const
+{
+    if (isDegenerate()) return 0;
+    if (_cross_start) {
+        if (_reverse) {
+            return _path_size - _to.curve_index + _from.curve_index + 1;
+        } else {
+            return _path_size - _from.curve_index + _to.curve_index + 1;
+        }
+    } else {
+        if (_reverse) {
+            return _from.curve_index - _to.curve_index + 1;
+        } else {
+            return _to.curve_index - _from.curve_index + 1;
+        }
+    }
+}
+
 PathTime PathInterval::inside(Coord min_dist) const
 {
     // If there is some node further than min_dist (in time coord) from the ends,
@@ -220,25 +238,25 @@ PathInterval PathInterval::from_direction(PathTime const &from, PathTime const &
 
 
 Path::Path(Rect const &r)
-    : _curves(new Sequence())
+    : _data(new PathData())
     , _closing_seg(new ClosingSegment(r.corner(3), r.corner(0)))
     , _closed(true)
     , _exception_on_stitch(true)
 {
     for (unsigned i = 0; i < 3; ++i) {
-        _curves->push_back(new LineSegment(r.corner(i), r.corner(i+1)));
+        _data->curves.push_back(new LineSegment(r.corner(i), r.corner(i+1)));
     }
-    _curves->push_back(_closing_seg);
+    _data->curves.push_back(_closing_seg);
 }
 
 Path::Path(ConvexHull const &ch)
-    : _curves(new Sequence())
+    : _data(new PathData())
     , _closing_seg(new ClosingSegment(Point(), Point()))
     , _closed(true)
     , _exception_on_stitch(true)
 {
     if (ch.empty()) {
-        _curves->push_back(_closing_seg);
+        _data->curves.push_back(_closing_seg);
         return;
     }
 
@@ -248,52 +266,52 @@ Path::Path(ConvexHull const &ch)
     Point last = ch.front();
 
     for (std::size_t i = 1; i < ch.size(); ++i) {
-        _curves->push_back(new LineSegment(last, ch[i]));
+        _data->curves.push_back(new LineSegment(last, ch[i]));
         last = ch[i];
     }
 
-    _curves->push_back(_closing_seg);
+    _data->curves.push_back(_closing_seg);
     _closed = true;
 }
 
 Path::Path(Circle const &c)
-    : _curves(new Sequence())
+    : _data(new PathData())
     , _closing_seg(NULL)
     , _closed(true)
     , _exception_on_stitch(true)
 {
     Point p1 = c.pointAt(0);
     Point p2 = c.pointAt(M_PI);
-    _curves->push_back(new EllipticalArc(p1, c.radius(), c.radius(), 0, false, true, p2));
-    _curves->push_back(new EllipticalArc(p2, c.radius(), c.radius(), 0, false, true, p1));
+    _data->curves.push_back(new EllipticalArc(p1, c.radius(), c.radius(), 0, false, true, p2));
+    _data->curves.push_back(new EllipticalArc(p2, c.radius(), c.radius(), 0, false, true, p1));
     _closing_seg = new ClosingSegment(p1, p1);
-    _curves->push_back(_closing_seg);
+    _data->curves.push_back(_closing_seg);
 }
 
 Path::Path(Ellipse const &e)
-    : _curves(new Sequence())
+    : _data(new PathData())
     , _closing_seg(NULL)
     , _closed(true)
     , _exception_on_stitch(true)
 {
     Point p1 = e.pointAt(0);
     Point p2 = e.pointAt(M_PI);
-    _curves->push_back(new EllipticalArc(p1, e.rays(), e.rotationAngle(), false, true, p2));
-    _curves->push_back(new EllipticalArc(p2, e.rays(), e.rotationAngle(), false, true, p1));
+    _data->curves.push_back(new EllipticalArc(p1, e.rays(), e.rotationAngle(), false, true, p2));
+    _data->curves.push_back(new EllipticalArc(p2, e.rays(), e.rotationAngle(), false, true, p1));
     _closing_seg = new ClosingSegment(p1, p1);
-    _curves->push_back(_closing_seg);
+    _data->curves.push_back(_closing_seg);
 }
 
 void Path::close(bool c)
 {
     if (c == _closed) return;
-    if (c && _curves->size() >= 2) {
+    if (c && _data->curves.size() >= 2) {
         // when closing, if last segment is linear and ends at initial point,
         // replace it with the closing segment
-        Sequence::iterator last = _curves->end() - 2;
+        Sequence::iterator last = _data->curves.end() - 2;
         if (last->isLineSegment() && last->finalPoint() == initialPoint()) {
             _closing_seg->setInitial(last->initialPoint());
-            _curves->erase(last);
+            _data->curves.erase(last);
         }
     }
     _closed = c;
@@ -302,27 +320,35 @@ void Path::close(bool c)
 void Path::clear()
 {
     _unshare();
-    _curves->pop_back().release();
-    _curves->clear();
+    _data->curves.pop_back().release();
+    _data->curves.clear();
     _closing_seg->setInitial(Point(0, 0));
     _closing_seg->setFinal(Point(0, 0));
-    _curves->push_back(_closing_seg);
+    _data->curves.push_back(_closing_seg);
     _closed = false;
 }
 
 OptRect Path::boundsFast() const
 {
     OptRect bounds;
-    if (empty())
+    if (empty()) {
         return bounds;
+    }
+    // if the path is not empty, we look for cached bounds
+    if (_data->fast_bounds) {
+        return _data->fast_bounds;
+    }
+
     bounds = front().boundsFast();
     const_iterator iter = begin();
-    // the closing path segment can be ignored, because it will always lie within the bbox of the rest of the path
+    // the closing path segment can be ignored, because it will always
+    // lie within the bbox of the rest of the path
     if (iter != end()) {
         for (++iter; iter != end(); ++iter) {
             bounds.unionWith(iter->boundsFast());
         }
     }
+    _data->fast_bounds = bounds;
     return bounds;
 }
 
@@ -377,11 +403,11 @@ bool Path::operator==(Path const &other) const
         return true;
     if (_closed != other._closed)
         return false;
-    return *_curves == *other._curves;
+    return _data->curves == other._data->curves;
 }
 
 void Path::start(Point const &p) {
-    if (_curves->size() > 1) {
+    if (_data->curves.size() > 1) {
         clear();
     }
     _closing_seg->setInitial(p);
@@ -442,86 +468,106 @@ std::vector<PathTime> Path::roots(Coord v, Dim2 d) const
 // Instead of O(N^2), this takes O(N + X), where X is the number of overlaps
 // between the bounding boxes of curves.
 
-struct CurveSweepTraits {
-    struct Bound {
-        Rect r;
+struct CurveIntersectionSweepSet
+{
+public:
+    struct CurveRecord {
+        boost::intrusive::list_member_hook<> _hook;
+        Curve const *curve;
+        Rect bounds;
         std::size_t index;
-        int which;
+        unsigned which;
+
+        CurveRecord(Curve const *pc, std::size_t idx, unsigned w)
+            : curve(pc)
+            , bounds(curve->boundsFast())
+            , index(idx)
+            , which(w)
+        {}
     };
-    typedef std::less<Coord> Compare;
-    inline static Coord entry_value(Bound const &b) { return b.r[X].min(); }
-    inline static Coord exit_value(Bound const &b) { return b.r[X].max(); }
-};
 
-class CurveSweeper
-    : public Sweeper<Curve const *, CurveSweepTraits>
-{
-public:
-    CurveSweeper(Path const &a, Path const &b, std::vector<PathIntersection> &result, Coord prec)
+    typedef std::vector<CurveRecord>::const_iterator ItemIterator;
+
+    CurveIntersectionSweepSet(std::vector<PathIntersection> &result,
+                              Path const &a, Path const &b, Coord precision)
         : _result(result)
-        , _precision(prec)
+        , _precision(precision)
+        , _sweep_dir(X)
     {
-        for (std::size_t i = 0; i < a.size(); ++i) {
-            Bound bound;
-            bound.r = a[i].boundsFast();
-            bound.index = i;
-            bound.which = 0;
-            insert(bound, &a[i]);
+        std::size_t asz = a.size(), bsz = b.size();
+        _records.reserve(asz + bsz);
+
+        for (std::size_t i = 0; i < asz; ++i) {
+            _records.push_back(CurveRecord(&a[i], i, 0));
         }
-        for (std::size_t i = 0; i < b.size(); ++i) {
-            Bound bound;
-            bound.r = b[i].boundsFast();
-            bound.index = i;
-            bound.which = 1;
-            insert(bound, &b[i]);
+        for (std::size_t i = 0; i < bsz; ++i) {
+            _records.push_back(CurveRecord(&b[i], i, 1));
         }
-    }
 
-protected:
-    void _enter(Record const &record) {
-        int which = record.bound.which;
-
-        for (RecordList::iterator i = _active_items.begin(); i != _active_items.end(); ++i) {
-            // do not intersect in the same path
-            if (i->bound.which == which) continue;
-            // do not intersect if boxes do not overlap in Y
-            if (!record.bound.r[Y].intersects(i->bound.r[Y])) continue;
-
-            std::vector<CurveIntersection> cx;
-            int ia = record.bound.index;
-            int ib = i->bound.index;
+        OptRect abb = a.boundsFast() | b.boundsFast();
+        if (abb && abb->height() > abb->width()) {
+            _sweep_dir = Y;
+        }
+    }
 
-            if (which == 0) {
-                cx = record.item->intersect(*i->item, _precision);
-            } else {
-                cx = i->item->intersect(*record.item, _precision);
-                std::swap(ia, ib);
-            }
+    std::vector<CurveRecord> const &items() { return _records; }
+    Interval itemBounds(ItemIterator ii) {
+        return ii->bounds[_sweep_dir];
+    }
 
-            for (std::size_t ci = 0; ci < cx.size(); ++ci) {
-                PathTime a(ia, cx[ci].first), b(ib, cx[ci].second);
-                PathIntersection px(a, b, cx[ci].point());
-                _result.push_back(px);
+    void addActiveItem(ItemIterator ii) {
+        unsigned w = ii->which;
+        unsigned ow = (w+1) % 2;
+
+        _active[w].push_back(const_cast<CurveRecord&>(*ii));
+
+        for (ActiveCurveList::iterator i = _active[ow].begin(); i != _active[ow].end(); ++i) {
+            if (!ii->bounds.intersects(i->bounds)) continue;
+            std::vector<CurveIntersection> cx = ii->curve->intersect(*i->curve, _precision);
+            for (std::size_t k = 0; k < cx.size(); ++k) {
+                PathTime tw(ii->index, cx[k].first), tow(i->index, cx[k].second);
+                _result.push_back(PathIntersection(
+                    w == 0 ? tw : tow,
+                    w == 0 ? tow : tw,
+                    cx[k].point()));
             }
         }
     }
+    void removeActiveItem(ItemIterator ii) {
+        ActiveCurveList &acl = _active[ii->which];
+        acl.erase(acl.iterator_to(*ii));
+    }
 
 private:
+    typedef boost::intrusive::list
+        < CurveRecord
+        , boost::intrusive::member_hook
+            < CurveRecord
+            , boost::intrusive::list_member_hook<>
+            , &CurveRecord::_hook
+            >
+        > ActiveCurveList;
+
+    std::vector<CurveRecord> _records;
     std::vector<PathIntersection> &_result;
+    ActiveCurveList _active[2];
     Coord _precision;
+    Dim2 _sweep_dir;
 };
 
 std::vector<PathIntersection> Path::intersect(Path const &other, Coord precision) const
 {
     std::vector<PathIntersection> result;
 
-    CurveSweeper sweeper(*this, other, result, precision);
+    CurveIntersectionSweepSet cisset(result, *this, other, precision);
+    Sweeper<CurveIntersectionSweepSet> sweeper(cisset);
     sweeper.process();
 
     // preprocessing to remove duplicate intersections at endpoints
+    std::size_t asz = size(), bsz = other.size();
     for (std::size_t i = 0; i < result.size(); ++i) {
-        result[i].first.normalizeForward(size());
-        result[i].second.normalizeForward(other.size());
+        result[i].first.normalizeForward(asz);
+        result[i].second.normalizeForward(bsz);
     }
     std::sort(result.begin(), result.end());
     result.erase(std::unique(result.begin(), result.end()), result.end());
@@ -672,7 +718,7 @@ PathTime Path::nearestTime(Point const &p, Coord *dist) const
     Coord mindist = std::numeric_limits<Coord>::max();
     PathTime ret;
 
-    if (_curves->size() == 1) {
+    if (_data->curves.size() == 1) {
         // naked moveto
         ret.curve_index = 0;
         ret.t = 0;
@@ -701,6 +747,16 @@ PathTime Path::nearestTime(Point const &p, Coord *dist) const
     return ret;
 }
 
+std::vector<Point> Path::nodes() const
+{
+    std::vector<Point> result;
+    size_type path_size = size_closed();
+    for (size_type i = 0; i < path_size; ++i) {
+        result.push_back(_data->curves[i].initialPoint());
+    }
+    return result;
+}
+
 void Path::appendPortionTo(Path &ret, double from, double to) const
 {
     if (!(from >= 0 && to >= 0)) {
@@ -797,16 +853,16 @@ Path Path::reversed() const
     Path ret(finalPoint());
     if (empty()) return ret;
 
-    ret._curves->pop_back(); // this also deletes the closing segment from ret
+    ret._data->curves.pop_back(); // this also deletes the closing segment from ret
 
-    RIter iter(_includesClosingSegment() ? _curves->end() : _curves->end() - 1);
-    RIter rend(_curves->begin());
+    RIter iter(_includesClosingSegment() ? _data->curves.end() : _data->curves.end() - 1);
+    RIter rend(_data->curves.begin());
 
     if (_closed) {
         // when the path is closed, there are two cases:
         if (front().isLineSegment()) {
             // 1. initial segment is linear: it becomes the new closing segment.
-            rend = RIter(_curves->begin() + 1);
+            rend = RIter(_data->curves.begin() + 1);
             ret._closing_seg = new ClosingSegment(front().finalPoint(), front().initialPoint());
         } else {
             // 2. initial segment is not linear: the closing segment becomes degenerate.
@@ -820,9 +876,9 @@ Path Path::reversed() const
     }
 
     for (; iter != rend; ++iter) {
-        ret._curves->push_back(iter->reverse());
+        ret._data->curves.push_back(iter->reverse());
     }
-    ret._curves->push_back(ret._closing_seg);
+    ret._data->curves.push_back(ret._closing_seg);
     ret._closed = _closed;
     return ret;
 }
@@ -896,10 +952,10 @@ void Path::replace(iterator first, iterator last, Path const &path)
 void Path::snapEnds(Coord precision)
 {
     if (!_closed) return;
-    if (_curves->size() > 1 && are_near(_closing_seg->length(precision), 0, precision)) {
+    if (_data->curves.size() > 1 && are_near(_closing_seg->length(precision), 0, precision)) {
         _unshare();
         _closing_seg->setInitial(_closing_seg->finalPoint());
-        (_curves->end() - 1)->setFinal(_closing_seg->finalPoint());
+        (_data->curves.end() - 1)->setFinal(_closing_seg->finalPoint());
     }
 }
 
@@ -908,7 +964,7 @@ void Path::snapEnds(Coord precision)
 void Path::do_update(Sequence::iterator first, Sequence::iterator last, Sequence &source)
 {
     // TODO: handle cases where first > last in closed paths?
-    bool last_beyond_closing_segment = (last == _curves->end());
+    bool last_beyond_closing_segment = (last == _data->curves.end());
 
     // special case:
     // if do_update replaces the closing segment, we have to regenerate it
@@ -916,7 +972,7 @@ void Path::do_update(Sequence::iterator first, Sequence::iterator last, Sequence
         if (first == last) return; // nothing to do
 
         // only removing some segments
-        if ((!_closed && first == _curves->begin()) || (!_closed && last == _curves->end() - 1) || last_beyond_closing_segment) {
+        if ((!_closed && first == _data->curves.begin()) || (!_closed && last == _data->curves.end() - 1) || last_beyond_closing_segment) {
             // just adjust the closing segment
             // do nothing
         } else if (first->initialPoint() != (last - 1)->finalPoint()) {
@@ -927,17 +983,17 @@ void Path::do_update(Sequence::iterator first, Sequence::iterator last, Sequence
         }
     } else {
         // replacing
-        if (first == _curves->begin() && last == _curves->end()) {
+        if (first == _data->curves.begin() && last == _data->curves.end()) {
             // special case: replacing everything should work the same in open and closed curves
-            _curves->erase(_curves->begin(), _curves->end() - 1);
+            _data->curves.erase(_data->curves.begin(), _data->curves.end() - 1);
             _closing_seg->setFinal(source.front().initialPoint());
             _closing_seg->setInitial(source.back().finalPoint());
-            _curves->transfer(_curves->begin(), source.begin(), source.end(), source);
+            _data->curves.transfer(_data->curves.begin(), source.begin(), source.end(), source);
             return;
         }
 
         // stitch in front
-        if (!_closed && first == _curves->begin()) {
+        if (!_closed && first == _data->curves.begin()) {
             // not necessary to stitch in front
         } else if (first->initialPoint() != source.front().initialPoint()) {
             if (_exception_on_stitch) {
@@ -947,7 +1003,7 @@ void Path::do_update(Sequence::iterator first, Sequence::iterator last, Sequence
         }
 
         // stitch at the end
-        if ((!_closed && last == _curves->end() - 1) || last_beyond_closing_segment) {
+        if ((!_closed && last == _data->curves.end() - 1) || last_beyond_closing_segment) {
             // repurpose the closing segment as the stitch segment
             // do nothing
         } else if (source.back().finalPoint() != (last - 1)->finalPoint()) {
@@ -962,8 +1018,8 @@ void Path::do_update(Sequence::iterator first, Sequence::iterator last, Sequence
     if (last_beyond_closing_segment) {
         --last;
     }
-    _curves->erase(first, last);
-    _curves->transfer(first, source.begin(), source.end(), source);
+    _data->curves.erase(first, last);
+    _data->curves.transfer(first, source.begin(), source.end(), source);
 
     // adjust closing segment
     if (size_open() == 0) {
@@ -978,7 +1034,7 @@ void Path::do_update(Sequence::iterator first, Sequence::iterator last, Sequence
 
 void Path::do_append(Curve *c)
 {
-    if (&_curves->front() == _closing_seg) {
+    if (&_data->curves.front() == _closing_seg) {
         _closing_seg->setFinal(c->initialPoint());
     } else {
         // if we can't freely move the closing segment, we check whether
@@ -994,20 +1050,20 @@ void Path::do_append(Curve *c)
             return;
         }
     }
-    _curves->insert(_curves->end() - 1, c);
+    _data->curves.insert(_data->curves.end() - 1, c);
     _closing_seg->setInitial(c->finalPoint());
 }
 
 void Path::checkContinuity() const
 {
-    Sequence::const_iterator i = _curves->begin(), j = _curves->begin();
+    Sequence::const_iterator i = _data->curves.begin(), j = _data->curves.begin();
     ++j;
-    for (; j != _curves->end(); ++i, ++j) {
+    for (; j != _data->curves.end(); ++i, ++j) {
         if (i->finalPoint() != j->initialPoint()) {
             THROW_CONTINUITYERROR();
         }
     }
-    if (_curves->front().initialPoint() != _curves->back().finalPoint()) {
+    if (_data->curves.front().initialPoint() != _data->curves.back().finalPoint()) {
         THROW_CONTINUITYERROR();
     }
 }
@@ -1040,6 +1096,16 @@ Piecewise<D2<SBasis> > paths_to_pw(PathVector const &paths)
     return ret;
 }
 
+bool are_near(Path const &a, Path const &b, Coord precision)
+{
+    if (a.size() != b.size()) return false;
+
+    for (unsigned i = 0; i < a.size(); ++i) {
+        if (!a[i].isNear(b[i], precision)) return false;
+    }
+    return true;
+}
+
 std::ostream &operator<<(std::ostream &out, Path const &path)
 {
     SVGPathWriter pw;
diff --git a/src/2geom/path.h b/src/2geom/path.h
index 3ca43e0e5..3552ec43e 100644
--- a/src/2geom/path.h
+++ b/src/2geom/path.h
@@ -55,6 +55,11 @@ namespace PathInternal {
 
 typedef boost::ptr_vector<Curve> Sequence;
 
+struct PathData {
+    Sequence curves;
+    OptRect fast_bounds;
+};
+
 template <typename P>
 class BaseIterator
     : public boost::random_access_iterator_helper
@@ -221,6 +226,7 @@ public:
     }
 
     size_type pathSize() const { return _path_size; }
+    size_type curveCount() const;
 
 private:
     PathTime _from, _to;
@@ -316,6 +322,7 @@ class Path
     : boost::equality_comparable< Path >
 {
 public:
+    typedef PathInternal::PathData PathData;
     typedef PathInternal::Sequence Sequence;
     typedef PathInternal::BaseIterator<Path> iterator;
     typedef PathInternal::BaseIterator<Path const> const_iterator;
@@ -342,31 +349,31 @@ public:
 
     /// Construct an empty path starting at the specified point.
     explicit Path(Point const &p = Point())
-        : _curves(new Sequence())
+        : _data(new PathData())
         , _closing_seg(new ClosingSegment(p, p))
         , _closed(false)
         , _exception_on_stitch(true)
     {
-        _curves->push_back(_closing_seg);
+        _data->curves.push_back(_closing_seg);
     }
 
     /// Construct a path containing a range of curves.
     template <typename Iter>
     Path(Iter first, Iter last, bool closed = false, bool stitch = false)
-        : _curves(new Sequence())
+        : _data(new PathData())
         , _closed(closed)
         , _exception_on_stitch(!stitch)
     {
         for (Iter i = first; i != last; ++i) {
-            _curves->push_back(i->duplicate());
+            _data->curves.push_back(i->duplicate());
         }
-        if (!_curves->empty()) {
-            _closing_seg = new ClosingSegment(_curves->back().finalPoint(),
-                                              _curves->front().initialPoint());
+        if (!_data->curves.empty()) {
+            _closing_seg = new ClosingSegment(_data->curves.back().finalPoint(),
+                                              _data->curves.front().initialPoint());
         } else {
             _closing_seg = new ClosingSegment();
         }
-        _curves->push_back(_closing_seg);
+        _data->curves.push_back(_closing_seg);
     }
 
     /// Create a path from a rectangle.
@@ -386,7 +393,7 @@ public:
      * @todo Add noexcept specifiers for C++11 */
     void swap(Path &other) throw() {
         using std::swap;
-        swap(other._curves, _curves);
+        swap(other._data, _data);
         swap(other._closing_seg, _closing_seg);
         swap(other._closed, _closed);
         swap(other._exception_on_stitch, _exception_on_stitch);
@@ -396,26 +403,26 @@ public:
     friend inline void swap(Path &a, Path &b) throw() { a.swap(b); }
 
     /** @brief Access a curve by index */
-    Curve const &operator[](size_type i) const { return (*_curves)[i]; }
+    Curve const &operator[](size_type i) const { return _data->curves[i]; }
     /** @brief Access a curve by index */
-    Curve const &at(size_type i) const { return _curves->at(i); }
+    Curve const &at(size_type i) const { return _data->curves.at(i); }
 
     /** @brief Access the first curve in the path.
      * Since the curve always contains at least a degenerate closing segment,
      * it is always safe to use this method. */
-    Curve const &front() const { return _curves->front(); }
+    Curve const &front() const { return _data->curves.front(); }
     /// Alias for front().
-    Curve const &initialCurve() const { return _curves->front(); }
+    Curve const &initialCurve() const { return _data->curves.front(); }
     /** @brief Access the last curve in the path. */
     Curve const &back() const { return back_default(); }
     Curve const &back_open() const {
-        if (empty()) return _curves->back();
-        return (*_curves)[_curves->size() - 2];
+        if (empty()) return _data->curves.back();
+        return _data->curves[_data->curves.size() - 2];
     }
     Curve const &back_closed() const {
         return _closing_seg->isDegenerate()
-            ? (*_curves)[_curves->size() - 2]
-            : (*_curves)[_curves->size() - 1];
+            ? _data->curves[_data->curves.size() - 2]
+            : _data->curves[_data->curves.size() - 1];
     }
     Curve const &back_default() const {
         return _includesClosingSegment()
@@ -436,13 +443,13 @@ public:
     iterator end_closed() { return iterator(*this, size_closed()); }
 
     /// Size without the closing segment, even if the path is closed.
-    size_type size_open() const { return _curves->size() - 1; }
+    size_type size_open() const { return _data->curves.size() - 1; }
 
     /** @brief Size with the closing segment, if it makes a difference.
      * If the closing segment is degenerate, i.e. its initial and final points
      * are exactly equal, then it is not included in this size. */
     size_type size_closed() const {
-        return _closing_seg->isDegenerate() ? _curves->size() - 1 : _curves->size();
+        return _closing_seg->isDegenerate() ? _data->curves.size() - 1 : _data->curves.size();
     }
 
     /// Natural size of the path.
@@ -452,7 +459,7 @@ public:
     /// Natural size of the path.
     size_type size() const { return size_default(); }
 
-    size_type max_size() const { return _curves->max_size() - 1; }
+    size_type max_size() const { return _data->curves.max_size() - 1; }
 
     /** @brief Check whether path is empty.
      * The path is empty if it contains only the closing segment, which according
@@ -460,7 +467,7 @@ public:
      * containers, two empty paths are not necessarily identical, because the
      * degenerate closing segment may be at a different point, affecting the operation
      * of methods such as appendNew(). */
-    bool empty() const { return (_curves->size() == 1); }
+    bool empty() const { return (_data->curves.size() == 1); }
 
     /// Check whether the path is closed.
     bool closed() const { return _closed; }
@@ -497,8 +504,8 @@ public:
     Path &operator*=(T const &tr) {
         BOOST_CONCEPT_ASSERT((TransformConcept<T>));
         _unshare();
-        for (std::size_t i = 0; i < _curves->size(); ++i) {
-            (*_curves)[i] *= tr;
+        for (std::size_t i = 0; i < _data->curves.size(); ++i) {
+            _data->curves[i] *= tr;
         }
         return *this;
     }
@@ -570,6 +577,8 @@ public:
     PathTime nearestTime(Point const &p, Coord *dist = NULL) const;
     std::vector<Coord> nearestTimePerCurve(Point const &p) const;
 
+    std::vector<Point> nodes() const;
+
     void appendPortionTo(Path &p, Coord f, Coord t) const;
 
     /** @brief Append a subset of this path to another path.
@@ -662,13 +671,13 @@ public:
     void setInitial(Point const &p) {
         _unshare();
         _closed = false;
-        _curves->front().setInitial(p);
+        _data->curves.front().setInitial(p);
         _closing_seg->setFinal(p);
     }
     void setFinal(Point const &p) {
         _unshare();
         _closed = false;
-        (*_curves)[size_open() - 1].setFinal(p);
+        _data->curves[size_open() - 1].setFinal(p);
         _closing_seg->setInitial(p);
     }
 
@@ -804,10 +813,10 @@ public:
 
 private:
     static Sequence::iterator seq_iter(iterator const &iter) {
-        return iter.path->_curves->begin() + iter.index;
+        return iter.path->_data->curves.begin() + iter.index;
     }
     static Sequence::const_iterator seq_iter(const_iterator const &iter) {
-        return iter.path->_curves->begin() + iter.index;
+        return iter.path->_data->curves.begin() + iter.index;
     }
 
     // whether the closing segment is part of the path
@@ -815,10 +824,13 @@ private:
         return _closed && !_closing_seg->isDegenerate();
     }
     void _unshare() {
-        if (!_curves.unique()) {
-            _curves.reset(new Sequence(*_curves));
-            _closing_seg = static_cast<ClosingSegment*>(&_curves->back());
+        // Called before every mutation.
+        // Ensure we have our own copy of curve data and reset cached values
+        if (!_data.unique()) {
+            _data.reset(new PathData(*_data));
+            _closing_seg = static_cast<ClosingSegment*>(&_data->curves.back());
         }
+        _data->fast_bounds = OptRect();
     }
     PathTime _factorTime(Coord t) const;
 
@@ -828,7 +840,7 @@ private:
     // n.b. takes ownership of curve object
     void do_append(Curve *curve);
 
-    boost::shared_ptr<Sequence> _curves;
+    boost::shared_ptr<PathData> _data;
     ClosingSegment *_closing_seg;
     bool _closed;
     bool _exception_on_stitch;
@@ -841,6 +853,8 @@ inline Coord nearest_time(Point const &p, Path const &c) {
     return pt.curve_index + pt.t;
 }
 
+bool are_near(Path const &a, Path const &b, Coord precision = EPSILON);
+
 std::ostream &operator<<(std::ostream &out, Path const &path);
 
 } // end namespace Geom
diff --git a/src/2geom/pathvector.cpp b/src/2geom/pathvector.cpp
index bff201c71..28ab26237 100644
--- a/src/2geom/pathvector.cpp
+++ b/src/2geom/pathvector.cpp
@@ -35,6 +35,7 @@
 #include <2geom/path.h>
 #include <2geom/pathvector.h>
 #include <2geom/svg-path-writer.h>
+#include <2geom/sweeper.h>
 
 namespace Geom {
 
@@ -133,19 +134,96 @@ void PathVector::snapEnds(Coord precision)
     }
 }
 
-std::vector<PVIntersection> PathVector::intersect(PathVector const &other, Coord precision) const
-{
-    typedef PathVectorTime PVPos;
-    std::vector<PVIntersection> result;
-    for (std::size_t i = 0; i < size(); ++i) {
-        for (std::size_t j = 0; j < other.size(); ++j) {
-            std::vector<PathIntersection> xs = (*this)[i].intersect(other[j], precision);
-            for (std::size_t k = 0; k < xs.size(); ++k) {
-                PVIntersection pvx(PVPos(i, xs[k].first), PVPos(j, xs[k].second), xs[k].point());
-                result.push_back(pvx);
+// sweepline optimization
+// this is very similar to CurveIntersectionSweepSet in path.cpp
+// should probably be merged
+class PathIntersectionSweepSet {
+public:
+    struct PathRecord {
+        boost::intrusive::list_member_hook<> _hook;
+        Path const *path;
+        std::size_t index;
+        unsigned which;
+
+        PathRecord(Path const &p, std::size_t i, unsigned w)
+            : path(&p)
+            , index(i)
+            , which(w)
+        {}
+    };
+
+    typedef std::vector<PathRecord>::iterator ItemIterator;
+
+    PathIntersectionSweepSet(std::vector<PVIntersection> &result,
+                             PathVector const &a, PathVector const &b, Coord precision)
+        : _result(result)
+        , _precision(precision)
+    {
+        _records.reserve(a.size() + b.size());
+        for (std::size_t i = 0; i < a.size(); ++i) {
+            _records.push_back(PathRecord(a[i], i, 0));
+        }
+        for (std::size_t i = 0; i < b.size(); ++i) {
+            _records.push_back(PathRecord(b[i], i, 1));
+        }
+    }
+
+    std::vector<PathRecord> &items() { return _records; }
+
+    Interval itemBounds(ItemIterator ii) {
+        OptRect r = ii->path->boundsFast();
+        return (*r)[X];
+    }
+
+    void addActiveItem(ItemIterator ii) {
+        unsigned w = ii->which;
+        unsigned ow = (ii->which + 1) % 2;
+
+        for (ActivePathList::iterator i = _active[ow].begin(); i != _active[ow].end(); ++i) {
+            if (!ii->path->boundsFast().intersects(i->path->boundsFast())) continue;
+            std::vector<PathIntersection> px = ii->path->intersect(*i->path, _precision);
+            for (std::size_t k = 0; k < px.size(); ++k) {
+                PathVectorTime tw(ii->index, px[k].first), tow(i->index, px[k].second);
+                _result.push_back(PVIntersection(
+                    w == 0 ? tw : tow,
+                    w == 0 ? tow : tw,
+                    px[k].point()));
             }
         }
+        _active[w].push_back(*ii);
+    }
+
+    void removeActiveItem(ItemIterator ii) {
+        ActivePathList &apl = _active[ii->which];
+        apl.erase(apl.iterator_to(*ii));
     }
+
+private:
+    typedef boost::intrusive::list
+        < PathRecord
+        , boost::intrusive::member_hook
+            < PathRecord
+            , boost::intrusive::list_member_hook<>
+            , &PathRecord::_hook
+            >
+        > ActivePathList;
+
+    std::vector<PVIntersection> &_result;
+    std::vector<PathRecord> _records;
+    ActivePathList _active[2];
+    Coord _precision;
+};
+
+std::vector<PVIntersection> PathVector::intersect(PathVector const &other, Coord precision) const
+{
+    std::vector<PVIntersection> result;
+
+    PathIntersectionSweepSet pisset(result, *this, other, precision);
+    Sweeper<PathIntersectionSweepSet> sweeper(pisset);
+    sweeper.process();
+
+    std::sort(result.begin(), result.end());
+
     return result;
 }
 
@@ -153,6 +231,7 @@ int PathVector::winding(Point const &p) const
 {
     int wind = 0;
     for (const_iterator i = begin(); i != end(); ++i) {
+        if (!i->boundsFast().contains(p)) continue;
         wind += i->winding(p);
     }
     return wind;
@@ -201,6 +280,18 @@ std::vector<PathVectorTime> PathVector::allNearestTimes(Point const &p, Coord *d
     return retval;
 }
 
+std::vector<Point> PathVector::nodes() const
+{
+    std::vector<Point> result;
+    for (size_type i = 0; i < size(); ++i) {
+        size_type path_size = (*this)[i].size_closed();
+        for (size_type j = 0; j < path_size; ++j) {
+            result.push_back((*this)[i][j].initialPoint());
+        }
+    }
+    return result;
+}
+
 PathVectorTime PathVector::_factorTime(Coord t) const
 {
     PathVectorTime ret;
diff --git a/src/2geom/pathvector.h b/src/2geom/pathvector.h
index 108f2aa05..0cbe5bf52 100644
--- a/src/2geom/pathvector.h
+++ b/src/2geom/pathvector.h
@@ -272,6 +272,8 @@ public:
     boost::optional<PathVectorTime> nearestTime(Point const &p, Coord *dist = NULL) const;
     std::vector<PathVectorTime> allNearestTimes(Point const &p, Coord *dist = NULL) const;
 
+    std::vector<Point> nodes() const;
+
 private:
     PathVectorTime _factorTime(Coord t) const;
 
diff --git a/src/2geom/piecewise.h b/src/2geom/piecewise.h
index a5a65a9fe..3a12671e0 100644
--- a/src/2geom/piecewise.h
+++ b/src/2geom/piecewise.h
@@ -222,7 +222,7 @@ class Piecewise {
     inline void setDomain(Interval dom) {
         if(empty()) return;
         /* dom can not be empty
-        if(dom.isEmpty()) {
+        if(dom.empty()) {
             cuts.clear(); segs.clear();
             return;
         }*/
diff --git a/src/2geom/ray.h b/src/2geom/ray.h
index a336e3132..4e60fd814 100644
--- a/src/2geom/ray.h
+++ b/src/2geom/ray.h
@@ -53,61 +53,62 @@ namespace Geom
 class Ray {
 private:
     Point _origin;
-    Point _versor;
+    Point _vector;
 
 public:
-    Ray() : _origin(0,0), _versor(1,0) {}
+    Ray() : _origin(0,0), _vector(1,0) {}
     Ray(Point const& origin, Coord angle)
         : _origin(origin)
     {
-        sincos(angle, _versor[Y], _versor[X]);
+        sincos(angle, _vector[Y], _vector[X]);
     }
     Ray(Point const& A, Point const& B) {
         setPoints(A, B);
     }
     Point origin() const { return _origin; }
-    Point versor() const { return _versor; }
+    Point vector() const { return _vector; }
+    Point versor() const { return _vector.normalized(); }
     void setOrigin(Point const &o) { _origin = o; }
-    void setVersor(Point const& v) { _versor = v; }
-    Coord angle() const { return std::atan2(_versor[Y], _versor[X]); }
-    void setAngle(Coord a) { sincos(a, _versor[Y], _versor[X]); }
+    void setVector(Point const& v) { _vector = v; }
+    Coord angle() const { return std::atan2(_vector[Y], _vector[X]); }
+    void setAngle(Coord a) { sincos(a, _vector[Y], _vector[X]); }
     void setPoints(Point const &a, Point const &b) {
         _origin = a;
-        _versor = b - a;
-        if (are_near(_versor, Point(0,0)) )
-            _versor = Point(0,0);
+        _vector = b - a;
+        if (are_near(_vector, Point(0,0)) )
+            _vector = Point(0,0);
         else
-            _versor.normalize();
+            _vector.normalize();
     }
     bool isDegenerate() const {
-        return ( _versor[X] == 0 && _versor[Y] == 0 );
+        return ( _vector[X] == 0 && _vector[Y] == 0 );
     }
     Point pointAt(Coord t) const {
-        return _origin + _versor * t;
+        return _origin + _vector * t;
     }
     Coord valueAt(Coord t, Dim2 d) const {
-        return _origin[d] + _versor[d] * t;
+        return _origin[d] + _vector[d] * t;
     }
     std::vector<Coord> roots(Coord v, Dim2 d) const {
         std::vector<Coord> result;
-        if ( _versor[d] != 0 ) {
-            double t = (v - _origin[d]) / _versor[d];
+        if ( _vector[d] != 0 ) {
+            double t = (v - _origin[d]) / _vector[d];
             if (t >= 0)	result.push_back(t);
-        } else if (_versor[(d+1)%2] == v) {
+        } else if (_vector[(d+1)%2] == v) {
             THROW_INFINITESOLUTIONS();
         }
         return result;
     }
     Coord nearestTime(Point const& point) const {
         if ( isDegenerate() ) return 0;
-        double t = dot(point - _origin, _versor);
+        double t = dot(point - _origin, _vector);
         if (t < 0) t = 0;
         return t;
     }
     Ray reverse() const {
         Ray result;
         result.setOrigin(_origin);
-        result.setVersor(-_versor);
+        result.setVector(-_vector);
         return result;
     }
     Curve *portion(Coord f, Coord t) const {
@@ -117,7 +118,7 @@ public:
         return LineSegment(pointAt(f), pointAt(t));
     }
     Ray transformed(Affine const& m) const {
-        return Ray(_origin * m, (_origin + _versor) * m);
+        return Ray(_origin * m, (_origin + _vector) * m);
     }
 }; // end class Ray
 
@@ -134,7 +135,7 @@ bool are_near(Point const& _point, Ray const& _ray, double eps = EPSILON) {
 
 inline
 bool are_same(Ray const& r1, Ray const& r2, double eps = EPSILON) {
-	return are_near(r1.versor(), r2.versor(), eps)
+	return are_near(r1.vector(), r2.vector(), eps)
 			&& are_near(r1.origin(), r2.origin(), eps);
 }
 
@@ -142,7 +143,7 @@ bool are_same(Ray const& r1, Ray const& r2, double eps = EPSILON) {
 // the returned value is an angle in the interval [0, 2PI[
 inline
 double angle_between(Ray const& r1, Ray const& r2, bool cw = true) {
-	double angle = angle_between(r1.versor(), r2.versor());
+	double angle = angle_between(r1.vector(), r2.vector());
 	if (angle < 0) angle += 2*M_PI;
 	if (!cw) angle = 2*M_PI - angle;
 	return angle;
@@ -170,7 +171,7 @@ Ray make_angle_bisector_ray(Ray const& r1, Ray const& r2, bool cw = true)
         THROW_RANGEERROR("passed rays do not have the same origin");
     }
 
-    Ray bisector(r1.origin(), r1.origin() +  r1.versor() * Rotate(angle_between(r1, r2) / 2.0));
+    Ray bisector(r1.origin(), r1.origin() +  r1.vector() * Rotate(angle_between(r1, r2) / 2.0));
     
     return (cw ? bisector : bisector.reverse());
 }
diff --git a/src/2geom/rect.cpp b/src/2geom/rect.cpp
index 383e72c8e..5a821e9f9 100644
--- a/src/2geom/rect.cpp
+++ b/src/2geom/rect.cpp
@@ -30,9 +30,56 @@
  */
 
 #include <2geom/rect.h>
+#include <2geom/transforms.h>
 
 namespace Geom {
 
+Point align_factors(Align g) {
+    Point p;
+    switch (g) {
+    case ALIGN_XMIN_YMIN:
+        p[X] = 0.0;
+        p[Y] = 0.0;
+        break;
+    case ALIGN_XMID_YMIN:
+        p[X] = 0.5;
+        p[Y] = 0.0;
+        break;
+    case ALIGN_XMAX_YMIN:
+        p[X] = 1.0;
+        p[Y] = 0.0;
+        break;
+    case ALIGN_XMIN_YMID:
+        p[X] = 0.0;
+        p[Y] = 0.5;
+        break;
+    case ALIGN_XMID_YMID:
+        p[X] = 0.5;
+        p[Y] = 0.5;
+        break;
+    case ALIGN_XMAX_YMID:
+        p[X] = 1.0;
+        p[Y] = 0.5;
+        break;
+    case ALIGN_XMIN_YMAX:
+        p[X] = 0.0;
+        p[Y] = 1.0;
+        break;
+    case ALIGN_XMID_YMAX:
+        p[X] = 0.5;
+        p[Y] = 1.0;
+        break;
+    case ALIGN_XMAX_YMAX:
+        p[X] = 1.0;
+        p[Y] = 1.0;
+        break;
+    default:
+        break;
+    }
+    return p;
+}
+
+
 /** @brief Transform the rectangle by an affine.
  * The result of the transformation might not be axis-aligned. The return value
  * of this operation will be the smallest axis-aligned rectangle containing
@@ -49,6 +96,37 @@ Rect &Rect::operator*=(Affine const &m) {
     return *this;
 }
 
+Affine Rect::transformTo(Rect const &viewport, Aspect const &aspect) const
+{
+    // 1. translate viewbox to origin
+    Geom::Affine total = Translate(-min());
+
+    // 2. compute scale
+    Geom::Point vdims = viewport.dimensions();
+    Geom::Point dims = dimensions();
+    Geom::Scale scale(vdims[X] / dims[X], vdims[Y] / dims[Y]);
+
+    if (aspect.align == ALIGN_NONE) {
+        // apply non-uniform scale
+        total *= scale * Translate(viewport.min());
+    } else {
+        double uscale = 0;
+        if (aspect.expansion == EXPANSION_MEET) {
+            uscale = std::min(scale[X], scale[Y]);
+        } else {
+            uscale = std::max(scale[X], scale[Y]);
+        }
+        scale = Scale(uscale);
+
+        // compute offset for align
+        Geom::Point offset = vdims - dims * scale;
+        offset *= Scale(align_factors(aspect.align));
+        total *= scale * Translate(viewport.min() + offset);
+    }
+
+    return total;
+}
+
 Coord distanceSq(Point const &p, Rect const &rect)
 {
     double dx = 0, dy = 0;
diff --git a/src/2geom/rect.h b/src/2geom/rect.h
index 51daf6b5a..d86741476 100644
--- a/src/2geom/rect.h
+++ b/src/2geom/rect.h
@@ -47,6 +47,43 @@
 
 namespace Geom {
 
+/** Values for the <align> parameter of preserveAspectRatio.
+ * See: http://www.w3.org/TR/SVG/coords.html#PreserveAspectRatioAttribute */
+enum Align {
+    ALIGN_NONE,
+    ALIGN_XMIN_YMIN,
+    ALIGN_XMID_YMIN,
+    ALIGN_XMAX_YMIN,
+    ALIGN_XMIN_YMID,
+    ALIGN_XMID_YMID,
+    ALIGN_XMAX_YMID,
+    ALIGN_XMIN_YMAX,
+    ALIGN_XMID_YMAX,
+    ALIGN_XMAX_YMAX
+};
+
+/** Values for the <meetOrSlice> parameter of preserveAspectRatio.
+ * See: http://www.w3.org/TR/SVG/coords.html#PreserveAspectRatioAttribute */
+enum Expansion {
+    EXPANSION_MEET,
+    EXPANSION_SLICE
+};
+
+/// Convert an align specification to coordinate fractions.
+Point align_factors(Align align);
+
+/** @brief Structure that specifies placement of within a viewport.
+ * Use this to create transformations that preserve aspect. */
+struct Aspect {
+    Align align;
+    Expansion expansion;
+    bool deferred; ///< for SVG compatibility
+
+    Aspect(Align a = ALIGN_NONE, Expansion ex = EXPANSION_MEET)
+        : align(a), expansion(ex), deferred(false)
+    {}
+};
+
 /**
  * @brief Axis aligned, non-empty rectangle.
  * @ingroup Primitives
@@ -113,6 +150,16 @@ public:
     }
     /// @}
 
+    /// @name SVG viewbox functionality.
+    /// @{
+    /** @brief Transform contents to viewport.
+     * Computes an affine that transforms the contents of this rectangle
+     * to the specified viewport. The aspect parameter specifies how to
+     * to the transformation (whether the aspect ratio of content
+     * should be kept and where it should be placed in the viewport). */
+    Affine transformTo(Rect const &viewport, Aspect const &aspect = Aspect()) const;
+    /// @}
+
     /// @name Operators
     /// @{
     Rect &operator*=(Affine const &m);
@@ -145,8 +192,14 @@ public:
     OptRect(OptIntRect const &r) : Base() {
         if (r) *this = Rect(*r);
     }
-    // actually, the only reason we have this class, instead of typedefing
-    // to GenericOptRect<Coord>, are the above constructors
+
+    Affine transformTo(Rect const &viewport, Aspect const &aspect = Aspect()) {
+        Affine ret = Affine::identity();
+        if (empty()) return ret;
+        ret = (*this)->transformTo(viewport, aspect);
+        return ret;
+    }
+
     bool operator==(OptRect const &other) const {
         return Base::operator==(other);
     }
diff --git a/src/2geom/sbasis-curve.h b/src/2geom/sbasis-curve.h
index affe7edc0..cfc4ee9a8 100644
--- a/src/2geom/sbasis-curve.h
+++ b/src/2geom/sbasis-curve.h
@@ -37,6 +37,7 @@
 #define LIB2GEOM_SEEN_SBASIS_CURVE_H
 
 #include <2geom/curve.h>
+#include <2geom/exception.h>
 #include <2geom/nearest-time.h>
 #include <2geom/sbasis-geometric.h>
 #include <2geom/transforms.h>
@@ -131,6 +132,10 @@ public:
         if (!other) return false;
         return inner == other->inner;
     }
+    virtual bool isNear(Curve const &/*c*/, Coord /*eps*/) const {
+        THROW_NOTIMPLEMENTED();
+        return false;
+    }
     virtual int degreesOfFreedom() const {
         return inner[0].degreesOfFreedom() + inner[1].degreesOfFreedom();
     }
diff --git a/src/2geom/sbasis-geometric.cpp b/src/2geom/sbasis-geometric.cpp
index 4c474f7f0..8aaa15144 100644
--- a/src/2geom/sbasis-geometric.cpp
+++ b/src/2geom/sbasis-geometric.cpp
@@ -228,7 +228,7 @@ Geom::unitVector(D2<SBasis> const &V_in, double tol, unsigned order){
     // -This done, unitVector will have jumps at zeros: fill the gaps with arcs of circles.
     D2<SBasis> V = RescaleForNonVanishingEnds(V_in);
 
-    if (V[0].empty() && V[1].empty())
+    if (V[0].isZero(0) && V[1].isZero(0))
         return Piecewise<D2<SBasis> >(D2<SBasis>(Linear(1),SBasis()));
     SBasis x = V[0], y = V[1];
     SBasis r_eqn1, r_eqn2;
diff --git a/src/2geom/sbasis-math.cpp b/src/2geom/sbasis-math.cpp
index 896eb18a7..e9ee917de 100644
--- a/src/2geom/sbasis-math.cpp
+++ b/src/2geom/sbasis-math.cpp
@@ -34,9 +34,8 @@
 //TODO: define a truncated compose(sb,sb, order) and extend it to pw<sb>.
 //TODO: in all these functions, compute 'order' according to 'tol'.
 
+#include <2geom/d2.h>
 #include <2geom/sbasis-math.h>
-
-#include <2geom/d2-sbasis.h>
 #include <stdio.h>
 #include <math.h>
 //#define ZERO 1e-3
diff --git a/src/2geom/sbasis-roots.cpp b/src/2geom/sbasis-roots.cpp
index 57bef4c0f..e3e5e4441 100644
--- a/src/2geom/sbasis-roots.cpp
+++ b/src/2geom/sbasis-roots.cpp
@@ -222,7 +222,7 @@ static void multi_roots_internal(SBasis const &f,
 				 double b,
 				 double fb){
 
-    if (f.size()==0){
+    if (f.isZero(0)){
         int idx;
         idx=upper_level(levels,0,vtol);
         if (idx<(int)levels.size()&&fabs(levels.at(idx))<=vtol){
@@ -414,7 +414,7 @@ static void level_sets_internal(SBasis const &f,
 				 double fb,
 				 double tol=1e-5){
 
-    if (f.size()==0){
+    if (f.isZero(0)){
         unsigned idx;
         idx=upper_level( levels, 0. );
         if (idx<levels.size() && levels[idx].contains(0.)){
@@ -614,6 +614,7 @@ std::vector<double> roots1(SBasis const & s, Interval const ivl) {
 std::vector<double> roots(SBasis const & s) {
     switch(s.size()) {
         case 0:
+            assert(false);
             return std::vector<double>();
         case 1:
             return roots1(s);
@@ -628,6 +629,7 @@ std::vector<double> roots(SBasis const & s) {
 std::vector<double> roots(SBasis const & s, Interval const ivl) {
     switch(s.size()) {
         case 0:
+        	assert(false);
             return std::vector<double>();
         case 1:
             return roots1(s, ivl);
diff --git a/src/2geom/sbasis-to-bezier.cpp b/src/2geom/sbasis-to-bezier.cpp
index 09fbb03ef..8a18cfd4a 100644
--- a/src/2geom/sbasis-to-bezier.cpp
+++ b/src/2geom/sbasis-to-bezier.cpp
@@ -100,9 +100,7 @@ int sgn(unsigned int j, unsigned int k)
 */
 void sbasis_to_bezier (Bezier & bz, SBasis const& sb, size_t sz)
 {
-    if (sb.size() == 0) {
-        THROW_RANGEERROR("size of sb is too small");
-    }
+    assert(sb.size() > 0);
 
     size_t q, n;
     bool even;
@@ -207,11 +205,11 @@ void sbasis_to_cubic_bezier (std::vector<Point> & bz, D2<SBasis> const& sb)
         xprime[i] = sb[X][0][1] - sb[X][0][0];
         yprime[i] = sb[Y][0][1] - sb[Y][0][0];
     }
-    if (sb[X].size() > 0) {
+    if (sb[X].size() > 1) {
         xprime[0] += sb[X][1][0];
         xprime[1] -= sb[X][1][1];
     }
-    if (sb[Y].size() > 0) {
+    if (sb[Y].size() > 1) {
         yprime[0] += sb[Y][1][0];
         yprime[1] -= sb[Y][1][1];
     }
diff --git a/src/2geom/sbasis.cpp b/src/2geom/sbasis.cpp
index 42d92d7b8..0f4159e94 100644
--- a/src/2geom/sbasis.cpp
+++ b/src/2geom/sbasis.cpp
@@ -279,9 +279,11 @@ SBasis multiply_add(SBasis const &a, SBasis const &b, SBasis c) {
 
 */
 SBasis multiply(SBasis const &a, SBasis const &b) {
-    SBasis c(a.size() + b.size(), Linear(0,0));
-    if(a.isZero() || b.isZero())
+    if(a.isZero() || b.isZero()) {
+        SBasis c(1, Linear(0,0));
         return c;
+    }
+    SBasis c(a.size() + b.size(), Linear(0,0));
     return multiply_add(a, b, c);
 }
 #endif 
diff --git a/src/2geom/sbasis.h b/src/2geom/sbasis.h
index 787e8b722..6923017be 100644
--- a/src/2geom/sbasis.h
+++ b/src/2geom/sbasis.h
@@ -83,49 +83,52 @@ public:
     const_iterator end() const { return d.end();}
     iterator begin() { return d.begin();}
     iterator end() { return d.end();}
-    bool empty() const {return d.empty();}
+    bool empty() const { return d.size() == 1 && d[0][0] == 0 && d[0][1] == 0; }
     Linear &back() {return d.back();}
     Linear const &back() const {return d.back();}
-    void pop_back() { d.pop_back();}
-    void resize(unsigned n) { d.resize(n);}
-    void resize(unsigned n, Linear const& l) { d.resize(n, l);}
+    void pop_back() {
+		if (d.size() > 1) {
+			d.pop_back();
+		} else {
+			d[0][0] = 0;
+			d[0][1] = 0;
+		}
+	}
+    void resize(unsigned n) { d.resize(std::max<unsigned>(n, 1));}
+    void resize(unsigned n, Linear const& l) { d.resize(std::max<unsigned>(n, 1), l);}
     void reserve(unsigned n) { d.reserve(n);}
-    void clear() {d.clear();}
+    void clear() {
+    	d.resize(1);
+    	d[0][0] = 0;
+    	d[0][1] = 0;
+    }
     void insert(iterator before, const_iterator src_begin, const_iterator src_end) { d.insert(before, src_begin, src_end);}
     Linear& at(unsigned i) { return d.at(i);}
     //void insert(Linear* before, int& n, Linear const &l) { d.insert(std::vector<Linear>::iterator(before), n, l);}
     bool operator==(SBasis const&B) const { return d == B.d;}
     bool operator!=(SBasis const&B) const { return d != B.d;}
-    operator std::vector<Linear>() { return d;}
-
     
-    SBasis() {}
+    SBasis()
+		: d(1, Linear(0, 0))
+	{}
     explicit SBasis(double a)
-        : d(1)
-    {
-        d[0][0] = a;
-        d[0][1] = a;
-    }
+        : d(1, Linear(a, a))
+    {}
     explicit SBasis(double a, double b)
-        : d(1)
-    {
-        d[0][0] = a;
-        d[0][1] = b;
-    }
-    SBasis(SBasis const & a) :
-        d(a.d)
+        : d(1, Linear(a, b))
     {}
-    SBasis(std::vector<Linear> const & ls) :
-        d(ls)
+    SBasis(SBasis const &a)
+        : d(a.d)
+    {}
+    SBasis(std::vector<Linear> const &ls)
+        : d(ls)
+    {}
+    SBasis(Linear const &bo)
+		: d(1, bo)
+	{}
+    SBasis(Linear* bo)
+		: d(1, bo ? *bo : Linear(0, 0))
     {}
-    SBasis(Linear const & bo) {
-        push_back(bo);
-    }
-    SBasis(Linear* bo) {
-        if (bo) {
-            push_back(*bo);
-        }
-    }
     explicit SBasis(size_t n, Linear const&l) : d(n, l) {}
 
     SBasis(Coord c0, Coord c1, Coord c2, Coord c3)
@@ -179,6 +182,7 @@ public:
     template <typename Iter>
     SBasis(Iter first, Iter last) {
         assert(std::distance(first, last) % 2 == 0);
+        assert(std::distance(first, last) >= 2);
         for (; first != last; ++first) {
             --last;
             push_back(Linear(*first, *last));
@@ -188,14 +192,14 @@ public:
     //IMPL: FragmentConcept
     typedef double output_type;
     inline bool isZero(double eps=EPSILON) const {
-        if(empty()) return true;
+    	assert(size() > 0);
         for(unsigned i = 0; i < size(); i++) {
             if(!(*this)[i].isZero(eps)) return false;
         }
         return true;
     }
     inline bool isConstant(double eps=EPSILON) const {
-        if (empty()) return true;
+    	assert(size() > 0);
         if(!(*this)[0].isConstant(eps)) return false;
         for (unsigned i = 1; i < size(); i++) {
             if(!(*this)[i].isZero(eps)) return false;
@@ -212,6 +216,7 @@ public:
     int degreesOfFreedom() const { return size()*2;}
 
     double valueAt(double t) const {
+    	assert(size() > 0);
         double s = t*(1-t);
         double p0 = 0, p1 = 0;
         for(unsigned k = size(); k > 0; k--) {
@@ -239,11 +244,11 @@ public:
 //MUTATOR PRISON
     //remove extra zeros
     void normalize() {
-        while(!empty() && 0 == back()[0] && 0 == back()[1])
+        while(size() > 1 && back().isZero(0))
             pop_back();
     }
 
-    void truncate(unsigned k) { if(k < size()) resize(k); }
+    void truncate(unsigned k) { if(k < size()) resize(std::max<size_t>(k, 1)); }
 private:
     void derive(); // in place version
 };
diff --git a/src/2geom/svg-path-parser.cpp b/src/2geom/svg-path-parser.cpp
index b6e6da869..93694156d 100644
--- a/src/2geom/svg-path-parser.cpp
+++ b/src/2geom/svg-path-parser.cpp
@@ -1424,7 +1424,7 @@ _match:
             Point point = _pop_point();
             bool sweep = _pop_flag();
             bool large_arc = _pop_flag();
-            double angle = deg_to_rad(_pop());
+            double angle = rad_from_deg(_pop());
             double ry = _pop();
             double rx = _pop();
 
@@ -1530,7 +1530,7 @@ _again:
             Point point = _pop_point();
             bool sweep = _pop_flag();
             bool large_arc = _pop_flag();
-            double angle = deg_to_rad(_pop());
+            double angle = rad_from_deg(_pop());
             double ry = _pop();
             double rx = _pop();
 
diff --git a/src/2geom/svg-path-writer.cpp b/src/2geom/svg-path-writer.cpp
index 1c40ba64e..c484a172b 100644
--- a/src/2geom/svg-path-writer.cpp
+++ b/src/2geom/svg-path-writer.cpp
@@ -150,7 +150,7 @@ void SVGPathWriter::arcTo(double rx, double ry, double angle,
     _setCommand('A');
     _current_pars.push_back(rx);
     _current_pars.push_back(ry);
-    _current_pars.push_back(rad_to_deg(angle));
+    _current_pars.push_back(deg_from_rad(angle));
     _current_pars.push_back(large_arc ? 1. : 0.);
     _current_pars.push_back(sweep ? 1. : 0.);
     _current_pars.push_back(p[X]);
diff --git a/src/2geom/sweeper.h b/src/2geom/sweeper.h
index 8c4e182a6..469108465 100644
--- a/src/2geom/sweeper.h
+++ b/src/2geom/sweeper.h
@@ -41,32 +41,24 @@
 
 namespace Geom {
 
-/** @brief Sweep traits class for interval bounds.
- * @relates Sweeper
- * @ingroup Utilities */
-struct IntervalSweepTraits {
-    typedef Interval Bound;
-    typedef std::less<Coord> Compare;
-    inline static Coord entry_value(Bound const &b) { return b.min(); }
-    inline static Coord exit_value(Bound const &b) { return b.max(); }
-};
+// exposition only
+template <typename Item>
+class SweepVector {
+public:
+    typedef typename std::vector<Item>::const_iterator ItemIterator;
 
-/** @brief Sweep traits class for rectangle bounds.
- * @tparam D Which axis to use for sweeping
- * @ingroup Utilities */
-template <Dim2 D>
-struct RectSweepTraits {
-    typedef Rect Bound;
-    typedef std::less<Coord> Compare;
-    inline static Coord entry_value(Bound const &b) { return b[D].min(); }
-    inline static Coord exit_value(Bound const &b) { return b[D].max(); }
-};
+    SweepVector(std::vector<Item> const &v)
+        : _items(v)
+    {}
 
-template <typename T>
-struct BoundsFast {
-    Rect operator()(T const &item) const {
-        return item.boundsFast();
-    }
+    std::vector<Item> const &items() { return _items; }
+    Interval itemBounds(ItemIterator /*ii*/) { return Interval(); }
+
+    void addActiveItem(ItemIterator /*ii*/) {}
+    void removeActiveItem(ItemIterator /*ii*/) {}
+
+private:
+    std::vector<Item> const &_items;
 };
 
 /** @brief Generic sweepline algorithm.
@@ -77,13 +69,19 @@ struct BoundsFast {
  * the line starts intersecting their bounds, and removed when it completely
  * passes over them.
  *
- * To use this, create a derived class and reimplement the _enter()
- * and/or _leave() virtual functions, insert all the objects,
- * and finally call process(). Inside _enter() and _leave(), the items that have
- * their bounds intersected by the sweepline are available in a list called
- * _active_items. This is an intrusive linked list, so you should access it using
- * iterators. Do not add or remove items from it. You can specify the bound type
- * and how it should be accessed by defining a custom SweepTraits class.
+ * To use this, create a class that exposes the following methods:
+ * - Range items() - returns a forward iterable range of items that will be swept.
+ * - Interval itemBounds(iterator i) - given an iterator from the above range,
+ *   compute the bounding interval of the referenced item in the direction of sweep.
+ * - void addActiveItem(iterator i) - add an item to the active list.
+ * - void removeActiveItem(iterator i) - remove an item from the active list.
+ *
+ * Create the object, then instantiate this template with the above class
+ * as the template parameter, pass it the constructed object of the class,
+ * and call the process() method.
+ *
+ * A good choice for the active list is a Boost intrusive list, which allows
+ * you to get an iterator from a value in constant time.
  *
  * Look in path.cpp for example usage.
  *
@@ -92,148 +90,88 @@ struct BoundsFast {
  *   how to interpret them and how to sort the events
  * @ingroup Utilities
  */
-template <typename Item, typename SweepTraits = IntervalSweepTraits>
+template <typename SweepSet>
 class Sweeper {
 public:
-    /// Type of the item's boundary - usually this will be an Interval or Rect.
-    typedef typename SweepTraits::Bound Bound;
-
-    Sweeper() {}
-
-    /** @brief Insert a single item for sweeping.
-     * @param bound Boundary of the item, as defined in sweep traits
-     * @param item The item itself */
-    void insert(Bound const &bound, Item const &item) {
-        assert(!(typename SweepTraits::Compare()(
-            SweepTraits::exit_value(bound),
-            SweepTraits::entry_value(bound))));
-        _items.push_back(Record(bound, item));
-    }
-
-    /** @brief Insert a range of items using the supplied bounds functor.
-     * The bounds are computed from items using the supplied bounds functor.
-     * @param first Start of range
-     * @param last End of range (one-past-the-end iterator)
-     * @param f Bounds functor */
-    template <typename Iter, typename BoundFunc>
-    void insert(Iter first, Iter last, BoundFunc f = BoundFunc()) {
-        for (; first != last; ++first) {
-            Bound b = f(*first);
-            assert(!(typename SweepTraits::Compare()(
-                SweepTraits::exit_value(b),
-                SweepTraits::entry_value(b))));
-            _items.push_back(Record(b, *first));
-        }
+    typedef typename SweepSet::ItemIterator Iter;
+
+    explicit Sweeper(SweepSet &set)
+        : _set(set)
+    {
+        std::size_t sz = std::distance(set.items().begin(), set.items().end());
+        _entry_events.reserve(sz);
+        _exit_events.reserve(sz);
     }
 
     /** @brief Process entry and exit events.
-     * This will iterate over all inserted items, calling the virtual protected
-     * functions _enter() and _leave() according to the order of the boundaries
-     * of each item. */
+     * This will iterate over all inserted items, calling the methods
+     * addActiveItem and removeActiveItem on the SweepSet passed at construction
+     * according to the order of the boundaries of each item. */
     void process() {
-        if (_items.empty()) return;
+        if (_set.items().empty()) return;
+
+        Iter last = _set.items().end();
+        for (Iter i = _set.items().begin(); i != last; ++i) {
+            Interval b = _set.itemBounds(i);
+            // guard against NANs
+            assert(b.min() == b.min() && b.max() == b.max());
+            _entry_events.push_back(Event(b.max(), i));
+            _exit_events.push_back(Event(b.min(), i));
+        }
 
-        typename SweepTraits::Compare cmp;
+        boost::make_heap(_entry_events);
+        boost::make_heap(_exit_events);
 
-        // we store the events in heaps, which is slightly more efficient
-        // than sorting them, since a heap requires linear time to construct
-        for (RecordIter i = _items.begin(); i != _items.end(); ++i) {
-            _entry_events.push_back(i);
-            _exit_events.push_back(i);
-        }
-        boost::make_heap(_entry_events, _entry_heap_compare);
-        boost::make_heap(_exit_events, _exit_heap_compare);
-        boost::pop_heap(_entry_events, _entry_heap_compare);
-        boost::pop_heap(_exit_events, _exit_heap_compare);
-
-        RecordIter next_entry = _entry_events.back();
-        RecordIter next_exit = _exit_events.back();
-        _entry_events.pop_back();
-        _exit_events.pop_back();
-
-        while (next_entry != _items.end() || next_exit != _items.end()) {
-            assert(next_exit != _items.end());
-
-            if (next_entry == _items.end() ||
-                cmp(SweepTraits::exit_value(next_exit->bound),
-                    SweepTraits::entry_value(next_entry->bound)))
-            {
+        Event next_entry = _get_next(_entry_events);
+        Event next_exit = _get_next(_exit_events);
+
+        while (next_entry || next_exit) {
+            assert(next_exit);
+
+            if (!next_entry || next_exit > next_entry) {
                 // exit event - remove record from active list
-                _leave(*next_exit);
-                _active_items.erase(_active_items.iterator_to(*next_exit));
-                if (!_exit_events.empty()) {
-                    boost::pop_heap(_exit_events, _exit_heap_compare);
-                    next_exit = _exit_events.back();
-                    _exit_events.pop_back();
-                } else {
-                    next_exit = _items.end();
-                    // we should end the loop after this happens
-                }
+                _set.removeActiveItem(next_exit.item);
+                next_exit = _get_next(_exit_events);
             } else {
                 // entry event - add record to active list
-                _enter(*next_entry);
-                _active_items.push_back(*next_entry);
-                if (!_entry_events.empty()) {
-                    boost::pop_heap(_entry_events, _entry_heap_compare);
-                    next_entry = _entry_events.back();
-                    _entry_events.pop_back();
-                } else {
-                    next_entry = _items.end();
-                }
+                _set.addActiveItem(next_entry.item);
+                next_entry = _get_next(_entry_events);
             }
         }
-
-        assert(_active_items.empty());
     }
 
-protected:
-    /// The item and its sweepline boundary.
-    struct Record {
-        boost::intrusive::list_member_hook<> _hook;
-        Bound bound;
-        Item item;
-
-        Record(Bound const &b, Item const &i)
-            : bound(b), item(i)
+private:
+    struct Event
+        : boost::totally_ordered<Event>
+    {
+        Coord coord;
+        Iter item;
+
+        Event(Coord c, Iter const &i)
+            : coord(c), item(i)
+        {}
+        Event()
+            : coord(nan("")), item()
         {}
+        bool operator<(Event const &other) const { return coord < other.coord; }
+        bool operator==(Event const &other) const { return coord == other.coord; }
+        operator bool() const { return !IS_NAN(coord); }
     };
-    typedef typename std::vector<Record>::iterator RecordIter;
-
-    typedef boost::intrusive::list
-    < Record
-    , boost::intrusive::member_hook
-        < Record
-        , boost::intrusive::list_member_hook<>
-        , &Record::_hook
-        >
-    > RecordList;
-
-    /** @brief Enter an item record.
-     * Override this to process an item as it is about to enter the active list.
-     * When called, the passed record will not be part of the active list. */
-    virtual void _enter(Record const &) {}
-    /** @brief Leave an item record.
-     * Override this to process an item as it is about to leave the active list.
-     * When called, the passed record will be part of the active list. */
-    virtual void _leave(Record const &) {}
-
-    /// The list of all item records undergoing sweeping.
-    std::vector<Record> _items;
-    /// The list of active item records.
-    RecordList _active_items;
 
-private:
-    inline static bool _entry_heap_compare(RecordIter a, RecordIter b) {
-        typename SweepTraits::Compare cmp;
-        return cmp(SweepTraits::entry_value(b->bound), SweepTraits::entry_value(a->bound));
-    }
-    inline static bool _exit_heap_compare(RecordIter a, RecordIter b) {
-        typename SweepTraits::Compare cmp;
-        return cmp(SweepTraits::exit_value(b->bound), SweepTraits::exit_value(a->bound));
+    static Event _get_next(std::vector<Event> &heap) {
+        if (heap.empty()) {
+            Event e;
+            return e;
+        }
+        boost::pop_heap(heap);
+        Event ret = heap.back();
+        heap.pop_back();
+        return ret;
     }
 
-    std::vector<RecordIter> _entry_events;
-    std::vector<RecordIter> _exit_events;
+    SweepSet &_set;
+    std::vector<Event> _entry_events;
+    std::vector<Event> _exit_events;
 };
 
 } // namespace Geom
diff --git a/src/2geom/utils.h b/src/2geom/utils.h
index bc0ad74b8..01579d81b 100644
--- a/src/2geom/utils.h
+++ b/src/2geom/utils.h
@@ -68,9 +68,33 @@ struct NullIterator
     NullIterator() {}
 
     template <typename T>
-    void operator=(T const &v) {}
+    void operator=(T const &) {}
 };
 
+/** @brief Get the next iterator in the container with wrap-around.
+ * If the iterator would become the end iterator after incrementing,
+ * return the begin iterator instead. */
+template <typename Iter, typename Container>
+Iter cyclic_next(Iter i, Container &c) {
+    ++i;
+    if (i == c.end()) {
+        i = c.begin();
+    }
+    return i;
+}
+
+/** @brief Get the previous iterator in the container with wrap-around.
+ * If the passed iterator is the begin iterator, return the iterator
+ * just before the end iterator instead. */
+template <typename Iter, typename Container>
+Iter cyclic_prior(Iter i, Container &c) {
+    if (i == c.begin()) {
+        i = c.end();
+    }
+    --i;
+    return i;
+}
+
 } // end namespace Geom
 
 #endif // LIB2GEOM_SEEN_UTILS_H
diff --git a/src/2geom/viewbox.cpp b/src/2geom/viewbox.cpp
deleted file mode 100644
index 69bd0c487..000000000
--- a/src/2geom/viewbox.cpp
+++ /dev/null
@@ -1,133 +0,0 @@
-/**
- * \file
- * \brief Convenience class for SVG viewBox handling
- *//*
- *  Authors:
- *    Krzysztof Kosiński <tweenk.pl@gmail.com>
- *
- * Copyright (C) 2013 Authors
- *
- * 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 <2geom/transforms.h>
-#include <2geom/viewbox.h>
-
-namespace Geom {
-
-/** Convert an align specification to coordinate fractions. */
-Point align_factors(Align g) {
-    Point p;
-    switch (g) {
-    case ALIGN_XMIN_YMIN:
-        p[X] = 0.0;
-        p[Y] = 0.0;
-        break;
-    case ALIGN_XMID_YMIN:
-        p[X] = 0.5;
-        p[Y] = 0.0;
-        break;
-    case ALIGN_XMAX_YMIN:
-        p[X] = 1.0;
-        p[Y] = 0.0;
-        break;
-    case ALIGN_XMIN_YMID:
-        p[X] = 0.0;
-        p[Y] = 0.5;
-        break;
-    case ALIGN_XMID_YMID:
-        p[X] = 0.5;
-        p[Y] = 0.5;
-        break;
-    case ALIGN_XMAX_YMID:
-        p[X] = 1.0;
-        p[Y] = 0.5;
-        break;
-    case ALIGN_XMIN_YMAX:
-        p[X] = 0.0;
-        p[Y] = 1.0;
-        break;
-    case ALIGN_XMID_YMAX:
-        p[X] = 0.5;
-        p[Y] = 1.0;
-        break;
-    case ALIGN_XMAX_YMAX:
-        p[X] = 1.0;
-        p[Y] = 1.0;
-        break;
-    default:
-        break;
-    }
-    return p;
-}
-
-/** Obtain transformation from the viewbox to the specified viewport. */
-Affine ViewBox::transformTo(Geom::Rect const &viewport) const
-{
-    if (!_box) {
-        return Geom::Affine::identity();
-    }
-
-    // 1. translate viewbox to origin
-    Geom::Affine total = Translate(-_box->min());
-
-    // 2. compute scale
-    Geom::Point vdims = viewport.dimensions();
-    Geom::Point bdims = _box->dimensions();
-    Geom::Scale scale(vdims[X] / bdims[X], vdims[Y] / bdims[Y]);
-
-    if (_align == ALIGN_NONE) {
-        // apply non-uniform scale
-        // = Scale(_box->dimensions()).inverse() * Scale(viewport.dimensions())
-        total *= scale * Translate(viewport.min());
-    } else {
-        double uscale = 0;
-        if (_expansion == EXPANSION_MEET) {
-            uscale = std::min(scale[X], scale[Y]);
-        } else {
-            uscale = std::max(scale[X], scale[Y]);
-        }
-        scale = Scale(uscale);
-
-        // compute offset for align
-        Geom::Point offset = bdims * scale - vdims;
-        offset *= Scale(align_factors(_align));
-        total *= Translate(-offset);
-    }
-
-    return total;
-}
-
-} // namespace Geom
-
-/*
-  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:fileencoding=utf-8:textwidth=99 :
diff --git a/src/2geom/viewbox.h b/src/2geom/viewbox.h
deleted file mode 100644
index 81f59ee36..000000000
--- a/src/2geom/viewbox.h
+++ /dev/null
@@ -1,102 +0,0 @@
-/**
- * \file
- * \brief Convenience class for SVG viewBox handling
- *//*
- *  Authors:
- *    Krzysztof Kosiński <tweenk.pl@gmail.com>
- *
- * Copyright (C) 2013 Authors
- *
- * 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.
- */
-
-#ifndef LIB2GEOM_SEEN_VIEWBOX_H
-#define LIB2GEOM_SEEN_VIEWBOX_H
-
-#include <2geom/affine.h>
-#include <2geom/rect.h>
-
-namespace Geom {
-
-/** Values for the <align> parameter of preserveAspectRatio.
- * See: http://www.w3.org/TR/SVG/coords.html#PreserveAspectRatioAttribute */
-enum Align {
-    ALIGN_NONE,
-    ALIGN_XMIN_YMIN,
-    ALIGN_XMID_YMIN,
-    ALIGN_XMAX_YMIN,
-    ALIGN_XMIN_YMID,
-    ALIGN_XMID_YMID,
-    ALIGN_XMAX_YMID,
-    ALIGN_XMIN_YMAX,
-    ALIGN_XMID_YMAX,
-    ALIGN_XMAX_YMAX
-};
-
-/** Values for the <meetOrSlice> parameter of preserveAspectRatio.
- * See: http://www.w3.org/TR/SVG/coords.html#PreserveAspectRatioAttribute */
-enum Expansion {
-    EXPANSION_MEET,
-    EXPANSION_SLICE
-};
-
-Point align_factors(Align align);
-
-class ViewBox {
-    OptRect _box;
-    Align _align;
-    Expansion _expansion;
-
-public:
-    explicit ViewBox(OptRect const &r = OptRect(), Align a = ALIGN_XMID_YMID, Expansion ex = EXPANSION_MEET)
-        : _box(r)
-        , _align(a)
-        , _expansion(ex)
-    {}
-
-    void setBox(OptRect const &r) { _box = r; }
-    void setAlign(Align a) { _align = a; }
-    void setExpansion(Expansion ex) { _expansion = ex; }
-    OptRect const &box() const { return _box; }
-    Align align() const { return _align; }
-    Expansion expansion() const { return _expansion; }
-
-    /** Obtain transformation from the viewbox to the specified viewport. */
-    Affine transformTo(Geom::Rect const &viewport) const;
-};
-
-} // namespace Geom
-
-#endif // !LIB2GEOM_SEEN_VIEWBOX_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:fileencoding=utf-8:textwidth=99 :