Update 2Geom to pull in integer rectangle class

(bzr r10347.1.1)

84 files changed, 2751 insertions, 1265 deletions

diff --git a/src/2geom/Makefile_insert b/src/2geom/Makefile_insert
index 4f7c3b6ef..a668a2b3b 100644
--- a/src/2geom/Makefile_insert
+++ b/src/2geom/Makefile_insert
@@ -78,6 +78,7 @@
 	2geom/quadtree.h	\
 	2geom/ray.h	\
 	2geom/rect.h	\
+	2geom/rect.cpp	\
 	2geom/region.cpp	\
 	2geom/region.h	\
 	2geom/sbasis-2d.cpp	\
diff --git a/src/2geom/affine.cpp b/src/2geom/affine.cpp
index 925f43820..2a1f18d77 100644
--- a/src/2geom/affine.cpp
+++ b/src/2geom/affine.cpp
@@ -1,9 +1,3 @@
-#define __Geom_MATRIX_C__
-
-/** \file
- * Various matrix routines.  Currently includes some Geom::Rotate etc. routines too.
- */
-
 /*
  * Authors:
  *   Lauris Kaplinski <lauris@kaplinski.com>
@@ -387,10 +381,10 @@ Coord Affine::descrim2() const {
 }
 
 /** @brief Calculate the descriminant.
- * If the matrix doesn't contain a non-uniform scaling or shearing component, this value says
- * how will the length any line segment change after applying this transformation
- * to arbitrary objects on a plane (the new length will be
- * @code line_seg.length() * m.descrim()) @endcode.
+ * If the matrix doesn't contain a shearing or non-uniform scaling component, this value says
+ * how will the length of any line segment change after applying this transformation
+ * to arbitrary objects on a plane. The new length will be
+ * @code line_seg.length() * m.descrim()) @endcode
  * @return \f$\sqrt{|\det A|}\f$. */
 Coord Affine::descrim() const {
     return sqrt(descrim2());
diff --git a/src/2geom/affine.h b/src/2geom/affine.h
index 277d8b4ee..b07fba0f7 100644
--- a/src/2geom/affine.h
+++ b/src/2geom/affine.h
@@ -1,7 +1,8 @@
-/** \file
- *  \brief 3x3 affine transformation matrix.
+/** 
+ * \file
+ * \brief 3x3 affine transformation matrix.
  *//*
- * Main authors:
+ * Authors:
  *   Lauris Kaplinski <lauris@kaplinski.com> (Original NRAffine definition and related macros)
  *   Nathan Hurst <njh@mail.csse.monash.edu.au> (Geom::Affine class version of the above)
  *   Michael G. Sloan <mgsloan@gmail.com> (reorganization and additions)
@@ -10,8 +11,8 @@
  * This code is in public domain.
  */
 
-#ifndef SEEN_LIB2GEOM_MATRIX_H
-#define SEEN_LIB2GEOM_MATRIX_H
+#ifndef SEEN_LIB2GEOM_AFFINE_H
+#define SEEN_LIB2GEOM_AFFINE_H
 
 #include <boost/operators.hpp>
 #include <2geom/forward.h>
@@ -236,9 +237,9 @@ inline Affine Affine::identity() {
     return ret; // allow NRVO
 }
 
-} /* namespace Geom */
+} // end namespace Geom
 
-#endif /* !__Geom_MATRIX_H__ */
+#endif // LIB2GEOM_SEEN_AFFINE_H
 
 /*
   Local Variables:
diff --git a/src/2geom/angle.h b/src/2geom/angle.h
index 42e3531f3..bdf546989 100644
--- a/src/2geom/angle.h
+++ b/src/2geom/angle.h
@@ -107,11 +107,18 @@ public:
         if (ret < 0) ret += 360;
         return ret;
     }
-
+    /** @brief Create an angle from its measure in radians. */
+    static Angle from_radians(Coord d) {
+        Angle a(d);
+        return a;
+    }
+    /** @brief Create an angle from its measure in degrees. */
     static Angle from_degrees(Coord d) {
         Angle a(d * M_PI / 180);
         return a;
     }
+    /** @brief Create an angle from its measure in degrees in clock convention.
+     * @see Angle::degreesClock() */
     static Angle from_degrees_clock(Coord d) {
         // first make sure d is in [0, 360)
         d = std::fmod(d, 360.0);
@@ -208,9 +215,12 @@ protected:
     bool _sweep;
 };
 
-inline double deg_to_rad(double deg) { return deg*M_PI/180.0;}
-
-inline double rad_to_deg(double rad) { return rad*180.0/M_PI;}
+/** @brief Given an angle in degrees, return radians
+ * @relates Angle */
+inline Coord deg_to_rad(Coord deg) { return deg*M_PI/180.0;}
+/** @brief Given an angle in radians, return degrees
+ * @relates Angle */
+inline Coord rad_to_deg(Coord rad) { return rad*180.0/M_PI;}
 
 /*
  *  start_angle and angle must belong to [0, 2PI[
@@ -294,8 +304,7 @@ bool arc_contains (double a, double sa, double ia, double ea)
 
 } // end namespace Geom
 
-#endif
-
+#endif // LIB2GEOM_SEEN_ANGLE_H
 
 /*
   Local Variables:
diff --git a/src/2geom/basic-intersection.h b/src/2geom/basic-intersection.h
index b07052449..5a813ae99 100644
--- a/src/2geom/basic-intersection.h
+++ b/src/2geom/basic-intersection.h
@@ -1,6 +1,6 @@
 /**
  * \file
- * \brief  \todo brief description
+ * \brief Basic intersection routines
  *
  * Authors:
  *      ? <?@?.?>
diff --git a/src/2geom/bezier-curve.cpp b/src/2geom/bezier-curve.cpp
index bde0e3ef1..46aff8b49 100644
--- a/src/2geom/bezier-curve.cpp
+++ b/src/2geom/bezier-curve.cpp
@@ -1,8 +1,5 @@
-/**
- * \file
- * \brief Bezier curve
+/* Bezier curve implementation
  *
- *//*
  * Authors:
  *   MenTaLguY <mental@rydia.net>
  *   Marco Cecchetti <mrcekets at gmail.com>
@@ -41,7 +38,7 @@ namespace Geom
 
 /**
  * @class BezierCurve
- * @brief Two-dimensional Bezier curve of arbitrary order. (this is an abstract class)
+ * @brief Two-dimensional Bezier curve of arbitrary order.
  *
  * Bezier curves are an expansion of the concept of linear interpolation to n points.
  * Linear segments in 2Geom are in fact Bezier curves of order 1.
@@ -82,28 +79,40 @@ namespace Geom
  * have an intutive geometric interpretation. Because of this, they are frequently used
  * in vector graphics editors.
  *
- * Every bezier curve is contained in its control polygon (the convex polygon composed
+ * Every Bezier curve is contained in its control polygon (the convex polygon composed
  * of its control points). This fact is useful for sweepline algorithms and intersection.
  *
- * Bezier curves of order 1, 2 and 3 are common enough to have their own more specific subclasses:
- * LineSegment, QuadraticBezier, and CubicBezier.
- * Note that you cannot create a generic BezierCurve, you can only create a BezierCurve of a
- * specific order, by creating a BezierCurveN.
+ * @par Implementation notes
+ * The order of a Bezier curve is immuable once it has been created. Normally, you should
+ * know the order at compile time and use the BezierCurveN template. If you need to determine
+ * the order at runtime, use the BezierCurve::create() function. It will create a BezierCurveN
+ * for orders 1, 2 and 3 (up to cubic Beziers), so you can later <tt>dynamic_cast</tt>
+ * to those types, and for higher orders it will create an instance of BezierCurve.
  *
  * @relates BezierCurveN
  * @ingroup Curves
  */
 
- /**
+/**
  * @class BezierCurveN
+ * @brief Bezier curve with compile-time specified order.
+ *
  * @tparam degree unsigned value indicating the order of the bezier curve
- * @brief Two-dimensional Bezier curve of specific order.
  * 
  * @relates BezierCurve 
  * @ingroup Curves
  */
 
- 
+
+BezierCurve::BezierCurve(std::vector<Point> const &pts)
+{
+    inner = D2<Bezier>(Bezier::Order(pts.size()-1), Bezier::Order(pts.size()-1));
+    for (unsigned d = 0; d < 2; ++d) {
+        for(unsigned i = 0; i <= pts.size(); i++)
+            inner[d][i] = pts[i][d];
+    }
+}
+
 Coord BezierCurve::length(Coord tolerance) const
 {
     switch (order())
@@ -127,6 +136,48 @@ Coord BezierCurve::length(Coord tolerance) const
     }
 }
 
+BezierCurve *BezierCurve::create(std::vector<Point> const &pts)
+{
+    switch (pts.size()) {
+    case 0:
+    case 1:
+        THROW_LOGICALERROR("BezierCurve::create: too few points in vector");
+        return NULL;
+    case 2:
+        return new LineSegment(pts[0], pts[1]);
+    case 3:
+        return new QuadraticBezier(pts[0], pts[1], pts[2]);
+    case 4:
+        return new CubicBezier(pts[0], pts[1], pts[2], pts[3]);
+    default:
+        return new BezierCurve(pts);
+    }
+}
+
+// optimized specializations for LineSegment
+
+template <>
+Curve *BezierCurveN<1>::derivative() const {
+    double dx = inner[X][1] - inner[X][0], dy = inner[Y][1] - inner[Y][0];
+    return new BezierCurveN<1>(Point(dx,dy),Point(dx,dy));
+}
+
+template<>
+Coord BezierCurveN<1>::nearestPoint(Point const& p, Coord from, Coord to) const
+{
+    if ( from > to ) std::swap(from, to);
+    Point ip = pointAt(from);
+    Point fp = pointAt(to);
+    Point v = fp - ip;
+    Coord l2v = L2sq(v);
+    if (l2v == 0) return 0;
+    Coord t = dot( p - ip, v ) / l2v;
+    if ( t <= 0 )  		return from;
+    else if ( t >= 1 )  return to;
+    else return from + t*(to-from);
+}
+
+
 static Coord bezier_length_internal(std::vector<Point> &v1, Coord tolerance)
 {
     /* The Bezier length algorithm used in 2Geom utilizes a simple fact:
@@ -178,7 +229,7 @@ static Coord bezier_length_internal(std::vector<Point> &v1, Coord tolerance)
      * After loop with i==2
      * # # 2 3 4
      *  # 1 ? ?
-     *   0 ? ? -> wirte 0 to v2[2]
+     *   0 ? ? -> write 0 to v2[2]
      *    ? ?
      *     ?
      *
@@ -204,7 +255,7 @@ static Coord bezier_length_internal(std::vector<Point> &v1, Coord tolerance)
 }
 
 /** @brief Compute the length of a bezier curve given by a vector of its control points
- * @relates BezierCurve */
+ * @relatesalso BezierCurve */
 Coord bezier_length(std::vector<Point> const &points, Coord tolerance)
 {
     if (points.size() < 2) return 0.0;
@@ -213,7 +264,7 @@ Coord bezier_length(std::vector<Point> const &points, Coord tolerance)
 }
 
 /** @brief Compute the length of a quadratic bezier curve given by its control points
- * @relates QuadraticBezier */
+ * @relatesalso QuadraticBezier */
 Coord bezier_length(Point a0, Point a1, Point a2, Coord tolerance)
 {
     Coord lower = distance(a0, a2);
@@ -231,7 +282,7 @@ Coord bezier_length(Point a0, Point a1, Point a2, Coord tolerance)
 }
 
 /** @brief Compute the length of a cubic bezier curve given by its control points
- * @relates CubicBezier */
+ * @relatesalso CubicBezier */
 Coord bezier_length(Point a0, Point a1, Point a2, Point a3, Coord tolerance)
 {
     Coord lower = distance(a0, a3);
diff --git a/src/2geom/bezier-curve.h b/src/2geom/bezier-curve.h
index 40da6f366..d13ff8321 100644
--- a/src/2geom/bezier-curve.h
+++ b/src/2geom/bezier-curve.h
@@ -1,14 +1,13 @@
 /**
  * \file
  * \brief Bezier curve
- *
  *//*
  * Authors:
  *   MenTaLguY <mental@rydia.net>
  *   Marco Cecchetti <mrcekets at gmail.com>
  *   Krzysztof Kosiński <tweenk.pl@gmail.com>
  * 
- * Copyright 2007-2009 Authors
+ * Copyright 2007-2011 Authors
  *
  * This library is free software; you can redistribute it and/or
  * modify it either under the terms of the GNU Lesser General Public
@@ -47,10 +46,13 @@ namespace Geom
 class BezierCurve : public Curve {
 protected:
     D2<Bezier> inner;
+    BezierCurve() {}
+    BezierCurve(BezierCurve const &b) : inner(b.inner) {}
+    BezierCurve(D2<Bezier> const &b) : inner(b) {}
+    BezierCurve(Bezier const &x, Bezier const &y) : inner(x, y) {}
+    BezierCurve(std::vector<Point> const &pts);
 
 public:
-   /// No constructors allowed!
-
     /// @name Access and modify control points
     /// @{
     /** @brief Get the order of the Bezier curve.
@@ -66,8 +68,10 @@ public:
         inner[X].setPoint(ix, v[X]);
         inner[Y].setPoint(ix, v[Y]);
     }
-    /** @brief Set new control points for this curve.
-     * @param ps Vector which must contain order() + 1 points.  Note that the caller is responsible for checking the size of this vector. */
+    /** @brief Set new control points.
+     * @param ps Vector which must contain order() + 1 points.
+     *           Note that the caller is responsible for checking the size of this vector.
+     * @throws LogicalError Thrown when the size of the vector does not match the order. */
     virtual void setPoints(std::vector<Point> const &ps) {
         // must be virtual, because HLineSegment will need to redefine it
         if (ps.size() != order() + 1)
@@ -76,12 +80,18 @@ public:
             setPoint(i, ps[i]);
         }
     }
-    /** Access control points of the curve.
-     * @param ix The (zero-based) index of the control point.  Note that the caller is responsible for checking that this value is <= order().
+    /** @brief Access control points of the curve.
+     * @param ix The (zero-based) index of the control point. Note that the caller is responsible for checking that this value is <= order().
      * @return The control point. No-reference return, use setPoint() to modify control points. */
     Point const operator[](unsigned ix) const { return Point(inner[X][ix], inner[Y][ix]); }
     /// @}
 
+    /// @name Construct a Bezier curve with runtime-determined order.
+    /// @{
+    /** @brief Construct a curve from a vector of control points. */
+    static BezierCurve *create(std::vector<Point> const &pts);
+    /// @}
+
     // implementation of virtual methods goes here
 #ifndef DOXYGEN_SHOULD_SKIP_THIS
     virtual Point initialPoint() const { return inner.at0(); }
@@ -100,6 +110,27 @@ public:
                                                    bounds_local(Geom::derivative(inner[Y]), i));
         return OptRect();
     }
+    virtual Curve *duplicate() const {
+        return new BezierCurve(*this);
+    }
+    virtual Curve *portion(Coord f, Coord t) const {
+        return new BezierCurve(Geom::portion(inner, f, t));
+    }
+    virtual Curve *reverse() const {
+        return new BezierCurve(Geom::reverse(inner));
+    }
+    virtual Curve *transformed(Affine const &m) const {
+        BezierCurve *ret = new BezierCurve();
+        std::vector<Point> ps = points();
+        for (unsigned i = 0;  i <= order(); i++) {
+            ps[i] = ps[i] * m;
+        }
+        ret->setPoints(ps);
+        return ret;
+    }
+    virtual Curve *derivative() const {
+        return new BezierCurve(Geom::derivative(inner[X]), Geom::derivative(inner[Y]));
+    }
     virtual int degreesOfFreedom() const {
         return 2 * (order() + 1);
     }
@@ -121,7 +152,7 @@ public:
   template <unsigned required_degree>
   static void assert_degree(BezierCurveN<required_degree> const *) {}
 
-    /// @name Construct the curve
+    /// @name Construct Bezier curves
     /// @{
     /** @brief Construct a Bezier curve of the specified order with all points zero. */
     BezierCurveN() {
@@ -175,7 +206,19 @@ public:
 
     /// @}
 
+    /** @brief Divide a Bezier curve into two curves
+     * @param t Time value
+     * @return Pair of Bezier curves \f$(\mathbf{D}, \mathbf{E})\f$ such that
+     *         \f$\mathbf{D}[ [0,1] ] = \mathbf{C}[ [0,t] ]\f$ and
+     *         \f$\mathbf{E}[ [0,1] ] = \mathbf{C}[ [t,1] ]\f$ */
+    std::pair<BezierCurveN, BezierCurveN> subdivide(Coord t) const {
+        std::pair<Bezier, Bezier> sx = inner[X].subdivide(t), sy = inner[Y].subdivide(t);
+        return std::make_pair(
+                   BezierCurveN(sx.first, sy.first),
+                   BezierCurveN(sx.second, sy.second));
+    }
 
+#ifndef DOXYGEN_SHOULD_SKIP_THIS
     virtual Curve *duplicate() const {
         return new BezierCurveN(*this);
     }
@@ -207,31 +250,29 @@ public:
         }
     }
     virtual Curve *derivative() const;
-
-    /** @brief Divide a Bezier curve into two curves
-     * @param t Time value
-     * @return Pair of Bezier curves \f$(\mathbf{D}, \mathbf{E})\f$ such that
-     *         \f$\mathbf{D}[ [0,1] ] = \mathbf{C}[ [0,t] ]\f$ and
-     *         \f$\mathbf{E}[ [0,1] ] = \mathbf{C}[ [t,1] ]\f$ */
-    std::pair<BezierCurveN, BezierCurveN> subdivide(Coord t) const {
-        std::pair<Bezier, Bezier> sx = inner[X].subdivide(t), sy = inner[Y].subdivide(t);
-        return std::make_pair(
-                   BezierCurveN(sx.first, sy.first),
-                   BezierCurveN(sx.second, sy.second));
-    }
-
-    double nearestPoint( Point const& p, double from = 0, double to = 1 ) const {
+    
+    // the method below is defined so that LineSegment can specialize it
+    virtual Coord nearestPoint(Point const& p, Coord from = 0, Coord to = 1) const {
         return Curve::nearestPoint(p, from, to);
     }
-
+#endif
 };
 
 // BezierCurveN<0> is meaningless; specialize it out
-template<> class BezierCurveN<0> : public BezierCurveN<1> { public: BezierCurveN();};
+template<> class BezierCurveN<0> : public BezierCurveN<1> { private: BezierCurveN();};
 
-// provide convenient names for common degree bezier curves
+/** @brief Line segment.
+ * Line segments are Bezier curves of order 1. They have only two control points,
+ * the starting point and the ending point.
+ * @ingroup Curves */
 typedef BezierCurveN<1> LineSegment;
+
+/** @brief Quadratic (order 2) Bezier curve.
+ * @ingroup Curves */
 typedef BezierCurveN<2> QuadraticBezier;
+
+/** @brief Cubic (order 3) Bezier curve.
+ * @ingroup Curves */
 typedef BezierCurveN<3> CubicBezier;
 
 template <unsigned degree>
@@ -239,28 +280,10 @@ inline
 Curve *BezierCurveN<degree>::derivative() const {
     return new BezierCurveN<degree-1>(Geom::derivative(inner[X]), Geom::derivative(inner[Y]));
 }
-template <>
-inline
-Curve *BezierCurveN<1>::derivative() const {
-    double dx = inner[X][1] - inner[X][0], dy = inner[Y][1] - inner[Y][0];
-    return new BezierCurveN<1>(Point(dx,dy),Point(dx,dy));
-}
 
-template<>
-inline
-double LineSegment::nearestPoint(Point const& p, double from, double to) const
-{
-    if ( from > to ) std::swap(from, to);
-    Point ip = pointAt(from);
-    Point fp = pointAt(to);
-    Point v = fp - ip;
-    Coord l2v = L2sq(v);
-    if (l2v == 0) return 0;
-    Coord t = dot( p - ip, v ) / l2v;
-    if ( t <= 0 )  		return from;
-    else if ( t >= 1 )  return to;
-    else return from + t*(to-from);
-}
+// optimized specializations for LineSegment
+template <> Curve *BezierCurveN<1>::derivative() const;
+template <> Coord BezierCurveN<1>::nearestPoint(Point const &, Coord, Coord) const;
 
 inline Point middle_point(LineSegment const& _segment) {
     return ( _segment.initialPoint() + _segment.finalPoint() ) / 2;
diff --git a/src/2geom/bezier-to-sbasis.h b/src/2geom/bezier-to-sbasis.h
index ba98a8a34..8cd4bf444 100644
--- a/src/2geom/bezier-to-sbasis.h
+++ b/src/2geom/bezier-to-sbasis.h
@@ -1,7 +1,7 @@
 /**
- * \file bezier-to-sbasis.h
- * \brief  \todo brief description
- *
+ * \file
+ * \brief Conversion between Bezier control points and SBasis curves
+ *//*
  * Copyright 2006 Nathan Hurst <njh@mail.csse.monash.edu.au>
  *
  * This library is free software; you can redistribute it and/or
diff --git a/src/2geom/bezier-utils.cpp b/src/2geom/bezier-utils.cpp
index eb317940f..af07db707 100644
--- a/src/2geom/bezier-utils.cpp
+++ b/src/2geom/bezier-utils.cpp
@@ -1,9 +1,5 @@
-#define __SP_BEZIER_UTILS_C__
-
-/** \file
- * Bezier interpolation for inkscape drawing code.
- */
-/*
+/* Bezier interpolation for inkscape drawing code.
+ *
  * Original code published in:
  *   An Algorithm for Automatically Fitting Digitized Curves
  *   by Philip J. Schneider
@@ -52,8 +48,7 @@
 #endif
 
 #include <2geom/bezier-utils.h>
-
-#include <2geom/isnan.h>
+#include <2geom/math-utils.h>
 #include <assert.h>
 
 namespace Geom {
diff --git a/src/2geom/bezier-utils.h b/src/2geom/bezier-utils.h
index 9689db82d..3e56e6e25 100644
--- a/src/2geom/bezier-utils.h
+++ b/src/2geom/bezier-utils.h
@@ -1,10 +1,7 @@
-#ifndef SEEN_GEOM_BEZIER_UTILS_H
-#define SEEN_GEOM_BEZIER_UTILS_H
-
 /**
  * \file
- * \brief  \todo brief description
- *
+ * \brief Bezier fitting algorithms
+ *//*
  * An Algorithm for Automatically Fitting Digitized Curves
  * by Philip J. Schneider
  * from "Graphics Gems", Academic Press, 1990
@@ -41,11 +38,13 @@
  *
  */
 
+#ifndef LIB2GEOM_SEEN_BEZIER_UTILS_H
+#define LIB2GEOM_SEEN_BEZIER_UTILS_H
+
 #include <2geom/point.h>
 
-namespace Geom{
+namespace Geom {
 
-/* Bezier approximation utils */
 Point bezier_pt(unsigned degree, Point const V[], double t);
 
 int bezier_fit_cubic(Point bezier[], Point const data[], int len, double error);
@@ -84,8 +83,9 @@ cubic_bezier_poly_coeff(iterator b, Point *pc) {
 	}
 }
 
-}
-#endif /* !SEEN_GEOM_BEZIER_UTILS_H */
+} // end namespace Geom
+
+#endif // LIB2GEOM_SEEN_BEZIER_UTILS_H
 
 /*
   Local Variables:
diff --git a/src/2geom/bezier.h b/src/2geom/bezier.h
index a7d75da45..48a1dc750 100644
--- a/src/2geom/bezier.h
+++ b/src/2geom/bezier.h
@@ -1,10 +1,13 @@
 /**
- * \file bezier.h
- * \brief \todo brief description
+ * @file
+ * @brief Bezier polynomial
+ *//*
+ * Authors:
+ *   MenTaLguY <mental@rydia.net>
+ *   Michael Sloan <mgsloan@gmail.com>
+ *   Nathan Hurst <njh@njhurst.com>
  *
- * Copyright 2007  MenTaLguY <mental@rydia.net>
- * Copyright 2007  Michael Sloan <mgsloan@gmail.com>
- * Copyright 2007  Nathan Hurst <njh@njhurst.com>
+ * Copyright 2007 Authors
  *
  * This library is free software; you can redistribute it and/or
  * modify it either under the terms of the GNU Lesser General Public
@@ -31,15 +34,15 @@
  *
  */
 
-#ifndef SEEN_BEZIER_H
-#define SEEN_BEZIER_H
+#ifndef LIB2GEOM_SEEN_BEZIER_H
+#define LIB2GEOM_SEEN_BEZIER_H
 
-#include <2geom/coord.h>
 #include <valarray>
-#include <2geom/isnan.h>
+#include <boost/optional.hpp>
+#include <2geom/coord.h>
+#include <2geom/math-utils.h>
 #include <2geom/d2.h>
 #include <2geom/solver.h>
-#include <boost/optional/optional.hpp>
 
 namespace Geom {
 
@@ -280,7 +283,7 @@ public:
     }
     std::vector<double> roots(Interval const ivl) const {
         std::vector<double> solutions;
-        find_bernstein_roots(&const_cast<std::valarray<Coord>&>(c_)[0], order(), solutions, 0, ivl[0], ivl[1]);
+        find_bernstein_roots(&const_cast<std::valarray<Coord>&>(c_)[0], order(), solutions, 0, ivl.min(), ivl.max());
         return solutions;
     }
 };
@@ -407,7 +410,7 @@ inline std::ostream &operator<< (std::ostream &out_file, const Bezier & b) {
 }
 
 }
-#endif //SEEN_BEZIER_H
+#endif // LIB2GEOM_SEEN_BEZIER_H
 
 /*
   Local Variables:
diff --git a/src/2geom/choose.h b/src/2geom/choose.h
index 3fecf1ba2..64ce76f39 100644
--- a/src/2geom/choose.h
+++ b/src/2geom/choose.h
@@ -1,7 +1,7 @@
 /**
- * \file choose.h
- * \brief  \todo brief description
- *
+ * \file
+ * \brief Calculation of binomial cefficients
+ *//*
  * Copyright 2006 Nathan Hurst <njh@mail.csse.monash.edu.au>
  *
  * This library is free software; you can redistribute it and/or
@@ -29,10 +29,13 @@
  *
  */
 
-#ifndef _CHOOSE_H
-#define _CHOOSE_H
+#ifndef LIB2GEOM_SEEN_CHOOSE_H
+#define LIB2GEOM_SEEN_CHOOSE_H
+
 #include <vector>
 
+namespace Geom {
+
 // XXX: Can we keep only the left terms easily?
 // this would more than halve the array
 // row index becomes n2 = n/2, row2 = n2*(n2+1)/2, row = row2*2+(n&1)?n2:0
@@ -121,13 +124,9 @@ class BinomialCoefficient
     container_type coefficients;
 };
 
+} // end namespace Geom
 
-
-
-
-
-
-#endif
+#endif // LIB2GEOM_SEEN_CHOOSE_H
 
 /*
   Local Variables:
diff --git a/src/2geom/circle.h b/src/2geom/circle.h
index ec58e163a..67a638437 100644
--- a/src/2geom/circle.h
+++ b/src/2geom/circle.h
@@ -1,7 +1,7 @@
 /**
  * \file
- * \brief Circle Curve
- *
+ * \brief Circles
+ *//*
  * Authors:
  *      Marco Cecchetti <mrcekets at gmail.com>
  *
@@ -31,18 +31,15 @@
  * the specific language governing rights and limitations.
  */
 
+#ifndef LIB2GEOM_SEEN_CIRCLE_H
+#define LIB2GEOM_SEEN_CIRCLE_H
 
-#ifndef _2GEOM_CIRCLE_H_
-#define _2GEOM_CIRCLE_H_
-
-
+#include <vector>
 #include <2geom/point.h>
 #include <2geom/exception.h>
 #include <2geom/path.h>
-#include <vector>
 
-namespace Geom
-{
+namespace Geom {
 
 class EllipticalArc;
 
@@ -115,13 +112,9 @@ class Circle
     Coord m_ray;
 };
 
-
 } // end namespace Geom
 
-
-
-#endif // _2GEOM_CIRCLE_H_
-
+#endif // LIB2GEOM_SEEN_CIRCLE_H
 
 /*
   Local Variables:
diff --git a/src/2geom/circulator.h b/src/2geom/circulator.h
index 1a70dc4d3..9671ce4a9 100644
--- a/src/2geom/circulator.h
+++ b/src/2geom/circulator.h
@@ -1,7 +1,7 @@
 /**
- * \file circulator.h
- * \brief  \todo brief description
- *
+ * @file circulator.h
+ * @brief Circular iterator adapter
+ *//*
  * Copyright 2006 MenTaLguY <mental@rydia.net>
  *
  * This library is free software; you can redistribute it and/or
@@ -36,6 +36,9 @@
 
 namespace Geom {
 
+/** @brief Circular iterator adapter
+ * This iterator adapter will loop indefinitely over a set of values
+ * from a random access container. */
 template <typename Iterator>
 class Circulator {
 public:
diff --git a/src/2geom/concepts.h b/src/2geom/concepts.h
index a03538d42..c89c3a224 100644
--- a/src/2geom/concepts.h
+++ b/src/2geom/concepts.h
@@ -1,7 +1,7 @@
 /**
  * \file
- * \brief Declares various mathematical concepts, for restriction of template parameters
- *
+ * \brief Template concepts used by 2Geom
+ *//*
  * Copyright 2007 Michael Sloan <mgsloan@gmail.com>
  *
  * This library is free software; you can redistribute it and/or
@@ -26,15 +26,15 @@
  * 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 SEEN_CONCEPTS_H
-#define SEEN_CONCEPTS_H
+#ifndef LIB2GEOM_SEEN_CONCEPTS_H
+#define LIB2GEOM_SEEN_CONCEPTS_H
 
 #include <2geom/sbasis.h>
 #include <2geom/interval.h>
 #include <2geom/point.h>
+#include <2geom/rect.h>
 #include <vector>
 #include <boost/concept_check.hpp>
 #include <2geom/forward.h>
@@ -49,7 +49,7 @@ template <> struct ResultTraits<double> {
   typedef SBasis sb_type;
 };
 
-template <> struct ResultTraits<Point > {
+template <> struct ResultTraits<Point> {
   typedef OptRect bounds_type;
   typedef D2<SBasis> sb_type;
 };
@@ -130,7 +130,7 @@ struct ScalableConcept {
     }
 };
 
-template <class T>
+template <typename T>
 struct AddableConcept {
     T i, j;
     void constraints() {
@@ -139,7 +139,7 @@ struct AddableConcept {
     }
 };
 
-template <class T>
+template <typename T>
 struct MultiplicableConcept {
     T i, j;
     void constraints() {
@@ -147,9 +147,9 @@ struct MultiplicableConcept {
     }
 };
 
-};
+} // end namespace Geom
 
-#endif //SEEN_CONCEPTS_H
+#endif // LIB2GEOM_SEEN_CONCEPTS_H
 
 /*
   Local Variables:
diff --git a/src/2geom/conic_section_clipper_impl.cpp b/src/2geom/conic_section_clipper_impl.cpp
index edfafb11c..33a218a8c 100644
--- a/src/2geom/conic_section_clipper_impl.cpp
+++ b/src/2geom/conic_section_clipper_impl.cpp
@@ -1,6 +1,4 @@
-/**
- * \file
- * \brief Conic section clipping with respect to a rectangle
+/* Conic section clipping with respect to a rectangle
  *
  * Authors:
  *      Marco Cecchetti <mrcekets at gmail>
@@ -31,30 +29,13 @@
  * the specific language governing rights and limitations.
  */
 
-
-
-
 #ifndef CLIP_WITH_CAIRO_SUPPORT
     #include <2geom/conic_section_clipper.h>
 #endif
 
-
-
-
 namespace Geom
 {
 
-struct lex_lesser
-{
-    bool operator() (const Point & P, const Point & Q) const
-    {
-        if (P[X] < Q[X])  return true;
-        if (P[X] == Q[X] && P[Y] < Q[Y]) return true;
-        return false;
-    }
-};
-
-
 /*
  *  Find rectangle-conic crossing points. They are returned in the
  *  "crossing_points" parameter.
@@ -192,7 +173,7 @@ bool CLIPPER_CLASS::intersect (std::vector<Point> & crossing_points) const
               cpts.size())
 
     // remove duplicates
-    std::sort (cpts.begin(), cpts.end(), lex_lesser());
+    std::sort (cpts.begin(), cpts.end(), Point::LexOrder<X>());
     cpts.erase (std::unique (cpts.begin(), cpts.end()), cpts.end());
 
 
diff --git a/src/2geom/conicsec.cpp b/src/2geom/conicsec.cpp
index 2a537a1f0..a7e8e0ad8 100644
--- a/src/2geom/conicsec.cpp
+++ b/src/2geom/conicsec.cpp
@@ -1,7 +1,4 @@
-/**
- * \file
- * \brief Circle Curve
- *
+/*
  * Authors:
  *      Nathan Hurst <njh@njhurst.com
  *
diff --git a/src/2geom/conjugate_gradient.h b/src/2geom/conjugate_gradient.h
index 8ea1b83b4..a34307d4b 100644
--- a/src/2geom/conjugate_gradient.h
+++ b/src/2geom/conjugate_gradient.h
@@ -1,7 +1,7 @@
 /**
- * \file
- * \brief \todo brief description
- *
+ * @file
+ * @brief Routines for solving a system of linear equations using the conjugate gradient method
+ *//*
  * Copyright 2006 Nathan Hurst <njh@mail.csse.monash.edu.au>
  *
  * This library is free software; you can redistribute it and/or
diff --git a/src/2geom/convex-cover.h b/src/2geom/convex-cover.h
index d5e2dee44..e4b5de200 100644
--- a/src/2geom/convex-cover.h
+++ b/src/2geom/convex-cover.h
@@ -1,9 +1,6 @@
-#ifndef GEOM_CONVEX_COVER_H
-#define GEOM_CONVEX_COVER_H
-
 /**
  * \file
- * \brief \todo brief description
+ * \brief Dynamic convex hull structure
  *
  * Copyright 2006 Nathan Hurst <njh@mail.csse.monash.edu.au>
  * Copyright 2006 Michael G. Sloan <mgsloan@gmail.com>
@@ -33,15 +30,18 @@
  *
  */
 
-/** A convex cover is a sequence of convex polygons that completely cover the path.  For now a
- * convex hull class is included here (the convex-hull header is wrong)
- */
+#ifndef GEOM_CONVEX_COVER_H
+#define GEOM_CONVEX_COVER_H
 
 #include <2geom/point.h>
 #include <vector>
 
 namespace Geom{
 
+/* A convex cover is a sequence of convex polygons that completely cover the path.  For now a
+ * convex hull class is included here (the convex-hull header is wrong)
+ */
+
 /** ConvexHull
  * A convexhull is a convex region - every point between two points in the convex hull is also in
  * the convex hull.  It is defined by a set of points travelling in a clockwise direction.  We require the first point to be top most, and of the topmost, leftmost.
diff --git a/src/2geom/coord.h b/src/2geom/coord.h
index 9c42f6bfc..c7bbcdcd4 100644
--- a/src/2geom/coord.h
+++ b/src/2geom/coord.h
@@ -1,7 +1,7 @@
 /**
  *  \file
  *  \brief Defines the Coord "real" type with sufficient precision for coordinates.
- *
+ *//*
  * Copyright 2006 Nathan Hurst <njh@mail.csse.monash.edu.au>
  *
  * This library is free software; you can redistribute it and/or
@@ -29,15 +29,16 @@
  *
  */
 
-#ifndef SEEN_Geom_COORD_H
-#define SEEN_Geom_COORD_H
+#ifndef LIB2GEOM_SEEN_COORD_H
+#define LIB2GEOM_SEEN_COORD_H
 
 #include <cmath>
 #include <limits>
+#include <2geom/forward.h>
 
 namespace Geom {
 
-/** @brief Axis enum (X or Y). */
+/** @brief 2D axis enumeration (X or Y). */
 enum Dim2 { X=0, Y=1 };
 
 /**
@@ -48,6 +49,7 @@ enum Dim2 { X=0, Y=1 };
  * differences of on-canvas points.
  */
 typedef double Coord;
+typedef int IntCoord;
 
 const Coord EPSILON = 1e-5; //1e-18;
 
@@ -57,13 +59,36 @@ inline Coord infinity() {  return std::numeric_limits<Coord>::infinity();  }
 inline bool are_near(Coord a, Coord b, double eps=EPSILON) { return a-b <= eps && a-b >= -eps; }
 inline bool rel_error_bound(Coord a, Coord b, double eps=EPSILON) { return a <= eps*b && a >= -eps*b; }
 
+template <typename C>
+struct CoordTraits {};
 
-typedef long IntCoord;
+template<>
+struct CoordTraits<IntCoord> {
+    typedef IntPoint PointType;
+    typedef IntInterval IntervalType;
+    typedef OptIntInterval OptIntervalType;
+    typedef IntRect RectType;
+    typedef OptIntRect OptRectType;
+    inline static bool contains(IntCoord low, IntCoord high, IntCoord testlow, IntCoord testhigh) {
+        return low <= testlow && testhigh < high;
+    }
+};
 
-} /* namespace Geom */
+template<>
+struct CoordTraits<Coord> {
+    typedef Point PointType;
+    typedef Interval IntervalType;
+    typedef OptInterval OptIntervalType;
+    typedef Rect RectType;
+    typedef OptRect OptRectType;
+    inline static bool contains(Coord low, Coord high, Coord testlow, Coord testhigh) {
+        return low <= testlow && testhigh <= high;
+    }
+};
 
+} // end namespace Geom
 
-#endif /* !SEEN_Geom_COORD_H */
+#endif // LIB2GEOM_SEEN_COORD_H
 
 /*
   Local Variables:
diff --git a/src/2geom/crossing.h b/src/2geom/crossing.h
index 62e447450..75c75fc24 100644
--- a/src/2geom/crossing.h
+++ b/src/2geom/crossing.h
@@ -1,10 +1,10 @@
 /**
- * \file
- * \brief  \todo brief description
- *
+ * @file
+ * @brief Structure representing the intersection of two curves
+ *//*
  * Authors:
- *      Michael Sloane <?@?.?>
- *      Marco <?@?.?>
+ *   Michael Sloan <mgsloan@gmail.com>
+ *   Marco Cecchetti <mrcekets at gmail.com>
  * 
  * Copyright 2006-2008  authors
  *
diff --git a/src/2geom/curve.cpp b/src/2geom/curve.cpp
index 49e011a8b..fe9d607d8 100644
--- a/src/2geom/curve.cpp
+++ b/src/2geom/curve.cpp
@@ -1,8 +1,5 @@
-/**
- * \file
- * \brief Abstract curve type - implementation of default methods
+/* Abstract curve type - implementation of default methods
  *
- *//*
  * Authors:
  *   MenTaLguY <mental@rydia.net>
  *   Marco Cecchetti <mrcekets at gmail.com>
diff --git a/src/2geom/curves.h b/src/2geom/curves.h
index 64cf3d4fb..319b1924d 100644
--- a/src/2geom/curves.h
+++ b/src/2geom/curves.h
@@ -32,9 +32,8 @@
  * the specific language governing rights and limitations.
  */
 
-#ifndef _2GEOM_CURVES_H_
-#define _2GEOM_CURVES_H_
-
+#ifndef LIB2GEOM_SEEN_CURVES_H
+#define LIB2GEOM_SEEN_CURVES_H
 
 #include <2geom/curve.h>
 #include <2geom/sbasis-curve.h>
@@ -43,7 +42,7 @@
 #include <2geom/elliptical-arc.h>
 #include <2geom/svg-elliptical-arc.h>
 
-#endif // _2GEOM_CURVES_H_
+#endif // LIB2GEOM_SEEN_CURVES_H
 
 /*
   Local Variables:
diff --git a/src/2geom/d2-sbasis.h b/src/2geom/d2-sbasis.h
index bd6c35805..95c0da4ed 100644
--- a/src/2geom/d2-sbasis.h
+++ b/src/2geom/d2-sbasis.h
@@ -1,11 +1,10 @@
 /**
  * \file
- * \brief  Do not include this file \todo brief description
+ * \brief  Do not include this file
  *
  * We don't actually want anyone to
- * include this, other than D2.h.  If somone else tries, D2
- * won't be defined.  If it is, this will already be included.
- *
+ * include this, other than D2.h.
+ *//*
  * Authors:
  *      ? <?@?.?>
  * 
diff --git a/src/2geom/d2.h b/src/2geom/d2.h
index 3e4de430e..73330295b 100644
--- a/src/2geom/d2.h
+++ b/src/2geom/d2.h
@@ -1,7 +1,7 @@
 /**
  * \file
  * \brief   Lifts one dimensional objects into 2d 
- *
+ *//*
  * Copyright 2007 Michael Sloan <mgsloan@gmail.com>
  *
  * This library is free software; you can redistribute it and/or
@@ -426,7 +426,6 @@ inline std::ostream &operator<< (std::ostream &out_file, const Geom::D2<T> &in_d
 
 } //end namespace Geom
 
-#include <2geom/rect.h>
 #include <2geom/d2-sbasis.h>
 
 namespace Geom{
diff --git a/src/2geom/forward.h b/src/2geom/forward.h
index 399344dda..b1cad6f1f 100644
--- a/src/2geom/forward.h
+++ b/src/2geom/forward.h
@@ -41,11 +41,23 @@ namespace Geom {
 
 // basic types
 typedef double Coord;
+typedef int IntCoord;
 class Point;
-class Interval;
-class OptInterval;
+class IntPoint;
 class Line;
 class Ray;
+template <typename> class GenericInterval;
+template <typename> class GenericOptInterval;
+class Interval;
+typedef GenericOptInterval<Coord> OptInterval;
+typedef GenericInterval<IntCoord> IntInterval;
+typedef GenericOptInterval<IntCoord> OptIntInterval;
+template <typename> class GenericRect;
+template <typename> class GenericOptRect;
+class Rect;
+typedef GenericOptRect<Coord> OptRect;
+typedef GenericRect<IntCoord> IntRect;
+typedef GenericOptRect<IntCoord> OptIntRect;
 
 // fragments
 class Linear;
@@ -90,9 +102,6 @@ class VShear;
 template <typename> class D2;
 template <typename> class Piecewise;
 
-typedef D2<Interval> Rect;
-class OptRect;
-
 class Shape;
 class Region;
 class Hat;
diff --git a/src/2geom/generic-interval.h b/src/2geom/generic-interval.h
new file mode 100644
index 000000000..d719c16c8
--- /dev/null
+++ b/src/2geom/generic-interval.h
@@ -0,0 +1,342 @@
+/**
+ *  @file
+ *  @brief Closed interval of generic values
+ *//*
+ * Copyright 2011 Krzysztof Kosiński <tweenk.pl@gmail.com>
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, 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_GENERIC_INTERVAL_H
+#define LIB2GEOM_SEEN_GENERIC_INTERVAL_H
+
+#include <cassert>
+#include <boost/none.hpp>
+#include <boost/optional.hpp>
+#include <boost/operators.hpp>
+
+namespace Geom {
+
+template <typename C>
+class GenericOptInterval;
+
+/**
+ * @brief A range of numbers which is never empty.
+ * @ingroup Primitives
+ */
+template <typename C>
+class GenericInterval
+    : boost::equality_comparable< GenericInterval<C>
+    , boost::additive< GenericInterval<C>
+    , boost::additive< GenericInterval<C>, C
+    , boost::orable< GenericInterval<C>
+      > > > >
+{
+    typedef GenericInterval<C> Self;
+protected:
+    C _b[2];
+public:
+    /// @name Create intervals.
+    /// @{
+    /** @brief Create an interval that contains only zero. */
+    GenericInterval() { _b[0] = 0;  _b[1] = 0; }
+    /** @brief Create an interval that contains a single point. */
+    explicit GenericInterval(C u) { _b[0] = _b[1] = u; }
+    /** @brief Create an interval that contains all points between @c u and @c v. */
+    GenericInterval(C u, C v) {
+        if (u <= v) {
+            _b[0] = u; _b[1] = v;
+        } else {
+            _b[0] = v; _b[1] = u;
+        }
+    }
+
+    /** @brief Create an interval containing a range of values.
+     * The resulting interval will contain all values from the given range.
+     * The return type of iterators must be convertible to C. The given range
+     * must not be empty. For potentially empty ranges, see GenericOptInterval.
+     * @param start Beginning of the range
+     * @param end   End of the range
+     * @return Interval that contains all values from [start, end). */
+    template <typename InputIterator>
+    static Self from_range(InputIterator start, InputIterator end) {
+        assert(start != end);
+        Self result(*start++);
+        for (; start != end; ++start) result.expandTo(*start);
+        return result;
+    }
+    /** @brief Create an interval from a C-style array of values it should contain. */
+    static Self from_array(C const *c, unsigned n) {
+        Self result = from_range(c, c+n);
+        return result;
+    }
+    /// @}
+
+    /// @name Inspect endpoints.
+    /// @{
+    C min() const { return _b[0]; }
+    C max() const { return _b[1]; }
+    C extent() const { return max() - min(); }
+    C middle() const { return (max() + min()) * 0.5; }
+    bool isSingular() const { return min() == max(); }
+    /// @}
+
+    /// @name Test coordinates and other intervals for inclusion.
+    /// @{
+    /** @brief Check whether the interval includes this number. */
+    bool contains(C val) const {
+        return CoordTraits<C>::contains(min(), max(), val, val);
+    }
+    /** @brief Check whether the interval includes the given interval. */
+    bool contains(Self const &val) const {
+        return CoordTraits<C>::contains(min(), max(), val.min(), val.max());
+    }
+    /** @brief Check whether the intervals have any common elements. */
+    bool intersects(Self const &val) const {
+        return contains(val.min()) || contains(val.max()) || val.contains(*this);
+    }
+    /// @}
+
+    /// @name Modify the interval.
+    /// @{
+    //TODO: NaN handleage for the next two?
+    /** @brief Set the lower boundary of the interval.
+     * When the given number is larger than the interval's largest element,
+     * it will be reduced to the single number @c val. */
+    void setMin(C val) {
+        if(val > _b[1]) {
+            _b[0] = _b[1] = val;
+        } else {
+            _b[0] = val;
+        }
+    }
+    /** @brief Set the upper boundary of the interval.
+     * When the given number is smaller than the interval's smallest element,
+     * it will be reduced to the single number @c val. */
+    void setMax(C val) {
+        if(val < _b[0]) {
+            _b[1] = _b[0] = val;
+        } else {
+            _b[1] = val;
+        }
+    }
+    /** @brief Extend the interval to include the given number. */
+    void expandTo(C val) {
+       if(val < _b[0]) _b[0] = val;
+       if(val > _b[1]) _b[1] = val;  //no else, as we want to handle NaN
+    }
+    /** @brief Expand or shrink the interval in both directions by the given amount.
+     * After this method, the interval's length (extent) will be increased by
+     * <code>amount * 2</code>. Negative values can be given; they will shrink the interval.
+     * Shrinking by a value larger than half the interval's length will create a degenerate
+     * interval containing only the midpoint of the original. */
+    void expandBy(C amount) {
+        _b[0] -= amount;
+        _b[1] += amount;
+        if (_b[0] > _b[1]) {
+            C halfway = (_b[0]+_b[1])/2;
+            _b[0] = _b[1] = halfway;
+        }
+    }
+    /** @brief Union the interval with another one.
+     * The resulting interval will contain all points of both intervals.
+     * It might also contain some points which didn't belong to either - this happens
+     * when the intervals did not have any common elements. */
+    void unionWith(Self const &a) {
+        if(a._b[0] < _b[0]) _b[0] = a._b[0];
+        if(a._b[1] > _b[1]) _b[1] = a._b[1];
+    }
+    /// @}
+
+    /// @name Operators
+    /// @{
+    //IMPL: OffsetableConcept
+    //TODO: rename output_type to something else in the concept
+    typedef C output_type;
+    /** @brief Offset the interval by a specified amount */
+    Self &operator+=(C amnt) {
+        _b[0] += amnt; _b[1] += amnt;
+        return *this;
+    }
+    /** @brief Offset the interval by the negation of the specified amount */
+    Self &operator-=(C amnt) {
+        _b[0] -= amnt; _b[1] -= amnt;
+        return *this;
+    }
+    
+    /** @brief Return an interval mirrored about 0 */
+    Self operator-() const { Self r(-_b[1], -_b[0]); return r; }
+    // IMPL: AddableConcept
+    /** @brief Add two intervals.
+     * Sum is defined as the set of points that can be obtained by adding any two values
+     * from both operands: \f$S = \{x \in A, y \in B: x + y\}\f$ */
+    Self &operator+=(Self const &o) {
+        _b[0] += o._b[0];
+        _b[1] += o._b[1];
+        return *this;
+    }
+    /** @brief Subtract two intervals.
+     * Difference is defined as the set of points that can be obtained by subtracting
+     * any value from the second operand from any value from the first operand:
+     * \f$S = \{x \in A, y \in B: x - y\}\f$ */
+    Self &operator-=(Self const &o) {
+        // equal to *this += -o
+        _b[0] -= o._b[1];
+        _b[1] -= o._b[0];
+        return *this;
+    }
+    /** @brief Union two intervals.
+     * Note that the intersection-and-assignment operator is not defined,
+     * because the result of an intersection can be empty, while Interval cannot. */
+    Self &operator|=(Self const &o) {
+        unionWith(o);
+        return *this;
+    }
+    /** @brief Test for interval equality. */
+    bool operator==(Self const &other) const {
+        return min() == other.min() && max() == other.max();
+    }
+    /// @}
+};
+
+/** @brief Union two intervals
+ * @relates GenericInterval */
+template <typename C>
+inline GenericInterval<C> unify(GenericInterval<C> const &a, GenericInterval<C> const &b) {
+    return a | b;
+}
+
+/**
+ * @brief A range of numbers that can be empty.
+ * @ingroup Primitives
+ */
+template <typename C>
+class GenericOptInterval
+    : public boost::optional<typename CoordTraits<C>::IntervalType>
+    , boost::orable< GenericOptInterval<C>, typename CoordTraits<C>::OptIntervalType
+    , boost::andable< GenericOptInterval<C>, typename CoordTraits<C>::OptIntervalType
+      > >
+{
+    typedef typename CoordTraits<C>::IntervalType CInterval;
+    typedef typename CoordTraits<C>::OptIntervalType OptCInterval;
+    typedef boost::optional<CInterval> Base;
+public:
+    /// @name Create optionally empty intervals of integers.
+    /// @{
+    /** @brief Create an empty interval. */
+    GenericOptInterval() : Base() {}
+    /** @brief Wrap an existing interval. */
+    GenericOptInterval(GenericInterval<C> const &a) : Base(CInterval(a)) {}
+    /** @brief Create an interval containing a single point. */
+    GenericOptInterval(C u) : Base(CInterval(u)) {}
+    /** @brief Create an interval containing a range of numbers. */
+    GenericOptInterval(C u, C v) : Base(CInterval(u,v)) {}
+
+    /** @brief Create a possibly empty interval containing a range of values.
+     * The resulting interval will contain all values from the given range.
+     * The return type of iterators must be convertible to C. The given range
+     * may be empty.
+     * @param start Beginning of the range
+     * @param end   End of the range
+     * @return Interval that contains all values from [start, end), or nothing if the range
+     *         is empty. */
+    template <typename InputIterator>
+    static GenericOptInterval<C> from_range(InputIterator start, InputIterator end) {
+        if (start == end) {
+            GenericOptInterval<C> ret;
+            return ret;
+        }
+        GenericOptInterval<C> ret(CInterval::from_range(start, end));
+        return ret;
+    }
+    /// @}
+
+    /** @brief Check whether this interval is empty. */
+    bool isEmpty() { return !*this; };
+
+    /** @brief Union with another interval, gracefully handling empty ones. */
+    void unionWith(GenericOptInterval<C> const &a) {
+        if (*this) { // check that we are not empty
+            (*this)->unionWith(*a);
+        } else if (a) {
+            *this = *a;
+        }
+    }
+    void intersectWith(GenericOptInterval<C> const &o) {
+        if (o && *this) {
+            if (!*this) return;
+            C u = std::max((*this)->min(), o->min());
+            C v = std::min((*this)->max(), o->max());
+            if (u <= v) {
+                *this = CInterval(u, v);
+                return;
+            }
+        }
+        (*static_cast<Base*>(this)) = boost::none;
+    }
+    GenericOptInterval<C> &operator|=(OptCInterval const &o) {
+        unionWith(o);
+        return *this;
+    }
+    GenericOptInterval<C> &operator&=(OptCInterval const &o) {
+        intersectWith(o);
+        return *this;
+    }
+};
+
+/** @brief Intersect two intervals and return a possibly empty range of numbers
+ * @relates GenericOptInterval */
+template <typename C>
+inline GenericOptInterval<C> intersect(GenericInterval<C> const &a, GenericInterval<C> const &b) {
+    return GenericOptInterval<C>(a) & GenericOptInterval<C>(b);
+}
+/** @brief Intersect two intervals and return a possibly empty range of numbers
+ * @relates GenericOptInterval */
+template <typename C>
+inline GenericOptInterval<C> operator&(GenericInterval<C> const &a, GenericInterval<C> const &b) {
+    return GenericOptInterval<C>(a) & GenericOptInterval<C>(b);
+}
+
+#ifdef _GLIBCXX_IOSTREAM
+template <typename C>
+inline std::ostream &operator<< (std::ostream &os, 
+                                 Geom::GenericInterval<C> const &I) {
+    os << "Interval("<<I.min() << ", "<<I.max() << ")";
+    return os;
+}
+#endif
+
+} // namespace Geom
+#endif // !LIB2GEOM_SEEN_GENERIC_INTERVAL_H
+
+/*
+  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/generic-rect.h b/src/2geom/generic-rect.h
new file mode 100644
index 000000000..cc0d7d42e
--- /dev/null
+++ b/src/2geom/generic-rect.h
@@ -0,0 +1,363 @@
+/**
+ * \file
+ * \brief Axis-aligned rectangle
+ *//*
+ * Authors:
+ *   Michael Sloan <mgsloan@gmail.com>
+ *   Krzysztof Kosiński <tweenk.pl@gmail.com>
+ * Copyright 2007-2011 Authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, output to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ *
+ * Authors of original rect class:
+ *   Lauris Kaplinski <lauris@kaplinski.com>
+ *   Nathan Hurst <njh@mail.csse.monash.edu.au>
+ *   bulia byak <buliabyak@users.sf.net>
+ *   MenTaLguY <mental@rydia.net>
+ */
+
+#ifndef LIB2GEOM_SEEN_GENERIC_RECT_H
+#define LIB2GEOM_SEEN_GENERIC_RECT_H
+
+#include <boost/optional.hpp>
+
+namespace Geom {
+
+template <typename C>
+class GenericOptRect;
+
+/**
+ * @brief Axis aligned, non-empty, generic rectangle.
+ * @ingroup Primitives
+ */
+template <typename C>
+class GenericRect
+    : boost::additive< GenericRect<C>, typename CoordTraits<C>::PointType
+    , boost::equality_comparable< GenericRect<C>
+    , boost::orable< GenericRect<C>
+    , boost::orable< GenericRect<C>, typename CoordTraits<C>::OptRectType
+      > > > >
+{
+    typedef typename CoordTraits<C>::IntervalType CInterval;
+    typedef typename CoordTraits<C>::PointType CPoint;
+    typedef typename CoordTraits<C>::RectType CRect;
+    typedef typename CoordTraits<C>::OptRectType OptCRect;
+protected:
+    CInterval f[2];
+public:
+    /// @name Create rectangles.
+    /// @{
+    /** @brief Create a rectangle that contains only the point at (0,0). */
+    GenericRect() { f[X] = f[Y] = Interval(); }
+    /** @brief Create a rectangle from X and Y intervals. */
+    GenericRect(CInterval const &a, CInterval const &b) {
+        f[X] = a;
+        f[Y] = b;
+    }
+    /** @brief Create a rectangle from two points. */
+    GenericRect(CPoint const &a, CPoint const &b) {
+        f[X] = Interval(a[X], b[X]);
+        f[Y] = Interval(a[Y], b[Y]);
+    }
+    /** @brief Create a rectangle from a range of points.
+     * The resulting rectangle will contain all ponts from the range.
+     * The return type of iterators must be convertible to Point.
+     * The range must not be empty. For possibly empty ranges, see OptRect.
+     * @param start Beginning of the range
+     * @param end   End of the range
+     * @return Rectangle that contains all points from [start, end). */
+    template <typename InputIterator>
+    static GenericRect<C> from_range(InputIterator start, InputIterator end) {
+        assert(start != end);
+        CPoint p1 = *start++;
+        GenericRect<C> result(p1, p1);
+        for (; start != end; ++start) {
+            result.expandTo(*start);
+        }
+        return result;
+    }
+    /** @brief Create a rectangle from a C-style array of points it should contain. */
+    static GenericRect<C> from_array(CPoint const *c, unsigned n) {
+        GenericRect<C> result = GenericRect<C>::from_range(c, c+n);
+        return result;
+    }
+    /// @}
+
+    /// @name Inspect dimensions.
+    /// @{
+    CInterval &operator[](unsigned i)              { return f[i]; }
+    CInterval const &operator[](unsigned i) const { return f[i]; }
+
+    CPoint min() const { return CPoint(f[X].min(), f[Y].min()); }
+    CPoint max() const { return CPoint(f[X].max(), f[Y].max()); }
+    /** @brief Return the n-th corner of the rectangle.
+     * If the Y axis grows upwards, this returns corners in clockwise order
+     * starting from the lower left. If Y grows downwards, it returns the corners
+     * in counter-clockwise order starting from the upper left. */
+    CPoint corner(unsigned i) const {
+        switch(i % 4) {
+            case 0:  return CPoint(f[X].min(), f[Y].min());
+            case 1:  return CPoint(f[X].max(), f[Y].min());
+            case 2:  return CPoint(f[X].max(), f[Y].max());
+            default: return CPoint(f[X].min(), f[Y].max());
+        }
+    }
+        
+    //We should probably remove these - they're coord sys gnostic
+    /** @brief Return top coordinate of the rectangle (+Y is downwards). */
+    C top() const { return f[Y].min(); }
+    /** @brief Return bottom coordinate of the rectangle (+Y is downwards). */
+    C bottom() const { return f[Y].max(); }
+    /** @brief Return leftmost coordinate of the rectangle (+X is to the right). */
+    C left() const { return f[X].min(); }
+    /** @brief Return rightmost coordinate of the rectangle (+X is to the right). */
+    C right() const { return f[X].max(); }
+
+    /** @brief Get the horizontal extent of the rectangle. */
+    C width() const { return f[X].extent(); }
+    /** @brief Get the vertical extent of the rectangle. */
+    C height() const { return f[Y].extent(); }
+
+    /** @brief Get rectangle's width and height as a point.
+     * @return Point with X coordinate corresponding to the width and the Y coordinate
+     *         corresponding to the height of the rectangle. */
+    CPoint dimensions() const { return CPoint(f[X].extent(), f[Y].extent()); }
+    /** @brief Get the point in the geometric center of the rectangle. */
+    CPoint midpoint() const { return CPoint(f[X].middle(), f[Y].middle()); }
+
+    /** @brief Compute rectangle's area. */
+    C area() const { return f[X].extent() * f[Y].extent(); }
+    /** @brief Check whether the rectangle has zero area. */
+    bool hasZeroArea() const { return (area() == 0); }
+
+    /** @brief Get the larger extent (width or height) of the rectangle. */
+    C maxExtent() const { return std::max(f[X].extent(), f[Y].extent()); }
+    /** @brief Get the smaller extent (width or height) of the rectangle. */
+    C minExtent() const { return std::min(f[X].extent(), f[Y].extent()); }
+    /// @}
+
+    /// @name Test other rectangles and points for inclusion.
+    /// @{
+    /** @brief Check whether the rectangles have any common points. */
+    bool intersects(GenericRect<C> const &r) const { 
+        return f[X].intersects(r[X]) && f[Y].intersects(r[Y]);
+    }
+    /** @brief Check whether the rectangle includes all points in the given rectangle. */
+    bool contains(GenericRect<C> const &r) const { 
+        return f[X].contains(r[X]) && f[Y].contains(r[Y]);
+    }
+
+    /** @brief Check whether the rectangles have any common points.
+     * A non-empty rectangle will not intersect empty rectangles. */
+    inline bool intersects(OptCRect const &r) const;
+    /** @brief Check whether the rectangle includes all points in the given rectangle.
+     * A non-empty rectangle will contain any empty rectangle. */
+    inline bool contains(OptCRect const &r) const;
+
+    /** @brief Check whether the given point is within the rectangle. */
+    bool contains(CPoint const &p) const {
+        return f[X].contains(p[X]) && f[Y].contains(p[Y]);
+    }
+    /// @}
+
+    /// @name Modify the rectangle.
+    /// @{
+    /** @brief Enlarge the rectangle to contain the given point. */
+    void expandTo(CPoint const &p)        { 
+        f[X].expandTo(p[X]);  f[Y].expandTo(p[Y]);
+    }
+    /** @brief Enlarge the rectangle to contain the given rectangle. */
+    void unionWith(GenericRect<C> const &b) { 
+        f[X].unionWith(b[X]); f[Y].unionWith(b[Y]);
+    }
+    /** @brief Enlarge the rectangle to contain the given rectangle.
+     * Unioning with an empty rectangle results in no changes. */
+    void unionWith(OptCRect const &b);
+
+    /** @brief Expand the rectangle in both directions by the specified amount.
+     * Note that this is different from scaling. Negative values wil shrink the
+     * rectangle. If <code>-amount</code> is larger than
+     * half of the width, the X interval will contain only the X coordinate
+     * of the midpoint; same for height. */
+    void expandBy(C amount) {
+        f[X].expandBy(amount);  f[Y].expandBy(amount);
+    }
+    /** @brief Expand the rectangle by the coordinates of the given point.
+     * This will expand the width by the X coordinate of the point in both directions
+     * and the height by Y coordinate of the point. Negative coordinate values will
+     * shrink the rectangle. If <code>-p[X]</code> is larger than half of the width,
+     * the X interval will contain only the X coordinate of the midpoint; same for height. */
+    void expandBy(CPoint const &p) { 
+        f[X].expandBy(p[X]);  f[Y].expandBy(p[Y]);
+    }
+    /// @}
+
+    /// @name Operators
+    /// @{
+    /** @brief Offset the rectangle by a vector. */
+    GenericRect<C> &operator+=(CPoint const &p) {
+        f[X] += p[X];
+        f[Y] += p[Y];
+        return *this;
+    }
+    /** @brief Offset the rectangle by the negation of a vector. */
+    GenericRect<C> &operator-=(CPoint const &p) {
+        f[X] -= p[X];
+        f[Y] -= p[Y];
+        return *this;
+    }
+    /** @brief Union two rectangles. */
+    GenericRect<C> &operator|=(GenericRect<C> const &o) {
+        unionWith(o);
+        return *this;
+    }
+    GenericRect<C> &operator|=(OptCRect const &o) {
+        unionWith(o);
+        return *this;
+    }
+    /** @brief Test for equality of rectangles. */
+    bool operator==(GenericRect<C> const &o) const { return f[X] == o[X] && f[Y] == o[Y]; }
+    /// @}
+};
+
+/**
+ * @brief Axis-aligned generic rectangle that can be empty.
+ * @ingroup Primitives
+ */
+template <typename C>
+class GenericOptRect
+    : public boost::optional<typename CoordTraits<C>::RectType>
+    , boost::orable< GenericOptRect<C>
+    , boost::andable< GenericOptRect<C>
+    , boost::andable< GenericOptRect<C>, typename CoordTraits<C>::RectType
+      > > >
+{
+    typedef typename CoordTraits<C>::IntervalType CInterval;
+    typedef typename CoordTraits<C>::OptIntervalType OptCInterval;
+    typedef typename CoordTraits<C>::PointType CPoint;
+    typedef typename CoordTraits<C>::RectType CRect;
+    typedef typename CoordTraits<C>::OptRectType OptCRect;
+    typedef boost::optional<CRect> Base;
+public:
+    GenericOptRect() : Base() {}
+    GenericOptRect(GenericRect<C> const &a) : Base(CRect(a)) {}
+    GenericOptRect(CPoint const &a, CPoint const &b) : Base(CRect(a, b)) {}
+    /**
+     * Creates an empty OptRect when one of the argument intervals is empty.
+     */
+    GenericOptRect(OptCInterval const &x_int, OptCInterval const &y_int) {
+        if (x_int && y_int) {
+            *this = CRect(*x_int, *y_int);
+        }
+        // else, stay empty.
+    }
+
+    /** @brief Check for emptiness. */
+    inline bool isEmpty() const { return !*this; };
+
+    bool intersects(CRect const &r) const { return r.intersects(*this); }
+    bool contains(CRect const &r) const { return *this && (*this)->contains(r); }
+
+    bool intersects(OptCRect const &r) const { return *this && (*this)->intersects(r); }
+    bool contains(OptCRect const &r) const { return *this && (*this)->contains(r); }
+
+    bool contains(CPoint const &p) const { return *this && (*this)->contains(p); }
+
+    void unionWith(CRect const &b) {
+        if (*this) {
+            (*this)->unionWith(b);
+        } else {
+            *this = b;
+        }
+    }
+    void unionWith(OptCRect const &b) {
+        if (b) unionWith(*b);
+    }
+    void intersectWith(CRect const &b) {
+        if (!*this) return;
+        OptCInterval x = (**this)[X] & b[X], y = (**this)[Y] & b[Y];
+        if (x && y) {
+            *this = CRect(*x, *y);
+        } else {
+            *(static_cast<Base*>(this)) = boost::none;
+        }
+    }
+    void intersectWith(OptCRect const &b) {
+        if (b) {
+            intersectWith(*b);
+        } else {
+            *(static_cast<Base*>(this)) = boost::none;
+        }
+    }
+    GenericOptRect<C> &operator|=(OptCRect const &b) {
+        unionWith(b);
+        return *this;
+    }
+    GenericOptRect<C> &operator&=(CRect const &b) {
+        intersectWith(b);
+        return *this;
+    }
+    GenericOptRect<C> &operator&=(OptCRect const &b) {
+        intersectWith(b);
+        return *this;
+    }
+};
+
+template <typename C>
+inline void GenericRect<C>::unionWith(OptCRect const &b) { 
+    if (b) {
+        unionWith(*b);
+    }
+}
+template <typename C>
+inline bool GenericRect<C>::intersects(OptCRect const &r) const {
+    return r && intersects(*r);
+}
+template <typename C>
+inline bool GenericRect<C>::contains(OptCRect const &r) const {
+    return !r || contains(*r);
+}
+
+#ifdef _GLIBCXX_IOSTREAM
+template <typename C>
+inline std::ostream &operator<<(std::ostream &out, GenericRect<C> const &r) {
+    out << "X: " << r[X] << "  Y: " << r[Y];
+    return out;
+}
+#endif
+
+} // end namespace Geom
+
+#endif // LIB2GEOM_SEEN_RECT_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/hvlinesegment.h b/src/2geom/hvlinesegment.h
index 9419be8f6..05252468e 100644
--- a/src/2geom/hvlinesegment.h
+++ b/src/2geom/hvlinesegment.h
@@ -1,8 +1,11 @@
 /**
  * \file
- * \brief  Horizontal and Vertical Line Segment
- *
- * Copyright 2008  Marco Cecchetti <mrcekets at gmail.com>
+ * \brief  Horizontal and vertical line segment
+ *//*
+ * Authors:
+ *   Marco Cecchetti <mrcekets at gmail.com>
+ *   Krzysztof Kosiński <tweenk.pl@gmail.com>
+ * Copyright 2008-2011 Authors
  *
  * This library is free software; you can redistribute it and/or
  * modify it either under the terms of the GNU Lesser General Public
@@ -28,14 +31,11 @@
  * the specific language governing rights and limitations.
  */
 
-
-#ifndef _2GEOM_HVLINESEGMENT_H_
-#define _2GEOM_HVLINESEGMENT_H_
-
+#ifndef LIB2GEOM_SEEN_HVLINESEGMENT_H
+#define LIB2GEOM_SEEN_HVLINESEGMENT_H
 
 #include <2geom/bezier-curve.h>
 
-
 namespace Geom
 {
 
@@ -44,6 +44,7 @@ class AxisLineSegment : public LineSegment
 {
 public:
     static const Dim2 other_axis = static_cast<Dim2>((axis + 1) % 2);
+#ifndef DOXYGEN_SHOULD_SKIP_THIS
     virtual void setInitial(Point const &p) {
         Point f = finalPoint();
         f[axis] = p[axis];
@@ -106,10 +107,6 @@ public:
         if (d != axis) return initialPoint()[other_axis];
         return initialPoint()[axis] + t * (finalPoint()[axis] - initialPoint()[axis]);
     }
-
-    /**
-    *  The size of the returned vector equals n+1.
-    */
     virtual std::vector<Point> pointAndDerivatives(Coord t, unsigned n) const {
         std::vector<Point> result;
         result.push_back(pointAt(t));
@@ -124,6 +121,7 @@ public:
         }
         return result;
     }
+#endif
 protected:
     AxisLineSegment(Point const &p0, Point const &p1) : LineSegment(p0, p1) {}
     AxisLineSegment() {}
@@ -171,6 +169,7 @@ public:
         return result;
     }
 
+#ifndef DOXYGEN_SHOULD_SKIP_THIS
     virtual Curve* duplicate() const { return new HLineSegment(*this); }
     virtual Curve *portion(Coord f, Coord t) const {
         Point ip = pointAt(f);
@@ -195,6 +194,7 @@ public:
         Coord x = finalPoint()[X] - initialPoint()[X];
         return new HLineSegment(x, x, 0);
     }
+#endif
 };  // end class HLineSegment
 
 
@@ -241,6 +241,7 @@ public:
         return result;
     }
 
+#ifndef DOXYGEN_SHOULD_SKIP_THIS
     virtual Curve *duplicate() const { return new VLineSegment(*this); }
     virtual Curve *portion(Coord f, Coord t) const {
         Point ip = pointAt(f);
@@ -264,13 +265,12 @@ public:
         Coord y = finalPoint()[Y] - initialPoint()[Y];
         return new VLineSegment(0, y, y);
     }
+#endif
 }; // end class VLineSegment
 
 }  // end namespace Geom
 
-
-#endif // _2GEOM_HVLINESEGMENT_H_
-
+#endif // LIB2GEOM_SEEN_HVLINESEGMENT_H
 
 /*
   Local Variables:
diff --git a/src/2geom/int-interval.h b/src/2geom/int-interval.h
new file mode 100644
index 000000000..0faf48d80
--- /dev/null
+++ b/src/2geom/int-interval.h
@@ -0,0 +1,63 @@
+/**
+ *  \file
+ *  \brief Closed interval of integer values
+ *//*
+ * Copyright 2011 Krzysztof Kosiński <tweenk.pl@gmail.com>
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, 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_INT_INTERVAL_H
+#define LIB2GEOM_SEEN_INT_INTERVAL_H
+
+#include <2geom/coord.h>
+#include <2geom/generic-interval.h>
+
+namespace Geom {
+
+/**
+ * @brief Range of integers that is never empty.
+ * @ingroup Primitives
+ */
+typedef GenericInterval<IntCoord> IntInterval;
+
+/**
+ * @brief Range of integers that can be empty.
+ * @ingroup Primitives
+ */
+typedef GenericOptInterval<IntCoord> OptIntInterval;
+
+} // namespace Geom
+#endif // !LIB2GEOM_SEEN_INT_INTERVAL_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/int-point.h b/src/2geom/int-point.h
new file mode 100644
index 000000000..cf2fe720f
--- /dev/null
+++ b/src/2geom/int-point.h
@@ -0,0 +1,157 @@
+/**
+ *  \file
+ *  \brief Cartesian point / 2D vector with integer coordinates
+ *//*
+ * Copyright 2011 Krzysztof Kosiński <tweenk.pl@gmail.com>
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, 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_INT_POINT_H
+#define LIB2GEOM_SEEN_INT_POINT_H
+
+#include <stdexcept>
+#include <boost/operators.hpp>
+#include <2geom/coord.h>
+
+namespace Geom {
+
+/**
+ * @brief Two-dimensional point with integer coordinates.
+ *
+ * This class is an exact equivalent of Point, except it stores integer coordinates.
+ * Integer points are useful in contexts related to rasterized graphics, for example
+ * for bounding boxes when rendering SVG.
+ *
+ * @see Point
+ * @ingroup Primitives */
+class IntPoint
+    : boost::additive< IntPoint
+    , boost::totally_ordered< IntPoint
+    > >
+{
+    IntCoord _pt[2];
+public:
+    /// @name Creating integer points
+    /// @{
+    IntPoint() { }
+    IntPoint(IntCoord x, IntCoord y) {
+        _pt[X] = x;
+        _pt[Y] = y;
+    }
+    IntPoint(IntPoint const &p) {
+        _pt[X] = p._pt[X];
+        _pt[Y] = p._pt[Y];
+    }
+    IntPoint &operator=(IntPoint const &p) {
+        _pt[X] = p._pt[X];
+        _pt[Y] = p._pt[Y];
+        return *this;
+    }
+    /// @}
+
+    /// @name Access the coordinates of a point
+    /// @{
+    IntCoord operator[](unsigned i) const {
+        if ( i > Y ) throw std::out_of_range("index out of range");
+        return _pt[i];
+    }
+    IntCoord &operator[](unsigned i) {
+        if ( i > Y ) throw std::out_of_range("index out of range");
+        return _pt[i];
+    }
+    IntCoord operator[](Dim2 d) const { return _pt[d]; }
+    IntCoord &operator[](Dim2 d) { return _pt[d]; }
+    /// @}
+
+    /// @name Vector-like arithmetic operations
+    /// @{
+    IntPoint &operator+=(IntPoint const &o) {
+        _pt[X] += o._pt[X];
+        _pt[Y] += o._pt[Y];
+        return *this;
+    }
+    IntPoint &operator-=(IntPoint const &o) {
+        _pt[X] -= o._pt[X];
+        _pt[Y] -= o._pt[Y];
+        return *this;
+    }
+    /// @}
+    
+    /// @name Various utilities
+    /// @{
+    /** @brief Equality operator. */
+    bool operator==(IntPoint const &in_pnt) const {
+        return ((_pt[X] == in_pnt[X]) && (_pt[Y] == in_pnt[Y]));
+    }
+    /** @brief Lexicographical ordering for points.
+     * Y coordinate is regarded as more significant. When sorting according to this
+     * ordering, the points will be sorted according to the Y coordinate, and within
+     * points with the same Y coordinate according to the X coordinate. */
+    bool operator<(IntPoint const &p) const {
+        return ( ( _pt[Y] < p[Y] ) ||
+             (( _pt[Y] == p[Y] ) && ( _pt[X] < p[X] )));
+    }
+    /// @}
+    
+    /** @brief Lexicographical ordering functor. */
+    template <Dim2 d> struct LexOrder;
+    /** @brief Lexicographical ordering functor with runtime dimension. */
+    class LexOrderRt {
+    public:
+        LexOrderRt(Dim2 d) : dim(d) {}
+        inline bool operator()(IntPoint const &a, IntPoint const &b);
+    private:
+        Dim2 dim;
+    };
+};
+
+template<> struct IntPoint::LexOrder<X> {
+    bool operator()(IntPoint const &a, IntPoint const &b) {
+        return a[X] < b[X] || (a[X] == b[X] && a[Y] < b[Y]);
+    }
+};
+template<> struct IntPoint::LexOrder<Y> {
+    bool operator()(IntPoint const &a, IntPoint const &b) {
+        return a[Y] < b[Y] || (a[Y] == b[Y] && a[X] < b[X]);
+    }
+};
+inline bool IntPoint::LexOrderRt::operator()(IntPoint const &a, IntPoint const &b) {
+    return dim ? IntPoint::LexOrder<Y>()(a, b) : IntPoint::LexOrder<X>()(a, b);
+}
+
+}  // namespace Geom
+
+#endif // !SEEN_GEOM_INT_POINT_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/int-rect.h b/src/2geom/int-rect.h
new file mode 100644
index 000000000..27fb06dfe
--- /dev/null
+++ b/src/2geom/int-rect.h
@@ -0,0 +1,74 @@
+/**
+ *  \file
+ *  \brief Axis-aligned rectangle with integer coordinates
+ *//*
+ * Copyright 2011 Krzysztof Kosiński <tweenk.pl@gmail.com>
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, 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_INT_RECT_H
+#define LIB2GEOM_SEEN_INT_RECT_H
+
+#include <2geom/coord.h>
+#include <2geom/generic-rect.h>
+
+namespace Geom {
+
+typedef GenericRect<IntCoord> IntRect;
+typedef GenericOptRect<IntCoord> OptIntRect;
+
+// the functions below do not work when defined generically
+inline OptIntRect operator&(IntRect const &a, IntRect const &b) {
+    OptIntRect ret(a);
+    ret.intersectWith(b);
+    return ret;
+}
+inline OptIntRect intersect(IntRect const &a, IntRect const &b) {
+    return a & b;
+}
+inline OptIntRect intersect(OptIntRect const &a, OptIntRect const &b) {
+    return a & b;
+}
+inline IntRect unify(IntRect const &a, IntRect const &b) {
+    return a | b;
+}
+inline OptIntRect unify(OptIntRect const &a, OptIntRect const &b) {
+    return a | b;
+}
+
+} // end namespace Geom
+
+#endif // !LIB2GEOM_SEEN_INT_RECT_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/interval.h b/src/2geom/interval.h
index a790a6c3b..ee6d674d2 100644
--- a/src/2geom/interval.h
+++ b/src/2geom/interval.h
@@ -37,19 +37,24 @@
 #ifndef LIB2GEOM_SEEN_INTERVAL_H
 #define LIB2GEOM_SEEN_INTERVAL_H
 
-#include <assert.h>
 #include <boost/none.hpp>
 #include <boost/optional.hpp>
 #include <boost/operators.hpp>
 #include <2geom/coord.h>
-#include <2geom/isnan.h>
+#include <2geom/math-utils.h>
+#include <2geom/generic-interval.h>
+#include <2geom/int-interval.h>
 
 namespace Geom {
 
-class OptInterval;
+/**
+ * @brief Range of real numbers that can be empty.
+ * @ingroup Primitives
+ */
+typedef GenericOptInterval<Coord> OptInterval;
 
 /** 
- * @brief Range of numbers that is never empty.
+ * @brief Range of real numbers that is never empty.
  *
  * Intervals are closed ranges \f$[a, b]\f$, which means they include their endpoints.
  * To use them as open ranges, you can use the interiorContains() methods.
@@ -57,32 +62,24 @@ class OptInterval;
  * @ingroup Primitives
  */
 class Interval
-    : boost::equality_comparable< Interval
-    , boost::additive< Interval
+    : public GenericInterval<Coord>
     , boost::multipliable< Interval
-    , boost::arithmetic< Interval, Coord
-    , boost::orable< Interval
-      > > > > >
+    , boost::multipliable< Interval, Coord
+      > >
 {
-private:
-    /// @invariant _b[0] <= _b[1]
-    Coord _b[2];
-
+    typedef GenericInterval<Coord> Base;
 public:
     /// @name Create intervals.
     /// @{
     /** @brief Create an interval that contains only zero. */
-    explicit Interval() { _b[0] = 0;  _b[1] = 0; }
+    Interval() {}
     /** @brief Create an interval that contains a single point. */
-    explicit Interval(Coord u) { _b[0] = _b[1] = u; }
+    explicit Interval(Coord u) : Base(u) {}
     /** @brief Create an interval that contains all points between @c u and @c v. */
-    Interval(Coord u, Coord v) {
-        if (u <= v) {
-            _b[0] = u; _b[1] = v;
-        } else {
-            _b[0] = v; _b[1] = u;
-        }
-    }
+    Interval(Coord u, Coord v) : Base(u,v) {}
+    /** @brief Convert from integer interval */
+    Interval(IntInterval const &i) : Base(i.min(), i.max()) {}
+    Interval(Base const &b) : Base(b) {}
 
     /** @brief Create an interval containing a range of values.
      * The resulting interval will contain all values from the given range.
@@ -93,9 +90,7 @@ public:
      * @return Interval that contains all values from [start, end). */
     template <typename InputIterator>
     static Interval from_range(InputIterator start, InputIterator end) {
-        assert(start != end);
-        Interval result(*start++);
-        for (; start != end; ++start) result.expandTo(*start);
+        Interval result = Base::from_range(start, end);
         return result;
     }
     /** @brief Create an interval from a C-style array of values it should contain. */
@@ -107,114 +102,38 @@ public:
 
     /// @name Inspect endpoints.
     /// @{
+    /** @brief Access endpoints by value.
+     * @deprecated Use min() and max() instead */
     Coord operator[](unsigned i) const { return _b[i]; }
+    /** @brief Access endpoints by reference.
+     * @deprecated Use min() and max() instead
+     * @todo Remove Interval index operator, which can be used to break the invariant */
     Coord& operator[](unsigned i) { return _b[i]; }
 
-    Coord min() const { return _b[0]; }
-    Coord max() const { return _b[1]; }
-    Coord extent() const { return _b[1] - _b[0]; }
-    Coord middle() const { return (_b[1] + _b[0]) * 0.5; }
-    bool isSingular() const { return _b[0] == _b[1]; }
     bool isFinite() const {
-        return IS_FINITE(_b[0]) && IS_FINITE(_b[1]);
+        return IS_FINITE(min()) && IS_FINITE(max());
     }
     /// @}
 
     /// @name Test coordinates and other intervals for inclusion.
     /// @{
-    /** @brief Check whether the interval includes this number. */
-    bool contains(Coord val) const { return _b[0] <= val && val <= _b[1]; }
     /** @brief Check whether the interior of the interval includes this number.
      * Interior means all numbers in the interval except its ends. */
-    bool interiorContains(Coord val) const { return _b[0] < val && val < _b[1]; }
-    /** @brief Check whether the interval includes the given interval. */
-    bool contains(Interval const &val) const { return _b[0] <= val._b[0] && val._b[1] <= _b[1]; }
+    bool interiorContains(Coord val) const { return min() < val && val < max(); }
     /** @brief Check whether the interior of the interval includes the given interval.
      * Interior means all numbers in the interval except its ends. */
-    bool interiorContains(Interval const &val) const { return _b[0] < val._b[0] && val._b[1] < _b[1]; }
-    /** @brief Check whether the intervals have any common elements. */
-    bool intersects(Interval const &val) const {
-        return contains(val._b[0]) || contains(val._b[1]) || val.contains(*this);
-    }
+    bool interiorContains(Interval const &val) const { return min() < val.min() && val.max() < max(); }
     /** @brief Check whether the interiors of the intervals have any common elements. */
     bool interiorIntersects(Interval const &val) const {
-        return interiorContains(val._b[0]) || interiorContains(val._b[1]) || val.interiorContains(*this);
-    }
-    /// @}
-
-    /// @name Modify the interval.
-    /// @{
-    //TODO: NaN handleage for the next two?
-    /** @brief Set the lower boundary of the interval.
-     * When the given number is larger than the interval's largest element,
-     * it will be reduced to the single number @c val. */
-    void setMin(Coord val) {
-        if(val > _b[1]) {
-            _b[0] = _b[1] = val;
-        } else {
-            _b[0] = val;
-        }
-    }
-    /** @brief Set the upper boundary of the interval.
-     * When the given number is smaller than the interval's smallest element,
-     * it will be reduced to the single number @c val. */
-    void setMax(Coord val) {
-        if(val < _b[0]) {
-            _b[1] = _b[0] = val;
-        } else {
-            _b[1] = val;
-        }
-    }
-    /** @brief Extend the interval to include the given number. */
-    void expandTo(Coord val) {
-       if(val < _b[0]) _b[0] = val;
-       if(val > _b[1]) _b[1] = val;  //no else, as we want to handle NaN
-    }
-    /** @brief Expand or shrink the interval in both directions by the given amount.
-     * After this method, the interval's length (extent) will be increased by
-     * <code>amount * 2</code>. Negative values can be given; they will shrink the interval.
-     * Shrinking by a value larger than half the interval's length will create a degenerate
-     * interval containing only the midpoint of the original. */
-    void expandBy(double amount) {
-        _b[0] -= amount;
-        _b[1] += amount;
-        if (_b[0] > _b[1]) {
-            Coord halfway = (_b[0]+_b[1])/2;
-            _b[0] = _b[1] = halfway;
-        }
-    }
-    /** @brief Union the interval with another one.
-     * The resulting interval will contain all points of both intervals.
-     * It might also contain some points which didn't belong to either - this happens
-     * when the intervals did not have any common elements. */
-    void unionWith(const Interval & a) {
-        if(a._b[0] < _b[0]) _b[0] = a._b[0];
-        if(a._b[1] > _b[1]) _b[1] = a._b[1];
+        return interiorContains(val.min()) || interiorContains(val.max()) || val.interiorContains(*this);
     }
     /// @}
 
     /// @name Operators
     /// @{
-    inline operator OptInterval();
-    bool operator==(Interval const &other) const { return _b[0] == other._b[0] && _b[1] == other._b[1]; }
-    
-    //IMPL: OffsetableConcept
-    //TODO: rename output_type to something else in the concept
-    typedef Coord output_type;
-    /** @brief Offset the interval by a specified amount */
-    Interval &operator+=(Coord amnt) {
-        _b[0] += amnt; _b[1] += amnt;
-        return *this;
-    }
-    /** @brief Offset the interval by the negation of the specified amount */
-    Interval &operator-=(Coord amnt) {
-        _b[0] -= amnt; _b[1] -= amnt;
-        return *this;
-    }
+    inline operator OptInterval() { return OptInterval(*this); }
     
     // IMPL: ScalableConcept
-    /** @brief Return an interval mirrored about 0 */
-    Interval operator-() const { return Interval(-_b[1], -_b[0]); }
     /** @brief Scale an interval */
     Interval &operator*=(Coord s) {
         _b[0] *= s;
@@ -229,25 +148,6 @@ public:
         if(s < 0) std::swap(_b[0], _b[1]);
         return *this;
     }
-    // IMPL: AddableConcept
-    /** @brief Add two intervals.
-     * Sum is defined as the set of points that can be obtained by adding any two values
-     * from both operands: \f$S = \{x \in A, y \in B: x + y\}\f$ */
-    Interval &operator+=(Interval const &o) {
-        _b[0] += o._b[0];
-        _b[1] += o._b[1];
-        return *this;
-    }
-    /** @brief Subtract two intervals.
-     * Difference is defined as the set of points that can be obtained by subtracting
-     * any value from the second operand from any value from the first operand:
-     * \f$S = \{x \in A, y \in B: x - y\}\f$ */
-    Interval &operator-=(Interval const &o) {
-        // equal to *this += -o
-        _b[0] -= o._b[1];
-        _b[1] -= o._b[0];
-        return *this;
-    }
     /** @brief Multiply two intervals.
      * Product is defined as the set of points that can be obtained by multiplying
      * any value from the second operand by any value from the first operand:
@@ -261,121 +161,25 @@ public:
         expandTo(mx * o.max());
         return *this;
     }
-    /** @brief Union two intervals.
-     * Note that intersection is only defined for OptIntervals, because the result
-     * of an intersection can be empty, while an Interval cannot. */
-    Interval &operator|=(Interval const &o) {
-        unionWith(o);
-        return *this;
-    }
     /// @}
-};
-
-/** @brief Union two intervals
- * @relates Interval */
-inline Interval unify(Interval const &a, Interval const &b) {
-    return a | b;
-}
-
-/**
- * @brief A range of numbers that can be empty.
- * @ingroup Primitives
- */
-class OptInterval
-    : public boost::optional<Interval>
-    , boost::orable< OptInterval
-    , boost::andable< OptInterval
-      > >
-{
-public:
-    /// @name Create optionally empty intervals.
+    
+    /// @name Rounding to integer values
     /// @{
-    /** @brief Create an empty interval. */
-    OptInterval() : boost::optional<Interval>() {};
-    /** @brief Wrap an existing interval. */
-    OptInterval(Interval const &a) : boost::optional<Interval>(a) {};
-    /** @brief Create an interval containing a single point. */
-    OptInterval(Coord u) : boost::optional<Interval>(Interval(u)) {};
-    /** @brief Create an interval containing a range of numbers. */
-    OptInterval(Coord u, Coord v) : boost::optional<Interval>(Interval(u,v)) {};
-
-    /** @brief Create a possibly empty interval containing a range of values.
-     * The resulting interval will contain all values from the given range.
-     * The return type of iterators must be convertible to double. The given range
-     * may be empty.
-     * @param start Beginning of the range
-     * @param end   End of the range
-     * @return Interval that contains all values from [start, end), or nothing if the range
-     *         is empty. */
-    template <typename InputIterator>
-    static OptInterval from_range(InputIterator start, InputIterator end) {
-        if (start == end) {
-            OptInterval ret;
-            return ret;
-        }
-        OptInterval ret(Interval::from_range(start, end));
+    /** @brief Return the smallest integer interval which contains this one. */
+    IntInterval roundOutwards() const {
+        IntInterval ret(floor(min()), ceil(max()));
         return ret;
     }
-    /// @}
-
-    /** @brief Check whether this OptInterval is empty. */
-    bool isEmpty() { return !*this; };
-
-    /** @brief Union with another interval, gracefully handling empty ones. */
-    inline void unionWith(OptInterval const &a) {
-        if (a) {
-            if (*this) { // check that we are not empty
-                (*this)->unionWith(*a);
-            } else {
-                *this = a;
-            }
-        }
-    }
-    inline void intersectWith(OptInterval const &o) {
-        if (o && *this) {
-            Coord u, v;
-            u = std::max((*this)->min(), o->min());
-            v = std::min((*this)->max(), o->max());
-            if (u <= v) {
-                *this = Interval(u, v);
-                return;
-            }
-        }
-        (*static_cast<boost::optional<Interval>*>(this)) = boost::none;
-    }
-    OptInterval &operator|=(OptInterval const &o) {
-        unionWith(o);
-        return *this;
-    }
-    OptInterval &operator&=(OptInterval const &o) {
-        intersectWith(o);
-        return *this;
+    /** @brief Return the largest integer interval which is contained in this one. */
+    OptIntInterval roundInwards() const {
+        IntCoord u = ceil(min()), v = floor(max());
+        if (u > v) { OptIntInterval e; return e; }
+        IntInterval ret(u, v);
+        return ret;
     }
+    /// @}
 };
 
-/** @brief Intersect two intervals and return a possibly empty range of numbers
- * @relates OptInterval */
-inline OptInterval intersect(Interval const &a, Interval const &b) {
-    return OptInterval(a) & OptInterval(b);
-}
-/** @brief Intersect two intervals and return a possibly empty range of numbers
- * @relates OptInterval */
-inline OptInterval operator&(Interval const &a, Interval const &b) {
-    return OptInterval(a) & OptInterval(b);
-}
-
-inline Interval::operator OptInterval() {
-    return OptInterval(*this);
-}
-
-#ifdef _GLIBCXX_IOSTREAM
-inline std::ostream &operator<< (std::ostream &os, 
-                                 const Geom::Interval &I) {
-    os << "Interval("<<I[0] << ", "<<I[1] << ")";
-    return os;
-}
-#endif
-
 }
 #endif //SEEN_INTERVAL_H
 
diff --git a/src/2geom/isnan.h b/src/2geom/isnan.h
deleted file mode 100644
index e20ab7f87..000000000
--- a/src/2geom/isnan.h
+++ /dev/null
@@ -1,116 +0,0 @@
-/**
- * \file
- * \brief  \todo brief description
- *
- * 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 _2GEOM_ISNAN_H__
-#define _2GEOM_ISNAN_H__
-
-/*
- * Temporary fix for various misdefinitions of isnan().
- * isnan() is becoming undef'd in some .h files. 
- * #include this last in your .cpp file to get it right.
- *
- * The problem is that isnan and isfinite are part of C99 but aren't part of
- * the C++ standard (which predates C99).
- *
- * Authors:
- *   Inkscape groupies and obsessive-compulsives
- *
- * Copyright (C) 2004 authors
- *
- * Released under GNU GPL, read the file 'COPYING' for more information
- *
- * 2005 modification hereby placed in public domain.  Probably supercedes 
- * the 2004 copyright for the code itself.
- */
-
-#include <math.h>
-/* You might try changing the above to <cmath> if you have problems.
- * Whether you use math.h or cmath, you may need to edit the .cpp file
- * and/or other .h files to use the same header file.
- */
-
-#if defined(__isnan)
-# define IS_NAN(_a) (__isnan(_a))
-#elif defined(__APPLE__) && __GNUC__ == 3
-# define IS_NAN(_a) (__isnan(_a))	/* MacOSX/Darwin definition < 10.4 */
-#elif defined(WIN32) || defined(_isnan)
-# define IS_NAN(_a) (_isnan(_a)) 	/* Win32 definition */
-#elif defined(isnan) || defined(__FreeBSD__) || defined(__osf__)
-# define IS_NAN(_a) (isnan(_a))		/* GNU definition */
-#elif defined (SOLARIS_2_8) && __GNUC__ == 3 && __GNUC_MINOR__ == 2
-# define IS_NAN(_a) (isnan(_a))		/* GNU definition */
-#else
-# define IS_NAN(_a) (std::isnan(_a))
-#endif
-/* If the above doesn't work, then try (a != a).
- * Also, please report a bug as per http://www.inkscape.org/report_bugs.php,
- * giving information about what platform and compiler version you're using.
- */
-
-
-#if defined(__isfinite)
-# define IS_FINITE(_a) (__isfinite(_a))
-#elif defined(__APPLE__) && __GNUC__ == 3
-# define IS_FINITE(_a) (__isfinite(_a))	/* MacOSX/Darwin definition < 10.4 */
-#elif defined(__sgi)
-# define IS_FINITE(_a) (_isfinite(_a))
-#elif defined(isfinite)
-# define IS_FINITE(_a) (isfinite(_a))
-#elif defined(__osf__)
-# define IS_FINITE(_a) (finite(_a) && !IS_NAN(_a))
-#elif defined (SOLARIS_2_8) && __GNUC__ == 3 && __GNUC_MINOR__ == 2
-#include  <ieeefp.h>
-#define IS_FINITE(_a) (finite(_a) && !IS_NAN(_a))
-#else
-# define IS_FINITE(_a) (std::isfinite(_a))
-#endif
-/* If the above doesn't work, then try (finite(_a) && !IS_NAN(_a)) or 
- * (!IS_NAN((_a) - (_a))).
- * Also, please report a bug as per http://www.inkscape.org/report_bugs.php,
- * giving information about what platform and compiler version you're using.
- */
-
-
-#endif /* _2GEOM_ISNAN_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/line.h b/src/2geom/line.h
index ccb0ae6c5..f2d31ecc6 100644
--- a/src/2geom/line.h
+++ b/src/2geom/line.h
@@ -1,8 +1,11 @@
 /**
  * \file
- * \brief  Infinite Straight Line
- *
- * Copyright 2008  Marco Cecchetti <mrcekets at gmail.com>
+ * \brief  Infinite straight line
+ *//*
+ * Authors:
+ *   Marco Cecchetti <mrcekets at gmail.com>
+ *   Krzysztof Kosiński <tweenk.pl@gmail.com>
+ * Copyright 2008-2011 Authors
  *
  * This library is free software; you can redistribute it and/or
  * modify it either under the terms of the GNU Lesser General Public
@@ -28,22 +31,17 @@
  * the specific language governing rights and limitations.
  */
 
-#ifndef _2GEOM_LINE_H_
-#define _2GEOM_LINE_H_
-
+#ifndef LIB2GEOM_SEEN_LINE_H
+#define LIB2GEOM_SEEN_LINE_H
 
 #include <cmath>
-
+#include <boost/optional.hpp>
 #include <2geom/bezier-curve.h> // for LineSegment
 #include <2geom/rect.h>
 #include <2geom/crossing.h>
 #include <2geom/exception.h>
-
 #include <2geom/ray.h>
 
-#include <boost/optional.hpp>
-
-
 namespace Geom
 {
 
@@ -226,8 +224,8 @@ public:
      * @return Ray starting at t and going in the direction of the versor */
     Ray ray(Coord t) {
         Ray result;
-        result.origin(pointAt(t));
-        result.versor(m_versor);
+        result.setOrigin(pointAt(t));
+        result.setVersor(m_versor);
         return result;
     }
 
@@ -448,17 +446,16 @@ OptCrossing intersection(LineSegment const& ls1, LineSegment const& ls2);
 } // end namespace Geom
 
 
-#endif // _2GEOM_LINE_H_
+#endif // LIB2GEOM_SEEN_LINE_H
 
 
 /*
   Local Variables:
   mode:c++
   c-file-style:"stroustrup"
-  c-file-offsets:((innamespace . 0)(substatement-open . 0))
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
   indent-tabs-mode:nil
-  c-brace-offset:0
   fill-column:99
   End:
-  vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4 :
 */
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
diff --git a/src/2geom/linear.h b/src/2geom/linear.h
index 1b6cca071..df6dd9904 100644
--- a/src/2geom/linear.h
+++ b/src/2geom/linear.h
@@ -35,7 +35,7 @@
 #ifndef SEEN_LINEAR_H
 #define SEEN_LINEAR_H
 #include <2geom/interval.h>
-#include <2geom/isnan.h>
+#include <2geom/math-utils.h>
 
 
 //#define USE_SBASIS_OF
diff --git a/src/2geom/math-utils.h b/src/2geom/math-utils.h
index 2c348f54b..77280aa50 100644
--- a/src/2geom/math-utils.h
+++ b/src/2geom/math-utils.h
@@ -1,6 +1,3 @@
-#ifndef LIB2GEOM_MATH_UTILS_HEADER
-#define LIB2GEOM_MATH_UTILS_HEADER
-
 /**
  * \file
  * \brief Low level math functions and compatibility wrappers
@@ -36,6 +33,9 @@
  *
  */
 
+#ifndef LIB2GEOM_SEEN_MATH_UTILS_H
+#define LIB2GEOM_SEEN_MATH_UTILS_H
+
 #include "config.h"
 #include <math.h> // sincos is usually only available in math.h
 #include <cmath>
@@ -92,10 +92,51 @@ inline void sincos(double angle, double &sin_, double &cos_) {
 #endif
 }
 
-}
+/* Temporary fix for various misdefinitions of isnan().
+ * isnan() is becoming undef'd in some .h files. 
+ * #include this last in your .cpp file to get it right.
+ *
+ * The problem is that isnan and isfinite are part of C99 but aren't part of
+ * the C++ standard (which predates C99).
+ */
+
+#if defined(__isnan)
+# define IS_NAN(_a) (__isnan(_a))
+#elif defined(__APPLE__) && __GNUC__ == 3
+# define IS_NAN(_a) (__isnan(_a))	/* MacOSX/Darwin definition < 10.4 */
+#elif defined(WIN32) || defined(_isnan)
+# define IS_NAN(_a) (_isnan(_a)) 	/* Win32 definition */
+#elif defined(isnan) || defined(__FreeBSD__) || defined(__osf__)
+# define IS_NAN(_a) (isnan(_a))		/* GNU definition */
+#elif defined (SOLARIS_2_8) && __GNUC__ == 3 && __GNUC_MINOR__ == 2
+# define IS_NAN(_a) (isnan(_a))		/* GNU definition */
+#else
+# define IS_NAN(_a) (std::isnan(_a))
+#endif
+/* If the above doesn't work, then try (a != a). */
+
 
+#if defined(__isfinite)
+# define IS_FINITE(_a) (__isfinite(_a))
+#elif defined(__APPLE__) && __GNUC__ == 3
+# define IS_FINITE(_a) (__isfinite(_a))	/* MacOSX/Darwin definition < 10.4 */
+#elif defined(__sgi)
+# define IS_FINITE(_a) (_isfinite(_a))
+#elif defined(isfinite)
+# define IS_FINITE(_a) (isfinite(_a))
+#elif defined(__osf__)
+# define IS_FINITE(_a) (finite(_a) && !IS_NAN(_a))
+#elif defined (SOLARIS_2_8) && __GNUC__ == 3 && __GNUC_MINOR__ == 2
+#include  <ieeefp.h>
+#define IS_FINITE(_a) (finite(_a) && !IS_NAN(_a))
+#else
+# define IS_FINITE(_a) (std::isfinite(_a))
 #endif
 
+} // end namespace Geom
+
+#endif // LIB2GEOM_SEEN_MATH_UTILS_H
+
 /*
   Local Variables:
   mode:c++
diff --git a/src/2geom/ord.h b/src/2geom/ord.h
index ca91af579..ce524ebf7 100644
--- a/src/2geom/ord.h
+++ b/src/2geom/ord.h
@@ -1,7 +1,7 @@
 /**
  * \file
- * \brief  \todo brief description
- *
+ * \brief Comparator template
+ *//*
  * Authors:
  *      ? <?@?.?>
  * 
diff --git a/src/2geom/path-intersection.h b/src/2geom/path-intersection.h
index de2a5b02c..2470e44fb 100644
--- a/src/2geom/path-intersection.h
+++ b/src/2geom/path-intersection.h
@@ -1,7 +1,7 @@
 /**
  * \file
- * \brief  \todo brief description
- *
+ * \brief Path intersection
+ *//*
  * Authors:
  *      ? <?@?.?>
  * 
diff --git a/src/2geom/path.h b/src/2geom/path.h
index cbd449248..48d7acaaf 100644
--- a/src/2geom/path.h
+++ b/src/2geom/path.h
@@ -1,12 +1,12 @@
 /**
  * \file
  * \brief  Path - Series of continuous curves
- *
+ *//*
  * Authors:
- * 		MenTaLguY <mental@rydia.net>
- * 		Marco Cecchetti <mrcekets at gmail.com>
+ *   MenTaLguY <mental@rydia.net>
+ *   Marco Cecchetti <mrcekets at gmail.com>
  * 
- * Copyright 2007-2008  authors
+ * Copyright 2007-2008 Authors
  *
  * This library is free software; you can redistribute it and/or
  * modify it either under the terms of the GNU Lesser General Public
@@ -32,23 +32,16 @@
  * the specific language governing rights and limitations.
  */
 
+#ifndef LIB2GEOM_SEEN_PATH_H
+#define LIB2GEOM_SEEN_PATH_H
 
-
-
-#ifndef SEEN_GEOM_PATH_H
-#define SEEN_GEOM_PATH_H
-
-
+#include <iterator>
+#include <algorithm>
 #include <boost/shared_ptr.hpp>
 #include <2geom/curve.h>
 #include <2geom/bezier-curve.h>
 
-#include <iterator>
-#include <algorithm>
-
-
-namespace Geom
-{
+namespace Geom {
 
 class Path;
 
@@ -696,17 +689,13 @@ Coord nearest_point(Point const& p, Path const& c)
 namespace std {
 
 template <>
-inline void swap<Geom::Path>(Geom::Path &a, Geom::Path &b)
-{
+inline void swap<Geom::Path>(Geom::Path &a, Geom::Path &b) {
   a.swap(b);
 }
 
 }  // end namespace std
 
-#endif // SEEN_GEOM_PATH_H
-
-
-
+#endif // LIB2GEOM_SEEN_PATH_H
 
 /*
   Local Variables:
diff --git a/src/2geom/pathvector.h b/src/2geom/pathvector.h
index 2f45b9d86..2b690a005 100644
--- a/src/2geom/pathvector.h
+++ b/src/2geom/pathvector.h
@@ -1,10 +1,9 @@
 /**
  * \file
- * \brief PathVector - std::vector containing Geom::Path
+ * \brief PathVector - std::vector containing Geom::Path.
  * This file provides a set of operations that can be performed on PathVector,
  * e.g. an affine transform.
- */
-/*
+ *//*
  * Authors:
  *  Johan Engelen <goejendaagh@zonnet.nl>
  * 
@@ -34,8 +33,8 @@
  * the specific language governing rights and limitations.
  */
 
-#ifndef SEEN_GEOM_PATHVECTOR_H
-#define SEEN_GEOM_PATHVECTOR_H
+#ifndef LIB2GEOM_SEEN_PATHVECTOR_H
+#define LIB2GEOM_SEEN_PATHVECTOR_H
 
 #include <2geom/forward.h>
 #include <2geom/path.h>
@@ -122,11 +121,9 @@ Point pointAt(PathVector const & path_in, PathVectorPosition const pvp) {
     return path_in[pvp.path_nr].pointAt(pvp.t);
 }
 
-
-
 } // end namespace Geom
 
-#endif // SEEN_GEOM_PATHVECTOR_H
+#endif // LIB2GEOM_SEEN_PATHVECTOR_H
 
 /*
   Local Variables:
diff --git a/src/2geom/piecewise.h b/src/2geom/piecewise.h
index 19c66d8f0..837f33ea7 100644
--- a/src/2geom/piecewise.h
+++ b/src/2geom/piecewise.h
@@ -32,13 +32,13 @@
 #ifndef SEEN_GEOM_PW_SB_H
 #define SEEN_GEOM_PW_SB_H
 
-#include <2geom/sbasis.h>
 #include <vector>
 #include <map>
-
-#include <2geom/concepts.h>
-#include <2geom/isnan.h>
 #include <boost/concept_check.hpp>
+#include <2geom/concepts.h>
+#include <2geom/math-utils.h>
+#include <2geom/sbasis.h>
+
 
 namespace Geom {
 /**
diff --git a/src/2geom/point.cpp b/src/2geom/point.cpp
index a9005ef61..cafc0fdba 100644
--- a/src/2geom/point.cpp
+++ b/src/2geom/point.cpp
@@ -1,3 +1,38 @@
+/**
+ *  \file
+ *  \brief Cartesian point / 2D vector and related operations
+ *//*
+ *  Authors:
+ *    Michael G. Sloan <mgsloan@gmail.com>
+ *    Nathan Hurst <njh@njhurst.com>
+ *    Krzysztof Kosiński <tweenk.pl@gmail.com>
+ *
+ * Copyright (C) 2006-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.
+ */
+
 #include <assert.h>
 #include <math.h>
 #include <2geom/point.h>
@@ -14,8 +49,8 @@ namespace Geom {
  * from the origin (point at 0,0) to the stored coordinates,
  * and has methods implementing several vector operations (like length()).
  *
- * \par Operator note
- * \par
+ * @par Operator note
+ * @par
  * Most operators are provided by Boost operator helpers, so they are not visible in this class.
  * If @a p, @a q, @a r denote points, @a s a floating-point scalar, and @a m a transformation matrix,
  * then the following operations are available:
@@ -149,12 +184,11 @@ Point unit_vector(Point const &a)
  * that the origin (0, 0), its negation is returned. You can check whether
  * the points' vectors have the same direction (e.g. lie
  * on the same line passing through the origin) using
- * @code abs(a).normalize() == abs(b).normalize() @endcode.
+ * @code abs(a).normalize() == abs(b).normalize() @endcode
  * To check with some margin of error, use
- * @code are_near(abs(a).normalize(), abs(b).normalize()) @endcode.
+ * @code are_near(abs(a).normalize(), abs(b).normalize()) @endcode
  * Although naively this should take the absolute value of each coordinate, such an operation
  * is not very useful.
- * @return \f$p' = (p_X, -p_Y)\f$
  * @relates Point */
 Point abs(Point const &b)
 {
@@ -178,8 +212,8 @@ Point &Point::operator*=(Affine const &m) {
     return *this;
 }
 
-/** @brief Snap the angle B - A - dir to miltiples of \f$2\pi/n\f$.
- * The 'dir' argument must be normalized (have an unit length), otherwise the result
+/** @brief Snap the angle B - A - dir to multiples of \f$2\pi/n\f$.
+ * The 'dir' argument must be normalized (have unit length), otherwise the result
  * is undefined.
  * @return Point with the same distance from A as B, with a snapped angle.
  * @post distance(A, B) == distance(A, result)
diff --git a/src/2geom/point.h b/src/2geom/point.h
index 3c6e12eff..69da8a4ae 100644
--- a/src/2geom/point.h
+++ b/src/2geom/point.h
@@ -42,7 +42,7 @@
 #include <boost/operators.hpp>
 #include <2geom/forward.h>
 #include <2geom/coord.h>
-#include <2geom/isnan.h> //temporary fix for isnan()
+#include <2geom/int-point.h>
 #include <2geom/math-utils.h>
 #include <2geom/utils.h>
 
@@ -61,8 +61,9 @@ class Point
       > > > > > > > > > // this uses chaining so it looks weird, but works
 {
     Coord _pt[2];
-
 public:
+    /// @name Create points
+    /// @{
     /** Construct a point on the origin. */
     Point()
     { _pt[X] = _pt[Y] = 0; }
@@ -71,6 +72,11 @@ public:
     Point(Coord x, Coord y) {
         _pt[X] = x; _pt[Y] = y;
     }
+    /** Construct from integer point. */
+    Point(IntPoint const &p) {
+        _pt[X] = p[X];
+        _pt[Y] = p[Y];
+    }
     Point(Point const &p) {
         for (unsigned i = 0; i < 2; ++i)
             _pt[i] = p._pt[i];
@@ -80,6 +86,23 @@ public:
             _pt[i] = p._pt[i];
         return *this;
     }
+    /** @brief Construct a point from its polar coordinates.
+     * The angle is specified in radians, in the mathematical convention (increasing
+     * counter-clockwise from +X). */
+    static Point polar(Coord angle, Coord radius) {
+        Point ret(polar(angle));
+        ret *= radius;
+        return ret;
+    }
+    /** @brief Construct an unit vector from its angle.
+     * The angle is specified in radians, in the mathematical convention (increasing
+     * counter-clockwise from +X). */
+    static Point polar(Coord angle) {
+        Point ret;
+        sincos(angle, ret[Y], ret[X]);
+        return ret;
+    }
+    /// @}
 
     /// @name Access the coordinates of a point
     /// @{
@@ -157,57 +180,54 @@ public:
     }
     /// @}
 
-    /// @name Various utilities
+    /// @name Conversion to integer points
     /// @{
-    /** @brief Lower the precision of the point.
-     * This will round both coordinates to multiples of \f$10^p\f$. */
-    void round (int p = 0) {
-        _pt[X] = (Coord)(decimal_round((double)_pt[X], p));
-        _pt[Y] = (Coord)(decimal_round((double)_pt[Y], p));
-        return;
+    /** @brief Round to nearest integer coordinates. */
+    IntPoint round() const {
+        IntPoint ret(::round(_pt[X]), ::round(_pt[Y]));
+        return ret;
+    }
+    /** @brief Round coordinates downwards. */
+    IntPoint floor() const {
+        IntPoint ret(::floor(_pt[X]), ::floor(_pt[Y]));
+        return ret;
+    }
+    /** @brief Round coordinates upwards. */
+    IntPoint ceil() const {
+        IntPoint ret(::ceil(_pt[X]), ::ceil(_pt[Y]));
+        return ret;
     }
+    /// @}
 
-    /** @brief Check whether both coordinates are finite.
-     * @return True if neither coordinate is infinite. */
+    /// @name Various utilities
+    /// @{
+    /** @brief Check whether both coordinates are finite. */
     bool isFinite() const {
         for ( unsigned i = 0 ; i < 2 ; ++i ) {
             if(!IS_FINITE(_pt[i])) return false;
         }
         return true;
     }
+    /** @brief Check whether both coordinates are zero. */
+    bool isZero() const {
+        return _pt[X] == 0 && _pt[Y] == 0;
+    }
+    /** @brief Check whether the length of the vector is close to 1. */
+    bool isNormalized(Coord eps=EPSILON) const {
+        return are_near(length(), 1.0, eps);
+    }
     /** @brief Equality operator.
      * This tests for exact identity (as opposed to are_near()). Note that due to numerical
      * errors, this test might return false even if the points should be identical. */
     bool operator==(const Point &in_pnt) const {
-        return ((_pt[X] == in_pnt[X]) && (_pt[Y] == in_pnt[Y]));
+        return (_pt[X] == in_pnt[X]) && (_pt[Y] == in_pnt[Y]);
     }
     /** @brief Lexicographical ordering for points.
      * Y coordinate is regarded as more significant. When sorting according to this
      * ordering, the points will be sorted according to the Y coordinate, and within
      * points with the same Y coordinate according to the X coordinate. */
     bool operator<(const Point &p) const {
-        return ( ( _pt[Y] < p[Y] ) ||
-             (( _pt[Y] == p[Y] ) && ( _pt[X] < p[X] )));
-    }
-    /// @}
-
-    /// @name Point factories
-    /// @{
-    /** @brief Construct a point from its polar coordinates.
-     * The angle is specified in radians, in the mathematical convention (increasing
-     * counter-clockwise from +X). */
-    static Point polar(Coord angle, Coord radius) {
-        Point ret(polar(angle));
-        ret *= radius;
-        return ret;
-    }
-    /** @brief Construct an unit vector from its angle.
-     * The angle is specified in radians, in the mathematical convention (increasing
-     * counter-clockwise from +X). */
-    static Point polar(Coord angle) {
-        Point ret;
-        sincos(angle, ret[Y], ret[X]);
-        return ret;
+        return _pt[Y] < p[Y] || (_pt[Y] == p[Y] && _pt[X] < p[X]);
     }
     /// @}
 
@@ -315,7 +335,7 @@ inline Point lerp(double const t, Point const &a, Point const &b)
  * For perpendicular vectors, it is zero. For parallel ones, its absolute value is highest,
  * and the sign depends on whether they point in the same direction (+) or opposite ones (-).
  * @return \f$a \cdot b = a_X b_X + a_Y b_Y\f$.
- * @relates Point*/
+ * @relates Point */
 inline Coord dot(Point const &a, Point const &b)
 {
     return a[0] * b[0] + a[1] * b[1];
@@ -349,8 +369,8 @@ Coord L1(Point const &p);
 Coord LInfty(Point const &p);
 bool is_zero(Point const &p);
 bool is_unit_vector(Point const &p);
-extern double atan2(Point const &p);
-extern double angle_between(Point const &a, Point const &b);
+double atan2(Point const &p);
+double angle_between(Point const &a, Point const &b);
 Point abs(Point const &b);
 Point constrain_angle(Point const &A, Point const &B, unsigned int n = 4, Geom::Point const &dir = Geom::Point(1,0));
 
diff --git a/src/2geom/poly.h b/src/2geom/poly.h
index 3567bda6d..7d93d0a85 100644
--- a/src/2geom/poly.h
+++ b/src/2geom/poly.h
@@ -1,7 +1,7 @@
 /**
  * \file
- * \brief  \todo brief description
- *
+ * \brief Polynomial in canonical (monomial) basis
+ *//*
  * Authors:
  *      ? <?@?.?>
  * 
@@ -34,7 +34,6 @@
 
 #ifndef LIB2GEOM_SEEN_POLY_H
 #define LIB2GEOM_SEEN_POLY_H
-
 #include <assert.h>
 #include <vector>
 #include <iostream>
diff --git a/src/2geom/quadtree.h b/src/2geom/quadtree.h
index 01ea33ed7..949a9b898 100644
--- a/src/2geom/quadtree.h
+++ b/src/2geom/quadtree.h
@@ -1,7 +1,7 @@
 /**
  * \file
- * \brief  \todo brief description
- *
+ * \brief Quad tree data structure
+ *//*
  * Authors:
  *      ? <?@?.?>
  * 
diff --git a/src/2geom/ray.h b/src/2geom/ray.h
index 638b86195..75cc72005 100644
--- a/src/2geom/ray.h
+++ b/src/2geom/ray.h
@@ -1,7 +1,7 @@
 /**
  * \file
- * \brief  Infinite Straight Ray
- *
+ * \brief  Infinite straight ray
+ *//*
  * Copyright 2008  Marco Cecchetti <mrcekets at gmail.com>
  *
  * This library is free software; you can redistribute it and/or
@@ -35,6 +35,7 @@
 #include <2geom/point.h>
 #include <2geom/bezier-curve.h> // for LineSegment
 #include <2geom/exception.h>
+#include <2geom/math-utils.h>
 
 namespace Geom
 {
@@ -49,171 +50,88 @@ namespace Geom
  */
 class Ray {
 private:
-    Point m_origin;
-    Point m_versor;
+    Point _origin;
+    Point _versor;
 
 public:
-	Ray()
-		: m_origin(0,0), m_versor(1,0)
-	{
-	}
-
-	Ray(Point const& _origin, Coord angle )
-		: m_origin(_origin), m_versor(std::cos(angle), std::sin(angle))
-	{
-	}
-
-	Ray(Point const& A, Point const& B)
-	{
-		setBy2Points(A, B);
-	}
-
-	Point origin() const
-	{
-		return m_origin;
-	}
-
-	Point versor() const
-	{
-		return m_versor;
-	}
-
-	void origin(Point const& _point)
-	{
-		m_origin = _point;
-	}
-
-	void versor(Point const& _versor)
-	{
-		m_versor = _versor;
-	}
-
-	Coord angle() const
-	{
-		double a = std::atan2(m_versor[Y], m_versor[X]);
-		if (a < 0) a += 2*M_PI;
-		return a;
-	}
-
-	void angle(Coord _angle)
-	{
-		m_versor[X] = std::cos(_angle);
-		m_versor[Y] = std::sin(_angle);
-	}
-
-	void setBy2Points(Point const& A, Point const& B)
-	{
-		m_origin = A;
-		m_versor = B - A;
-		if ( are_near(m_versor, Point(0,0)) )
-			m_versor = Point(0,0);
-		else
-			m_versor.normalize();
-	}
-
-	bool isDegenerate() const
-	{
-		return ( m_versor[X] == 0 && m_versor[Y] == 0 );
-	}
-
-	Point pointAt(Coord t) const
-	{
-		if (t < 0)	THROW_RANGEERROR("Ray::pointAt, negative t value passed");
-		return m_origin + m_versor * t;
-	}
-
-	Coord valueAt(Coord t, Dim2 d) const
-	{
-		if (t < 0)
-			THROW_RANGEERROR("Ray::valueAt, negative t value passed");
-		if (d < 0 || d > 1)
-			THROW_RANGEERROR("Ray::valueAt, dimension argument out of range");
-		return m_origin[d] + m_versor[d] * t;
-	}
-
-	std::vector<Coord> roots(Coord v, Dim2 d) const
-	{
-		if (d < 0 || d > 1)
-			THROW_RANGEERROR("Ray::roots, dimension argument out of range");
-		std::vector<Coord> result;
-		if ( m_versor[d] != 0 )
-		{
-			double t = (v - m_origin[d]) / m_versor[d];
-			if (t >= 0)	result.push_back(t);
-		}
-		// TODO: else ?
-		return result;
-	}
-
-	// require are_near(_point, *this)
-	// on the contrary the result value is meaningless
-	Coord timeAt(Point const& _point) const
-	{
-		Coord t;
-		if ( m_versor[X] != 0 )
-		{
-			t = (_point[X] - m_origin[X]) / m_versor[X];
-		}
-		else if ( m_versor[Y] != 0 )
-		{
-			t = (_point[Y] - m_origin[Y]) / m_versor[Y];
-		}
-		else // degenerate case
-		{
-			t = 0;
-		}
-		return t;
-	}
-
-	Coord nearestPoint(Point const& _point) const
-	{
-		if ( isDegenerate() ) return 0;
-		double t = dot( _point - m_origin, m_versor );
-		if (t < 0) t = 0;
-		return t;
-	}
-
-	Ray reverse() const
-	{
-		Ray result;
-		result.origin(m_origin);
-		result.versor(-m_versor);
-		return result;
-	}
-
-	Curve* portion(Coord  f, Coord t) const
-	{
-		LineSegment* seg = new LineSegment(pointAt(f), pointAt(t));
-		return seg;
-	}
-
-	LineSegment segment(Coord  f, Coord t) const
-	{
-		return LineSegment(pointAt(f), pointAt(t));
-	}
-
-	Ray transformed(Affine const& m) const
-	{
-		return Ray(m_origin * m, (m_origin + m_versor) * m);
-	}
-};  // end class ray
+    Ray() : _origin(0,0), _versor(1,0) {}
+    Ray(Point const& origin, Coord angle)
+        : _origin(origin)
+    {
+        sincos(angle, _versor[Y], _versor[X]);
+    }
+    Ray(Point const& A, Point const& B) {
+        setPoints(A, B);
+    }
+    Point origin() const { return _origin; }
+    Point versor() const { return _versor; }
+    void setOrigin(Point const &o) { _origin = o; }
+    void setVersor(Point const& v) { _versor = v; }
+    Coord angle() const { return std::atan2(_versor[Y], _versor[X]); }
+    void setAngle(Coord a) { sincos(a, _versor[Y], _versor[X]); }
+    void setPoints(Point const &a, Point const &b) {
+        _origin = a;
+        _versor = b - a;
+        if (are_near(_versor, Point(0,0)) )
+            _versor = Point(0,0);
+        else
+            _versor.normalize();
+    }
+    bool isDegenerate() const {
+        return ( _versor[X] == 0 && _versor[Y] == 0 );
+    }
+    Point pointAt(Coord t) const {
+        return _origin + _versor * t;
+    }
+    Coord valueAt(Coord t, Dim2 d) const {
+        return _origin[d] + _versor[d] * t;
+    }
+    std::vector<Coord> roots(Coord v, Dim2 d) const {
+        std::vector<Coord> result;
+        if ( _versor[d] != 0 ) {
+            double t = (v - _origin[d]) / _versor[d];
+            if (t >= 0)	result.push_back(t);
+        } else if (_versor[(d+1)%2] == v) {
+            THROW_INFINITESOLUTIONS();
+        }
+        return result;
+    }
+    Coord nearestPoint(Point const& point) const {
+        if ( isDegenerate() ) return 0;
+        double t = dot(point - _origin, _versor);
+        if (t < 0) t = 0;
+        return t;
+    }
+    Ray reverse() const {
+        Ray result;
+        result.setOrigin(_origin);
+        result.setVersor(-_versor);
+        return result;
+    }
+    Curve *portion(Coord f, Coord t) const {
+        return new LineSegment(pointAt(f), pointAt(t));
+    }
+    LineSegment segment(Coord f, Coord t) const {
+        return LineSegment(pointAt(f), pointAt(t));
+    }
+    Ray transformed(Affine const& m) const {
+        return Ray(_origin * m, (_origin + _versor) * m);
+    }
+}; // end class Ray
 
 inline
-double distance(Point const& _point, Ray const& _ray)
-{
+double distance(Point const& _point, Ray const& _ray) {
 	double t = _ray.nearestPoint(_point);
 	return ::Geom::distance(_point, _ray.pointAt(t));
 }
 
 inline
-bool are_near(Point const& _point, Ray const& _ray, double eps = EPSILON)
-{
+bool are_near(Point const& _point, Ray const& _ray, double eps = EPSILON) {
 	return are_near(distance(_point, _ray), 0, eps);
 }
 
 inline
-bool are_same(Ray const& r1, Ray const& r2, double eps = EPSILON)
-{
+bool are_same(Ray const& r1, Ray const& r2, double eps = EPSILON) {
 	return are_near(r1.versor(), r2.versor(), eps)
 			&& are_near(r1.origin(), r2.origin(), eps);
 }
@@ -221,15 +139,13 @@ bool are_same(Ray const& r1, Ray const& r2, double eps = EPSILON)
 // evaluate the angle between r1 and r2 rotating r1 in cw or ccw direction on r2
 // the returned value is an angle in the interval [0, 2PI[
 inline
-double angle_between(Ray const& r1, Ray const& r2, bool cw = true)
-{
+double angle_between(Ray const& r1, Ray const& r2, bool cw = true) {
 	double angle = angle_between(r1.versor(), r2.versor());
 	if (angle < 0) angle += 2*M_PI;
 	if (!cw) angle = 2*M_PI - angle;
 	return angle;
 }
 
-
 inline
 Ray make_angle_bisector_ray(Ray const& r1, Ray const& r2)
 {
@@ -243,22 +159,17 @@ Ray make_angle_bisector_ray(Ray const& r1, Ray const& r2)
     return Ray(r1.origin(), M);
 }
 
-
 }  // end namespace Geom
 
-
-
-#endif /*_2GEOM_RAY_H_*/
-
+#endif // LIB2GEOM_SEEN_RAY_H
 
 /*
   Local Variables:
   mode:c++
   c-file-style:"stroustrup"
-  c-file-offsets:((innamespace . 0)(substatement-open . 0))
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
   indent-tabs-mode:nil
-  c-brace-offset:0
   fill-column:99
   End:
-  vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4 :
 */
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
diff --git a/src/2geom/rect.cpp b/src/2geom/rect.cpp
new file mode 100644
index 000000000..0cb842d29
--- /dev/null
+++ b/src/2geom/rect.cpp
@@ -0,0 +1,98 @@
+/* Axis-aligned rectangle
+ *
+ * Authors:
+ *   Michael Sloan <mgsloan@gmail.com>
+ *   Krzysztof Kosiński <tweenk.pl@gmail.com>
+ * Copyright 2007-2011 Authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+#include <2geom/rect.h>
+
+namespace Geom {
+
+/** @brief Transform the rectangle by an affine.
+ * The result of the transformation might not be axis-aligned. The return value
+ * of this operation will be the smallest axis-aligned rectangle containing
+ * all points of the true result. */
+Rect &Rect::operator*=(Affine const &m) {
+    Point pts[4];
+    for (unsigned i=0; i<4; ++i) pts[i] = corner(i) * m;
+    Coord minx = std::min(std::min(pts[0][X], pts[1][X]), std::min(pts[2][X], pts[3][X]));
+    Coord miny = std::min(std::min(pts[0][Y], pts[1][Y]), std::min(pts[2][Y], pts[3][Y]));
+    Coord maxx = std::max(std::max(pts[0][X], pts[1][X]), std::max(pts[2][X], pts[3][X]));
+    Coord maxy = std::max(std::max(pts[0][Y], pts[1][Y]), std::max(pts[2][Y], pts[3][Y]));
+    f[X].setMin(minx); f[X].setMax(maxx);
+    f[Y].setMin(miny); f[Y].setMax(maxy);
+    return *this;
+}
+
+Coord distanceSq(Point const &p, Rect const &rect)
+{
+    double dx = 0, dy = 0;
+    if ( p[X] < rect.left() ) {
+        dx = p[X] - rect.left();
+    } else if ( p[X] > rect.right() ) {
+        dx = rect.right() - p[X];
+    }
+    if (p[Y] < rect.top() ) {
+        dy = rect.top() - p[Y];
+    } else if (  p[Y] > rect.bottom() ) {
+        dy = p[Y] - rect.bottom();
+    }
+    return dx*dx+dy*dy;
+}
+
+/** @brief Returns the smallest distance between p and rect.
+ * @relates Rect */
+Coord distance(Point const &p, Rect const &rect)
+{
+    // copy of distanceSq, because we need to use hypot()
+    double dx = 0, dy = 0;
+    if ( p[X] < rect.left() ) {
+        dx = p[X] - rect.left();
+    } else if ( p[X] > rect.right() ) {
+        dx = rect.right() - p[X];
+    }
+    if (p[Y] < rect.top() ) {
+        dy = rect.top() - p[Y];
+    } else if (  p[Y] > rect.bottom() ) {
+        dy = p[Y] - rect.bottom();
+    }
+    return hypot(dx, dy);
+}
+
+} // namespace Geom
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
diff --git a/src/2geom/rect.h b/src/2geom/rect.h
index 65bb1bb76..72b659a81 100644
--- a/src/2geom/rect.h
+++ b/src/2geom/rect.h
@@ -2,7 +2,10 @@
  * \file
  * \brief Axis-aligned rectangle
  *//*
- * Copyright 2007 Michael Sloan <mgsloan@gmail.com>
+ * Authors:
+ *   Michael Sloan <mgsloan@gmail.com>
+ *   Krzysztof Kosiński <tweenk.pl@gmail.com>
+ * Copyright 2007-2011 Authors
  *
  * This library is free software; you can redistribute it and/or
  * modify it either under the terms of the GNU Lesser General Public
@@ -34,48 +37,41 @@
  *   MenTaLguY <mental@rydia.net>
  */
 
-#include <2geom/d2.h>
-
-#ifndef LIB2GEOM_RECT_H
-#define LIB2GEOM_RECT_H
+#ifndef LIB2GEOM_SEEN_RECT_H
+#define LIB2GEOM_SEEN_RECT_H
 
+#include <boost/optional.hpp>
 #include <2geom/affine.h>
-#include <boost/optional/optional.hpp>
+#include <2geom/interval.h>
+#include <2geom/int-rect.h>
 
 namespace Geom {
 
 /**
- * @brief Axis-aligned, non-empty rectangle - convenience typedef
+ * @brief Axis-aligned rectangle that can be empty.
  * @ingroup Primitives
  */
-typedef D2<Interval> Rect;
-class OptRect;
-
-inline Rect unify(Rect const &, Rect const &);
+typedef GenericOptRect<Coord> OptRect;
 
 /**
  * @brief Axis aligned, non-empty rectangle.
  * @ingroup Primitives
  */
-template<>
-class D2<Interval> {
-private:
-    Interval f[2];
+class Rect
+    : public GenericRect<Coord>
+    , boost::multipliable< Rect, Affine >
+{
+    typedef GenericRect<Coord> Base;
 public:
     /// @name Create rectangles.
     /// @{
     /** @brief Create a rectangle that contains only the point at (0,0). */
-    D2<Interval>() { f[X] = f[Y] = Interval(); }
+    Rect() {}
     /** @brief Create a rectangle from X and Y intervals. */
-    D2<Interval>(Interval const &a, Interval const &b) {
-        f[X] = a;
-        f[Y] = b;
-    }
+    Rect(Interval const &a, Interval const &b) : Base(a,b) {}
     /** @brief Create a rectangle from two points. */
-    D2<Interval>(Point const & a, Point const & b) {
-        f[X] = Interval(a[X], b[X]);
-        f[Y] = Interval(a[Y], b[Y]);
-    }
+    Rect(Point const &a, Point const &b) : Base(a,b) {}
+    Rect(Base const &b) : Base(b) {}
     /** @brief Create a rectangle from a range of points.
      * The resulting rectangle will contain all ponts from the range.
      * The return type of iterators must be convertible to Point.
@@ -85,12 +81,7 @@ public:
      * @return Rectangle that contains all points from [start, end). */
     template <typename InputIterator>
     static Rect from_range(InputIterator start, InputIterator end) {
-        assert(start != end);
-        Point p1 = *start++;
-        Rect result(p1, p1);
-        for (; start != end; ++start) {
-            result.expandTo(*start);
-        }
+        Rect result = Base::from_range(start, end);
         return result;
     }
     /** @brief Create a rectangle from a C-style array of points it should contain. */
@@ -102,260 +93,85 @@ public:
 
     /// @name Inspect dimensions.
     /// @{
-    Interval& operator[](unsigned i)              { return f[i]; }
-    Interval const & operator[](unsigned i) const { return f[i]; }
-
-    Point min() const { return Point(f[X].min(), f[Y].min()); }
-    Point max() const { return Point(f[X].max(), f[Y].max()); }
-    /** @brief Return the n-th corner of the rectangle.
-     * If the Y axis grows upwards, this returns corners in clockwise order
-     * starting from the lower left. If Y grows downwards, it returns the corners
-     * in counter-clockwise order starting from the upper left. */
-    Point corner(unsigned i) const {
-        switch(i % 4) {
-            case 0:  return Point(f[X].min(), f[Y].min());
-            case 1:  return Point(f[X].max(), f[Y].min());
-            case 2:  return Point(f[X].max(), f[Y].max());
-            default: return Point(f[X].min(), f[Y].max());
-        }
-    }
-        
-    //We should probably remove these - they're coord sys gnostic
-    /** @brief Return top coordinate of the rectangle (+Y is downwards). */
-    Coord top() const { return f[Y].min(); }
-    /** @brief Return bottom coordinate of the rectangle (+Y is downwards). */
-    Coord bottom() const { return f[Y].max(); }
-    /** @brief Return leftmost coordinate of the rectangle (+X is to the right). */
-    Coord left() const { return f[X].min(); }
-    /** @brief Return rightmost coordinate of the rectangle (+X is to the right). */
-    Coord right() const { return f[X].max(); }
-
-    Coord width() const { return f[X].extent(); }
-    Coord height() const { return f[Y].extent(); }
-
-    /** @brief Get rectangle's width and height as a point.
-     * @return Point with X coordinate corresponding to the width and the Y coordinate
-     *         corresponding to the height of the rectangle. */
-    Point dimensions() const { return Point(f[X].extent(), f[Y].extent()); }
-    Point midpoint() const { return Point(f[X].middle(), f[Y].middle()); }
-
-/**
- * \brief Compute the area of this rectangle.
- *
- * Note that a zero area rectangle is not empty - just as the interval [0,0] contains one point, the rectangle [0,0] x [0,0] contains 1 point and no area.
- * \retval For a valid return value, the rect must be tested for emptyness first.
- */
-    /** @brief Compute rectangle's area. */
-    Coord area() const { return f[X].extent() * f[Y].extent(); }
     /** @brief Check whether the rectangle has zero area up to specified tolerance.
      * @param eps Maximum value of the area to consider empty
      * @return True if rectangle has an area smaller than tolerance, false otherwise */
-    bool hasZeroArea(double eps = EPSILON) const { return (area() <= eps); }
-
-    /** @brief Get the larger extent (width or height) of the rectangle. */
-    Coord maxExtent() const { return std::max(f[X].extent(), f[Y].extent()); }
-    /** @brief Get the smaller extent (width or height) of the rectangle. */
-    Coord minExtent() const { return std::min(f[X].extent(), f[Y].extent()); }
+    bool hasZeroArea(Coord eps = EPSILON) const { return (area() <= eps); }
     /// @}
 
     /// @name Test other rectangles and points for inclusion.
     /// @{
-    /** @brief Check whether the rectangles have any common points. */
-    bool intersects(Rect const &r) const { 
-        return f[X].intersects(r[X]) && f[Y].intersects(r[Y]);
-    }
     /** @brief Check whether the interiors of the rectangles have any common points. */
     bool interiorIntersects(Rect const &r) const {
         return f[X].interiorIntersects(r[X]) && f[Y].interiorIntersects(r[Y]);
     }
-    /** @brief Check whether the rectangle includes all points in the given rectangle. */
-    bool contains(Rect const &r) const { 
-        return f[X].contains(r[X]) && f[Y].contains(r[Y]);
-    }
     /** @brief Check whether the interior includes all points in the given rectangle.
      * Interior of the rectangle is the entire rectangle without its borders. */
     bool interiorContains(Rect const &r) const { 
         return f[X].interiorContains(r[X]) && f[Y].interiorContains(r[Y]);
     }
-
-    /** @brief Check whether the rectangles have any common points.
-     * A non-empty rectangle will not intersect empty rectangles. */
-    inline bool intersects(OptRect const &r) const;
-    /** @brief Check whether the rectangle includes all points in the given rectangle.
-     * A non-empty rectangle will contain any empty rectangle. */
-    inline bool contains(OptRect const &r) const;
-    /** @brief Check whether the interior includes all points in the given rectangle.
-     * The interior of a non-empty rectangle will contain any empty rectangle. */
     inline bool interiorContains(OptRect const &r) const;
+    /// @}
 
-    /** @brief Check whether the given point is within the rectangle. */
-    bool contains(Point const &p) const {
-        return f[X].contains(p[X]) && f[Y].contains(p[Y]);
+    /// @name Rounding to integer coordinates
+    /// @{
+    /** @brief Return the smallest integer rectangle which contains this one. */
+    IntRect roundOutwards() const {
+        IntRect ir(f[X].roundOutwards(), f[Y].roundOutwards());
+        return ir;
     }
-    /** @brief Check whether the given point is in the rectangle's interior.
-     * This means the point must lie within the rectangle but not on its border. */
-    bool interiorContains(Point const &p) const {
-        return f[X].interiorContains(p[X]) && f[Y].interiorContains(p[Y]);
+    /** @brief Return the largest integer rectangle which is contained in this one. */
+    OptIntRect roundInwards() const {
+        OptIntRect oir(f[X].roundInwards(), f[Y].roundInwards());
+        return oir;
     }
     /// @}
 
-    /// @name Modify the rectangle.
+    /// @name Operators
     /// @{
-    /** @brief Enlarge the rectangle to contain the given point. */
-    void expandTo(Point p)        { 
-        f[X].expandTo(p[X]);  f[Y].expandTo(p[Y]);
-    }
-    /** @brief Enlarge the rectangle to contain the given rectangle. */
-    void unionWith(Rect const &b) { 
-        f[X].unionWith(b[X]); f[Y].unionWith(b[Y]);
-    }
-    /** @brief Enlarge the rectangle to contain the given rectangle.
-     * Unioning with an empty rectangle results in no changes. */
-    void unionWith(OptRect const &b);
-    
-    //TODO: figure out how these work with negative values and OptRect
-    /** @brief Expand the rectangle in both directions by the specified amount.
-     * Note that this is different from scaling. Negative values wil shrink the
-     * rectangle. If <code>-amount</code> is larger than
-     * half of the width, the X interval will contain only the X coordinate
-     * of the midpoint; same for height. */
-    void expandBy(Coord amount) {
-        f[X].expandBy(amount);  f[Y].expandBy(amount);
-    }
-    /** @brief Expand the rectangle by the coordinates of the given point.
-     * This will expand the width by the X coordinate of the point in both directions
-     * and the height by Y coordinate of the point. Negative coordinate values will
-     * shrink the rectangle. If <code>-p[X]</code> is larger than half of the width,
-     * the X interval will contain only the X coordinate of the midpoint; same for height. */
-    void expandBy(Point const p) { 
-        f[X].expandBy(p[X]);  f[Y].expandBy(p[Y]);
-    }
+    Rect &operator*=(Affine const &m);
     /// @}
 };
 
-inline Rect unify(Rect const & a, Rect const & b) {
-    return Rect(unify(a[X], b[X]), unify(a[Y], b[Y]));
+Coord distanceSq(Point const &p, Rect const &rect);
+Coord distance(Point const &p, Rect const &rect);
+
+inline bool Rect::interiorContains(OptRect const &r) const {
+    return !r || interiorContains(static_cast<Rect const &>(*r));
 }
 
-inline Rect union_list(std::vector<Rect> const &r) {
-    if(r.empty()) return Rect(Interval(0,0), Interval(0,0));
-    Rect ret = r[0];
-    for(unsigned i = 1; i < r.size(); i++)
-        ret.unionWith(r[i]);
+// the functions below do not work when defined generically
+inline OptRect operator&(Rect const &a, Rect const &b) {
+    OptRect ret(a);
+    ret.intersectWith(b);
     return ret;
 }
-
-inline
-Coord distanceSq( Point const& p, Rect const& rect )
-{
-    double dx = 0, dy = 0;
-    if ( p[X] < rect.left() )
-    {
-        dx = p[X] - rect.left();
-    }
-    else if ( p[X] > rect.right() )
-    {
-        dx = rect.right() - p[X];
-    }
-    if ( p[Y] < rect.top() )
-    {
-        dy = rect.top() - p[Y];
-    }
-    else if (  p[Y] > rect.bottom() )
-    {
-        dy = p[Y] - rect.bottom();
-    }
-    return dx*dx + dy*dy;
+inline OptRect intersect(Rect const &a, Rect const &b) {
+    return a & b;
 }
-
-/**
- * Returns the smallest distance between p and rect.
- */
-inline 
-Coord distance( Point const& p, Rect const& rect )
-{
-    return std::sqrt(distanceSq(p, rect));
+inline OptRect intersect(OptRect const &a, OptRect const &b) {
+    return a & b;
 }
-
-/**
- * @brief Axis-aligned rectangle that can be empty.
- * @ingroup Primitives
- */
-class OptRect : public boost::optional<Rect> {
-public:
-    OptRect() : boost::optional<Rect>() {};
-    OptRect(Rect const &a) : boost::optional<Rect>(a) {};
-
-    /**
-     * Creates an empty OptRect when one of the argument intervals is empty.
-     */
-    OptRect(OptInterval const &x_int, OptInterval const &y_int) {
-        if (x_int && y_int) {
-            *this = Rect(*x_int, *y_int);
-        }
-        // else, stay empty.
-    }
-
-    /** @brief Check for emptiness. */
-    inline bool isEmpty() const { return (*this == false); };
-
-    bool intersects(Rect const &r) const { return r.intersects(*this); }
-    bool contains(Rect const &r) const { return *this && (*this)->contains(r); }
-    bool interiorContains(Rect const &r) const { return *this && (*this)->interiorContains(r); }
-
-    bool intersects(OptRect const &r) const { return *this && (*this)->intersects(r); }
-    bool contains(OptRect const &r) const { return *this && (*this)->contains(r); }
-    bool interiorContains(OptRect const &r) const { return *this && (*this)->interiorContains(r); }
-
-    bool contains(Point const &p) const { return *this && (*this)->contains(p); }
-    bool interiorContains(Point const &p) const { return *this && (*this)->contains(p); }
-
-    inline void unionWith(OptRect const &b) {
-        if (*this) { // check that we are not empty
-            (*this)->unionWith(b);
-        } else {
-            *this = b;
-        }
-    }
-};
-
-
-/** 
- * Returns the smallest rectangle that encloses both rectangles.
- * An empty argument is assumed to be an empty rectangle
- */
-inline OptRect unify(OptRect const & a, OptRect const & b) {
-    if (!a) {
-        return b;
-    } else if (!b) {
-        return a;
-    } else {
-        return unify(*a, *b);
-    }
+inline Rect unify(Rect const &a, Rect const &b) {
+    return a | b;
 }
-
-inline OptRect intersect(Rect const & a, Rect const & b) {
-    return OptRect(intersect(a[X], b[X]), intersect(a[Y], b[Y]));
+inline OptRect unify(OptRect const &a, OptRect const &b) {
+    return a | b;
 }
 
-inline void Rect::unionWith(OptRect const &b) { 
-    if (b) {
-        unionWith(*b);
-    }
-}
-inline bool Rect::intersects(OptRect const &r) const {
-    return r && intersects(*r);
-}
-inline bool Rect::contains(OptRect const &r) const {
-    return !r || contains(*r);
-}
-inline bool Rect::interiorContains(OptRect const &r) const {
-    return !r || interiorContains(*r);
+/** @brief Union a list of rectangles
+ * @deprecated Use OptRect::from_range instead */
+inline Rect union_list(std::vector<Rect> const &r) {
+    if(r.empty()) return Rect(Interval(0,0), Interval(0,0));
+    Rect ret = r[0];
+    for(unsigned i = 1; i < r.size(); i++)
+        ret.unionWith(r[i]);
+    return ret;
 }
 
 } // end namespace Geom
 
-#endif //_2GEOM_RECT
+#endif // LIB2GEOM_SEEN_RECT_H
 
 /*
   Local Variables:
diff --git a/src/2geom/region.h b/src/2geom/region.h
index e23d6a158..06a4f63e9 100644
--- a/src/2geom/region.h
+++ b/src/2geom/region.h
@@ -1,7 +1,7 @@
 /**
  * \file
- * \brief  \todo brief description
- *
+ * \brief Uncrossed path for boolean algorithms
+ *//*
  * Authors:
  *      ? <?@?.?>
  * 
diff --git a/src/2geom/sbasis-2d.h b/src/2geom/sbasis-2d.h
index f1218b028..00429e259 100644
--- a/src/2geom/sbasis-2d.h
+++ b/src/2geom/sbasis-2d.h
@@ -1,7 +1,7 @@
 /**
  * \file
- * \brief  \todo brief description
- *
+ * \brief Obsolete 2D SBasis function class
+ *//*
  * Authors:
  *      Nathan Hurst <?@?.?>
  *      JFBarraud <?@?.?>
diff --git a/src/2geom/sbasis-curve.h b/src/2geom/sbasis-curve.h
index 22fe4fc1f..554b702e6 100644
--- a/src/2geom/sbasis-curve.h
+++ b/src/2geom/sbasis-curve.h
@@ -33,8 +33,8 @@
  * the specific language governing rights and limitations.
  */
 
-#ifndef _2GEOM_SBASIS_CURVE_H_
-#define _2GEOM_SBASIS_CURVE_H_
+#ifndef LIB2GEOM_SEEN_SBASIS_CURVE_H
+#define LIB2GEOM_SEEN_SBASIS_CURVE_H
 
 #include <2geom/curve.h>
 #include <2geom/nearest-point.h>
@@ -45,24 +45,45 @@ namespace Geom
 
 /** @brief Symmetric power basis curve.
  *
- * Symmetric power basis (S-basis for short) polynomials are a versatile numeric representation
- * of arbitrary continuous curves. They combine the properties of Bezier curves
- * (geometric interpretation of parameters, numerical stability near ends of the curve)
- * and the monomial basis (fast evaluation). They are the main representation of curves
+ * Symmetric power basis (S-basis for short) polynomials are a versatile numeric
+ * representation of arbitrary continuous curves. They are the main representation of curves
  * in 2Geom.
  *
+ * S-basis is defined for odd degrees and composed of the following polynomials:
+ * \f{align*}{
+         P_k^0(t) &= t^k (1-t)^{k+1} \\
+         P_k^1(t) &= t^{k+1} (1-t)^k \f}
+ * This can be understood more easily with the help of the chart below. Each square
+ * represents a product of a specific number of \f$t\f$ and \f$(1-t)\f$ terms. Red dots
+ * are the canonical (monomial) basis, the green dots are the Bezier basis, and the blue
+ * dots are the S-basis, all of them of degree 7.
+ *
+ * @image html sbasis.png "Illustration of the monomial, Bezier and symmetric power bases"
+ *
+ * The S-Basis has several important properties:
+ * - S-basis polynomials are closed under multiplication.
+ * - Evaluation is fast, using a modified Horner scheme.
+ * - Degree change is as trivial as in the monomial basis. To elevate, just add extra
+ *   zero coefficients. To reduce the degree, truncate the terms in the highest powers.
+ *   Compare this with Bezier curves, where degree change is complicated.
+ * - Conversion between S-basis and Bezier basis is numerically stable.
+ *
+ * More in-depth information can be found in the following paper:
+ * J Sanchez-Reyes, "The symmetric analogue of the polynomial power basis".
+ * ACM Transactions on Graphics, Vol. 16, No. 3, July 1997, pages 319--357.
+ * http://portal.acm.org/citation.cfm?id=256162
+ *
  * @ingroup Curves
  */
 class SBasisCurve : public Curve {
-
 private:
-    SBasisCurve();
     D2<SBasis> inner;
   
 public:
     explicit SBasisCurve(D2<SBasis> const &sb) : inner(sb) {}
     explicit SBasisCurve(Curve const &other) : inner(other.toSBasis()) {}
 
+#ifndef DOXYGEN_SHOULD_SKIP_THIS
     virtual Curve *duplicate() const { return new SBasisCurve(*this); }
     virtual Point initialPoint() const    { return inner.at0(); }
     virtual Point finalPoint() const      { return inner.at1(); }
@@ -104,16 +125,12 @@ public:
     virtual int degreesOfFreedom() const {
         return inner[0].degreesOfFreedom() + inner[1].degreesOfFreedom();
     }
+#endif
 };
 
-
 } // end namespace Geom
 
-
-#endif // _2GEOM_SBASIS_CURVE_H_
-
-
-
+#endif // LIB2GEOM_SEEN_SBASIS_CURVE_H
 
 /*
   Local Variables:
diff --git a/src/2geom/sbasis-to-bezier.h b/src/2geom/sbasis-to-bezier.h
index 5b88a40fa..819aa87d6 100644
--- a/src/2geom/sbasis-to-bezier.h
+++ b/src/2geom/sbasis-to-bezier.h
@@ -1,7 +1,7 @@
 /**
  * \file
- * \brief  \todo brief description
- *
+ * \brief Conversion between SBasis and Bezier basis polynomials
+ *//*
  * Authors:
  *      ? <?@?.?>
  * 
diff --git a/src/2geom/sbasis.cpp b/src/2geom/sbasis.cpp
index e313ad08d..89640af5c 100644
--- a/src/2geom/sbasis.cpp
+++ b/src/2geom/sbasis.cpp
@@ -34,7 +34,7 @@
 #include <cmath>
 
 #include <2geom/sbasis.h>
-#include <2geom/isnan.h>
+#include <2geom/math-utils.h>
 
 namespace Geom{
 
diff --git a/src/2geom/sbasis.h b/src/2geom/sbasis.h
index b1b0b6c2a..7a7e33fe4 100644
--- a/src/2geom/sbasis.h
+++ b/src/2geom/sbasis.h
@@ -345,7 +345,7 @@ SBasis compose_inverse(SBasis const &f, SBasis const &g, unsigned order=2, doubl
  \relates SBasis
 */
 inline SBasis portion(const SBasis &t, double from, double to) { return compose(t, Linear(from, to)); }
-inline SBasis portion(const SBasis &t, Interval ivl) { return compose(t, Linear(ivl[0], ivl[1])); }
+inline SBasis portion(const SBasis &t, Interval ivl) { return compose(t, Linear(ivl.min(), ivl.max())); }
 
 // compute f(g)
 inline SBasis
diff --git a/src/2geom/sturm.h b/src/2geom/sturm.h
deleted file mode 100644
index 4fef1b954..000000000
--- a/src/2geom/sturm.h
+++ /dev/null
@@ -1,70 +0,0 @@
-#ifndef LIB2GEOM_STURM_HEADER
-#define LIB2GEOM_STURM_HEADER
-
-#include <2geom/poly.h>
-#include <2geom/utils.h>
-
-namespace Geom {
-
-class sturm : public std::vector<Poly>{
-public:
-    sturm(Poly const &X) {
-        push_back(X);
-        push_back(derivative(X));
-        Poly Xi = back();
-        Poly Xim1 = X;
-        std::cout << "sturm:\n" << Xim1 << std::endl;
-        std::cout << Xi << std::endl;
-        while(Xi.size() > 1) {
-            Poly r;
-            divide(Xim1, Xi, r);
-            std::cout << r << std::endl;
-            assert(r.size() < Xi.size());
-            Xim1 = Xi;
-            Xi = -r;
-            assert(Xim1.size() > Xi.size());
-            push_back(Xi);
-        }
-    }
-    
-    unsigned count_signs(double t) {
-        unsigned n_signs = 0;/*  Number of sign-changes */
-        const double big = 1e20; // a number such that practical polys would overflow on evaluation
-        if(t >= big) {
-            int old_sign = sgn((*this)[0].back());
-            for (unsigned i = 1; i < size(); i++) {
-                int sign = sgn((*this)[i].back());
-                if (sign != old_sign)
-                    n_signs++;
-                old_sign = sign;
-            }
-        } else {
-            int old_sign = sgn((*this)[0].eval(t));
-            for (unsigned i = 1; i < size(); i++) {
-                int sign = sgn((*this)[i].eval(t));
-                if (sign != old_sign)
-                    n_signs++;
-                old_sign = sign;
-            }
-        }
-        return n_signs;
-    }
-    
-    unsigned n_roots_between(double l, double r) {
-        return count_signs(l) - count_signs(r);
-    }
-};
-
-} //namespace Geom
-
-#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/svg-elliptical-arc.h b/src/2geom/svg-elliptical-arc.h
index 79497cdb3..ba0a18257 100644
--- a/src/2geom/svg-elliptical-arc.h
+++ b/src/2geom/svg-elliptical-arc.h
@@ -1,7 +1,6 @@
 /**
  * \file
  * \brief SVG 1.1-compliant elliptical arc curve
- *
  *//*
  * Authors:
  *    MenTaLguY <mental@rydia.net>
@@ -35,8 +34,8 @@
  */
 
 
-#ifndef _2GEOM_SVG_ELLIPTICAL_ARC_H_
-#define _2GEOM_SVG_ELLIPTICAL_ARC_H_
+#ifndef LIB2GEOM_SEEN_SVG_ELLIPTICAL_ARC_H
+#define LIB2GEOM_SEEN_SVG_ELLIPTICAL_ARC_H
 
 #include <2geom/curve.h>
 #include <2geom/angle.h>
@@ -49,8 +48,7 @@
 #include <2geom/numeric/fitting-model.h>
 #include <algorithm>
 
-namespace Geom
-{
+namespace Geom {
 
 class SVGEllipticalArc : public EllipticalArc {
 public:
@@ -267,13 +265,9 @@ class make_elliptical_arc
       bool svg_compliant;
 };
 
-
 } // end namespace Geom
 
-
-
-
-#endif /* _2GEOM_SVG_ELLIPTICAL_ARC_H_ */
+#endif // LIB2GEOM_SEEN_SVG_ELLIPTICAL_ARC_H
 
 /*
   Local Variables:
diff --git a/src/2geom/sweep.h b/src/2geom/sweep.h
index 1c73efee0..91371e6fb 100644
--- a/src/2geom/sweep.h
+++ b/src/2geom/sweep.h
@@ -1,7 +1,7 @@
 /**
  * \file
- * \brief  \todo brief description
- *
+ * \brief Sweepline intersection of groups of rectangles
+ *//*
  * Authors:
  *      ? <?@?.?>
  * 
diff --git a/src/2geom/toposweep.cpp b/src/2geom/toposweep.cpp
new file mode 100644
index 000000000..cfb91857c
--- /dev/null
+++ b/src/2geom/toposweep.cpp
@@ -0,0 +1,663 @@
+#include <2geom/toposweep.h>
+
+#include <2geom/path-intersection.h>
+#include <2geom/basic-intersection.h>
+
+//using namespace Geom;
+
+namespace Geom {
+
+TopoGraph::Edge &TopoGraph::Vertex::operator[](unsigned ix) {
+    ix %= degree();
+    return ix < enters.size() ? enters[ix] : exits[ix - enters.size()];
+}
+
+TopoGraph::Edge TopoGraph::Vertex::operator[](unsigned ix) const {
+    ix %= degree();
+    return ix < enters.size() ? enters[ix] : exits[ix - enters.size()];
+}
+
+void TopoGraph::Vertex::erase(unsigned ix) {
+    ix %= degree();
+    if(ix < enters.size())
+        enters.erase(enters.begin() + ix);
+    else
+        exits.erase(exits.begin() + (ix - enters.size()));
+}
+
+void TopoGraph::Vertex::insert(unsigned ix, Edge v) {
+    ix %= degree();
+    if(ix < enters.size())
+        enters.insert(enters.begin() + ix, v);
+    else
+        exits.insert(exits.begin() + (ix - enters.size()), v);
+}
+
+unsigned TopoGraph::Vertex::find_section(boost::shared_ptr<Section> section) const {
+    unsigned i = 0;
+    for(; i < degree(); i++)
+        if((*this)[i].section == section) return i;
+    return i;
+}
+
+TopoGraph::Edge TopoGraph::remove_edge(unsigned ix, unsigned jx) {
+    Vertex &v = vertices[ix];
+    if(v.degree()) {
+        jx %= v.degree();
+        Edge &ret = v[jx];
+        v.erase(jx);
+        v = vertices[ret.other];
+        if(v.degree()) {
+            v.erase(v.find_section(ret.section));
+            return ret;
+        }
+    }
+    assert(0);
+}
+
+void TopoGraph::cannonize() {
+    std::vector<unsigned> vix;
+    unsigned ix = 0;
+    for(unsigned i = 0; i < vertices.size(); i++) {
+        vix.push_back(ix);
+        if(vertices[i].degree() != 0) vertices[ix++] = vertices[i];
+    }
+    
+    for(unsigned i = 0; i < ix; i++)
+        for(unsigned j = 0; j < vertices[i].degree(); j++)
+            vertices[i][j].other = vix[vertices[i][j].other];
+}
+
+
+void TopoGraph::assert_invariants() const {
+    for(unsigned i = 0; i < vertices.size(); i++) {
+        for(unsigned j = 0; j < vertices[i].degree(); j++) {
+            Edge e = vertices[i][j];
+            assert(e.other != i);
+            assert(are_near(e.section->fp, vertices[i].avg, tol) || are_near(e.section->tp, vertices[i].avg, tol));
+            assert(!are_near(e.section->fp, e.section->tp, tol));
+            assert(e.section.get());
+            unsigned oix = vertices[e.other].find_section(e.section);
+            assert(oix != vertices[e.other].degree());
+        }
+    }
+}
+
+//near predicate utilized in process_splits
+template<typename T>
+struct NearPredicate { bool operator()(T x, T y) { return are_near(x, y); } };
+
+// ensures that f and t are elements of a vector, sorts and uniqueifies
+// also asserts that no values fall outside of f and t
+// if f is greater than t, the sort is in reverse
+void process_splits(std::vector<double> &splits, double f, double t) {
+    splits.push_back(f);
+    std::sort(splits.begin(), splits.end());
+    while(are_near(splits.back(), t)) splits.erase(splits.end() - 1);
+    splits.push_back(t);
+    if(f > t) std::reverse(splits.begin(), splits.end());
+
+    //remove any splits which fall outside t / f
+    while(!splits.empty() && splits.front() != f) splits.erase(splits.begin());
+    while(!splits.empty() && splits.back() != t) splits.erase(splits.end() - 1);
+    
+    std::vector<double>::iterator end = std::unique(splits.begin(), splits.end(), NearPredicate<double>());
+    splits.resize(end - splits.begin());
+}
+
+// A little sugar for appending a list to another
+template<typename T>
+void concatenate(T &a, T const &b) { a.insert(a.end(), b.begin(), b.end()); }
+
+//returns a list of monotonic sections of a path
+//TODO: handle saddle points
+std::vector<boost::shared_ptr<Section> > mono_sections(PathVector const &ps, Dim2 d) {
+    std::vector<boost::shared_ptr<Section> > monos;
+    for(unsigned i = 0; i < ps.size(); i++) {
+        //TODO: necessary? can we have empty paths?
+        if(ps[i].size()) {
+            for(unsigned j = 0; j < ps[i].size(); j++) {
+                //find the points of 0 derivative
+                Curve* deriv = ps[i][j].derivative();
+                std::vector<double> splits = deriv->roots(0, X);
+                concatenate(splits, deriv->roots(0, Y));
+                delete deriv;
+                process_splits(splits, 0, 1);
+                //split on points of 0 derivative
+                for(unsigned k = 1; k < splits.size(); k++)
+                    monos.push_back(boost::shared_ptr<Section>(new Section(CurveIx(i,j), splits[k-1], splits[k], ps, d)));
+            }
+        }
+    }
+    return monos;
+}
+
+//finds the t-value on a section, which corresponds to a particular horizontal or vertical line
+//d indicates the dimension along which the roots is performed.
+//-1 is returned if no root is found
+double section_root(Section const &s, PathVector const &ps, double v, Dim2 d) {
+    std::vector<double> roots = s.curve.get(ps).roots(v, d);
+    for(unsigned j = 0; j < roots.size(); j++)
+        if(Interval(s.f, s.t).contains(roots[j])) return roots[j];
+    return -1;
+}
+
+bool SectionSorter::section_order(Section const &a, double at, Section const &b, double bt) const {
+    Point ap = a.curve.get(ps).pointAt(at);
+    Point bp = b.curve.get(ps).pointAt(bt);
+    if(are_near(ap[dim], bp[dim], tol)) {
+        // since the sections are monotonic, if the endpoints are on opposite sides of this
+        // coincidence, the order is determinable
+        if(a.tp[dim] < ap[dim] && b.tp[dim] > bp[dim]) return true;
+        if(a.tp[dim] > ap[dim] && b.tp[dim] < bp[dim]) return false;
+        //TODO: sampling / higher derivatives when unit tangents match
+        Point ad = a.curve.get(ps).unitTangentAt(a.f);
+        Point bd = b.curve.get(ps).unitTangentAt(b.f);
+        // tangent can point backwards
+        if(ad[1-dim] < 0) ad = -ad;
+        if(bd[1-dim] < 0) bd = -bd;
+        return ad[dim] < bd[dim];
+    }
+    return ap[dim] < bp[dim];
+}
+
+bool SectionSorter::operator()(Section const &a, Section const &b) const {
+    if(&a == &b) return false;
+    Rect ra = a.bbox(), rb = b.bbox();
+    //TODO: should we use tol in these conditions?
+    if(ra[dim].max() <= rb[dim].min()) return true;
+    if(rb[dim].max() <= ra[dim].min()) return false;
+    //we know that the rects intersect on dim
+    //by referencing f / t we are assuming that the section was constructed with 1-dim
+    if(ra[1-dim].intersects(rb[1-dim])) {
+        if(are_near(a.fp[1-dim], b.fp[1-dim], tol)) {
+            return section_order(a, a.f > a.t ? a.f - 0.01 : a.f + 0.01,
+                                 b, b.f > b.t ? b.f - 0.01 : b.f + 0.01);
+        } else if(a.fp[1-dim] < b.fp[1-dim]) {
+            //b inside a
+            double ta = section_root(a, ps, b.fp[1-dim], Dim2(1-dim));
+            //TODO: fix bug that necessitates this
+            if(ta == -1) ta = (a.t + a.f) / 2;
+            return section_order(a, ta, b, b.f);
+        } else {
+            //a inside b
+            double tb = section_root(b, ps, a.fp[1-dim], Dim2(1-dim));
+            //TODO: fix bug that necessitates this
+            if(tb == -1) tb = (b.t + b.f) / 2;
+            return section_order(a, a.f, b, tb);
+        }
+    }
+    
+    return Point::LexOrderRt(dim)(a.fp, b.fp);
+}
+
+// splits a section into pieces, as specified by an array of doubles, mutating the section to
+// represent the first part, and returning the rest
+//TODO: output iterator?
+std::vector<boost::shared_ptr<Section> > split_section(boost::shared_ptr<Section> s, PathVector const &ps, std::vector<double> &cuts, Dim2 d) {
+    std::vector<boost::shared_ptr<Section> > ret;
+    
+    process_splits(cuts, s->f, s->t);
+    if(cuts.size() <= 2) return ret;
+    
+    s->t = cuts[1];
+    s->tp = s->curve.get(ps)(cuts[1]);
+    assert(Point::LexOrderRt(d)(s->fp, s->tp));
+    
+    ret.reserve(cuts.size() - 2);
+    for(int i = cuts.size() - 1; i > 1; i--) ret.push_back(boost::shared_ptr<Section>(new Section(s->curve, cuts[i-1], cuts[i], ps, d)));
+    return ret;
+}
+
+//merges the sorted lists a and b according to comparison z
+template<typename X, typename Z>
+void merge(X &a, X const &b, Z const &z) {
+     a.reserve(a.size() + b.size());
+     unsigned start = a.size();
+     concatenate(a, b);
+     std::inplace_merge(a.begin(), a.begin() + start, a.end(), z);
+}
+
+//TODO: faster than linear
+unsigned find_vertex(std::vector<TopoGraph::Vertex> const &vertices, Point p, double tol) {
+    for(unsigned i = 0; i < vertices.size(); i++)
+        if(are_near(vertices[i].avg, p, tol)) return i;
+    return vertices.size();
+}
+
+//takes a vector of T pointers, and returns a vector of T with copies
+template<typename T>
+std::vector<T> deref_vector(std::vector<boost::shared_ptr<T> > const &xs, unsigned from = 0) {
+    std::vector<T> ret;
+    ret.reserve(xs.size() - from);
+    for(unsigned i = from; i < xs.size(); i++)
+        ret.push_back(T(*xs[i]));
+    return ret;
+}
+
+//used to create reversed sorting predicates
+template<typename C>
+struct ReverseAdapter {
+    typedef typename C::second_argument_type first_argument_type;
+    typedef typename C::first_argument_type second_argument_type;
+    typedef typename C::result_type result_type;
+    const C &comp;
+    ReverseAdapter(const C &c) : comp(c) {}
+    result_type operator()(const first_argument_type &a, const second_argument_type &b) const { return comp(b, a); }
+};
+
+//used to sort std::vector<Section*>
+template<typename C>
+struct DerefAdapter {
+    typedef typename boost::shared_ptr<typename C::first_argument_type> first_argument_type;
+    typedef typename boost::shared_ptr<typename C::second_argument_type> second_argument_type;
+    typedef typename C::result_type result_type;
+    const C &comp;
+    DerefAdapter(const C &c) : comp(c) {}
+    result_type operator()(const first_argument_type a, const second_argument_type b) const {
+        if(!a) return false;
+        if(!b) return true;
+        return comp(*a, *b);
+    }
+};
+
+struct EdgeSorter {
+    typedef TopoGraph::Edge first_argument_type;
+    typedef TopoGraph::Edge second_argument_type;
+    typedef bool result_type;
+    SectionSorter s;
+    EdgeSorter(const PathVector &rs, Dim2 d, double t) : s(rs, d, t) {}
+    bool operator()(TopoGraph::Edge const &e1, TopoGraph::Edge const &e2) const { return s(*e1.section, *e2.section); }
+};
+
+#ifdef SWEEP_GRAPH_DEBUG
+//used for debugging purposes - each element represents a subsequent iteration of the algorithm.
+std::vector<std::vector<Section> > monoss;
+std::vector<std::vector<Section> > chopss;
+std::vector<std::vector<Section> > contexts;
+#endif
+
+/*
+ 1) take item off sweep sorted todo
+ 2) find all of the to-values before the beginning of this section
+ 3) sort these lexicographically, process them in order, grouping other sections in the context, and constructing a vertex in one fell swoop.
+ 4) add our section into context, splitting on intersections
+ 
+ 3 is novel, we perform it by storing 
+ */
+
+template<typename A, typename B, typename Z>
+struct MergeIterator {
+    A const &a;
+    B &b;
+    Z const &z;
+    unsigned ai;
+    bool on_a;
+    MergeIterator(A const &av, B &bv, Z const &zv) : a(av), b(bv), z(zv), ai(0), on_a(b.empty() || z(a[0], b.back())) {}
+    MergeIterator &operator++() {
+        if(!done()) {
+            on_a = b.empty() ? true : (ai >= a.size() ? false : z(a[ai], b.back()));
+            if(on_a) {
+                ++ai;
+                if(ai >= a.size()) on_a = false;
+            } else {
+                b.erase(b.end());
+                if(b.empty()) on_a = true;
+            }
+        }
+        return *this;
+    }
+    typename A::value_type operator*() {
+        assert(!done());
+        return on_a ? a[ai] : b.back();
+    }
+    bool done() { return b.empty() && ai >= a.size() - 1; }
+    typename A::value_type operator->() { assert(!done()); return on_a ? a[ai] : b.back(); }
+};
+
+void modify_windings(std::vector<int> &windings, boost::shared_ptr<Section> sec, Dim2 d) {
+    unsigned k = sec->curve.path;
+    if(k >= windings.size() || sec->fp[d] == sec->tp[d]) return;
+    if(sec->f < sec->t) windings[k]++;
+    if(sec->f > sec->t) windings[k]--;
+}
+
+struct Context {
+    boost::shared_ptr<Section> section;
+    int from_vert;
+    int to_vert;
+    Context(boost::shared_ptr<Section> sect, int from) : section(sect), from_vert(from), to_vert(-1) {}
+};
+
+template<typename C>
+struct ContextAdapter {
+    typedef Context first_argument_type;
+    typedef typename C::second_argument_type second_argument_type;
+    typedef typename C::result_type result_type;
+    const C &comp;
+    ContextAdapter(const C &c) : comp(c) {}
+    result_type operator()(const Context &a, const second_argument_type &b) const { return comp(a.section, b); }
+};
+
+#define DINF std::numeric_limits<double>::infinity()
+
+TopoGraph::TopoGraph(PathVector const &ps, Dim2 d, double t) : dim(d), tol(t) {
+    //s_sort = vertical section order
+    ContextAdapter<DerefAdapter<SectionSorter> > s_sort = DerefAdapter<SectionSorter>(SectionSorter(ps, (Dim2)(1-d), tol));
+    //sweep_sort = horizontal sweep order
+    DerefAdapter<SweepSorter> sweep_sort = DerefAdapter<SweepSorter>(SweepSorter(d));
+    //heap_sort = reverse horizontal sweep order
+    ReverseAdapter<DerefAdapter<SweepSorter> > heap_sort = ReverseAdapter<DerefAdapter<SweepSorter> >(sweep_sort);
+    //edge_sort = sorter for edges
+    EdgeSorter edge_sort = EdgeSorter(ps, (Dim2)(1-d), tol);
+    
+    std::vector<boost::shared_ptr<Section> > input_sections = mono_sections(ps, d), chops;
+    std::sort(input_sections.begin(), input_sections.end(), sweep_sort);
+    
+    std::vector<Context> context;
+    
+    vertices.reserve(input_sections.size());
+    
+    //std::vector<unsigned> to_process;
+    
+    std::vector<int> windings(ps.size(), 0);
+    for(MergeIterator<Area, Area, DerefAdapter<SweepSorter> > iter(input_sections, chops, sweep_sort); ; ++iter) {
+        //represents our position in the sweep, which controls what we finalize
+        //if we have no more to process, finish the rest by setting our position to infinity
+        Point lim;
+        if(iter.done()) lim[X] = lim[Y] = DINF; else lim = iter->fp;
+        
+        /*
+        //finalize vertices
+        for(unsigned i = 0; i < to_process.size(); i++)  {
+            if(vertices[to_process[i]].avg[d] + tol < lim[d])
+            for(unsigned j = 0; j < context.size(); j++) {
+                
+            }
+        } */
+        
+        //find all sections to remove
+        for(int i = context.size() - 1; i >= 0; i--) {
+            boost::shared_ptr<Section> sec = context[i].section;
+            if(Point::LexOrderRt(d)(lim, sec->tp)) {
+                //sec->tp is less than or equal to lim
+                if(context[i].to_vert == -1) {
+                    //we need to create a new vertex; add everything that enters it
+                    //Point avg;
+                    //unsigned cnt;
+                    std::vector<Edge> enters;
+                    std::fill(windings.begin(), windings.end(), 0);
+                    for(unsigned j = 0; j < context.size(); j++) {
+                        modify_windings(windings, context[j].section, d);
+                        if(are_near(sec->tp, context[j].section->tp, tol)) {
+                            assert(-1 == context[j].to_vert);
+                            context[j].section->windings = windings;
+                            context[j].to_vert = vertices.size();
+                            enters.push_back(Edge(context[j].section, context[j].from_vert));
+                            //avg += context[j].section->tp;
+                            //cnt++;
+                        }
+                    }
+                    //Vertex &v(avg / (double)cnt);
+                    Vertex v(context[i].section->tp);
+                    v.enters = enters;
+                    vertices.push_back(v);
+                    //to_process.push_back(vertices.size() - 1);
+                }
+                context.erase(context.begin() + i);
+            }
+        }
+        
+        if(!iter.done()) {
+            boost::shared_ptr<Section> s = *iter;
+            
+            //create a new context, associate a beginning vertex, and insert it in the proper location
+            unsigned ix = find_vertex(vertices, s->fp, tol);
+            if(ix == vertices.size()) {
+                vertices.push_back(Vertex(s->fp));
+                //to_process.push_back(vertices.size() - 1);
+            }
+            unsigned context_ix = std::lower_bound(context.begin(), context.end(), s, s_sort) - context.begin();
+
+            context.insert(context.begin() + context_ix, Context(s, ix));
+            
+            Interval si = Interval(s->fp[1-d], s->tp[1-d]);
+            
+            // Now we intersect with neighbors - do a sweep!
+            std::vector<double> this_splits;
+            for(unsigned i = 0; i < context.size(); i++) {
+                if(context[i].section == context[context_ix].section) continue;
+                
+                boost::shared_ptr<Section> sec = context[i].section;
+                
+                if(!si.intersects(Interval(sec->fp[1-d], sec->tp[1-d]))) continue;
+                
+                std::vector<double> other_splits;
+                Crossings xs = mono_intersect(s->curve.get(ps), Interval(s->f, s->t),
+                                              sec->curve.get(ps), Interval(sec->f, sec->t));
+                if(xs.empty()) continue;
+                
+                for(unsigned j = 0; j < xs.size(); j++) {
+                    this_splits.push_back(xs[j].ta);
+                    other_splits.push_back(xs[j].tb);
+                }
+                merge(chops, split_section(sec, ps, other_splits, d), heap_sort);
+            }
+            if(!this_splits.empty())
+                merge(chops, split_section(context[context_ix].section, ps, this_splits, d), heap_sort);
+            
+            std::sort(chops.begin(), chops.end(), heap_sort);
+            
+            if(context[context_ix].section->tp[d] - context[context_ix].section->fp[d] <= tol) {
+                if(!are_near(context[context_ix].section->tp, context[context_ix].section->fp, tol)) {
+                    ix = find_vertex(vertices, context[context_ix].section->tp, tol);
+                    if(ix != vertices.size()) {
+                        boost::shared_ptr<Section> sec = context[context_ix].section;
+                        Edge e(sec, context[context_ix].from_vert);
+                        
+                        std::vector<Edge>::iterator it = std::lower_bound(vertices[ix].enters.begin(), vertices[ix].enters.end(), e, edge_sort);
+                        
+                        if(vertices[ix].enters.empty()) {                        
+                            std::fill(windings.begin(), windings.end(), 0);
+                            for(unsigned j = 0; j <= context_ix; j++) modify_windings(windings, context[j].section, d);
+                        } else if(it == vertices[ix].enters.end()) {
+                            windings = (it-1)->section->windings;
+                            modify_windings(windings, (it-1)->section, d);
+                        } else {
+                            windings = it->section->windings;
+                        }
+                                   
+                        sec->windings = windings;
+                        modify_windings(windings, sec, d);
+                        
+                        for(std::vector<Edge>::iterator it2 = it; it2 != vertices[ix].enters.end(); ++it2) {
+                            it2->section->windings = windings;
+                            modify_windings(windings, it2->section, d);
+                        }
+                        
+                        vertices[ix].enters.insert(it, e);
+                        context.erase(context.begin() + context_ix);
+                    }
+                } else context.erase(context.begin() + context_ix);
+            }
+        }
+        
+        #ifdef SWEEP_GRAPH_DEBUG
+        std::vector<Section> rem;
+        for(unsigned i = iter.ai + 1; i < iter.a.size(); i++) rem.push_back(*iter.a[i]);
+        monoss.push_back(rem);
+        chopss.push_back(deref_vector(iter.b));
+        rem.clear();
+        for(unsigned i = 0; i < context.size(); i++) rem.push_back(*context[i].section);
+        contexts.push_back(rem);
+        #endif
+        
+        if(iter.done() && context.empty()) return;
+    }
+}
+
+void trim_whiskers(TopoGraph &g) {
+    std::vector<unsigned> affected;
+    
+    for(unsigned i = 0; i < g.size(); i++)
+        if(g[i].degree() == 1) affected.push_back(i);
+    
+    while(!affected.empty()) {
+        unsigned j = 0;
+        for(unsigned i = 0; i < affected.size(); i++)
+            if(g[affected[i]].degree() == 1)
+                affected[j++] = g.remove_edge(affected[i], 0).other;
+        affected.resize(j);
+    }
+}
+
+void add_edge_at(TopoGraph &g, unsigned ix, boost::shared_ptr<Section> s, TopoGraph::Edge jx, bool before = true) {
+    TopoGraph::Vertex &v = g[ix];
+    for(unsigned i = 0; i < v.enters.size(); i++) {
+        if(v.enters[i].section == s) {
+            v.enters.insert(v.enters.begin() + (before ? i : i + 1), jx);
+            return;
+        }
+    }
+    for(unsigned i = 0; i < v.exits.size(); i++) {
+        if(v.exits[i].section == s) {
+            v.exits.insert(v.exits.begin() + (before ? i : i + 1), jx);
+            return;
+        }
+    }
+    //TODO: fix the fall through to here
+    //assert(false);
+}
+
+void double_whiskers(TopoGraph &g) {
+    for(unsigned i = 0; i < g.size(); i++) {
+        if(g[i].degree() == 1) {
+            unsigned j = i;
+            TopoGraph::Edge e = g[i][0];
+            while(true) {
+                TopoGraph::Edge next_edge = g[j][1 - g[j].find_section(e.section)];
+                boost::shared_ptr<Section> new_section = boost::shared_ptr<Section>(new Section(*e.section));
+                add_edge_at(g, j,       e.section, TopoGraph::Edge(new_section, e.other), false);
+                add_edge_at(g, e.other, e.section, TopoGraph::Edge(new_section, j), true);
+                
+                if(g[e.other].degree() == 3) {
+                    j = e.other;
+                    e = next_edge;
+                } else break;
+            }
+        }
+    }
+}
+
+/*
+void remove_degenerate(TopoGraph &g) {
+    for(unsigned i = 0; i < g.size(); i++) {
+        for(int j = g[i].degree(); j >= 0; j--) {
+            if(g[i][j].other == i) 
+        }
+    }
+}*/
+
+/*
+void remove_vestigial(TopoGraph &g) {
+    for(unsigned i = 0; i < g.size(); i++) {
+        if(g[i].enters.size() == 1 && g[i].exits.size() == 1) {
+            TopoGraph::Edge &e1 = g[i][0], &e2 = g[i][1];
+            if(e1.section == e2.section) {
+                //vestigial vert
+                Section *new_section = new Section(e1.section->curve,
+                                         e1.section->f,  e2.section->t,
+                                         e1.section->fp, e2.section->tp);
+                
+                e1.other
+                
+                Vertex *v1 = e1.other, *v2 = e2.other;
+                v1->lookup_section(e1.section) = Edge(new_section, v2);
+                v2->lookup_section(e2.section) = Edge(new_section, v1);
+                g.erase(g.begin() + i);
+            }    
+        }
+    }
+}*/
+
+//planar area finding
+//linear on number of edges
+Areas traverse_areas(TopoGraph const &g) {
+    Areas ret;
+    
+    //stores which edges we've visited
+    std::vector<std::vector<bool> > visited;
+    for(unsigned i = 0; i < g.size(); i++) visited.push_back(std::vector<bool>(g[i].degree(), false));
+    
+    for(unsigned vix = 0; vix < g.size(); vix++) {
+        while(true) {
+            //find an unvisited edge to start on
+            
+            unsigned e_ix = std::find(visited[vix].begin(), visited[vix].end(), false) - visited[vix].begin();
+            if(e_ix == g[vix].degree()) break;
+            
+            unsigned start = e_ix;
+            unsigned cur = vix;
+            
+            Area area;
+            //std::vector<std::vector<bool> > before(visited);
+            while(cur < g.size() && !visited[cur][e_ix]) {
+                visited[cur][e_ix] = true;
+                
+                TopoGraph::Edge e = g[cur][e_ix];
+                
+                area.push_back(e.section);
+                
+                //go to clockwise edge
+                cur = e.other;
+                unsigned deg = g[cur].degree();
+                e_ix = g[cur].find_section(e.section);
+                
+                if(deg == 1 || e_ix == deg) {
+                   visited[cur][e_ix] = true;
+                   break;
+                }
+                
+                e_ix = (e_ix + 1) % deg;
+                
+                if(cur == vix && start == e_ix) break;
+            }
+            //if(vix == cur && start == e_ix) {
+                ret.push_back(area);
+            //} else visited = before;
+        }
+    }
+    return ret;
+}
+
+void remove_area_whiskers(Areas &areas) {
+    for(int i = areas.size() - 1; i >= 0; i--)
+        if(areas[i].size() == 2 && *areas[i][0] == *areas[i][1]) 
+            areas.erase(areas.begin() + i);
+}
+
+Path area_to_path(PathVector const &ps, Area const &area) {
+    Path ret;
+    if(area.size() == 0) return ret;
+    Point prev = area[0]->fp;
+    for(unsigned i = 0; i < area.size(); i++) {
+        bool forward = are_near(area[i]->fp, prev, 0.01);
+        Curve *curv = area[i]->curve.get(ps).portion(
+                          forward ? area[i]->f : area[i]->t,
+                          forward ? area[i]->t : area[i]->f);
+        ret.append(*curv, Path::STITCH_DISCONTINUOUS);
+        delete curv;
+        prev = forward ? area[i]->tp : area[i]->fp;
+    }
+    return ret;
+}
+
+PathVector areas_to_paths(PathVector const &ps, Areas const &areas) {
+    std::vector<Path> ret;
+    ret.reserve(areas.size());
+    for(unsigned i = 0; i < areas.size(); i++)
+        ret.push_back(area_to_path(ps, areas[i]));
+    return ret;
+}
+
+} // end namespace Geom
diff --git a/src/2geom/toposweep.h b/src/2geom/toposweep.h
new file mode 100644
index 000000000..428115dd3
--- /dev/null
+++ b/src/2geom/toposweep.h
@@ -0,0 +1,222 @@
+
+/**
+ * \file
+ * \brief  TopoSweep - topology / graph representation of a PathVector, for boolean operations and related tasks
+ *
+ * Authors:
+ * 		Michael Sloan <mgsloan at gmail.com>
+ * 		Nathan Hurst <njhurst at njhurst.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 SEEN_GEOM_TOPOSWEEP_H
+#define SEEN_GEOM_TOPOSWEEP_H
+
+#include <2geom/coord.h>
+#include <2geom/point.h>
+#include <2geom/pathvector.h>
+#include <2geom/rect.h>
+#include <2geom/path.h>
+#include <2geom/curve.h>
+
+#include <boost/shared_ptr.hpp>
+
+namespace Geom {
+
+// indicates a particular curve in a pathvector
+struct CurveIx {
+    unsigned path, ix;
+    CurveIx(unsigned p, unsigned i) : path(p), ix(i) {}
+    // retrieves the indicated curve from the pathvector
+    Curve const &get(PathVector const &ps) const {
+        return ps[path][ix];
+    }
+    bool operator==(CurveIx const &other) const {
+        return other.path == path && other.ix == ix;
+    }
+};
+
+// represents a monotonic section of a path
+struct Section {
+    CurveIx curve;
+    double f, t;
+    Point fp, tp;
+    std::vector<int> windings;
+    Section(CurveIx cix, double fd, double td, Point fdp, Point tdp) : curve(cix), f(fd), t(td), fp(fdp), tp(tdp) { }
+    Section(CurveIx cix, double fd, double td, PathVector ps, Dim2 d) : curve(cix), f(fd), t(td) {
+        fp = curve.get(ps).pointAt(f), tp = curve.get(ps).pointAt(t);
+        if (Point::LexOrderRt(d)(tp, fp)) {
+            //swap from and to, since tp is left or above fp
+            std::swap(f, t);
+            std::swap(fp, tp);
+        }
+    }
+    Rect bbox() const { return Rect(fp, tp); }
+    bool operator==(Section const &other) const {
+        return (curve == other.curve) && (f == other.f) && (t == other.t);
+    }
+};
+
+class TopoGraph {
+  public:
+
+    // Represents an e    double tol;dge on a vertex
+    class Edge {
+      public:
+        boost::shared_ptr<Section> section;    // section associated with this edge
+        unsigned other; // index of the vertex this edge points to
+        Edge(boost::shared_ptr<Section> s, unsigned o) : section(s), other(o) {}
+    };
+    
+    // Represents a vertex in the graph, in terms of a point and edges which enter and exit.
+    // A vertex has an "avg" point, which is a representative point for the vertex.  All
+    // edges have an endpoint tol away.
+    class Vertex {
+      public:
+        std::vector<Edge> enters, exits; // indexes of the enter / exit edges
+        Point avg;
+        Vertex(Point p) : avg(p) {}
+        inline unsigned degree() const { return enters.size() + exits.size(); }
+        Edge operator[](unsigned ix) const;
+        Edge &operator[](unsigned ix);
+        void erase(unsigned ix);
+        void insert(unsigned ix, Edge e);
+        unsigned find_section(boost::shared_ptr<Section> section) const;
+    };
+    
+    TopoGraph(PathVector const &ps, Dim2 d, double t);
+    
+    unsigned size() const { return vertices.size(); }
+    
+    Vertex       &operator[](unsigned ix)       { return vertices[ix]; }
+    Vertex const &operator[](unsigned ix) const { return vertices[ix]; }
+    
+    //removes both edges, and returns the vertices[ix][jx] one
+    Edge remove_edge(unsigned ix, unsigned jx);
+    
+    //returns a graph with all zero degree vertices and unused edges removed
+    void cannonize();
+    
+    //checks invariants
+    void assert_invariants() const;
+    
+    std::vector<Vertex> vertices;
+    Dim2 dim;
+    double tol;
+};
+
+//TODO: convert to classes
+typedef std::vector<boost::shared_ptr<Section> > Area;
+typedef std::vector<Area> Areas;
+
+//TopoGraph sweep_graph(PathVector const &ps, Dim2 d = X, double tol = 0.00001);
+
+void trim_whiskers(TopoGraph &g);
+void double_whiskers(TopoGraph &g);
+//void remove_degenerate(TopoGraph &g);
+//void remove_vestigial(TopoGraph &g);
+//Areas traverse_areas(TopoGraph const &g);
+
+
+void remove_area_whiskers(Areas &areas);
+PathVector areas_to_paths(PathVector const &ps, Areas const &areas);
+
+class SectionSorter {
+    const PathVector &ps;
+    Dim2 dim;
+    double tol;
+    bool section_order(Section const &a, double at, Section const &b, double bt) const;
+  public:
+    typedef Section first_argument_type;
+    typedef Section second_argument_type;
+    typedef bool result_type;
+  
+    SectionSorter(const PathVector &rs, Dim2 d, double t = 0.00001) : ps(rs), dim(d), tol(t) {}
+    bool operator()(Section const &a, Section const &b) const;
+};
+
+//sorter used to create the initial sweep of sections, such that they are dealt with in order
+struct SweepSorter {
+    typedef Section first_argument_type;
+    typedef Section second_argument_type;
+    typedef bool result_type;
+    Dim2 dim;
+    SweepSorter(Dim2 d) : dim(d) {}
+    bool operator()(const Section &a, const Section &b) const {
+        return Point::LexOrderRt(dim)(a.fp, b.fp);
+    }
+};
+
+struct UnionOp {
+    unsigned ix;
+    bool nz1, nz2;
+    UnionOp(unsigned i, bool a, bool b) : ix(i), nz1(a), nz2(b) {}
+    bool operator()(std::vector<int> const &windings) const {
+        int w1 = 0, w2 = 0;
+        for(unsigned j = 0; j < ix; j++) w1 += windings[j];
+        for(unsigned j = ix; j < windings.size(); j++) w2 += windings[j];
+        return (nz1 ? w1 : w1 % 2) != 0 || (nz2 ? w2 : w2 % 2) != 0;
+    }
+};
+
+//returns all areas for which the winding -> bool function yields true
+template<class Z>
+Areas filter_areas(PathVector const &ps, Areas const & areas, Z const &z) {
+    Areas ret;
+    SweepSorter sorty = SweepSorter(Y);
+    SectionSorter sortx = SectionSorter(ps, X);
+    for(unsigned i = 0; i < areas.size(); i++) {
+        if(areas[i].size() < 2) continue;
+        //find a representative section
+        unsigned rj = 0;
+        bool rev = are_near(areas[i][0]->fp, areas[i][1]->tp);
+        for(unsigned j = 1; j < areas[i].size(); j++)
+            if(sorty(*areas[i][rj], *areas[i][j])) rj = j;
+        if(sortx(*areas[i][rj], *areas[i][(rj+areas[i].size() - 1) % areas[i].size()])) {
+            rj = 0;
+            for(unsigned j = 1; j < areas[i].size(); j++)
+                if(sorty(*areas[i][j], *areas[i][rj])) rj = j;
+        }
+        if(z(areas[i][rj]->windings)) ret.push_back(areas[i]);
+    }
+    return ret;
+}
+
+} // end namespace Geom
+
+#endif // SEEN_GEOM_TOPOSWEEP_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 3a1866c13..2658719c4 100644
--- a/src/2geom/transforms.cpp
+++ b/src/2geom/transforms.cpp
@@ -60,12 +60,12 @@ Point &Point::operator*=(Rotate const &r)
 }
 Point &Point::operator*=(HShear const &h)
 {
-    _pt[X] += h.f * _pt[Y];
+    _pt[X] += h.f * _pt[X];
     return *this;
 }
 Point &Point::operator*=(VShear const &v)
 {
-    _pt[Y] += v.f * _pt[X];
+    _pt[Y] += v.f * _pt[Y];
     return *this;
 }
 
diff --git a/src/2geom/transforms.h b/src/2geom/transforms.h
index 48d4b1dba..9623bed26 100644
--- a/src/2geom/transforms.h
+++ b/src/2geom/transforms.h
@@ -32,12 +32,12 @@
  * the specific language governing rights and limitations.
  */
 
-#ifndef SEEN_Geom_TRANSFORMS_H
-#define SEEN_Geom_TRANSFORMS_H
+#ifndef LIB2GEOM_SEEN_TRANSFORMS_H
+#define LIB2GEOM_SEEN_TRANSFORMS_H
 
+#include <cmath>
 #include <2geom/forward.h>
 #include <2geom/affine.h>
-#include <cmath>
 
 namespace Geom {
 
@@ -66,7 +66,8 @@ struct TransformConcept {
     }
 };
 
-/** @brief Base template for transforms. */
+/** @brief Base template for transforms.
+ * This class is an implementation detail and should not be used directly. */
 template <typename T>
 class TransformOperations
     : boost::equality_comparable< T
@@ -198,6 +199,7 @@ public:
 };
 
 /** @brief Common base for shearing transforms.
+ * This class is an implementation detail and should not be used directly.
  * @ingroup Transforms */
 template <typename S>
 class ShearBase
@@ -259,10 +261,9 @@ inline Translate pow(Translate const &t, int n) {
 
 //TODO: matrix to trans/scale/rotate
 
-} /* namespace Geom */
-
+} // end namespace Geom
 
-#endif /* !SEEN_Geom_TRANSFORMS_H */
+#endif // LIB2GEOM_SEEN_TRANSFORMS_H
 
 /*
   Local Variables:
diff --git a/src/2geom/utils.h b/src/2geom/utils.h
index e90a4623b..6a72d42c4 100644
--- a/src/2geom/utils.h
+++ b/src/2geom/utils.h
@@ -1,10 +1,7 @@
-#ifndef LIB2GEOM_UTILS_HEADER
-#define LIB2GEOM_UTILS_HEADER
-
 /**
  * \file
  * \brief  Various utility functions.
- *
+ *//*
  * Copyright 2007 Johan Engelen <goejendaagh@zonnet.nl>
  * Copyright 2006 Michael G. Sloan <mgsloan@gmail.com>
  *
@@ -33,6 +30,9 @@
  *
  */
 
+#ifndef SEEN_LIB2GEOM_UTILS_H
+#define SEEN_LIB2GEOM_UTILS_H
+
 #include <cstddef>
 #include <vector>
 
@@ -59,9 +59,9 @@ struct MultipliableNoncommutative : B
     }
 };
 
-}
+} // end namespace Geom
 
-#endif
+#endif // SEEN_LIB2GEOM_UTILS_H
 
 /*
   Local Variables:
diff --git a/src/connector-context.cpp b/src/connector-context.cpp
index 251b41066..2aa9c41ee 100644
--- a/src/connector-context.cpp
+++ b/src/connector-context.cpp
@@ -1306,12 +1306,12 @@ cc_connector_rerouting_finish(SPConnectorContext *const cc, Geom::Point *const p
 
         if (found) {
             if (cc->clickedhandle == cc->endpt_handle[0]) {
-                cc->clickeditem->setAttribute("inkscape:connection-start", shape_label, false);
-                cc->clickeditem->setAttribute("inkscape:connection-start-point", cpid, false);
+                cc->clickeditem->setAttribute("inkscape:connection-start", shape_label, NULL);
+                cc->clickeditem->setAttribute("inkscape:connection-start-point", cpid, NULL);
             }
             else {
-                cc->clickeditem->setAttribute("inkscape:connection-end", shape_label, false);
-                cc->clickeditem->setAttribute("inkscape:connection-end-point", cpid, false);
+                cc->clickeditem->setAttribute("inkscape:connection-end", shape_label, NULL);
+                cc->clickeditem->setAttribute("inkscape:connection-end-point", cpid, NULL);
             }
             g_free(shape_label);
         }
@@ -1451,23 +1451,23 @@ spcc_flush_white(SPConnectorContext *cc, SPCurve *gc)
 
         bool connection = false;
         cc->newconn->setAttribute( "inkscape:connector-type",
-                                   cc->isOrthogonal ? "orthogonal" : "polyline", false );
+                                   cc->isOrthogonal ? "orthogonal" : "polyline", NULL );
         cc->newconn->setAttribute( "inkscape:connector-curvature",
-                                   Glib::Ascii::dtostr(cc->curvature).c_str(), false );
+                                   Glib::Ascii::dtostr(cc->curvature).c_str(), NULL );
         if (cc->shref)
         {
-            cc->newconn->setAttribute( "inkscape:connection-start", cc->shref, false);
+            cc->newconn->setAttribute( "inkscape:connection-start", cc->shref, NULL);
             if (cc->scpid) {
-                cc->newconn->setAttribute( "inkscape:connection-start-point", cc->scpid, false);
+                cc->newconn->setAttribute( "inkscape:connection-start-point", cc->scpid, NULL);
             }
             connection = true;
         }
 
         if (cc->ehref)
         {
-            cc->newconn->setAttribute( "inkscape:connection-end", cc->ehref, false);
+            cc->newconn->setAttribute( "inkscape:connection-end", cc->ehref, NULL);
             if (cc->ecpid) {
-                cc->newconn->setAttribute( "inkscape:connection-end-point", cc->ecpid, false);
+                cc->newconn->setAttribute( "inkscape:connection-end-point", cc->ecpid, NULL);
             }
             connection = true;
         }
@@ -1950,7 +1950,7 @@ void cc_selection_set_avoid(bool const set_avoid)
         char const *value = (set_avoid) ? "true" : NULL;
 
         if (cc_item_is_shape(item)) {
-            item->setAttribute("inkscape:connector-avoid", value, false);
+            item->setAttribute("inkscape:connector-avoid", value, NULL);
             item->avoidRef->handleSettingChange();
             changes++;
         }
diff --git a/src/display/nr-arena-image.cpp b/src/display/nr-arena-image.cpp
index 36d733eb8..a943a6214 100644
--- a/src/display/nr-arena-image.cpp
+++ b/src/display/nr-arena-image.cpp
@@ -301,7 +301,7 @@ nr_arena_image_rect (NRArenaImage *image)
         Geom::Point p(image->ox, image->oy);
         Geom::Point wh(vw, vh);
         Geom::Rect view(p, p+wh);
-        Geom::OptRect res = Geom::intersect(r, view);
+        Geom::OptRect res = r & view;
         r = res ? *res : r;
     }
 
diff --git a/src/display/nr-filter-composite.cpp b/src/display/nr-filter-composite.cpp
index d4cf47af4..694ccaec5 100644
--- a/src/display/nr-filter-composite.cpp
+++ b/src/display/nr-filter-composite.cpp
@@ -11,7 +11,6 @@
 
 #include <cmath>
 
-#include "2geom/isnan.h"
 #include "display/cairo-templates.h"
 #include "display/cairo-utils.h"
 #include "display/nr-filter-composite.h"
diff --git a/src/display/nr-filter-gaussian.cpp b/src/display/nr-filter-gaussian.cpp
index 326c37160..884e832ef 100644
--- a/src/display/nr-filter-gaussian.cpp
+++ b/src/display/nr-filter-gaussian.cpp
@@ -25,8 +25,6 @@
 #include <omp.h>
 #endif //HAVE_OPENMP
 
-#include "2geom/isnan.h"
-
 #include "display/cairo-utils.h"
 #include "display/nr-filter-primitive.h"
 #include "display/nr-filter-gaussian.h"
diff --git a/src/dyna-draw-context.cpp b/src/dyna-draw-context.cpp
index aa7d840bc..a3a665b1c 100644
--- a/src/dyna-draw-context.cpp
+++ b/src/dyna-draw-context.cpp
@@ -34,7 +34,7 @@
 #include "svg/svg.h"
 #include "display/canvas-bpath.h"
 #include "display/cairo-utils.h"
-#include <2geom/isnan.h>
+#include <2geom/math-utils.h>
 #include <2geom/pathvector.h>
 #include <2geom/bezier-utils.h>
 #include "display/curve.h"
diff --git a/src/eraser-context.cpp b/src/eraser-context.cpp
index 8ac765b9e..de6c7d86f 100644
--- a/src/eraser-context.cpp
+++ b/src/eraser-context.cpp
@@ -62,7 +62,7 @@
 #include "display/canvas-bpath.h"
 #include "display/canvas-arena.h"
 #include "livarot/Shape.h"
-#include <2geom/isnan.h>
+#include <2geom/math-utils.h>
 #include <2geom/pathvector.h>
 
 #include "eraser-context.h"
diff --git a/src/helper/recthull.h b/src/helper/recthull.h
index a9cad4466..d82450ce8 100644
--- a/src/helper/recthull.h
+++ b/src/helper/recthull.h
@@ -38,7 +38,7 @@ public:
 	void add(Rect const &r) {
 		// Note that this is a hack.  when convexhull actually works
 		// you will need to add all four points.
-        _bounds = unify(_bounds, r);
+        _bounds.unionWith(r);
 	}
 	void add(RectHull const &h) {
         if (h._bounds) {
diff --git a/src/libcola/cola.cpp b/src/libcola/cola.cpp
index 2a3b525a7..62771ece2 100644
--- a/src/libcola/cola.cpp
+++ b/src/libcola/cola.cpp
@@ -2,7 +2,7 @@
 #include "conjugate_gradient.h"
 #include "straightener.h"
 #include "shortest_paths.h"
-#include "2geom/isnan.h"
+#include <2geom/math-utils.h>
 
 namespace cola {
 
diff --git a/src/libcola/gradient_projection.cpp b/src/libcola/gradient_projection.cpp
index fb8702ec7..47109a4b0 100644
--- a/src/libcola/gradient_projection.cpp
+++ b/src/libcola/gradient_projection.cpp
@@ -17,7 +17,7 @@
 #include <libvpsc/constraint.h>
 #include "gradient_projection.h"
 #include <iostream>
-#include "2geom/isnan.h"
+#include <2geom/math-utils.h>
 #include "isinf.h"
 #include <math.h>
 
diff --git a/src/libnr/nr-point-fns.cpp b/src/libnr/nr-point-fns.cpp
index cd6d6927b..ac58eddb7 100644
--- a/src/libnr/nr-point-fns.cpp
+++ b/src/libnr/nr-point-fns.cpp
@@ -1,5 +1,5 @@
 #include <libnr/nr-point-fns.h>
-#include <2geom/isnan.h>
+#include <2geom/math-utils.h>
 
 using NR::Point;
 
diff --git a/src/libnr/nr-types.cpp b/src/libnr/nr-types.cpp
index 0231c91d5..5da5d5cf6 100644
--- a/src/libnr/nr-types.cpp
+++ b/src/libnr/nr-types.cpp
@@ -3,8 +3,7 @@
  */
 
 #include <libnr/nr-types.h>
-
-#include "2geom/isnan.h"
+#include <2geom/math-utils.h>
 
 /** Scales this vector to make it a unit vector (within rounding error).
  *
diff --git a/src/libvpsc/generate-constraints.cpp b/src/libvpsc/generate-constraints.cpp
index c57966e26..0c35ab51c 100644
--- a/src/libvpsc/generate-constraints.cpp
+++ b/src/libvpsc/generate-constraints.cpp
@@ -16,7 +16,7 @@
 #include "generate-constraints.h"
 #include "constraint.h"
 
-#include "2geom/isnan.h" /* Include last */
+#include <2geom/math-utils.h>
 
 using std::set;
 using std::vector;
diff --git a/src/live_effects/lpe-spiro.cpp b/src/live_effects/lpe-spiro.cpp
index 54554ebb2..22974fe13 100644
--- a/src/live_effects/lpe-spiro.cpp
+++ b/src/live_effects/lpe-spiro.cpp
@@ -12,7 +12,6 @@
 #include <2geom/affine.h>
 #include <2geom/bezier-curve.h>
 #include <2geom/hvlinesegment.h>
-#include <2geom/isnan.h>
 #include "helper/geom-nodetype.h"
 #include "helper/geom-curves.h"
 
diff --git a/src/object-edit.cpp b/src/object-edit.cpp
index 553c125a3..743ef573a 100644
--- a/src/object-edit.cpp
+++ b/src/object-edit.cpp
@@ -35,7 +35,7 @@
 #include <glibmm/i18n.h>
 #include "object-edit.h"
 #include "xml/repr.h"
-#include "2geom/isnan.h"
+#include <2geom/math-utils.h>
 
 #define sp_round(v,m) (((v) < 0.0) ? ((ceil((v) / (m) - 0.5)) * (m)) : ((floor((v) / (m) + 0.5)) * (m)))
 
diff --git a/src/selection-chemistry.cpp b/src/selection-chemistry.cpp
index 6f385b8f5..9b88077e7 100644
--- a/src/selection-chemistry.cpp
+++ b/src/selection-chemistry.cpp
@@ -777,7 +777,7 @@ enclose_items(GSList const *items)
 
     Geom::OptRect r;
     for (GSList const *i = items; i; i = i->next) {
-        r = Geom::unify(r, ((SPItem *) i->data)->getBboxDesktop());
+        r.unionWith(((SPItem *) i->data)->getBboxDesktop());
     }
     return r;
 }
diff --git a/src/selection.cpp b/src/selection.cpp
index 7564fad3a..3c4ccccf2 100644
--- a/src/selection.cpp
+++ b/src/selection.cpp
@@ -375,7 +375,7 @@ Geom::OptRect Selection::bounds(SPItem::BBoxType type) const
 
     Geom::OptRect bbox;
     for ( GSList const *i = items ; i != NULL ; i = i->next ) {
-        bbox = unify(bbox, SP_ITEM(i->data)->getBboxDesktop(type));
+        bbox.unionWith(SP_ITEM(i->data)->getBboxDesktop(type));
     }
     return bbox;
 }
diff --git a/src/sp-item.cpp b/src/sp-item.cpp
index 7e5f5f96a..424107426 100644
--- a/src/sp-item.cpp
+++ b/src/sp-item.cpp
@@ -825,7 +825,7 @@ void SPItem::invoke_bbox_full( Geom::OptRect &bbox, Geom::Affine const &transfor
     // would therefore be translated into empty Geom::OptRect() (see bug https://bugs.launchpad.net/inkscape/+bug/168684)
     Geom::OptRect temp_bbox_new = Geom::Rect(Geom::Point(temp_bbox.x0, temp_bbox.y0), Geom::Point(temp_bbox.x1, temp_bbox.y1));
 
-    bbox = Geom::unify(bbox, temp_bbox_new);
+    bbox.unionWith(temp_bbox_new);
 }
 
 // DEPRECATED to phase out the use of NRRect in favor of Geom::OptRect
diff --git a/src/spray-context.cpp b/src/spray-context.cpp
index 36c135924..aa14e6ee5 100644
--- a/src/spray-context.cpp
+++ b/src/spray-context.cpp
@@ -56,7 +56,6 @@
 #include "display/canvas-arena.h"
 #include "display/curve.h"
 #include "livarot/Shape.h"
-#include <2geom/isnan.h>
 #include <2geom/transforms.h>
 #include "preferences.h"
 #include "style.h"
diff --git a/src/style.cpp b/src/style.cpp
index bb25a5f46..37a784e2a 100644
--- a/src/style.cpp
+++ b/src/style.cpp
@@ -42,7 +42,6 @@
 #include "xml/repr.h"
 #include "xml/simple-document.h"
 #include "unit-constants.h"
-#include "2geom/isnan.h"
 #include "macros.h"
 #include "preferences.h"
 
diff --git a/src/tweak-context.cpp b/src/tweak-context.cpp
index faa08ee91..022869c69 100644
--- a/src/tweak-context.cpp
+++ b/src/tweak-context.cpp
@@ -63,7 +63,6 @@
 #include "display/canvas-arena.h"
 #include "display/curve.h"
 #include "livarot/Shape.h"
-#include <2geom/isnan.h>
 #include <2geom/transforms.h>
 #include "preferences.h"
 #include "style.h"
diff --git a/src/widgets/desktop-widget.cpp b/src/widgets/desktop-widget.cpp
index 0d890fa86..6f3b4dcb9 100644
--- a/src/widgets/desktop-widget.cpp
+++ b/src/widgets/desktop-widget.cpp
@@ -1161,7 +1161,7 @@ bool SPDesktopWidget::showInfoDialog( Glib::ustring const &message )
                 GTK_BUTTONS_OK,
                 "%s", message.c_str());
         gtk_window_set_title( GTK_WINDOW(dialog), _("Note:")); // probably want to take this as a parameter.
-        gint response = gtk_dialog_run(GTK_DIALOG(dialog));
+        gtk_dialog_run(GTK_DIALOG(dialog));
         gtk_widget_destroy(dialog);
     }
     return result;
@@ -1887,7 +1887,7 @@ sp_desktop_widget_update_scrollbars (SPDesktopWidget *dtw, double scale)
     Geom::Rect darea ( Geom::Point(-doc->getWidth(), -doc->getHeight()),
                      Geom::Point(2 * doc->getWidth(), 2 * doc->getHeight())  );
 
-    Geom::OptRect deskarea = Geom::unify(darea, doc->getRoot()->getBboxDesktop());
+    Geom::OptRect deskarea = darea | doc->getRoot()->getBboxDesktop();
 
     /* Canvas region we always show unconditionally */
     Geom::Rect carea( Geom::Point(deskarea->min()[Geom::X] * scale - 64, deskarea->max()[Geom::Y] * -scale - 64),