update 2geom to r2049. fixes bugs!

(bzr r10837)

38 files changed, 368 insertions, 151 deletions

diff --git a/src/2geom/CMakeLists.txt b/src/2geom/CMakeLists.txt
index 4aeb7ffac..b1c30f678 100644
--- a/src/2geom/CMakeLists.txt
+++ b/src/2geom/CMakeLists.txt
@@ -38,6 +38,7 @@ set(2geom_SRC
 	sbasis-to-bezier.cpp
 	sbasis.cpp
 	shape.cpp
+	solve-bezier.cpp
 	solve-bezier-one-d.cpp
 	solve-bezier-parametric.cpp
 	svg-elliptical-arc.cpp
diff --git a/src/2geom/basic-intersection.h b/src/2geom/basic-intersection.h
index 8e29b4617..5a813ae99 100644
--- a/src/2geom/basic-intersection.h
+++ b/src/2geom/basic-intersection.h
@@ -1,8 +1,7 @@
 /**
- * @file
- * Basic intersection routines.
- */
-/*
+ * \file
+ * \brief Basic intersection routines
+ *
  * Authors:
  *      ? <?@?.?>
  * 
diff --git a/src/2geom/bezier.h b/src/2geom/bezier.h
index 48a1dc750..1e1b88587 100644
--- a/src/2geom/bezier.h
+++ b/src/2geom/bezier.h
@@ -40,6 +40,8 @@
 #include <valarray>
 #include <boost/optional.hpp>
 #include <2geom/coord.h>
+#include <2geom/choose.h>
+#include <valarray>
 #include <2geom/math-utils.h>
 #include <2geom/d2.h>
 #include <2geom/solver.h>
@@ -111,7 +113,6 @@ inline T bernsteinValueAt(double t, T const *c_, unsigned n) {
     return (tmp + tn*t*c_[n]);
 }
 
-
 class Bezier {
 private:
     std::valarray<Coord> c_;
@@ -122,6 +123,12 @@ private:
 
     friend Bezier derivative(const Bezier & a);
 
+    friend class Bernstein;
+
+    void
+    find_bezier_roots(std::vector<double> & solutions,
+                      double l, double r) const;
+
 protected:
     Bezier(Coord const c[], unsigned ord) : c_(c, ord+1){
         //std::copy(c, c+order()+1, &c_[0]);
@@ -221,17 +228,10 @@ public:
             tmp = (tmp + tn*bc*c_[i])*u;
         }
         return (tmp + tn*t*c_[n]);
-        //return subdivideArr(t, &c_[0], NULL, NULL, order());
     }
     inline Coord operator()(double t) const { return valueAt(t); }
 
     SBasis toSBasis() const;
-//    inline SBasis toSBasis() const {
-//        SBasis sb;
-//        bezier_to_sbasis(sb, (*this));
-//        return sb;
-//        //return bezier_to_sbasis(&c_[0], order());
-//    }
 
     //Only mutator
     inline Coord &operator[](unsigned ix) { return c_[ix]; }
@@ -278,7 +278,7 @@ public:
 
     std::vector<double> roots() const {
         std::vector<double> solutions;
-        find_bernstein_roots(&const_cast<std::valarray<Coord>&>(c_)[0], order(), solutions, 0, 0.0, 1.0);
+        find_bezier_roots(solutions, 0, 1);
         return solutions;
     }
     std::vector<double> roots(Interval const ivl) const {
@@ -286,11 +286,87 @@ public:
         find_bernstein_roots(&const_cast<std::valarray<Coord>&>(c_)[0], order(), solutions, 0, ivl.min(), ivl.max());
         return solutions;
     }
+
+    Bezier forward_difference(unsigned k) {
+        Bezier fd(Order(order()-k));
+        unsigned n = fd.size();
+        
+        for(unsigned i = 0; i < n; i++) {
+            fd[i] = 0;
+            for(unsigned j = i; j < n; j++) {
+                fd[i] += (((j)&1)?-c_[j]:c_[j])*choose<double>(n, j-i);
+            }
+        }
+        return fd;
+    }
+  
+    Bezier elevate_degree() const {
+        Bezier ed(Order(order()+1));
+        unsigned n = size();
+        ed[0] = c_[0];
+        ed[n] = c_[n-1];
+        for(unsigned i = 1; i < n; i++) {
+            ed[i] = (i*c_[i-1] + (n - i)*c_[i])/(n);
+        }
+        return ed;
+    }
+
+    Bezier reduce_degree() const {
+        if(order() == 0) return *this;
+        Bezier ed(Order(order()-1));
+        unsigned n = size();
+        ed[0] = c_[0];
+        ed[n-1] = c_[n]; // ensure exact endpoints
+        unsigned middle = n/2;
+        for(unsigned i = 1; i < middle; i++) {
+            ed[i] = (n*c_[i] - i*ed[i-1])/(n-i);
+        }
+        for(unsigned i = n-1; i >= middle; i--) {
+            ed[i] = (n*c_[i] - i*ed[n-i])/(i);
+        }
+        return ed;
+    }
+
+    Bezier elevate_to_degree(unsigned newDegree) const {
+        Bezier ed = *this;
+        for(unsigned i = degree(); i < newDegree; i++) {
+            ed = ed.elevate_degree();
+        }
+        return ed;
+    }
+
+    Bezier deflate() {
+        if(order() == 0) return *this;
+        unsigned n = order();
+        Bezier b(Order(n-1));
+        for(unsigned i = 0; i < n; i++) {
+            b[i] = (n*c_[i+1])/(i+1);
+        }
+        return b;
+    }
 };
 
 
 void bezier_to_sbasis (SBasis & sb, Bezier const& bz);
 
+inline
+Bezier multiply(Bezier const& f, Bezier const& g) {
+    unsigned m = f.order();
+    unsigned n = g.order();
+    Bezier h(Bezier::Order(m+n));
+    // h_k = sum_(i+j=k) (m i)f_i (n j)g_j / (m+n k)
+    
+    for(unsigned i = 0; i <= m; i++) {
+        const double fi = choose<double>(m,i)*f[i];
+        for(unsigned j = 0; j <= n; j++) {
+            h[i+j] += fi * choose<double>(n,j)*g[j];
+        }
+    }
+    for(unsigned k = 0; k <= m+n; k++) {
+        h[k] /= choose<double>(m+n, k);
+    }
+    return h;
+}
 
 inline
 SBasis Bezier::toSBasis() const {
@@ -403,10 +479,11 @@ inline OptInterval bounds_local(Bezier const & b, OptInterval i) {
 }
 
 inline std::ostream &operator<< (std::ostream &out_file, const Bezier & b) {
+    out_file << "Bezier(";
     for(unsigned i = 0; i < b.size(); i++) {
         out_file << b[i] << ", ";
     }
-    return out_file;
+    return out_file << ")";
 }
 
 }
diff --git a/src/2geom/conic_section_clipper.h b/src/2geom/conic_section_clipper.h
index 3e4ac8429..a02cda4d3 100644
--- a/src/2geom/conic_section_clipper.h
+++ b/src/2geom/conic_section_clipper.h
@@ -1,8 +1,7 @@
 /**
- * @file
- * Conic section clipping with respect to a rectangle.
- */
-/*
+ * \file
+ * \brief Conic section clipping with respect to a rectangle
+ *
  * Authors:
  *      Marco Cecchetti <mrcekets at gmail>
  *
diff --git a/src/2geom/conic_section_clipper_cr.h b/src/2geom/conic_section_clipper_cr.h
index 687fa182d..31f5a4269 100644
--- a/src/2geom/conic_section_clipper_cr.h
+++ b/src/2geom/conic_section_clipper_cr.h
@@ -1,8 +1,7 @@
 /**
- * @file
- * Conic section clipping with respect to a rectangle.
- */
-/*
+ * \file
+ * \brief Conic section clipping with respect to a rectangle
+ *
  * Authors:
  *      Marco Cecchetti <mrcekets at gmail>
  *
diff --git a/src/2geom/conic_section_clipper_impl.cpp b/src/2geom/conic_section_clipper_impl.cpp
index 33a218a8c..e091d1e05 100644
--- a/src/2geom/conic_section_clipper_impl.cpp
+++ b/src/2geom/conic_section_clipper_impl.cpp
@@ -438,7 +438,7 @@ bool CLIPPER_CLASS::clip (std::vector<RatQuad> & arcs)
     // if the conic does not cross any line passing through a rectangle edge or
     // it is tangent to only one edge then it is an ellipse
     if (no_crossing
-            || (crossing_points.size() == 1 && single_points.size() == 0))
+            || (crossing_points.size() == 1 && single_points.empty()))
     {
         // if the ellipse centre is inside the rectangle
         // then so it is the ellipse
@@ -477,7 +477,7 @@ bool CLIPPER_CLASS::clip (std::vector<RatQuad> & arcs)
         // in case the conic section intersects any of the four lines passing
         // through the rectangle edges but it does not cross any rectangle edge
         // then the conic is all outer of the rectangle
-        if (crossing_points.size() == 0) return false;
+        if (crossing_points.empty()) return false;
         // else we need to pair crossing points, and to find an arc inner point
         // in order to generate a RatQuad object
         pairing (paired_points, inner_points, crossing_points);
diff --git a/src/2geom/conic_section_clipper_impl.h b/src/2geom/conic_section_clipper_impl.h
index ba213b8d5..37415df97 100644
--- a/src/2geom/conic_section_clipper_impl.h
+++ b/src/2geom/conic_section_clipper_impl.h
@@ -1,8 +1,7 @@
 /**
- * @file
- * Conic section clipping with respect to a rectangle.
- */
-/*
+ * \file
+ * \brief Conic section clipping with respect to a rectangle
+ *
  * Authors:
  *      Marco Cecchetti <mrcekets at gmail>
  *
@@ -116,7 +115,7 @@ class CLIPPER_CLASS
 
     bool found_any_isolated_point() const
     {
-        return (single_points.size() != 0);
+        return ( !single_points.empty() );
     }
 
     const std::vector<Point> & isolated_points() const
diff --git a/src/2geom/conicsec.h b/src/2geom/conicsec.h
index d9c5e7bc5..e48d149ec 100644
--- a/src/2geom/conicsec.h
+++ b/src/2geom/conicsec.h
@@ -1,8 +1,7 @@
 /**
- * @file
- * Conic Section.
- */
-/*
+ * \file
+ * \brief Conic Section
+ *
  * Authors:
  *      Nathan Hurst <njh@njhurst.com>
  *
@@ -485,7 +484,7 @@ public:
     Point nearestPoint (const Point P) const
     {
         std::vector<Point> points = allNearestPoints (P);
-        if (points.size() != 0)
+        if ( !points.empty() )
         {
             return points.front();
         }
diff --git a/src/2geom/convex-cover.cpp b/src/2geom/convex-cover.cpp
index d50accadf..550779f65 100644
--- a/src/2geom/convex-cover.cpp
+++ b/src/2geom/convex-cover.cpp
@@ -139,13 +139,22 @@ ConvexHull::angle_sort() {
 
 void
 ConvexHull::graham_scan() {
-    if (boundary.size() < 4) return;
+    // prune out equal points.  points are sorted, so equals are adjacent
+    std::vector<Point>::iterator e = 
+        std::unique(boundary.begin(), boundary.end());
+    boundary.resize(e - boundary.begin());
+    for(unsigned int i = 2; i < boundary.size(); i++) {
+        
+    }
+    if (boundary.size() < 4) {
+        return;
+    }
     unsigned stac = 2;
     for(unsigned int i = 2; i < boundary.size(); i++) {
         double o = SignedTriangleArea(boundary[stac-2],
                                       boundary[stac-1],
                                       boundary[i]);
-        if(fabs(o) < 1e-8) { // colinear - dangerous...
+        if(o == 0) { // colinear - dangerous...
             stac--;
         } else if(o < 0) { // anticlockwise
         } else { // remove concavity
@@ -161,13 +170,93 @@ ConvexHull::graham_scan() {
     boundary.resize(stac);
 }
 
+// following code is from marco.
+
+/*
+ * return true in case the oriented polyline p0, p1, p2 is a right turn
+ */
+inline
+bool is_a_right_turn (Point const& p0, Point const& p1, Point const& p2)
+{
+    if (p1 == p2) return false;
+    Point q1 = p1 - p0;
+    Point q2 = p2 - p0;
+    if (q1 == -q2) return false;
+    return (cross (q1, q2) < 0);
+}
+
+/*
+ * return true if p < q wrt the lexicographyc order induced by the coordinates
+ */
+struct lex_less
+{
+    bool operator() (Point const& p, Point const& q)
+    {
+      return ((p[Y] < q[Y]) || (p[Y] == q[Y] && p[X] < q[X]));
+    }
+};
+
+/*
+ * return true if p > q wrt the lexicographyc order induced by the coordinates
+ */
+struct lex_greater
+{
+    bool operator() (Point const& p, Point const& q)
+    {
+        return ((p[Y] > q[Y]) || (p[Y] == q[Y] && p[X] > q[X]));
+    }
+};
+
+/*
+ * Compute the convex hull of a set of points.
+ * The implementation is based on the Andrew's scan algorithm
+ * note: in the Bezier clipping for collinear normals it seems
+ * to be more stable wrt the Graham's scan algorithm and in general
+ * a bit quikier
+ */
+void ConvexHull::andrew_scan ()
+{
+    vector<Point> & P = boundary;
+    size_t n = P.size();
+    if (n < 2)  return;
+    std::sort(P.begin(), P.end(), lex_less());
+    if (n < 4) return;
+    // upper hull
+    size_t u = 2;
+    for (size_t i = 2; i < n; ++i)
+    {
+        while (u > 1 && !is_a_right_turn(P[u-2], P[u-1], P[i]))
+        {
+            --u;
+        }
+        std::swap(P[u], P[i]);
+        ++u;
+    }
+    std::sort(P.begin() + u, P.end(), lex_greater());
+    std::rotate(P.begin(), P.begin() + 1, P.end());
+    // lower hull
+    size_t l = u;
+    size_t k = u - 1;
+    for (size_t i = l; i < n; ++i)
+    {
+        while (l > k && !is_a_right_turn(P[l-2], P[l-1], P[i]))
+        {
+            --l;
+        }
+        std::swap(P[l], P[i]);
+        ++l;
+    }
+    P.resize(l);
+}
+
 void
 ConvexHull::graham() {
-    if(is_degenerate()) // nothing to do
-        return;
-    find_pivot();
-    angle_sort();
-    graham_scan();
+    /*if(is_degenerate()) // nothing to do
+      return;*/
+    //find_pivot();
+    //angle_sort();
+    //graham_scan();
+    andrew_scan();
 }
 
 //Mathematically incorrect mod, but more useful.
@@ -181,17 +270,17 @@ int mod(int i, int l) {
  * Tests if a point is left (outside) of a particular segment, n. */
 bool
 ConvexHull::is_left(Point p, int n) {
-    return SignedTriangleArea((*this)[n], (*this)[n+1], p) >= 0;
+    return SignedTriangleArea((*this)[n], (*this)[n+1], p) > 0;
 }
 
 /*** ConvexHull::strict_left
  * Tests if a point is left (outside) of a particular segment, n. */
 bool
 ConvexHull::is_strict_left(Point p, int n) {
-    return SignedTriangleArea((*this)[n], (*this)[n+1], p) > 0;
+    return SignedTriangleArea((*this)[n], (*this)[n+1], p) >= 0;
 }
 
-/*** ConvexHull::find_positive
+/*** ConvexHull::find_left
  * May return any number n where the segment n -> n + 1 (possibly looped around) in the hull such
  * that the point is on the wrong side to be within the hull.  Returns -1 if it is within the hull.*/
 int
@@ -223,6 +312,7 @@ ConvexHull::find_strict_left(Point p) {
  * sure that each triangle made from an edge and the point has positive area. */
 bool
 ConvexHull::contains_point(Point p) {
+    if(size() == 0) return false;
     return find_left(p) == -1;
 }
 
@@ -231,6 +321,7 @@ ConvexHull::contains_point(Point p) {
  * sure that each triangle made from an edge and the point has positive area. */
 bool
 ConvexHull::strict_contains_point(Point p) {
+    if(size() == 0) return false;
     return find_strict_left(p) == -1;
 }
 
@@ -244,15 +335,16 @@ ConvexHull::merge(Point p) {
     int l = boundary.size();
 
     if(l < 2) {
-        boundary.push_back(p);
+        if(boundary.empty() || boundary[0] != p)
+            boundary.push_back(p);
         return;
     }
 
     bool pushed = false;
 
-    bool pre = is_strict_left(p, -1);
+    bool pre = is_left(p, -1);
     for(int i = 0; i < l; i++) {
-        bool cur = is_strict_left(p, i);
+        bool cur = is_left(p, i);
         if(pre) {
             if(cur) {
                 if(!pushed) {
@@ -277,6 +369,7 @@ ConvexHull::merge(Point p) {
 
 /*** ConvexHull::is_clockwise
  * We require that successive pairs of edges always turn right.
+ * We return false on collinear points
  * proposed algorithm: walk successive edges and require triangle area is positive.
  */
 bool
@@ -302,6 +395,7 @@ ConvexHull::is_clockwise() const {
  */
 bool
 ConvexHull::top_point_first() const {
+    if(size() <= 1) return true;
     std::vector<Point>::const_iterator pivot = boundary.begin();
     for(std::vector<Point>::const_iterator it(boundary.begin()+1),
             e(boundary.end());
@@ -316,19 +410,9 @@ ConvexHull::top_point_first() const {
 }
 //OPT: since the Y values are orderly there should be something like a binary search to do this.
 
-/*** ConvexHull::no_colinear_points
- * We require that no three vertices are colinear.
-proposed algorithm:  We must be very careful about rounding here.
-*/
-bool
-ConvexHull::no_colinear_points() const {
-    // XXX: implement me!
-    THROW_NOTIMPLEMENTED();
-}
-
 bool
 ConvexHull::meets_invariants() const {
-    return is_clockwise() && top_point_first() && no_colinear_points();
+    return is_clockwise() && top_point_first();
 }
 
 /*** ConvexHull::is_degenerate
@@ -394,7 +478,7 @@ std::vector<pair<int, int> > bridges(ConvexHull a, ConvexHull b) {
             bi++;
             bp_angle += angle_between(b[bi] - b[bi-1], b[bi+1] - b[bi]);
             L[1] = b[bi];
-            std::cout << "parallel\n";
+            //std::cout << "parallel\n";
         } else if(ap_angle < bp_angle) {
             ai++;
             ap_angle += angle_between(a[ai] - a[ai-1], a[ai+1] - a[ai]);
@@ -475,7 +559,7 @@ T idx_to_pair(pair<T, T> p, int idx) {
 ConvexHull merge(ConvexHull a, ConvexHull b) {
     ConvexHull ret;
 
-    std::cout << "---\n";
+    //std::cout << "---\n";
     std::vector<pair<int, int> > bpair = bridges(a, b);
     
     // Given our list of bridges {(pb1, qb1), ..., (pbk, qbk)}
@@ -490,10 +574,10 @@ ConvexHull merge(ConvexHull a, ConvexHull b) {
     
     for(unsigned k = 0; k < bpair.size(); k++) {
         unsigned limit = idx_to_pair(bpair[k], state);
-        std::cout << bpair[k].first << " , " << bpair[k].second << "; "
+        /*std::cout << bpair[k].first << " , " << bpair[k].second << "; "
                   << idx << ", " << limit << ", s: "
                   << state
-                  << " \n";
+                  << " \n";*/
         while(idx <= limit) {
             ret.boundary.push_back(chs[state][idx++]);
         }
@@ -551,6 +635,26 @@ ConvexHull graham_merge(ConvexHull a, ConvexHull b) {
 
     return result;
 }
+
+ConvexHull andrew_merge(ConvexHull a, ConvexHull b) {
+    ConvexHull result;
+
+    // we can avoid the find pivot step because of top_point_first
+    if(b.boundary[0] <= a.boundary[0])
+        std::swap(a, b);
+
+    result.boundary = a.boundary;
+    result.boundary.insert(result.boundary.end(),
+                           b.boundary.begin(), b.boundary.end());
+
+/** if we modified graham scan to work top to bottom as proposed in lect754.pdf we could replace the
+ angle sort with a simple merge sort type algorithm. furthermore, we could do the graham scan
+ online, avoiding a bunch of memory copies.  That would probably be linear. -- njh*/
+    result.andrew_scan();
+
+    return result;
+}
+
 //TODO: reinstate
 /*ConvexCover::ConvexCover(Path const &sp) : path(&sp) {
     cc.reserve(sp.size());
@@ -570,7 +674,7 @@ double ConvexHull::centroid_and_area(Geom::Point& centroid) const {
     Geom::Point centroid_tmp(0,0);
     double atmp = 0;
     for (unsigned i = n-1, j = 0; j < n; i = j, j++) {
-        const double ai = -cross(boundary[j], boundary[i]);
+        const double ai = cross(boundary[j], boundary[i]);
         atmp += ai;
         centroid_tmp += (boundary[j] + boundary[i])*ai; // first moment.
     }
diff --git a/src/2geom/convex-cover.h b/src/2geom/convex-cover.h
index 4c1b59d17..a0443e9dc 100644
--- a/src/2geom/convex-cover.h
+++ b/src/2geom/convex-cover.h
@@ -1,8 +1,7 @@
 /**
- * @file
- * Dynamic convex hull structure.
- */
-/*
+ * \file
+ * \brief Dynamic convex hull structure
+ *
  * Copyright 2006 Nathan Hurst <njh@mail.csse.monash.edu.au>
  * Copyright 2006 Michael G. Sloan <mgsloan@gmail.com>
  *
@@ -59,6 +58,7 @@ public: // XXX: should be private :)
     void find_pivot();
     void angle_sort();
     void graham_scan();
+    void andrew_scan();
     void graham();
 public:
     std::vector<Point> boundary;
@@ -100,7 +100,6 @@ public:
     /** Is the convex hull clockwise?  We use the definition of clockwise from point.h
     **/
     bool is_clockwise() const;
-    bool no_colinear_points() const;
     bool top_point_first() const;
     bool meets_invariants() const;
 
@@ -131,6 +130,16 @@ public:
     double narrowest_diameter(Point &a, Point &b, Point &c);
 
 };
+/** @brief Output operator for points.
+ * Prints out all the coordinates. */
+inline std::ostream &operator<< (std::ostream &out_file, const Geom::ConvexHull &in_cvx) {
+    out_file << "ConvexHull(";
+    for(unsigned i = 0; i < in_cvx.size(); i++) {
+        out_file << in_cvx.boundary[i] << ", ";
+    }
+    out_file << ")";
+    return out_file;
+}
 
 // do two convex hulls intersect?
 bool intersectp(ConvexHull a, ConvexHull b);
@@ -147,6 +156,9 @@ ConvexHull merge(ConvexHull a, ConvexHull b);
 // naive approach
 ConvexHull graham_merge(ConvexHull a, ConvexHull b);
 
+// naive approach
+ConvexHull andrew_merge(ConvexHull a, ConvexHull b);
+
 unsigned find_bottom_right(ConvexHull const &a);
 
 /*** Arbitrary transform operator.
diff --git a/src/2geom/crossing.cpp b/src/2geom/crossing.cpp
index 13affa8e9..379eb4b00 100644
--- a/src/2geom/crossing.cpp
+++ b/src/2geom/crossing.cpp
@@ -154,7 +154,7 @@ Crossings reverse_tb(Crossings const &cr, unsigned split, std::vector<double> ma
     Crossings ret;
     for(Crossings::const_iterator i = cr.begin(); i != cr.end(); ++i) {
         double mx = max[i->b - split];
-        std::cout << i->b << "\n";
+        //std::cout << i->b << "\n";
         ret.push_back(Crossing(i->ta, i->tb > mx+0.01 ? (1 - (i->tb - mx) + mx) : mx - i->tb,
                                !i->dir));
     }
diff --git a/src/2geom/ellipse.h b/src/2geom/ellipse.h
index 2d6ba399f..297254366 100644
--- a/src/2geom/ellipse.h
+++ b/src/2geom/ellipse.h
@@ -1,8 +1,7 @@
 /**
- * @file
- * Ellipse Curve.
- */
-/*
+ * \file
+ * \brief  Ellipse Curve
+ *
  * Authors:
  *      Marco Cecchetti <mrcekets at gmail.com>
  *
diff --git a/src/2geom/generic-rect.h b/src/2geom/generic-rect.h
index 719b37385..f0863f439 100644
--- a/src/2geom/generic-rect.h
+++ b/src/2geom/generic-rect.h
@@ -214,6 +214,22 @@ public:
 
     /// @name Modify the rectangle.
     /// @{
+    /** @brief Set the minimum X coordinate of the rectangle. */
+    void setLeft(C val) {
+        f[X].setMin(val);
+    }
+    /** @brief Set the maximum X coordinate of the rectangle. */
+    void setRight(C val) {
+        f[X].setMax(val);
+    }
+    /** @brief Set the minimum Y coordinate of the rectangle. */
+    void setTop(C val) {
+        f[Y].setMin(val);
+    }
+    /** @brief Set the maximum Y coordinate of the rectangle. */
+    void setBottom(C val) {
+        f[Y].setMax(val);
+    }
     /** @brief Set the upper left point of the rectangle. */
     void setMin(CPoint const &p) {
         f[X].setMin(p[X]);
diff --git a/src/2geom/geom.h b/src/2geom/geom.h
index b9d910e2a..5aeded23d 100644
--- a/src/2geom/geom.h
+++ b/src/2geom/geom.h
@@ -1,8 +1,7 @@
 /**
- * @file
- * Various geometrical calculations.
- */
-/*
+ *  \file
+ *  \brief Various geometrical calculations
+ *
  *  Authors:
  *   Nathan Hurst <njh@mail.csse.monash.edu.au>
  *
diff --git a/src/2geom/line.cpp b/src/2geom/line.cpp
index 5ef8dfddc..471ce256c 100644
--- a/src/2geom/line.cpp
+++ b/src/2geom/line.cpp
@@ -464,7 +464,7 @@ boost::optional<LineSegment> clip (Line const& l, Rect const& r)
             return opt_linesegment(ls);
         }
     }
-    if (points.size() != 0)
+    if ( !points.empty() )
     {
         result.setInitial (points[0]);
         result.setFinal (points[0]);
diff --git a/src/2geom/linear.h b/src/2geom/linear.h
index 6302d810e..df6dd9904 100644
--- a/src/2geom/linear.h
+++ b/src/2geom/linear.h
@@ -1,8 +1,7 @@
 /**
- * @file
- * Linear fragment function class.
- */
-/*
+ * \file
+ * \brief  Linear fragment function class
+ *
  *  Authors:
  *   Nathan Hurst <njh@mail.csse.monash.edu.au>
  *   Michael Sloan <mgsloan@gmail.com>
diff --git a/src/2geom/nearest-point.cpp b/src/2geom/nearest-point.cpp
index 71beab9cb..68a863434 100644
--- a/src/2geom/nearest-point.cpp
+++ b/src/2geom/nearest-point.cpp
@@ -61,6 +61,7 @@ double nearest_point( Point const& p,
     }
     if (c.isConstant()) return from;
     SBasis dd = dot(c - p, dc);
+    //std::cout << dd << std::endl;
     std::vector<double> zeros = Geom::roots(dd);
 
     double closest = from;
diff --git a/src/2geom/nearest-point.h b/src/2geom/nearest-point.h
index 484c47afc..19485242c 100644
--- a/src/2geom/nearest-point.h
+++ b/src/2geom/nearest-point.h
@@ -1,8 +1,7 @@
 /**
- * @file
- * nearest point routines for D2<SBasis> and Piecewise<D2<SBasis>>.
- */
-/*
+ * \file
+ * \brief  nearest point routines for D2<SBasis> and Piecewise<D2<SBasis>>
+ *
  * Authors:
  *
  * 		Marco Cecchetti <mrcekets at gmail.com>
diff --git a/src/2geom/path-intersection.cpp b/src/2geom/path-intersection.cpp
index c38776304..c680f3a31 100644
--- a/src/2geom/path-intersection.cpp
+++ b/src/2geom/path-intersection.cpp
@@ -460,8 +460,8 @@ void mono_pair(Path const &A, double Al, double Ah,
 /** This returns the times when the x or y derivative is 0 in the curve. */
 std::vector<double> curve_mono_splits(Curve const &d) {
     Curve* deriv = d.derivative();
-    std::vector<double> rs = d.roots(0, X);
-    append(rs, d.roots(0, Y));
+    std::vector<double> rs = deriv->roots(0, X);
+    append(rs, deriv->roots(0, Y));
     delete deriv;
     std::sort(rs.begin(), rs.end());
     return rs;
diff --git a/src/2geom/path-intersection.h b/src/2geom/path-intersection.h
index 2470e44fb..512c31167 100644
--- a/src/2geom/path-intersection.h
+++ b/src/2geom/path-intersection.h
@@ -57,7 +57,7 @@ Crossings curve_sweep(Path const &a, Path const &b) {
     std::vector<Rect> bounds_a = bounds(a), bounds_b = bounds(b);
     std::vector<std::vector<unsigned> > ixs = sweep_bounds(bounds_a, bounds_b);
     for(unsigned i = 0; i < a.size(); i++) {
-        for(std::vector<unsigned>::iterator jp = ixs[i].begin(); jp != ixs[i].end(); jp++) {
+        for(std::vector<unsigned>::iterator jp = ixs[i].begin(); jp != ixs[i].end(); ++jp) {
             Crossings cc = t.crossings(a[i], b[*jp]);
             offset_crossings(cc, i, *jp);
             ret.insert(ret.end(), cc.begin(), cc.end());
diff --git a/src/2geom/path.cpp b/src/2geom/path.cpp
index 58d6b9b5e..266694390 100644
--- a/src/2geom/path.cpp
+++ b/src/2geom/path.cpp
@@ -255,10 +255,9 @@ 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(Geom::Point(0,0),Geom::Point(0,0));
 	for ( unsigned int i = si + 1; i < ei; ++i )
 	{
-		bb = (_path[i].boundsFast());
+            Rect bb = (_path[i].boundsFast());
 		dsq = distanceSq(_point, bb);
 		if ( mindistsq <= dsq ) continue;
 		t = _path[i].nearestPoint(_point);
@@ -270,7 +269,7 @@ double Path::nearestPoint(Point const &_point, double from, double to, double *d
 			mindistsq = dsq;
 		}
 	}
-	bb = (_path[ei].boundsFast());
+	Rect bb = (_path[ei].boundsFast());
 	dsq = distanceSq(_point, bb);
 	if ( mindistsq > dsq )
 	{
diff --git a/src/2geom/piecewise.h b/src/2geom/piecewise.h
index 310de23b1..837f33ea7 100644
--- a/src/2geom/piecewise.h
+++ b/src/2geom/piecewise.h
@@ -1,4 +1,7 @@
-/*
+/**
+ * \file
+ * \brief Piecewise function class
+ *
  * Copyright 2007 Michael Sloan <mgsloan@gmail.com>
  *
  * This library is free software; you can redistribute it and/or
diff --git a/src/2geom/quadtree.cpp b/src/2geom/quadtree.cpp
index 08e6dd7e2..f30475415 100644
--- a/src/2geom/quadtree.cpp
+++ b/src/2geom/quadtree.cpp
@@ -131,7 +131,7 @@ void QuadTree::insert(double x0, double y0, double x1, double y1, int shape) {
     // loop until a quad would break the box.
 
     // empty root => empty QuadTree. Create initial bounding box (0,0), (1,1)
-    if(root == 0) {
+    if(root == NULL) {
         root = new Quad;
             
         bx0 = 0;
@@ -179,7 +179,7 @@ void QuadTree::insert(double x0, double y0, double x1, double y1, int shape) {
         by1 = byy1;
     }
 
-    while(q) {
+    while(true) {
         // Find the center of the temp bounding box
         double cx = (bxx0 + bxx1)/2;
         double cy = (byy0 + byy1)/2;
@@ -225,7 +225,7 @@ void QuadTree::insert(double x0, double y0, double x1, double y1, int shape) {
         */
         assert(i < 4);
         Quad *qq = q->children[i];
-        if(qq == 0) {
+        if(qq == NULL) {
             qq = new Quad;
             q->children[i] = qq;
         }
@@ -233,6 +233,8 @@ void QuadTree::insert(double x0, double y0, double x1, double y1, int shape) {
     }
     q->data.push_back(shape);
 }
+
+
 void QuadTree::erase(Quad *q, int shape) {
     for(Quad::iterator i = q->data.begin();  i != q->data.end(); i++) {
         if(*i == shape) {
diff --git a/src/2geom/sbasis-geometric.cpp b/src/2geom/sbasis-geometric.cpp
index 1c180e143..017b377c9 100644
--- a/src/2geom/sbasis-geometric.cpp
+++ b/src/2geom/sbasis-geometric.cpp
@@ -4,8 +4,7 @@
 //#include <2geom/solver.h>
 #include <2geom/sbasis-geometric.h>
 
-/*
- * Geometric operators on D2<SBasis> (1D->2D).
+/** Geometric operators on D2<SBasis> (1D->2D).
  * Copyright 2007 JF Barraud
  * Copyright 2007 N Hurst
  *
@@ -746,18 +745,26 @@ Geom::cubics_with_prescribed_curvature(Point const &M0,   Point const &M1,
 }
 
 
+namespace Geom {
 /**
 * \brief returns all the parameter values of A whose tangent passes through P.
 * \relates D2
 */
-std::vector<double> Geom::find_tangents(Point P, D2<SBasis> const &A) {
+std::vector<double> find_tangents(Point P, D2<SBasis> const &A) {
     SBasis crs (cross(A - P, derivative(A)));
-    crs = shift(crs*Linear(-1, 0)*Linear(-1, 0), -2); // We know that there is a double root at t=0 so we divide out t^2
-// JFB points out that this is equivalent to (t-1)^2 followed by a divide by s^2 (shift)
     return roots(crs);
 }
 
+/**
+* \brief returns all the parameter values of A whose normal passes through P.
+* \relates D2
+*/
+std::vector<double> find_normals(Point P, D2<SBasis> const &A) {
+    SBasis crs (dot(A - P, derivative(A)));
+    return roots(crs);
+}
 
+}
 //}; // namespace
 
 
diff --git a/src/2geom/sbasis-geometric.h b/src/2geom/sbasis-geometric.h
index 841d75a28..00b64b430 100644
--- a/src/2geom/sbasis-geometric.h
+++ b/src/2geom/sbasis-geometric.h
@@ -1,14 +1,13 @@
-#ifndef SEEN_SBASIS_GEOMETRIC_
-#define SEEN_SBASIS_GEOMETRIC_
+#ifndef _SBASIS_GEOMETRIC
+#define _SBASIS_GEOMETRIC
 #include <2geom/d2.h>
 #include <2geom/piecewise.h>
 #include <vector>
 
 /**
- * @file
- * two-dimensional geometric operators.
- */
-/*
+ * \file
+ * \brief two-dimensional geometric operators.  
+ *
  * Copyright 2007, JFBarraud
  * Copyright 2007, njh
  * 
@@ -102,10 +101,11 @@ cubics_with_prescribed_curvature(Point const &M0,   Point const &M1,
 
 
 std::vector<double> find_tangents(Point P, D2<SBasis> const &A);
+std::vector<double> find_normals(Point P, D2<SBasis> const &A);
 
 };
 
-#endif // SEEN_SBASIS_GEOMETRIC_
+#endif
 
 /*
   Local Variables:
diff --git a/src/2geom/sbasis-math.h b/src/2geom/sbasis-math.h
index c3f7518eb..e6d40a3de 100644
--- a/src/2geom/sbasis-math.h
+++ b/src/2geom/sbasis-math.h
@@ -1,8 +1,7 @@
 /**
- * @file
- * some std functions to work with (pw)s-basis.
- */
-/*
+ * \file
+ * \brief some std functions to work with (pw)s-basis
+ *
  *  Authors:
  *   Jean-Francois Barraud
  *
diff --git a/src/2geom/sbasis-poly.h b/src/2geom/sbasis-poly.h
index 70abfeea1..e0bef9333 100644
--- a/src/2geom/sbasis-poly.h
+++ b/src/2geom/sbasis-poly.h
@@ -5,10 +5,9 @@
 #include <2geom/sbasis.h>
 
 /**
- * @file
- * Conversion between SBasis and Poly.  Not recommended for general use due to instability.
- */
-/*
+ * \file
+ * \brief Conversion between SBasis and Poly.  Not recommended for general use due to instability.
+ *
  * Authors:
  *      ? <?@?.?>
  * 
diff --git a/src/2geom/sbasis-roots.cpp b/src/2geom/sbasis-roots.cpp
index 4602fced9..1b870d88f 100644
--- a/src/2geom/sbasis-roots.cpp
+++ b/src/2geom/sbasis-roots.cpp
@@ -1,5 +1,4 @@
-/*
- * root finding for sbasis functions.
+/** root finding for sbasis functions.
  * Copyright 2006 N Hurst
  * Copyright 2007 JF Barraud
  *
diff --git a/src/2geom/sbasis-to-bezier.cpp b/src/2geom/sbasis-to-bezier.cpp
index aabafabea..8e47580ce 100644
--- a/src/2geom/sbasis-to-bezier.cpp
+++ b/src/2geom/sbasis-to-bezier.cpp
@@ -390,20 +390,24 @@ path_from_piecewise(Geom::Piecewise<Geom::D2<Geom::SBasis> > const &B, double to
     Geom::Point start = B[0].at0();
     pb.moveTo(start);
     for(unsigned i = 0; ; i++) {
-        if(i+1 == B.size() || !are_near(B[i+1].at0(), B[i].at1(), tol)) {
+        if ( (i+1 == B.size()) 
+             || !are_near(B[i+1].at0(), B[i].at1(), tol) )
+        {
             //start of a new path
-            if(are_near(start, B[i].at1()) && sbasis_size(B[i]) <= 1) {
+            if (are_near(start, B[i].at1()) && sbasis_size(B[i]) <= 1) {
                 pb.closePath();
                 //last line seg already there (because of .closePath())
                 goto no_add;
             }
             build_from_sbasis(pb, B[i], tol, only_cubicbeziers);
-            if(are_near(start, B[i].at1())) {
+            if (are_near(start, B[i].at1())) {
                 //it's closed, the last closing segment was not a straight line so it needed to be added, but still make it closed here with degenerate straight line.
                 pb.closePath();
             }
           no_add:
-            if(i+1 >= B.size()) break;
+            if (i+1 >= B.size()) {
+                break;
+            }
             start = B[i+1].at0();
             pb.moveTo(start);
         } else {
diff --git a/src/2geom/sbasis.h b/src/2geom/sbasis.h
index f3598ccc8..dd2df55b9 100644
--- a/src/2geom/sbasis.h
+++ b/src/2geom/sbasis.h
@@ -1,8 +1,7 @@
 /**
- * @file
- * Defines S-power basis function class.
- */
-/*
+ * \file
+ * \brief Defines S-power basis function class
+ *
  *  Authors:
  *   Nathan Hurst <njh@mail.csse.monash.edu.au>
  *   Michael Sloan <mgsloan@gmail.com>
@@ -110,6 +109,9 @@ public:
     SBasis(SBasis const & a) :
         d(a.d)
     {}
+    SBasis(std::vector<Linear> const & ls) :
+        d(ls)
+    {}
     SBasis(Linear const & bo) {
         push_back(bo);
     }
diff --git a/src/2geom/shape.cpp b/src/2geom/shape.cpp
index c389b9c6e..e9f5e55dc 100644
--- a/src/2geom/shape.cpp
+++ b/src/2geom/shape.cpp
@@ -1,5 +1,5 @@
-/*
- * Shapes are special paths on which boolops can be performed.
+/**
+ * \brief  Shapes are special paths on which boolops can be performed
  *
  * Authors:
  *      Michael G. Sloan <mgsloan@gmail.com>
diff --git a/src/2geom/shape.h b/src/2geom/shape.h
index 9f7ead4aa..0a7ee9709 100644
--- a/src/2geom/shape.h
+++ b/src/2geom/shape.h
@@ -1,4 +1,6 @@
-/*
+/**
+ * \brief  Shapes are special paths on which boolops can be performed
+ *
  * Authors:
  *      Michael G. Sloan <mgsloan@gmail.com>
  *      Nathan Hurst <njh@mail.csse.monash.edu.au>
@@ -61,9 +63,6 @@ enum {
   BOOLOP_UNION        = BOOLOP_JUST_A | BOOLOP_JUST_B | BOOLOP_BOTH
 };
 
-/**
- * Shapes are special paths on which boolops can be performed.
- */
 class Shape {
     Regions content;
     mutable bool fill;
diff --git a/src/2geom/solve-bezier-one-d.cpp b/src/2geom/solve-bezier-one-d.cpp
index c74e3243c..f14f701d0 100644
--- a/src/2geom/solve-bezier-one-d.cpp
+++ b/src/2geom/solve-bezier-one-d.cpp
@@ -71,7 +71,7 @@ find_bernstein_roots(double const *w, /* The control points  */
 
 void
 find_bernstein_roots(std::vector<double> &solutions, /* RETURN candidate t-values */
-                     Geom::Bezier const& bz, /* The control points  */
+                     Geom::Bezier const &bz, /* The control points  */
                      double left_t, double right_t)
 {
     Bernsteins B(bz.degree(), solutions);
diff --git a/src/2geom/svg-path-parser.cpp b/src/2geom/svg-path-parser.cpp
index 24f9c81ee..60a734ece 100644
--- a/src/2geom/svg-path-parser.cpp
+++ b/src/2geom/svg-path-parser.cpp
@@ -1,6 +1,7 @@
 #line 1 "/opt/shared/work/programming/eclipse/eclipse_3.4/lib2geom/src/2geom/svg-path-parser.rl"
-/*
- * parse SVG path specifications.
+/**
+ * \file
+ * \brief parse SVG path specifications
  *
  * Copyright 2007 MenTaLguY <mental@rydia.net>
  * Copyright 2007 Aaron Spike <aaron@ekips.org>
@@ -45,7 +46,7 @@ namespace {
 
 class Parser {
 public:
-    Parser(SVGPathSink &sink) : _sink(sink) {}
+    Parser(SVGPathSink &sink) : _absolute(false), _sink(sink) {}
 
     void parse(char const *str) throw(SVGPathParseError);
 
diff --git a/src/2geom/svg-path-parser.h b/src/2geom/svg-path-parser.h
index 1aab8bb36..365287d8c 100644
--- a/src/2geom/svg-path-parser.h
+++ b/src/2geom/svg-path-parser.h
@@ -1,8 +1,7 @@
 /**
- * @file
- * parse SVG path specifications.
- */
-/*
+ * \file
+ * \brief  parse SVG path specifications
+ *
  * Copyright 2007 MenTaLguY <mental@rydia.net>
  * Copyright 2007 Aaron Spike <aaron@ekips.org>
  *
@@ -57,13 +56,18 @@ inline std::vector<Path> parse_svg_path(char const *str) throw(SVGPathParseError
     return subpaths;
 }
 
+inline std::vector<Path> read_svgd_f(FILE * fi) throw(SVGPathParseError) {
+    char input[1024 * 10];
+    fgets(input, 1024 * 10, fi);
+    return parse_svg_path(input);
+}
+
 inline std::vector<Path> read_svgd(char const * name) throw(SVGPathParseError) {
     FILE* fi = fopen(name, "r");
     if(fi == NULL) throw(std::runtime_error("Error opening file"));
-    char input[1024 * 10];
-    fgets(input, 1024 * 10, fi);
+    std::vector<Path> out = read_svgd_f(fi);
     fclose(fi);
-    return parse_svg_path(input);
+    return out;
 }
 
 }
diff --git a/src/2geom/svg-path.h b/src/2geom/svg-path.h
index 89dcb4759..89192fb72 100644
--- a/src/2geom/svg-path.h
+++ b/src/2geom/svg-path.h
@@ -1,8 +1,7 @@
 /**
- * @file
- * callback interface for SVG path data.
- */
-/*
+ * \file
+ * \brief  callback interface for SVG path data
+ *
  * Copyright 2007 MenTaLguY <mental@rydia.net>
  *
  * This library is free software; you can redistribute it and/or
diff --git a/src/2geom/toposweep.h b/src/2geom/toposweep.h
index b8c1bdcf4..428115dd3 100644
--- a/src/2geom/toposweep.h
+++ b/src/2geom/toposweep.h
@@ -2,8 +2,7 @@
 /**
  * \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>
diff --git a/src/2geom/utils.cpp b/src/2geom/utils.cpp
index 4c9629deb..a40b7253d 100644
--- a/src/2geom/utils.cpp
+++ b/src/2geom/utils.cpp
@@ -1,5 +1,4 @@
-/*
- * Various utility functions.
+/** Various utility functions.
  *
  * Copyright 2008 Marco Cecchetti <mrcekets at gmail.com>
  * Copyright 2007 Johan Engelen <goejendaagh@zonnet.nl>