update 2geom to r2049. fixes bugs!

(bzr r10837)

1 files changed, 133 insertions, 29 deletions

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.
     }