Merge from fe-moved

(bzr r6891)

39 files changed, 696 insertions, 406 deletions

diff --git a/src/2geom/bezier-curve.h b/src/2geom/bezier-curve.h
index 135aa4f71..589335d22 100644
--- a/src/2geom/bezier-curve.h
+++ b/src/2geom/bezier-curve.h
@@ -104,15 +104,16 @@ public:
   void setPoint(unsigned ix, Point v) { inner[X].setPoint(ix, v[X]); inner[Y].setPoint(ix, v[Y]); }
   Point const operator[](unsigned ix) const { return Point(inner[X][ix], inner[Y][ix]); }
 
-  Rect boundsFast() const { return bounds_fast(inner); }
-  Rect boundsExact() const { return bounds_exact(inner); }
-  Rect boundsLocal(Interval i, unsigned deg) const {
-      if(i.min() == 0 && i.max() == 1) return boundsFast();
+  virtual OptRect boundsFast() const { return bounds_fast(inner); }
+  virtual OptRect boundsExact() const { return bounds_exact(inner); }
+  virtual OptRect boundsLocal(OptInterval i, unsigned deg) const {
+      if (!i) return OptRect();
+      if(i->min() == 0 && i->max() == 1) return boundsFast();
       if(deg == 0) return bounds_local(inner, i);
       // TODO: UUUUUUGGGLLY
-      if(deg == 1 && order > 1) return Rect(bounds_local(Geom::derivative(inner[X]), i),
-                                            bounds_local(Geom::derivative(inner[Y]), i));
-      return Rect(Interval(0,0), Interval(0,0));
+      if(deg == 1 && order > 1) return OptRect(bounds_local(Geom::derivative(inner[X]), i),
+                                               bounds_local(Geom::derivative(inner[Y]), i));
+      return OptRect();
   }
 //TODO: local
 
diff --git a/src/2geom/bezier.h b/src/2geom/bezier.h
index 94dd909ca..889dde9ed 100644
--- a/src/2geom/bezier.h
+++ b/src/2geom/bezier.h
@@ -72,7 +72,7 @@ private:
 
     friend Bezier portion(const Bezier & a, Coord from, Coord to);
 
-    friend Interval bounds_fast(Bezier const & b);
+    friend OptInterval bounds_fast(Bezier const & b);
 
     friend Bezier derivative(const Bezier & a);
 
@@ -185,10 +185,9 @@ public:
     std::vector<Coord> valueAndDerivatives(Coord t, unsigned n_derivs) const {
         std::vector<Coord> val_n_der;
         Coord d_[order()+1];
-        unsigned nn = n_derivs + 1; // the size of the result vector equals n_derivs+1
+        unsigned nn = n_derivs + 1; 	// the size of the result vector equals n_derivs+1 ...
         if(nn > order())
-            //nn = order();
-            nn = order()+1;
+            nn = order()+1; 			// .. but with a maximum of order() + 1!
         for(unsigned i = 0; i < size(); i++)
             d_[i] = c_[i];
         for(unsigned di = 0; di < nn; di++) {
@@ -309,18 +308,22 @@ inline Bezier integral(const Bezier & a) {
     return inte;
 }
 
-inline Interval bounds_fast(Bezier const & b) {
+inline OptInterval bounds_fast(Bezier const & b) {
     return Interval::fromArray(&b.c_[0], b.size());
 }
 
 //TODO: better bounds exact
-inline Interval bounds_exact(Bezier const & b) {
+inline OptInterval bounds_exact(Bezier const & b) {
     return bounds_exact(b.toSBasis());
 }
 
-inline Interval bounds_local(Bezier const & b, Interval i) {
-    return bounds_fast(portion(b, i.min(), i.max()));
+inline OptInterval bounds_local(Bezier const & b, OptInterval i) {
     //return bounds_local(b.toSBasis(), i);
+    if (i) {
+        return bounds_fast(portion(b, i->min(), i->max()));
+    } else {
+        return OptInterval();
+    }
 }
 
 inline std::ostream &operator<< (std::ostream &out_file, const Bezier & b) {
diff --git a/src/2geom/concepts.h b/src/2geom/concepts.h
index 8f4d98ef2..9c57db44d 100644
--- a/src/2geom/concepts.h
+++ b/src/2geom/concepts.h
@@ -37,21 +37,20 @@
 #include <2geom/point.h>
 #include <vector>
 #include <boost/concept_check.hpp>
+#include <2geom/forward.h>
 
 namespace Geom {
 
 //forward decls
-template <typename T> class D2;
-
 template <typename T> struct ResultTraits;
 
 template <> struct ResultTraits<double> {
-  typedef Interval bounds_type;
+  typedef OptInterval bounds_type;
   typedef SBasis sb_type;
 };
 
 template <> struct ResultTraits<Point > {
-  typedef D2<Interval> bounds_type;
+  typedef OptRect bounds_type;
   typedef D2<SBasis> sb_type;
 };
 
diff --git a/src/2geom/crossing.h b/src/2geom/crossing.h
index 291f69382..3ad034ac9 100644
--- a/src/2geom/crossing.h
+++ b/src/2geom/crossing.h
@@ -120,7 +120,12 @@ typedef std::vector<Crossings> CrossingSet;
 template<typename C>
 std::vector<Rect> bounds(C const &a) {
     std::vector<Rect> rs;
-    for(unsigned i = 0; i < a.size(); i++) rs.push_back(a[i].boundsFast());
+    for (unsigned i = 0; i < a.size(); i++) {
+        OptRect bb = a[i].boundsFast();
+        if (bb) {
+            rs.push_back(*bb);
+        }
+    }
     return rs;
 }
 
diff --git a/src/2geom/curve.h b/src/2geom/curve.h
index b81548ba8..af02d6edb 100644
--- a/src/2geom/curve.h
+++ b/src/2geom/curve.h
@@ -74,10 +74,10 @@ public:
 
   virtual Curve *duplicate() const = 0;
 
-  virtual Rect boundsFast() const = 0;
-  virtual Rect boundsExact() const = 0;
-  virtual Rect boundsLocal(Interval i, unsigned deg) const = 0;
-  Rect boundsLocal(Interval i) const { return boundsLocal(i, 0); }
+  virtual OptRect boundsFast() const = 0;
+  virtual OptRect boundsExact() const = 0;
+  virtual OptRect boundsLocal(OptInterval i, unsigned deg) const = 0;
+  OptRect boundsLocal(OptInterval i) const { return boundsLocal(i, 0); }
 
   virtual std::vector<double> roots(double v, Dim2 d) const = 0;
 
@@ -133,12 +133,12 @@ public:
    */
   virtual Point unitTangentAt(Coord t, unsigned n = 3) const
   {
-    for (unsigned deriv_n = 1; deriv_n <= n; deriv_n++) {
-      Point deriv = pointAndDerivatives(t, deriv_n)[deriv_n];
-      Coord length = deriv.length();
+	std::vector<Point> derivs = pointAndDerivatives(t, n);
+	for (unsigned deriv_n = 1; deriv_n < derivs.size(); deriv_n++) {
+      Coord length = derivs[deriv_n].length();
       if ( ! are_near(length, 0) ) {
          // length of derivative is non-zero, so return unit vector
-         return deriv / length;
+    	  derivs[deriv_n] / length;
       }
     }
     return Point (0,0);
diff --git a/src/2geom/d2-sbasis.cpp b/src/2geom/d2-sbasis.cpp
index 01f83bf5e..55c7ef55e 100644
--- a/src/2geom/d2-sbasis.cpp
+++ b/src/2geom/d2-sbasis.cpp
@@ -150,6 +150,81 @@ split_at_discontinuities (Geom::Piecewise<Geom::D2<Geom::SBasis> > const & pwsbi
     return ret;
 }
 
+static void set_first_point(Piecewise<D2<SBasis> > &f, Point a){
+    if ( f.empty() ){
+        f.concat(Piecewise<D2<SBasis> >(D2<SBasis>(Linear(a[X]),Linear(a[Y]))));
+        return;
+    }
+    for (unsigned dim=0; dim<2; dim++){
+        if (f.segs.front()[dim].size() == 0){
+            f.segs.front()[dim].push_back(Linear(a[dim],0));
+        }else{
+            f.segs.front()[dim][0][0] = a[dim];
+        }
+    }
+}
+static void set_last_point(Piecewise<D2<SBasis> > &f, Point a){
+    if ( f.empty() ){
+        f.concat(Piecewise<D2<SBasis> >(D2<SBasis>(Linear(a[X]),Linear(a[Y]))));
+        return;
+    }
+    for (unsigned dim=0; dim<2; dim++){
+        if (f.segs.back()[dim].size() == 0){
+            f.segs.back()[dim].push_back(Linear(0,a[dim]));
+        }else{
+            f.segs.back()[dim][0][1] = a[dim];
+        }
+    }
+}
+
+std::vector<Piecewise<D2<SBasis> > > fuse_nearby_ends(std::vector<Piecewise<D2<SBasis> > > const &f, double tol){
+
+    if ( f.size()==0 ) return f;
+    std::vector<Piecewise<D2<SBasis> > > result;
+    std::vector<std::vector<unsigned> > pre_result;
+    for (unsigned i=0; i<f.size(); i++){
+        bool inserted = false;
+        Point a = f[i].firstValue();
+        Point b = f[i].lastValue();
+        for (unsigned j=0; j<pre_result.size(); j++){
+            Point aj = f.at(pre_result[j].back()).lastValue();
+            Point bj = f.at(pre_result[j].front()).firstValue();
+            if ( L2(a-aj) < tol ) {
+                pre_result[j].push_back(i);
+                inserted = true;
+                break;
+            }
+            if ( L2(b-bj) < tol ) {
+                pre_result[j].insert(pre_result[j].begin(),i);
+                inserted = true;
+                break;
+            }
+        }
+        if (!inserted) {
+            pre_result.push_back(std::vector<unsigned>());
+            pre_result.back().push_back(i);
+        }
+    }
+    for (unsigned i=0; i<pre_result.size(); i++){
+        Piecewise<D2<SBasis> > comp;
+        for (unsigned j=0; j<pre_result[i].size(); j++){
+            Piecewise<D2<SBasis> > new_comp = f.at(pre_result[i][j]);
+            if ( j>0 ){
+                set_first_point( new_comp, comp.segs.back().at1() );
+            }
+            comp.concat(new_comp);
+        }
+        if ( L2(comp.firstValue()-comp.lastValue()) < tol ){
+            //TODO: check sizes!!!
+            set_last_point( comp, comp.segs.front().at0() ); 
+        }
+        result.push_back(comp);
+    }
+    return result;
+    return f;
+}
+
+
 }  // namespace Geom
 
 
diff --git a/src/2geom/d2-sbasis.h b/src/2geom/d2-sbasis.h
index dd1a8e11c..c61052da7 100644
--- a/src/2geom/d2-sbasis.h
+++ b/src/2geom/d2-sbasis.h
@@ -79,6 +79,8 @@ Piecewise<D2<SBasis> > operator*(Piecewise<D2<SBasis> > const &a, Matrix const &
 
 Piecewise<D2<SBasis> > force_continuity(Piecewise<D2<SBasis> > const &f, double tol=0, bool closed=false);
 
+std::vector<Piecewise<D2<SBasis> > > fuse_nearby_ends(std::vector<Piecewise<D2<SBasis> > > const &f, double tol=0);
+
 std::vector<Geom::Piecewise<Geom::D2<Geom::SBasis> > > split_at_discontinuities (Geom::Piecewise<Geom::D2<Geom::SBasis> > const & pwsbin, double tol = .0001);
 
 class CoordIterator
@@ -114,13 +116,25 @@ inline CoordIterator iterateCoord(Piecewise<D2<SBasis> > const &a, unsigned d) {
 }
 
 //bounds specializations with order
-inline Rect bounds_fast(D2<SBasis> const & s, unsigned order=0) {
-    return Rect(bounds_fast(s[X], order),
-                bounds_fast(s[Y], order));
+inline OptRect bounds_fast(D2<SBasis> const & s, unsigned order=0) {
+    OptRect retval;
+    OptInterval xint = bounds_fast(s[X], order);
+    if (xint) {
+        OptInterval yint = bounds_fast(s[Y], order);
+        if (yint) {
+            retval = Rect(*xint, *yint);
+        }
+    }
+    return retval;
 }
-inline Rect bounds_local(D2<SBasis> const & s, Interval i, unsigned order=0) {
-    return Rect(bounds_local(s[X], i, order),
-                bounds_local(s[Y], i, order));
+inline OptRect bounds_local(D2<SBasis> const & s, OptInterval i, unsigned order=0) {
+    OptRect retval;
+    OptInterval xint = bounds_local(s[X], i, order);
+    OptInterval yint = bounds_local(s[Y], i, order);
+    if (xint && yint) {
+        retval = Rect(*xint, *yint);
+    }
+    return retval;
 }
 
 }
diff --git a/src/2geom/d2.h b/src/2geom/d2.h
index 7e2bbae53..a14e3b0eb 100644
--- a/src/2geom/d2.h
+++ b/src/2geom/d2.h
@@ -400,19 +400,19 @@ namespace Geom{
 
 //Some D2 Fragment implementation which requires rect:
 template <typename T>
-Rect bounds_fast(const D2<T> &a) {
-    boost::function_requires<FragmentConcept<T> >();        
-    return Rect(bounds_fast(a[X]), bounds_fast(a[Y]));
+OptRect bounds_fast(const D2<T> &a) {
+    boost::function_requires<FragmentConcept<T> >();
+    return OptRect(bounds_fast(a[X]), bounds_fast(a[Y]));
 }
 template <typename T>
-Rect bounds_exact(const D2<T> &a) {
-    boost::function_requires<FragmentConcept<T> >();        
-    return Rect(bounds_exact(a[X]), bounds_exact(a[Y]));
+OptRect bounds_exact(const D2<T> &a) {
+    boost::function_requires<FragmentConcept<T> >();
+    return OptRect(bounds_exact(a[X]), bounds_exact(a[Y]));
 }
 template <typename T>
-Rect bounds_local(const D2<T> &a, const Interval &t) {
-    boost::function_requires<FragmentConcept<T> >();        
-    return Rect(bounds_local(a[X], t), bounds_local(a[Y], t));
+OptRect bounds_local(const D2<T> &a, const OptInterval &t) {
+    boost::function_requires<FragmentConcept<T> >();
+    return OptRect(bounds_local(a[X], t), bounds_local(a[Y], t));
 }
 };
 
diff --git a/src/2geom/elliptical-arc.cpp b/src/2geom/elliptical-arc.cpp
index bf2b72b63..f2b6b6be2 100644
--- a/src/2geom/elliptical-arc.cpp
+++ b/src/2geom/elliptical-arc.cpp
@@ -31,6 +31,7 @@
 #include <2geom/elliptical-arc.h>
 #include <2geom/bezier-curve.h>
 #include <2geom/poly.h>
+#include <2geom/sbasis-math.h>
 
 #include <cfloat>
 #include <limits>
@@ -42,7 +43,7 @@ namespace Geom
 {
 
 
-Rect EllipticalArc::boundsExact() const
+OptRect EllipticalArc::boundsExact() const
 {
 	std::vector<double> extremes(4);
 	double cosrot = std::cos(rotation_angle());
diff --git a/src/2geom/elliptical-arc.h b/src/2geom/elliptical-arc.h
index 24b4fcf46..25c79a363 100644
--- a/src/2geom/elliptical-arc.h
+++ b/src/2geom/elliptical-arc.h
@@ -168,17 +168,17 @@ class EllipticalArc : public Curve
     }
     
     
-    Rect boundsFast() const
+    virtual OptRect boundsFast() const
     {
-    	return boundsExact();
+        return boundsExact();
     }
   
-    Rect boundsExact() const;
+    virtual OptRect boundsExact() const;
     
     // TODO: native implementation of the following methods
-    Rect boundsLocal(Interval i, unsigned int deg) const
+    virtual OptRect boundsLocal(OptInterval i, unsigned int deg) const
     {
-    	return SBasisCurve(toSBasis()).boundsLocal(i, deg);
+        return SBasisCurve(toSBasis()).boundsLocal(i, deg);
     }
     
     std::vector<double> roots(double v, Dim2 d) const;
diff --git a/src/2geom/forward.h b/src/2geom/forward.h
index 184801ef4..15740faf0 100644
--- a/src/2geom/forward.h
+++ b/src/2geom/forward.h
@@ -62,6 +62,7 @@ class NotInvertible;
 class ContinuityError;
 
 class Interval;
+class OptInterval;
 class Linear;
 class Hat;
 class Tri;
@@ -81,6 +82,7 @@ class SBasis;
 class SBasisCurve;
 
 typedef D2<Interval> Rect;
+class OptRect;
 
 class Shape;
 class Region;
diff --git a/src/2geom/geom.cpp b/src/2geom/geom.cpp
index d0689981a..f9b1a664b 100644
--- a/src/2geom/geom.cpp
+++ b/src/2geom/geom.cpp
@@ -9,6 +9,7 @@
 #include <2geom/geom.h>
 #include <2geom/point.h>
 #include <algorithm>
+#include <2geom/rect.h>
 
 namespace Geom {
 
@@ -313,6 +314,14 @@ rect_line_intersect(Geom::Point const &c0, Geom::Point const &c1,
     return results;
 }
 
+std::vector<Geom::Point>
+rect_line_intersect(Geom::Rect &r,
+                    Geom::Point const &p0, Geom::Point const &p1)
+{
+    return rect_line_intersect(r.min(), r.max(), p0, p1);
+}
+
+
 /**
  * polyCentroid: Calculates the centroid (xCentroid, yCentroid) and area of a polygon, given its
  * vertices (x[0], y[0]) ... (x[n-1], y[n-1]). It is assumed that the contour is closed, i.e., that
diff --git a/src/2geom/geom.h b/src/2geom/geom.h
index aeb40f7e1..d0af7d7d2 100644
--- a/src/2geom/geom.h
+++ b/src/2geom/geom.h
@@ -38,6 +38,7 @@
 
 #include <vector>
 #include <2geom/point.h>
+#include <2geom/rect.h>
 
 namespace Geom {
 
@@ -73,6 +74,11 @@ std::vector<Geom::Point>
 rect_line_intersect(Geom::Point const &E, Geom::Point const &F,
                     Geom::Point const &p0, Geom::Point const &p1);
 
+
+std::vector<Geom::Point>
+rect_line_intersect(Geom::Rect &r,
+                    Geom::Point const &p0, Geom::Point const &p1);
+
 int centroid(std::vector<Geom::Point> p, Geom::Point& centroid, double &area);
 
 }
diff --git a/src/2geom/hvlinesegment.h b/src/2geom/hvlinesegment.h
index 3e1287682..ccb95afdb 100644
--- a/src/2geom/hvlinesegment.h
+++ b/src/2geom/hvlinesegment.h
@@ -116,17 +116,17 @@ class HLineSegment : public Curve
         m_line_seg.setFinal( Point(finalPoint()[X], _y) );
     }
 
-    Rect boundsFast() const
+    virtual OptRect boundsFast() const
     {
         return boundsExact();
     }
 
-    Rect boundsExact() const
+    virtual OptRect boundsExact() const
     {
         return Rect( initialPoint(), finalPoint() );
     }
 
-    Rect boundsLocal(Interval i, unsigned deg) const
+    virtual OptRect boundsLocal(OptInterval i, unsigned deg) const
     {
         return m_line_seg.boundsLocal(i, deg);
     }
@@ -355,17 +355,17 @@ class VLineSegment : public Curve
         m_line_seg.setFinal( Point(_x, finalPoint()[Y]) );
     }
 
-    Rect boundsFast() const
+    virtual OptRect boundsFast() const
     {
         return boundsExact();
     }
 
-    Rect boundsExact() const
+    virtual OptRect boundsExact() const
     {
         return Rect( initialPoint(), finalPoint() );
     }
 
-    Rect boundsLocal(Interval i, unsigned deg) const
+    virtual OptRect boundsLocal(OptInterval i, unsigned deg) const
     {
         return m_line_seg.boundsLocal(i, deg);
     }
diff --git a/src/2geom/interval.h b/src/2geom/interval.h
index 8f7a0b2a1..009528586 100644
--- a/src/2geom/interval.h
+++ b/src/2geom/interval.h
@@ -44,17 +44,20 @@
 
 namespace Geom {
 
-/* Although an Interval where _b[0] > _b[1] is considered empty, for proper functioning of other methods,
- * a proper empty Interval is [+infinity, -infinity]. Then, expandTo(p) will set the interval to [p,p].
+class Interval;
+
+/** 
+ * \brief This class represents a range of numbers that is never empty.
+ *
+ * The endpoints are included in the range.
  */
 class Interval {
 private:
     Coord _b[2];
 
 public:
-    // The default constructor creates an empty interval, that ranges from +infinity to -infinity.
-    // Doing an expandTo(p) on this empty interval will correctly set the whole interval to [p,p].
-    explicit Interval() { _b[0] = +infinity();  _b[1] = -infinity(); }
+    /// The default constructor creates an interval [0,0]  DO NOT RELY ON THIS, BEST NOT TO USE THIS CONSTRUCTOR
+    explicit Interval() { _b[0] = 0;  _b[1] = 0; }
     explicit Interval(Coord u) { _b[0] = _b[1] = u; }
     /* When creating an Interval using the constructor specifying the exact range, the created interval
      * will be [u,v] when u<=v ; and will be [v,u] when v < u !!!
@@ -78,7 +81,8 @@ public:
     inline Coord extent() const { return _b[1] - _b[0]; }
     inline Coord middle() const { return (_b[1] + _b[0]) * 0.5; }
     
-    inline bool isEmpty() const { return _b[0] > _b[1]; }
+//    inline bool isEmpty() const { return _b[0] > _b[1]; }
+    inline bool isSingular() const { return _b[0] == _b[1]; }
     inline bool contains(Coord val) const { return _b[0] <= val && val <= _b[1]; }
     bool contains(const Interval & val) const { return _b[0] <= val._b[0] && val._b[1] <= _b[1]; }
     bool intersects(const Interval & val) const {
@@ -167,9 +171,15 @@ public:
         return result;
     }
     
+    /** When this would create an empty interval, the interval will be the centerpoint of the old range only.
+     */
     inline void expandBy(double amnt) {
         _b[0] -= amnt;
         _b[1] += amnt;
+        if (_b[0] > _b[1]) {
+            Coord halfway = (_b[0]+_b[1])/2;
+            _b[0] = _b[1] = halfway;
+        }
     }
     
     inline void unionWith(const Interval & a) {
@@ -214,12 +224,42 @@ inline Interval unify(const Interval & a, const Interval & b) {
     return Interval(std::min(a.min(), b.min()),
                     std::max(a.max(), b.max()));
 }
-inline boost::optional<Interval> intersect(const Interval & a, const Interval & b) {
+
+/**
+ * \brief OptInterval is an Interval that can be empty.
+ */
+class OptInterval : public boost::optional<Interval> {
+public:
+    OptInterval() : boost::optional<Interval>() {};
+    OptInterval(Interval const &a) : boost::optional<Interval>(a) {};
+    OptInterval(Coord u) : boost::optional<Interval>(Interval(u)) {};
+    OptInterval(Coord u, Coord v) : boost::optional<Interval>(Interval(u,v)) {};
+
+    /**
+     * Check whether this OptInterval is empty or not.
+     */
+    inline bool isEmpty() { return (*this == false); };
+    
+    /**
+     * If \c this is empty, copy argument \c a. Otherwise, union with it (and do nothing when \c a is empty)
+     */
+    inline void unionWith(const OptInterval & a) {
+        if (a) {
+            if (*this) { // check that we are not empty
+                (*this)->unionWith(*a);
+            } else {
+                *this = a;
+            }
+        }
+    }
+};
+
+inline OptInterval intersect(const Interval & a, const Interval & b) {
     Coord u = std::max(a.min(), b.min()),
           v = std::min(a.max(), b.max());
     //technically >= might be incorrect, but singulars suck
-    return u >= v ? boost::optional<Interval>()
-                  : boost::optional<Interval>(Interval(u, v));
+    return u >= v ? OptInterval()
+                  : OptInterval(Interval(u, v));
 }
 
 }
diff --git a/src/2geom/linear.h b/src/2geom/linear.h
index 4594e17be..e8828b283 100644
--- a/src/2geom/linear.h
+++ b/src/2geom/linear.h
@@ -37,6 +37,15 @@
 #include <2geom/interval.h>
 #include <2geom/isnan.h>
 
+
+//#define USE_SBASIS_OF
+
+#ifdef USE_SBASIS_OF
+
+#include "linear-of.h"
+
+#else
+
 namespace Geom{
 
 inline double lerp(double t, double a, double b) { return a*(1-t) + b*t; }
@@ -100,9 +109,9 @@ public:
     //defined in sbasis.h
     inline SBasis toSBasis() const;
 
-    inline Interval bounds_exact() const { return Interval(a[0], a[1]); }
-    inline Interval bounds_fast() const { return bounds_exact(); }
-    inline Interval bounds_local(double u, double v) const { return Interval(valueAt(u), valueAt(v)); }
+    inline OptInterval bounds_exact() const { return Interval(a[0], a[1]); }
+    inline OptInterval bounds_fast() const { return bounds_exact(); }
+    inline OptInterval bounds_local(double u, double v) const { return Interval(valueAt(u), valueAt(v)); }
 
     operator Tri() const {
         return a[1] - a[0];
@@ -169,7 +178,9 @@ inline Linear operator/=(Linear & a, double b) {
     a[0] /= b; a[1] /= b;
     return a;
 }
-};
+
+}
+#endif
 
 #endif //SEEN_LINEAR_H
 
diff --git a/src/2geom/nearest-point.cpp b/src/2geom/nearest-point.cpp
index 705767f13..71beab9cb 100644
--- a/src/2geom/nearest-point.cpp
+++ b/src/2geom/nearest-point.cpp
@@ -50,34 +50,34 @@ namespace Geom
  */
 
 double nearest_point( Point const& p,
-					  D2<SBasis> const& c,
-					  D2<SBasis> const& dc,
-		              double from, double to )
+                      D2<SBasis> const& c,
+                      D2<SBasis> const& dc,
+                      double from, double to )
 {
-	if ( from > to ) std::swap(from, to);
-	if ( from < 0 || to > 1 )
-	{
-		THROW_RANGEERROR("[from,to] interval out of bounds");
-	}
-	if (c.isConstant()) return from;
-	SBasis dd = dot(c - p, dc);
-	std::vector<double> zeros = Geom::roots(dd);
+    if ( from > to ) std::swap(from, to);
+    if ( from < 0 || to > 1 )
+    {
+        THROW_RANGEERROR("[from,to] interval out of bounds");
+    }
+    if (c.isConstant()) return from;
+    SBasis dd = dot(c - p, dc);
+    std::vector<double> zeros = Geom::roots(dd);
 
-	double closest = from;
-	double min_dist_sq = L2sq(c(from) - p);
-	double distsq;
-	for ( unsigned int i = 0; i < zeros.size(); ++i )
-	{
-		distsq = L2sq(c(zeros[i]) - p);
-		if ( min_dist_sq > L2sq(c(zeros[i]) - p) )
-		{
-			closest = zeros[i];
-			min_dist_sq = distsq;
-		}
-	}
-	if ( min_dist_sq > L2sq( c(to) - p ) )
-		closest = to;
-	return closest;
+    double closest = from;
+    double min_dist_sq = L2sq(c(from) - p);
+    double distsq;
+    for ( unsigned int i = 0; i < zeros.size(); ++i )
+    {
+        distsq = L2sq(c(zeros[i]) - p);
+        if ( min_dist_sq > L2sq(c(zeros[i]) - p) )
+        {
+            closest = zeros[i];
+            min_dist_sq = distsq;
+        }
+    }
+    if ( min_dist_sq > L2sq( c(to) - p ) )
+        closest = to;
+    return closest;
 
 }
 
@@ -89,54 +89,54 @@ double nearest_point( Point const& p,
 
 std::vector<double>
 all_nearest_points( Point const& p,
-		    D2<SBasis> const& c,
-		    D2<SBasis> const& dc,
-		    double from, double to )
+            D2<SBasis> const& c,
+            D2<SBasis> const& dc,
+            double from, double to )
 {
-	std::swap(from, to);
-	if ( from > to ) std::swap(from, to);
-	if ( from < 0 || to > 1 )
-	{
-		THROW_RANGEERROR("[from,to] interval out of bounds");
-	}
+    std::swap(from, to);
+    if ( from > to ) std::swap(from, to);
+    if ( from < 0 || to > 1 )
+    {
+        THROW_RANGEERROR("[from,to] interval out of bounds");
+    }
 
-	std::vector<double> result;
-	if (c.isConstant())
-	{
-	    result.push_back(from);
-	    return result;
-	}
-	SBasis dd = dot(c - p, dc);
+    std::vector<double> result;
+    if (c.isConstant())
+    {
+        result.push_back(from);
+        return result;
+    }
+    SBasis dd = dot(c - p, dc);
 
-	std::vector<double> zeros = Geom::roots(dd);
-	std::vector<double> candidates;
-	candidates.push_back(from);
-	candidates.insert(candidates.end(), zeros.begin(), zeros.end());
-	candidates.push_back(to);
-	std::vector<double> distsq;
-	distsq.reserve(candidates.size());
-	for ( unsigned int i = 0; i < candidates.size(); ++i )
-	{
-		distsq.push_back( L2sq(c(candidates[i]) - p) );
-	}
-	unsigned int closest = 0;
-	double dsq = distsq[0];
-	for ( unsigned int i = 1; i < candidates.size(); ++i )
-	{
-		if ( dsq > distsq[i] )
-		{
-			closest = i;
-			dsq = distsq[i];
-		}
-	}
-	for ( unsigned int i = 0; i < candidates.size(); ++i )
-	{
-		if( distsq[closest] == distsq[i] )
-		{
-			result.push_back(candidates[i]);
-		}
-	}
-	return result;
+    std::vector<double> zeros = Geom::roots(dd);
+    std::vector<double> candidates;
+    candidates.push_back(from);
+    candidates.insert(candidates.end(), zeros.begin(), zeros.end());
+    candidates.push_back(to);
+    std::vector<double> distsq;
+    distsq.reserve(candidates.size());
+    for ( unsigned int i = 0; i < candidates.size(); ++i )
+    {
+        distsq.push_back( L2sq(c(candidates[i]) - p) );
+    }
+    unsigned int closest = 0;
+    double dsq = distsq[0];
+    for ( unsigned int i = 1; i < candidates.size(); ++i )
+    {
+        if ( dsq > distsq[i] )
+        {
+            closest = i;
+            dsq = distsq[i];
+        }
+    }
+    for ( unsigned int i = 0; i < candidates.size(); ++i )
+    {
+        if( distsq[closest] == distsq[i] )
+        {
+            result.push_back(candidates[i]);
+        }
+    }
+    return result;
 }
 
 
@@ -145,141 +145,141 @@ all_nearest_points( Point const& p,
 
 
 double nearest_point( Point const& p,
-					  Piecewise< D2<SBasis> > const& c,
-		              double from, double to )
+                      Piecewise< D2<SBasis> > const& c,
+                      double from, double to )
 {
-	if ( from > to ) std::swap(from, to);
-	if ( from < c.cuts[0] || to > c.cuts[c.size()] )
-	{
-		THROW_RANGEERROR("[from,to] interval out of bounds");
-	}
+    if ( from > to ) std::swap(from, to);
+    if ( from < c.cuts[0] || to > c.cuts[c.size()] )
+    {
+        THROW_RANGEERROR("[from,to] interval out of bounds");
+    }
 
-	unsigned int si = c.segN(from);
-	unsigned int ei = c.segN(to);
-	if ( si == ei )
-	{
-		double nearest=
-			nearest_point(p, c[si], c.segT(from, si), c.segT(to, si));
-		return c.mapToDomain(nearest, si);
-	}
-	double t;
-	double nearest = nearest_point(p, c[si], c.segT(from, si));
-	unsigned int ni = si;
-	double dsq;
-	double mindistsq = distanceSq(p, c[si](nearest));
-	Rect bb;
-	for ( unsigned int i = si + 1; i < ei; ++i )
-	{
-		bb = bounds_fast(c[i]);
-		dsq = distanceSq(p, bb);
-		if ( mindistsq <= dsq ) continue;
-		t = nearest_point(p, c[i]);
-		dsq = distanceSq(p, c[i](t));
-		if ( mindistsq > dsq )
-		{
-			nearest = t;
-			ni = i;
-			mindistsq = dsq;
-		}
-	}
-	bb = bounds_fast(c[ei]);
-	dsq = distanceSq(p, bb);
-	if ( mindistsq > dsq )
-	{
-		t = nearest_point(p, c[ei], 0, c.segT(to, ei));
-		dsq = distanceSq(p, c[ei](t));
-		if ( mindistsq > dsq )
-		{
-			nearest = t;
-			ni = ei;
-		}
-	}
-	return c.mapToDomain(nearest, ni);
+    unsigned int si = c.segN(from);
+    unsigned int ei = c.segN(to);
+    if ( si == ei )
+    {
+        double nearest=
+            nearest_point(p, c[si], c.segT(from, si), c.segT(to, si));
+        return c.mapToDomain(nearest, si);
+    }
+    double t;
+    double nearest = nearest_point(p, c[si], c.segT(from, si));
+    unsigned int ni = si;
+    double dsq;
+    double mindistsq = distanceSq(p, c[si](nearest));
+    Rect bb(Geom::Point(0,0),Geom::Point(0,0));
+    for ( unsigned int i = si + 1; i < ei; ++i )
+    {
+        bb = *bounds_fast(c[i]);
+        dsq = distanceSq(p, bb);
+        if ( mindistsq <= dsq ) continue;
+        t = nearest_point(p, c[i]);
+        dsq = distanceSq(p, c[i](t));
+        if ( mindistsq > dsq )
+        {
+            nearest = t;
+            ni = i;
+            mindistsq = dsq;
+        }
+    }
+    bb = *bounds_fast(c[ei]);
+    dsq = distanceSq(p, bb);
+    if ( mindistsq > dsq )
+    {
+        t = nearest_point(p, c[ei], 0, c.segT(to, ei));
+        dsq = distanceSq(p, c[ei](t));
+        if ( mindistsq > dsq )
+        {
+            nearest = t;
+            ni = ei;
+        }
+    }
+    return c.mapToDomain(nearest, ni);
 }
 
 std::vector<double>
 all_nearest_points( Point const& p,
                     Piecewise< D2<SBasis> > const& c,
-		            double from, double to )
+                    double from, double to )
 {
-	if ( from > to ) std::swap(from, to);
-	if ( from < c.cuts[0] || to > c.cuts[c.size()] )
-	{
-		THROW_RANGEERROR("[from,to] interval out of bounds");
-	}
+    if ( from > to ) std::swap(from, to);
+    if ( from < c.cuts[0] || to > c.cuts[c.size()] )
+    {
+        THROW_RANGEERROR("[from,to] interval out of bounds");
+    }
 
-	unsigned int si = c.segN(from);
-	unsigned int ei = c.segN(to);
-	if ( si == ei )
-	{
-		std::vector<double>	all_nearest =
-			all_nearest_points(p, c[si], c.segT(from, si), c.segT(to, si));
-		for ( unsigned int i = 0; i < all_nearest.size(); ++i )
-		{
-			all_nearest[i] = c.mapToDomain(all_nearest[i], si);
-		}
-		return all_nearest;
-	}
-	std::vector<double> all_t;
-	std::vector< std::vector<double> > all_np;
-	all_np.push_back( all_nearest_points(p, c[si], c.segT(from, si)) );
-	std::vector<unsigned int> ni;
-	ni.push_back(si);
-	double dsq;
-	double mindistsq = distanceSq( p, c[si](all_np.front().front()) );
-	Rect bb;
-	for ( unsigned int i = si + 1; i < ei; ++i )
-	{
-		bb = bounds_fast(c[i]);
-		dsq = distanceSq(p, bb);
-		if ( mindistsq < dsq ) continue;
-		all_t = all_nearest_points(p, c[i]);
-		dsq = distanceSq( p, c[i](all_t.front()) );
-		if ( mindistsq > dsq )
-		{
-			all_np.clear();
-			all_np.push_back(all_t);
-			ni.clear();
-			ni.push_back(i);
-			mindistsq = dsq;
-		}
-		else if ( mindistsq == dsq )
-		{
-			all_np.push_back(all_t);
-			ni.push_back(i);
-		}
-	}
-	bb = bounds_fast(c[ei]);
-	dsq = distanceSq(p, bb);
-	if ( mindistsq >= dsq )
-	{
-		all_t = all_nearest_points(p, c[ei], 0, c.segT(to, ei));
-		dsq = distanceSq( p, c[ei](all_t.front()) );
-		if ( mindistsq > dsq )
-		{
-			for ( unsigned int i = 0; i < all_t.size(); ++i )
-			{
-				all_t[i] = c.mapToDomain(all_t[i], ei);
-			}
-			return all_t;
-		}
-		else if ( mindistsq == dsq )
-		{
-			all_np.push_back(all_t);
-			ni.push_back(ei);
-		}
-	}
-	std::vector<double> all_nearest;
-	for ( unsigned int i = 0; i < all_np.size(); ++i )
-	{
-		for ( unsigned int j = 0; j < all_np[i].size(); ++j )
-		{
-			all_nearest.push_back( c.mapToDomain(all_np[i][j], ni[i]) );
-		}
-	}
-	all_nearest.erase(std::unique(all_nearest.begin(), all_nearest.end()),
-	                  all_nearest.end());
-	return all_nearest;
+    unsigned int si = c.segN(from);
+    unsigned int ei = c.segN(to);
+    if ( si == ei )
+    {
+        std::vector<double>	all_nearest =
+            all_nearest_points(p, c[si], c.segT(from, si), c.segT(to, si));
+        for ( unsigned int i = 0; i < all_nearest.size(); ++i )
+        {
+            all_nearest[i] = c.mapToDomain(all_nearest[i], si);
+        }
+        return all_nearest;
+    }
+    std::vector<double> all_t;
+    std::vector< std::vector<double> > all_np;
+    all_np.push_back( all_nearest_points(p, c[si], c.segT(from, si)) );
+    std::vector<unsigned int> ni;
+    ni.push_back(si);
+    double dsq;
+    double mindistsq = distanceSq( p, c[si](all_np.front().front()) );
+    Rect bb(Geom::Point(0,0),Geom::Point(0,0));
+    for ( unsigned int i = si + 1; i < ei; ++i )
+    {
+        bb = *bounds_fast(c[i]);
+        dsq = distanceSq(p, bb);
+        if ( mindistsq < dsq ) continue;
+        all_t = all_nearest_points(p, c[i]);
+        dsq = distanceSq( p, c[i](all_t.front()) );
+        if ( mindistsq > dsq )
+        {
+            all_np.clear();
+            all_np.push_back(all_t);
+            ni.clear();
+            ni.push_back(i);
+            mindistsq = dsq;
+        }
+        else if ( mindistsq == dsq )
+        {
+            all_np.push_back(all_t);
+            ni.push_back(i);
+        }
+    }
+    bb = *bounds_fast(c[ei]);
+    dsq = distanceSq(p, bb);
+    if ( mindistsq >= dsq )
+    {
+        all_t = all_nearest_points(p, c[ei], 0, c.segT(to, ei));
+        dsq = distanceSq( p, c[ei](all_t.front()) );
+        if ( mindistsq > dsq )
+        {
+            for ( unsigned int i = 0; i < all_t.size(); ++i )
+            {
+                all_t[i] = c.mapToDomain(all_t[i], ei);
+            }
+            return all_t;
+        }
+        else if ( mindistsq == dsq )
+        {
+            all_np.push_back(all_t);
+            ni.push_back(ei);
+        }
+    }
+    std::vector<double> all_nearest;
+    for ( unsigned int i = 0; i < all_np.size(); ++i )
+    {
+        for ( unsigned int j = 0; j < all_np[i].size(); ++j )
+        {
+            all_nearest.push_back( c.mapToDomain(all_np[i][j], ni[i]) );
+        }
+    }
+    all_nearest.erase(std::unique(all_nearest.begin(), all_nearest.end()),
+                      all_nearest.end());
+    return all_nearest;
 }
 
 } // end namespace Geom
diff --git a/src/2geom/path-intersection.cpp b/src/2geom/path-intersection.cpp
index 715c67c23..5a4e76020 100644
--- a/src/2geom/path-intersection.cpp
+++ b/src/2geom/path-intersection.cpp
@@ -23,9 +23,9 @@ int winding(Path const &path, Point p) {
   Path::const_iterator start;
   for(Path::const_iterator iter = path.begin(); ; ++iter) {
     if(iter == path.end_closed()) { return 0; }
-    if(iter->initialPoint()[Y]!=p[Y])  { start = iter; break; }
-    if(iter->finalPoint()[Y]!=p[Y])    { start = iter; break; }
-    if(iter->boundsFast().height()!=0.){ start = iter; break; }
+    if(iter->initialPoint()[Y]!=p[Y])   { start = iter; break; }
+    if(iter->finalPoint()[Y]!=p[Y])     { start = iter; break; }
+    if(iter->boundsFast()->height()!=0.){ start = iter; break; }
   }
   int wind = 0;
   unsigned cnt = 0;
@@ -36,7 +36,7 @@ int winding(Path const &path, Point p) {
     cnt++;
     if(cnt > path.size()) return wind;  //some bug makes this required
     starting = false;
-    Rect bounds = iter->boundsFast();
+    Rect bounds = *(iter->boundsFast());
     Coord x = p[X], y = p[Y];
     
     if(x > bounds.right() || !bounds[Y].contains(y)) continue; //ray doesn't intersect box
@@ -70,7 +70,7 @@ int winding(Path const &path, Point p) {
         next++;
         for(; ; next++) {
             if(next == path.end_closed()) next = path.begin();
-            Rect bnds = next->boundsFast();
+            Rect bnds = *(next->boundsFast());
             //TODO: X considerations
             if(bnds.height() > 0) {
                 //It has diverged
@@ -267,13 +267,13 @@ void pair_intersect(Curve const & A, double Al, double Ah,
                     Curve const & B, double Bl, double Bh,
                     Crossings &ret,  unsigned depth=0) {
    // std::cout << depth << "(" << Al << ", " << Ah << ")\n";
-    Rect Ar = A.boundsLocal(Interval(Al, Ah));
-    if(Ar.isEmpty()) return;
+    OptRect Ar = A.boundsLocal(Interval(Al, Ah));
+    if (!Ar) return;
 
-    Rect Br = B.boundsLocal(Interval(Bl, Bh));
-    if(Br.isEmpty()) return;
+    OptRect Br = B.boundsLocal(Interval(Bl, Bh));
+    if (!Br) return;
     
-    if(!Ar.intersects(Br)) return;
+    if(! Ar->intersects(*Br)) return;
     
     //Checks the general linearity of the function
     if((depth > 12)) { // || (A.boundsLocal(Interval(Al, Ah), 1).maxExtent() < 0.1 
diff --git a/src/2geom/path.cpp b/src/2geom/path.cpp
index c04d9d08d..136e6d481 100644
--- a/src/2geom/path.cpp
+++ b/src/2geom/path.cpp
@@ -43,21 +43,21 @@ using namespace Geom::PathInternal;
 namespace Geom
 {
 
-Rect Path::boundsFast() const {
-  Rect bounds;
+OptRect Path::boundsFast() const {
+  OptRect bounds;
   if (empty()) return bounds;
   bounds=front().boundsFast();
   const_iterator iter = begin();
   if ( iter != end() ) {
-	  for ( ++iter; iter != end() ; ++iter ) {
-	    bounds.unionWith(iter->boundsFast());
-	  }
+    for ( ++iter; iter != end() ; ++iter ) {
+      bounds.unionWith(iter->boundsFast());
+    }
   }
   return bounds;
 }
 
-Rect Path::boundsExact() const {
-  Rect bounds;
+OptRect Path::boundsExact() const {
+  OptRect bounds;
   if (empty()) return bounds;
   bounds=front().boundsExact();
   const_iterator iter = begin();
@@ -143,10 +143,10 @@ Path::allNearestPoints(Point const& _point, double from, double to) const
 	double dsq;
 	double mindistsq
 		= distanceSq( _point, _path[si].pointAt( all_np.front().front() ) );
-	Rect bb;
+	Rect bb(Geom::Point(0,0),Geom::Point(0,0));
 	for ( unsigned int i = si + 1; i < ei; ++i )
 	{
-		bb = _path[i].boundsFast();
+		bb = *(_path[i].boundsFast());
 		dsq = distanceSq(_point, bb);
 		if ( mindistsq < dsq ) continue;
 		all_t = _path[i].allNearestPoints(_point);
@@ -165,7 +165,7 @@ Path::allNearestPoints(Point const& _point, double from, double to) const
 			ni.push_back(i);
 		}
 	}
-	bb = _path[ei].boundsFast();
+	bb = *(_path[ei].boundsFast());
 	dsq = distanceSq(_point, bb);
 	if ( mindistsq >= dsq )
 	{
@@ -254,10 +254,10 @@ double Path::nearestPoint(Point const &_point, double from, double to, double *d
 	unsigned int ni = si;
 	double dsq;
 	double mindistsq = distanceSq(_point, _path[si].pointAt(nearest));
-	Rect bb;
+	Rect bb(Geom::Point(0,0),Geom::Point(0,0));
 	for ( unsigned int i = si + 1; i < ei; ++i )
 	{
-		bb = _path[i].boundsFast();
+		bb = *(_path[i].boundsFast());
 		dsq = distanceSq(_point, bb);
 		if ( mindistsq <= dsq ) continue;
 		t = _path[i].nearestPoint(_point);
@@ -269,7 +269,7 @@ double Path::nearestPoint(Point const &_point, double from, double to, double *d
 			mindistsq = dsq;
 		}
 	}
-	bb = _path[ei].boundsFast();
+	bb = *(_path[ei].boundsFast());
 	dsq = distanceSq(_point, bb);
 	if ( mindistsq > dsq )
 	{
diff --git a/src/2geom/path.h b/src/2geom/path.h
index e02b749e7..b95a54eaa 100644
--- a/src/2geom/path.h
+++ b/src/2geom/path.h
@@ -259,8 +259,8 @@ public:
   bool closed() const { return closed_; }
   void close(bool closed=true) { closed_ = closed; }
 
-  Rect boundsFast() const;
-  Rect boundsExact() const;
+  OptRect boundsFast() const;
+  OptRect boundsExact() const;
 
   Piecewise<D2<SBasis> > toPwSb() const {
     Piecewise<D2<SBasis> > ret;
@@ -566,8 +566,9 @@ public:
   }
 
 
-  /*
+  /**
    * It is important to note that the coordinates passed to appendNew should be finite!
+   * If one of the coordinates is infinite, 2geom will throw a ContinuityError exception.
    */
 
   template <typename CurveType, typename A>
diff --git a/src/2geom/pathvector.cpp b/src/2geom/pathvector.cpp
index 39a00d8dc..790265c76 100644
--- a/src/2geom/pathvector.cpp
+++ b/src/2geom/pathvector.cpp
@@ -56,12 +56,11 @@ PathVector reverse_paths_and_order (PathVector const & path_in)
     return path_out;
 }
 
-// FIXME: this function does not work for empty paths, because Rect cannot be initialized empty yet. check rect.h
-Rect bounds_fast( PathVector const& pv )
+OptRect bounds_fast( PathVector const& pv )
 {
     typedef PathVector::const_iterator const_iterator;
     
-    Rect bound;
+    OptRect bound;
     if (pv.empty()) return bound;
     
     bound = (pv.begin())->boundsFast();
@@ -72,12 +71,11 @@ Rect bounds_fast( PathVector const& pv )
     return bound;
 }
 
-// FIXME: this function does not work for empty paths, because Rect cannot be initialized empty yet. check rect.h
-Rect bounds_exact( PathVector const& pv )
+OptRect bounds_exact( PathVector const& pv )
 {
     typedef PathVector::const_iterator const_iterator;
     
-    Rect bound;
+    OptRect bound;
     if (pv.empty()) return bound;
     
     bound = (pv.begin())->boundsExact();
diff --git a/src/2geom/pathvector.h b/src/2geom/pathvector.h
index 5df3ad00c..9efae7c73 100644
--- a/src/2geom/pathvector.h
+++ b/src/2geom/pathvector.h
@@ -102,8 +102,8 @@ Geom::Point finalPoint(PathVector const &path_in)
 
 PathVector reverse_paths_and_order (PathVector const & path_in);
 
-Rect bounds_fast( PathVector const & pv );
-Rect bounds_exact( PathVector const & pv );
+OptRect bounds_fast( PathVector const & pv );
+OptRect bounds_exact( PathVector const & pv );
 
 struct PathVectorPosition {
     // pathvector[path_nr].pointAt(t) is the position
diff --git a/src/2geom/piecewise.h b/src/2geom/piecewise.h
index ec3ba26ce..bbd1f054a 100644
--- a/src/2geom/piecewise.h
+++ b/src/2geom/piecewise.h
@@ -80,6 +80,8 @@ class Piecewise {
         push_cut(1.);
     }
 
+    unsigned input_dim(){return 1;}
+
     typedef typename T::output_type output_type;
 
     explicit Piecewise(const output_type & v) {
@@ -198,14 +200,18 @@ class Piecewise {
     //Transforms the domain into another interval
     inline void setDomain(Interval dom) {
         if(empty()) return;
+        /* dom can not be empty
         if(dom.isEmpty()) {
             cuts.clear(); segs.clear();
             return;
-        }
+        }*/
         double cf = cuts.front();
         double o = dom.min() - cf, s = dom.extent() / (cuts.back() - cf);
         for(unsigned i = 0; i <= size(); i++)
             cuts[i] = (cuts[i] - cf) * s + o;
+        //fix floating point precision errors.
+        cuts[0] = dom.min();
+        cuts[size()] = dom.max();
     }
 
     //Concatenates this Piecewise function with another, offseting time of the other to match the end.
@@ -227,15 +233,15 @@ class Piecewise {
     inline void continuousConcat(const Piecewise<T> &other) {
         boost::function_requires<AddableConcept<typename T::output_type> >();
         if(other.empty()) return;
-        typename T::output_type y = segs.back().at1() - other.segs.front().at0();
 
         if(empty()) {
             for(unsigned i = 0; i < other.size(); i++)
-                push_seg(other[i] + y);
+                push_seg(other[i]);
             cuts = other.cuts;
             return;
         }
 
+        typename T::output_type y = segs.back().at1() - other.segs.front().at0();
         double t = cuts.back() - other.cuts.front();
         for(unsigned i = 0; i < other.size(); i++)
             push(other[i] + y, other.cuts[i + 1] + t);
@@ -278,11 +284,12 @@ inline typename FragmentConcept<T>::BoundsType bounds_exact(const Piecewise<T> &
 }
 
 template<typename T>
-inline typename FragmentConcept<T>::BoundsType bounds_local(const Piecewise<T> &f, const Interval &m) {
+inline typename FragmentConcept<T>::BoundsType bounds_local(const Piecewise<T> &f, const OptInterval &_m) {
     boost::function_requires<FragmentConcept<T> >();
 
-    if(f.empty()) return typename FragmentConcept<T>::BoundsType();
-    if(m.isEmpty()) return typename FragmentConcept<T>::BoundsType(f(m.min()));
+    if(f.empty() || !_m) return typename FragmentConcept<T>::BoundsType();
+    Interval const &m = *_m;
+    if(m.isSingular()) return typename FragmentConcept<T>::BoundsType(f(m.min()));
 
     unsigned fi = f.segN(m.min()), ti = f.segN(m.max());
     double ft = f.segT(m.min(), fi), tt = f.segT(m.max(), ti);
@@ -640,7 +647,7 @@ Piecewise<T> compose(Piecewise<T> const &f, SBasis const &g){
     }
 
     //first check bounds...
-    Interval bs = bounds_fast(g);
+    Interval bs = *bounds_fast(g);
     if (f.cuts.front() > bs.max()  || bs.min() > f.cuts.back()){
         int idx = (bs.max() < f.cuts[1]) ? 0 : f.cuts.size()-2;
         double t0 = f.cuts[idx], width = f.cuts[idx+1] - t0;
diff --git a/src/2geom/point.h b/src/2geom/point.h
index e6e74242d..d89b53f83 100644
--- a/src/2geom/point.h
+++ b/src/2geom/point.h
@@ -9,6 +9,7 @@
 #include <iostream>
 
 #include <2geom/coord.h>
+#include <2geom/isnan.h>
 #include <2geom/utils.h>
 
 namespace Geom {
@@ -22,6 +23,7 @@ class Point {
     Coord _pt[2];
 
   public:
+    /// The default constructor creates an Point(0,0)  DO NOT RELY ON THIS, BEST NOT TO USE THIS CONSTRUCTOR
     inline Point()
     { _pt[X] = _pt[Y] = 0; }
 
@@ -81,6 +83,13 @@ class Point {
 
     void normalize();
 
+    inline bool isFinite() const {
+        for ( unsigned i = 0 ; i < 2 ; ++i ) {
+            if(!IS_FINITE(_pt[i])) return false;
+        }
+        return true;
+    }
+
     inline Point operator+(Point const &o) const {
         return Point(_pt[X] + o._pt[X], _pt[Y] + o._pt[Y]);
     }
diff --git a/src/2geom/rect.h b/src/2geom/rect.h
index 563c8a058..fb42ff92d 100644
--- a/src/2geom/rect.h
+++ b/src/2geom/rect.h
@@ -48,6 +48,7 @@
 namespace Geom {
 /** D2<Interval> specialization to Rect */
 typedef D2<Interval> Rect;
+class OptRect;
 
 Rect unify(const Rect &, const Rect &);
 /**
@@ -60,10 +61,12 @@ class D2<Interval> {
   private:
     Interval f[2];
   public:
-    /* The default constructor creates an empty rect, constructed of two empty Intervals. (users rely on this!)
+    /** Best not to use this constructor, do not rely on what it initializes the object to.
+     *The default constructor creates a rect of default intervals.
      */
     D2<Interval>() { f[X] = f[Y] = Interval(); }
     
+    public:
     D2<Interval>(Interval const &a, Interval const &b) {
         f[X] = a;
         f[Y] = b;
@@ -105,14 +108,20 @@ class D2<Interval> {
     inline Point midpoint() const { return Point(f[X].middle(), f[Y].middle()); }
 
 /**
- * Compute the area of this rectangle.  Note that a zero area rectangle is not necessarily empty - just as the interval [0,0] contains one point, the rectangle [0,0] x [0,0] contains 1 point and no area.
+ * \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.
  */
     inline double area() const { return f[X].extent() * f[Y].extent(); }
+    inline bool hasZeroArea(double eps = EPSILON) const { return (area() <= eps); }
+
     inline double maxExtent() const { return std::max(f[X].extent(), f[Y].extent()); }
+    inline double minExtent() const { return std::min(f[X].extent(), f[Y].extent()); }
 
-    inline bool isEmpty()                 const { 
-        return f[X].isEmpty()        || f[Y].isEmpty(); 
-    }
+//    inline bool isEmpty()                 const { 
+//        return f[X].isEmpty()        || f[Y].isEmpty(); 
+//    }
     inline bool intersects(Rect const &r) const { 
         return f[X].intersects(r[X]) && f[Y].intersects(r[Y]); 
     }
@@ -129,40 +138,20 @@ class D2<Interval> {
     inline void unionWith(Rect const &b) { 
         f[X].unionWith(b[X]); f[Y].unionWith(b[Y]); 
     }
+    void unionWith(OptRect const &b);
 
-    inline void expandBy(double amnt)    { 
+    inline void expandBy(double amnt)    {
         f[X].expandBy(amnt);  f[Y].expandBy(amnt); 
     }
     inline void expandBy(Point const p)  { 
         f[X].expandBy(p[X]);  f[Y].expandBy(p[Y]); 
     }
-    
-    /** Transforms the rect by m. Note that it gives correct results only for scales and translates,
-        in the case of rotations, the area of the rect will grow as it cannot rotate. */
-    inline Rect operator*(Matrix const m) const { 
-        return unify(Rect(corner(0) * m, corner(2) * m),
-                     Rect(corner(1) * m, corner(3) * m));
-    }
 };
 
 inline Rect unify(Rect const & a, Rect const & b) {
     return Rect(unify(a[X], b[X]), unify(a[Y], b[Y]));
 }
 
-/** 
- * Returns the smallest rectangle that encloses both rectangles.
- * An empty argument is assumed to be an empty rectangle
- */
-inline boost::optional<Rect> unify(boost::optional<Rect> const & a, boost::optional<Rect> const & b) {
-    if (!a) {
-        return b;
-    } else if (!b) {
-        return a;
-    } else {
-        return unify(*a, *b);
-    }
-}
-
 inline Rect union_list(std::vector<Rect> const &r) {
     if(r.empty()) return Rect(Interval(0,0), Interval(0,0));
     Rect ret = r[0];
@@ -171,12 +160,6 @@ inline Rect union_list(std::vector<Rect> const &r) {
     return ret;
 }
 
-inline boost::optional<Rect> intersect(Rect const & a, Rect const & b) {
-    boost::optional<Interval> x = intersect(a[X], b[X]);
-    boost::optional<Interval> y = intersect(a[Y], b[Y]);
-    return x && y ? boost::optional<Rect>(Rect(*x, *y)) : boost::optional<Rect>();
-}
-
 inline
 double distanceSq( Point const& p, Rect const& rect )
 {
@@ -209,6 +192,70 @@ double distance( Point const& p, Rect const& rect )
     return std::sqrt(distanceSq(p, rect));
 }
 
+/**
+ * The OptRect class can represent and empty Rect and non-empty Rects.
+ * If OptRect is not empty, it means that both X and Y intervals are not empty.
+ * 
+ */
+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.
+    }
+
+    /**
+     * Check whether this OptRect is empty or not.
+     */
+    inline bool isEmpty() { return (*this == false); };
+
+    /**
+     * If \c this is empty, copy argument \c b. Otherwise, union with it (and do nothing when \c b is empty)
+     */
+    inline void unionWith(OptRect const &b) {
+        if (b) {
+            if (*this) { // check that we are not empty
+                (**this)[X].unionWith((*b)[X]);
+                (**this)[Y].unionWith((*b)[Y]);
+            } 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 OptRect intersect(Rect const & a, Rect const & b) {
+    return OptRect(intersect(a[X], b[X]), intersect(a[Y], b[Y]));
+}
+
+inline void Rect::unionWith(OptRect const &b) { 
+    if (b) {
+        unionWith(*b);
+    }
+}
 
 } // end namespace Geom
 
diff --git a/src/2geom/region.h b/src/2geom/region.h
index 7b2f5763d..fe2517e23 100644
--- a/src/2geom/region.h
+++ b/src/2geom/region.h
@@ -48,14 +48,14 @@ class Region {
     friend Shape shape_boolean(bool rev, Shape const & a, Shape const & b, CrossingSet const & crs);
 
     Path boundary;
-    mutable boost::optional<Rect> box;
+    mutable OptRect box;
     bool fill;
   public:
     Region() : fill(true) {}
     explicit Region(Path const &p) : boundary(p) { fill = path_direction(p); }
     Region(Path const &p, bool dir) : boundary(p), fill(dir) {}
-    Region(Path const &p, boost::optional<Rect> const &b) : boundary(p), box(b) { fill = path_direction(p); }
-    Region(Path const &p, boost::optional<Rect> const &b, bool dir) : boundary(p), box(b), fill(dir) {}
+    Region(Path const &p, OptRect const &b) : boundary(p), box(b) { fill = path_direction(p); }
+    Region(Path const &p, OptRect const &b, bool dir) : boundary(p), box(b), fill(dir) {}
     
     unsigned size() const { return boundary.size(); }
     
@@ -65,7 +65,7 @@ class Region {
     
     operator Path() const { return boundary; }
     Rect boundsFast() const {
-        if(!box) box = boost::optional<Rect>(boundary.boundsFast());
+        if(!box) box = boundary.boundsFast();  /// \todo this doesn't look right at all...
         return *box;
     }
     
diff --git a/src/2geom/sbasis-2d.cpp b/src/2geom/sbasis-2d.cpp
index 9566e0a19..e3a407094 100644
--- a/src/2geom/sbasis-2d.cpp
+++ b/src/2geom/sbasis-2d.cpp
@@ -173,7 +173,7 @@ sb2d_cubic_solve(SBasis2d const &f, Geom::Point const &A, Geom::Point const &B){
     double error = -1;
     unsigned best = 0;
     for (unsigned i=0; i<candidates.size(); i++){
-        Interval bounds = bounds_fast(compose(f,candidates[i]));
+        Interval bounds = *bounds_fast(compose(f,candidates[i]));
         double new_error = (fabs(bounds.max())>fabs(bounds.min()) ? fabs(bounds.max()) : fabs(bounds.min()) );
         if ( new_error < error || error < 0 ){
             error = new_error;
diff --git a/src/2geom/sbasis-curve.h b/src/2geom/sbasis-curve.h
index e11919401..b45e63eb8 100644
--- a/src/2geom/sbasis-curve.h
+++ b/src/2geom/sbasis-curve.h
@@ -68,9 +68,9 @@ public:
   void setInitial(Point v) { for(unsigned d = 0; d < 2; d++) { inner[d][0][0] = v[d]; } }
   void setFinal(Point v)   { for(unsigned d = 0; d < 2; d++) { inner[d][0][1] = v[d]; } }
 
-  Rect boundsFast() const  { return bounds_fast(inner); }
-  Rect boundsExact() const { return bounds_exact(inner); }
-  Rect boundsLocal(Interval i, unsigned deg) const { return bounds_local(inner, i, deg); }
+  virtual OptRect boundsFast() const  { return bounds_fast(inner); }
+  virtual OptRect boundsExact() const { return bounds_exact(inner); }
+  virtual OptRect boundsLocal(OptInterval i, unsigned deg) const { return bounds_local(inner, i, deg); }
 
   std::vector<double> roots(double v, Dim2 d) const { return Geom::roots(inner[d] - v); }
   
diff --git a/src/2geom/sbasis-geometric.cpp b/src/2geom/sbasis-geometric.cpp
index f0170ec6b..96a5ed0ce 100644
--- a/src/2geom/sbasis-geometric.cpp
+++ b/src/2geom/sbasis-geometric.cpp
@@ -528,7 +528,7 @@ unsigned Geom::centroid(Piecewise<D2<SBasis> > const &p, Point& centroid, double
  * Below are basic functions dedicated to solving this assuming a0 and a1 !=0.
  */
 
-static Interval
+static OptInterval
 find_bounds_for_lambda0(double aa0,double aa1,double cc0,double cc1,
     int insist_on_speeds_signs){
 
@@ -541,7 +541,7 @@ find_bounds_for_lambda0(double aa0,double aa1,double cc0,double cc1,
     double c = (c0<c1 ? c0 : c1);
     double delta = 1-4*a*c;
     if ( delta < 0 )
-        return Interval();//return empty interval
+        return OptInterval();//return empty interval
     double lambda_max = (1+std::sqrt(delta))/2/a;
     
     result = Interval(c,lambda_max);
@@ -549,9 +549,10 @@ find_bounds_for_lambda0(double aa0,double aa1,double cc0,double cc1,
         result *= -1;
     if (insist_on_speeds_signs == 1){
         if (result.max() < 0)//Caution: setMin with max<new min...
-            return Interval();//return empty interval
+            return OptInterval();//return empty interval
         result.setMin(0);
     }
+    result = Interval(result.min()-.1,result.max()+.1);//just in case all our approx. were exact...
     return result;
 }
 
@@ -565,13 +566,13 @@ solve_lambda0(double a0,double a1,double c0,double c1,
     p.push_back(Linear( -a1*a0*(a0+2*c0), -a1*a0*(3*a0+2*c0) ));
     p.push_back(Linear( a1*a0*a0 ));
 
-    Interval domain = find_bounds_for_lambda0(a0,a1,c0,c1,insist_on_speeds_signs);
-    if ( domain.isEmpty() ) 
+    OptInterval domain = find_bounds_for_lambda0(a0,a1,c0,c1,insist_on_speeds_signs);
+    if ( !domain ) 
         return std::vector<double>();
-    p = compose(p,Linear(domain.min(),domain.max()));
+    p = compose(p,Linear(domain->min(),domain->max()));
     std::vector<double>rts = roots(p);
     for (unsigned i=0; i<rts.size(); i++){
-        rts[i] = domain.min()+rts[i]*domain.extent();
+        rts[i] = domain->min() + rts[i] * domain->extent();
     }
     return rts;
 }
diff --git a/src/2geom/sbasis-math.cpp b/src/2geom/sbasis-math.cpp
index 1d179a563..e08023e84 100644
--- a/src/2geom/sbasis-math.cpp
+++ b/src/2geom/sbasis-math.cpp
@@ -41,6 +41,7 @@
 namespace Geom {
 
 
+#include <2geom/d2-sbasis.h>
 #include <stdio.h>
 #include <math.h>
 
@@ -327,18 +328,44 @@ Piecewise<SBasis> reciprocalOnDomain(Interval range, double tol){
 }
 
 Piecewise<SBasis> reciprocal(SBasis const &f, double tol, int order){
-    Piecewise<SBasis> reciprocal_fn=reciprocalOnDomain(bounds_fast(f), tol);
+    Piecewise<SBasis> reciprocal_fn=reciprocalOnDomain(*bounds_fast(f), tol);
     Piecewise<SBasis> result=compose(reciprocal_fn,f);
     truncateResult(result,order);
     return(result);
 }
 Piecewise<SBasis> reciprocal(Piecewise<SBasis> const &f, double tol, int order){
-    Piecewise<SBasis> reciprocal_fn=reciprocalOnDomain(bounds_fast(f), tol);
+    Piecewise<SBasis> reciprocal_fn=reciprocalOnDomain(*bounds_fast(f), tol);
     Piecewise<SBasis> result=compose(reciprocal_fn,f);
     truncateResult(result,order);
     return(result);
 }
 
+/**
+ * \brief Retruns a Piecewise SBasis with prescribed values at prescribed times.
+ * 
+ * \param times: vector of times at which the values are given. Should be sorted in increasing order.
+ * \param values: vector of prescribed values. Should have the same size as times and be sorted accordingly.
+ * \param smoothness: (defaults to 1) regularity class of the result: 0=piecewise linear, 1=continuous derivative, etc...
+ */
+Piecewise<SBasis> interpolate(std::vector<double> times, std::vector<double> values, unsigned smoothness){
+    assert ( values.size() == times.size() );
+    if ( values.size() == 0 ) return Piecewise<SBasis>();
+    if ( values.size() == 1 ) return Piecewise<SBasis>(values[0]);//what about time??
+
+    SBasis sk = shift(Linear(1.),smoothness);
+    SBasis bump_in = integral(sk);
+    bump_in -= bump_in.at0();
+    bump_in /= bump_in.at1();
+    SBasis bump_out = reverse( bump_in );
+    
+    Piecewise<SBasis> result;
+    result.cuts.push_back(times[0]);
+    for (unsigned i = 0; i<values.size()-1; i++){
+        result.push(bump_out*values[i]+bump_in*values[i+1],times[i+1]);
+    }
+    return result;
+}
+
 }
 
 /*
diff --git a/src/2geom/sbasis-math.h b/src/2geom/sbasis-math.h
index 053c2d285..49ad965d4 100644
--- a/src/2geom/sbasis-math.h
+++ b/src/2geom/sbasis-math.h
@@ -42,6 +42,7 @@
 
 #include <2geom/sbasis.h>
 #include <2geom/piecewise.h>
+#include <2geom/d2.h>
 
 namespace Geom{
 //-|x|---------------------------------------------------------------
@@ -81,6 +82,8 @@ Piecewise<SBasis> reciprocalOnDomain(Interval range, double tol=1e-3);
 Piecewise<SBasis> reciprocal(          SBasis const &f, double tol=1e-3, int order=3);
 Piecewise<SBasis> reciprocal(Piecewise<SBasis>const &f, double tol=1e-3, int order=3);
 
+//--interpolate------------------------------------------------------------
+Piecewise<SBasis> interpolate( std::vector<double> times, std::vector<double> values, unsigned smoothness = 1);
 }
 
 #endif //SEEN_GEOM_PW_SB_CALCULUS_H
diff --git a/src/2geom/sbasis-roots.cpp b/src/2geom/sbasis-roots.cpp
index 2eb1e4cbb..5249053fa 100644
--- a/src/2geom/sbasis-roots.cpp
+++ b/src/2geom/sbasis-roots.cpp
@@ -57,7 +57,19 @@ namespace Geom{
  \returns inteval
 
 */
-Interval bounds_exact(SBasis const &a) {
+
+#ifdef USE_SBASIS_OF
+OptInterval bounds_exact(SBasisOf<double> const &a) {
+    Interval result = Interval(a.at0(), a.at1());
+    SBasisOf<double> df = derivative(a);
+    vector<double>extrema = roots(df);
+    for (unsigned i=0; i<extrema.size(); i++){
+        result.extendTo(a(extrema[i]));
+    }
+    return result;
+}
+#else
+OptInterval bounds_exact(SBasis const &a) {
     Interval result = Interval(a.at0(), a.at1());
     SBasis df = derivative(a);
     vector<double>extrema = roots(df);
@@ -66,6 +78,7 @@ Interval bounds_exact(SBasis const &a) {
     }
     return result;
 }
+#endif
 
 /** Find a small interval that bounds a
  \param a sbasis function
@@ -73,7 +86,11 @@ Interval bounds_exact(SBasis const &a) {
 
 */
 // I have no idea how this works, some clever bounding argument by jfb.
-Interval bounds_fast(const SBasis &sb, int order) {
+#ifdef USE_SBASIS_OF
+OptInterval bounds_fast(const SBasisOf<double> &sb, int order) {
+#else
+OptInterval bounds_fast(const SBasis &sb, int order) {
+#endif
     Interval res(0,0); // an empty sbasis is 0.
 
     for(int j = sb.size()-1; j>=order; j--) {
@@ -105,11 +122,15 @@ Interval bounds_fast(const SBasis &sb, int order) {
  \param sb sbasis function
  \param i domain interval
  \param order number of terms
- \returns inteval
+ \return interval
 
 */
-Interval bounds_local(const SBasis &sb, const Interval &i, int order) {
-    double t0=i.min(), t1=i.max(), lo=0., hi=0.;
+#ifdef USE_SBASIS_OF
+OptInterval bounds_local(const SBasisOf<double> &sb, const OptInterval &i, int order) {
+#else
+OptInterval bounds_local(const SBasis &sb, const OptInterval &i, int order) {
+#endif
+    double t0=i->min(), t1=i->max(), lo=0., hi=0.;
     for(int j = sb.size()-1; j>=order; j--) {
         double a=sb[j][0];
         double b=sb[j][1];
@@ -154,8 +175,13 @@ static int upper_level(vector<double> const &levels,double x,double tol=0.){
     return(upper_bound(levels.begin(),levels.end(),x-tol)-levels.begin());
 }
 
+#ifdef USE_SBASIS_OF
+static void multi_roots_internal(SBasis const &f,
+				 SBasis const &df,
+#else
 static void multi_roots_internal(SBasis const &f,
 				 SBasis const &df,
+#endif
 				 std::vector<double> const &levels,
 				 std::vector<std::vector<double> > &roots,
 				 double htol,
@@ -209,7 +235,7 @@ static void multi_roots_internal(SBasis const &f,
     int idxa=upper_level(levels,fa,vtol);
     int idxb=upper_level(levels,fb,vtol);
 
-    Interval bs = bounds_local(df,Interval(a,b));
+    Interval bs = *bounds_local(df,Interval(a,b));
 
     //first times when a level (higher or lower) can be reached from a or b.
     double ta_hi,tb_hi,ta_lo,tb_lo;
@@ -306,8 +332,8 @@ std::vector<std::vector<double> > multi_roots(SBasis const &f,
 void subdiv_sbasis(SBasis const & s,
                    std::vector<double> & roots,
                    double left, double right) {
-    Interval bs = bounds_fast(s);
-    if(bs.min() > 0 || bs.max() < 0)
+    OptInterval bs = bounds_fast(s);
+    if(!bs || bs->min() > 0 || bs->max() < 0)
         return; // no roots here
     if(s.tailError(1) < 1e-7) {
         double t = s[0][0] / (s[0][0] - s[0][1]);
diff --git a/src/2geom/sbasis.cpp b/src/2geom/sbasis.cpp
index 2619da594..0bd672c15 100644
--- a/src/2geom/sbasis.cpp
+++ b/src/2geom/sbasis.cpp
@@ -43,7 +43,7 @@ namespace Geom{
  \returns the largest possible error this truncation could give
 */
 double SBasis::tailError(unsigned tail) const {
-  Interval bs = bounds_fast(*this, tail);
+  Interval bs = *bounds_fast(*this, tail);
   return std::max(fabs(bs.min()),fabs(bs.max()));
 }
 
@@ -212,7 +212,7 @@ SBasis shift(SBasis const &a, int sh) {
 */
 SBasis shift(Linear const &a, int sh) {
     SBasis c;
-    if(sh > 0) {
+    if(sh >= 0) {
         c.insert(c.begin(), sh, Linear(0,0));
         c.push_back(a);
     }
diff --git a/src/2geom/sbasis.h b/src/2geom/sbasis.h
index 72bf422e7..f2681a722 100644
--- a/src/2geom/sbasis.h
+++ b/src/2geom/sbasis.h
@@ -43,7 +43,14 @@
 #include <2geom/utils.h>
 #include <2geom/exception.h>
 
-namespace Geom {
+
+#ifdef USE_SBASIS_OF
+
+#include "sbasis-of.h"
+
+#else
+
+namespace Geom{
 
 /*** An empty SBasis is identically 0. */
 class SBasis : public std::vector<Linear>{
@@ -135,9 +142,9 @@ private:
 inline SBasis Linear::toSBasis() const { return SBasis(*this); }
 
 //implemented in sbasis-roots.cpp
-Interval bounds_exact(SBasis const &a);
-Interval bounds_fast(SBasis const &a, int order = 0);
-Interval bounds_local(SBasis const &a, const Interval &t, int order = 0);
+OptInterval bounds_exact(SBasis const &a);
+OptInterval bounds_fast(SBasis const &a, int order = 0);
+OptInterval bounds_local(SBasis const &a, const OptInterval &t, int order = 0);
 
 /** Returns a function which reverses the domain of a.
  \param a sbasis function
@@ -325,6 +332,7 @@ std::vector<std::vector<double> > multi_roots(SBasis const &f,
                                  double b=1);
     
 }
+#endif
 
 /*
   Local Variables:
diff --git a/src/2geom/shape.cpp b/src/2geom/shape.cpp
index 814c6c68c..64abf48ff 100644
--- a/src/2geom/shape.cpp
+++ b/src/2geom/shape.cpp
@@ -349,7 +349,7 @@ void crossing_dual(unsigned &i, unsigned &j, CrossingSet const & crs) {
 
 //locate a crossing on the outside, by casting a ray through the middle of the bbox
 void outer_crossing(unsigned &ix, unsigned &jx, bool & dir, std::vector<Path> const & ps, CrossingSet const & crs) {
-    Rect bounds = ps[ix].boundsFast();
+    Rect bounds = *(ps[ix].boundsFast());
     double ry = bounds[Y].middle();
     double max_val = bounds.left(), max_t = 0;
     ix = ps.size();
diff --git a/src/2geom/svg-elliptical-arc.cpp b/src/2geom/svg-elliptical-arc.cpp
index 23645fe11..42c787eca 100644
--- a/src/2geom/svg-elliptical-arc.cpp
+++ b/src/2geom/svg-elliptical-arc.cpp
@@ -47,7 +47,7 @@ namespace Geom
 {
 
 
-Rect SVGEllipticalArc::boundsExact() const
+OptRect SVGEllipticalArc::boundsExact() const
 {
     if (isDegenerate() && is_svg_compliant())
         return chord().boundsExact();
diff --git a/src/2geom/svg-elliptical-arc.h b/src/2geom/svg-elliptical-arc.h
index 1d6d2d74d..92ec51b49 100644
--- a/src/2geom/svg-elliptical-arc.h
+++ b/src/2geom/svg-elliptical-arc.h
@@ -200,15 +200,15 @@ class SVGEllipticalArc : public Curve
         return m_svg_compliant;
     }
 
-    Rect boundsFast() const
+    virtual OptRect boundsFast() const
     {
         return boundsExact();
     }
 
-    Rect boundsExact() const;
+    virtual OptRect boundsExact() const;
 
     // TODO: native implementation of the following methods
-    Rect boundsLocal(Interval i, unsigned int deg) const
+    virtual OptRect boundsLocal(OptInterval i, unsigned int deg) const
     {
         if (isDegenerate() && is_svg_compliant())
             return chord().boundsLocal(i, deg);
diff --git a/src/2geom/svg-path-parser.cpp b/src/2geom/svg-path-parser.cpp
index 2f26870a5..071b171b3 100644
--- a/src/2geom/svg-path-parser.cpp
+++ b/src/2geom/svg-path-parser.cpp
@@ -38,6 +38,7 @@
 
 #include <2geom/point.h>
 #include <2geom/svg-path-parser.h>
+#include <2geom/angle.h>
 
 namespace Geom {
 
@@ -139,7 +140,7 @@ private:
 };
 
 
-#line 143 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.cpp"
+#line 144 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.cpp"
 static const char _svg_path_actions[] = {
 	0, 1, 0, 1, 1, 1, 2, 1, 
 	3, 1, 4, 1, 5, 1, 15, 1, 
@@ -1144,7 +1145,7 @@ static const int svg_path_first_final = 270;
 
 static const int svg_path_en_main = 1;
 
-#line 143 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
+#line 144 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
 
 
 void Parser::parse(char const *str)
@@ -1157,12 +1158,12 @@ throw(SVGPathParseError)
     _reset();
 
     
-#line 1161 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.cpp"
+#line 1162 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.cpp"
 	{
 	cs = svg_path_start;
 	}
 
-#line 1166 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.cpp"
+#line 1167 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.cpp"
 	{
 	int _klen;
 	unsigned int _trans;
@@ -1235,13 +1236,13 @@ _match:
 		switch ( *_acts++ )
 		{
 	case 0:
-#line 155 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
+#line 156 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
 	{
             start = p;
         }
 	break;
 	case 1:
-#line 159 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
+#line 160 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
 	{
             char const *end=p;
             std::string buf(start, end);
@@ -1250,55 +1251,55 @@ _match:
         }
 	break;
 	case 2:
-#line 166 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
+#line 167 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
 	{
             _push(1.0);
         }
 	break;
 	case 3:
-#line 170 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
+#line 171 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
 	{
             _push(0.0);
         }
 	break;
 	case 4:
-#line 174 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
+#line 175 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
 	{
             _absolute = true;
         }
 	break;
 	case 5:
-#line 178 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
+#line 179 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
 	{
             _absolute = false;
         }
 	break;
 	case 6:
-#line 182 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
+#line 183 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
 	{
             _moveTo(_pop_point());
         }
 	break;
 	case 7:
-#line 186 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
+#line 187 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
 	{
             _lineTo(_pop_point());
         }
 	break;
 	case 8:
-#line 190 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
+#line 191 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
 	{
             _hlineTo(Point(_pop_coord(X), _current[Y]));
         }
 	break;
 	case 9:
-#line 194 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
+#line 195 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
 	{
             _vlineTo(Point(_current[X], _pop_coord(Y)));
         }
 	break;
 	case 10:
-#line 198 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
+#line 199 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
 	{
             Point p = _pop_point();
             Point c1 = _pop_point();
@@ -1307,7 +1308,7 @@ _match:
         }
 	break;
 	case 11:
-#line 205 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
+#line 206 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
 	{
             Point p = _pop_point();
             Point c1 = _pop_point();
@@ -1315,7 +1316,7 @@ _match:
         }
 	break;
 	case 12:
-#line 211 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
+#line 212 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
 	{
             Point p = _pop_point();
             Point c = _pop_point();
@@ -1323,14 +1324,14 @@ _match:
         }
 	break;
 	case 13:
-#line 217 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
+#line 218 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
 	{
             Point p = _pop_point();
             _quadTo(_quad_tangent, p);
         }
 	break;
 	case 14:
-#line 222 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
+#line 223 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
 	{
             Point point = _pop_point();
             bool sweep = _pop_flag();
@@ -1343,16 +1344,16 @@ _match:
         }
 	break;
 	case 15:
-#line 233 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
+#line 234 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
 	{
             _closePath();
         }
 	break;
 	case 16:
-#line 369 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
+#line 370 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
 	{goto _out;}
 	break;
-#line 1356 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.cpp"
+#line 1357 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.cpp"
 		}
 	}
 
@@ -1363,7 +1364,7 @@ _again:
 	goto _resume;
 	_out: {}
 	}
-#line 379 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
+#line 380 "/home/njh/svn/lib2geom/src/2geom/svg-path-parser.rl"
 
 
     if ( cs < svg_path_first_final ) {
diff --git a/src/2geom/sweep.cpp b/src/2geom/sweep.cpp
index 227c822bd..53b953255 100644
--- a/src/2geom/sweep.cpp
+++ b/src/2geom/sweep.cpp
@@ -9,10 +9,8 @@ std::vector<std::vector<unsigned> > sweep_bounds(std::vector<Rect> rs) {
     std::vector<std::vector<unsigned> > pairs(rs.size());
 
     for(unsigned i = 0; i < rs.size(); i++) {
-        if(!rs[i].isEmpty()) {
-            events.push_back(Event(rs[i].left(), i, false));
-            events.push_back(Event(rs[i].right(), i, true));
-        }
+        events.push_back(Event(rs[i].left(), i, false));
+        events.push_back(Event(rs[i].right(), i, true));
     }
     std::sort(events.begin(), events.end());
 
@@ -48,10 +46,8 @@ std::vector<std::vector<unsigned> > sweep_bounds(std::vector<Rect> a, std::vecto
         events[n].reserve(sz*2);
         for(unsigned i = 0; i < sz; i++) {
             Rect r = n ? b[i] : a[i];
-            if(!r.isEmpty()) {
-                events[n].push_back(Event(r.left(), i, false));
-                events[n].push_back(Event(r.right(), i, true));
-            }
+            events[n].push_back(Event(r.left(), i, false));
+            events[n].push_back(Event(r.right(), i, true));
         }
         std::sort(events[n].begin(), events[n].end());
     }