Upgrade to 2Geom r2413

(bzr r14059.2.18)

32 files changed, 817 insertions, 553 deletions

diff --git a/src/2geom/affine.cpp b/src/2geom/affine.cpp
index e9ca5748e..48179e864 100644
--- a/src/2geom/affine.cpp
+++ b/src/2geom/affine.cpp
@@ -6,9 +6,10 @@
  * This code is in public domain
  */
 
-#include <2geom/utils.h>
 #include <2geom/affine.h>
 #include <2geom/point.h>
+#include <2geom/polynomial.h>
+#include <2geom/utils.h>
 
 namespace Geom {
 
@@ -36,8 +37,7 @@ Point Affine::yAxis() const {
     return Point(_c[2], _c[3]);
 }
 
-/** Gets the translation imparted by the Affine.
- */
+/// Gets the translation imparted by the Affine.
 Point Affine::translation() const {
     return Point(_c[4], _c[5]);
 }
@@ -52,8 +52,7 @@ void Affine::setYAxis(Point const &vec) {
         _c[i + 2] = vec[i];
 }
 
-/** Sets the translation imparted by the Affine.
- */
+/// Sets the translation imparted by the Affine.
 void Affine::setTranslation(Point const &loc) {
     for(int i = 0; i < 2; i++)
         _c[i + 4] = loc[i];
@@ -61,33 +60,35 @@ void Affine::setTranslation(Point const &loc) {
 
 /** Calculates the amount of x-scaling imparted by the Affine.  This is the scaling applied to
  *  the original x-axis region.  It is \emph{not} the overall x-scaling of the transformation.
- *  Equivalent to L2(m.xAxis())
- */
+ *  Equivalent to L2(m.xAxis()). */
 double Affine::expansionX() const {
     return sqrt(_c[0] * _c[0] + _c[1] * _c[1]);
 }
 
 /** Calculates the amount of y-scaling imparted by the Affine.  This is the scaling applied before
  *  the other transformations.  It is \emph{not} the overall y-scaling of the transformation. 
- *  Equivalent to L2(m.yAxis())
- */
+ *  Equivalent to L2(m.yAxis()). */
 double Affine::expansionY() const {
     return sqrt(_c[2] * _c[2] + _c[3] * _c[3]);
 }
 
 void Affine::setExpansionX(double val) {
     double exp_x = expansionX();
-    if(!are_near(exp_x, 0.0)) {  //TODO: best way to deal with it is to skip op?
+    if (exp_x != 0.0) {  //TODO: best way to deal with it is to skip op?
         double coef = val / expansionX();
-        for(unsigned i=0;i<2;i++) _c[i] *= coef;
+        for (unsigned i = 0; i < 2; ++i) {
+            _c[i] *= coef;
+        }
     }
 }
 
 void Affine::setExpansionY(double val) {
     double exp_y = expansionY();
-    if(!are_near(exp_y, 0.0)) {  //TODO: best way to deal with it is to skip op?
+    if (exp_y != 0.0) {  //TODO: best way to deal with it is to skip op?
         double coef = val / expansionY();
-        for(unsigned i=2; i<4; i++) _c[i] *= coef;
+        for (unsigned i = 2; i < 4; ++i) {
+            _c[i] *= coef;
+        }
     }
 }
 
@@ -468,55 +469,38 @@ Affine elliptic_quadratic_form(Affine const &m) {
 Eigen::Eigen(Affine const &m) {
     double const B = -m[0] - m[3];
     double const C = m[0]*m[3] - m[1]*m[2];
-    double const center = -B/2.0;
-    double const delta = sqrt(B*B-4*C)/2.0;
-    values[0] = center + delta; values[1] = center - delta;
-    for (int i = 0; i < 2; i++) {
-        vectors[i] = unit_vector(rot90(Point(m[0]-values[i], m[1])));
-    }
-}
 
-static void quadratic_roots(const double q0, const double q1, const double q2, int &n, double&r0, double&r1) {
-    if(q2 == 0) {
-        if(q1 == 0) { // zero or infinite roots
-            n = 0;
-        } else {
-            n = 1;
-            r0 = -q0/q1;
-        }
-    } else {
-        double desc = q1*q1 - 4*q2*q0;
-        if (desc < 0)
-            n = 0;
-        else if (desc == 0) {
-            n = 1;
-            r0 = -q1/(2*q2);
-        } else {
-            n = 2;
-            desc = std::sqrt(desc);
-            double t = -0.5*(q1+sgn(q1)*desc);
-            r0 = t/q2;
-            r1 = q0/t;
-        }
+    std::vector<double> v = solve_quadratic(1, B, C);
+
+    for (unsigned i = 0; i < v.size(); ++i) {
+        values[i] = v[i];
+        vectors[i] = unit_vector(rot90(Point(m[0] - values[i], m[1])));
+    }
+    for (unsigned i = v.size(); i < 2; ++i) {
+        values[i] = 0;
+        vectors[i] = Point(0,0);
     }
 }
 
 Eigen::Eigen(double m[2][2]) {
     double const B = -m[0][0] - m[1][1];
     double const C = m[0][0]*m[1][1] - m[1][0]*m[0][1];
-    //double const desc = B*B-4*C;
-    //double t = -0.5*(B+sgn(B)*desc);
-    int n;
-    values[0] = values[1] = 0;
-    quadratic_roots(C, B, 1, n, values[0], values[1]);
-    for (int i = 0; i < n; i++)
-        vectors[i] = unit_vector(rot90(Point(m[0][0]-values[i], m[0][1])));
-    for (int i = n; i < 2; i++) 
+
+    std::vector<double> v = solve_quadratic(1, B, C);
+
+    for (unsigned i = 0; i < v.size(); ++i) {
+        values[i] = v[i];
+        vectors[i] = unit_vector(rot90(Point(m[0][0] - values[i], m[0][1])));
+    }
+    for (unsigned i = v.size(); i < 2; ++i) {
+        values[i] = 0;
         vectors[i] = Point(0,0);
+    }
 }
 
-/** @brief Nearness predicate for affine transforms
- * @returns True if all entries of matrices are within eps of each other */
+/** @brief Nearness predicate for affine transforms.
+ * @returns True if all entries of matrices are within eps of each other.
+ * @relates Affine */
 bool are_near(Affine const &a, Affine const &b, Coord eps)
 {
     return are_near(a[0], b[0], eps) && are_near(a[1], b[1], eps) &&
diff --git a/src/2geom/angle.h b/src/2geom/angle.h
index 9ece3b9fe..af4442e88 100644
--- a/src/2geom/angle.h
+++ b/src/2geom/angle.h
@@ -59,6 +59,9 @@ namespace Geom {
  * to <tt>double</tt> is in the range \f$[-\pi, \pi)\f$ - the convention used by C's
  * math library.
  *
+ * This class holds only a single floating point value, so passing it by value will generally
+ * be faster than passing it by const reference.
+ *
  * @ingroup Primitives
  */
 class Angle
@@ -74,12 +77,12 @@ public:
     explicit Angle(Point const &p) : _angle(atan2(p)) { _normalize(); }
     Angle(Point const &a, Point const &b) : _angle(angle_between(a, b)) { _normalize(); }
     operator Coord() const { return radians(); }
-    Angle &operator+=(Angle const &o) {
+    Angle &operator+=(Angle o) {
         _angle += o._angle;
         _normalize();
         return *this;
     }
-    Angle &operator-=(Angle const &o) {
+    Angle &operator-=(Angle o) {
         _angle -= o._angle;
         _normalize();
         return *this;
@@ -92,7 +95,7 @@ public:
         *this -= Angle(a);
         return *this;
     }
-    bool operator==(Angle const &o) const {
+    bool operator==(Angle o) const {
         return _angle == o._angle;
     }
     bool operator==(Coord c) const {
@@ -174,36 +177,112 @@ inline Angle distance(Angle const &a, Angle const &b) {
  * the end angle for 1, and interpolate linearly for other values. Note that such functions
  * are not continuous if the interval crosses the angle \f$\pi\f$.
  *
- * It is currently not possible to represent the full angle with this class.
- * If you specify the same start and end angle, the interval will be treated as empty
- * except for that value.
- *
- * This class is immutable - you cannot change the values of start and end angles
- * without creating a new instance of this class.
+ * This class can represent all directed angular intervals, including empty ones.
+ * However, not all possible intervals can be created with the constructors.
+ * For full control, use the setInitial(), setFinal() and setAngles() methods.
  *
  * @ingroup Primitives
  */
-class AngleInterval {
+class AngleInterval
+    : boost::equality_comparable< AngleInterval >
+{
 public:
-    /** @brief Create an angular interval.
+    AngleInterval() {}
+    /** @brief Create an angular interval from two angles and direction.
+     * If the initial and final angle are the same, a degenerate interval
+     * (containing only one angle) will be created.
      * @param s Starting angle
      * @param e Ending angle
      * @param cw Which direction the interval goes. True means that it goes
      *   in the direction of increasing angles, while false means in the direction
      *   of decreasing angles. */
-    AngleInterval(Angle const &s, Angle const &e, bool cw = false)
-        : _start_angle(s), _end_angle(e), _sweep(cw)
+    AngleInterval(Angle s, Angle e, bool cw = false)
+        : _start_angle(s), _end_angle(e), _sweep(cw), _full(false)
     {}
     AngleInterval(double s, double e, bool cw = false)
-        : _start_angle(s), _end_angle(e), _sweep(cw)
+        : _start_angle(s), _end_angle(e), _sweep(cw), _full(false)
     {}
+    /** @brief Create an angular interval from three angles.
+     * If the inner angle is exactly equal to initial or final angle,
+     * the sweep flag will be set to true, i.e. the interval will go
+     * in the direction of increasing angles.
+     *
+     * If the initial and final angle are the same, but the inner angle
+     * is different, a full angle in the direction of increasing angles
+     * will be created.
+     *
+     * @param s Initial angle
+     * @param inner Angle contained in the interval
+     * @param e Final angle */
+    AngleInterval(Angle s, Angle inner, Angle e)
+        : _start_angle(s)
+        , _end_angle(e)
+        , _sweep((inner-s).radians0() <= (e-s).radians0())
+        , _full(s == e && s != inner)
+    {
+        if (_full) {
+            _sweep = true;
+        }
+    }
 
     /// Get the start angle.
-    Angle const &initialAngle() const { return _start_angle; }
+    Angle initialAngle() const { return _start_angle; }
     /// Get the end angle.
-    Angle const &finalAngle() const { return _end_angle; }
+    Angle finalAngle() const { return _end_angle; }
+    /// Check whether the interval goes in the direction of increasing angles.
+    bool sweep() const { return _sweep; }
     /// Check whether the interval contains only a single angle.
-    bool isDegenerate() const { return initialAngle() == finalAngle(); }
+    bool isDegenerate() const {
+        return _start_angle == _end_angle && !_full;
+    }
+    /// Check whether the interval contains all angles.
+    bool isFull() const {
+        return _start_angle == _end_angle && _full;
+    }
+
+    /** @brief Set the initial angle.
+     * @param a Angle to set
+     * @param prefer_full Whether to set a full angular interval when
+     * the initial angle is set to the final angle */
+    void setInitial(Angle a, bool prefer_full = false) {
+        _start_angle = a;
+        _full = prefer_full && a == _end_angle;
+    }
+
+    /** @brief Set the final angle.
+     * @param a Angle to set
+     * @param prefer_full Whether to set a full angular interval when
+     * the initial angle is set to the final angle */
+    void setFinal(Angle a, bool prefer_full = false) {
+        _end_angle = a;
+        _full = prefer_full && a == _start_angle;
+    }
+    /** @brief Set both angles at once.
+     * The direction (sweep flag) is left unchanged.
+     * @param s Initial angle
+     * @param e Final angle
+     * @param prefer_full Whether to set a full interval when the passed
+     * initial and final angle are the same */
+    void setAngles(Angle s, Angle e, bool prefer_full = false) {
+        _start_angle = s;
+        _end_angle = e;
+        _full = prefer_full && s == e;
+    }
+    /// Set whether the interval goes in the direction of increasing angles.
+    void setSweep(bool s) { _sweep = s; }
+
+    /// Reverse the direction of the interval while keeping contained values the same.
+    void reverse() {
+        using std::swap;
+        swap(_start_angle, _end_angle);
+        _sweep = !_sweep;
+    }
+    /// Get a new interval with reversed direction.
+    AngleInterval reversed() const {
+        AngleInterval result(*this);
+        result.reverse();
+        return result;
+    }
 
     /// Get an angle corresponding to the specified time value.
     Angle angleAt(Coord t) const {
@@ -214,8 +293,14 @@ public:
     Angle operator()(Coord t) const { return angleAt(t); }
 
     /** @brief Compute a time value that would evaluate to the given angle.
-     * If the start and end angle are exactly the same, NaN will be returned. */
-    Coord timeAtAngle(Angle const &a) const {
+     * If the start and end angle are exactly the same, NaN will be returned.
+     * Negative values will be returned for angles between the initial angle
+     * and the angle exactly opposite the midpoint of the interval. */
+    Coord timeAtAngle(Angle a) const {
+        if (_full) {
+            Angle ta = _sweep ? a - _start_angle : _start_angle - a;
+            return ta.radians0() / (2*M_PI);
+        }
         Coord ex = extent();
         Coord outex = 2*M_PI - ex;
         if (_sweep) {
@@ -237,8 +322,9 @@ public:
         }
     }
 
-    /** @brief Check whether the interval includes the given angle. */
-    bool contains(Angle const &a) const {
+    /// Check whether the interval includes the given angle.
+    bool contains(Angle a) const {
+        if (_full) return true;
         Coord s = _start_angle.radians0();
         Coord e = _end_angle.radians0();
         Coord x = a.radians0();
@@ -254,6 +340,7 @@ public:
      * Equivalent to the absolute value of the sweep angle.
      * @return Extent in range \f$[0, 2\pi)\f$. */
     Coord extent() const {
+        if (_full) return 2*M_PI;
         return _sweep
             ? (_end_angle - _start_angle).radians0()
             : (_start_angle - _end_angle).radians0();
@@ -263,16 +350,35 @@ public:
      * It is positive when sweep is true. Denoted as \f$\Delta\theta\f$ in the SVG
      * elliptical arc implementation notes. */
     Coord sweepAngle() const {
+        if (_full) return _sweep ? 2*M_PI : -2*M_PI;
         Coord sa = _end_angle.radians0() - _start_angle.radians0();
         if (_sweep && sa < 0) sa += 2*M_PI;
         if (!_sweep && sa > 0) sa -= 2*M_PI;
         return sa;
     }
-protected:
-    AngleInterval() {}
+
+    /// Check another interval for equality.
+    bool operator==(AngleInterval const &other) const {
+        if (_start_angle != other._start_angle) return false;
+        if (_end_angle != other._end_angle) return false;
+        if (_sweep != other._sweep) return false;
+        if (_full != other._full) return false;
+        return true;
+    }
+
+    static AngleInterval create_full(Angle start, bool sweep = true) {
+        AngleInterval result;
+        result._start_angle = result._end_angle = start;
+        result._sweep = sweep;
+        result._full = true;
+        return result;
+    }
+
+private:
     Angle _start_angle;
     Angle _end_angle;
     bool _sweep;
+    bool _full;
 };
 
 /** @brief Given an angle in degrees, return radians
@@ -282,86 +388,6 @@ inline Coord deg_to_rad(Coord deg) { return deg*M_PI/180.0;}
  * @relates Angle */
 inline Coord rad_to_deg(Coord rad) { return rad*180.0/M_PI;}
 
-/*
- *  start_angle and angle must belong to [0, 2PI[
- *  and angle must belong to the cirsular arc defined by
- *  start_angle, end_angle and with rotation direction cw
- */
-inline
-double map_circular_arc_on_unit_interval( double angle, double start_angle, double end_angle, bool cw = true )
-{
-    double d = end_angle - start_angle;
-    double t = angle - start_angle;
-    if ( !cw )
-    {
-    	d = -d;
-    	t = -t;
-    }
-    d = std::fmod(d, 2*M_PI);
-    t = std::fmod(t, 2*M_PI);
-    if ( d < 0 ) d += 2*M_PI;
-    if ( t < 0 ) t += 2*M_PI;
-    return t / d;
-}
-
-inline
-Coord map_unit_interval_on_circular_arc(Coord t, double start_angle, double end_angle, bool cw = true)
-{
-	double sweep_angle = end_angle - start_angle;
-	if ( !cw ) sweep_angle = -sweep_angle;
-	sweep_angle = std::fmod(sweep_angle, 2*M_PI);
-	if ( sweep_angle < 0 ) sweep_angle += 2*M_PI;
-
-	Coord angle = start_angle;
-    if ( cw )
-    {
-        angle += sweep_angle * t;
-    }
-    else
-    {
-        angle -= sweep_angle * t;
-    }
-    angle = std::fmod(angle, 2*M_PI);
-    if (angle < 0) angle += 2*M_PI;
-    return angle;
-}
-
-/*
- *  Return true if the given angle is contained in the circular arc determined
- *  by the passed angles.
- *
- *  a:     the angle to be tested
- *  sa:    the angle the arc start from
- *  ia:    an angle strictly inner to the arc
- *  ea:    the angle the arc end to
- *
- *  prerequisite: the inner angle has to be not equal (mod 2PI) to the start
- *                or the end angle, except when they are equal each other, too.
- *  warning: when sa == ea (mod 2PI) they define a whole circle
- *           if ia != sa (mod 2PI), on the contrary if ia == sa == ea (mod 2PI)
- *           they define a single point.
- */
-inline
-bool arc_contains (double a, double sa, double ia, double ea)
-{
-    a -= sa;
-    a = std::fmod(a, 2*M_PI);
-    if (a < 0) a += 2*M_PI;
-    ia -= sa;
-    ia = std::fmod(ia, 2*M_PI);
-    if (ia < 0) ia += 2*M_PI;
-    ea -= sa;
-    ea = std::fmod(ea, 2*M_PI);
-    if (ea < 0) ea += 2*M_PI;
-
-    if (ia == 0 && ea == 0)  return (a == 0);
-    if (ia == 0 || ia == ea)
-    {
-        THROW_RANGEERROR ("arc_contains: passed angles do not define an arc");
-    }
-    return (ia < ea && a <= ea) || (ia > ea && (a >= ea || a == 0));
-}
-
 } // end namespace Geom
 
 namespace std {
diff --git a/src/2geom/bezier-curve.h b/src/2geom/bezier-curve.h
index 636a263a7..9ac4d7b4d 100644
--- a/src/2geom/bezier-curve.h
+++ b/src/2geom/bezier-curve.h
@@ -129,8 +129,21 @@ public:
         return new BezierCurve(Geom::reverse(inner));
     }
 
-    virtual void transform(Affine const &m) {
-        for (unsigned i = 0;  i < size(); ++i) {
+    using Curve::operator*=;
+    virtual void operator*=(Translate const &tr) {
+        for (unsigned i = 0; i < size(); ++i) {
+            inner[X][i] += tr[X];
+            inner[Y][i] += tr[Y];
+        }
+    }
+    virtual void operator*=(Scale const &s) {
+        for (unsigned i = 0; i < size(); ++i) {
+            inner[X][i] *= s[X];
+            inner[Y][i] *= s[Y];
+        }
+    }
+    virtual void operator*=(Affine const &m) {
+        for (unsigned i = 0; i < size(); ++i) {
             setPoint(i, controlPoint(i) * m);
         }
     }
diff --git a/src/2geom/circle.cpp b/src/2geom/circle.cpp
index ec59bbe4a..553981a72 100644
--- a/src/2geom/circle.cpp
+++ b/src/2geom/circle.cpp
@@ -101,6 +101,12 @@ Zoom Circle::inverseUnitCircleTransform() const
     return ret;
 }
 
+Point Circle::initialPoint() const
+{
+    Point p(_center);
+    p[X] += _radius;
+    return p;
+}
 
 Point Circle::pointAt(Coord t) const {
     return _center + Point::polar(t) * _radius;
@@ -112,11 +118,11 @@ Coord Circle::valueAt(Coord t, Dim2 d) const {
 }
 
 Coord Circle::timeAt(Point const &p) const {
+    if (_center == p) return 0;
     return atan2(p - _center);
 }
 
 Coord Circle::nearestTime(Point const &p) const {
-    if (_center == p) return 0;
     return timeAt(p);
 }
 
diff --git a/src/2geom/circle.h b/src/2geom/circle.h
index c42cb7f80..a4d5f2097 100644
--- a/src/2geom/circle.h
+++ b/src/2geom/circle.h
@@ -76,6 +76,7 @@ public:
     Coord center(Dim2 d) const { return _center[d]; }
     Coord radius() const { return _radius; }
     Coord area() const { return M_PI * _radius * _radius; }
+    bool isDegenerate() const { return _radius == 0; }
 
     void setCenter(Point const &p) { _center = p; }
     void setRadius(Coord c) { _radius = c; }
@@ -83,6 +84,8 @@ public:
     Rect boundsFast() const;
     Rect boundsExact() const { return boundsFast(); }
 
+    Point initialPoint() const;
+    Point finalPoint() const { return initialPoint(); }
     Point pointAt(Coord t) const;
     Coord valueAt(Coord t, Dim2 d) const;
     Coord timeAt(Point const &p) const;
@@ -138,7 +141,13 @@ bool are_near(Circle const &a, Circle const &b, Coord eps=EPSILON);
 
 std::ostream &operator<<(std::ostream &out, Circle const &c);
 
-//bool are_near(Circle const &a, Circle const &b, Coord eps = EPSILON);
+template <>
+struct ShapeTraits<Circle> {
+    typedef Coord TimeType;
+    typedef Interval IntervalType;
+    typedef Ellipse AffineClosureType;
+    typedef Intersection<> IntersectionType;
+};
 
 } // end namespace Geom
 
diff --git a/src/2geom/concepts.h b/src/2geom/concepts.h
index 7bc1bc0fd..c331f078b 100644
--- a/src/2geom/concepts.h
+++ b/src/2geom/concepts.h
@@ -39,6 +39,7 @@
 #include <vector>
 #include <boost/concept/assert.hpp>
 #include <2geom/forward.h>
+#include <2geom/transforms.h>
 
 namespace Geom {
 
@@ -82,8 +83,8 @@ struct FragmentConcept {
         o = t.valueAt(d);
         o = t(d);
         v = t.valueAndDerivatives(d, u-1);
-		//Is a pure derivative (ignoring others) accessor ever much faster?
-		//u = number of values returned. first val is value.
+        //Is a pure derivative (ignoring others) accessor ever much faster?
+        //u = number of values returned. first val is value.
         sb = t.toSBasis();
         t = reverse(t);
         i = bounds_fast(t);
@@ -104,26 +105,30 @@ struct ShapeConcept {
     typedef typename ShapeTraits<T>::TimeType Time;
     typedef typename ShapeTraits<T>::IntervalType Interval;
     typedef typename ShapeTraits<T>::AffineClosureType AffineClosure;
-    //typedef typename ShapeTraits<T>::IntersectionType Isect;
+    typedef typename ShapeTraits<T>::IntersectionType Isect;
 
-    T shape;
+    T shape, other;
     Time t;
     Point p;
     AffineClosure ac;
     Affine m;
+    Translate tr;
     Coord c;
     bool bool_;
+    std::vector<Isect> ivec;
 
     void constraints() {
         p = shape.pointAt(t);
         c = shape.valueAt(t, X);
-        //ivec = shape.intersect(other);
+        ivec = shape.intersect(other);
         t = shape.nearestTime(p);
+        shape *= tr;
         ac = shape;
         ac *= m;
         bool_ = (shape == shape);
-        bool_ = (shape != shape);
+        bool_ = (shape != other);
         bool_ = shape.isDegenerate();
+        //bool_ = are_near(shape, other, c);
     }
 };
 
diff --git a/src/2geom/conicsec.cpp b/src/2geom/conicsec.cpp
index 889797de3..3b36137be 100644
--- a/src/2geom/conicsec.cpp
+++ b/src/2geom/conicsec.cpp
@@ -80,41 +80,6 @@ static double boxprod(Point a, Point b, Point c) {
     return det(a,b) - det(a,c) + det(b,c);
 }
 
-
-/**
- * Find the roots of (q2x + q1)x+q0 = 0
- * Tries to be numerically robust.
- */
-template <typename T>
-static std::vector<T> quadratic_roots(T q0, T q1, T q2) {
-    std::vector<double> r;
-    if(q2 == 0) {
-        if(q1 == 0) { // zero or infinite roots
-            return r;
-        }
-        r.push_back(-q0/q1);
-    } else {
-        double desc = q1*q1 - 4*q2*q0;
-        /*cout << q2 << ", "
-          << q1 << ", "
-          << q0 << "; "
-          << desc << "\n";*/
-        if (desc < 0)
-            return r;
-        else if (desc == 0)
-            r.push_back(-q1/(2*q2));
-        else {
-            desc = std::sqrt(desc);
-            double t = -0.5*(q1+sgn(q1)*desc);
-            r.push_back(t/q2);
-            r.push_back(q0/t);
-        }
-    }
-    return r;
-}
-
-
-
 class BadConversion : public std::runtime_error {
 public:
     BadConversion(const std::string& s)
diff --git a/src/2geom/conicsec.h b/src/2geom/conicsec.h
index e9c466978..dbe564872 100644
--- a/src/2geom/conicsec.h
+++ b/src/2geom/conicsec.h
@@ -462,13 +462,8 @@ public:
     bool arc_contains (const Point & _point, const Point & _initial,
                        const Point & _inner, const Point & _final) const
     {
-        double pa = angle_at (_point);
-        double sa = angle_at (_initial);
-        double ia = angle_at (_inner);
-        double ea = angle_at (_final);
-        // we test if _point and _inner have the same position
-        // wrt _initial and _final
-        return Geom::arc_contains (pa, sa, ia, ea);
+        AngleInterval ai(angle_at(_initial), angle_at(_inner), angle_at(_final));
+        return ai.contains(angle_at(_point));
     }
 
     Rect arc_bound (const Point & P1, const Point & Q, const Point & P2) const;
diff --git a/src/2geom/coord.cpp b/src/2geom/coord.cpp
index 9ee8066f2..8b5e28586 100644
--- a/src/2geom/coord.cpp
+++ b/src/2geom/coord.cpp
@@ -3656,7 +3656,7 @@ std::string format_coord_nice(Coord x)
 {
     static DoubleToStringConverter conv(
         DoubleToStringConverter::UNIQUE_ZERO,
-        "Inf", "NaN", 'e', -6, 21, 0, 0);
+        "inf", "NaN", 'e', -6, 21, 0, 0);
     std::string ret;
     ret.reserve(32);
     conv.ToShortest(x, ret);
@@ -3669,7 +3669,7 @@ Coord parse_coord(std::string const &s)
         StringToDoubleConverter::ALLOW_LEADING_SPACES |
         StringToDoubleConverter::ALLOW_TRAILING_SPACES |
         StringToDoubleConverter::ALLOW_SPACES_AFTER_SIGN,
-        0.0, nan(""), "Inf", "NaN");
+        0.0, nan(""), "inf", "NaN");
     int dummy;
     return conv.StringToDouble(s.c_str(), s.length(), &dummy);
 }
diff --git a/src/2geom/coord.h b/src/2geom/coord.h
index 416e6443d..9cc220db7 100644
--- a/src/2geom/coord.h
+++ b/src/2geom/coord.h
@@ -1,7 +1,10 @@
 /** @file
  * @brief Integral and real coordinate types and some basic utilities
  *//*
- * Copyright 2006 Nathan Hurst <njh@mail.csse.monash.edu.au>
+ * Authors:
+ *   Nathan Hurst <njh@mail.csse.monash.edu.au>
+ *   Krzysztof Kosiński <tweenk.pl@gmail.com>
+ * Copyright 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
@@ -40,10 +43,12 @@
 
 namespace Geom {
 
-/// 2D axis enumeration (X or Y).
+/** @brief 2D axis enumeration (X or Y).
+ * @ingroup Primitives */
 enum Dim2 { X=0, Y=1 };
 
-/// Get the other (perpendicular) dimension.
+/** @brief Get the other (perpendicular) dimension.
+ * @ingroup Primitives */
 inline Dim2 other_dimension(Dim2 d) { return d == Y ? X : Y; }
 
 // TODO: make a smarter implementation with C++11
@@ -54,47 +59,48 @@ struct D2Traits {
     typedef typename T::D1ConstReference D1ConstReference;
 };
 
-// for use with things such as transform_iterator
-template <typename T>
-struct GetX {
-    typedef typename D2Traits<T>::D1Value result_type;
-    typedef T argument_type;
-    typename D2Traits<T>::D1Value operator()(T const &a) const {
-        return a[X];
-    }
-};
-template <typename T>
-struct GetY {
+/** @brief Axis extraction functor.
+ * For use with things such as Boost's transform_iterator.
+ * @ingroup Utilities */
+template <Dim2 D, typename T>
+struct GetAxis {
     typedef typename D2Traits<T>::D1Value result_type;
     typedef T argument_type;
     typename D2Traits<T>::D1Value operator()(T const &a) const {
-        return a[Y];
+        return a[D];
     }
 };
 
-/**
- * @brief Floating point type used to store coordinates.
- *
- * You may safely assume that double (or even float) provides enough precision for storing
- * on-canvas points, and hence that double provides enough precision for dot products of
- * differences of on-canvas points.
- */
+/** @brief Floating point type used to store coordinates.
+ * @ingroup Primitives */
 typedef double Coord;
+
+/** @brief Type used for integral coordinates.
+ * @ingroup Primitives */
 typedef int IntCoord;
 
-const Coord EPSILON = 1e-5; //1e-18;
+/** @brief Default "acceptably small" value.
+ * @ingroup Primitives */
+const Coord EPSILON = 1e-6; //1e-18;
 
+/** @brief Get a value representing infinity.
+ * @ingroup Primitives */
 inline Coord infinity() {  return std::numeric_limits<Coord>::infinity();  }
 
-//IMPL: NearConcept
+/** @brief Nearness predicate for values.
+ * @ingroup Primitives */
 inline bool are_near(Coord a, Coord b, double eps=EPSILON) { return a-b <= eps && a-b >= -eps; }
 inline bool rel_error_bound(Coord a, Coord b, double eps=EPSILON) { return a <= eps*b && a >= -eps*b; }
 
-/// Numerically stable linear interpolation.
+/** @brief Numerically stable linear interpolation.
+ * @ingroup Primitives */
 inline Coord lerp(Coord t, Coord a, Coord b) {
     return (1 - t) * a + t * b;
 }
 
+/** @brief Traits class used with coordinate types.
+ * Defines point, interval and rectangle types for the given coordinate type.
+ * @ingroup Utilities */
 template <typename C>
 struct CoordTraits {
     typedef D2<C> PointType;
@@ -174,19 +180,22 @@ struct CoordTraits<Coord> {
         RectOps;
 };
 
-/** @brief Convert coordinate to shortest possible string
- * @return The shortest string that parses back to the original value. */
+/** @brief Convert coordinate to shortest possible string.
+ * @return The shortest string that parses back to the original value.
+ * @relates Coord */
 std::string format_coord_shortest(Coord x);
 
-/** @brief Convert coordinate to human-readable string
+/** @brief Convert coordinate to human-readable string.
  * Unlike format_coord_shortest, this function will not omit a leading zero
  * before a decimal point or use small negative exponents. The output format
- * is similar to Javascript functions. */
+ * is similar to Javascript functions.
+ * @relates Coord */
 std::string format_coord_nice(Coord x);
 
-/** @brief Parse coordinate string
+/** @brief Parse coordinate string.
  * When using this function in conjunction with format_coord_shortest()
- * or format_coord_nice(), the value is guaranteed to be preserved exactly. */
+ * or format_coord_nice(), the value is guaranteed to be preserved exactly.
+ * @relates Coord */
 Coord parse_coord(std::string const &s);
 
 } // end namespace Geom
diff --git a/src/2geom/curve.h b/src/2geom/curve.h
index 7da0d17a0..abbdb1100 100644
--- a/src/2geom/curve.h
+++ b/src/2geom/curve.h
@@ -184,7 +184,17 @@ public:
      * Because of this method, all curve types must be closed under affine
      * transformations.
      * @param m Affine describing the affine transformation */
-    virtual void transform(Affine const &m) = 0;
+    void transform(Affine const &m) {
+         *this *= m;
+    }
+
+    virtual void operator*=(Translate const &tr) { *this *= Affine(tr); }
+    virtual void operator*=(Scale const &s) { *this *= Affine(s); }
+    virtual void operator*=(Rotate const &r) { *this *= Affine(r); }
+    virtual void operator*=(HShear const &hs) { *this *= Affine(hs); }
+    virtual void operator*=(VShear const &vs) { *this *= Affine(vs); }
+    virtual void operator*=(Zoom const &z) { *this *= Affine(z); }
+    virtual void operator*=(Affine const &m) = 0;
 
     /** @brief Create a curve transformed by an affine transformation.
      * This method returns a new curve instead modifying the existing one.
diff --git a/src/2geom/d2.h b/src/2geom/d2.h
index bf764539d..dca6fa614 100644
--- a/src/2geom/d2.h
+++ b/src/2geom/d2.h
@@ -43,20 +43,16 @@
 
 namespace Geom {
 /**
- * The D2 class takes two instances of a scalar data type and treats them
- * like a point. All operations which make sense on a point are defined for D2.
- * A D2<double> is a Point. A D2<Interval> is a standard axis aligned rectangle.
- * D2<SBasis> provides a 2d parametric function which maps t to a point
- * x(t), y(t)
+ * @brief Adaptor that creates 2D functions from 1D ones.
+ * @ingroup Fragments
  */
-template <class T>
-class D2{
-    //BOOST_CLASS_REQUIRE(T, boost, AssignableConcept);
-  private:
+template <typename T>
+class D2
+{
+private:
     T f[2];
 
-  public:
-
+public:
     typedef T D1Value;
     typedef T &D1Reference;
     typedef T const &D1ConstReference;
@@ -74,24 +70,24 @@ class D2{
     template <typename Iter>
     D2(Iter first, Iter last) {
         typedef typename std::iterator_traits<Iter>::value_type V;
-        typedef typename boost::transform_iterator<GetX<V>, Iter> XIter;
-        typedef typename boost::transform_iterator<GetY<V>, Iter> YIter;
+        typedef typename boost::transform_iterator<GetAxis<X,V>, Iter> XIter;
+        typedef typename boost::transform_iterator<GetAxis<Y,V>, Iter> YIter;
 
-        XIter xfirst(first, GetX<V>()), xlast(last, GetX<V>());
+        XIter xfirst(first, GetAxis<X,V>()), xlast(last, GetAxis<X,V>());
         f[X] = T(xfirst, xlast);
-        YIter yfirst(first, GetY<V>()), ylast(last, GetY<V>());
+        YIter yfirst(first, GetAxis<Y,V>()), ylast(last, GetAxis<Y,V>());
         f[Y] = T(yfirst, ylast);
     }
 
     D2(std::vector<Point> const &vec) {
         typedef Point V;
         typedef std::vector<Point>::const_iterator Iter;
-        typedef boost::transform_iterator<GetX<V>, Iter> XIter;
-        typedef boost::transform_iterator<GetY<V>, Iter> YIter;
+        typedef boost::transform_iterator<GetAxis<X,V>, Iter> XIter;
+        typedef boost::transform_iterator<GetAxis<Y,V>, Iter> YIter;
 
-        XIter xfirst(vec.begin(), GetX<V>()), xlast(vec.end(), GetX<V>());
+        XIter xfirst(vec.begin(), GetAxis<X,V>()), xlast(vec.end(), GetAxis<X,V>());
         f[X] = T(xfirst, xlast);
-        YIter yfirst(vec.begin(), GetY<V>()), ylast(vec.end(), GetY<V>());
+        YIter yfirst(vec.begin(), GetAxis<Y,V>()), ylast(vec.end(), GetAxis<Y,V>());
         f[Y] = T(yfirst, ylast);
     }
 
diff --git a/src/2geom/ellipse.cpp b/src/2geom/ellipse.cpp
index 46c60d85d..0264ab4ae 100644
--- a/src/2geom/ellipse.cpp
+++ b/src/2geom/ellipse.cpp
@@ -97,6 +97,14 @@ void Ellipse::setCoefficients(double A, double B, double C, double D, double E,
     makeCanonical();
 }
 
+Point Ellipse::initialPoint() const
+{
+    Coord sinrot, cosrot;
+    sincos(_angle, sinrot, cosrot);
+    Point p(ray(X) * cosrot + center(X), ray(X) * sinrot + center(Y));
+    return p;
+}
+
 
 Affine Ellipse::unitCircleTransform() const
 {
@@ -299,8 +307,19 @@ Point Ellipse::pointAt(Coord t) const
 
 Coord Ellipse::valueAt(Coord t, Dim2 d) const
 {
-    // TODO: more efficient version.
-    return pointAt(t)[d];
+    Coord sinrot, cosrot, cost, sint;
+    sincos(rotationAngle(), sinrot, cosrot);
+    sincos(t, sint, cost);
+
+    if ( d == X ) {
+        return    ray(X) * cosrot * cost
+                - ray(Y) * sinrot * sint
+                + center(X);
+    } else {
+        return    ray(X) * sinrot * cost
+                + ray(Y) * cosrot * sint
+                + center(Y);
+    }
 }
 
 Coord Ellipse::timeAt(Point const &p) const
@@ -319,6 +338,20 @@ Coord Ellipse::timeAt(Point const &p) const
     return Angle(atan2(p * iuct)).radians0(); // return a value in [0, 2pi)
 }
 
+Point Ellipse::unitTangentAt(Coord t) const
+{
+    Point p = Point::polar(t + M_PI/2);
+    p *= unitCircleTransform().withoutTranslation();
+    p.normalize();
+    return p;
+}
+
+bool Ellipse::contains(Point const &p) const
+{
+    Point tp = p * inverseUnitCircleTransform();
+    return tp.length() <= 1;
+}
+
 std::vector<ShapeIntersection> Ellipse::intersect(Line const &line) const
 {
 
diff --git a/src/2geom/ellipse.h b/src/2geom/ellipse.h
index 6fb5ed254..4b1701cce 100644
--- a/src/2geom/ellipse.h
+++ b/src/2geom/ellipse.h
@@ -124,6 +124,10 @@ public:
     Coord ray(Dim2 d) const { return _rays[d]; }
     /// Get the angle the X ray makes with the +X axis.
     Angle rotationAngle() const { return _angle; }
+    /// Get the point corresponding to the +X ray of the ellipse.
+    Point initialPoint() const;
+    /// Get the point corresponding to the +X ray of the ellipse.
+    Point finalPoint() const { return initialPoint(); }
 
     /** @brief Create an ellipse passing through the specified points
      * At least five points have to be specified. */
@@ -184,6 +188,12 @@ public:
      * with the ellipse. Note that this is NOT the nearest point on the ellipse. */
     Coord timeAt(Point const &p) const;
 
+    /// Get the value of the derivative at time t normalized to unit length.
+    Point unitTangentAt(Coord t) const;
+
+    /// Check whether the ellipse contains the given point.
+    bool contains(Point const &p) const;
+
     /// Compute intersections with an infinite line.
     std::vector<ShapeIntersection> intersect(Line const &line) const;
     /// Compute intersections with a line segment.
diff --git a/src/2geom/elliptical-arc.cpp b/src/2geom/elliptical-arc.cpp
index a04b730b3..25a78b4ad 100644
--- a/src/2geom/elliptical-arc.cpp
+++ b/src/2geom/elliptical-arc.cpp
@@ -83,9 +83,11 @@ namespace Geom
  * and the ellipse's rotation. Each set of those parameters corresponds to four different arcs,
  * with two of them larger than half an ellipse and two of them turning clockwise while traveling
  * from initial to final point. The two flags disambiguate between them: "large arc flag" selects
- * the bigger arc, while the "sweep flag" selects the clockwise arc.
+ * the bigger arc, while the "sweep flag" selects the arc going in the direction of positive
+ * angles. Angles always increase when going from the +X axis in the direction of the +Y axis,
+ * so if Y grows downwards, this means clockwise.
  *
- * @image html elliptical-arc-flags.png "The four possible arcs and the meaning of flags"
+ * @image html elliptical-arc-flags.png "Meaning of arc flags (Y grows downwards)"
  *
  * @ingroup Curves
  */
@@ -142,31 +144,7 @@ Point EllipticalArc::pointAtAngle(Coord t) const
 
 Coord EllipticalArc::valueAtAngle(Coord t, Dim2 d) const
 {
-    Coord sinrot, cosrot, cost, sint;
-    sincos(rotationAngle(), sinrot, cosrot);
-    sincos(t, sint, cost);
-
-    if ( d == X ) {
-        return    ray(X) * cosrot * cost
-                - ray(Y) * sinrot * sint
-                + center(X);
-    } else {
-        return    ray(X) * sinrot * cost
-                + ray(Y) * cosrot * sint
-                + center(Y);
-    }
-}
-
-Affine EllipticalArc::unitCircleTransform() const
-{
-    Affine ret = _ellipse.unitCircleTransform();
-    return ret;
-}
-
-Affine EllipticalArc::inverseUnitCircleTransform() const
-{
-    Affine ret = _ellipse.inverseUnitCircleTransform();
-    return ret;
+    return _ellipse.valueAt(t, d);
 }
 
 std::vector<Coord> EllipticalArc::roots(Coord v, Dim2 d) const
@@ -176,6 +154,7 @@ std::vector<Coord> EllipticalArc::roots(Coord v, Dim2 d) const
     }
 
     std::vector<Coord> sol;
+    Interval unit_interval(0, 1);
 
     if ( are_near(ray(X), 0) && are_near(ray(Y), 0) ) {
         if ( center(d) == v )
@@ -251,7 +230,7 @@ std::vector<Coord> EllipticalArc::roots(Coord v, Dim2 d) const
             {
                 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 ( logical_xor( sweep(), are_near(initialAngle(), M_PI/2) )  )
                     {
                         if ( s < 0 ) s += 2*M_PI;
                     }
@@ -263,7 +242,7 @@ std::vector<Coord> EllipticalArc::roots(Coord v, Dim2 d) const
                     break;
                 case Y:
                     s = std::acos(s); // return a value in [0,PI]
-                    if ( logical_xor( _sweep, are_near(initialAngle(), 0) ) )
+                    if ( logical_xor( sweep(), are_near(initialAngle(), 0) ) )
                     {
                         s = 2*M_PI - s;
                         if ( !(s < 2*M_PI) ) s -= 2*M_PI;
@@ -272,10 +251,11 @@ std::vector<Coord> EllipticalArc::roots(Coord v, Dim2 d) const
             }
 
             //std::cerr << "s = " << rad_to_deg(s);
-            s = map_to_01(s);
+            s = timeAtAngle(s);
             //std::cerr << " -> t: " << s << std::endl;
-            if ( !(s < 0 || s > 1) )
+            if (unit_interval.contains(s)) {
                 sol.push_back(s);
+            }
             return sol;
         }
     }
@@ -331,13 +311,13 @@ std::vector<Coord> EllipticalArc::roots(Coord v, Dim2 d) const
     }
 
     std::vector<double> arc_sol;
-    for (unsigned int i = 0; i < sol.size(); ++i )
-    {
+    for (unsigned int i = 0; i < sol.size(); ++i ) {
         //std::cerr << "s = " << rad_to_deg(sol[i]);
-        sol[i] = map_to_01(sol[i]);
+        sol[i] = timeAtAngle(sol[i]);
         //std::cerr << " -> t: " << sol[i] << std::endl;
-        if ( !(sol[i] < 0 || sol[i] > 1) )
+        if (unit_interval.contains(sol[i])) {
             arc_sol.push_back(sol[i]);
+        }
     }
     return arc_sol;
 }
@@ -355,16 +335,8 @@ Curve *EllipticalArc::derivative() const
 
     EllipticalArc *result = static_cast<EllipticalArc*>(duplicate());
     result->_ellipse.setCenter(0, 0);
-    result->_start_angle += M_PI/2;
-    if( !( result->_start_angle < 2*M_PI ) )
-    {
-        result->_start_angle -= 2*M_PI;
-    }
-    result->_end_angle += M_PI/2;
-    if( !( result->_end_angle < 2*M_PI ) )
-    {
-        result->_end_angle -= 2*M_PI;
-    }
+    result->_angles.setInitial(result->_angles.initialAngle() + M_PI/2);
+    result->_angles.setFinal(result->_angles.finalAngle() + M_PI/2);
     result->_initial_point = result->pointAtAngle( result->initialAngle() );
     result->_final_point = result->pointAtAngle( result->finalAngle() );
     return result;
@@ -381,15 +353,14 @@ EllipticalArc::pointAndDerivatives(Coord t, unsigned int n) const
     unsigned int nn = n+1; // nn represents the size of the result vector.
     std::vector<Point> result;
     result.reserve(nn);
-    double angle = map_unit_interval_on_circular_arc(t, initialAngle(),
-                                                     finalAngle(), _sweep);
+    double angle = angleAt(t);
     std::auto_ptr<EllipticalArc> ea( static_cast<EllipticalArc*>(duplicate()) );
     ea->_ellipse.setCenter(0, 0);
     unsigned int m = std::min(nn, 4u);
     for ( unsigned int i = 0; i < m; ++i )
     {
         result.push_back( ea->pointAtAngle(angle) );
-        angle += (_sweep ? M_PI/2 : -M_PI/2);
+        angle += (sweep() ? M_PI/2 : -M_PI/2);
         if ( !(angle < 2*M_PI) ) angle -= 2*M_PI;
     }
     m = nn / 4;
@@ -408,17 +379,18 @@ EllipticalArc::pointAndDerivatives(Coord t, unsigned int n) const
     return result;
 }
 
-bool EllipticalArc::containsAngle(Coord angle) const
-{
-    return AngleInterval::contains(angle);
-}
-
 Point EllipticalArc::pointAt(Coord t) const
 {
     if (isChord()) return chord().pointAt(t);
     return _ellipse.pointAt(angleAt(t));
 }
 
+Coord EllipticalArc::valueAt(Coord t, Dim2 d) const
+{
+    if (isChord()) return chord().valueAt(t, d);
+    return valueAtAngle(angleAt(t), d);
+}
+
 Curve* EllipticalArc::portion(double f, double t) const
 {
     // fix input arguments
@@ -427,39 +399,30 @@ Curve* EllipticalArc::portion(double f, double t) const
     if (t < 0) t = 0;
     if (t > 1) t = 1;
 
-    if (f == t)
-    {
-        EllipticalArc *arc = static_cast<EllipticalArc*>(duplicate());
+    if (f == t) {
+        EllipticalArc *arc = new EllipticalArc();
         arc->_initial_point = arc->_final_point = pointAt(f);
-        arc->_start_angle = arc->_end_angle = 0;
-        arc->_ellipse = Ellipse();
-        arc->_sweep = false;
-        arc->_large_arc = false;
         return arc;
     }
 
     EllipticalArc *arc = static_cast<EllipticalArc*>(duplicate());
     arc->_initial_point = pointAt(f);
     arc->_final_point = pointAt(t);
-    if ( f > t ) arc->_sweep = !_sweep;
-    if ( _large_arc && fabs(sweepAngle() * (t-f)) < M_PI) {
+    arc->_angles.setAngles(angleAt(f), angleAt(t));
+    if (f > t) arc->_angles.setSweep(!sweep());
+    if ( _large_arc && fabs(angularExtent() * (t-f)) <= M_PI) {
         arc->_large_arc = false;
     }
-    //arc->_start_angle = angleAt(f);
-    //arc->_end_angle = angleAt(t);
-    arc->_updateCenterAndAngles(); //TODO: be more clever
     return arc;
 }
 
 // the arc is the same but traversed in the opposite direction
 Curve *EllipticalArc::reverse() const
 {
+    using std::swap;
     EllipticalArc *rarc = static_cast<EllipticalArc*>(duplicate());
-    rarc->_sweep = !_sweep;
-    rarc->_initial_point = _final_point;
-    rarc->_final_point = _initial_point;
-    rarc->_start_angle = _end_angle;
-    rarc->_end_angle = _start_angle;
+    rarc->_angles.reverse();
+    swap(rarc->_initial_point, rarc->_final_point);
     return rarc;
 }
 
@@ -625,14 +588,14 @@ std::vector<double> EllipticalArc::allNearestTimes( Point const& p, double from,
         }
     }
 
-    double t = map_to_01( real_sol[mi1] );
+    double t = timeAtAngle(real_sol[mi1]);
     if ( !(t < from || t > to) )
     {
         result.push_back(t);
     }
 
     bool second_sol = false;
-    t = map_to_01( real_sol[mi2] );
+    t = timeAtAngle(real_sol[mi2]);
     if ( real_sol.size() == 4 && !(t < from || t > to) )
     {
         if ( result.empty() || are_near(mindistsq1, mindistsq2) )
@@ -757,25 +720,23 @@ void EllipticalArc::_updateCenterAndAngles()
 
     // if ip = sp, the arc contains no other points
     if (initialPoint() == finalPoint()) {
-        _start_angle = _end_angle = 0;
-        _ellipse.setRotationAngle(0);
+        _ellipse = Ellipse();
         _ellipse.setCenter(initialPoint());
-        _ellipse.setRays(0, 0);
-        _large_arc = _sweep = false;
+        _angles = AngleInterval();
+        _large_arc = false;
         return;
     }
 
     // rays should be positive
     _ellipse.setRays(std::fabs(ray(X)), std::fabs(ray(Y)));
 
-    if (ray(X) == 0 || ray(Y) == 0) {
-        _start_angle = 0;
-        _end_angle = M_PI;
+    if (isChord()) {
         _ellipse.setRays(L2(d) / 2, 0);
         _ellipse.setRotationAngle(atan2(d));
         _ellipse.setCenter(mid);
+        _angles.setAngles(0, M_PI);
+        _angles.setSweep(false);
         _large_arc = false;
-        _sweep = false;
         return;
     }
 
@@ -801,7 +762,7 @@ void EllipticalArc::_updateCenterAndAngles()
         // normally rad should never be negative, but numerical inaccuracy may cause this
         if (rad > 0) {
             rad = std::sqrt(rad);
-            if (_sweep == _large_arc) {
+            if (sweep() == _large_arc) {
                 rad = -rad;
             }
             c = rad * Point(r[X]*p[Y]/r[Y], -r[Y]*p[X]/r[X]);
@@ -816,12 +777,13 @@ void EllipticalArc::_updateCenterAndAngles()
     Point sp((p[X] - c[X]) / r[X], (p[Y] - c[Y]) / r[Y]);
     Point ep((-p[X] - c[X]) / r[X], (-p[Y] - c[Y]) / r[Y]);
     Point v(1, 0);
-    _start_angle = angle_between(v, sp);
-    double sweep_angle = angle_between(sp, ep);
-    if (!_sweep && sweep_angle > 0) sweep_angle -= 2*M_PI;
-    if (_sweep && sweep_angle < 0) sweep_angle += 2*M_PI;
 
-    _end_angle = _start_angle + sweep_angle;
+    _angles.setInitial(angle_between(v, sp));
+    _angles.setFinal(angle_between(v, ep));
+
+    /*double sweep_angle = angle_between(sp, ep);
+    if (!sweep() && sweep_angle > 0) sweep_angle -= 2*M_PI;
+    if (sweep() && sweep_angle < 0) sweep_angle += 2*M_PI;*/
 }
 
 D2<SBasis> EllipticalArc::toSBasis() const
@@ -852,7 +814,39 @@ D2<SBasis> EllipticalArc::toSBasis() const
     return arc;
 }
 
-void EllipticalArc::transform(Affine const& m)
+// All operations that do not contain skew can be evaulated
+// without passing through the implicit form of the ellipse,
+// which preserves precision.
+
+void EllipticalArc::operator*=(Translate const &tr)
+{
+    _initial_point *= tr;
+    _final_point *= tr;
+    _ellipse *= tr;
+}
+
+void EllipticalArc::operator*=(Scale const &s)
+{
+    _initial_point *= s;
+    _final_point *= s;
+    _ellipse *= s;
+}
+
+void EllipticalArc::operator*=(Rotate const &r)
+{
+    _initial_point *= r;
+    _final_point *= r;
+    _ellipse *= r;
+}
+
+void EllipticalArc::operator*=(Zoom const &z)
+{
+    _initial_point *= z;
+    _final_point *= z;
+    _ellipse *= z;
+}
+
+void EllipticalArc::operator*=(Affine const& m)
 {
     if (isChord()) {
         _initial_point *= m;
@@ -867,14 +861,14 @@ void EllipticalArc::transform(Affine const& m)
     _final_point *= m;
     _ellipse *= m;
     if (m.det() < 0) {
-        _sweep = !_sweep;
+        _angles.setSweep(!sweep());
     }
 
     // ellipse transformation does not preserve its functional form,
     // i.e. e.pointAt(0.5)*m and (e*m).pointAt(0.5) can be different.
     // We need to recompute start / end angles.
-    _start_angle = _ellipse.timeAt(_initial_point);
-    _end_angle = _ellipse.timeAt(_final_point);
+    _angles.setInitial(_ellipse.timeAt(_initial_point));
+    _angles.setFinal(_ellipse.timeAt(_final_point));
 }
 
 bool EllipticalArc::operator==(Curve const &c) const
@@ -888,7 +882,7 @@ bool EllipticalArc::operator==(Curve const &c) const
     if (rays() != other->rays()) return false;
     if (rotationAngle() != other->rotationAngle()) return false;
     if (_large_arc != other->_large_arc) return false;
-    if (_sweep != other->_sweep) return false;
+    if (sweep() != other->sweep()) return false;
     return true;
 }
 
@@ -897,16 +891,9 @@ void EllipticalArc::feed(PathSink &sink, bool moveto_initial) const
     if (moveto_initial) {
         sink.moveTo(_initial_point);
     }
-    sink.arcTo(ray(X), ray(Y), rotationAngle(), _large_arc, _sweep, _final_point);
-}
-
-Coord EllipticalArc::map_to_01(Coord angle) const
-{
-    return map_circular_arc_on_unit_interval(angle, initialAngle(),
-                                             finalAngle(), _sweep);
+    sink.arcTo(ray(X), ray(Y), rotationAngle(), _large_arc, sweep(), _final_point);
 }
 
-
 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 83909370e..fbb290dca 100644
--- a/src/2geom/elliptical-arc.h
+++ b/src/2geom/elliptical-arc.h
@@ -49,15 +49,14 @@
 namespace Geom 
 {
 
-class EllipticalArc : public Curve, public AngleInterval
+class EllipticalArc : public Curve
 {
 public:
-    /** @brief Creates an arc with all variables set to zero, and both flags to true. */
+    /** @brief Creates an arc with all variables set to zero. */
     EllipticalArc()
-        : AngleInterval(0, 0, true)
-        , _initial_point(0,0)
+        : _initial_point(0,0)
         , _final_point(0,0)
-        , _large_arc(true)
+        , _large_arc(false)
     {}
     /** @brief Create a new elliptical arc.
      * @param ip Initial point of the arc
@@ -72,61 +71,73 @@ public:
                    Coord rot_angle, bool large_arc, bool sweep,
                    Point const &fp
                  )
-        : AngleInterval(0,0,sweep)
-        , _initial_point(ip)
+        : _initial_point(ip)
         , _final_point(fp)
         , _ellipse(0, 0, r[X], r[Y], rot_angle)
+        , _angles(0, 0, sweep)
         , _large_arc(large_arc)
     {
         _updateCenterAndAngles();
     }
 
+    /// Create a new elliptical arc, giving the ellipse's rays as separate coordinates.
     EllipticalArc( Point const &ip, Coord rx, Coord ry,
                    Coord rot_angle, bool large_arc, bool sweep,
                    Point const &fp
                  )
-        : AngleInterval(0,0,sweep)
-        , _initial_point(ip)
+        : _initial_point(ip)
         , _final_point(fp)
         , _ellipse(0, 0, rx, ry, rot_angle)
+        , _angles(0, 0, sweep)
         , _large_arc(large_arc)
     {
         _updateCenterAndAngles();
     }
 
-    // methods new to EllipticalArc go here
-
-    /// @name Retrieve and modify parameters
+    /// @name Retrieve basic information
     /// @{
-    /** @brief Get the interval of angles the arc contains
-     * @return The interval between the final and initial angles of the arc */
-    Interval angleInterval() const { return Interval(initialAngle(), finalAngle()); }
+
     /** @brief Get a coordinate of the elliptical arc's center.
      * @param d The dimension to retrieve
      * @return The selected coordinate of the center */
-    /** @brief Get the defining ellipse's rotation
-     * @return Angle between the +X ray of the ellipse and the +X axis */
-    Angle rotationAngle() const {
-        return _ellipse.rotationAngle();
-    }
+    Coord center(Dim2 d) const { return _ellipse.center(d); }
+
+    /** @brief Get the arc's center
+     * @return The arc's center, situated on the intersection of the ellipse's rays */
+    Point center() const { return _ellipse.center(); }
+
     /** @brief Get one of the ellipse's rays
      * @param d Dimension to retrieve
      * @return The selected ray of the ellipse */
     Coord ray(Dim2 d) const { return _ellipse.ray(d); }
+
     /** @brief Get both rays as a point
      * @return Point with X equal to the X ray and Y to Y ray */
     Point rays() const { return _ellipse.rays(); }
+
+    /** @brief Get the defining ellipse's rotation
+     * @return Angle between the +X ray of the ellipse and the +X axis */
+    Angle rotationAngle() const {
+        return _ellipse.rotationAngle();
+    }
+
     /** @brief Whether the arc is larger than half an ellipse.
      * @return True if the arc is larger than \f$\pi\f$, false otherwise */
     bool largeArc() const { return _large_arc; }
+
     /** @brief Whether the arc turns clockwise
      * @return True if the arc makes a clockwise turn when going from initial to final
      *         point, false otherwise */
-    bool sweep() const { return _sweep; }
-    /** @brief Get the line segment connecting the arc's endpoints.
-     * @return A linear segment with initial and final point correspoding to those of the arc. */
-    LineSegment chord() const { return LineSegment(_initial_point, _final_point); }
-    /** @brief Change the arc's parameters. */ 
+    bool sweep() const { return _angles.sweep(); }
+
+    Angle initialAngle() const { return _angles.initialAngle(); }
+    Angle finalAngle() const { return _angles.finalAngle(); }
+    /// @}
+
+    /// @name Modify parameters
+    /// @{
+
+    /// Change all of the arc's parameters.
     void set( Point const &ip, double rx, double ry,
               double rot_angle, bool large_arc, bool sweep,
               Point const &fp
@@ -136,10 +147,26 @@ public:
         _final_point = fp;
         _ellipse.setRays(rx, ry);
         _ellipse.setRotationAngle(rot_angle);
+        _angles.setSweep(sweep);
+        _large_arc = large_arc;
+        _updateCenterAndAngles();
+    }
+
+    /// Change all of the arc's parameters.
+    void set( Point const &ip, Point const &r,
+              Angle rot_angle, bool large_arc, bool sweep,
+              Point const &fp
+            )
+    {
+        _initial_point = ip;
+        _final_point = fp;
+        _ellipse.setRays(r);
+        _ellipse.setRotationAngle(rot_angle);
+        _angles.setSweep(sweep);
         _large_arc = large_arc;
-        _sweep = sweep;
         _updateCenterAndAngles();
     }
+
     /** @brief Change the initial and final point in one operation.
      * This method exists because modifying any of the endpoints causes rather costly
      * recalculations of the center and extreme angles.
@@ -152,52 +179,76 @@ public:
     }
     /// @}
 
-    /// @name Access computed parameters of the arc
-    /// @{
-    Coord center(Dim2 d) const { return _ellipse.center(d); }
-    /** @brief Get the arc's center
-     * @return The arc's center, situated on the intersection of the ellipse's rays */
-    Point center() const { return _ellipse.center(); }
-    /// @}
-    
-    /// @name Angular evaluation
+    /// @name Evaluate the arc as a function
     /// @{
     /** Check whether the arc contains the given angle
      * @param t The angle to check
      * @return True if the arc contains the angle, false otherwise */
-    bool containsAngle(Coord angle) const;
+    bool containsAngle(Angle angle) const { return _angles.contains(angle); }
+
     /** @brief Evaluate the arc at the specified angular coordinate
      * @param t Angle
      * @return Point corresponding to the given angle */
     Point pointAtAngle(Coord t) const;
+
     /** @brief Evaluate one of the arc's coordinates at the specified angle
      * @param t Angle
      * @param d The dimension to retrieve
      * @return Selected coordinate of the arc at the specified angle */
     Coord valueAtAngle(Coord t, Dim2 d) const;
-    /** @brief Retrieve the unit circle transform.
+
+    /// Compute the curve time value corresponding to the given angular value.
+    Coord timeAtAngle(Angle a) const { return _angles.timeAtAngle(a); }
+
+    /// Compute the angular domain value corresponding to the given time value.
+    Angle angleAt(Coord t) const { return _angles.angleAt(t); }
+
+    /** @brief Compute the amount by which the angle parameter changes going from start to end.
+     * This has range \f$(-2\pi, 2\pi)\f$ and thus cannot be represented as instance
+     * of the class Angle. Add this to the initial angle to obtain the final angle. */
+    Coord sweepAngle() const { return _angles.sweepAngle(); }
+
+    /** @brief Get the elliptical angle spanned by the arc.
+     * This is basically the absolute value of sweepAngle(). */
+    Coord angularExtent() const { return _angles.extent(); }
+
+    /// Get the angular interval of the arc.
+    AngleInterval angularInterval() const { return _angles; }
+
+    /// Evaluate the arc in the curve domain, i.e. \f$[0, 1]\$.
+    virtual Point pointAt(Coord t) const;
+
+    /// Evaluate a single coordinate on the arc in the curve domain.
+    virtual Coord valueAt(Coord t, Dim2 d) const;
+
+    /** @brief Compute a transform that maps the unit circle to the arc's ellipse.
      * Each ellipse can be interpreted as a translated, scaled and rotate unit circle.
      * This function returns the transform that maps the unit circle to the arc's ellipse.
      * @return Transform from unit circle to the arc's ellipse */
-    Affine unitCircleTransform() const;
-    Affine inverseUnitCircleTransform() const;
+    Affine unitCircleTransform() const {
+        Affine result = _ellipse.unitCircleTransform();
+        return result;
+    }
+
+    /** @brief Compute a transform that maps the arc's ellipse to the unit circle. */
+    Affine inverseUnitCircleTransform() const {
+        Affine result = _ellipse.inverseUnitCircleTransform();
+        return result;
+    }
     /// @}
 
+    /// @name Deal with degenerate ellipses.
+    /// @{
     /** @brief Check whether both rays are nonzero.
      * If they are not, the arc is represented as a line segment instead. */
     bool isChord() const {
         return ray(X) == 0 || ray(Y) == 0;
     }
 
-    std::pair<EllipticalArc, EllipticalArc> subdivide(Coord t) const {
-        EllipticalArc* arc1 = static_cast<EllipticalArc*>(portion(0, t));
-        EllipticalArc* arc2 = static_cast<EllipticalArc*>(portion(t, 1));
-        assert( arc1 != NULL && arc2 != NULL);
-        std::pair<EllipticalArc, EllipticalArc> arc_pair(*arc1, *arc2);        
-        delete arc1;
-        delete arc2;
-        return arc_pair;
-    }
+    /** @brief Get the line segment connecting the arc's endpoints.
+     * @return A linear segment with initial and final point correspoding to those of the arc. */
+    LineSegment chord() const { return LineSegment(_initial_point, _final_point); }
+    /// @}
 
     // implementation of overloads goes here
     virtual Point initialPoint() const { return _initial_point; }
@@ -226,58 +277,50 @@ public:
     virtual std::vector<double> roots(double v, Dim2 d) const;
 #ifdef HAVE_GSL
     virtual std::vector<double> allNearestTimes( Point const& p, double from = 0, double to = 1 ) const;
-#endif
     virtual double nearestTime( Point const& p, double from = 0, double to = 1 ) const {
         if ( are_near(ray(X), ray(Y)) && are_near(center(), p) ) {
             return from;
         }
         return allNearestTimes(p, from, to).front();
     }
+#endif
     virtual std::vector<CurveIntersection> intersect(Curve const &other, Coord eps=EPSILON) const;
     virtual int degreesOfFreedom() const { return 7; }
     virtual Curve *derivative() const;
-    virtual void transform(Affine const &m);
-    virtual Curve &operator*=(Translate const &m) {
-        _initial_point *= m;
-        _final_point *= m;
-        _ellipse *= m;
-        return *this;
-    }
 
+    using Curve::operator*=;
+    virtual void operator*=(Translate const &tr);
+    virtual void operator*=(Scale const &s);
+    virtual void operator*=(Rotate const &r);
+    virtual void operator*=(Zoom const &z);
+    virtual void operator*=(Affine const &m);
 
-    /**
-    *  The size of the returned vector equals n+1.
-    */
     virtual std::vector<Point> pointAndDerivatives(Coord t, unsigned int n) const;
-
     virtual D2<SBasis> toSBasis() const;
-    virtual double valueAt(Coord t, Dim2 d) const {
-        if (isChord()) return chord().valueAt(t, d);
-        return valueAtAngle(angleAt(t), d);
-    }
-    virtual Point pointAt(Coord t) const;
-    virtual Curve* portion(double f, double t) const;
-    virtual Curve* reverse() const;
+    virtual Curve *portion(double f, double t) const;
+    virtual Curve *reverse() const;
     virtual bool operator==(Curve const &c) const;
     virtual void feed(PathSink &sink, bool moveto_initial) const;
 
-protected:
+private:
     void _updateCenterAndAngles();
     void _filterIntersections(std::vector<ShapeIntersection> &xs, bool is_first) const;
 
     Point _initial_point, _final_point;
     Ellipse _ellipse;
+    AngleInterval _angles;
     bool _large_arc;
-
-private:
-    Coord map_to_01(Coord angle) const; 
 }; // end class EllipticalArc
 
 
 // implemented in elliptical-arc-from-sbasis.cpp
+/** @brief Fit an elliptical arc to an SBasis fragment.
+ * @relates EllipticalArc */
 bool arc_from_sbasis(EllipticalArc &ea, D2<SBasis> const &in,
                      double tolerance = EPSILON, unsigned num_samples = 20);
 
+/** @brief Debug output for elliptical arcs.
+ * @relates EllipticalArc */
 std::ostream &operator<<(std::ostream &out, EllipticalArc const &ea);
 
 } // end namespace Geom
diff --git a/src/2geom/intersection-graph.cpp b/src/2geom/intersection-graph.cpp
index ff96c28a8..e18561f67 100644
--- a/src/2geom/intersection-graph.cpp
+++ b/src/2geom/intersection-graph.cpp
@@ -39,7 +39,7 @@
 
 namespace Geom {
 
-struct IntersectionVertexLess {
+struct PathIntersectionGraph::IntersectionVertexLess {
     bool operator()(IntersectionVertex const &a, IntersectionVertex const &b) const {
         return a.pos < b.pos;
     }
@@ -49,7 +49,8 @@ struct IntersectionVertexLess {
  * @brief Intermediate data for computing Boolean operations on paths.
  *
  * This class implements the Greiner-Hormann clipping algorithm,
- * with improvements by Foster and Overfelt.
+ * with improvements inspired by Foster and Overfelt as well as some
+ * original contributions.
  *
  * @ingroup Paths
  */
@@ -69,15 +70,16 @@ PathIntersectionGraph::PathIntersectionGraph(PathVector const &a, PathVector con
     }
 
     std::vector<PVIntersection> pxs = _a.intersect(_b, 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.size() % 2) return;
 
     // prepare intersection lists for each path component
     for (std::size_t i = 0; i < _a.size(); ++i) {
-        _xalists.push_back(new IntersectionList());
+        _apaths.push_back(new PathData());
     }
     for (std::size_t i = 0; i < _b.size(); ++i) {
-        _xblists.push_back(new IntersectionList());
+        _bpaths.push_back(new PathData());
     }
 
     for (std::size_t i = 0; i < pxs.size(); ++i) {
@@ -92,30 +94,32 @@ PathIntersectionGraph::PathIntersectionGraph(PathVector const &a, PathVector con
         xb->neighbor = xa;
         _xs.push_back(xa);
         _xs.push_back(xb);
-        _xalists[xa->pos.path_index].push_back(*xa);
-        _xblists[xb->pos.path_index].push_back(*xb);
+        _apaths[xa->pos.path_index].xlist.push_back(*xa);
+        _bpaths[xb->pos.path_index].xlist.push_back(*xb);
     }
 
-    for (std::size_t i = 0; i < _xalists.size(); ++i) {
-        _xalists[i].sort(IntersectionVertexLess());
+    for (std::size_t i = 0; i < _apaths.size(); ++i) {
+        _apaths[i].xlist.sort(IntersectionVertexLess());
     }
-    for (std::size_t i = 0; i < _xblists.size(); ++i) {
-        _xblists[i].sort(IntersectionVertexLess());
+    for (std::size_t i = 0; i < _bpaths.size(); ++i) {
+        _bpaths[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<IntersectionList> &ls = npv ? _xblists : _xalists;
+        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) {
-            for (Iter i = ls[li].begin(); i != ls[li].end(); ++i) {
+            IntersectionList &xl = ls[li].xlist;
+            for (Iter i = xl.begin(); i != xl.end(); ++i) {
                 Iter n = boost::next(i);
-                if (n == ls[li].end()) {
-                    n = ls[li].begin();
+                if (n == xl.end()) {
+                    n = xl.begin();
                 }
                 std::size_t pi = i->pos.path_index;
 
@@ -123,7 +127,6 @@ PathIntersectionGraph::PathIntersectionGraph(PathVector const &a, PathVector con
                 PathTime mid = ival.inside(precision);
 
                 // TODO check for degenerate cases
-                // requires changes in the winding routine
                 int w = other.winding(pv[pi].pointAt(mid));
                 if (w % 2) {
                     i->next = POINT_INSIDE;
@@ -134,9 +137,22 @@ PathIntersectionGraph::PathIntersectionGraph(PathVector const &a, PathVector con
                 }
             }
 
-            // assign exit / entry flags
-            for (Iter i = ls[li].begin(); i != ls[li].end(); ++i) {
-                i->entry = ((i->next == POINT_INSIDE) && (i->previous == POINT_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;
+                    }
+                } else {
+                    i->entry = ((i->next == POINT_INSIDE) && (i->previous == POINT_OUTSIDE));
+                    ++i;
+                }
             }
         }
     }
@@ -162,6 +178,7 @@ PathVector PathIntersectionGraph::getAminusB()
 {
     PathVector result = _getResult(false, true);
     _handleNonintersectingPaths(result, 0, false);
+    _handleNonintersectingPaths(result, 1, true);
     return result;
 }
 
@@ -169,6 +186,7 @@ PathVector PathIntersectionGraph::getBminusA()
 {
     PathVector result = _getResult(true, false);
     _handleNonintersectingPaths(result, 1, false);
+    _handleNonintersectingPaths(result, 0, true);
     return result;
 }
 
@@ -180,13 +198,22 @@ PathVector PathIntersectionGraph::getXOR()
     return r1;
 }
 
+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();
+    }
+    return result;
+}
+
 std::vector<Point> PathIntersectionGraph::intersectionPoints() const
 {
     std::vector<Point> result;
 
     typedef IntersectionList::const_iterator Iter;
-    for (std::size_t i = 0; i < _xalists.size(); ++i) {
-        for (Iter j = _xalists[i].begin(); j != _xalists[i].end(); ++j) {
+    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);
         }
     }
@@ -201,9 +228,9 @@ PathVector PathIntersectionGraph::_getResult(bool enter_a, bool enter_b)
 
     // reset processed status
     for (unsigned npv = 0; npv < 2; ++npv) {
-        boost::ptr_vector<IntersectionList> &ls = npv ? _xblists : _xalists;
+        boost::ptr_vector<PathData> &ls = npv ? _bpaths : _apaths;
         for (std::size_t li = 0; li < ls.size(); ++li) {
-            for (Iter k = ls[li].begin(); k != ls[li].end(); ++k) {
+            for (Iter k = ls[li].xlist.begin(); k != ls[li].xlist.end(); ++k) {
                 k->processed = false;
             }
         }
@@ -213,7 +240,7 @@ PathVector PathIntersectionGraph::_getResult(bool enter_a, bool enter_b)
 
     while (true) {
         PathVector const *cur = &_a, *other = &_b;
-        boost::ptr_vector<IntersectionList> *lscur = &_xalists, *lsother = &_xblists;
+        boost::ptr_vector<PathData> *lscur = &_apaths, *lsother = &_bpaths;
 
         // find unprocessed intersection
         Iter i;
@@ -239,14 +266,14 @@ PathVector PathIntersectionGraph::_getResult(bool enter_a, bool enter_b)
 
             // get next intersection
             if (reverse) {
-                if (i == (*lscur)[pi].begin()) {
-                    i = (*lscur)[pi].end();
+                if (i == (*lscur)[pi].xlist.begin()) {
+                    i = (*lscur)[pi].xlist.end();
                 }
                 --i;
             } else {
                 ++i;
-                if (i == (*lscur)[pi].end()) {
-                    i = (*lscur)[pi].begin();
+                if (i == (*lscur)[pi].xlist.end()) {
+                    i = (*lscur)[pi].xlist.begin();
                 }
             }
 
@@ -263,7 +290,7 @@ PathVector PathIntersectionGraph::_getResult(bool enter_a, bool enter_b)
             n_processed += 2;
 
             // switch to the other path
-            i = (*lsother)[i->neighbor->pos.path_index].iterator_to(*i->neighbor);
+            i = (*lsother)[i->neighbor->pos.path_index].xlist.iterator_to(*i->neighbor);
             std::swap(lscur, lsother);
             std::swap(cur, other);
         }
@@ -272,7 +299,7 @@ PathVector PathIntersectionGraph::_getResult(bool enter_a, bool enter_b)
         assert(!result.back().empty());
     }
 
-    assert(n_processed == _xs.size());
+    assert(n_processed == size() * 2);
 
     return result;
 }
@@ -285,15 +312,29 @@ 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<IntersectionList> const &ls = which ? _xblists : _xalists;
+    boost::ptr_vector<PathData> const &ls = which ? _bpaths : _apaths;
     PathVector const &cur = which ? _b : _a;
     PathVector const &other = which ? _a : _b;
 
     for (std::size_t i = 0; i < cur.size(); ++i) {
-        if (!ls.empty() && !ls[i].empty()) continue;
+        // the path data vector might have been left empty if there were no intersections at all
+        bool has_path_data = !ls.empty();
+        // Skip if the path has intersections
+        if (has_path_data && !ls[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) {
+            path_inside = true;
+        } else if (has_path_data && ls[i].removed_out != 0) {
+            path_inside = false;
+        } else {
+            int w = other.winding(cur[i].initialPoint());
+            path_inside = w % 2 != 0;
+        }
 
-        int w = other.winding(cur[i].initialPoint());
-        bool path_inside = w % 2 != 0;
         if (path_inside == inside) {
             result.push_back(cur[i]);
         }
@@ -304,11 +345,16 @@ bool PathIntersectionGraph::_findUnprocessed(IntersectionList::iterator &result)
 {
     typedef IntersectionList::iterator Iter;
 
-    Iter it;
+    Iter it, last;
 
-    for (std::size_t k = 0; k < _xalists.size(); ++k) {
-        it = _xalists[k].begin();
-        for (; it != _xalists[k].end(); ++it) {
+    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;
@@ -319,6 +365,21 @@ bool PathIntersectionGraph::_findUnprocessed(IntersectionList::iterator &result)
     return false;
 }
 
+
+std::ostream &operator<<(std::ostream &os, PathIntersectionGraph const &pig)
+{
+    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) {
+            os << j->pos << " - " << j->neighbor->pos << " @ " << j->p << "\n";
+        }
+    }
+    return os;
+}
+
 } // namespace Geom
 
 /*
diff --git a/src/2geom/intersection-graph.h b/src/2geom/intersection-graph.h
index 44beaf46b..bd9aaee81 100644
--- a/src/2geom/intersection-graph.h
+++ b/src/2geom/intersection-graph.h
@@ -43,33 +43,6 @@
 
 namespace Geom {
 
-enum InOutFlag {
-    POINT_INSIDE,
-    POINT_OUTSIDE,
-    POINT_ON_EDGE
-};
-
-struct IntersectionVertex {
-    boost::intrusive::list_member_hook<> _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
-};
-
-typedef boost::intrusive::list
-    < IntersectionVertex
-    , boost::intrusive::member_hook
-        < IntersectionVertex
-        , boost::intrusive::list_member_hook<>
-        , &IntersectionVertex::_hook
-        >
-    > IntersectionList;
-
-
 class PathIntersectionGraph
 {
     // this is called PathIntersectionGraph so that we can also have a class for polygons,
@@ -83,19 +56,63 @@ public:
     PathVector getBminusA();
     PathVector getXOR();
 
+    /// Returns the number of intersections used when computing Boolean operations.
+    std::size_t size() const;
     std::vector<Point> intersectionPoints() const;
 
 private:
+    enum InOutFlag {
+        POINT_INSIDE,
+        POINT_OUTSIDE
+    };
+
+    struct IntersectionVertex {
+        boost::intrusive::list_member_hook<> _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
+    };
+
+    typedef boost::intrusive::list
+        < IntersectionVertex
+        , boost::intrusive::member_hook
+            < IntersectionVertex
+            , boost::intrusive::list_member_hook<>
+            , &IntersectionVertex::_hook
+            >
+        > IntersectionList;
+
+    struct PathData {
+        IntersectionList xlist;
+        unsigned removed_in;
+        unsigned removed_out;
+
+        PathData()
+            : removed_in(0)
+            , removed_out(0)
+        {}
+    };
+
+    struct IntersectionVertexLess;
+
     PathVector _getResult(bool enter_a, bool enter_b);
     void _handleNonintersectingPaths(PathVector &result, int which, bool inside);
     bool _findUnprocessed(IntersectionList::iterator &result);
 
     PathVector _a, _b;
     boost::ptr_vector<IntersectionVertex> _xs;
-    boost::ptr_vector<IntersectionList> _xalists;
-    boost::ptr_vector<IntersectionList> _xblists;
+    boost::ptr_vector<PathData> _apaths;
+    boost::ptr_vector<PathData> _bpaths;
+
+    friend std::ostream &operator<<(std::ostream &, PathIntersectionGraph const &);
 };
 
+std::ostream &operator<<(std::ostream &os, PathIntersectionGraph const &pig);
+
 } // namespace Geom
 
 #endif // SEEN_LIB2GEOM_PATH_GRAPH_H
diff --git a/src/2geom/intersection.h b/src/2geom/intersection.h
index ae4b50dd1..bbce19947 100644
--- a/src/2geom/intersection.h
+++ b/src/2geom/intersection.h
@@ -44,6 +44,7 @@ namespace Geom {
  */
 template <typename TimeA = Coord, typename TimeB = TimeA>
 class Intersection
+    : boost::totally_ordered< Intersection<TimeA, TimeB> >
 {
 public:
     /** @brief Construct from shape references and time values.
@@ -79,6 +80,17 @@ public:
         swap(a._point, b._point);
     }
 
+    bool operator==(Intersection const &other) const {
+        if (first != other.first) return false;
+        if (second != other.second) return false;
+        return true;
+    }
+    bool operator<(Intersection const &other) const {
+        if (first < other.first) return true;
+        if (first == other.first && second < other.second) return true;
+        return false;
+    }
+
 public:
     /// First shape and time value.
     TimeA first;
diff --git a/src/2geom/line.cpp b/src/2geom/line.cpp
index 8bea33638..bada8ef38 100644
--- a/src/2geom/line.cpp
+++ b/src/2geom/line.cpp
@@ -60,6 +60,7 @@ namespace Geom
  * satisfy the line equation with the given coefficients. */
 void Line::setCoefficients (Coord a, Coord b, Coord c)
 {
+    // degenerate case
     if (a == 0 && b == 0) {
         if (c != 0) {
             THROW_LOGICALERROR("the passed coefficients give the empty set");
@@ -68,22 +69,36 @@ void Line::setCoefficients (Coord a, Coord b, Coord c)
         _final = Point(0,0);
         return;
     }
+
+    // The way final / initial points are set based on coefficients is somewhat unusual.
+    // This is done to make sure that calling coefficients() will give back
+    // (almost) the same values.
+
+    // vertical case
     if (a == 0) {
         // b must be nonzero
-        _initial = Point(0, -c / b);
+        _initial = Point(-b/2, -c / b);
         _final = _initial;
-        _final[X] = 1;
+        _final[X] = b/2;
         return;
     }
+
+    // horizontal case
     if (b == 0) {
-        _initial = Point(-c / a, 0);
+        _initial = Point(-c / a, a/2);
         _final = _initial;
-        _final[Y] = 1;
+        _final[Y] = -a/2;
         return;
     }
 
-    _initial = Point(-c / a, 0);
-    _final = Point(0, -c / b);
+    // This gives reasonable results regardless of the magnitudes of a, b and c.
+    _initial = Point(-b/2,a/2);
+    _final = Point(b/2,-a/2);
+
+    Point offset(-c/(2*a), -c/(2*b));
+
+    _initial += offset;
+    _final += offset;
 }
 
 void Line::coefficients(Coord &a, Coord &b, Coord &c) const
@@ -94,15 +109,18 @@ void Line::coefficients(Coord &a, Coord &b, Coord &c) const
     c = cross(_initial, _final);
 }
 
-/** @brief Get the line equation coefficients of this line.
+/** @brief Get the implicit line equation coefficients.
+ * Note that conversion to implicit form always causes loss of
+ * precision when dealing with lines that start far from the origin
+ * and end very close to it. It is recommended to normalize the line
+ * before converting it to implicit form.
  * @return Vector with three values corresponding to the A, B and C
  *         coefficients of the line equation for this line. */
 std::vector<Coord> Line::coefficients() const
 {
-    Coord c[3];
+    std::vector<Coord> c(3);
     coefficients(c[0], c[1], c[2]);
-    std::vector<Coord> coeff(c, c+3);
-    return coeff;
+    return c;
 }
 
 /** @brief Find intersection with an axis-aligned line.
diff --git a/src/2geom/line.h b/src/2geom/line.h
index 5d92c436d..7d4766e12 100644
--- a/src/2geom/line.h
+++ b/src/2geom/line.h
@@ -49,15 +49,7 @@ namespace Geom
 
 // class docs in cpp file
 class Line
-    : boost::equality_comparable1< Line
-    , MultipliableNoncommutative< Line, Translate
-    , MultipliableNoncommutative< Line, Scale
-    , MultipliableNoncommutative< Line, Rotate
-    , MultipliableNoncommutative< Line, HShear
-    , MultipliableNoncommutative< Line, VShear
-    , MultipliableNoncommutative< Line, Zoom
-    , MultipliableNoncommutative< Line, Affine
-      > > > > > > > >
+    : boost::equality_comparable< Line >
 {
 private:
     Point _initial;
@@ -203,8 +195,14 @@ public:
         return _initial[X] == _final[X];
     }
 
-    /** @brief Reparametrize the line so that it has unit speed. */
+    /** @brief Reparametrize the line so that it has unit speed.
+     * Note that the direction of the line may also change. */
     void normalize() {
+        // this helps with the nasty case of a line that starts somewhere far
+        // and ends very close to the origin
+        if (L2sq(_final) < L2sq(_initial)) {
+            std::swap(_initial, _final);
+        }
         Point v = _final - _initial;
         v.normalize();
         _final = _initial + v;
@@ -380,6 +378,14 @@ public:
         if (distance(pointAt(nearestTime(other._final)), other._final) != 0) return false;
         return true;
     }
+
+    template <typename T>
+    friend Line operator*(Line const &l, T const &tr) {
+        BOOST_CONCEPT_ASSERT((TransformConcept<T>));
+        Line result(l);
+        result *= tr;
+        return result;
+    }
 }; // end class Line
 
 /** @brief Removes intersections outside of the unit interval.
diff --git a/src/2geom/linear.h b/src/2geom/linear.h
index cd7108835..b0306fb9f 100644
--- a/src/2geom/linear.h
+++ b/src/2geom/linear.h
@@ -51,7 +51,7 @@ namespace Geom {
 class SBasis;
 
 /**
- * @brief Linear function fragment
+ * @brief Function that interpolates linearly between two values.
  * @ingroup Fragments
  */
 class Linear {
diff --git a/src/2geom/path.cpp b/src/2geom/path.cpp
index 5fbc65b3e..836e65013 100644
--- a/src/2geom/path.cpp
+++ b/src/2geom/path.cpp
@@ -380,20 +380,6 @@ bool Path::operator==(Path const &other) const
     return *_curves == *other._curves;
 }
 
-Path &Path::operator*=(Affine const &m)
-{
-    _unshare();
-    Sequence::iterator last = _curves->end() - 1;
-    Sequence::iterator it;
-
-    for (it = _curves->begin(); it != last; ++it) {
-        it->transform(m);
-    }
-    _closing_seg->transform(m);
-    checkContinuity();
-    return *this;
-}
-
 void Path::start(Point const &p) {
     if (_curves->size() > 1) {
         clear();
@@ -507,7 +493,7 @@ protected:
             int ib = i->bound.index;
 
             if (which == 0) {
-                cx = record.item->intersect(*i->item);
+                cx = record.item->intersect(*i->item, _precision);
             } else {
                 cx = i->item->intersect(*record.item, _precision);
                 std::swap(ia, ib);
@@ -535,7 +521,14 @@ std::vector<PathIntersection> Path::intersect(Path const &other, Coord precision
     CurveSweeper sweeper(*this, other, result, precision);
     sweeper.process();
 
-    // TODO: remove multiple intersections within precision of each other?
+    // preprocessing to remove duplicate intersections at endpoints
+    for (std::size_t i = 0; i < result.size(); ++i) {
+        result[i].first.normalizeForward(size());
+        result[i].second.normalizeForward(other.size());
+    }
+    std::sort(result.begin(), result.end());
+    result.erase(std::unique(result.begin(), result.end()), result.end());
+
     return result;
 }
 
diff --git a/src/2geom/path.h b/src/2geom/path.h
index a6473e0d2..a2d1e751e 100644
--- a/src/2geom/path.h
+++ b/src/2geom/path.h
@@ -159,7 +159,7 @@ struct PathTime
 };
 
 inline std::ostream &operator<<(std::ostream &os, PathTime const &pos) {
-    os << pos.curve_index << ": " << pos.t;
+    os << pos.curve_index << ": " << format_coord_nice(pos.t);
     return os;
 }
 
@@ -313,9 +313,7 @@ struct ShapeTraits<Path> {
  *
  * @ingroup Paths */
 class Path
-    : boost::equality_comparable1< Path
-    , MultipliableNoncommutative< Path, Affine
-      > >
+    : boost::equality_comparable< Path >
 {
 public:
     typedef PathInternal::Sequence Sequence;
@@ -493,8 +491,24 @@ public:
     /// Test paths for exact equality.
     bool operator==(Path const &other) const;
 
-    /// Apply an affine transform.
-    Path &operator*=(Affine const &m);
+    /// Apply a transform to each curve.
+    template <typename T>
+    Path &operator*=(T const &tr) {
+        BOOST_CONCEPT_ASSERT((TransformConcept<T>));
+        _unshare();
+        for (std::size_t i = 0; i < _curves->size(); ++i) {
+            (*_curves)[i] *= tr;
+        }
+        return *this;
+    }
+
+    template <typename T>
+    friend Path operator*(Path const &path, T const &tr) {
+        BOOST_CONCEPT_ASSERT((TransformConcept<T>));
+        Path result(path);
+        result *= tr;
+        return result;
+    }
 
     /** @brief Get the allowed range of time values.
      * @return Values for which pointAt() and valueAt() yield valid results. */
diff --git a/src/2geom/pathvector.h b/src/2geom/pathvector.h
index 6636cbf2e..108f2aa05 100644
--- a/src/2geom/pathvector.h
+++ b/src/2geom/pathvector.h
@@ -81,6 +81,11 @@ struct PathVectorTime
     }
 };
 
+inline std::ostream &operator<<(std::ostream &os, PathVectorTime const &pvt) {
+    os << pvt.path_index << ": " << pvt.asPathTime();
+    return os;
+}
+
 typedef Intersection<PathVectorTime> PathVectorIntersection;
 typedef PathVectorIntersection PVIntersection; ///< Alias to save typing
 
diff --git a/src/2geom/piecewise.h b/src/2geom/piecewise.h
index 981d67144..a5a65a9fe 100644
--- a/src/2geom/piecewise.h
+++ b/src/2geom/piecewise.h
@@ -42,12 +42,14 @@
 
 namespace Geom {
 /**
- * %Piecewise function class.
+ * @brief Function defined as discrete pieces.
+ *
  * The Piecewise class manages a sequence of elements of a type as segments and
  * the ’cuts’ between them. These cuts are time values which separate the pieces.
  * This function representation allows for more interesting functions, as it provides
  * a viable output for operations such as inversion, which may require multiple
  * SBasis to properly invert the original.
+ *
  * As for technical details, while the actual SBasis segments begin on the first
  * cut and end on the last, the function is defined throughout all inputs by ex-
  * tending the first and last segments. The exact switching between segments is
@@ -62,6 +64,8 @@ namespace Geom {
  *      s_n,& c_n <= t
  *      \end{array}\right.
  * \f]
+ *
+ * @ingroup Fragments
  */
 template <typename T>
 class Piecewise {
diff --git a/src/2geom/polynomial.cpp b/src/2geom/polynomial.cpp
index a0689f0c5..ca2389f80 100644
--- a/src/2geom/polynomial.cpp
+++ b/src/2geom/polynomial.cpp
@@ -34,6 +34,7 @@
 
 #include <algorithm>
 #include <2geom/polynomial.h>
+#include <2geom/math-utils.h>
 #include <math.h>
 
 #ifdef HAVE_GSL
@@ -235,12 +236,17 @@ std::vector<Coord> solve_quadratic(Coord a, Coord b, Coord c)
 
     if (delta == 0) {
         // one root
-        result.push_back(-0.5 * b / a);
+        result.push_back(-b / (2*a));
     } else if (delta > 0) {
         // two roots
         Coord delta_sqrt = sqrt(delta);
-        result.push_back((-b + delta_sqrt)/(2*a));
-        result.push_back((-b - delta_sqrt)/(2*a));
+
+        // Use different formulas depending on sign of b to preserve
+        // numerical stability. See e.g.:
+        // http://people.csail.mit.edu/bkph/articles/Quadratics.pdf
+        Coord t = -0.5 * (b + sgn(b) * delta_sqrt);
+        result.push_back(t / a);
+        result.push_back(c / t);
     }
     // no roots otherwise
 
diff --git a/src/2geom/polynomial.h b/src/2geom/polynomial.h
index aeac2f3fa..5ab2aa4c8 100644
--- a/src/2geom/polynomial.h
+++ b/src/2geom/polynomial.h
@@ -44,6 +44,8 @@
 
 namespace Geom {
 
+/** @brief Polynomial in canonical (monomial) basis.
+ * @ingroup Fragments */
 class Poly : public std::vector<double>{
 public:
     // coeff; // sum x^i*coeff[i]
diff --git a/src/2geom/sbasis-curve.h b/src/2geom/sbasis-curve.h
index 17ee8f4c9..affe7edc0 100644
--- a/src/2geom/sbasis-curve.h
+++ b/src/2geom/sbasis-curve.h
@@ -119,7 +119,8 @@ public:
         return new SBasisCurve(Geom::portion(inner, f, t));
     }
 
-    virtual void transform(Affine const &m) { inner = inner * m; }
+    using Curve::operator*=;
+    virtual void operator*=(Affine const &m) { inner = inner * m; }
 
     virtual Curve *derivative() const {
         return new SBasisCurve(Geom::derivative(inner));
diff --git a/src/2geom/sbasis-math.cpp b/src/2geom/sbasis-math.cpp
index 97cdf45ce..896eb18a7 100644
--- a/src/2geom/sbasis-math.cpp
+++ b/src/2geom/sbasis-math.cpp
@@ -297,11 +297,11 @@ Piecewise<SBasis> reciprocalOnDomain(Interval range, double tol){
     if (a<=tol){
         reciprocal_fn.push_cut(0);
         int i0=(int) floor(std::log(tol)/std::log(R));
-        a=pow(R,i0);
+        a = std::pow(R,i0);
         reciprocal_fn.push(Linear(1/a),a);
     }else{
         int i0=(int) floor(std::log(a)/std::log(R));
-        a=pow(R,i0);
+        a = std::pow(R,i0);
         reciprocal_fn.cuts.push_back(a);
     }  
 
diff --git a/src/2geom/sweeper.h b/src/2geom/sweeper.h
index 1cbce52b0..8c4e182a6 100644
--- a/src/2geom/sweeper.h
+++ b/src/2geom/sweeper.h
@@ -41,6 +41,9 @@
 
 namespace Geom {
 
+/** @brief Sweep traits class for interval bounds.
+ * @relates Sweeper
+ * @ingroup Utilities */
 struct IntervalSweepTraits {
     typedef Interval Bound;
     typedef std::less<Coord> Compare;
@@ -48,12 +51,22 @@ struct IntervalSweepTraits {
     inline static Coord exit_value(Bound const &b) { return b.max(); }
 };
 
-template <Dim2 d>
+/** @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(); }
+    inline static Coord entry_value(Bound const &b) { return b[D].min(); }
+    inline static Coord exit_value(Bound const &b) { return b[D].max(); }
+};
+
+template <typename T>
+struct BoundsFast {
+    Rect operator()(T const &item) const {
+        return item.boundsFast();
+    }
 };
 
 /** @brief Generic sweepline algorithm.
@@ -66,18 +79,30 @@ struct RectSweepTraits {
  *
  * To use this, create a derived class and reimplement the _enter()
  * and/or _leave() virtual functions, insert all the objects,
- * and finally call process(). You can specify the bound type
+ * 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.
  *
  * Look in path.cpp for example usage.
+ *
+ * @tparam Item The type of items to sweep
+ * @tparam SweepTraits Traits class that defines the items' bounds,
+ *   how to interpret them and how to sort the events
+ * @ingroup Utilities
  */
 template <typename Item, typename SweepTraits = IntervalSweepTraits>
 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),
@@ -85,6 +110,11 @@ public:
         _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) {
@@ -105,6 +135,8 @@ public:
 
         typename SweepTraits::Compare cmp;
 
+        // 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);
diff --git a/src/2geom/transforms.h b/src/2geom/transforms.h
index 7c03c5226..de4e6871f 100644
--- a/src/2geom/transforms.h
+++ b/src/2geom/transforms.h
@@ -40,6 +40,7 @@
 #include <2geom/forward.h>
 #include <2geom/affine.h>
 #include <2geom/angle.h>
+#include <boost/concept/assert.hpp>
 
 namespace Geom {
 
@@ -95,6 +96,7 @@ public:
  * @ingroup Transforms */
 template <typename T>
 T pow(T const &t, int n) {
+    BOOST_CONCEPT_ASSERT((TransformConcept<T>));
     if (n == 0) return T::identity();
     T result(T::identity());
     T x(n < 0 ? t.inverse() : t);