Update to 2Geom revision 2396

(bzr r14059.2.16)

70 files changed, 2203 insertions, 1669 deletions

diff --git a/src/2geom/CMakeLists.txt b/src/2geom/CMakeLists.txt
index eb25074ef..d0c196f56 100644
--- a/src/2geom/CMakeLists.txt
+++ b/src/2geom/CMakeLists.txt
@@ -7,7 +7,6 @@ set(2geom_SRC
 	bezier-curve.cpp
 	bezier-utils.cpp
 	cairo-path-sink.cpp
-	circle-circle.cpp
 	circle.cpp
 	# conic_section_clipper_impl.cpp
 	# conicsec.cpp
@@ -18,6 +17,7 @@ set(2geom_SRC
 	d2-sbasis.cpp
 	ellipse.cpp
 	elliptical-arc.cpp
+	elliptical-arc-from-sbasis.cpp
 	geom.cpp
 	intersection-graph.cpp
 	line.cpp
@@ -29,8 +29,7 @@ set(2geom_SRC
 	pathvector.cpp
 	piecewise.cpp
 	point.cpp
-	poly.cpp
-	quadtree.cpp
+	polynomial.cpp
 	rect.cpp
 	# recursive-bezier-intersection.cpp
 	sbasis-2d.cpp
@@ -43,8 +42,8 @@ set(2geom_SRC
 	solve-bezier.cpp
 	solve-bezier-one-d.cpp
 	solve-bezier-parametric.cpp
-	svg-elliptical-arc.cpp
 	svg-path-parser.cpp
+	svg-path-writer.cpp
 	sweep.cpp
 	toposweep.cpp
 	transforms.cpp
@@ -101,8 +100,7 @@ set(2geom_SRC
 	piecewise.h
 	point-ops.h
 	point.h
-	poly.h
-	quadtree.h
+	polynomial.h
 	ray.h
 	rect.h
 	region.h
@@ -115,8 +113,8 @@ set(2geom_SRC
 	sbasis.h
 	shape.h
 	solver.h
-	svg-elliptical-arc.h
 	svg-path-parser.h
+	svg-path-writer.h
 	sweep.h
 	toposweep.h
 	transforms.h
diff --git a/src/2geom/Makefile_insert b/src/2geom/Makefile_insert
index aeb0b20a1..01f48ab39 100644
--- a/src/2geom/Makefile_insert
+++ b/src/2geom/Makefile_insert
@@ -23,7 +23,6 @@
 	2geom/cairo-path-sink.cpp \
 	2geom/cairo-path-sink.h \
 	2geom/choose.h \
-	2geom/circle-circle.cpp \
 	2geom/circle.cpp \
 	2geom/circle.h \
 	2geom/circulator.h \
@@ -51,6 +50,7 @@
 	2geom/ellipse.h \
 	2geom/elliptical-arc.cpp \
 	2geom/elliptical-arc.h \
+	2geom/elliptical-arc-from-sbasis.cpp \
 	2geom/exception.h \
 	2geom/forward.h \
 	2geom/generic-interval.h \
@@ -84,10 +84,8 @@
 	2geom/piecewise.h \
 	2geom/point.cpp \
 	2geom/point.h \
-	2geom/poly.cpp \
-	2geom/poly.h \
-	2geom/quadtree.cpp \
-	2geom/quadtree.h \
+	2geom/polynomial.cpp \
+	2geom/polynomial.h \
 	2geom/ray.h \
 	2geom/rect.cpp \
 	2geom/rect.h \
@@ -110,14 +108,13 @@
 	2geom/solve-bezier-one-d.cpp \
 	2geom/solve-bezier-parametric.cpp \
 	2geom/solver.h \
-	2geom/svg-elliptical-arc.cpp \
-	2geom/svg-elliptical-arc.h \
 	2geom/svg-path-parser.cpp \
 	2geom/svg-path-parser.h \
 	2geom/svg-path-writer.cpp \
 	2geom/svg-path-writer.h \
-	2geom/sweep.cpp \
-	2geom/sweep.h \
+	2geom/sweep-bounds.cpp \
+	2geom/sweep-bounds.h \
+	2geom/sweeper.h \
 	2geom/toposweep.cpp \
 	2geom/toposweep.h \
 	2geom/transforms.cpp \
diff --git a/src/2geom/angle.h b/src/2geom/angle.h
index 77ecd5eb4..9ece3b9fe 100644
--- a/src/2geom/angle.h
+++ b/src/2geom/angle.h
@@ -63,11 +63,13 @@ namespace Geom {
  */
 class Angle
     : boost::additive< Angle
+    , boost::additive< Angle, Coord
     , boost::equality_comparable< Angle
-      > >
+    , boost::equality_comparable< Angle, Coord
+      > > > >
 {
 public:
-    Angle() : _angle(0) {} //added default constructor because of cython
+    Angle() : _angle(0) {}
     Angle(Coord v) : _angle(v) { _normalize(); } // this can be called implicitly
     explicit Angle(Point const &p) : _angle(atan2(p)) { _normalize(); }
     Angle(Point const &a, Point const &b) : _angle(angle_between(a, b)) { _normalize(); }
@@ -82,9 +84,20 @@ public:
         _normalize();
         return *this;
     }
+    Angle &operator+=(Coord a) {
+        *this += Angle(a);
+        return *this;
+    }
+    Angle &operator-=(Coord a) {
+        *this -= Angle(a);
+        return *this;
+    }
     bool operator==(Angle const &o) const {
         return _angle == o._angle;
     }
+    bool operator==(Coord c) const {
+        return _angle == Angle(c)._angle;
+    }
 
     /** @brief Get the angle as radians.
      * @return Number in range \f$[-\pi, \pi)\f$. */
@@ -136,12 +149,22 @@ public:
 private:
 
     void _normalize() {
-        _angle -= floor(_angle * (1.0/(2*M_PI))) * 2*M_PI;
+        _angle = std::fmod(_angle, 2*M_PI);
+        if (_angle < 0) _angle += 2*M_PI;
+        //_angle -= floor(_angle * (1.0/(2*M_PI))) * 2*M_PI;
     }
     Coord _angle; // this is always in [0, 2pi)
     friend class AngleInterval;
 };
 
+inline Angle distance(Angle const &a, Angle const &b) {
+    // the distance cannot be larger than M_PI.
+    Coord ac = a.radians0();
+    Coord bc = b.radians0();
+    Coord d = fabs(ac - bc);
+    return Angle(d > M_PI ? 2*M_PI - d : d);
+}
+
 /** @brief Directed angular interval.
  *
  * Wrapper for directed angles with defined start and end values. Useful e.g. for representing
@@ -149,7 +172,11 @@ private:
  * in the interval (it is a closed interval). Angular intervals can also be interptered
  * as functions \f$f: [0, 1] \to [-\pi, \pi)\f$, which return the start angle for 0,
  * the end angle for 1, and interpolate linearly for other values. Note that such functions
- * are not continuous if the interval contains the zero angle.
+ * 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.
@@ -158,39 +185,58 @@ private:
  */
 class AngleInterval {
 public:
+    /** @brief Create an angular interval.
+     * @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(double s, double e, bool cw = false)
         : _start_angle(s), _end_angle(e), _sweep(cw)
     {}
-    /** @brief Get the angular coordinate of the interval's initial point
-     * @return Angle in range \f$[0,2\pi)\f$ corresponding to the start of arc */
+
+    /// Get the start angle.
     Angle const &initialAngle() const { return _start_angle; }
-    /** @brief Get the angular coordinate of the interval's final point
-     * @return Angle in range \f$[0,2\pi)\f$ corresponding to the end of arc */
+    /// Get the end angle.
     Angle const &finalAngle() const { return _end_angle; }
+    /// Check whether the interval contains only a single angle.
     bool isDegenerate() const { return initialAngle() == finalAngle(); }
-    /** @brief Get an angle corresponding to the specified time value. */
+
+    /// Get an angle corresponding to the specified time value.
     Angle angleAt(Coord t) const {
         Coord span = extent();
         Angle ret = _start_angle.radians0() + span * (_sweep ? t : -t);
         return ret;
     }
     Angle operator()(Coord t) const { return angleAt(t); }
-#if 0
-    /** @brief Find an angle nearest to the specified time value.
-     * @param a Query angle
-     * @return If the interval contains the query angle, a number from the range [0, 1]
-     *         such that a = angleAt(t); otherwise, 0 or 1, depending on which extreme
-     *         angle is nearer. */
-    Coord nearestAngle(Angle const &a) const {
-        Coord dist = distance(_start_angle, a, _sweep).radians0();
-        Coord span = distance(_start_angle, _end_angle, _sweep).radians0();
-        if (dist <= span) return dist / span;
-        else return distance_abs(_start_angle, a).radians0() > distance_abs(_end_angle, a).radians0() ? 1.0 : 0.0;
+
+    /** @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 {
+        Coord ex = extent();
+        Coord outex = 2*M_PI - ex;
+        if (_sweep) {
+            Angle midout = _start_angle - outex / 2;
+            Angle acmp = a - midout, scmp = _start_angle - midout;
+            if (acmp.radians0() >= scmp.radians0()) {
+                return (a - _start_angle).radians0() / ex;
+            } else {
+                return -(_start_angle - a).radians0() / ex;
+            }
+        } else {
+            Angle midout = _start_angle + outex / 2;
+            Angle acmp = a - midout, scmp = _start_angle - midout;
+            if (acmp.radians0() <= scmp.radians0()) {
+                return (_start_angle - a).radians0() / ex;
+            } else {
+                return -(a - _start_angle).radians0() / ex;
+            }
+        }
     }
-#endif
+
     /** @brief Check whether the interval includes the given angle. */
     bool contains(Angle const &a) const {
         Coord s = _start_angle.radians0();
@@ -205,12 +251,22 @@ public:
         }
     }
     /** @brief Extent of the angle interval.
-     * @return Extent in range \f$[0, 2\pi)\f$ */
+     * Equivalent to the absolute value of the sweep angle.
+     * @return Extent in range \f$[0, 2\pi)\f$. */
     Coord extent() const {
-        Coord d = _end_angle - _start_angle;
-        if (!_sweep) d = -d;
-        if (d < 0) d += 2*M_PI;
-        return d;
+        return _sweep
+            ? (_end_angle - _start_angle).radians0()
+            : (_start_angle - _end_angle).radians0();
+    }
+    /** @brief Get the sweep angle of the interval.
+     * This is the value you need to add to the initial angle to get the final angle.
+     * It is positive when sweep is true. Denoted as \f$\Delta\theta\f$ in the SVG
+     * elliptical arc implementation notes. */
+    Coord sweepAngle() const {
+        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() {}
@@ -308,6 +364,10 @@ bool arc_contains (double a, double sa, double ia, double ea)
 
 } // end namespace Geom
 
+namespace std {
+template <> class iterator_traits<Geom::Angle> {};
+}
+
 #endif // LIB2GEOM_SEEN_ANGLE_H
 
 /*
diff --git a/src/2geom/bezier-curve.cpp b/src/2geom/bezier-curve.cpp
index b81041f29..17221264b 100644
--- a/src/2geom/bezier-curve.cpp
+++ b/src/2geom/bezier-curve.cpp
@@ -143,9 +143,15 @@ BezierCurve::intersect(Curve const &other, Coord eps) const
 {
     std::vector<CurveIntersection> result;
 
-    // optimization for the common case of no intersections
-    if (!boundsFast().intersects(other.boundsFast())) return result;
+    // in case we encounter an order-1 curve created from a vector
+    // or a degenerate elliptical arc
+    if (isLineSegment()) {
+        LineSegment ls(initialPoint(), finalPoint());
+        result = ls.intersect(other);
+        return result;
+    }
 
+    // here we are sure that this curve is at least a quadratic Bezier
     BezierCurve const *bez = dynamic_cast<BezierCurve const *>(&other);
     if (bez) {
         std::vector<std::pair<double, double> > xs;
@@ -154,10 +160,12 @@ BezierCurve::intersect(Curve const &other, Coord eps) const
             CurveIntersection x(*this, other, xs[i].first, xs[i].second);
             result.push_back(x);
         }
-    } else {
-        THROW_NOTIMPLEMENTED();
+        return result;
     }
 
+    // pass other intersection types to the other curve
+    result = other.intersect(*this, eps);
+    transpose_in_place(result);
     return result;
 }
 
@@ -253,6 +261,26 @@ Coord BezierCurveN<1>::nearestTime(Point const& p, Coord from, Coord to) const
 }
 
 template <>
+std::vector<CurveIntersection> BezierCurveN<1>::intersect(Curve const &other, Coord eps) const
+{
+    std::vector<CurveIntersection> result;
+
+    // only handle intersections with other LineSegments here
+    if (other.isLineSegment()) {
+        Line this_line(initialPoint(), finalPoint());
+        Line other_line(other.initialPoint(), other.finalPoint());
+        result = this_line.intersect(other_line);
+        filter_line_segment_intersections(result, true, true);
+        return result;
+    }
+
+    // pass all other types to the other curve
+    result = other.intersect(*this, eps);
+    transpose_in_place(result);
+    return result;
+}
+
+template <>
 void BezierCurveN<1>::feed(PathSink &sink, bool moveto_initial) const
 {
     if (moveto_initial) {
diff --git a/src/2geom/bezier-curve.h b/src/2geom/bezier-curve.h
index 9b08466f8..636a263a7 100644
--- a/src/2geom/bezier-curve.h
+++ b/src/2geom/bezier-curve.h
@@ -69,6 +69,7 @@ public:
     /** @brief Get the control points.
      * @return Vector with order() + 1 control points. */
     std::vector<Point> controlPoints() const { return bezier_points(inner); }
+    D2<Bezier> const &fragment() const { return inner; }
 
     /** @brief Modify a control point.
      * @param ix The zero-based index of the point to modify. Note that the caller is responsible for checking that this value is <= order().
@@ -104,6 +105,7 @@ public:
     virtual Point initialPoint() const { return inner.at0(); }
     virtual Point finalPoint() const { return inner.at1(); }
     virtual bool isDegenerate() const { return inner.isConstant(0); }
+    virtual bool isLineSegment() const { return size() == 2; }
     virtual void setInitial(Point const &v) { setPoint(0, v); }
     virtual void setFinal(Point const &v) { setPoint(order(), v); }
     virtual Rect boundsFast() const { return *bounds_fast(inner); }
@@ -145,7 +147,7 @@ public:
     virtual Coord nearestTime(Point const &p, Coord from = 0, Coord to = 1) const;
     virtual Coord length(Coord tolerance) const;
     virtual std::vector<CurveIntersection> intersect(Curve const &other, Coord eps = EPSILON) const;
-    virtual Point pointAt(Coord t) const { return inner.valueAt(t); }
+    virtual Point pointAt(Coord t) const { return inner.pointAt(t); }
     virtual std::vector<Point> pointAndDerivatives(Coord t, unsigned n) const {
         return inner.valueAndDerivatives(t, n);
     }
@@ -229,6 +231,10 @@ public:
                    BezierCurveN(sx.second, sy.second));
     }
 
+    virtual bool isLineSegment() const {
+        return size() == 2;
+    }
+
     virtual Curve *duplicate() const {
         return new BezierCurveN(*this);
     }
@@ -251,6 +257,10 @@ public:
     virtual Coord nearestTime(Point const &p, Coord from = 0, Coord to = 1) const {
         return BezierCurve::nearestTime(p, from, to);
     }
+    virtual std::vector<CurveIntersection> intersect(Curve const &other, Coord eps = EPSILON) const {
+        // call super. this is implemented only to allow specializations
+        return BezierCurve::intersect(other, eps);
+    }
     virtual void feed(PathSink &sink, bool moveto_initial) const {
         // call super. this is implemented only to allow specializations
         BezierCurve::feed(sink, moveto_initial);
@@ -281,8 +291,10 @@ Curve *BezierCurveN<degree>::derivative() const {
 }
 
 // optimized specializations
+template <> inline bool BezierCurveN<1>::isLineSegment() const { return true; }
 template <> Curve *BezierCurveN<1>::derivative() const;
 template <> Coord BezierCurveN<1>::nearestTime(Point const &, Coord, Coord) const;
+template <> std::vector<CurveIntersection> BezierCurveN<1>::intersect(Curve const &, Coord) const;
 template <> void BezierCurveN<1>::feed(PathSink &sink, bool moveto_initial) const;
 template <> void BezierCurveN<2>::feed(PathSink &sink, bool moveto_initial) const;
 template <> void BezierCurveN<3>::feed(PathSink &sink, bool moveto_initial) const;
diff --git a/src/2geom/bezier.cpp b/src/2geom/bezier.cpp
index 45496b75c..0c9d12c3b 100644
--- a/src/2geom/bezier.cpp
+++ b/src/2geom/bezier.cpp
@@ -210,7 +210,7 @@ Bezier &Bezier::operator-=(Bezier const &other)
 
 
 
-Bezier multiply(Bezier const &f, Bezier const &g)
+Bezier operator*(Bezier const &f, Bezier const &g)
 {
     unsigned m = f.order();
     unsigned n = g.order();
diff --git a/src/2geom/bezier.h b/src/2geom/bezier.h
index be7df1a6b..c41c2b3a7 100644
--- a/src/2geom/bezier.h
+++ b/src/2geom/bezier.h
@@ -308,7 +308,11 @@ public:
 
 void bezier_to_sbasis (SBasis &sb, Bezier const &bz);
 
-Bezier multiply(Bezier const &f, Bezier const &g);
+Bezier operator*(Bezier const &f, Bezier const &g);
+inline Bezier multiply(Bezier const &f, Bezier const &g) {
+    Bezier result = f * g;
+    return result;
+}
 
 inline Bezier reverse(const Bezier & a) {
     Bezier result = Bezier(Bezier::Order(a));
diff --git a/src/2geom/cairo-path-sink.cpp b/src/2geom/cairo-path-sink.cpp
index f327bf04d..244a08ba4 100644
--- a/src/2geom/cairo-path-sink.cpp
+++ b/src/2geom/cairo-path-sink.cpp
@@ -31,7 +31,7 @@
 
 #include <cairo.h>
 #include <2geom/cairo-path-sink.h>
-#include <2geom/svg-elliptical-arc.h>
+#include <2geom/elliptical-arc.h>
 
 namespace Geom {
 
@@ -71,7 +71,7 @@ void CairoPathSink::quadTo(Point const &p1, Point const &p2)
 void CairoPathSink::arcTo(double rx, double ry, double angle,
                           bool large_arc, bool sweep, Point const &p)
 {
-    SVGEllipticalArc arc(_current_point, rx, ry, angle, large_arc, sweep, p);
+    EllipticalArc arc(_current_point, rx, ry, angle, large_arc, sweep, p);
     // Cairo only does circular arcs.
     // To do elliptical arcs, we must use a temporary transform.
     Affine uct = arc.unitCircleTransform();
diff --git a/src/2geom/circle-circle.cpp b/src/2geom/circle-circle.cpp
deleted file mode 100644
index 134fa33a2..000000000
--- a/src/2geom/circle-circle.cpp
+++ /dev/null
@@ -1,139 +0,0 @@
-/* circle_circle_intersection() *
- * Determine the points where 2 circles in a common plane intersect.
- *
- * int circle_circle_intersection(
- *                                // center and radius of 1st circle
- *                                double x0, double y0, double r0,
- *                                // center and radius of 2nd circle
- *                                double x1, double y1, double r1,
- *                                // 1st intersection point
- *                                double *xi, double *yi,              
- *                                // 2nd intersection point
- *                                double *xi_prime, double *yi_prime)
- *
- * This is a public domain work. 3/26/2005 Tim Voght
- * Ported to lib2geom, 2006 Nathan Hurst
- *
- * Copyright 2006 Nathan Hurst <njh@mail.csse.monash.edu.au>
- *
- * This library is free software; you can redistribute it and/or
- * modify it either under the terms of the GNU Lesser General Public
- * License version 2.1 as published by the Free Software Foundation
- * (the "LGPL") or, at your option, under the terms of the Mozilla
- * Public License Version 1.1 (the "MPL"). If you do not alter this
- * notice, a recipient may use your version of this file under either
- * the MPL or the LGPL.
- *
- * You should have received a copy of the LGPL along with this library
- * in the file COPYING-LGPL-2.1; if not, write to the Free Software
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
- * You should have received a copy of the MPL along with this library
- * in the file COPYING-MPL-1.1
- *
- * The contents of this file are subject to the Mozilla Public License
- * Version 1.1 (the "License"); you may not use this file except in
- * compliance with the License. You may obtain a copy of the License at
- * http://www.mozilla.org/MPL/
- *
- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
- * OF ANY KIND, either express or implied. See the LGPL or the MPL for
- * the specific language governing rights and limitations.
- * 
- */
-#include <stdio.h>
-#include <math.h>
-#include <2geom/point.h>
-
-namespace Geom{
-
-int circle_circle_intersection(Point X0, double r0,
-                               Point X1, double r1,
-                               Point & p0, Point & p1)
-{
-  /* dx and dy are the vertical and horizontal distances between
-   * the circle centers.
-   */
-  Point D = X1 - X0;
-
-  /* Determine the straight-line distance between the centers. */
-  double d = L2(D);
-
-  /* Check for solvability. */
-  if (d > (r0 + r1))
-  {
-    /* no solution. circles do not intersect. */
-    return 0;
-  }
-  if (d <= fabs(r0 - r1))
-  {
-    /* no solution. one circle is contained in the other */
-    return 1;
-  }
-  
-  /* 'point 2' is the point where the line through the circle
-   * intersection points crosses the line between the circle
-   * centers.  
-   */
-
-  /* Determine the distance from point 0 to point 2. */
-  double a = ((r0*r0) - (r1*r1) + (d*d)) / (2.0 * d) ;
-
-  /* Determine the coordinates of point 2. */
-  Point p2 = X0 + D * (a/d);
-
-  /* Determine the distance from point 2 to either of the
-   * intersection points.
-   */
-  double h = sqrt((r0*r0) - (a*a));
-
-  /* Now determine the offsets of the intersection points from
-   * point 2.
-   */
-  Point r = (h/d)*rot90(D);
-
-  /* Determine the absolute intersection points. */
-  p0 = p2 + r;
-  p1 = p2 - r;
-
-  return 2;
-}
-
-};
-
-
-#ifdef TEST
-
-void run_test(double x0, double y0, double r0,
-              double x1, double y1, double r1)
-{
-  printf("x0=%F, y0=%F, r0=%F, x1=%F, y1=%F, r1=%F :\n",
-          x0, y0, r0, x1, y1, r1);
-  Geom::Point p0, p1;
-  Geom::circle_circle_intersection(Geom::Point(x0, y0), r0, 
-				   Geom::Point(x1, y1), r1,
-				   p0, p1);
-  printf("  x3=%F, y3=%F, x3_prime=%F, y3_prime=%F\n",
-            p0[0], p0[1], p1[0], p1[1]);
-}
-
-int main(void)
-{
-  /* Add more! */    
-  run_test(-1.0, -1.0, 1.5, 1.0, 1.0, 2.0);
-  run_test(1.0, -1.0, 1.5, -1.0, 1.0, 2.0);
-  run_test(-1.0, 1.0, 1.5, 1.0, -1.0, 2.0);
-  run_test(1.0, 1.0, 1.5, -1.0, -1.0, 2.0);
-  exit(0);
-}
-#endif
-
-/*
-  Local Variables:
-  mode:c++
-  c-file-style:"stroustrup"
-  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
-  indent-tabs-mode:nil
-  fill-column:99
-  End:
-*/
-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
diff --git a/src/2geom/circle.cpp b/src/2geom/circle.cpp
index 0b1dddc8e..ec59bbe4a 100644
--- a/src/2geom/circle.cpp
+++ b/src/2geom/circle.cpp
@@ -33,12 +33,19 @@
 
 #include <2geom/circle.h>
 #include <2geom/ellipse.h>
-#include <2geom/svg-elliptical-arc.h>
+#include <2geom/elliptical-arc.h>
 #include <2geom/numeric/fitting-tool.h>
 #include <2geom/numeric/fitting-model.h>
 
 namespace Geom {
 
+Rect Circle::boundsFast() const
+{
+    Point rr(_radius, _radius);
+    Rect bbox(_center - rr, _center + rr);
+    return bbox;
+}
+
 void Circle::setCoefficients(Coord A, Coord B, Coord C, Coord D)
 {
     if (A == 0) {
@@ -62,44 +69,197 @@ void Circle::setCoefficients(Coord A, Coord B, Coord C, Coord D)
     _radius = std::sqrt(r2);
 }
 
+void Circle::coefficients(Coord &A, Coord &B, Coord &C, Coord &D) const
+{
+    A = 1;
+    B = -2 * _center[X];
+    C = -2 * _center[Y];
+    D = _center[X] * _center[X] + _center[Y] * _center[Y] - _radius * _radius;
+}
 
-void Circle::fit(std::vector<Point> const& points)
+std::vector<Coord> Circle::coefficients() const
 {
-    size_t sz = points.size();
-    if (sz < 2) {
-        THROW_RANGEERROR("fitting error: too few points passed");
+    std::vector<Coord> c(4);
+    coefficients(c[0], c[1], c[2], c[3]);
+    return c;
+}
+
+
+Zoom Circle::unitCircleTransform() const
+{
+    Zoom ret(_radius, _center / _radius);
+    return ret;
+}
+
+Zoom Circle::inverseUnitCircleTransform() const
+{
+    if (_radius == 0) {
+        THROW_RANGEERROR("degenerate circle does not have an inverse unit circle transform");
     }
-    if (sz == 2) {
-        _center = points[0] * 0.5 + points[1] * 0.5;
-        _radius = distance(points[0], points[1]) / 2;
-        return;
+
+    Zoom ret(1/_radius, Translate(-_center));
+    return ret;
+}
+
+
+Point Circle::pointAt(Coord t) const {
+    return _center + Point::polar(t) * _radius;
+}
+
+Coord Circle::valueAt(Coord t, Dim2 d) const {
+    Coord delta = (d == X ? std::cos(t) : std::sin(t));
+    return _center[d] + delta * _radius;
+}
+
+Coord Circle::timeAt(Point const &p) const {
+    return atan2(p - _center);
+}
+
+Coord Circle::nearestTime(Point const &p) const {
+    if (_center == p) return 0;
+    return timeAt(p);
+}
+
+bool Circle::contains(Rect const &r) const
+{
+    for (unsigned i = 0; i < 4; ++i) {
+        if (!contains(r.corner(i))) return false;
     }
+    return true;
+}
 
-    NL::LFMCircle model;
-    NL::least_squeares_fitter<NL::LFMCircle> fitter(model, sz);
+bool Circle::contains(Circle const &other) const
+{
+    Coord cdist = distance(_center, other._center);
+    Coord rdist = fabs(_radius - other._radius);
+    return cdist <= rdist;
+}
 
-    for (size_t i = 0; i < sz; ++i) {
-        fitter.append(points[i]);
+bool Circle::intersects(Line const &l) const
+{
+    // http://mathworld.wolfram.com/Circle-LineIntersection.html
+    Coord dr = l.versor().length();
+    Coord r = _radius;
+    Coord D = cross(l.initialPoint(), l.finalPoint());
+    Coord delta = r*r * dr*dr - D*D;
+    if (delta >= 0) return true;
+    return false;
+}
+
+bool Circle::intersects(Circle const &other) const
+{
+    Coord cdist = distance(_center, other._center);
+    Coord rsum = _radius + other._radius;
+    return cdist <= rsum;
+}
+
+
+std::vector<ShapeIntersection> Circle::intersect(Line const &l) const
+{
+    // http://mathworld.wolfram.com/Circle-LineIntersection.html
+    Coord dr = l.versor().length();
+    Coord dx = l.versor().x();
+    Coord dy = l.versor().y();
+    Coord D = cross(l.initialPoint() - _center, l.finalPoint() - _center);
+    Coord delta = _radius*_radius * dr*dr - D*D;
+
+    std::vector<ShapeIntersection> result;
+    if (delta < 0) return result;
+    if (delta == 0) {
+        Coord ix = (D*dy) / (dr*dr);
+        Coord iy = (-D*dx) / (dr*dr);
+        Point ip(ix, iy); ip += _center;
+        result.push_back(ShapeIntersection(timeAt(ip), l.timeAt(ip), ip));
+        return result;
     }
-    fitter.update();
 
-    NL::Vector z(sz, 0.0);
-    model.instance(*this, fitter.result(z));
+    Coord sqrt_delta = std::sqrt(delta);
+    Coord signmod = dy < 0 ? -1 : 1;
+
+    Coord i1x = (D*dy + signmod * dx * sqrt_delta) / (dr*dr);
+    Coord i1y = (-D*dx + fabs(dy) * sqrt_delta) / (dr*dr);
+    Point i1p(i1x, i1y); i1p += _center;
+
+    Coord i2x = (D*dy - signmod * dx * sqrt_delta) / (dr*dr);
+    Coord i2y = (-D*dx - fabs(dy) * sqrt_delta) / (dr*dr);
+    Point i2p(i2x, i2y); i2p += _center;
+
+    result.push_back(ShapeIntersection(timeAt(i1p), l.timeAt(i1p), i1p));
+    result.push_back(ShapeIntersection(timeAt(i2p), l.timeAt(i2p), i2p));
+    return result;
+}
+
+std::vector<ShapeIntersection> Circle::intersect(LineSegment const &l) const
+{
+    std::vector<ShapeIntersection> result = intersect(Line(l));
+    filter_line_segment_intersections(result);
+    return result;
+}
+
+std::vector<ShapeIntersection> Circle::intersect(Circle const &other) const
+{
+    std::vector<ShapeIntersection> result;
+
+    if (*this == other) {
+        THROW_INFINITESOLUTIONS();
+    }
+    if (contains(other)) return result;
+    if (!intersects(other)) return result;
+
+    // See e.g. http://mathworld.wolfram.com/Circle-CircleIntersection.html
+    // Basically, we figure out where is the third point of a triangle
+    // with two points in the centers and with edge lengths equal to radii
+    Point cv = other._center - _center;
+    Coord d = cv.length();
+    Coord R = radius(), r = other.radius();
+
+    if (d == R + r) {
+        Point px = lerp(R / d, _center, other._center);
+        Coord T = timeAt(px), t = other.timeAt(px);
+        result.push_back(ShapeIntersection(T, t, px));
+        return result;
+    }
+
+    // q is the distance along the line between centers to the perpendicular line
+    // that goes through both intersections.
+    Coord q = (d*d - r*r + R*R) / (2*d);
+    Point qp = lerp(q/d, _center, other._center);
+
+    // The triangle given by the points:
+    // _center, qp, intersection
+    // is a right triangle. Determine the distance between qp and intersection
+    // using the Pythagorean theorem.
+    Coord h = std::sqrt(R*R - q*q);
+    Point qd = (h/d) * cv.cw();
+
+    // now compute the intersection points
+    Point x1 = qp + qd;
+    Point x2 = qp - qd;
+
+    result.push_back(ShapeIntersection(timeAt(x1), other.timeAt(x1), x1));
+    result.push_back(ShapeIntersection(timeAt(x2), other.timeAt(x2), x2));
+    return result;
 }
 
 /**
  @param inner a point whose angle with the circle center is inside the angle that the arc spans
  */
 EllipticalArc *
-Circle::arc(Point const& initial, Point const& inner, Point const& final,
-             bool svg_compliant)
+Circle::arc(Point const& initial, Point const& inner, Point const& final) const
 {
     // TODO native implementation!
     Ellipse e(_center[X], _center[Y], _radius, _radius, 0);
-    return e.arc(initial, inner, final, svg_compliant);
+    return e.arc(initial, inner, final);
+}
+
+bool Circle::operator==(Circle const &other) const
+{
+    if (_center != other._center) return false;
+    if (_radius != other._radius) return false;
+    return true;
 }
 
-D2<SBasis> Circle::toSBasis()
+D2<SBasis> Circle::toSBasis() const
 {
     D2<SBasis> B;
     Linear bo = Linear(0, 2 * M_PI);
@@ -112,22 +272,52 @@ D2<SBasis> Circle::toSBasis()
     return B;
 }
 
-void
-Circle::getPath(PathVector &path_out) {
-    Path pb;
 
-    D2<SBasis> B = toSBasis();
+void Circle::fit(std::vector<Point> const& points)
+{
+    size_t sz = points.size();
+    if (sz < 2) {
+        THROW_RANGEERROR("fitting error: too few points passed");
+    }
+    if (sz == 2) {
+        _center = points[0] * 0.5 + points[1] * 0.5;
+        _radius = distance(points[0], points[1]) / 2;
+        return;
+    }
+
+    NL::LFMCircle model;
+    NL::least_squeares_fitter<NL::LFMCircle> fitter(model, sz);
 
-    pb.append(SBasisCurve(B));
+    for (size_t i = 0; i < sz; ++i) {
+        fitter.append(points[i]);
+    }
+    fitter.update();
 
-    path_out.push_back(pb);
+    NL::Vector z(sz, 0.0);
+    model.instance(*this, fitter.result(z));
 }
 
 
-}  // end namespace Geom
-
+bool are_near(Circle const &a, Circle const &b, Coord eps)
+{
+    // to check whether no point on a is further than eps from b,
+    // we check two things:
+    // 1. if radii differ by more than eps, there is definitely a point that fails
+    // 2. if they differ by less, we check the centers. They have to be closer
+    //    together if the radius differs, since the maximum distance will be
+    //    equal to sum of radius difference and distance between centers.
+    if (!are_near(a.radius(), b.radius(), eps)) return false;
+    Coord adjusted_eps = eps - fabs(a.radius() - b.radius());
+    return are_near(a.center(), b.center(), adjusted_eps);
+}
 
+std::ostream &operator<<(std::ostream &out, Circle const &c)
+{
+    out << "Circle(" << c.center() << ", " << format_coord_nice(c.radius()) << ")";
+    return out;
+}
 
+}  // end namespace Geom
 
 /*
   Local Variables:
diff --git a/src/2geom/circle.h b/src/2geom/circle.h
index 3c2115b12..c42cb7f80 100644
--- a/src/2geom/circle.h
+++ b/src/2geom/circle.h
@@ -34,8 +34,10 @@
 #ifndef LIB2GEOM_SEEN_CIRCLE_H
 #define LIB2GEOM_SEEN_CIRCLE_H
 
-#include <2geom/point.h>
 #include <2geom/forward.h>
+#include <2geom/intersection.h>
+#include <2geom/point.h>
+#include <2geom/rect.h>
 #include <2geom/transforms.h>
 
 namespace Geom {
@@ -45,10 +47,11 @@ class EllipticalArc;
 /** @brief Set of all points at a fixed distance from the center
  * @ingroup Shapes */
 class Circle
-    : MultipliableNoncommutative< Circle, Translate
+    : boost::equality_comparable1< Circle
+    , MultipliableNoncommutative< Circle, Translate
     , MultipliableNoncommutative< Circle, Rotate
     , MultipliableNoncommutative< Circle, Zoom
-      > > >
+      > > > >
 {
     Point _center;
     Coord _radius;
@@ -66,28 +69,51 @@ public:
         setCoefficients(A, B, C, D);
     }
 
-    // build a circle by its implicit equation:
-    // Ax^2 + Ay^2 + Bx + Cy + D = 0
-    void setCoefficients(Coord A, Coord B, Coord C, Coord D);
-
-    /** @brief Fit the circle to the passed points using the least squares method.
-     * @param points Samples at the perimeter of the circle */
-    void fit(std::vector<Point> const &points);
-
-    EllipticalArc *
-    arc(Point const& initial, Point const& inner, Point const& final,
-        bool svg_compliant = true);
-
-    D2<SBasis> toSBasis();
-    void getPath(PathVector &path_out);
+    // Construct the unique circle passing through three points.
+    //Circle(Point const &a, Point const &b, Point const &c);
 
     Point center() const { return _center; }
     Coord center(Dim2 d) const { return _center[d]; }
     Coord radius() const { return _radius; }
+    Coord area() const { return M_PI * _radius * _radius; }
 
     void setCenter(Point const &p) { _center = p; }
     void setRadius(Coord c) { _radius = c; }
 
+    Rect boundsFast() const;
+    Rect boundsExact() const { return boundsFast(); }
+
+    Point pointAt(Coord t) const;
+    Coord valueAt(Coord t, Dim2 d) const;
+    Coord timeAt(Point const &p) const;
+    Coord nearestTime(Point const &p) const;
+
+    bool contains(Point const &p) const { return distance(p, _center) <= _radius; }
+    bool contains(Rect const &other) const;
+    bool contains(Circle const &other) const;
+
+    bool intersects(Line const &l) const;
+    bool intersects(LineSegment const &l) const;
+    bool intersects(Circle const &other) const;
+
+    std::vector<ShapeIntersection> intersect(Line const &other) const;
+    std::vector<ShapeIntersection> intersect(LineSegment const &other) const;
+    std::vector<ShapeIntersection> intersect(Circle const &other) const;
+
+    // build a circle by its implicit equation:
+    // Ax^2 + Ay^2 + Bx + Cy + D = 0
+    void setCoefficients(Coord A, Coord B, Coord C, Coord D);
+    void coefficients(Coord &A, Coord &B, Coord &C, Coord &D) const;
+    std::vector<Coord> coefficients() const;
+
+    Zoom unitCircleTransform() const;
+    Zoom inverseUnitCircleTransform() const;
+
+    EllipticalArc *
+    arc(Point const& initial, Point const& inner, Point const& final) const;
+
+    D2<SBasis> toSBasis() const;
+
     Circle &operator*=(Translate const &t) {
         _center *= t;
         return *this;
@@ -100,8 +126,20 @@ public:
         _radius *= z.scale();
         return *this;
     }
+
+    bool operator==(Circle const &other) const;
+
+    /** @brief Fit the circle to the passed points using the least squares method.
+     * @param points Samples at the perimeter of the circle */
+    void fit(std::vector<Point> const &points);
 };
 
+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);
+
 } // end namespace Geom
 
 #endif // LIB2GEOM_SEEN_CIRCLE_H
diff --git a/src/2geom/concepts.h b/src/2geom/concepts.h
index f571ddc60..7bc1bc0fd 100644
--- a/src/2geom/concepts.h
+++ b/src/2geom/concepts.h
@@ -131,6 +131,16 @@ template <typename T>
 inline T portion(const T& t, const Interval& i) { return portion(t, i.min(), i.max()); }
 
 template <typename T>
+struct EqualityComparableConcept {
+    T a, b;
+    bool bool_;
+    void constaints() {
+        bool_ = (a == b);
+        bool_ = (a != b);
+    }
+};
+
+template <typename T>
 struct NearConcept {
     T a, b;
     double tol;
diff --git a/src/2geom/crossing.h b/src/2geom/crossing.h
index be41fe893..425fa58f5 100644
--- a/src/2geom/crossing.h
+++ b/src/2geom/crossing.h
@@ -38,7 +38,7 @@
 
 #include <vector>
 #include <2geom/rect.h>
-#include <2geom/sweep.h>
+#include <2geom/sweep-bounds.h>
 #include <boost/optional/optional.hpp>
 #include <2geom/pathvector.h>
 
diff --git a/src/2geom/curve.h b/src/2geom/curve.h
index 893dc6bdb..7da0d17a0 100644
--- a/src/2geom/curve.h
+++ b/src/2geom/curve.h
@@ -94,6 +94,9 @@ public:
      *         no other points (its value set contains only one element). */
     virtual bool isDegenerate() const = 0;
 
+    /// Check whether the curve is a line segment.
+    virtual bool isLineSegment() const { return false; }
+
     /** @brief Get the interval of allowed time values.
      * @return \f$[0, 1]\f$ */
     virtual Interval timeRange() const {
@@ -294,13 +297,7 @@ public:
      * derivative could be found.
      * @param t Time value
      * @param n The maximum order of derivative to consider
-     * @return Unit tangent vector \f$\mathbf{v}(t)\f$
-     * @bug This method might currently break for the case of t being exactly 1. A workaround
-     *      is to reverse the curve and use the negated unit tangent at 0 like this:
-     * @code
-        Curve *c_reverse = c.reverse();
-        Point tangent = - c_reverse->unitTangentAt(0);
-        delete c_reverse; @endcode */
+     * @return Unit tangent vector \f$\mathbf{v}(t)\f$ */
     virtual Point unitTangentAt(Coord t, unsigned n = 3) const;
 
     /** @brief Convert the curve to a symmetric power basis polynomial.
diff --git a/src/2geom/curves.h b/src/2geom/curves.h
index 6022c71d8..46fb6d973 100644
--- a/src/2geom/curves.h
+++ b/src/2geom/curves.h
@@ -38,7 +38,6 @@
 #include <2geom/sbasis-curve.h>
 #include <2geom/bezier-curve.h>
 #include <2geom/elliptical-arc.h>
-#include <2geom/svg-elliptical-arc.h>
 
 #endif // LIB2GEOM_SEEN_CURVES_H
 
diff --git a/src/2geom/d2.h b/src/2geom/d2.h
index 714319f99..bf764539d 100644
--- a/src/2geom/d2.h
+++ b/src/2geom/d2.h
@@ -33,7 +33,7 @@
 #define LIB2GEOM_SEEN_D2_H
 
 #include <iterator>
-#include <boost/concept_check.hpp>
+#include <boost/concept/assert.hpp>
 #include <boost/iterator/transform_iterator.hpp>
 #include <2geom/point.h>
 #include <2geom/interval.h>
@@ -107,30 +107,36 @@ class D2{
     //IMPL: FragmentConcept
     typedef Point output_type;
     bool isZero(double eps=EPSILON) const {
-        boost::function_requires<FragmentConcept<T> >();
+        BOOST_CONCEPT_ASSERT((FragmentConcept<T>));
         return f[X].isZero(eps) && f[Y].isZero(eps);
     }
     bool isConstant(double eps=EPSILON) const {
-        boost::function_requires<FragmentConcept<T> >();
+        BOOST_CONCEPT_ASSERT((FragmentConcept<T>));
         return f[X].isConstant(eps) && f[Y].isConstant(eps);
     }
     bool isFinite() const {
-        boost::function_requires<FragmentConcept<T> >();
+        BOOST_CONCEPT_ASSERT((FragmentConcept<T>));
         return f[X].isFinite() && f[Y].isFinite();
     }
     Point at0() const { 
-        boost::function_requires<FragmentConcept<T> >();
+        BOOST_CONCEPT_ASSERT((FragmentConcept<T>));
         return Point(f[X].at0(), f[Y].at0());
     }
     Point at1() const {
-        boost::function_requires<FragmentConcept<T> >();
+        BOOST_CONCEPT_ASSERT((FragmentConcept<T>));
         return Point(f[X].at1(), f[Y].at1());
     }
+    Point pointAt(double t) const {
+        BOOST_CONCEPT_ASSERT((FragmentConcept<T>));
+        return (*this)(t);
+    }
     Point valueAt(double t) const {
-        boost::function_requires<FragmentConcept<T> >();
+        // TODO: remove this alias
+        BOOST_CONCEPT_ASSERT((FragmentConcept<T>));
         return (*this)(t);
     }
     std::vector<Point > valueAndDerivatives(double t, unsigned n) const {
+        BOOST_CONCEPT_ASSERT((FragmentConcept<T>));
         std::vector<Coord> x = f[X].valueAndDerivatives(t, n),
                            y = f[Y].valueAndDerivatives(t, n); // always returns a vector of size n+1
         std::vector<Point> res(n+1);
@@ -140,7 +146,7 @@ class D2{
         return res;
     }
     D2<SBasis> toSBasis() const {
-        boost::function_requires<FragmentConcept<T> >();
+        BOOST_CONCEPT_ASSERT((FragmentConcept<T>));
         return D2<SBasis>(f[X].toSBasis(), f[Y].toSBasis());
     }
 
@@ -149,33 +155,33 @@ class D2{
 };
 template <typename T>
 inline D2<T> reverse(const D2<T> &a) {
-    boost::function_requires<FragmentConcept<T> >();
+    BOOST_CONCEPT_ASSERT((FragmentConcept<T>));
     return D2<T>(reverse(a[X]), reverse(a[Y]));
 }
 
 template <typename T>
 inline D2<T> portion(const D2<T> &a, Coord f, Coord t) {
-    boost::function_requires<FragmentConcept<T> >();
+    BOOST_CONCEPT_ASSERT((FragmentConcept<T>));
     return D2<T>(portion(a[X], f, t), portion(a[Y], f, t));
 }
 
 template <typename T>
 inline D2<T> portion(const D2<T> &a, Interval i) {
-    boost::function_requires<FragmentConcept<T> >();
+    BOOST_CONCEPT_ASSERT((FragmentConcept<T>));
     return D2<T>(portion(a[X], i), portion(a[Y], i));
 }
 
-//IMPL: boost::EqualityComparableConcept
+//IMPL: EqualityComparableConcept
 template <typename T>
 inline bool
 operator==(D2<T> const &a, D2<T> const &b) {
-    boost::function_requires<boost::EqualityComparableConcept<T> >();
+    BOOST_CONCEPT_ASSERT((EqualityComparableConcept<T>));
     return a[0]==b[0] && a[1]==b[1];
 }
 template <typename T>
 inline bool
 operator!=(D2<T> const &a, D2<T> const &b) {
-    boost::function_requires<boost::EqualityComparableConcept<T> >();
+    BOOST_CONCEPT_ASSERT((EqualityComparableConcept<T>));
     return a[0]!=b[0] || a[1]!=b[1];
 }
 
@@ -183,7 +189,7 @@ operator!=(D2<T> const &a, D2<T> const &b) {
 template <typename T>
 inline bool
 are_near(D2<T> const &a, D2<T> const &b, double tol) {
-    boost::function_requires<NearConcept<T> >();
+    BOOST_CONCEPT_ASSERT((NearConcept<T>));
     return are_near(a[0], b[0], tol) && are_near(a[1], b[1], tol);
 }
 
@@ -191,7 +197,7 @@ are_near(D2<T> const &a, D2<T> const &b, double tol) {
 template <typename T>
 inline D2<T>
 operator+(D2<T> const &a, D2<T> const &b) {
-    boost::function_requires<AddableConcept<T> >();
+    BOOST_CONCEPT_ASSERT((AddableConcept<T>));
 
     D2<T> r;
     for(unsigned i = 0; i < 2; i++)
@@ -201,7 +207,7 @@ operator+(D2<T> const &a, D2<T> const &b) {
 template <typename T>
 inline D2<T>
 operator-(D2<T> const &a, D2<T> const &b) {
-    boost::function_requires<AddableConcept<T> >();
+    BOOST_CONCEPT_ASSERT((AddableConcept<T>));
 
     D2<T> r;
     for(unsigned i = 0; i < 2; i++)
@@ -211,7 +217,7 @@ operator-(D2<T> const &a, D2<T> const &b) {
 template <typename T>
 inline D2<T>
 operator+=(D2<T> &a, D2<T> const &b) {
-    boost::function_requires<AddableConcept<T> >();
+    BOOST_CONCEPT_ASSERT((AddableConcept<T>));
 
     for(unsigned i = 0; i < 2; i++)
         a[i] += b[i];
@@ -220,7 +226,7 @@ operator+=(D2<T> &a, D2<T> const &b) {
 template <typename T>
 inline D2<T>
 operator-=(D2<T> &a, D2<T> const & b) {
-    boost::function_requires<AddableConcept<T> >();
+    BOOST_CONCEPT_ASSERT((AddableConcept<T>));
 
     for(unsigned i = 0; i < 2; i++)
         a[i] -= b[i];
@@ -231,7 +237,7 @@ operator-=(D2<T> &a, D2<T> const & b) {
 template <typename T>
 inline D2<T>
 operator-(D2<T> const & a) {
-    boost::function_requires<ScalableConcept<T> >();
+    BOOST_CONCEPT_ASSERT((ScalableConcept<T>));
     D2<T> r;
     for(unsigned i = 0; i < 2; i++)
         r[i] = -a[i];
@@ -240,7 +246,7 @@ operator-(D2<T> const & a) {
 template <typename T>
 inline D2<T>
 operator*(D2<T> const & a, Point const & b) {
-    boost::function_requires<ScalableConcept<T> >();
+    BOOST_CONCEPT_ASSERT((ScalableConcept<T>));
 
     D2<T> r;
     for(unsigned i = 0; i < 2; i++)
@@ -250,7 +256,7 @@ operator*(D2<T> const & a, Point const & b) {
 template <typename T>
 inline D2<T>
 operator/(D2<T> const & a, Point const & b) {
-    boost::function_requires<ScalableConcept<T> >();
+    BOOST_CONCEPT_ASSERT((ScalableConcept<T>));
     //TODO: b==0?
     D2<T> r;
     for(unsigned i = 0; i < 2; i++)
@@ -260,7 +266,7 @@ operator/(D2<T> const & a, Point const & b) {
 template <typename T>
 inline D2<T>
 operator*=(D2<T> &a, Point const & b) {
-    boost::function_requires<ScalableConcept<T> >();
+    BOOST_CONCEPT_ASSERT((ScalableConcept<T>));
 
     for(unsigned i = 0; i < 2; i++)
         a[i] *= b[i];
@@ -269,7 +275,7 @@ operator*=(D2<T> &a, Point const & b) {
 template <typename T>
 inline D2<T>
 operator/=(D2<T> &a, Point const & b) {
-    boost::function_requires<ScalableConcept<T> >();
+    BOOST_CONCEPT_ASSERT((ScalableConcept<T>));
     //TODO: b==0?
     for(unsigned i = 0; i < 2; i++)
         a[i] /= b[i];
@@ -287,8 +293,8 @@ inline D2<T> operator/=(D2<T> & a, double b) { a[0] /= b; a[1] /= b; return a; }
 
 template<typename T>
 D2<T> operator*(D2<T> const &v, Affine const &m) {
-    boost::function_requires<AddableConcept<T> >();
-    boost::function_requires<ScalableConcept<T> >();
+    BOOST_CONCEPT_ASSERT((AddableConcept<T>));
+    BOOST_CONCEPT_ASSERT((ScalableConcept<T>));
     D2<T> ret;
     for(unsigned i = 0; i < 2; i++)
         ret[i] = v[X] * m[i] + v[Y] * m[i + 2] + m[i + 4];
@@ -299,7 +305,7 @@ D2<T> operator*(D2<T> const &v, Affine const &m) {
 template <typename T>
 inline D2<T>
 operator*(D2<T> const & a, T const & b) {
-    boost::function_requires<MultiplicableConcept<T> >();
+    BOOST_CONCEPT_ASSERT((MultiplicableConcept<T>));
     D2<T> ret;
     for(unsigned i = 0; i < 2; i++)
         ret[i] = a[i] * b;
@@ -312,7 +318,7 @@ operator*(D2<T> const & a, T const & b) {
 template <typename T>
 inline D2<T>
 operator+(D2<T> const & a, Point b) {
-    boost::function_requires<OffsetableConcept<T> >();
+    BOOST_CONCEPT_ASSERT((OffsetableConcept<T>));
     D2<T> r;
     for(unsigned i = 0; i < 2; i++)
         r[i] = a[i] + b[i];
@@ -321,7 +327,7 @@ operator+(D2<T> const & a, Point b) {
 template <typename T>
 inline D2<T>
 operator-(D2<T> const & a, Point b) {
-    boost::function_requires<OffsetableConcept<T> >();
+    BOOST_CONCEPT_ASSERT((OffsetableConcept<T>));
     D2<T> r;
     for(unsigned i = 0; i < 2; i++)
         r[i] = a[i] - b[i];
@@ -330,7 +336,7 @@ operator-(D2<T> const & a, Point b) {
 template <typename T>
 inline D2<T>
 operator+=(D2<T> & a, Point b) {
-    boost::function_requires<OffsetableConcept<T> >();
+    BOOST_CONCEPT_ASSERT((OffsetableConcept<T>));
     for(unsigned i = 0; i < 2; i++)
         a[i] += b[i];
     return a;
@@ -338,7 +344,7 @@ operator+=(D2<T> & a, Point b) {
 template <typename T>
 inline D2<T>
 operator-=(D2<T> & a, Point b) {
-    boost::function_requires<OffsetableConcept<T> >();
+    BOOST_CONCEPT_ASSERT((OffsetableConcept<T>));
     for(unsigned i = 0; i < 2; i++)
         a[i] -= b[i];
     return a;
@@ -347,8 +353,8 @@ operator-=(D2<T> & a, Point b) {
 template <typename T>
 inline T
 dot(D2<T> const & a, D2<T> const & b) {
-    boost::function_requires<AddableConcept<T> >();
-    boost::function_requires<MultiplicableConcept<T> >();
+    BOOST_CONCEPT_ASSERT((AddableConcept<T>));
+    BOOST_CONCEPT_ASSERT((MultiplicableConcept<T>));
 
     T r;
     for(unsigned i = 0; i < 2; i++)
@@ -362,8 +368,8 @@ dot(D2<T> const & a, D2<T> const & b) {
 template <typename T>
 inline T
 dot(D2<T> const & a, Point const & b) {
-    boost::function_requires<AddableConcept<T> >();
-    boost::function_requires<ScalableConcept<T> >();
+    BOOST_CONCEPT_ASSERT((AddableConcept<T>));
+    BOOST_CONCEPT_ASSERT((ScalableConcept<T>));
 
     T r;
     for(unsigned i = 0; i < 2; i++) {
@@ -378,8 +384,8 @@ dot(D2<T> const & a, Point const & b) {
 template <typename T>
 inline T
 cross(D2<T> const & a, D2<T> const & b) {
-    boost::function_requires<ScalableConcept<T> >();
-    boost::function_requires<MultiplicableConcept<T> >();
+    BOOST_CONCEPT_ASSERT((ScalableConcept<T>));
+    BOOST_CONCEPT_ASSERT((MultiplicableConcept<T>));
 
     return a[1] * b[0] - a[0] * b[1];
 }
@@ -389,7 +395,7 @@ cross(D2<T> const & a, D2<T> const & b) {
 template <typename T>
 inline D2<T>
 rot90(D2<T> const & a) {
-    boost::function_requires<ScalableConcept<T> >();
+    BOOST_CONCEPT_ASSERT((ScalableConcept<T>));
     return D2<T>(-a[Y], a[X]);
 }
 
@@ -468,17 +474,17 @@ namespace Geom {
 //Some D2 Fragment implementation which requires rect:
 template <typename T>
 OptRect bounds_fast(const D2<T> &a) {
-    boost::function_requires<FragmentConcept<T> >();
+    BOOST_CONCEPT_ASSERT((FragmentConcept<T>));
     return OptRect(bounds_fast(a[X]), bounds_fast(a[Y]));
 }
 template <typename T>
 OptRect bounds_exact(const D2<T> &a) {
-    boost::function_requires<FragmentConcept<T> >();
+    BOOST_CONCEPT_ASSERT((FragmentConcept<T>));
     return OptRect(bounds_exact(a[X]), bounds_exact(a[Y]));
 }
 template <typename T>
 OptRect bounds_local(const D2<T> &a, const OptInterval &t) {
-    boost::function_requires<FragmentConcept<T> >();
+    BOOST_CONCEPT_ASSERT((FragmentConcept<T>));
     return OptRect(bounds_local(a[X], t), bounds_local(a[Y], t));
 }
 
diff --git a/src/2geom/ellipse.cpp b/src/2geom/ellipse.cpp
index f23c9b24e..46c60d85d 100644
--- a/src/2geom/ellipse.cpp
+++ b/src/2geom/ellipse.cpp
@@ -32,12 +32,18 @@
  */
 
 #include <2geom/ellipse.h>
-#include <2geom/svg-elliptical-arc.h>
+#include <2geom/elliptical-arc.h>
 #include <2geom/numeric/fitting-tool.h>
 #include <2geom/numeric/fitting-model.h>
 
 namespace Geom {
 
+Ellipse::Ellipse(Geom::Circle const &c)
+    : _center(c.center())
+    , _rays(c.radius(), c.radius())
+    , _angle(0)
+{}
+
 void Ellipse::setCoefficients(double A, double B, double C, double D, double E, double F)
 {
     double den = 4*A*C - B*B;
@@ -57,17 +63,11 @@ void Ellipse::setCoefficients(double A, double B, double C, double D, double E,
 
 
     //evaluate ellipse rotation angle
-    double rot = std::atan2( -B, -(A - C) )/2;
-//      std::cerr << "rot = " << rot << std::endl;
-    bool swap_axes = false;
-
-    if (rot >= M_PI/2 || rot < 0) {
-        swap_axes = true;
-    }
+    _angle = std::atan2( -B, -(A - C) )/2;
 
     // evaluate the length of the ellipse rays
     double sinrot, cosrot;
-    sincos(rot, sinrot, cosrot);
+    sincos(_angle, sinrot, cosrot);
     double cos2 = cosrot * cosrot;
     double sin2 = sinrot * sinrot;
     double cossin = cosrot * sinrot;
@@ -80,7 +80,7 @@ void Ellipse::setCoefficients(double A, double B, double C, double D, double E,
     if (rx2 < 0) {
         THROW_RANGEERROR("rx2 < 0, while computing 'rx' coefficient");
     }
-    double rx = std::sqrt(rx2);
+    _rays[X] = std::sqrt(rx2);
 
     den = C * cos2 - B * cossin + A * sin2;
     if (den == 0) {
@@ -90,24 +90,11 @@ void Ellipse::setCoefficients(double A, double B, double C, double D, double E,
     if (ry2 < 0) {
         THROW_RANGEERROR("ry2 < 0, while computing 'rx' coefficient");
     }
-    double ry = std::sqrt(ry2);
+    _rays[Y] = std::sqrt(ry2);
 
     // the solution is not unique so we choose always the ellipse
     // with a rotation angle between 0 and PI/2
-    if (swap_axes) {
-        std::swap(rx, ry);
-    }
-
-    if (rx == ry) {
-        rot = 0;
-    }
-    if (rot < 0) {
-        rot += M_PI/2;
-    }
-
-    _rays[X] = rx;
-    _rays[Y] = ry;
-    _angle = rot;
+    makeCanonical();
 }
 
 
@@ -118,14 +105,49 @@ Affine Ellipse::unitCircleTransform() const
     return ret;
 }
 
+Affine Ellipse::inverseUnitCircleTransform() const
+{
+    if (ray(X) == 0 || ray(Y) == 0) {
+        THROW_RANGEERROR("a degenerate ellipse doesn't have an inverse unit circle transform");
+    }
+    Affine ret = Translate(-center()) * Rotate(-_angle) * Scale(1/ray(X), 1/ray(Y));
+    return ret;
+}
+
+
+LineSegment Ellipse::axis(Dim2 d) const
+{
+    Point a(0, 0), b(0, 0);
+    a[d] = -1;
+    b[d] = 1;
+    LineSegment ls(a, b);
+    ls.transform(unitCircleTransform());
+    return ls;
+}
+
+LineSegment Ellipse::semiaxis(Dim2 d, int sign) const
+{
+    Point a(0, 0), b(0, 0);
+    b[d] = sgn(sign);
+    LineSegment ls(a, b);
+    ls.transform(unitCircleTransform());
+    return ls;
+}
+
 
 std::vector<double> Ellipse::coefficients() const
 {
+    std::vector<double> c(6);
+    coefficients(c[0], c[1], c[2], c[3], c[4], c[5]);
+    return c;
+}
+
+void Ellipse::coefficients(Coord &A, Coord &B, Coord &C, Coord &D, Coord &E, Coord &F) const
+{
     if (ray(X) == 0 || ray(Y) == 0) {
         THROW_RANGEERROR("a degenerate ellipse doesn't have an implicit form");
     }
 
-    std::vector<double> coeff(6);
     double cosrot, sinrot;
     sincos(_angle, sinrot, cosrot);
     double cos2 = cosrot * cosrot;
@@ -134,16 +156,15 @@ std::vector<double> Ellipse::coefficients() const
     double invrx2 = 1 / (ray(X) * ray(X));
     double invry2 = 1 / (ray(Y) * ray(Y));
 
-    coeff[0] = invrx2 * cos2 + invry2 * sin2;
-    coeff[1] = 2 * (invrx2 - invry2) * cossin;
-    coeff[2] = invrx2 * sin2 + invry2 * cos2;
-    coeff[3] = -(2 * coeff[0] * center(X) + coeff[1] * center(Y));
-    coeff[4] = -(2 * coeff[2] * center(Y) + coeff[1] * center(X));
-    coeff[5] = coeff[0] * center(X) * center(X)
-             + coeff[1] * center(X) * center(Y)
-             + coeff[2] * center(Y) * center(Y)
-             - 1;
-    return coeff;
+    A = invrx2 * cos2 + invry2 * sin2;
+    B = 2 * (invrx2 - invry2) * cossin;
+    C = invrx2 * sin2 + invry2 * cos2;
+    D = -2 * A * center(X) - B * center(Y);
+    E = -2 * C * center(Y) - B * center(X);
+    F =   A * center(X) * center(X)
+        + B * center(X) * center(Y)
+        + C * center(Y) * center(Y)
+        - 1;
 }
 
 
@@ -167,8 +188,7 @@ void Ellipse::fit(std::vector<Point> const &points)
 
 
 EllipticalArc *
-Ellipse::arc(Point const &ip, Point const &inner, Point const &fp,
-             bool _svg_compliant)
+Ellipse::arc(Point const &ip, Point const &inner, Point const &fp)
 {
     // This is resistant to degenerate ellipses:
     // both flags evaluate to false in that case.
@@ -196,28 +216,14 @@ Ellipse::arc(Point const &ip, Point const &inner, Point const &fp,
         sweep_flag = true;
     }
 
-    EllipticalArc *ret_arc;
-    if (_svg_compliant) {
-        ret_arc = new SVGEllipticalArc(ip, ray(X), ray(Y), rotationAngle(),
-                                       large_arc_flag, sweep_flag, fp);
-    } else {
-        ret_arc = new EllipticalArc(ip, ray(X), ray(Y), rotationAngle(),
-                                    large_arc_flag, sweep_flag, fp);
-    }
+    EllipticalArc *ret_arc = new EllipticalArc(ip, ray(X), ray(Y), rotationAngle(),
+                                               large_arc_flag, sweep_flag, fp);
     return ret_arc;
 }
 
 Ellipse &Ellipse::operator*=(Rotate const &r)
 {
     _angle += r.angle();
-    // keep the angle in the first quadrant
-    if (_angle < 0) {
-        _angle += M_PI;
-    }
-    if (_angle >= M_PI/2) {
-        std::swap(_rays[X], _rays[Y]);
-        _angle -= M_PI/2;
-    }
     _center *= r;
     return *this;
 }
@@ -261,10 +267,311 @@ Ellipse &Ellipse::operator*=(Affine const& m)
     return *this;
 }
 
-Ellipse::Ellipse(Geom::Circle const &c)
+Ellipse Ellipse::canonicalForm() const
+{
+    Ellipse result(*this);
+    result.makeCanonical();
+    return result;
+}
+
+void Ellipse::makeCanonical()
+{
+    if (_rays[X] == _rays[Y]) {
+        _angle = 0;
+        return;
+    }
+
+    if (_angle < 0) {
+        _angle += M_PI;
+    }
+    if (_angle >= M_PI/2) {
+        std::swap(_rays[X], _rays[Y]);
+        _angle -= M_PI/2;
+    }
+}
+
+Point Ellipse::pointAt(Coord t) const
+{
+    Point p = Point::polar(t);
+    p *= unitCircleTransform();
+    return p;
+}
+
+Coord Ellipse::valueAt(Coord t, Dim2 d) const
+{
+    // TODO: more efficient version.
+    return pointAt(t)[d];
+}
+
+Coord Ellipse::timeAt(Point const &p) const
+{
+    // degenerate ellipse is basically a reparametrized line segment
+    if (ray(X) == 0 || ray(Y) == 0) {
+        if (ray(X) != 0) {
+            return asin(Line(axis(X)).timeAt(p));
+        } else if (ray(Y) != 0) {
+            return acos(Line(axis(Y)).timeAt(p));
+        } else {
+            return 0;
+        }
+    }
+    Affine iuct = inverseUnitCircleTransform();
+    return Angle(atan2(p * iuct)).radians0(); // return a value in [0, 2pi)
+}
+
+std::vector<ShapeIntersection> Ellipse::intersect(Line const &line) const
+{
+
+    std::vector<ShapeIntersection> result;
+
+    if (line.isDegenerate()) return result;
+    if (ray(X) == 0 || ray(Y) == 0) {
+        // TODO intersect with line segment.
+        return result;
+    }
+
+    // Ax^2 + Bxy + Cy^2 + Dx + Ey + F
+    Coord A, B, C, D, E, F;
+    coefficients(A, B, C, D, E, F);
+    Affine iuct = inverseUnitCircleTransform();
+
+    if (line.isHorizontal()) {
+        // substitute y into the ellipse equation and solve
+        Coord y = line.initialPoint()[Y];
+        std::vector<Coord> xs = solve_quadratic(A, B*y + D, C*y*y + E*y + F);
+        for (unsigned i = 0; i < xs.size(); ++i) {
+            Point p(xs[i], y);
+            result.push_back(ShapeIntersection(atan2(p * iuct), line.timeAt(p), p));
+        }
+        return result;
+    }
+    if (line.isVertical()) {
+        // substitute y into the ellipse equation and solve
+        Coord x = line.initialPoint()[X];
+        std::vector<Coord> ys = solve_quadratic(C, B*x + E, A*x*x + D*x + F);
+        for (unsigned i = 0; i < ys.size(); ++i) {
+            Point p(x, ys[i]);
+            result.push_back(ShapeIntersection(atan2(p * iuct), line.timeAt(p), p));
+        }
+        return result;
+    }
+
+    // generic case
+    Coord a, b, c;
+    line.coefficients(a, b, c);
+
+    // y = -a/b x - C/B
+    // TODO: when is it better to substitute X?
+    Coord q = -a/b;
+    Coord r = -c/b;
+
+    // substitute that into the ellipse equation, making it quadratic in x
+    Coord I = A + B*q + C*q*q;          // x^2 terms
+    Coord J = B*r + C*2*q*r + D + E*q;  // x^1 terms
+    Coord K = C*r*r + E*r + F;          // x^0 terms
+    std::vector<Coord> xs = solve_quadratic(I, J, K);
+
+    for (unsigned i = 0; i < xs.size(); ++i) {
+        Point p(xs[i], q*xs[i] + r);
+        result.push_back(ShapeIntersection(atan2(p * iuct), line.timeAt(p), p));
+    }
+    return result;
+}
+
+std::vector<ShapeIntersection> Ellipse::intersect(LineSegment const &seg) const
+{
+    // we simply re-use the procedure for lines and filter out
+    // results where the line time value is outside of the unit interval.
+    std::vector<ShapeIntersection> result = intersect(Line(seg));
+    filter_line_segment_intersections(result);
+    return result;
+}
+
+std::vector<ShapeIntersection> Ellipse::intersect(Ellipse const &other) const
+{
+    // handle degenerate cases first
+    if (ray(X) == 0 || ray(Y) == 0) {
+        
+    }
+    // intersection of two ellipses can be solved analytically.
+    // http://maptools.home.comcast.net/~maptools/BivariateQuadratics.pdf
+
+    Coord A, B, C, D, E, F;
+    Coord a, b, c, d, e, f;
+
+    // NOTE: the order of coefficients is different to match the convention in the PDF above
+    // Ax^2 + Bx^2 + Cx + Dy + Exy + F
+    this->coefficients(A, E, B, C, D, F);
+    other.coefficients(a, e, b, c, d, f);
+
+    // Assume that Q is the ellipse equation given by uppercase letters
+    // and R is the equation given by lowercase ones. An intersection exists when
+    // there is a coefficient mu such that
+    // mu Q + R = 0
+    //
+    // This can be written in the following way:
+    //
+    //                    |  ff  cc/2 dd/2 | |1|
+    // mu Q + R = [1 x y] | cc/2  aa  ee/2 | |x| = 0
+    //                    | dd/2 ee/2  bb  | |y|
+    //
+    // where aa = mu A + a and so on. The determinant can be explicitly written out,
+    // giving an equation which is cubic in mu and can be solved analytically.
+
+    Coord I, J, K, L;
+    I = (-E*E*F + 4*A*B*F + C*D*E - A*D*D - B*C*C) / 4;
+    J = -((E*E - 4*A*B) * f + (2*E*F - C*D) * e + (2*A*D - C*E) * d +
+          (2*B*C - D*E) * c + (C*C - 4*A*F) * b + (D*D - 4*B*F) * a) / 4;
+    K = -((e*e - 4*a*b) * F + (2*e*f - c*d) * E + (2*a*d - c*e) * D +
+          (2*b*c - d*e) * C + (c*c - 4*a*f) * B + (d*d - 4*b*f) * A) / 4;
+    L = (-e*e*f + 4*a*b*f + c*d*e - a*d*d - b*c*c) / 4;
+
+    std::vector<Coord> mus = solve_cubic(I, J, K, L);
+    Coord mu = infinity();
+    std::vector<ShapeIntersection> result;
+
+    // Now that we have solved for mu, we need to check whether the conic
+    // determined by mu Q + R is reducible to a product of two lines. If it's not,
+    // it means that there are no intersections. If it is, the intersections of these
+    // lines with the original ellipses (if there are any) give the coordinates
+    // of intersections.
+
+    // Prefer middle root if there are three.
+    // Out of three possible pairs of lines that go through four points of intersection
+    // of two ellipses, this corresponds to cross-lines. These intersect the ellipses
+    // at less shallow angles than the other two options.
+    if (mus.size() == 3) {
+        std::swap(mus[1], mus[0]);
+    }
+    for (unsigned i = 0; i < mus.size(); ++i) {
+        Coord aa = mus[i] * A + a;
+        Coord bb = mus[i] * B + b;
+        Coord ee = mus[i] * E + e;
+        Coord delta = ee*ee - 4*aa*bb;
+        if (delta < 0) continue;
+        mu = mus[i];
+        break;
+    }
+
+    // if no suitable mu was found, there are no intersections
+    if (mu == infinity()) return result;
+
+    Coord aa = mu * A + a;
+    Coord bb = mu * B + b;
+    Coord cc = mu * C + c;
+    Coord dd = mu * D + d;
+    Coord ee = mu * E + e;
+    Coord ff = mu * F + f;
+
+    Line lines[2];
+
+    if (aa != 0) {
+        bb /= aa; cc /= aa; dd /= aa; ee /= aa; ff /= aa;
+        Coord s = (ee + std::sqrt(ee*ee - 4*bb)) / 2;
+        Coord q = ee - s;
+        Coord alpha = (dd - cc*q) / (s - q);
+        Coord beta = cc - alpha;
+
+        lines[0] = Line(1, q, alpha);
+        lines[1] = Line(1, s, beta);
+    } else if (bb != 0) {
+        cc /= bb; dd /= bb; ee /= bb; ff /= bb;
+        Coord s = ee;
+        Coord q = 0;
+        Coord alpha = cc / ee;
+        Coord beta = ff * ee / cc;
+
+        lines[0] = Line(q, 1, alpha);
+        lines[1] = Line(s, 1, beta);
+    } else {
+        lines[0] = Line(ee, 0, dd);
+        lines[1] = Line(0, 1, cc/ee);
+    }
+
+    // intersect with the obtained lines and report intersections
+    for (unsigned li = 0; li < 2; ++li) {
+        std::vector<ShapeIntersection> as = intersect(lines[li]);
+        std::vector<ShapeIntersection> bs = other.intersect(lines[li]);
+
+        if (!as.empty() && as.size() == bs.size()) {
+            for (unsigned i = 0; i < as.size(); ++i) {
+                ShapeIntersection ix(as[i].first, bs[i].first,
+                    middle_point(as[i].point(), bs[i].point()));
+                result.push_back(ix);
+            }
+        }
+    }
+    return result;
+}
+
+std::vector<ShapeIntersection> Ellipse::intersect(D2<Bezier> const &b) const
+{
+    Coord A, B, C, D, E, F;
+    coefficients(A, B, C, D, E, F);
+
+    Bezier x = A*b[X]*b[X] + B*b[X]*b[Y] + C*b[Y]*b[Y] + D*b[X] + E*b[Y] + F;
+    std::vector<Coord> r = x.roots();
+
+    std::vector<ShapeIntersection> result;
+    for (unsigned i = 0; i < r.size(); ++i) {
+        Point p = b.valueAt(r[i]);
+        result.push_back(ShapeIntersection(timeAt(p), r[i], p));
+    }
+    return result;
+}
+
+bool Ellipse::operator==(Ellipse const &other) const
+{
+    if (_center != other._center) return false;
+
+    Ellipse a = this->canonicalForm();
+    Ellipse b = other.canonicalForm();
+
+    if (a._rays != b._rays) return false;
+    if (a._angle != b._angle) return false;
+
+    return true;
+}
+
+
+bool are_near(Ellipse const &a, Ellipse const &b, Coord precision)
+{
+    // We want to know whether no point on ellipse a is further than precision
+    // from the corresponding point on ellipse b. To check this, we compute
+    // the four extreme points at the end of each ray for each ellipse
+    // and check whether they are sufficiently close.
+
+    // First, we need to correct the angles on the ellipses, so that they are
+    // no further than M_PI/4 apart. This can always be done by rotating
+    // and exchanging axes.
+    Ellipse ac = a, bc = b;
+    if (distance(ac.rotationAngle(), bc.rotationAngle()).radians0() >= M_PI/2) {
+        ac.setRotationAngle(ac.rotationAngle() + M_PI);
+    }
+    if (distance(ac.rotationAngle(), bc.rotationAngle()) >= M_PI/4) {
+        Angle d1 = distance(ac.rotationAngle() + M_PI/2, bc.rotationAngle());
+        Angle d2 = distance(ac.rotationAngle() - M_PI/2, bc.rotationAngle());
+        Coord adj = d1.radians0() < d2.radians0() ? M_PI/2 : -M_PI/2;
+        ac.setRotationAngle(ac.rotationAngle() + adj);
+        ac.setRays(ac.ray(Y), ac.ray(X));
+    }
+
+    // Do the actual comparison by computing four points on each ellipse.
+    Point tps[] = {Point(1,0), Point(0,1), Point(-1,0), Point(0,-1)};
+    for (unsigned i = 0; i < 4; ++i) {
+        if (!are_near(tps[i] * ac.unitCircleTransform(),
+                      tps[i] * bc.unitCircleTransform(),
+                      precision))
+            return false;
+    }
+    return true;
+}
+
+std::ostream &operator<<(std::ostream &out, Ellipse const &e)
 {
-    _center = c.center();
-    _rays[X] = _rays[Y] = c.radius();
+    out << "Ellipse(" << e.center() << ", " << e.rays()
+        << ", " << format_coord_nice(e.rotationAngle()) << ")";
+    return out;
 }
 
 }  // end namespace Geom
diff --git a/src/2geom/ellipse.h b/src/2geom/ellipse.h
index a67969933..6fb5ed254 100644
--- a/src/2geom/ellipse.h
+++ b/src/2geom/ellipse.h
@@ -37,8 +37,10 @@
 
 #include <vector>
 #include <2geom/angle.h>
+#include <2geom/bezier-curve.h>
 #include <2geom/exception.h>
-#include <2geom/point.h>
+#include <2geom/forward.h>
+#include <2geom/line.h>
 #include <2geom/transforms.h>
 
 namespace Geom {
@@ -46,7 +48,14 @@ namespace Geom {
 class EllipticalArc;
 class Circle;
 
-/** @brief Set of points with a constant sum of distances from two foci
+/** @brief Set of points with a constant sum of distances from two foci.
+ *
+ * An ellipse can be specified in several ways. Internally, 2Geom uses
+ * the SVG style representation: center, rays and angle between the +X ray
+ * and the +X axis. Another popular way is to use an implicit equation,
+ * which is as follows:
+ * \f$Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0\f$
+ *
  * @ingroup Shapes */
 class Ellipse
     : boost::multipliable< Ellipse, Translate
@@ -54,7 +63,8 @@ class Ellipse
     , boost::multipliable< Ellipse, Rotate
     , boost::multipliable< Ellipse, Zoom
     , boost::multipliable< Ellipse, Affine
-      > > > > >
+    , boost::equality_comparable< Ellipse
+      > > > > > >
 {
     Point _center;
     Point _rays;
@@ -74,13 +84,16 @@ public:
     Ellipse(double A, double B, double C, double D, double E, double F) {
         setCoefficients(A, B, C, D, E, F);
     }
+    /// Construct ellipse from a circle.
     Ellipse(Geom::Circle const &c);
 
+    /// Set center, rays and angle.
     void set(Point const &c, Point const &r, Coord angle) {
         _center = c;
         _rays = r;
         _angle = angle;
     }
+    /// Set center, rays and angle as constituent values.
     void set(Coord cx, Coord cy, Coord rx, Coord ry, Coord a) {
         _center[X] = cx;
         _center[Y] = cy;
@@ -88,31 +101,97 @@ public:
         _rays[Y] = ry;
         _angle = a;
     }
-
-    // build an ellipse by its implicit equation:
-    // Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0
+    /// Set an ellipse by solving its implicit equation.
     void setCoefficients(double A, double B, double C, double D, double E, double F);
-
-    // biuld up the best fitting ellipse wrt the passed points
-    // prerequisite: at least 5 points must be passed
-    void fit(std::vector<Point> const& points);
-
-    EllipticalArc *arc(Point const &ip, Point const &inner, Point const &fp,
-                       bool svg_compliant = true);
+    /// Set the center.
+    void setCenter(Point const &p) { _center = p; }
+    /// Set the center by coordinates.
+    void setCenter(Coord cx, Coord cy) { _center[X] = cx; _center[Y] = cy; }
+    /// Set both rays of the ellipse.
+    void setRays(Point const &p) { _rays = p; }
+    /// Set both rays of the ellipse as coordinates.
+    void setRays(Coord x, Coord y) { _rays[X] = x; _rays[Y] = y; }
+    /// Set one of the rays of the ellipse.
+    void setRay(Coord r, Dim2 d) { _rays[d] = r; }
+    /// Set the angle the X ray makes with the +X axis.
+    void setRotationAngle(Angle a) { _angle = a; }
 
     Point center() const { return _center; }
     Coord center(Dim2 d) const { return _center[d]; }
+    /// Get both rays as a point.
     Point rays() const { return _rays; }
+    /// Get one ray of the ellipse.
     Coord ray(Dim2 d) const { return _rays[d]; }
+    /// Get the angle the X ray makes with the +X axis.
     Angle rotationAngle() const { return _angle; }
 
+    /** @brief Create an ellipse passing through the specified points
+     * At least five points have to be specified. */
+    void fit(std::vector<Point> const& points);
+
+    /** @brief Create an elliptical arc from a section of the ellipse.
+     * This is mainly useful to determine the flags of the new arc.
+     * The passed points should lie on the ellipse, otherwise the results
+     * will be undefined.
+     * @param ip Initial point of the arc
+     * @param inner Point in the middle of the arc, used to pick one of two possibilities
+     * @param fp Final point of the arc
+     * @return Newly allocated arc, delete when no longer used */
+    EllipticalArc *arc(Point const &ip, Point const &inner, Point const &fp);
+
+    /** @brief Return an ellipse with less degrees of freedom.
+     * The canonical form always has the angle less than \f$\frac{\pi}{2}\f$,
+     * and zero if the rays are equal (i.e. the ellipse is a circle). */
+    Ellipse canonicalForm() const;
+    void makeCanonical();
+
     /** @brief Compute the transform that maps the unit circle to this 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 this ellipse.
      * @return Transform from unit circle to the ellipse */
     Affine unitCircleTransform() const;
+    /** @brief Compute the transform that maps this ellipse to the unit circle.
+     * This may be a little more precise and/or faster than simply using
+     * unitCircleTransform().inverse(). An exception will be thrown for
+     * degenerate ellipses. */
+    Affine inverseUnitCircleTransform() const;
+
+    LineSegment majorAxis() const { return ray(X) >= ray(Y) ? axis(X) : axis(Y); }
+    LineSegment minorAxis() const { return ray(X) < ray(Y) ? axis(X) : axis(Y); }
+    LineSegment semimajorAxis(int sign = 1) const {
+        return ray(X) >= ray(Y) ? semiaxis(X, sign) : semiaxis(Y, sign);
+    }
+    LineSegment semiminorAxis(int sign = 1) const {
+        return ray(X) < ray(Y) ? semiaxis(X, sign) : semiaxis(Y, sign);
+    }
+    LineSegment axis(Dim2 d) const;
+    LineSegment semiaxis(Dim2 d, int sign = 1) const;
 
+    /// Get the coefficients of the ellipse's implicit equation.
     std::vector<double> coefficients() const;
+    void coefficients(Coord &A, Coord &B, Coord &C, Coord &D, Coord &E, Coord &F) const;
+
+    /** @brief Evaluate a point on the ellipse.
+     * The parameter range is \f$[0, 2\pi)\f$; larger and smaller values
+     * wrap around. */
+    Point pointAt(Coord t) const;
+    /// Evaluate a single coordinate of a point on the ellipse.
+    Coord valueAt(Coord t, Dim2 d) const;
+
+    /** @brief Find the time value of a point on an ellipse.
+     * If the point is not on the ellipse, the returned time value will correspond
+     * to an intersection with a ray from the origin passing through the point
+     * with the ellipse. Note that this is NOT the nearest point on the ellipse. */
+    Coord timeAt(Point const &p) const;
+
+    /// Compute intersections with an infinite line.
+    std::vector<ShapeIntersection> intersect(Line const &line) const;
+    /// Compute intersections with a line segment.
+    std::vector<ShapeIntersection> intersect(LineSegment const &seg) const;
+    /// Compute intersections with another ellipse.
+    std::vector<ShapeIntersection> intersect(Ellipse const &other) const;
+    /// Compute intersections with a 2D Bezier polynomial.
+    std::vector<ShapeIntersection> intersect(D2<Bezier> const &other) const;
 
     Ellipse &operator*=(Translate const &t) {
         _center *= t;
@@ -130,8 +209,21 @@ public:
     }
     Ellipse &operator*=(Rotate const &r);
     Ellipse &operator*=(Affine const &m);
+
+    /// Compare ellipses for exact equality.
+    bool operator==(Ellipse const &other) const;
 };
 
+/** @brief Test whether two ellipses are approximately the same.
+ * This will check whether no point on ellipse a is further away from
+ * the corresponding point on ellipse b than precision.
+ * @relates Ellipse */
+bool are_near(Ellipse const &a, Ellipse const &b, Coord precision = EPSILON);
+
+/** @brief Outputs ellipse data, useful for debugging.
+ * @relates Ellipse */
+std::ostream &operator<<(std::ostream &out, Ellipse const &e);
+
 } // end namespace Geom
 
 #endif // LIB2GEOM_SEEN_ELLIPSE_H
diff --git a/src/2geom/svg-elliptical-arc.cpp b/src/2geom/elliptical-arc-from-sbasis.cpp
index 3a5154d08..54f995fb6 100644
--- a/src/2geom/svg-elliptical-arc.cpp
+++ b/src/2geom/elliptical-arc-from-sbasis.cpp
@@ -1,6 +1,8 @@
-/*
- * SVG Elliptical Arc Class
+/** @file
+ * @brief Fitting elliptical arc to SBasis
  *
+ * This file contains the implementation of the function arc_from_sbasis.
+ *//*
  * Copyright 2008  Marco Cecchetti <mrcekets at gmail.com>
  *
  * This library is free software; you can redistribute it and/or
@@ -27,33 +29,120 @@
  * the specific language governing rights and limitations.
  */
 
-
+#include <2geom/curve.h>
+#include <2geom/angle.h>
+#include <2geom/utils.h>
 #include <2geom/bezier-curve.h>
-#include <2geom/ellipse.h>
-#include <2geom/numeric/fitting-model.h>
-#include <2geom/numeric/fitting-tool.h>
+#include <2geom/elliptical-arc.h>
+#include <2geom/sbasis-curve.h>  // for non-native methods
 #include <2geom/numeric/vector.h>
-#include <2geom/poly.h>
-#include <2geom/sbasis-geometric.h>
-#include <2geom/svg-elliptical-arc.h>
-
-#include <cfloat>
-#include <limits>
-#include <memory>
+#include <2geom/numeric/fitting-tool.h>
+#include <2geom/numeric/fitting-model.h>
+#include <algorithm>
 
+namespace Geom {
 
-namespace Geom
+// forward declation
+namespace detail
 {
+    struct ellipse_equation;
+}
 
-/**
- * @class SVGEllipticalArc
- * @brief SVG 1.1-compliant elliptical arc.
+/*
+ * make_elliptical_arc
  *
- * This class is almost identical to the normal elliptical arc, but it differs slightly
- * in the handling of degenerate arcs to be compliant with SVG 1.1 implementation guidelines.
+ * convert a parametric polynomial curve given in symmetric power basis form
+ * into an EllipticalArc type; in order to be successfull the input curve
+ * has to look like an actual elliptical arc even if a certain tolerance
+ * is allowed through an ad-hoc parameter.
+ * The conversion is performed through an interpolation on a certain amount of
+ * sample points computed on the input curve;
+ * the interpolation computes the coefficients of the general implicit equation
+ * of an ellipse (A*X^2 + B*XY + C*Y^2 + D*X + E*Y + F = 0), then from the
+ * implicit equation we compute the parametric form.
  *
- * @ingroup Curves
  */
+class make_elliptical_arc
+{
+  public:
+    typedef D2<SBasis> curve_type;
+
+    /*
+     * constructor
+     *
+     * it doesn't execute the conversion but set the input and output parameters
+     *
+     * _ea:         the output EllipticalArc that will be generated;
+     * _curve:      the input curve to be converted;
+     * _total_samples: the amount of sample points to be taken
+     *                 on the input curve for performing the conversion
+     * _tolerance:     how much likelihood is required between the input curve
+     *                 and the generated elliptical arc; the smaller it is the
+     *                 the tolerance the higher it is the likelihood.
+     */
+    make_elliptical_arc( EllipticalArc& _ea,
+                         curve_type const& _curve,
+                         unsigned int _total_samples,
+                         double _tolerance );
+
+  private:
+    bool bound_exceeded( unsigned int k, detail::ellipse_equation const & ee,
+                         double e1x, double e1y, double e2 );
+
+    bool check_bound(double A, double B, double C, double D, double E, double F);
+
+    void fit();
+
+    bool make_elliptiarc();
+
+    void print_bound_error(unsigned int k)
+    {
+        std::cerr
+            << "tolerance error" << std::endl
+            << "at point: " << k << std::endl
+            << "error value: "<< dist_err << std::endl
+            << "bound: " << dist_bound << std::endl
+            << "angle error: " << angle_err
+            << " (" << angle_tol << ")" << std::endl;
+    }
+
+  public:
+    /*
+     * perform the actual conversion
+     * return true if the conversion is successfull, false on the contrary
+     */
+    bool operator()()
+    {
+        // initialize the reference
+        const NL::Vector & coeff = fitter.result();
+        fit();
+        if ( !check_bound(1, coeff[0], coeff[1], coeff[2], coeff[3], coeff[4]) )
+            return false;
+        if ( !(make_elliptiarc()) ) return false;
+        return true;
+    }
+
+  private:
+      EllipticalArc& ea;                 // output elliptical arc
+      const curve_type & curve;             // input curve
+      Piecewise<D2<SBasis> > dcurve;        // derivative of the input curve
+      NL::LFMEllipse model;                 // model used for fitting
+      // perform the actual fitting task
+      NL::least_squeares_fitter<NL::LFMEllipse> fitter;
+      // tolerance: the user-defined tolerance parameter;
+      // tol_at_extr: the tolerance at end-points automatically computed
+      // on the value of "tolerance", and usually more strict;
+      // tol_at_center: tolerance at the center of the ellipse
+      // angle_tol: tolerance for the angle btw the input curve tangent
+      // versor and the ellipse normal versor at the sample points
+      double tolerance, tol_at_extr, tol_at_center, angle_tol;
+      Point initial_point, final_point;     // initial and final end-points
+      unsigned int N;                       // total samples
+      unsigned int last; // N-1
+      double partitions; // N-1
+      std::vector<Point> p;                 // sample points
+      double dist_err, dist_bound, angle_err;
+};
 
 namespace detail
 {
@@ -97,7 +186,7 @@ struct ellipse_equation
     double A, B, C, D, E, F;
 };
 
-}
+} // end namespace detail
 
 make_elliptical_arc::
 make_elliptical_arc( EllipticalArc& _ea,
@@ -110,8 +199,7 @@ make_elliptical_arc( EllipticalArc& _ea,
       tolerance(_tolerance), tol_at_extr(tolerance/2),
       tol_at_center(0.1), angle_tol(0.1),
       initial_point(curve.at0()), final_point(curve.at1()),
-      N(_total_samples), last(N-1), partitions(N-1), p(N),
-      svg_compliant(true)
+      N(_total_samples), last(N-1), partitions(N-1), p(N)
 {
 }
 
@@ -216,33 +304,12 @@ bool make_elliptical_arc::make_elliptiarc()
 
     Point inner_point = curve(0.5);
 
-    if (svg_compliant_flag())
-    {
-#ifdef CPP11
-        std::unique_ptr<EllipticalArc> arc( e.arc(initial_point, inner_point, final_point, true) );
-#else
-        std::auto_ptr<EllipticalArc> arc( e.arc(initial_point, inner_point, final_point, true) );
-#endif
-        ea = *arc;
-    }
-    else
-    {
-        try
-        {
 #ifdef CPP11
-            std::unique_ptr<EllipticalArc>
+    std::unique_ptr<EllipticalArc> arc( e.arc(initial_point, inner_point, final_point) );
 #else
-            std::auto_ptr<EllipticalArc>
+    std::auto_ptr<EllipticalArc> arc( e.arc(initial_point, inner_point, final_point) );
 #endif
-                eap( e.arc(initial_point, inner_point, final_point, false) );
-            ea = *eap;
-        }
-        catch(RangeError const &exc)
-        {
-            return false;
-        }
-    }
-
+    ea = *arc;
 
     if ( !are_near( e.center(),
                     ea.center(),
@@ -257,10 +324,14 @@ bool make_elliptical_arc::make_elliptiarc()
 
 
 
-} // end namespace Geom
-
-
+bool arc_from_sbasis(EllipticalArc &ea, D2<SBasis> const &in,
+                     double tolerance, unsigned num_samples)
+{
+    make_elliptical_arc convert(ea, in, num_samples, tolerance);
+    return convert();
+}
 
+} // end namespace Geom
 
 /*
   Local Variables:
@@ -272,4 +343,3 @@ bool make_elliptical_arc::make_elliptiarc()
   End:
 */
 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
-
diff --git a/src/2geom/elliptical-arc.cpp b/src/2geom/elliptical-arc.cpp
index 601f0c8c8..a04b730b3 100644
--- a/src/2geom/elliptical-arc.cpp
+++ b/src/2geom/elliptical-arc.cpp
@@ -38,7 +38,6 @@
 #include <2geom/ellipse.h>
 #include <2geom/elliptical-arc.h>
 #include <2geom/path-sink.h>
-#include <2geom/poly.h>
 #include <2geom/sbasis-geometric.h>
 #include <2geom/transforms.h>
 #include <2geom/utils.h>
@@ -93,11 +92,16 @@ namespace Geom
 
 Rect EllipticalArc::boundsExact() const
 {
+    if (isChord()) {
+        return chord().boundsExact();
+    }
+
     using std::swap;
 
+    // TODO: simplify / document what is going on here.
     double extremes[4];
     double sinrot, cosrot;
-    sincos(_rot_angle, sinrot, cosrot);
+    sincos(rotationAngle(), sinrot, cosrot);
 
     extremes[0] = std::atan2( -ray(Y) * sinrot, ray(X) * cosrot );
     extremes[1] = extremes[0] + M_PI;
@@ -132,14 +136,14 @@ Rect EllipticalArc::boundsExact() const
 
 Point EllipticalArc::pointAtAngle(Coord t) const
 {
-    Point ret = Point::polar(t) * unitCircleTransform();
+    Point ret = _ellipse.pointAt(t);
     return ret;
 }
 
 Coord EllipticalArc::valueAtAngle(Coord t, Dim2 d) const
 {
     Coord sinrot, cosrot, cost, sint;
-    sincos(_rot_angle, sinrot, cosrot);
+    sincos(rotationAngle(), sinrot, cosrot);
     sincos(t, sint, cost);
 
     if ( d == X ) {
@@ -155,13 +159,22 @@ Coord EllipticalArc::valueAtAngle(Coord t, Dim2 d) const
 
 Affine EllipticalArc::unitCircleTransform() const
 {
-    Affine ret = Scale(ray(X), ray(Y)) * Rotate(_rot_angle);
-    ret.setTranslation(center());
+    Affine ret = _ellipse.unitCircleTransform();
+    return ret;
+}
+
+Affine EllipticalArc::inverseUnitCircleTransform() const
+{
+    Affine ret = _ellipse.inverseUnitCircleTransform();
     return ret;
 }
 
 std::vector<Coord> EllipticalArc::roots(Coord v, Dim2 d) const
 {
+    if (isChord()) {
+        return chord().roots(v, d);
+    }
+
     std::vector<Coord> sol;
 
     if ( are_near(ray(X), 0) && are_near(ray(Y), 0) ) {
@@ -190,7 +203,7 @@ std::vector<Coord> EllipticalArc::roots(Coord v, Dim2 d) const
 
     for ( unsigned int dim = 0; dim < 2; ++dim )
     {
-        if ( are_near(ray((Dim2) dim), 0) )
+        if (ray((Dim2) dim) == 0)
         {
             if ( initialPoint()[d] == v && finalPoint()[d] == v )
             {
@@ -212,18 +225,18 @@ std::vector<Coord> EllipticalArc::roots(Coord v, Dim2 d) const
                 case X:
                     switch(dim)
                     {
-                        case X: ray_prj = -ray(Y) * std::sin(_rot_angle);
+                        case X: ray_prj = -ray(Y) * std::sin(rotationAngle());
                                 break;
-                        case Y: ray_prj = ray(X) * std::cos(_rot_angle);
+                        case Y: ray_prj = ray(X) * std::cos(rotationAngle());
                                 break;
                     }
                     break;
                 case Y:
                     switch(dim)
                     {
-                        case X: ray_prj = ray(Y) * std::cos(_rot_angle);
+                        case X: ray_prj = ray(Y) * std::cos(rotationAngle());
                                 break;
-                        case Y: ray_prj = ray(X) * std::sin(_rot_angle);
+                        case Y: ray_prj = ray(X) * std::sin(rotationAngle());
                                 break;
                     }
                     break;
@@ -269,10 +282,10 @@ std::vector<Coord> EllipticalArc::roots(Coord v, Dim2 d) const
 
     double rotx, roty;
     if (d == X) {
-        sincos(_rot_angle, roty, rotx);
+        sincos(rotationAngle(), roty, rotx);
         roty = -roty;
     } else {
-        sincos(_rot_angle, rotx, roty);
+        sincos(rotationAngle(), rotx, roty);
     }
 
     double rxrotx = ray(X) * rotx;
@@ -336,8 +349,12 @@ std::vector<Coord> EllipticalArc::roots(Coord v, Dim2 d) const
 // of such an angle in the cw direction
 Curve *EllipticalArc::derivative() const
 {
+    if (isChord()) {
+        return chord().derivative();
+    }
+
     EllipticalArc *result = static_cast<EllipticalArc*>(duplicate());
-    result->_center[X] = result->_center[Y] = 0;
+    result->_ellipse.setCenter(0, 0);
     result->_start_angle += M_PI/2;
     if( !( result->_start_angle < 2*M_PI ) )
     {
@@ -357,13 +374,17 @@ Curve *EllipticalArc::derivative() const
 std::vector<Point>
 EllipticalArc::pointAndDerivatives(Coord t, unsigned int n) const
 {
+    if (isChord()) {
+        return chord().pointAndDerivatives(t, n);
+    }
+
     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);
     std::auto_ptr<EllipticalArc> ea( static_cast<EllipticalArc*>(duplicate()) );
-    ea->_center = Point(0,0);
+    ea->_ellipse.setCenter(0, 0);
     unsigned int m = std::min(nn, 4u);
     for ( unsigned int i = 0; i < m; ++i )
     {
@@ -392,6 +413,12 @@ 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));
+}
+
 Curve* EllipticalArc::portion(double f, double t) const
 {
     // fix input arguments
@@ -400,14 +427,14 @@ Curve* EllipticalArc::portion(double f, double t) const
     if (t < 0) t = 0;
     if (t > 1) t = 1;
 
-    if ( are_near(f, t) )
+    if (f == t)
     {
         EllipticalArc *arc = static_cast<EllipticalArc*>(duplicate());
-        arc->_center = arc->_initial_point = arc->_final_point = pointAt(f);
-        arc->_start_angle = arc->_end_angle = _start_angle;
-        arc->_rot_angle = _rot_angle;
-        arc->_sweep = _sweep;
-        arc->_large_arc = _large_arc;
+        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;
     }
 
@@ -415,21 +442,24 @@ Curve* EllipticalArc::portion(double f, double t) const
     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)
+    if ( _large_arc && fabs(sweepAngle() * (t-f)) < M_PI) {
         arc->_large_arc = false;
-    arc->_updateCenterAndAngles(arc->isSVGCompliant()); //TODO: be more clever
+    }
+    //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 {
+Curve *EllipticalArc::reverse() const
+{
     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->_updateCenterAndAngles(rarc->isSVGCompliant());
     return rarc;
 }
 
@@ -455,8 +485,8 @@ std::vector<double> EllipticalArc::allNearestTimes( Point const& p, double from,
         Point np = seg.pointAt( seg.nearestTime(p) );
         if ( are_near(ray(Y), 0) )
         {
-            if ( are_near(_rot_angle, M_PI/2)
-                 || are_near(_rot_angle, 3*M_PI/2) )
+            if ( are_near(rotationAngle(), M_PI/2)
+                 || are_near(rotationAngle(), 3*M_PI/2) )
             {
                 result = roots(np[Y], Y);
             }
@@ -467,8 +497,8 @@ std::vector<double> EllipticalArc::allNearestTimes( Point const& p, double from,
         }
         else
         {
-            if ( are_near(_rot_angle, M_PI/2)
-                 || are_near(_rot_angle, 3*M_PI/2) )
+            if ( are_near(rotationAngle(), M_PI/2)
+                 || are_near(rotationAngle(), 3*M_PI/2) )
             {
                 result = roots(np[X], X);
             }
@@ -527,7 +557,7 @@ std::vector<double> EllipticalArc::allNearestTimes( Point const& p, double from,
     Point p_c = p - center();
     double rx2_ry2 = (ray(X) - ray(Y)) * (ray(X) + ray(Y));
     double sinrot, cosrot;
-    sincos(_rot_angle, sinrot, cosrot);
+    sincos(rotationAngle(), sinrot, cosrot);
     double expr1 = ray(X) * (p_c[X] * cosrot + p_c[Y] * sinrot);
     Poly coeff;
     coeff.resize(5);
@@ -662,227 +692,156 @@ std::vector<double> EllipticalArc::allNearestTimes( Point const& p, double from,
 }
 #endif
 
-/*
- * NOTE: this implementation follows Standard SVG 1.1 implementation guidelines
- * for elliptical arc curves. See Appendix F.6.
- */
-void EllipticalArc::_updateCenterAndAngles(bool svg)
-{
-    Point d = initialPoint() - finalPoint();
 
-    // TODO move this to SVGElipticalArc?
-    if (svg)
-    {
-        if ( initialPoint() == finalPoint() )
-        {
-            _rot_angle = _start_angle = _end_angle = 0;
-            _center = initialPoint();
-            _rays = Geom::Point(0,0);
-            _large_arc = _sweep = false;
-            return;
+void EllipticalArc::_filterIntersections(std::vector<ShapeIntersection> &xs, bool is_first) const
+{
+    Interval unit(0, 1);
+    std::vector<ShapeIntersection>::reverse_iterator i = xs.rbegin(), last = xs.rend();
+    while (i != last) {
+        Coord &t = is_first ? i->first : i->second;
+        assert(are_near(_ellipse.pointAt(t), i->point(), 1e-6));
+        t = timeAtAngle(t);
+        if (!unit.contains(t)) {
+            xs.erase((++i).base());
+            continue;
+        } else {
+            assert(are_near(pointAt(t), i->point(), 1e-6));
+            ++i;
         }
+    }
+}
 
-        _rays[X] = std::fabs(_rays[X]);
-        _rays[Y] = std::fabs(_rays[Y]);
-
-        if ( are_near(ray(X), 0) || are_near(ray(Y), 0) )
-        {
-            _rays[X] = L2(d) / 2;
-            _rays[Y] = 0;
-            _rot_angle = std::atan2(d[Y], d[X]);
-            _start_angle = 0;
-            _end_angle = M_PI;
-            _center = middle_point(initialPoint(), finalPoint());
-            _large_arc = false;
-            _sweep = false;
-            return;
-        }
+std::vector<CurveIntersection> EllipticalArc::intersect(Curve const &other, Coord eps) const
+{
+    if (isLineSegment()) {
+        LineSegment ls(_initial_point, _final_point);
+        return ls.intersect(other, eps);
     }
-    else
-    {
-        if ( are_near(initialPoint(), finalPoint()) )
-        {
-            if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )
-            {
-                _start_angle = _end_angle = 0;
-                _center = initialPoint();
-                return;
-            }
-            else
-            {
-                THROW_RANGEERROR("initial and final point are the same");
-            }
-        }
-        if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )
-        { // but initialPoint != finalPoint
-            THROW_RANGEERROR(
-                "there is no ellipse that satisfies the given constraints: "
-                "ray(X) == 0 && ray(Y) == 0 but initialPoint != finalPoint"
-            );
-        }
-        if ( are_near(ray(Y), 0) )
-        {
-            Point v = initialPoint() - finalPoint();
-            if ( are_near(L2sq(v), 4*ray(X)*ray(X)) )
-            {
-                Angle angle(v);
-                if ( are_near( angle, _rot_angle ) )
-                {
-                    _start_angle = 0;
-                    _end_angle = M_PI;
-                    _center = v/2 + finalPoint();
-                    return;
-                }
-                angle -= M_PI;
-                if ( are_near( angle, _rot_angle ) )
-                {
-                    _start_angle = M_PI;
-                    _end_angle = 0;
-                    _center = v/2 + finalPoint();
-                    return;
-                }
-                THROW_RANGEERROR(
-                    "there is no ellipse that satisfies the given constraints: "
-                    "ray(Y) == 0 "
-                    "and slope(initialPoint - finalPoint) != rotation_angle "
-                    "and != rotation_angle + PI"
-                );
-            }
-            if ( L2sq(v) > 4*ray(X)*ray(X) )
-            {
-                THROW_RANGEERROR(
-                    "there is no ellipse that satisfies the given constraints: "
-                    "ray(Y) == 0 and distance(initialPoint, finalPoint) > 2*ray(X)"
-                );
-            }
-            else
-            {
-                THROW_RANGEERROR(
-                    "there is infinite ellipses that satisfy the given constraints: "
-                    "ray(Y) == 0  and distance(initialPoint, finalPoint) < 2*ray(X)"
-                );
-            }
 
-        }
+    std::vector<CurveIntersection> result;
 
-        if ( are_near(ray(X), 0) )
-        {
-            Point v = initialPoint() - finalPoint();
-            if ( are_near(L2sq(v), 4*ray(Y)*ray(Y)) )
-            {
-                double angle = std::atan2(v[Y], v[X]);
-                if (angle < 0) angle += 2*M_PI;
-                double rot_angle = _rot_angle.radians() + M_PI/2;
-                if ( !(rot_angle < 2*M_PI) ) rot_angle -= 2*M_PI;
-                if ( are_near( angle, rot_angle ) )
-                {
-                    _start_angle = M_PI/2;
-                    _end_angle = 3*M_PI/2;
-                    _center = v/2 + finalPoint();
-                    return;
-                }
-                angle -= M_PI;
-                if ( angle < 0 ) angle += 2*M_PI;
-                if ( are_near( angle, rot_angle ) )
-                {
-                    _start_angle = 3*M_PI/2;
-                    _end_angle = M_PI/2;
-                    _center = v/2 + finalPoint();
-                    return;
-                }
-                THROW_RANGEERROR(
-                    "there is no ellipse that satisfies the given constraints: "
-                    "ray(X) == 0 "
-                    "and slope(initialPoint - finalPoint) != rotation_angle + PI/2 "
-                    "and != rotation_angle + (3/2)*PI"
-                );
-            }
-            if ( L2sq(v) > 4*ray(Y)*ray(Y) )
-            {
-                THROW_RANGEERROR(
-                    "there is no ellipse that satisfies the given constraints: "
-                    "ray(X) == 0 and distance(initialPoint, finalPoint) > 2*ray(Y)"
-                );
-            }
-            else
-            {
-                THROW_RANGEERROR(
-                    "there is infinite ellipses that satisfy the given constraints: "
-                    "ray(X) == 0  and distance(initialPoint, finalPoint) < 2*ray(Y)"
-                );
-            }
+    if (other.isLineSegment()) {
+        LineSegment ls(other.initialPoint(), other.finalPoint());
+        result = _ellipse.intersect(ls);
+        _filterIntersections(result, true);
+        return result;
+    }
 
-        }
+    BezierCurve const *bez = dynamic_cast<BezierCurve const *>(&other);
+    if (bez) {
+        result = _ellipse.intersect(bez->fragment());
+        _filterIntersections(result, true);
+        return result;
+    }
 
+    EllipticalArc const *arc = dynamic_cast<EllipticalArc const *>(&other);
+    if (arc) {
+        result = _ellipse.intersect(arc->_ellipse);
+        _filterIntersections(result, true);
+        arc->_filterIntersections(result, false);
+        return result;
     }
 
-    Rotate rm(_rot_angle);
-    Affine m(rm);
-    m[1] = -m[1];
-    m[2] = -m[2];
-
-    Point p = (d / 2) * m;
-    double rx2 = _rays[X] * _rays[X];
-    double ry2 = _rays[Y] * _rays[Y];
-    double rxpy = _rays[X] * p[Y];
-    double rypx = _rays[Y] * p[X];
-    double rx2py2 = rxpy * rxpy;
-    double ry2px2 = rypx * rypx;
-    double num = rx2 * ry2;
-    double den = rx2py2 + ry2px2;
-    assert(den != 0);
-    double rad = num / den;
-    Point c(0,0);
-    if (rad > 1)
-    {
-        rad -= 1;
-        rad = std::sqrt(rad);
+    // in case someone wants to make a custom curve type
+    result = other.intersect(*this, eps);
+    transpose_in_place(result);
+    return result;
+}
 
-        if (_large_arc == _sweep) rad = -rad;
-        c = rad * Point(rxpy / ray(Y), -rypx / ray(X));
 
-        _center = c * rm + middle_point(initialPoint(), finalPoint());
+void EllipticalArc::_updateCenterAndAngles()
+{
+    // See: http://www.w3.org/TR/SVG/implnote.html#ArcImplementationNotes
+    Point d = initialPoint() - finalPoint();
+    Point mid = middle_point(initialPoint(), finalPoint());
+
+    // if ip = sp, the arc contains no other points
+    if (initialPoint() == finalPoint()) {
+        _start_angle = _end_angle = 0;
+        _ellipse.setRotationAngle(0);
+        _ellipse.setCenter(initialPoint());
+        _ellipse.setRays(0, 0);
+        _large_arc = _sweep = false;
+        return;
     }
-    else if (rad == 1 || svg)
-    {
-        double lamda = std::sqrt(1 / rad);
-        _rays[X] *= lamda;
-        _rays[Y] *= lamda;
-        _center = middle_point(initialPoint(), finalPoint());
+
+    // 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;
+        _ellipse.setRays(L2(d) / 2, 0);
+        _ellipse.setRotationAngle(atan2(d));
+        _ellipse.setCenter(mid);
+        _large_arc = false;
+        _sweep = false;
+        return;
     }
-    else
-    {
-        THROW_RANGEERROR(
-            "there is no ellipse that satisfies the given constraints"
-        );
+
+    Rotate rot(rotationAngle()); // the matrix in F.6.5.3
+    Rotate invrot = rot.inverse(); // the matrix in F.6.5.1
+
+    Point r = rays();
+    Point p = (initialPoint() - mid) * invrot; // x', y' in F.6.5.1
+    Point c(0,0); // cx', cy' in F.6.5.2
+
+    // Correct out-of-range radii
+    Coord lambda = hypot(p[X]/r[X], p[Y]/r[Y]);
+    if (lambda > 1) {
+        r *= lambda;
+        _ellipse.setRays(r);
+        _ellipse.setCenter(mid);
+    } else {
+        // evaluate F.6.5.2
+        Coord rxry = r[X]*r[X] * r[Y]*r[Y];
+        Coord pxry = p[X]*p[X] * r[Y]*r[Y];
+        Coord rxpy = r[X]*r[X] * p[Y]*p[Y];
+        Coord rad = (rxry - pxry - rxpy)/(rxpy + pxry);
+        // normally rad should never be negative, but numerical inaccuracy may cause this
+        if (rad > 0) {
+            rad = std::sqrt(rad);
+            if (_sweep == _large_arc) {
+                rad = -rad;
+            }
+            c = rad * Point(r[X]*p[Y]/r[Y], -r[Y]*p[X]/r[X]);
+            _ellipse.setCenter(c * rot + mid);
+        } else {
+            _ellipse.setCenter(mid);
+        }
     }
 
-    Point sp((p[X] - c[X]) / ray(X), (p[Y] - c[Y]) / ray(Y));
-    Point ep((-p[X] - c[X]) / ray(X), (-p[Y] - c[Y]) / ray(Y));
+    // Compute start and end angles.
+    // If the ellipse was enlarged, c will be zero - this is correct.
+    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;
-    _end_angle += sweep_angle;
+    _end_angle = _start_angle + sweep_angle;
 }
 
 D2<SBasis> EllipticalArc::toSBasis() const
 {
+    if (isChord()) {
+        return chord().toSBasis();
+    }
+
     D2<SBasis> arc;
     // the interval of parametrization has to be [0,1]
-    Coord et = initialAngle().radians() + ( _sweep ? sweepAngle() : -sweepAngle() );
-    Linear param(initialAngle(), et);
-    Coord cos_rot_angle, sin_rot_angle;
-    sincos(_rot_angle, sin_rot_angle, cos_rot_angle);
+    Coord et = initialAngle().radians() + sweepAngle();
+    Linear param(initialAngle().radians(), et);
+    Coord cosrot, sinrot;
+    sincos(rotationAngle(), sinrot, cosrot);
 
     // order = 4 seems to be enough to get a perfect looking elliptical arc
     SBasis arc_x = ray(X) * cos(param,4);
     SBasis arc_y = ray(Y) * sin(param,4);
-    arc[0] = arc_x * cos_rot_angle - arc_y * sin_rot_angle + Linear(center(X),center(X));
-    arc[1] = arc_x * sin_rot_angle + arc_y * cos_rot_angle + Linear(center(Y),center(Y));
+    arc[0] = arc_x * cosrot - arc_y * sinrot + Linear(center(X), center(X));
+    arc[1] = arc_x * sinrot + arc_y * cosrot + Linear(center(Y), center(Y));
 
     // ensure that endpoints remain exact
     for ( int d = 0 ; d < 2 ; d++ ) {
@@ -898,22 +857,24 @@ void EllipticalArc::transform(Affine const& m)
     if (isChord()) {
         _initial_point *= m;
         _final_point *= m;
-        _center = middle_point(_initial_point, _final_point);
-        _rays = Point(0,0);
-        _rot_angle = 0;
+        _ellipse.setCenter(middle_point(_initial_point, _final_point));
+        _ellipse.setRays(0, 0);
+        _ellipse.setRotationAngle(0);
         return;
     }
 
-    // TODO avoid allocating a new arc here
-    Ellipse e(center(X), center(Y), ray(X), ray(Y), _rot_angle);
-    e *= m;
-    Point inner_point = pointAt(0.5);
-    EllipticalArc *arc = e.arc(initialPoint() * m,
-                               inner_point * m,
-                               finalPoint() * m,
-                               isSVGCompliant());
-    *this = *arc;
-    delete arc;
+    _initial_point *= m;
+    _final_point *= m;
+    _ellipse *= m;
+    if (m.det() < 0) {
+        _sweep = !_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);
 }
 
 bool EllipticalArc::operator==(Curve const &c) const
@@ -924,8 +885,8 @@ bool EllipticalArc::operator==(Curve const &c) const
     if (_final_point != other->_final_point) return false;
     // TODO: all arcs with ellipse rays which are too small
     //       and fall back to a line should probably be equal
-    if (_rays != other->_rays) return false;
-    if (_rot_angle != other->_rot_angle) return false;
+    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;
     return true;
@@ -936,7 +897,7 @@ void EllipticalArc::feed(PathSink &sink, bool moveto_initial) const
     if (moveto_initial) {
         sink.moveTo(_initial_point);
     }
-    sink.arcTo(_rays[X], _rays[Y], _rot_angle, _large_arc, _sweep, _final_point);
+    sink.arcTo(ray(X), ray(Y), rotationAngle(), _large_arc, _sweep, _final_point);
 }
 
 Coord EllipticalArc::map_to_01(Coord angle) const
@@ -945,6 +906,19 @@ Coord EllipticalArc::map_to_01(Coord angle) const
                                              finalAngle(), _sweep);
 }
 
+
+std::ostream &operator<<(std::ostream &out, EllipticalArc const &ea)
+{
+    out << "EllipticalArc("
+        << ea.initialPoint() << ", "
+        << format_coord_nice(ea.ray(X)) << ", " << format_coord_nice(ea.ray(Y)) << ", "
+        << format_coord_nice(ea.rotationAngle()) << ", "
+        << "large_arc=" << (ea.largeArc() ? "true" : "false") << ", "
+        << "sweep=" << (ea.sweep() ? "true" : "false") << ", "
+        << ea.finalPoint() << ")";
+    return out;
+}
+
 } // end namespace Geom
 
 /*
diff --git a/src/2geom/elliptical-arc.h b/src/2geom/elliptical-arc.h
index 342b8c8c9..83909370e 100644
--- a/src/2geom/elliptical-arc.h
+++ b/src/2geom/elliptical-arc.h
@@ -38,10 +38,11 @@
 #define LIB2GEOM_SEEN_ELLIPTICAL_ARC_H
 
 #include <algorithm>
+#include <2geom/affine.h>
 #include <2geom/angle.h>
 #include <2geom/bezier-curve.h>
 #include <2geom/curve.h>
-#include <2geom/affine.h>
+#include <2geom/ellipse.h>
 #include <2geom/sbasis-curve.h>  // for non-native methods
 #include <2geom/utils.h>
 
@@ -56,21 +57,30 @@ public:
         : AngleInterval(0, 0, true)
         , _initial_point(0,0)
         , _final_point(0,0)
-        , _rays(0,0)
-        , _center(0,0)
-        , _rot_angle(0)
         , _large_arc(true)
     {}
     /** @brief Create a new elliptical arc.
      * @param ip Initial point of the arc
-     * @param rx First ray of the ellipse
-     * @param ry Second ray of the ellipse
+     * @param r Rays of the ellipse as a point
      * @param rot Angle of rotation of the X axis of the ellipse in radians
      * @param large If true, the large arc is chosen (always >= 180 degrees), otherwise
      *              the smaller arc is chosen
      * @param sweep If true, the clockwise arc is chosen, otherwise the counter-clockwise
      *              arc is chosen
      * @param fp Final point of the arc */
+    EllipticalArc( Point const &ip, Point const &r,
+                   Coord rot_angle, bool large_arc, bool sweep,
+                   Point const &fp
+                 )
+        : AngleInterval(0,0,sweep)
+        , _initial_point(ip)
+        , _final_point(fp)
+        , _ellipse(0, 0, r[X], r[Y], rot_angle)
+        , _large_arc(large_arc)
+    {
+        _updateCenterAndAngles();
+    }
+
     EllipticalArc( Point const &ip, Coord rx, Coord ry,
                    Coord rot_angle, bool large_arc, bool sweep,
                    Point const &fp
@@ -78,11 +88,10 @@ public:
         : AngleInterval(0,0,sweep)
         , _initial_point(ip)
         , _final_point(fp)
-        , _rays(rx, ry)
-        , _rot_angle(rot_angle)
+        , _ellipse(0, 0, rx, ry, rot_angle)
         , _large_arc(large_arc)
     {
-        _updateCenterAndAngles(false);
+        _updateCenterAndAngles();
     }
 
     // methods new to EllipticalArc go here
@@ -98,15 +107,15 @@ public:
     /** @brief Get the defining ellipse's rotation
      * @return Angle between the +X ray of the ellipse and the +X axis */
     Angle rotationAngle() const {
-        return _rot_angle;
+        return _ellipse.rotationAngle();
     }
     /** @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 _rays[d]; }
+    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 _rays; }
+    Point rays() const { return _ellipse.rays(); }
     /** @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; }
@@ -125,12 +134,11 @@ public:
     {
         _initial_point = ip;
         _final_point = fp;
-        _rays[X] = rx;
-        _rays[Y] = ry;
-        _rot_angle = Angle(rot_angle);
+        _ellipse.setRays(rx, ry);
+        _ellipse.setRotationAngle(rot_angle);
         _large_arc = large_arc;
         _sweep = sweep;
-        _updateCenterAndAngles(isSVGCompliant());
+        _updateCenterAndAngles();
     }
     /** @brief Change the initial and final point in one operation.
      * This method exists because modifying any of the endpoints causes rather costly
@@ -140,21 +148,16 @@ public:
     void setEndpoints(Point const &ip, Point const &fp) {
         _initial_point = ip;
         _final_point = fp;
-        _updateCenterAndAngles(isSVGCompliant());
+        _updateCenterAndAngles();
     }
     /// @}
 
     /// @name Access computed parameters of the arc
     /// @{
-    Coord center(Dim2 d) const { return _center[d]; }
+    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 _center; }
-    /** @brief Get the extent of the arc
-     * @return The angle between the initial and final point, in arc's angular coordinates */
-    Coord sweepAngle() const {
-        return extent();
-    }
+    Point center() const { return _ellipse.center(); }
     /// @}
     
     /// @name Angular evaluation
@@ -177,14 +180,13 @@ public:
      * 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;
     /// @}
 
-    /** @brief Check whether the arc adheres to SVG 1.1 implementation guidelines */
-    virtual bool isSVGCompliant() const { return false; }
-
-    /// Check whether both rays are nonzero
+    /** @brief Check whether both rays are nonzero.
+     * If they are not, the arc is represented as a line segment instead. */
     bool isChord() const {
-        return _rays[0] == 0 || _rays[Y] == 0;
+        return ray(X) == 0 || ray(Y) == 0;
     }
 
     std::pair<EllipticalArc, EllipticalArc> subdivide(Coord t) const {
@@ -203,15 +205,16 @@ public:
     virtual Curve* duplicate() const { return new EllipticalArc(*this); }
     virtual void setInitial(Point const &p) {
         _initial_point = p;
-        _updateCenterAndAngles(isSVGCompliant());
+        _updateCenterAndAngles();
     }
     virtual void setFinal(Point const &p) {
         _final_point = p;
-        _updateCenterAndAngles(isSVGCompliant());
+        _updateCenterAndAngles();
     }
     virtual bool isDegenerate() const {
         return _initial_point == _final_point;
     }
+    virtual bool isLineSegment() const { return isChord(); }
     virtual Rect boundsFast() const {
         return boundsExact();
     }
@@ -230,13 +233,14 @@ public:
         }
         return allNearestTimes(p, from, to).front();
     }
+    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.vector();
-        _final_point += m.vector();
-        _center += m.vector();
+        _initial_point *= m;
+        _final_point *= m;
+        _ellipse *= m;
         return *this;
     }
 
@@ -248,28 +252,34 @@ public:
 
     virtual D2<SBasis> toSBasis() const;
     virtual double valueAt(Coord t, Dim2 d) const {
-    	return valueAtAngle(angleAt(t), d);
-    }
-    virtual Point pointAt(Coord t) const {
-        return pointAtAngle(angleAt(t));
+        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 bool operator==(Curve const &c) const;
     virtual void feed(PathSink &sink, bool moveto_initial) const;
 
 protected:
-    void _updateCenterAndAngles(bool svg);
+    void _updateCenterAndAngles();
+    void _filterIntersections(std::vector<ShapeIntersection> &xs, bool is_first) const;
 
     Point _initial_point, _final_point;
-    Point _rays, _center;
-    Angle _rot_angle;
+    Ellipse _ellipse;
     bool _large_arc;
 
 private:
     Coord map_to_01(Coord angle) const; 
 }; // end class EllipticalArc
 
+
+// implemented in elliptical-arc-from-sbasis.cpp
+bool arc_from_sbasis(EllipticalArc &ea, D2<SBasis> const &in,
+                     double tolerance = EPSILON, unsigned num_samples = 20);
+
+std::ostream &operator<<(std::ostream &out, EllipticalArc const &ea);
+
 } // end namespace Geom
 
 #endif // LIB2GEOM_SEEN_ELLIPTICAL_ARC_H
diff --git a/src/2geom/forward.h b/src/2geom/forward.h
index a6e420122..cf9c75988 100644
--- a/src/2geom/forward.h
+++ b/src/2geom/forward.h
@@ -37,7 +37,7 @@
 
 namespace Geom {
 
-// basic types
+// primitives
 typedef double Coord;
 typedef int IntCoord;
 class Point;
@@ -61,6 +61,12 @@ typedef GenericOptRect<IntCoord> OptIntRect;
 class Linear;
 class Bezier;
 class SBasis;
+class Poly;
+
+// shapes
+class Circle;
+class Ellipse;
+class ConvexHull;
 
 // curves
 class Curve;
@@ -73,7 +79,6 @@ typedef BezierCurveN<1> LineSegment;
 typedef BezierCurveN<2> QuadraticBezier;
 typedef BezierCurveN<3> CubicBezier;
 class EllipticalArc;
-class SVGEllipticalArc;
 
 // paths and path sequences
 class Path;
diff --git a/src/2geom/intersection-graph.cpp b/src/2geom/intersection-graph.cpp
index b9e2feeed..ff96c28a8 100644
--- a/src/2geom/intersection-graph.cpp
+++ b/src/2geom/intersection-graph.cpp
@@ -35,6 +35,7 @@
 #include <2geom/path.h>
 #include <2geom/pathvector.h>
 #include <iostream>
+#include <iterator>
 
 namespace Geom {
 
@@ -157,6 +158,41 @@ PathVector PathIntersectionGraph::getIntersection()
     return result;
 }
 
+PathVector PathIntersectionGraph::getAminusB()
+{
+    PathVector result = _getResult(false, true);
+    _handleNonintersectingPaths(result, 0, false);
+    return result;
+}
+
+PathVector PathIntersectionGraph::getBminusA()
+{
+    PathVector result = _getResult(true, false);
+    _handleNonintersectingPaths(result, 1, false);
+    return result;
+}
+
+PathVector PathIntersectionGraph::getXOR()
+{
+    PathVector r1 = getAminusB();
+    PathVector r2 = getBminusA();
+    std::copy(r2.begin(), r2.end(), std::back_inserter(r1));
+    return r1;
+}
+
+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) {
+            result.push_back(j->p);
+        }
+    }
+    return result;
+}
+
 PathVector PathIntersectionGraph::_getResult(bool enter_a, bool enter_b)
 {
     typedef IntersectionList::iterator Iter;
@@ -231,6 +267,7 @@ PathVector PathIntersectionGraph::_getResult(bool enter_a, bool enter_b)
             std::swap(lscur, lsother);
             std::swap(cur, other);
         }
+        result.back().close(true);
 
         assert(!result.back().empty());
     }
diff --git a/src/2geom/intersection-graph.h b/src/2geom/intersection-graph.h
index dd70d3d5b..44beaf46b 100644
--- a/src/2geom/intersection-graph.h
+++ b/src/2geom/intersection-graph.h
@@ -79,6 +79,11 @@ public:
 
     PathVector getUnion();
     PathVector getIntersection();
+    PathVector getAminusB();
+    PathVector getBminusA();
+    PathVector getXOR();
+
+    std::vector<Point> intersectionPoints() const;
 
 private:
     PathVector _getResult(bool enter_a, bool enter_b);
diff --git a/src/2geom/intersection.h b/src/2geom/intersection.h
index 848797682..ae4b50dd1 100644
--- a/src/2geom/intersection.h
+++ b/src/2geom/intersection.h
@@ -92,7 +92,7 @@ private:
 };
 
 
-// TODO: move into new header
+// TODO: move into new header?
 template <typename T>
 struct ShapeTraits {
     typedef Coord TimeType;
@@ -101,6 +101,24 @@ struct ShapeTraits {
     typedef Intersection<> IntersectionType;
 };
 
+template <typename A, typename B> inline
+std::vector< Intersection<A, B> > transpose(std::vector< Intersection<B, A> > const &in) {
+    std::vector< Intersection<A, B> > result;
+    for (std::size_t i = 0; i < in.size(); ++i) {
+        result.push_back(Intersection<A, B>(in[i].second, in[i].first, in[i].point()));
+    }
+    return result;
+}
+
+template <typename T> inline
+void transpose_in_place(std::vector< Intersection<T, T> > &xs) {
+    for (std::size_t i = 0; i < xs.size(); ++i) {
+        std::swap(xs[i].first, xs[i].second);
+    }
+}
+
+typedef Intersection<> ShapeIntersection;
+
 
 } // namespace Geom
 
diff --git a/src/2geom/interval.h b/src/2geom/interval.h
index 09ea46923..fd446a1c6 100644
--- a/src/2geom/interval.h
+++ b/src/2geom/interval.h
@@ -93,7 +93,7 @@ public:
 
     /// @name Inspect contained values.
     /// @{
-    /** @brief Check whether both endpoints are finite. */
+    /// Check whether both endpoints are finite.
     bool isFinite() const {
         return IS_FINITE(min()) && IS_FINITE(max());
     }
@@ -102,11 +102,16 @@ public:
     Coord valueAt(Coord t) {
         return lerp(t, min(), max());
     }
-    /** @brief Find closest time in [0,1] that maps to the given value. */
-    Coord nearestTime(Coord t) {
-        if (t < min()) return 0;
-        if (t > max()) return 1;
-        return (t - min()) / extent();
+    /** @brief Compute a time value that maps to the given value.
+     * The supplied value does not need to be in the interval for this method to work. */
+    Coord timeAt(Coord v) {
+        return (v - min()) / extent();
+    }
+    /// Find closest time in [0,1] that maps to the given value. */
+    Coord nearestTime(Coord v) {
+        if (v <= min()) return 0;
+        if (v >= max()) return 1;
+        return timeAt(v);
     }
     /// @}
 
diff --git a/src/2geom/line.cpp b/src/2geom/line.cpp
index 35cf2d379..8bea33638 100644
--- a/src/2geom/line.cpp
+++ b/src/2geom/line.cpp
@@ -62,7 +62,7 @@ void Line::setCoefficients (Coord a, Coord b, Coord c)
 {
     if (a == 0 && b == 0) {
         if (c != 0) {
-            THROW_LOGICALERROR("the passed coefficients gives the empty set");
+            THROW_LOGICALERROR("the passed coefficients give the empty set");
         }
         _initial = Point(0,0);
         _final = Point(0,0);
@@ -70,20 +70,20 @@ void Line::setCoefficients (Coord a, Coord b, Coord c)
     }
     if (a == 0) {
         // b must be nonzero
-        _initial = Point(0, c / b);
+        _initial = Point(0, -c / b);
         _final = _initial;
         _final[X] = 1;
         return;
     }
     if (b == 0) {
-        _initial = Point(c / a, 0);
+        _initial = Point(-c / a, 0);
         _final = _initial;
         _final[Y] = 1;
         return;
     }
 
-    _initial = Point(c / a, 0);
-    _final = Point(0, c / b);
+    _initial = Point(-c / a, 0);
+    _final = Point(0, -c / b);
 }
 
 void Line::coefficients(Coord &a, Coord &b, Coord &c) const
@@ -217,6 +217,66 @@ Coord Line::timeAt(Point const &p) const
     }
 }
 
+std::vector<ShapeIntersection> Line::intersect(Line const &other) const
+{
+    std::vector<ShapeIntersection> result;
+
+    Point v1 = versor();
+    Point v2 = other.versor();
+    Coord cp = cross(v1, v2);
+    if (cp == 0) return result;
+
+    Point odiff = other.initialPoint() - initialPoint();
+    Coord t1 = cross(odiff, v2) / cp;
+    Coord t2 = cross(odiff, v1) / cp;
+    result.push_back(ShapeIntersection(*this, other, t1, t2));
+    return result;
+}
+
+std::vector<ShapeIntersection> Line::intersect(Ray const &r) const
+{
+    Line other(r);
+    std::vector<ShapeIntersection> result = intersect(other);
+    filter_ray_intersections(result, false, true);
+    return result;
+}
+
+std::vector<ShapeIntersection> Line::intersect(LineSegment const &ls) const
+{
+    Line other(ls);
+    std::vector<ShapeIntersection> result = intersect(other);
+    filter_line_segment_intersections(result, false, true);
+    return result;
+}
+
+
+
+void filter_line_segment_intersections(std::vector<ShapeIntersection> &xs, bool a, bool b)
+{
+    Interval unit(0, 1);
+    std::vector<ShapeIntersection>::reverse_iterator i = xs.rbegin(), last = xs.rend();
+    while (i != last) {
+        if ((a && !unit.contains(i->first)) || (b && !unit.contains(i->second))) {
+            xs.erase((++i).base());
+        } else {
+            ++i;
+        }
+    }
+}
+
+void filter_ray_intersections(std::vector<ShapeIntersection> &xs, bool a, bool b)
+{
+    Interval unit(0, 1);
+    std::vector<ShapeIntersection>::reverse_iterator i = xs.rbegin(), last = xs.rend();
+    while (i != last) {
+        if ((a && i->first < 0) || (b && i->second < 0)) {
+            xs.erase((++i).base());
+        } else {
+            ++i;
+        }
+    }
+}
+
 namespace detail
 {
 
diff --git a/src/2geom/line.h b/src/2geom/line.h
index 2c47b4486..5d92c436d 100644
--- a/src/2geom/line.h
+++ b/src/2geom/line.h
@@ -49,9 +49,15 @@ namespace Geom
 
 // class docs in cpp file
 class Line
-    : boost::equality_comparable1<Line
-    , MultipliableNoncommutative<Line, Affine
-      > >
+    : 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
+      > > > > > > > >
 {
 private:
     Point _initial;
@@ -188,6 +194,14 @@ public:
     bool isDegenerate() const {
         return _initial == _final;
     }
+    /// Check if the line is horizontal (y is constant).
+    bool isHorizontal() const {
+        return _initial[Y] == _final[Y];
+    }
+    /// Check if the line is vertical (x is constant).
+    bool isVertical() const {
+        return _initial[X] == _final[X];
+    }
 
     /** @brief Reparametrize the line so that it has unit speed. */
     void normalize() {
@@ -349,13 +363,18 @@ public:
     }
     /// @}
 
-    //std::vector<LineIntersection> intersect(Line const &other, Coord precision = EPSILON) const;
+    std::vector<ShapeIntersection> intersect(Line const &other) const;
+    std::vector<ShapeIntersection> intersect(Ray const &r) const;
+    std::vector<ShapeIntersection> intersect(LineSegment const &ls) const;
 
-    Line &operator*=(Affine const &m) {
-        _initial *= m;
-        _final *= m;
+    template <typename T>
+    Line &operator*=(T const &tr) {
+        BOOST_CONCEPT_ASSERT((TransformConcept<T>));
+        _initial *= tr;
+        _final *= tr;
         return *this;
     }
+
     bool operator==(Line const &other) const {
         if (distance(pointAt(nearestTime(other._initial)), other._initial) != 0) return false;
         if (distance(pointAt(nearestTime(other._final)), other._final) != 0) return false;
@@ -363,6 +382,15 @@ public:
     }
 }; // end class Line
 
+/** @brief Removes intersections outside of the unit interval.
+ * A helper used to implement line segment intersections.
+ * @param xs Line intersections
+ * @param a Whether the first time value has to be in the unit interval
+ * @param b Whether the second time value has to be in the unit interval
+ * @return Appropriately filtered intersections */
+void filter_line_segment_intersections(std::vector<ShapeIntersection> &xs, bool a=false, bool b=true);
+void filter_ray_intersections(std::vector<ShapeIntersection> &xs, bool a=false, bool b=true);
+
 /// @brief Compute distance from point to line.
 /// @relates Line
 inline
diff --git a/src/2geom/numeric/fitting-model.h b/src/2geom/numeric/fitting-model.h
index 5cc6dde7a..fb96d1d2a 100644
--- a/src/2geom/numeric/fitting-model.h
+++ b/src/2geom/numeric/fitting-model.h
@@ -39,7 +39,7 @@
 #include <2geom/sbasis.h>
 #include <2geom/bezier.h>
 #include <2geom/bezier-curve.h>
-#include <2geom/poly.h>
+#include <2geom/polynomial.h>
 #include <2geom/ellipse.h>
 #include <2geom/circle.h>
 #include <2geom/utils.h>
diff --git a/src/2geom/path-intersection.h b/src/2geom/path-intersection.h
index e3a8d1de3..f06eeaf94 100644
--- a/src/2geom/path-intersection.h
+++ b/src/2geom/path-intersection.h
@@ -35,11 +35,9 @@
 #ifndef LIB2GEOM_SEEN_PATH_INTERSECTION_H
 #define LIB2GEOM_SEEN_PATH_INTERSECTION_H
 
-#include <2geom/path.h>
-
 #include <2geom/crossing.h>
-
-#include <2geom/sweep.h>
+#include <2geom/path.h>
+#include <2geom/sweep-bounds.h>
 
 namespace Geom {
 
diff --git a/src/2geom/path-sink.cpp b/src/2geom/path-sink.cpp
index 6ccc21e7b..3b8d407f8 100644
--- a/src/2geom/path-sink.cpp
+++ b/src/2geom/path-sink.cpp
@@ -31,6 +31,8 @@
 #include <2geom/sbasis-to-bezier.h>
 #include <2geom/path-sink.h>
 #include <2geom/exception.h>
+#include <2geom/circle.h>
+#include <2geom/ellipse.h>
 
 namespace Geom {
 
@@ -68,6 +70,26 @@ void PathSink::feed(Rect const &r) {
     closePath();
 }
 
+void PathSink::feed(Circle const &e) {
+    Coord r = e.radius();
+    Point c = e.center();
+    Point a = c + Point(0, c[Y] + r);
+    Point b = c + Point(0, c[Y] - r);
+
+    moveTo(a);
+    arcTo(r, r, 0, false, false, b);
+    arcTo(r, r, 0, false, false, a);
+    closePath();
+}
+
+void PathSink::feed(Ellipse const &e) {
+    Point s = e.pointAt(0);
+    moveTo(s);
+    arcTo(e.ray(X), e.ray(Y), e.rotationAngle(), false, false, e.pointAt(M_PI));
+    arcTo(e.ray(X), e.ray(Y), e.rotationAngle(), false, false, s);
+    closePath();
+}
+
 }
 
 /*
diff --git a/src/2geom/path-sink.h b/src/2geom/path-sink.h
index d6806031d..17ede18a4 100644
--- a/src/2geom/path-sink.h
+++ b/src/2geom/path-sink.h
@@ -32,6 +32,7 @@
 #ifndef LIB2GEOM_SEEN_PATH_SINK_H
 #define LIB2GEOM_SEEN_PATH_SINK_H
 
+#include <2geom/forward.h>
 #include <2geom/pathvector.h>
 #include <2geom/curves.h>
 #include <iterator>
@@ -101,6 +102,10 @@ public:
     virtual void feed(PathVector const &v);
     /// Output an axis-aligned rectangle, using moveTo, lineTo and closePath.
     virtual void feed(Rect const &);
+    /// Output a circle as two elliptical arcs.
+    virtual void feed(Circle const &e);
+    /// Output an ellipse as two elliptical arcs.
+    virtual void feed(Ellipse const &e);
 
     virtual ~PathSink() {}
 };
@@ -152,8 +157,8 @@ public:
         if (!_in_path) {
             moveTo(_start_p);
         }
-        _path.template appendNew<SVGEllipticalArc>(rx, ry, angle,
-                                                 large_arc, sweep, p);
+        _path.template appendNew<EllipticalArc>(rx, ry, angle,
+                                                large_arc, sweep, p);
     }
 
     bool backspace()
diff --git a/src/2geom/path.cpp b/src/2geom/path.cpp
index 8cba316f5..5fbc65b3e 100644
--- a/src/2geom/path.cpp
+++ b/src/2geom/path.cpp
@@ -35,8 +35,11 @@
 #include <2geom/path.h>
 #include <2geom/pathvector.h>
 #include <2geom/transforms.h>
+#include <2geom/circle.h>
+#include <2geom/ellipse.h>
 #include <2geom/convex-hull.h>
 #include <2geom/svg-path-writer.h>
+#include <2geom/sweeper.h>
 #include <algorithm>
 #include <limits>
 
@@ -231,7 +234,7 @@ Path::Path(Rect const &r)
 Path::Path(ConvexHull const &ch)
     : _curves(new Sequence())
     , _closing_seg(new ClosingSegment(Point(), Point()))
-    , _closed(false)
+    , _closed(true)
     , _exception_on_stitch(true)
 {
     if (ch.empty()) {
@@ -253,6 +256,49 @@ Path::Path(ConvexHull const &ch)
     _closed = true;
 }
 
+Path::Path(Circle const &c)
+    : _curves(new Sequence())
+    , _closing_seg(NULL)
+    , _closed(true)
+    , _exception_on_stitch(true)
+{
+    Point p1 = c.pointAt(0);
+    Point p2 = c.pointAt(M_PI);
+    _curves->push_back(new EllipticalArc(p1, c.radius(), c.radius(), 0, false, true, p2));
+    _curves->push_back(new EllipticalArc(p2, c.radius(), c.radius(), 0, false, true, p1));
+    _closing_seg = new ClosingSegment(p1, p1);
+    _curves->push_back(_closing_seg);
+}
+
+Path::Path(Ellipse const &e)
+    : _curves(new Sequence())
+    , _closing_seg(NULL)
+    , _closed(true)
+    , _exception_on_stitch(true)
+{
+    Point p1 = e.pointAt(0);
+    Point p2 = e.pointAt(M_PI);
+    _curves->push_back(new EllipticalArc(p1, e.rays(), e.rotationAngle(), false, true, p2));
+    _curves->push_back(new EllipticalArc(p2, e.rays(), e.rotationAngle(), false, true, p1));
+    _closing_seg = new ClosingSegment(p1, p1);
+    _curves->push_back(_closing_seg);
+}
+
+void Path::close(bool c)
+{
+    if (c == _closed) return;
+    if (c && _curves->size() >= 2) {
+        // when closing, if last segment is linear and ends at initial point,
+        // replace it with the closing segment
+        Sequence::iterator last = _curves->end() - 2;
+        if (last->isLineSegment() && last->finalPoint() == initialPoint()) {
+            _closing_seg->setInitial(last->initialPoint());
+            _curves->erase(last);
+        }
+    }
+    _closed = c;
+}
+
 void Path::clear()
 {
     _unshare();
@@ -405,23 +451,91 @@ std::vector<PathTime> Path::roots(Coord v, Dim2 d) const
     return res;
 }
 
-std::vector<PathIntersection> Path::intersect(Path const &other, Coord precision) const
+
+// The class below implements sweepline optimization for curve intersection in paths.
+// Instead of O(N^2), this takes O(N + X), where X is the number of overlaps
+// between the bounding boxes of curves.
+namespace {
+
+struct CurveSweepTraits {
+    struct Bound {
+        Rect r;
+        std::size_t index;
+        int which;
+    };
+    typedef std::less<Coord> Compare;
+    inline static Coord entry_value(Bound const &b) { return b.r[X].min(); }
+    inline static Coord exit_value(Bound const &b) { return b.r[X].max(); }
+};
+
+class CurveSweeper
+    : public Sweeper<Curve const *, CurveSweepTraits>
 {
-    std::vector<PathIntersection> result;
+public:
+    CurveSweeper(Path const &a, Path const &b, std::vector<PathIntersection> &result, Coord prec)
+        : _result(result)
+        , _precision(prec)
+    {
+        for (std::size_t i = 0; i < a.size(); ++i) {
+            Bound bound;
+            bound.r = a[i].boundsFast();
+            bound.index = i;
+            bound.which = 0;
+            insert(bound, &a[i]);
+        }
+        for (std::size_t i = 0; i < b.size(); ++i) {
+            Bound bound;
+            bound.r = b[i].boundsFast();
+            bound.index = i;
+            bound.which = 1;
+            insert(bound, &b[i]);
+        }
+    }
+
+protected:
+    void _enter(Record const &record) {
+        int which = record.bound.which;
+
+        for (RecordList::iterator i = _active_items.begin(); i != _active_items.end(); ++i) {
+            // do not intersect in the same path
+            if (i->bound.which == which) continue;
+            // do not intersect if boxes do not overlap in Y
+            if (!record.bound.r[Y].intersects(i->bound.r[Y])) continue;
+
+            std::vector<CurveIntersection> cx;
+            int ia = record.bound.index;
+            int ib = i->bound.index;
+
+            if (which == 0) {
+                cx = record.item->intersect(*i->item);
+            } else {
+                cx = i->item->intersect(*record.item, _precision);
+                std::swap(ia, ib);
+            }
 
-    // TODO: sweepline optimization
-    // TODO: remove multiple intersections within precision of each other?
-    for (size_type i = 0; i < size(); ++i) {
-        for (size_type j = 0; j < other.size(); ++j) {
-            std::vector<CurveIntersection> cx = (*this)[i].intersect(other[j], precision);
             for (std::size_t ci = 0; ci < cx.size(); ++ci) {
-                PathTime a(i, cx[ci].first), b(j, cx[ci].second);
+                PathTime a(ia, cx[ci].first), b(ib, cx[ci].second);
                 PathIntersection px(a, b, cx[ci].point());
-                result.push_back(px);
+                _result.push_back(px);
             }
         }
     }
 
+private:
+    std::vector<PathIntersection> &_result;
+    Coord _precision;
+};
+
+} // end anonymous namespace
+
+std::vector<PathIntersection> Path::intersect(Path const &other, Coord precision) const
+{
+    std::vector<PathIntersection> result;
+
+    CurveSweeper sweeper(*this, other, result, precision);
+    sweeper.process();
+
+    // TODO: remove multiple intersections within precision of each other?
     return result;
 }
 
@@ -700,11 +814,10 @@ Path Path::reversed() const
 
     if (_closed) {
         // when the path is closed, there are two cases:
-        BezierCurve const *iseg = dynamic_cast<BezierCurve const *>(&front());
-        if (iseg && iseg->size() == 2) {
+        if (front().isLineSegment()) {
             // 1. initial segment is linear: it becomes the new closing segment.
             rend = RIter(_curves->begin() + 1);
-            ret._closing_seg = new ClosingSegment(iseg->finalPoint(), iseg->initialPoint());
+            ret._closing_seg = new ClosingSegment(front().finalPoint(), front().initialPoint());
         } else {
             // 2. initial segment is not linear: the closing segment becomes degenerate.
             // However, skip it if it's already degenerate.
@@ -883,6 +996,13 @@ void Path::do_append(Curve *c)
         if (c->initialPoint() != _closing_seg->initialPoint()) {
             THROW_CONTINUITYERROR();
         }
+        if (_closed && c->isLineSegment() &&
+            c->finalPoint() == _closing_seg->finalPoint())
+        {
+            // appending a curve that matches the closing segment has no effect
+            delete c;
+            return;
+        }
     }
     _curves->insert(_curves->end() - 1, c);
     _closing_seg->setInitial(c->finalPoint());
diff --git a/src/2geom/path.h b/src/2geom/path.h
index 5340df0c5..a6473e0d2 100644
--- a/src/2geom/path.h
+++ b/src/2geom/path.h
@@ -371,11 +371,14 @@ public:
         _curves->push_back(_closing_seg);
     }
 
-    /// Construct a path from a rectangle.
-    Path(Rect const &r);
-
-    /// Construct a path from a convex hull.
-    Path(ConvexHull const &);
+    /// Create a path from a rectangle.
+    explicit Path(Rect const &r);
+    /// Create a path from a convex hull.
+    explicit Path(ConvexHull const &);
+    /// Create a path from a circle, using two elliptical arcs.
+    explicit Path(Circle const &c);
+    /// Create a path from an ellipse, using two elliptical arcs.
+    explicit Path(Ellipse const &e);
 
     virtual ~Path() {}
 
@@ -463,8 +466,12 @@ public:
     /// Check whether the path is closed.
     bool closed() const { return _closed; }
 
-    /// Set whether the path is closed.
-    void close(bool closed = true) { _closed = closed; }
+    /** @brief Set whether the path is closed.
+     * When closing a path where the last segment can be represented as a closing
+     * segment, the last segment will be removed. When opening a path, the closing
+     * segment will be erased. This means that closing and then opening a path
+     * will not always give back the original path. */
+    void close(bool closed = true);
 
     /** @brief Remove all curves from the path.
      * The initial and final points of the closing segment are set to (0,0).
diff --git a/src/2geom/pathvector.cpp b/src/2geom/pathvector.cpp
index 720428e97..bff201c71 100644
--- a/src/2geom/pathvector.cpp
+++ b/src/2geom/pathvector.cpp
@@ -31,10 +31,10 @@
  * the specific language governing rights and limitations.
  */
 
-#include <2geom/pathvector.h>
-
-#include <2geom/path.h>
 #include <2geom/affine.h>
+#include <2geom/path.h>
+#include <2geom/pathvector.h>
+#include <2geom/svg-path-writer.h>
 
 namespace Geom {
 
@@ -221,6 +221,14 @@ PathVectorTime PathVector::_factorTime(Coord t) const
     return ret;
 }
 
+std::ostream &operator<<(std::ostream &out, PathVector const &pv)
+{
+    SVGPathWriter wr;
+    wr.feed(pv);
+    out << wr.str();
+    return out;
+}
+
 } // namespace Geom
 
 /*
diff --git a/src/2geom/pathvector.h b/src/2geom/pathvector.h
index 955cd1d37..6636cbf2e 100644
--- a/src/2geom/pathvector.h
+++ b/src/2geom/pathvector.h
@@ -276,6 +276,8 @@ private:
 inline OptRect bounds_fast(PathVector const &pv) { return pv.boundsFast(); }
 inline OptRect bounds_exact(PathVector const &pv) { return pv.boundsExact(); }
 
+std::ostream &operator<<(std::ostream &out, PathVector const &pv);
+
 } // end namespace Geom
 
 #endif // LIB2GEOM_SEEN_PATHVECTOR_H
diff --git a/src/2geom/point.cpp b/src/2geom/point.cpp
index e2662becd..bcdbab429 100644
--- a/src/2geom/point.cpp
+++ b/src/2geom/point.cpp
@@ -35,6 +35,7 @@
 
 #include <assert.h>
 #include <math.h>
+#include <2geom/coord.h>
 #include <2geom/point.h>
 #include <2geom/transforms.h>
 
@@ -230,6 +231,13 @@ Point constrain_angle(Point const &A, Point const &B, unsigned int n, Point cons
     return A + dir * Rotate(k * 2.0 * M_PI / (double)n) * L2(diff);
 }
 
+std::ostream &operator<<(std::ostream &out, const Geom::Point &p)
+{
+    out << "(" << format_coord_nice(p[X]) << ", "
+               << format_coord_nice(p[Y]) << ")";
+    return out;
+}
+
 } // end namespace Geom
 
 /*
diff --git a/src/2geom/point.h b/src/2geom/point.h
index 3e56ae358..f2659d351 100644
--- a/src/2geom/point.h
+++ b/src/2geom/point.h
@@ -250,17 +250,12 @@ public:
         Dim2 dim;
     };
     //template <typename First = std::less<Coord>, typename Second = std::less<Coord> > LexOrder
-
-    friend inline std::ostream &operator<< (std::ostream &out_file, const Geom::Point &in_pnt);
 };
 
 /** @brief Output operator for points.
  * Prints out the coordinates.
  * @relates Point */
-inline std::ostream &operator<< (std::ostream &out_file, const Geom::Point &in_pnt) {
-    out_file << "X: " << in_pnt[X] << "  Y: " << in_pnt[Y];
-    return out_file;
-}
+std::ostream &operator<<(std::ostream &out, const Geom::Point &p);
 
 template<> struct Point::LexLess<X> {
     typedef std::less<Coord> Primary;
diff --git a/src/2geom/poly.cpp b/src/2geom/polynomial.cpp
index de0229172..a0689f0c5 100644
--- a/src/2geom/poly.cpp
+++ b/src/2geom/polynomial.cpp
@@ -1,4 +1,40 @@
-#include <2geom/poly.h>
+/**
+ * \file
+ * \brief Polynomial in canonical (monomial) basis
+ *//*
+ * Authors:
+ *    MenTaLguY <mental@rydia.net>
+ *    Krzysztof Kosiński <tweenk.pl@gmail.com>
+ * 
+ * Copyright 2007-2015 Authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+#include <algorithm>
+#include <2geom/polynomial.h>
+#include <math.h>
 
 #ifdef HAVE_GSL
 #include <gsl/gsl_poly.h>
@@ -183,6 +219,95 @@ Poly gcd(Poly const &a, Poly const &b, const double /*tol*/) {
 
 
 
+
+std::vector<Coord> solve_quadratic(Coord a, Coord b, Coord c)
+{
+    std::vector<Coord> result;
+
+    if (a == 0) {
+        // linear equation
+        if (b == 0) return result;
+        result.push_back(-c/b);
+        return result;
+    }
+
+    Coord delta = b*b - 4*a*c;
+
+    if (delta == 0) {
+        // one root
+        result.push_back(-0.5 * b / 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));
+    }
+    // no roots otherwise
+
+    std::sort(result.begin(), result.end());
+    return result;
+}
+
+
+std::vector<Coord> solve_cubic(Coord a, Coord b, Coord c, Coord d)
+{
+    // based on:
+    // http://mathworld.wolfram.com/CubicFormula.html
+
+    if (a == 0) {
+        return solve_quadratic(b, c, d);
+    }
+    if (d == 0) {
+        // divide by x
+        std::vector<Coord> result = solve_quadratic(a, b, c);
+        result.push_back(0);
+        std::sort(result.begin(), result.end());
+        return result;
+    }
+
+    std::vector<Coord> result;
+
+    // 1. divide everything by a to bring to canonical form
+    b /= a;
+    c /= a;
+    d /= a;
+
+    // 2. eliminate x^2 term: x^3 + 3Qx - 2R = 0
+    Coord Q = (3*c - b*b) / 9;
+    Coord R = (-27 * d + b * (9*c - 2*b*b)) / 54;
+
+    // 3. compute polynomial discriminant
+    Coord D = Q*Q*Q + R*R;
+    Coord term1 = b/3;
+
+    if (D > 0) {
+        // only one real root
+        Coord S = cbrt(R + sqrt(D));
+        Coord T = cbrt(R - sqrt(D));
+        result.push_back(-b/3 + S + T);
+    } else if (D == 0) {
+        // 3 real roots, 2 of which are equal
+        Coord rroot = cbrt(R);
+        result.reserve(3);
+        result.push_back(-term1 + 2*rroot);
+        result.push_back(-term1 - rroot);
+        result.push_back(-term1 - rroot);
+    } else {
+        // 3 distinct real roots
+        assert(Q < 0);
+        Coord theta = acos(R / sqrt(-Q*Q*Q));
+        Coord rroot = 2 * sqrt(-Q);
+        result.reserve(3);
+        result.push_back(-term1 + rroot * cos(theta / 3));
+        result.push_back(-term1 + rroot * cos((theta + 2*M_PI) / 3));
+        result.push_back(-term1 + rroot * cos((theta + 4*M_PI) / 3));
+    }
+
+    std::sort(result.begin(), result.end());
+    return result;
+}
+
+
 /*Poly divide_out_root(Poly const & p, double x) {
     assert(1);
     }*/
diff --git a/src/2geom/poly.h b/src/2geom/polynomial.h
index 7d93d0a85..aeac2f3fa 100644
--- a/src/2geom/poly.h
+++ b/src/2geom/polynomial.h
@@ -3,9 +3,10 @@
  * \brief Polynomial in canonical (monomial) basis
  *//*
  * Authors:
- *      ? <?@?.?>
+ *    MenTaLguY <mental@rydia.net>
+ *    Krzysztof Kosiński <tweenk.pl@gmail.com>
  * 
- * Copyright ?-?  authors
+ * Copyright 2007-2015 Authors
  *
  * This library is free software; you can redistribute it and/or
  * modify it either under the terms of the GNU Lesser General Public
@@ -29,7 +30,6 @@
  * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
  * OF ANY KIND, either express or implied. See the LGPL or the MPL for
  * the specific language governing rights and limitations.
- *
  */
 
 #ifndef LIB2GEOM_SEEN_POLY_H
@@ -39,6 +39,7 @@
 #include <iostream>
 #include <algorithm>
 #include <complex>
+#include <2geom/coord.h>
 #include <2geom/utils.h>
 
 namespace Geom {
@@ -214,6 +215,18 @@ std::vector<double> solve_reals(const Poly & p);
 #endif
 double polish_root(Poly const & p, double guess, double tol);
 
+
+/** @brief Analytically solve quadratic equation.
+ * The equation is given in the standard form: ax^2 + bx + c = 0.
+ * Only real roots are returned. */
+std::vector<Coord> solve_quadratic(Coord a, Coord b, Coord c);
+
+/** @brief Analytically solve cubic equation.
+ * The equation is given in the standard form: ax^3 + bx^2 + cx + d = 0.
+ * Only real roots are returned. */
+std::vector<Coord> solve_cubic(Coord a, Coord b, Coord c, Coord d);
+
+
 inline std::ostream &operator<< (std::ostream &out_file, const Poly &in_poly) {
     if(in_poly.size() == 0)
         out_file << "0";
diff --git a/src/2geom/quadtree.cpp b/src/2geom/quadtree.cpp
deleted file mode 100644
index 76dc68774..000000000
--- a/src/2geom/quadtree.cpp
+++ /dev/null
@@ -1,286 +0,0 @@
-#include <2geom/quadtree.h>
-
-namespace Geom{
-Quad* QuadTree::search(Rect const &r) {
-    return search(r[0].min(), r[1].min(),
-                  r[0].max(), r[1].max());
-}
-
-void QuadTree::insert(Rect const &r, int shape) {
-    insert(r[0].min(), r[1].min(),
-           r[0].max(), r[1].max(), shape);
-}
-
-
-Quad* QuadTree::search(double x0, double y0, double x1, double y1) {
-    Quad *q = root;
-        
-    double bxx0 = bx1, bxx1 = bx1;
-    double byy0 = by1, byy1 = by1;
-    while(q) {
-        double cx = (bxx0 + bxx1)/2;
-        double cy = (byy0 + byy1)/2;
-        unsigned i = 0;
-        if(x0 >= cx) {
-            i += 1;
-            bxx0 = cx; // zoom in a quad
-        } else if(x1 <= cx) {
-            bxx1 = cx;
-        } else
-            break;
-        if(y0 >= cy) {
-            i += 2;
-            byy0 = cy;
-        } else if(y1 <= cy) {
-            byy1 = cy;
-        } else
-            break;
-        
-        assert(i < 4);
-        Quad *qq = q->children[i];
-        if(qq == 0) break; // last non-null
-        q = qq;
-    }
-    return q;
-}
-
-
-/*
-Comments by Vangelis (use with caution :P )
-
-Insert Rect (x0, y0), (x1, y1) in the QuadTree Q.
-
-===================================================================================
-* QuadTree Q has: Quadtree's Quad root R, QuadTree's bounding box B. 
-
-* Each Quad has a Quad::data where we store the id of the Rect that belong to 
-this Quad. (In reality we'll store a pointer to the shape).
-
-* Each Quad has 4 Quad children: 0, 1, 2, 3. Each child Quad represents one of the following quarters
-of the bounding box B:
-
-+---------------------+
-|          |          |
-|  NW=0    |  NE=1    |
-|          |          |
-|          |          |
-+---------------------+
-|          |          |
-|  SW=2    |  SE=3    |
-|          |          |
-|          |          |
-+---------------------+ 
-
-Each Quad can further be divided in 4 Quads as above and so on. Below there is an example 
- 
-       Root
-      / || \
-    /  /  \  \
-   0  1   2   3
-     /\
-  / | | \
-  0 1 2 3
-
-+---------------------+
-|          | 1-0 | 1-1|
-|    0     |     |    |
-|          |-----|----|
-|          | 1-2 | 1-3|
-|          |     |    |
-+---------------------+
-|          |          |
-|          |          |
-|     2    |     3    |
-|          |          |
-+---------------------+ 
-
-
-
-===================================================================================
-Insert Rect (x0, y0), (x1, y1) in the QuadTree Q. Algorithm:
-1) check if Rect is bigger than QuadTree's bounding box
-2) find in which Quad we should add the Rect:
-
-
-
------------------------------------------------------------------------------------
-How we find in which Quad we should add the Rect R:
-
-Q = Quadtree's Quad root
-B = QuadTree's bounding box B
-WHILE (Q) {
-    IF ( Rect cannot fit in one unique quarter of B ){
-        Q = current Quad ;
-        BREAK;
-    }
-    IF ( Rect can fit in the quarter I ) {
-        IF (Q.children[I] doesn't exist) {
-            create the Quad Q.children[I];
-        }
-        B = bounding box of the Quad Q.children[I] ;
-        Q = Q.children[I] ;
-        CHECK(R, B) ;
-    }
-}
-add Rect R to Q ;
-
-
-*/
-    
-void QuadTree::insert(double x0, double y0, double x1, double y1, int shape) {
-    // loop until a quad would break the box.
-
-    // empty root => empty QuadTree. Create initial bounding box (0,0), (1,1)
-    if(root == NULL) {
-        root = new Quad;
-            
-        bx0 = 0;
-        bx1 = 1;
-        by0 = 0;
-        by1 = 1;
-    }
-    Quad *q = root;
-
-    //A temp bounding box. Same as root's bounting box (ie of the whole QuadTree)
-    double bxx0 = bx0, bxx1 = bx1;
-    double byy0 = by0, byy1 = by1;
-
-    while((bxx0 > x0) ||
-          (bxx1 < x1) ||
-          (byy0 > y0) ||
-          (byy1 < y1)) { 
-        // QuadTree has small size, can't accomodate new rect. Double the size:
-        unsigned i = 0;
-
-        if(bxx0 > x0) {
-            bxx0 = 2*bxx0 - bxx1;
-            i += 1;
-        } else {
-            bxx1 = 2*bxx1 - bxx0;
-        }
-        if(byy0 > y0) {
-            byy0 = 2*byy0 - byy1;
-            i += 2;
-        } else {
-            byy1 = 2*byy1 - byy0;
-        }
-        q = new Quad;
-        //check if root is empty (no rects, no quad children)
-        if( clean_root() ){
-            root = q;
-        }
-        else{
-            q->children[i] = root;
-            root = q;
-        }
-        bx0 = bxx0;
-        bx1 = bxx1;
-        by0 = byy0;
-        by1 = byy1;
-    }
-
-    while(true) {
-        // Find the center of the temp bounding box
-        double cx = (bxx0 + bxx1)/2;
-        double cy = (byy0 + byy1)/2;
-        unsigned i = 0;
-        assert(x0 >= bxx0);
-        assert(x1 <= bxx1);
-        assert(y0 >= byy0);
-        assert(y1 <= byy1);
-
-        if(x0 >= cx) {
-            i += 1;
-            bxx0 = cx; // zoom in a quad
-        } else if(x1 <= cx) {
-            bxx1 = cx;
-        } else{
-            // rect does not fit in one unique quarter (in X axis) of the temp bounding box
-            break;
-        }
-        if(y0 >= cy) {
-            i += 2;
-            byy0 = cy;
-        } else if(y1 <= cy) {
-            byy1 = cy;
-        } else{
-            // rect does not fit in one unique quarter (in Y axis) of the temp bounding box
-            break;
-        }
-
-        // check if rect's bounding box has size 1x1. This means that rect is defined by 2 points
-        // that are in the same place.
-        if( ( fabs(bxx0 - bxx1) < 1.0 ) && ( fabs(byy0 - byy1) < 1.0 )){
-            break;
-        }
-
-        /*
-        1 rect does fit in one unique quarter of the temp bounding box. And we have found which.
-        2 temp bounding box = bounding box of this quarter. 
-        3 "Go in" this quarter (create if doesn't exist)
-        */
-        assert(i < 4);
-        Quad *qq = q->children[i];
-        if(qq == NULL) {
-            qq = new Quad;
-            q->children[i] = qq;
-        }
-        q = qq;
-    }
-    q->data.push_back(shape);
-}
-
-
-void QuadTree::erase(Quad *q, int shape) {
-    for(Quad::iterator i = q->data.begin();  i != q->data.end(); ++i) {
-        if(*i == shape) {
-            q->data.erase(i);
-            if(q->data.empty()) {
-
-            }
-        }
-    }
-    return;
-}
-
-/*
-Returns:
-false:  if root isn't empty
-true:   if root is empty it cleans root
-*/
-bool QuadTree::clean_root() { 
-    assert(root);
-
-    // false if root *has* rects assigned to it.
-    bool all_clean = root->data.empty(); 
-
-    // if root has children we get false
-    for(unsigned i = 0; i < 4; i++)
-    {
-        if(root->children[i])
-        {
-            all_clean = false;
-        }
-    }
-
-    if(all_clean)
-    {
-        delete root;
-        root=0;
-        return true;
-    }
-    return false;
-}
-
-};
-
-/*
-  Local Variables:
-  mode:c++
-  c-file-style:"stroustrup"
-  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
-  indent-tabs-mode:nil
-  fill-column:99
-  End:
-*/
-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
diff --git a/src/2geom/quadtree.h b/src/2geom/quadtree.h
deleted file mode 100644
index 62f2b8321..000000000
--- a/src/2geom/quadtree.h
+++ /dev/null
@@ -1,105 +0,0 @@
-/**
- * \file
- * \brief Quad tree data structure
- *//*
- * Authors:
- *      ? <?@?.?>
- * 
- * Copyright ?-?  authors
- *
- * This library is free software; you can redistribute it and/or
- * modify it either under the terms of the GNU Lesser General Public
- * License version 2.1 as published by the Free Software Foundation
- * (the "LGPL") or, at your option, under the terms of the Mozilla
- * Public License Version 1.1 (the "MPL"). If you do not alter this
- * notice, a recipient may use your version of this file under either
- * the MPL or the LGPL.
- *
- * You should have received a copy of the LGPL along with this library
- * in the file COPYING-LGPL-2.1; if not, write to the Free Software
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
- * You should have received a copy of the MPL along with this library
- * in the file COPYING-MPL-1.1
- *
- * The contents of this file are subject to the Mozilla Public License
- * Version 1.1 (the "License"); you may not use this file except in
- * compliance with the License. You may obtain a copy of the License at
- * http://www.mozilla.org/MPL/
- *
- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
- * OF ANY KIND, either express or implied. See the LGPL or the MPL for
- * the specific language governing rights and limitations.
- *
- */
-
-#ifndef LIB2GEOM_SEEN_QUADTREE_H
-#define LIB2GEOM_SEEN_QUADTREE_H
-
-#include <vector>
-#include <cassert>
-
-#include <2geom/d2.h>
-
-namespace Geom{
-
-class Quad{
-public:
-    Quad* children[4];
-    std::vector<int> data;
-    Quad() {
-        for(int i = 0; i < 4; i++)
-            children[i] = 0;
-    }
-    typedef std::vector<int>::iterator iterator;
-    Rect bounds(unsigned i, double x, double y, double d) {
-        double dd = d/2;
-        switch(i % 4) {
-            case 0:
-                return Rect(Interval(x, x+dd), Interval(y, y+dd));
-            case 1:
-                return Rect(Interval(x+dd, x+d), Interval(y, y+dd));
-            case 2:
-                return Rect(Interval(x, x+dd), Interval(y+dd, y+d));
-            case 3:
-                return Rect(Interval(x+dd, x+d), Interval(y+dd, y+d));
-            default: 
-                /* just to suppress warning message
-                 * this case should be never reached */
-                assert(false);
-        }        
-        assert(false);
-        return Rect();
-    }
-};
-
-class QuadTree{
-public:
-    Quad* root;
-    double scale;
-    double bx0, bx1;
-    double by0, by1;
-
-    QuadTree() : root(0), scale(1) {}
-
-    Quad* search(double x0, double y0, double x1, double y1);
-    void insert(double x0, double y0, double x1, double y1, int shape);
-    Quad* search(Geom::Rect const &r);
-    void insert(Geom::Rect const &r, int shape);
-    void erase(Quad *q, int shape);
-private:
-    bool clean_root();
-};
-
-};
-
-#endif
-/*
-  Local Variables:
-  mode:c++
-  c-file-style:"stroustrup"
-  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
-  indent-tabs-mode:nil
-  fill-column:99
-  End:
-*/
-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
diff --git a/src/2geom/sbasis-curve.h b/src/2geom/sbasis-curve.h
index a7b74dcee..17ee8f4c9 100644
--- a/src/2geom/sbasis-curve.h
+++ b/src/2geom/sbasis-curve.h
@@ -88,6 +88,7 @@ public:
     virtual Point initialPoint() const    { return inner.at0(); }
     virtual Point finalPoint() const      { return inner.at1(); }
     virtual bool isDegenerate() const     { return inner.isConstant(0); }
+    virtual bool isLineSegment() const    { return inner[X].size() == 1; }
     virtual Point pointAt(Coord t) const  { return inner.valueAt(t); }
     virtual std::vector<Point> pointAndDerivatives(Coord t, unsigned n) const {
         return inner.valueAndDerivatives(t, n);
diff --git a/src/2geom/sbasis-poly.h b/src/2geom/sbasis-poly.h
index 81eeca72c..d18bc369e 100644
--- a/src/2geom/sbasis-poly.h
+++ b/src/2geom/sbasis-poly.h
@@ -33,7 +33,7 @@
 #ifndef LIB2GEOM_SEEN_SBASIS_POLY_H
 #define LIB2GEOM_SEEN_SBASIS_POLY_H
 
-#include <2geom/poly.h>
+#include <2geom/polynomial.h>
 #include <2geom/sbasis.h>
 
 namespace Geom{
diff --git a/src/2geom/sbasis-to-bezier.cpp b/src/2geom/sbasis-to-bezier.cpp
index dfd07d84c..09fbb03ef 100644
--- a/src/2geom/sbasis-to-bezier.cpp
+++ b/src/2geom/sbasis-to-bezier.cpp
@@ -239,8 +239,6 @@ void sbasis_to_cubic_bezier (std::vector<Point> & bz, D2<SBasis> const& sb)
 
     midx = 8*midx - 4*bz[0][X] - 4*bz[3][X];  // re-define relative to center
     midy = 8*midy - 4*bz[0][Y] - 4*bz[3][Y];
-    if ((std::abs(midx) < EPSILON) && (std::abs(midy) < EPSILON))
-        return;
 
     if ((std::abs(xprime[0]) < EPSILON) && (std::abs(yprime[0]) < EPSILON)
     && ((std::abs(xprime[1]) > EPSILON) || (std::abs(yprime[1]) > EPSILON)))  { // degenerate handle at 0 : use distance of closest approach
@@ -264,11 +262,11 @@ void sbasis_to_cubic_bezier (std::vector<Point> & bz, D2<SBasis> const& sb)
         double test2 = (bz[2][Y] - bz[0][Y])*(bz[3][X] - bz[0][X]) - (bz[2][X] - bz[0][X])*(bz[3][Y] - bz[0][Y]);
         if (test1*test2 < 0) // reject anti-symmetric case, LP Bug 1428267 & Bug 1428683
             return;
-        denom = xprime[1]*yprime[0] - yprime[1]*xprime[0];
+        denom = 3.0*(xprime[1]*yprime[0] - yprime[1]*xprime[0]);
         for (int i = 0; i < 2; ++i) {
             numer = xprime[1 - i]*midy - yprime[1 - i]*midx;
-            delx[i] = xprime[i]*numer/denom/3;
-            dely[i] = yprime[i]*numer/denom/3;
+            delx[i] = xprime[i]*numer/denom;
+            dely[i] = yprime[i]*numer/denom;
         }
     } else if ((xprime[0]*xprime[1] < 0) || (yprime[0]*yprime[1] < 0)) { // symmetric case : use distance of closest approach
         numer = midx*xprime[0] + midy*yprime[0];
@@ -515,10 +513,10 @@ path_from_sbasis(D2<SBasis> const &B, double tol, bool only_cubicbeziers) {
  TODO: some of this logic should be lifted into svg-path
 */
 PathVector
-path_from_piecewise(Piecewise<D2<SBasis> > const &B, double tol, bool only_cubicbeziers) {
-    PathBuilder pb;
+path_from_piecewise(Geom::Piecewise<Geom::D2<Geom::SBasis> > const &B, double tol, bool only_cubicbeziers) {
+    Geom::PathBuilder pb;
     if(B.size() == 0) return pb.peek();
-    Point start = B[0].at0();
+    Geom::Point start = B[0].at0();
     pb.moveTo(start);
     for(unsigned i = 0; ; i++) {
         if ( (i+1 == B.size()) 
diff --git a/src/2geom/svg-elliptical-arc.h b/src/2geom/svg-elliptical-arc.h
deleted file mode 100644
index 90e16a999..000000000
--- a/src/2geom/svg-elliptical-arc.h
+++ /dev/null
@@ -1,280 +0,0 @@
-/**
- * \file
- * \brief SVG 1.1-compliant elliptical arc curve
- *//*
- * Authors:
- *    MenTaLguY <mental@rydia.net>
- *    Marco Cecchetti <mrcekets at gmail.com>
- *    Krzysztof Kosiński <tweenk.pl@gmail.com>
- * 
- * Copyright 2007-2009 Authors
- *
- * This library is free software; you can redistribute it and/or
- * modify it either under the terms of the GNU Lesser General Public
- * License version 2.1 as published by the Free Software Foundation
- * (the "LGPL") or, at your option, under the terms of the Mozilla
- * Public License Version 1.1 (the "MPL"). If you do not alter this
- * notice, a recipient may use your version of this file under either
- * the MPL or the LGPL.
- *
- * You should have received a copy of the LGPL along with this library
- * in the file COPYING-LGPL-2.1; if not, write to the Free Software
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
- * You should have received a copy of the MPL along with this library
- * in the file COPYING-MPL-1.1
- *
- * The contents of this file are subject to the Mozilla Public License
- * Version 1.1 (the "License"); you may not use this file except in
- * compliance with the License. You may obtain a copy of the License at
- * http://www.mozilla.org/MPL/
- *
- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
- * OF ANY KIND, either express or implied. See the LGPL or the MPL for
- * the specific language governing rights and limitations.
- */
-
-
-#ifndef LIB2GEOM_SEEN_SVG_ELLIPTICAL_ARC_H
-#define LIB2GEOM_SEEN_SVG_ELLIPTICAL_ARC_H
-
-#include <2geom/curve.h>
-#include <2geom/angle.h>
-#include <2geom/utils.h>
-#include <2geom/bezier-curve.h>
-#include <2geom/elliptical-arc.h>
-#include <2geom/sbasis-curve.h>  // for non-native methods
-#include <2geom/numeric/vector.h>
-#include <2geom/numeric/fitting-tool.h>
-#include <2geom/numeric/fitting-model.h>
-#include <algorithm>
-
-namespace Geom {
-
-class SVGEllipticalArc : public EllipticalArc {
-public:
-    SVGEllipticalArc()
-        : EllipticalArc()
-    {}
-    SVGEllipticalArc( Point ip, double rx, double ry,
-                      double rot_angle, bool large_arc, bool sweep,
-                      Point fp
-                    )
-        : EllipticalArc()
-    {
-        _initial_point = ip;
-        _final_point = fp;
-        _rays[X] = rx; _rays[Y] = ry;
-        _rot_angle = rot_angle;
-        _large_arc = large_arc;
-        _sweep = sweep;
-        _updateCenterAndAngles(true);
-    }
-
-    virtual Curve *duplicate() const {
-        return new SVGEllipticalArc(*this);
-    }
-    virtual Coord valueAt(Coord t, Dim2 d) const {
-        if (isChord()) return chord().valueAt(t, d);
-        return EllipticalArc::valueAt(t, d);
-    }
-    virtual Point pointAt(Coord t) const {
-        if (isChord()) return chord().pointAt(t);
-        return EllipticalArc::pointAt(t);
-    }
-    virtual std::vector<Point> pointAndDerivatives(Coord t, unsigned int n) const {
-        if (isChord()) return chord().pointAndDerivatives(t, n);
-        return EllipticalArc::pointAndDerivatives(t, n);
-    }
-    virtual Rect boundsExact() const {
-        if (isChord()) return chord().boundsExact();
-        return EllipticalArc::boundsExact();
-    }
-    virtual OptRect boundsLocal(OptInterval const &i, unsigned int deg) const {
-        if (isChord()) return chord().boundsLocal(i, deg);
-        return EllipticalArc::boundsLocal(i, deg);
-    }
-
-    virtual Curve *derivative() const {
-        if (isChord()) return chord().derivative();
-        return EllipticalArc::derivative();
-    }
-
-    virtual std::vector<Coord> roots(Coord v, Dim2 d) const {
-        if (isChord()) return chord().roots(v, d);
-        return EllipticalArc::roots(v, d);
-    }
-#ifdef HAVE_GSL
-    virtual std::vector<Coord> allNearestTimes( Point const& p, double from = 0, double to = 1 ) const {
-        if (isChord()) {
-            std::vector<Coord> result;
-            result.push_back(chord().nearestTime(p, from, to));
-            return result;
-        }
-        return EllipticalArc::allNearestTimes(p, from, to);
-    }
-#endif
-    virtual D2<SBasis> toSBasis() const {
-        if (isChord()) return chord().toSBasis();
-        return EllipticalArc::toSBasis();
-    }
-    virtual bool isSVGCompliant() const { return true; }
-    // TODO move SVG-specific behavior here.
-//protected:
-    //virtual void _updateCenterAndAngles();
-}; // end class SVGEllipticalArc
-
-/*
- * useful for testing and debugging
- */
-template< class charT >
-inline
-std::basic_ostream<charT> &
-operator<< (std::basic_ostream<charT> & os, const SVGEllipticalArc & ea)
-{
-    os << "{ cx: " << ea.center(X) << ", cy: " <<  ea.center(Y)
-       << ", rx: " << ea.ray(X) << ", ry: " << ea.ray(Y)
-       << ", rot angle: " << decimal_round(rad_to_deg(ea.rotationAngle()),2)
-       << ", start angle: " << decimal_round(rad_to_deg(ea.initialAngle()),2)
-       << ", end angle: " << decimal_round(rad_to_deg(ea.finalAngle()),2)
-       << " }";
-
-    return os;
-}
-
-
-
-
-// forward declation
-namespace detail
-{
-    struct ellipse_equation;
-}
-
-// TODO this needs to be rewritten and moved to EllipticalArc header
-/*
- * make_elliptical_arc
- *
- * convert a parametric polynomial curve given in symmetric power basis form
- * into an SVGEllipticalArc type; in order to be successfull the input curve
- * has to look like an actual elliptical arc even if a certain tolerance
- * is allowed through an ad-hoc parameter.
- * The conversion is performed through an interpolation on a certain amount of
- * sample points computed on the input curve;
- * the interpolation computes the coefficients of the general implicit equation
- * of an ellipse (A*X^2 + B*XY + C*Y^2 + D*X + E*Y + F = 0), then from the
- * implicit equation we compute the parametric form.
- *
- */
-class make_elliptical_arc
-{
-  public:
-    typedef D2<SBasis> curve_type;
-
-    /*
-     * constructor
-     *
-     * it doesn't execute the conversion but set the input and output parameters
-     *
-     * _ea:         the output SVGEllipticalArc that will be generated;
-     * _curve:      the input curve to be converted;
-     * _total_samples: the amount of sample points to be taken
-     *                 on the input curve for performing the conversion
-     * _tolerance:     how much likelihood is required between the input curve
-     *                 and the generated elliptical arc; the smaller it is the
-     *                 the tolerance the higher it is the likelihood.
-     */
-    make_elliptical_arc( EllipticalArc& _ea,
-                         curve_type const& _curve,
-                         unsigned int _total_samples,
-                         double _tolerance );
-
-  private:
-    bool bound_exceeded( unsigned int k, detail::ellipse_equation const & ee,
-                         double e1x, double e1y, double e2 );
-
-    bool check_bound(double A, double B, double C, double D, double E, double F);
-
-    void fit();
-
-    bool make_elliptiarc();
-
-    void print_bound_error(unsigned int k)
-    {
-        std::cerr
-            << "tolerance error" << std::endl
-            << "at point: " << k << std::endl
-            << "error value: "<< dist_err << std::endl
-            << "bound: " << dist_bound << std::endl
-            << "angle error: " << angle_err
-            << " (" << angle_tol << ")" << std::endl;
-    }
-
-  public:
-    /*
-     * perform the actual conversion
-     * return true if the conversion is successfull, false on the contrary
-     */
-    bool operator()()
-    {
-        // initialize the reference
-        const NL::Vector & coeff = fitter.result();
-        fit();
-        if ( !check_bound(1, coeff[0], coeff[1], coeff[2], coeff[3], coeff[4]) )
-            return false;
-        if ( !(make_elliptiarc()) ) return false;
-        return true;
-    }
-
-    /*
-     * you can set a boolean parameter to tell the conversion routine
-     * if the output elliptical arc has to be svg compliant or not;
-     * the default value is true
-     */
-    bool svg_compliant_flag() const
-    {
-        return svg_compliant;
-    }
-
-    void svg_compliant_flag(bool _svg_compliant)
-    {
-        svg_compliant = _svg_compliant;
-    }
-
-  private:
-      EllipticalArc& ea;                 // output elliptical arc
-      const curve_type & curve;             // input curve
-      Piecewise<D2<SBasis> > dcurve;        // derivative of the input curve
-      NL::LFMEllipse model;                 // model used for fitting
-      // perform the actual fitting task
-      NL::least_squeares_fitter<NL::LFMEllipse> fitter;
-      // tolerance: the user-defined tolerance parameter;
-      // tol_at_extr: the tolerance at end-points automatically computed
-      // on the value of "tolerance", and usually more strict;
-      // tol_at_center: tolerance at the center of the ellipse
-      // angle_tol: tolerance for the angle btw the input curve tangent
-      // versor and the ellipse normal versor at the sample points
-      double tolerance, tol_at_extr, tol_at_center, angle_tol;
-      Point initial_point, final_point;     // initial and final end-points
-      unsigned int N;                       // total samples
-      unsigned int last; // N-1
-      double partitions; // N-1
-      std::vector<Point> p;                 // sample points
-      double dist_err, dist_bound, angle_err;
-      bool svg_compliant;
-};
-
-} // end namespace Geom
-
-#endif // LIB2GEOM_SEEN_SVG_ELLIPTICAL_ARC_H
-
-/*
-  Local Variables:
-  mode:c++
-  c-file-style:"stroustrup"
-  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
-  indent-tabs-mode:nil
-  fill-column:99
-  End:
-*/
-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
-
diff --git a/src/2geom/svg-path-parser.cpp b/src/2geom/svg-path-parser.cpp
index 6ee55171e..b6e6da869 100644
--- a/src/2geom/svg-path-parser.cpp
+++ b/src/2geom/svg-path-parser.cpp
@@ -1209,7 +1209,7 @@ void SVGPathParser::_arcTo(Coord rx, Coord ry, Coord angle,
         return; // ignore invalid (ambiguous) arc segments where start and end point are the same (per SVG spec)
     }
 
-    _pushCurve(new SVGEllipticalArc(_current, rx, ry, angle, large_arc, sweep, p));
+    _pushCurve(new EllipticalArc(_current, rx, ry, angle, large_arc, sweep, p));
     _quad_tangent = _cubic_tangent = _current = p;
 }
 
diff --git a/src/2geom/svg-path-parser.h b/src/2geom/svg-path-parser.h
index 9f304b9ed..e25316cf8 100644
--- a/src/2geom/svg-path-parser.h
+++ b/src/2geom/svg-path-parser.h
@@ -37,6 +37,7 @@
 #include <iterator>
 #include <stdexcept>
 #include <vector>
+#include <cstdio>
 #include <2geom/exception.h>
 #include <2geom/point.h>
 #include <2geom/path-sink.h>
diff --git a/src/2geom/sweep.cpp b/src/2geom/sweep-bounds.cpp
index 6910d5d02..48f168b97 100644
--- a/src/2geom/sweep.cpp
+++ b/src/2geom/sweep-bounds.cpp
@@ -1,9 +1,27 @@
-#include <2geom/sweep.h>
+#include <2geom/sweep-bounds.h>
 
 #include <algorithm>
 
 namespace Geom {
 
+struct Event {
+    double x;
+    unsigned ix;
+    bool closing;
+    Event(double pos, unsigned i, bool c) : x(pos), ix(i), closing(c) {}
+// Lexicographic ordering by x then closing
+    bool operator<(Event const &other) const {
+        if(x < other.x) return true;
+        if(x > other.x) return false;
+        return closing < other.closing;
+    }
+    bool operator==(Event const &other) const {
+        return other.x == x && other.ix == ix && other.closing == closing;
+    }
+};
+
+std::vector<std::vector<unsigned> > fake_cull(unsigned a, unsigned b);
+
 /**
  * \brief Make a list of pairs of self intersections in a list of Rects.
  * 
diff --git a/src/2geom/sweep.h b/src/2geom/sweep-bounds.h
index 9c012958b..e0ebf2975 100644
--- a/src/2geom/sweep.h
+++ b/src/2geom/sweep-bounds.h
@@ -40,29 +40,12 @@
 
 namespace Geom {
 
-struct Event {
-    double x;
-    unsigned ix;
-    bool closing;
-    Event(double pos, unsigned i, bool c) : x(pos), ix(i), closing(c) {}
-// Lexicographic ordering by x then closing
-    bool operator<(Event const &other) const {
-        if(x < other.x) return true;
-        if(x > other.x) return false;
-        return closing < other.closing;
-    }
-    bool operator==(Event const &other) const {
-        return other.x == x && other.ix == ix && other.closing == closing;
-    }
-};
 std::vector<std::vector<unsigned> >
 sweep_bounds(std::vector<Rect>, Dim2 dim = X);
 
 std::vector<std::vector<unsigned> >
 sweep_bounds(std::vector<Rect>, std::vector<Rect>, Dim2 dim = X);
 
-std::vector<std::vector<unsigned> > fake_cull(unsigned a, unsigned b);
-
 }
 
 #endif
diff --git a/src/2geom/sweeper.h b/src/2geom/sweeper.h
new file mode 100644
index 000000000..1cbce52b0
--- /dev/null
+++ b/src/2geom/sweeper.h
@@ -0,0 +1,220 @@
+/** @file
+ * @brief Class for implementing sweepline algorithms
+ *//*
+  * Authors:
+  *   Krzysztof Kosiński <tweenk.pl@gmail.com>
+  *
+  * Copyright 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
+  * License version 2.1 as published by the Free Software Foundation
+  * (the "LGPL") or, at your option, under the terms of the Mozilla
+  * Public License Version 1.1 (the "MPL"). If you do not alter this
+  * notice, a recipient may use your version of this file under either
+  * the MPL or the LGPL.
+  *
+  * You should have received a copy of the LGPL along with this library
+  * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+  * You should have received a copy of the MPL along with this library
+  * in the file COPYING-MPL-1.1
+  *
+  * The contents of this file are subject to the Mozilla Public License
+  * Version 1.1 (the "License"); you may not use this file except in
+  * compliance with the License. You may obtain a copy of the License at
+  * http://www.mozilla.org/MPL/
+  *
+  * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+  * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+  * the specific language governing rights and limitations.
+  */
+
+#ifndef LIB2GEOM_SEEN_SWEEPER_H
+#define LIB2GEOM_SEEN_SWEEPER_H
+
+#include <2geom/coord.h>
+#include <algorithm>
+#include <vector>
+#include <boost/intrusive/list.hpp>
+#include <boost/range/algorithm/heap_algorithm.hpp>
+
+namespace Geom {
+
+struct IntervalSweepTraits {
+    typedef Interval Bound;
+    typedef std::less<Coord> Compare;
+    inline static Coord entry_value(Bound const &b) { return b.min(); }
+    inline static Coord exit_value(Bound const &b) { return b.max(); }
+};
+
+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(); }
+};
+
+/** @brief Generic sweepline algorithm.
+ *
+ * This class encapsulates an algorithm that sorts the objects according
+ * to their bounds, then moves an imaginary line (sweepline) over those
+ * bounds from left to right. Objects are added to the active list when
+ * the line starts intersecting their bounds, and removed when it completely
+ * passes over them.
+ *
+ * To use this, create a derived class and reimplement the _enter()
+ * and/or _leave() virtual functions, insert all the objects,
+ * and finally call process(). 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.
+ */
+template <typename Item, typename SweepTraits = IntervalSweepTraits>
+class Sweeper {
+public:
+    typedef typename SweepTraits::Bound Bound;
+
+    Sweeper() {}
+
+    void insert(Bound const &bound, Item const &item) {
+        assert(!(typename SweepTraits::Compare()(
+            SweepTraits::exit_value(bound),
+            SweepTraits::entry_value(bound))));
+        _items.push_back(Record(bound, item));
+    }
+
+    template <typename Iter, typename BoundFunc>
+    void insert(Iter first, Iter last, BoundFunc f = BoundFunc()) {
+        for (; first != last; ++first) {
+            Bound b = f(*first);
+            assert(!(typename SweepTraits::Compare()(
+                SweepTraits::exit_value(b),
+                SweepTraits::entry_value(b))));
+            _items.push_back(Record(b, *first));
+        }
+    }
+
+    /** @brief Process entry and exit events.
+     * This will iterate over all inserted items, calling the virtual protected
+     * functions _enter() and _leave() according to the order of the boundaries
+     * of each item. */
+    void process() {
+        if (_items.empty()) return;
+
+        typename SweepTraits::Compare cmp;
+
+        for (RecordIter i = _items.begin(); i != _items.end(); ++i) {
+            _entry_events.push_back(i);
+            _exit_events.push_back(i);
+        }
+        boost::make_heap(_entry_events, _entry_heap_compare);
+        boost::make_heap(_exit_events, _exit_heap_compare);
+        boost::pop_heap(_entry_events, _entry_heap_compare);
+        boost::pop_heap(_exit_events, _exit_heap_compare);
+
+        RecordIter next_entry = _entry_events.back();
+        RecordIter next_exit = _exit_events.back();
+        _entry_events.pop_back();
+        _exit_events.pop_back();
+
+        while (next_entry != _items.end() || next_exit != _items.end()) {
+            assert(next_exit != _items.end());
+
+            if (next_entry == _items.end() ||
+                cmp(SweepTraits::exit_value(next_exit->bound),
+                    SweepTraits::entry_value(next_entry->bound)))
+            {
+                // exit event - remove record from active list
+                _leave(*next_exit);
+                _active_items.erase(_active_items.iterator_to(*next_exit));
+                if (!_exit_events.empty()) {
+                    boost::pop_heap(_exit_events, _exit_heap_compare);
+                    next_exit = _exit_events.back();
+                    _exit_events.pop_back();
+                } else {
+                    next_exit = _items.end();
+                    // we should end the loop after this happens
+                }
+            } else {
+                // entry event - add record to active list
+                _enter(*next_entry);
+                _active_items.push_back(*next_entry);
+                if (!_entry_events.empty()) {
+                    boost::pop_heap(_entry_events, _entry_heap_compare);
+                    next_entry = _entry_events.back();
+                    _entry_events.pop_back();
+                } else {
+                    next_entry = _items.end();
+                }
+            }
+        }
+
+        assert(_active_items.empty());
+    }
+
+protected:
+    /// The item and its sweepline boundary.
+    struct Record {
+        boost::intrusive::list_member_hook<> _hook;
+        Bound bound;
+        Item item;
+
+        Record(Bound const &b, Item const &i)
+            : bound(b), item(i)
+        {}
+    };
+    typedef typename std::vector<Record>::iterator RecordIter;
+
+    typedef boost::intrusive::list
+    < Record
+    , boost::intrusive::member_hook
+        < Record
+        , boost::intrusive::list_member_hook<>
+        , &Record::_hook
+        >
+    > RecordList;
+
+    /** @brief Enter an item record.
+     * Override this to process an item as it is about to enter the active list.
+     * When called, the passed record will not be part of the active list. */
+    virtual void _enter(Record const &) {}
+    /** @brief Leave an item record.
+     * Override this to process an item as it is about to leave the active list.
+     * When called, the passed record will be part of the active list. */
+    virtual void _leave(Record const &) {}
+
+    /// The list of all item records undergoing sweeping.
+    std::vector<Record> _items;
+    /// The list of active item records.
+    RecordList _active_items;
+
+private:
+    inline static bool _entry_heap_compare(RecordIter a, RecordIter b) {
+        typename SweepTraits::Compare cmp;
+        return cmp(SweepTraits::entry_value(b->bound), SweepTraits::entry_value(a->bound));
+    }
+    inline static bool _exit_heap_compare(RecordIter a, RecordIter b) {
+        typename SweepTraits::Compare cmp;
+        return cmp(SweepTraits::exit_value(b->bound), SweepTraits::exit_value(a->bound));
+    }
+
+    std::vector<RecordIter> _entry_events;
+    std::vector<RecordIter> _exit_events;
+};
+
+} // namespace Geom
+
+#endif // !LIB2GEOM_SEEN_SWEEPER_H
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
diff --git a/src/2geom/transforms.cpp b/src/2geom/transforms.cpp
index 091079d5a..41d395297 100644
--- a/src/2geom/transforms.cpp
+++ b/src/2geom/transforms.cpp
@@ -139,6 +139,24 @@ Affine &Affine::operator*=(Zoom const &z) {
     return *this;
 }
 
+Affine Rotate::around(Point const &p, Coord angle)
+{
+    Affine result = Translate(-p) * Rotate(angle) * Translate(p);
+    return result;
+}
+
+Affine reflection(Point const & vector, Point const & origin)
+{
+    Geom::Point vn = unit_vector(vector);
+    Coord cx2 = vn[X] * vn[X];
+    Coord cy2 = vn[Y] * vn[Y];
+    Coord c2xy = 2 * vn[X] * vn[Y];
+    Affine mirror ( cx2 - cy2, c2xy,
+                    c2xy, cy2 - cx2,
+                    0, 0 );
+    return Translate(-origin) * mirror * Translate(origin);
+}
+
 // this checks whether the requirements of TransformConcept are satisfied for all transforms.
 // if you add a new transform type, include it here!
 void check_transforms()
@@ -173,14 +191,6 @@ void check_transforms()
     m = z * t; m = z * s; m = z * r; m = z * h; m = z * v; m = z * z;
 }
 
-Affine reflection(Point const & vector, Point const & origin) {
-    Geom::Point vec_norm = unit_vector(vector);
-    Affine mirror ( vec_norm[X]*vec_norm[X] - vec_norm[Y]*vec_norm[Y], 2 * vec_norm[X] * vec_norm[Y] ,
-                    2 * vec_norm[X] * vec_norm[Y], vec_norm[Y]*vec_norm[Y] - vec_norm[X]*vec_norm[X] ,
-                    0 ,0 );
-    return Translate(-origin) * mirror * Translate(origin);
-}
-
 } // namespace Geom
 
 /*
diff --git a/src/2geom/transforms.h b/src/2geom/transforms.h
index 2e805786d..7c03c5226 100644
--- a/src/2geom/transforms.h
+++ b/src/2geom/transforms.h
@@ -217,6 +217,7 @@ public:
         Coord rad = (deg / 180.0) * M_PI;
         return Rotate(rad);
     }
+    static Affine around(Point const &p, Coord angle);
 
     friend class Point;
 };
diff --git a/src/extension/internal/wmf-print.cpp b/src/extension/internal/wmf-print.cpp
index 271dec702..431053085 100644
--- a/src/extension/internal/wmf-print.cpp
+++ b/src/extension/internal/wmf-print.cpp
@@ -30,7 +30,7 @@
 
 
 #include <2geom/sbasis-to-bezier.h>
-#include <2geom/svg-elliptical-arc.h>
+#include <2geom/elliptical-arc.h>
 
 #include <2geom/path.h>
 #include <2geom/pathvector.h>
diff --git a/src/helper/geom-pathstroke.cpp b/src/helper/geom-pathstroke.cpp
index a152af62c..888e55180 100644
--- a/src/helper/geom-pathstroke.cpp
+++ b/src/helper/geom-pathstroke.cpp
@@ -10,83 +10,14 @@
 #include <2geom/path-sink.h>
 #include <2geom/point.h>
 #include <2geom/bezier-curve.h>
-#include <2geom/svg-elliptical-arc.h>
+#include <2geom/elliptical-arc.h>
 #include <2geom/sbasis-to-bezier.h> // cubicbezierpath_from_sbasis
 #include <2geom/path-intersection.h>
+#include <2geom/circle.h>
 
 #include "helper/geom-pathstroke.h"
 
 namespace Geom {
-// 2geom/circle-circle.cpp, no header
-int circle_circle_intersection(Point X0, double r0, Point X1, double r1, Point &p0, Point &p1);
-
-/**
- * Determine the intersection points between a circle C0 and a line defined
- * by two points, X0 and X1.
- *
- * Which intersection point is assigned to p0 or p1 is unspecified, and callers
- * should not depend on any particular intersection always being assigned to p0.
- *
- * Returns:
- *   If the line and circle do not cross, 0 is returned.
- *   If solution(s) exist, 2 is returned, and the results are written to p0 and p1.
- */
-static int circle_line_intersection(Circle C0, Point X0, Point X1, Point &p0, Point &p1)
-{
-    /* equation of a circle: (x - h)^2 + (y - k)^2 = r^2 */
-    Coord r = C0.radius();
-    Coord h = C0.center()[X];
-    Coord k = C0.center()[Y];
-
-    Coord x0, y0;
-    Coord x1, y1;
-
-    if (are_near(X1[X], X0[X])) {
-        /* slope is undefined (vertical line) */
-        Coord c = X0[X];
-        Coord det = r*r - (c-h)*(c-h);
-
-        /* no intersection */
-        if (det < 0)
-            return 0;
-
-        /* solve for y */
-        y0 = k + std::sqrt(det);
-        y1 = k - std::sqrt(det);
-
-        // x == c (always)
-        x0 = c;
-        x1 = c;
-    } else {
-        /* equation of a line: y = mx + b */
-        Coord m = (X1[Y] - X0[Y]) / (X1[X] - X0[X]);
-        Coord b = X0[Y] - m*X0[X];
-
-        /* obtain quadratic for x: */
-        Coord A = m*m + 1;
-        Coord B = 2*h - 2*b*m + 2*k*m;
-        Coord C = b*b + h*h + k*k - r*r - 2*b*k;
-
-        Coord det = B*B - 4*A*C;
-        
-        /* no intersection, circle and line do not cross */
-        if (det < 0)
-            return 0;
-
-        /* solve quadratic */
-        x0 = (B + std::sqrt(det)) / (2*A);
-        x1 = (B - std::sqrt(det)) / (2*A);
-
-        /* substitute the calculated x times to determine the y values */
-        y0 = m*x0 + b;
-        y1 = m*x1 + b;
-    }
-
-    p0 = Point(x0, y0);
-    p1 = Point(x1, y1);
-
-    return 2;
-}
 
 static Point intersection_point(Point origin_a, Point vector_a, Point origin_b, Point vector_b)
 {
@@ -160,7 +91,7 @@ void bevel_join(join_data jd)
 
 void round_join(join_data jd)
 {
-    jd.res.appendNew<Geom::SVGEllipticalArc>(jd.width, jd.width, 0, false, jd.width <= 0, jd.outgoing.initialPoint());
+    jd.res.appendNew<Geom::EllipticalArc>(jd.width, jd.width, 0, false, jd.width <= 0, jd.outgoing.initialPoint());
     jd.res.append(jd.outgoing);
 }
 
@@ -226,18 +157,20 @@ void miter_join_internal(join_data jd, bool clip)
 void miter_join(join_data jd) { miter_join_internal(jd, false); }
 void miter_clip_join(join_data jd) { miter_join_internal(jd, true); }
 
-Geom::Point pick_solution(Geom::Point points[2], Geom::Point tang2, Geom::Point endPt)
+Geom::Point pick_solution(std::vector<Geom::ShapeIntersection> points, Geom::Point tang2, Geom::Point endPt)
 {
+    assert(points.size() == 2);
     Geom::Point sol;
-    if ( dot(tang2,points[0]-endPt) > 0 ) {
+    if ( dot(tang2, points[0].point() - endPt) > 0 ) {
         // points[0] is bad, choose points[1]
         sol = points[1];
-    } else if ( dot(tang2,points[1]-endPt) > 0 ) { // points[0] could be good, now check points[1]
+    } else if ( dot(tang2, points[1].point() - endPt) > 0 ) { // points[0] could be good, now check points[1]
         // points[1] is bad, choose points[0]
         sol = points[0];
     } else {
         // both points are good, choose nearest
-        sol = ( distanceSq(endPt, points[0]) < distanceSq(endPt, points[1]) ) ? points[0] : points[1];
+        sol = ( distanceSq(endPt, points[0].point()) < distanceSq(endPt, points[1].point()) )
+            ? points[0].point() : points[1].point();
     }
     return sol;
 }
@@ -261,9 +194,8 @@ void extrapolate_join(join_data jd)
     bool inc_ls = !circle1.center().isFinite();
     bool out_ls = !circle2.center().isFinite();
 
-    Geom::Point points[2];
+    std::vector<Geom::ShapeIntersection> points;
 
-    int solutions = 0;
     Geom::EllipticalArc *arc1 = NULL;
     Geom::EllipticalArc *arc2 = NULL;
     Geom::Point sol;
@@ -272,33 +204,29 @@ void extrapolate_join(join_data jd)
 
     if (!inc_ls && !out_ls) {
         // Two circles
-        solutions = Geom::circle_circle_intersection(circle1.center(), circle1.radius(),
-                                                     circle2.center(), circle2.radius(),
-                                                     points[0], points[1]);
-        if (solutions == 2) {
+        points = circle1.intersect(circle2);
+        if (points.size() == 2) {
             sol = pick_solution(points, tang2, endPt);
-            arc1 = circle1.arc(startPt, 0.5*(startPt+sol), sol, true);
-            arc2 = circle2.arc(sol, 0.5*(sol+endPt), endPt, true);
+            arc1 = circle1.arc(startPt, 0.5*(startPt+sol), sol);
+            arc2 = circle2.arc(sol, 0.5*(sol+endPt), endPt);
         }
     } else if (inc_ls && !out_ls) {
         // Line and circle
-        solutions = Geom::circle_line_intersection(circle2, incoming.initialPoint(), incoming.finalPoint(), points[0], points[1]);
-
-        if (solutions == 2) {
+        points = circle2.intersect(Line(incoming.initialPoint(), incoming.finalPoint()));
+        if (points.size() == 2) {
             sol = pick_solution(points, tang2, endPt);
-            arc2 = circle2.arc(sol, 0.5*(sol+endPt), endPt, true);
+            arc2 = circle2.arc(sol, 0.5*(sol+endPt), endPt);
         }
     } else if (!inc_ls && out_ls) {
         // Circle and line
-        solutions = Geom::circle_line_intersection(circle1, outgoing.initialPoint(), outgoing.finalPoint(), points[0], points[1]);
-        
-        if (solutions == 2) {
+        points = circle1.intersect(Line(outgoing.initialPoint(), outgoing.finalPoint()));
+        if (points.size() == 2) {
             sol = pick_solution(points, tang2, endPt);
-            arc1 = circle1.arc(startPt, 0.5*(sol+startPt), sol, true);
+            arc1 = circle1.arc(startPt, 0.5*(sol+startPt), sol);
         }
     }
 
-    if (solutions != 2)
+    if (points.size() != 2)
         // no solutions available, fall back to miter
         return miter_clip_join(jd);
 
@@ -348,19 +276,18 @@ void extrapolate_join(join_data jd)
             limit_line.setAngle(start_ray.angle() + limit_angle);
         }
 
-        Geom::EllipticalArc *arc_center = circle_center.arc(point_on_path, 0.5*(point_on_path + sol), sol, true);
+        Geom::EllipticalArc *arc_center = circle_center.arc(point_on_path, 0.5*(point_on_path + sol), sol);
         if (arc_center && arc_center->sweepAngle() > limit_angle) {
             // We need to clip
             clipped = true;
 
             if (!inc_ls) {
                 // Incoming circular
-                solutions = Geom::circle_line_intersection(circle1, limit_line.pointAt(0), limit_line.pointAt(1), points[0], points[1]);
-        
-                if (solutions == 2) {
+                points = circle1.intersect(limit_line);
+                if (points.size() == 2) {
                     p1 = pick_solution(points, tang2, endPt);
                     delete arc1;
-                    arc1 = circle1.arc(startPt, 0.5*(p1+startPt), p1, true);
+                    arc1 = circle1.arc(startPt, 0.5*(p1+startPt), p1);
                 }
             } else {
                 p1 = Geom::intersection_point(startPt, tang1, limit_line.pointAt(0), limit_line.versor());
@@ -368,12 +295,11 @@ void extrapolate_join(join_data jd)
 
             if (!out_ls) {
                 // Outgoing circular
-                solutions = Geom::circle_line_intersection(circle2, limit_line.pointAt(0), limit_line.pointAt(1), points[0], points[1]);
-        
-                if (solutions == 2) {
+                points = circle2.intersect(limit_line);
+                if (points.size() == 2) {
                     p2 = pick_solution(points, tang1, endPt);
                     delete arc2;
-                    arc2 = circle2.arc(p2, 0.5*(p2+endPt), endPt, true);
+                    arc2 = circle2.arc(p2, 0.5*(p2+endPt), endPt);
                 }
             } else {
                 p2 = Geom::intersection_point(endPt, tang2, limit_line.pointAt(0), limit_line.versor());
diff --git a/src/livarot/PathCutting.cpp b/src/livarot/PathCutting.cpp
index d6ad6c94a..086b30557 100644
--- a/src/livarot/PathCutting.cpp
+++ b/src/livarot/PathCutting.cpp
@@ -308,7 +308,7 @@ Path::MakePathVector()
             {
                 /* TODO: add testcase for this descr_arcto case */
                 PathDescrArcTo *nData = dynamic_cast<PathDescrArcTo *>(descr_cmd[i]);
-                currentpath->appendNew<Geom::SVGEllipticalArc>( nData->rx, nData->ry, nData->angle*M_PI/180.0, nData->large, !nData->clockwise, nData->p );
+                currentpath->appendNew<Geom::EllipticalArc>( nData->rx, nData->ry, nData->angle*M_PI/180.0, nData->large, !nData->clockwise, nData->p );
                 lastP = nData->p;
             }
             break;
@@ -400,11 +400,11 @@ void  Path::AddCurve(Geom::Curve const &c)
         Geom::Point tme = 3 * ((*cubic_bezier)[3] - (*cubic_bezier)[2]);
         CubicTo (tmp, tms, tme);
     }
-    else if(Geom::SVGEllipticalArc const *svg_elliptical_arc = dynamic_cast<Geom::SVGEllipticalArc const *>(&c)) {
-        ArcTo( svg_elliptical_arc->finalPoint(),
-               svg_elliptical_arc->ray(Geom::X), svg_elliptical_arc->ray(Geom::Y),
-               svg_elliptical_arc->rotationAngle()*180.0/M_PI,  // convert from radians to degrees
-               svg_elliptical_arc->largeArc(), !svg_elliptical_arc->sweep() );
+    else if(Geom::EllipticalArc const *elliptical_arc = dynamic_cast<Geom::EllipticalArc const *>(&c)) {
+        ArcTo( elliptical_arc->finalPoint(),
+               elliptical_arc->ray(Geom::X), elliptical_arc->ray(Geom::Y),
+               elliptical_arc->rotationAngle()*180.0/M_PI,  // convert from radians to degrees
+               elliptical_arc->largeArc(), !elliptical_arc->sweep() );
     } else { 
         //this case handles sbasis as well as all other curve types
         Geom::Path sbasis_path = Geom::cubicbezierpath_from_sbasis(c.toSBasis(), 0.1);
diff --git a/src/live_effects/lpe-circle_3pts.cpp b/src/live_effects/lpe-circle_3pts.cpp
index cbabb5a6c..dbb1f4b6b 100644
--- a/src/live_effects/lpe-circle_3pts.cpp
+++ b/src/live_effects/lpe-circle_3pts.cpp
@@ -17,6 +17,7 @@
 // You might need to include other 2geom files. You can add them here:
 #include <2geom/path.h>
 #include <2geom/circle.h>
+#include <2geom/path-sink.h>
 
 namespace Inkscape {
 namespace LivePathEffect {
@@ -47,7 +48,7 @@ static void _circle3(Geom::Point const &A, Geom::Point const &B, Geom::Point con
     double radius = L2(M - A);
 
     Geom::Circle c(M, radius);
-    c.getPath(path_out);
+    path_out = Geom::Path(c);
 }
 
 Geom::PathVector
diff --git a/src/live_effects/lpe-circle_with_radius.cpp b/src/live_effects/lpe-circle_with_radius.cpp
index 32d0c0bf5..8f2156044 100644
--- a/src/live_effects/lpe-circle_with_radius.cpp
+++ b/src/live_effects/lpe-circle_with_radius.cpp
@@ -15,11 +15,9 @@
 #include "display/curve.h"
 
 // You might need to include other 2geom files. You can add them here:
-#include <2geom/path.h>
-#include <2geom/sbasis.h>
-#include <2geom/bezier-to-sbasis.h>
-#include <2geom/d2.h>
+#include <2geom/pathvector.h>
 #include <2geom/circle.h>
+#include <2geom/path-sink.h>
 
 using namespace Geom;
 
@@ -51,9 +49,7 @@ LPECircleWithRadius::doEffect_path (Geom::PathVector const & path_in)
     double radius = Geom::L2(pt - center);
 
     Geom::Circle c(center, radius);
-    c.getPath(path_out);
-
-    return path_out;
+    return Geom::Path(c);
 }
 
 /*
diff --git a/src/live_effects/lpe-fillet-chamfer.cpp b/src/live_effects/lpe-fillet-chamfer.cpp
index 992c45a43..a9a96864a 100644
--- a/src/live_effects/lpe-fillet-chamfer.cpp
+++ b/src/live_effects/lpe-fillet-chamfer.cpp
@@ -16,7 +16,7 @@
 #include "live_effects/lpe-fillet-chamfer.h"
  
 #include <2geom/sbasis-to-bezier.h>
-#include <2geom/svg-elliptical-arc.h>
+#include <2geom/elliptical-arc.h>
 #include <2geom/line.h>
 
 #include "desktop.h"
@@ -582,7 +582,7 @@ LPEFilletChamfer::doEffect_path(Geom::PathVector const &path_in)
                     Geom::Path path_chamfer;
                     path_chamfer.start(path_out.finalPoint());
                     if((is_straight_curve(*curve_it1) && is_straight_curve(*curve_it2Fixed) && method != FM_BEZIER )|| method == FM_ARC){ 
-                        path_chamfer.appendNew<SVGEllipticalArc>(rx, ry, angleArc, 0, ccwToggle, endArcPoint);
+                        path_chamfer.appendNew<EllipticalArc>(rx, ry, angleArc, 0, ccwToggle, endArcPoint);
                     } else {
                         path_chamfer.appendNew<Geom::CubicBezier>(handle1, handle2, endArcPoint);
                     }
@@ -598,7 +598,7 @@ LPEFilletChamfer::doEffect_path(Geom::PathVector const &path_in)
                     path_chamfer.start(path_out.finalPoint());
                     if((is_straight_curve(*curve_it1) && is_straight_curve(*curve_it2Fixed) && method != FM_BEZIER )|| method == FM_ARC){ 
                         ccwToggle = ccwToggle?0:1;
-                        path_chamfer.appendNew<SVGEllipticalArc>(rx, ry, angleArc, 0, ccwToggle, endArcPoint);
+                        path_chamfer.appendNew<EllipticalArc>(rx, ry, angleArc, 0, ccwToggle, endArcPoint);
                     }else{
                         path_chamfer.appendNew<Geom::CubicBezier>(inverseHandle1, inverseHandle2, endArcPoint);
                     }
@@ -611,13 +611,13 @@ LPEFilletChamfer::doEffect_path(Geom::PathVector const &path_in)
                 } else if (type == 2) {
                     if((is_straight_curve(*curve_it1) && is_straight_curve(*curve_it2Fixed) && method != FM_BEZIER )|| method == FM_ARC){ 
                         ccwToggle = ccwToggle?0:1;
-                        path_out.appendNew<SVGEllipticalArc>(rx, ry, angleArc, 0, ccwToggle, endArcPoint);
+                        path_out.appendNew<EllipticalArc>(rx, ry, angleArc, 0, ccwToggle, endArcPoint);
                     }else{
                         path_out.appendNew<Geom::CubicBezier>(inverseHandle1, inverseHandle2, endArcPoint);
                     }
                 } else if (type == 1){
                     if((is_straight_curve(*curve_it1) && is_straight_curve(*curve_it2Fixed) && method != FM_BEZIER )|| method == FM_ARC){ 
-                        path_out.appendNew<SVGEllipticalArc>(rx, ry, angleArc, 0, ccwToggle, endArcPoint);
+                        path_out.appendNew<EllipticalArc>(rx, ry, angleArc, 0, ccwToggle, endArcPoint);
                     } else {
                         path_out.appendNew<Geom::CubicBezier>(handle1, handle2, endArcPoint);
                     }
diff --git a/src/live_effects/lpe-jointype.cpp b/src/live_effects/lpe-jointype.cpp
index 893f3ab4e..cc5587194 100644
--- a/src/live_effects/lpe-jointype.cpp
+++ b/src/live_effects/lpe-jointype.cpp
@@ -20,7 +20,7 @@
 #include "display/curve.h"
 
 #include <2geom/path.h>
-#include <2geom/svg-elliptical-arc.h>
+#include <2geom/elliptical-arc.h>
 
 #include "lpe-jointype.h"
 
diff --git a/src/live_effects/lpe-offset.cpp b/src/live_effects/lpe-offset.cpp
index c841a5442..dd7af92a2 100644
--- a/src/live_effects/lpe-offset.cpp
+++ b/src/live_effects/lpe-offset.cpp
@@ -20,7 +20,7 @@
 #include <2geom/path.h>
 #include <2geom/piecewise.h>
 #include <2geom/sbasis-geometric.h>
-#include <2geom/svg-elliptical-arc.h>
+#include <2geom/elliptical-arc.h>
 #include <2geom/transforms.h>
 
 namespace Inkscape {
@@ -52,7 +52,7 @@ static void append_half_circle(Geom::Piecewise<Geom::D2<Geom::SBasis> > &pwd2,
     using namespace Geom;
 
     double r = L2(dir);
-    SVGEllipticalArc cap(center + dir, r, r, angle_between(Point(1,0), dir), false, false, center - dir);
+    EllipticalArc cap(center + dir, r, r, angle_between(Point(1,0), dir), false, false, center - dir);
     Piecewise<D2<SBasis> > cap_pwd2(cap.toSBasis());
     pwd2.continuousConcat(cap_pwd2);
 }
diff --git a/src/live_effects/lpe-powerstroke.cpp b/src/live_effects/lpe-powerstroke.cpp
index 6a823e943..6ce912bcb 100644
--- a/src/live_effects/lpe-powerstroke.cpp
+++ b/src/live_effects/lpe-powerstroke.cpp
@@ -27,12 +27,13 @@
 #include <2geom/sbasis-geometric.h>
 #include <2geom/transforms.h>
 #include <2geom/bezier-utils.h>
-#include <2geom/svg-elliptical-arc.h>
+#include <2geom/elliptical-arc.h>
 #include <2geom/sbasis-to-bezier.h>
 #include <2geom/path-sink.h>
 #include <2geom/path-intersection.h>
 #include <2geom/crossing.h>
 #include <2geom/ellipse.h>
+#include <2geom/circle.h>
 #include <2geom/math-utils.h>
 #include <math.h>
 
@@ -94,71 +95,6 @@ static Ellipse find_ellipse(Point P, Point Q, Point O)
 }
 
 /**
- * Refer to: Weisstein, Eric W. "Circle-Circle Intersection."
-             From MathWorld--A Wolfram Web Resource.
-             http://mathworld.wolfram.com/Circle-CircleIntersection.html
- *
- * @return 0 if no intersection
- * @return 1 if one circle is contained in the other
- * @return 2 if intersections are found (they are written to p0 and p1)
- */
-static int circle_circle_intersection(Circle const &circle0, Circle const &circle1,
-                                      Point & p0, Point & p1)
-{
-    Point X0 = circle0.center();
-    double r0 = circle0.radius();
-    Point X1 = circle1.center();
-    double r1 = circle1.radius();
-
-    /* dx and dy are the vertical and horizontal distances between
-    * the circle centers.
-    */
-    Point D = X1 - X0;
-
-    /* Determine the straight-line distance between the centers. */
-    double d = L2(D);
-
-    /* Check for solvability. */
-    if (d > (r0 + r1))
-    {
-        /* no solution. circles do not intersect. */
-        return 0;
-    }
-    if (d <= fabs(r0 - r1))
-    {
-        /* no solution. one circle is contained in the other */
-        return 1;
-    }
-
-    /* 'point 2' is the point where the line through the circle
-    * intersection points crosses the line between the circle
-    * centers.  
-    */
-
-    /* Determine the distance from point 0 to point 2. */
-    double a = ((r0*r0) - (r1*r1) + (d*d)) / (2.0 * d) ;
-
-    /* Determine the coordinates of point 2. */
-    Point p2 = X0 + D * (a/d);
-
-    /* Determine the distance from point 2 to either of the
-    * intersection points.
-    */
-    double h = std::sqrt((r0*r0) - (a*a));
-
-    /* Now determine the offsets of the intersection points from
-    * point 2.
-    */
-    Point r = (h/d)*rot90(D);
-
-    /* Determine the absolute intersection points. */
-    p0 = p2 + r;
-    p1 = p2 - r;
-
-    return 2;
-}
-
-/**
  * Find circle that touches inside of the curve, with radius matching the curvature, at time value \c t.
  * Because this method internally uses unitTangentAt, t should be smaller than 1.0 (see unitTangentAt).
  */
@@ -484,24 +420,24 @@ static Geom::Path path_from_piecewise_fix_cusps( Geom::Piecewise<Geom::D2<Geom::
                     // Extrapolate using the curvature at the end of the path segments to join
                     Geom::Circle circle1 = Geom::touching_circle(reverse(B[prev_i]), 0.0);
                     Geom::Circle circle2 = Geom::touching_circle(B[i], 0.0);
-                    Geom::Point points[2];
-                    int solutions = circle_circle_intersection(circle1, circle2, points[0], points[1]);
-                    if (solutions == 2) {
+                    std::vector<Geom::ShapeIntersection> solutions;
+                    solutions = circle1.intersect(circle2);
+                    if (solutions.size() == 2) {
                         Geom::Point sol(0.,0.);
-                        if ( dot(tang2,points[0]-B[i].at0()) > 0 ) {
+                        if ( dot(tang2, solutions[0].point() - B[i].at0()) > 0 ) {
                             // points[0] is bad, choose points[1]
-                            sol = points[1];
-                        } else if ( dot(tang2,points[1]-B[i].at0()) > 0 ) { // points[0] could be good, now check points[1]
+                            sol = solutions[1].point();
+                        } else if ( dot(tang2, solutions[1].point() - B[i].at0()) > 0 ) { // points[0] could be good, now check points[1]
                             // points[1] is bad, choose points[0]
-                            sol = points[0];
+                            sol = solutions[0].point();
                         } else {
                             // both points are good, choose nearest
-                            sol = ( distanceSq(B[i].at0(), points[0]) < distanceSq(B[i].at0(), points[1]) ) ?
-                                    points[0] : points[1];
+                            sol = ( distanceSq(B[i].at0(), solutions[0].point()) < distanceSq(B[i].at0(), solutions[1].point()) ) ?
+                                    solutions[0].point() : solutions[1].point();
                         }
 
-                        Geom::EllipticalArc *arc0 = circle1.arc(B[prev_i].at1(), 0.5*(B[prev_i].at1()+sol), sol, true);
-                        Geom::EllipticalArc *arc1 = circle2.arc(sol, 0.5*(sol+B[i].at0()), B[i].at0(), true);
+                        Geom::EllipticalArc *arc0 = circle1.arc(B[prev_i].at1(), 0.5*(B[prev_i].at1()+sol), sol);
+                        Geom::EllipticalArc *arc1 = circle2.arc(sol, 0.5*(sol+B[i].at0()), B[i].at0());
 
                         if (arc0) {
                             build_from_sbasis(pb,arc0->toSBasis(), tol, false);
@@ -739,7 +675,7 @@ LPEPowerStroke::doEffect_path (Geom::PathVector const & path_in)
             default:
             {
                 double radius1 = 0.5 * distance(pwd2_out.lastValue(), mirrorpath.firstValue());
-                fixed_path.appendNew<SVGEllipticalArc>( radius1, radius1, M_PI/2., false, y.lastValue() < 0, mirrorpath.firstValue() );
+                fixed_path.appendNew<EllipticalArc>( radius1, radius1, M_PI/2., false, y.lastValue() < 0, mirrorpath.firstValue() );
                 break;
             }
         }
@@ -777,7 +713,7 @@ LPEPowerStroke::doEffect_path (Geom::PathVector const & path_in)
             default:
             {
                 double radius2 = 0.5 * distance(pwd2_out.firstValue(), mirrorpath.lastValue());
-                fixed_path.appendNew<SVGEllipticalArc>( radius2, radius2, M_PI/2., false, y.firstValue() < 0, pwd2_out.firstValue() );
+                fixed_path.appendNew<EllipticalArc>( radius2, radius2, M_PI/2., false, y.firstValue() < 0, pwd2_out.firstValue() );
                 break;
             }
         }
diff --git a/src/live_effects/parameter/filletchamferpointarray.cpp b/src/live_effects/parameter/filletchamferpointarray.cpp
index a3ff4c96f..955fcd621 100644
--- a/src/live_effects/parameter/filletchamferpointarray.cpp
+++ b/src/live_effects/parameter/filletchamferpointarray.cpp
@@ -11,7 +11,6 @@
 #include <2geom/piecewise.h>
 #include <2geom/sbasis-to-bezier.h>
 #include <2geom/sbasis-geometric.h>
-#include <2geom/svg-elliptical-arc.h>
 #include <2geom/line.h>
 #include <2geom/path-intersection.h>
 
diff --git a/src/object-snapper.cpp b/src/object-snapper.cpp
index 1293d19aa..634d56aa6 100644
--- a/src/object-snapper.cpp
+++ b/src/object-snapper.cpp
@@ -19,6 +19,7 @@
 #include <2geom/rect.h>
 #include <2geom/line.h>
 #include <2geom/circle.h>
+#include <2geom/path-sink.h>
 #include "document.h"
 #include "sp-namedview.h"
 #include "sp-image.h"
@@ -626,7 +627,10 @@ void Inkscape::ObjectSnapper::_snapPathsConstrained(IntermSnapResults &isr,
     Geom::PathVector constraint_path;
     if (c.isCircular()) {
         Geom::Circle constraint_circle(dt->dt2doc(c.getPoint()), c.getRadius());
-        constraint_circle.getPath(constraint_path);
+        Geom::PathBuilder pb;
+        pb.feed(constraint_circle);
+        pb.flush();
+        constraint_path = pb.peek();
     } else {
         Geom::Path constraint_line;
         constraint_line.start(p_min_on_cl);
diff --git a/src/sp-ellipse.cpp b/src/sp-ellipse.cpp
index 932a3a1b7..0675b7781 100644
--- a/src/sp-ellipse.cpp
+++ b/src/sp-ellipse.cpp
@@ -18,6 +18,7 @@
 #include <glibmm/i18n.h>
 
 #include <2geom/angle.h>
+#include <2geom/circle.h>
 #include <2geom/ellipse.h>
 #include <2geom/path-sink.h>
 #include <2geom/pathvector.h>
diff --git a/src/svg/svg-path.cpp b/src/svg/svg-path.cpp
index ab51a5537..ba9e11452 100644
--- a/src/svg/svg-path.cpp
+++ b/src/svg/svg-path.cpp
@@ -82,11 +82,11 @@ static void sp_svg_write_curve(Inkscape::SVG::PathString & str, Geom::Curve cons
                      (*cubic_bezier)[2][0], (*cubic_bezier)[2][1],
                      (*cubic_bezier)[3][0], (*cubic_bezier)[3][1] );
     }
-    else if(Geom::SVGEllipticalArc const *svg_elliptical_arc = dynamic_cast<Geom::SVGEllipticalArc const *>(c)) {
-        str.arcTo( svg_elliptical_arc->ray(Geom::X), svg_elliptical_arc->ray(Geom::Y),
-                   Geom::rad_to_deg(svg_elliptical_arc->rotationAngle()),
-                   svg_elliptical_arc->largeArc(), svg_elliptical_arc->sweep(),
-                   svg_elliptical_arc->finalPoint() );
+    else if(Geom::EllipticalArc const *elliptical_arc = dynamic_cast<Geom::EllipticalArc const *>(c)) {
+        str.arcTo( elliptical_arc->ray(Geom::X), elliptical_arc->ray(Geom::Y),
+                   Geom::rad_to_deg(elliptical_arc->rotationAngle()),
+                   elliptical_arc->largeArc(), elliptical_arc->sweep(),
+                   elliptical_arc->finalPoint() );
     } else { 
         //this case handles sbasis as well as all other curve types
         Geom::Path sbasis_path = Geom::cubicbezierpath_from_sbasis(c->toSBasis(), 0.1);
diff --git a/src/ui/tools/calligraphic-tool.cpp b/src/ui/tools/calligraphic-tool.cpp
index 15e6527a3..28195eb75 100644
--- a/src/ui/tools/calligraphic-tool.cpp
+++ b/src/ui/tools/calligraphic-tool.cpp
@@ -140,8 +140,7 @@ void CalligraphicTool::setup() {
 
     {
         /* TODO: have a look at DropperTool::setup where the same is done.. generalize? */
-        Geom::PathVector path;
-        Geom::Circle(0, 0, 1).getPath(path);
+        Geom::PathVector path = Geom::Path(Geom::Circle(0,0,1));
 
         SPCurve *c = new SPCurve(path);
 
diff --git a/src/ui/tools/dropper-tool.cpp b/src/ui/tools/dropper-tool.cpp
index bda9d8e8a..c838c27d5 100644
--- a/src/ui/tools/dropper-tool.cpp
+++ b/src/ui/tools/dropper-tool.cpp
@@ -84,8 +84,7 @@ void DropperTool::setup() {
     ToolBase::setup();
 
     /* TODO: have a look at CalligraphicTool::setup where the same is done.. generalize? */
-    Geom::PathVector path;
-    Geom::Circle(0, 0, 1).getPath(path);
+    Geom::PathVector path = Geom::Path(Geom::Circle(0,0,1));
 
     SPCurve *c = new SPCurve(path);
 
diff --git a/src/ui/tools/spray-tool.cpp b/src/ui/tools/spray-tool.cpp
index 26d74733a..6e3a06345 100644
--- a/src/ui/tools/spray-tool.cpp
+++ b/src/ui/tools/spray-tool.cpp
@@ -207,8 +207,7 @@ void SprayTool::setup() {
 
     {
         /* TODO: have a look at sp_dyna_draw_context_setup where the same is done.. generalize? at least make it an arcto! */
-        Geom::PathVector path;
-        Geom::Circle(0, 0, 1).getPath(path);
+        Geom::PathVector path = Geom::Path(Geom::Circle(0,0,1));
 
         SPCurve *c = new SPCurve(path);
 
diff --git a/src/ui/tools/tweak-tool.cpp b/src/ui/tools/tweak-tool.cpp
index 76b52f9be..94f7aa135 100644
--- a/src/ui/tools/tweak-tool.cpp
+++ b/src/ui/tools/tweak-tool.cpp
@@ -259,8 +259,7 @@ void TweakTool::setup() {
 
     {
         /* TODO: have a look at sp_dyna_draw_context_setup where the same is done.. generalize? at least make it an arcto! */
-        Geom::PathVector path;
-        Geom::Circle(0, 0, 1).getPath(path);
+        Geom::PathVector path = Geom::Path(Geom::Circle(0,0,1));
 
         SPCurve *c = new SPCurve(path);