X-Git-Url: https://git.tokkee.org/?a=blobdiff_plain;f=src%2F2geom%2Fconvex-cover.cpp;h=e8ea2280d9d77ed808ef50e2e5f9a5efe8346d19;hb=698cd69991e23d81f6dccf30b2014d8cf0981d11;hp=1c704e0ede2839733d38e0d0fef07418ca0b8c77;hpb=981b809bc6ed10a21e89444d9447e5475801874f;p=inkscape.git diff --git a/src/2geom/convex-cover.cpp b/src/2geom/convex-cover.cpp index 1c704e0ed..e8ea2280d 100644 --- a/src/2geom/convex-cover.cpp +++ b/src/2geom/convex-cover.cpp @@ -29,9 +29,11 @@ * */ -#include "convex-cover.h" +#include <2geom/convex-cover.h> +#include <2geom/exception.h> #include #include + /** Todo: + modify graham scan to work top to bottom, rather than around angles + intersection @@ -61,12 +63,16 @@ class angle_cmp{ public: Point o; angle_cmp(Point o) : o(o) {} - + +#if 0 bool operator()(Point a, Point b) { + // not remove this check or std::sort could crash + if (a == b) return false; Point da = a - o; Point db = b - o; - + if (da == -db) return false; + #if 1 double aa = da[0]; double ab = db[0]; @@ -91,16 +97,33 @@ public: return L2sq(da) < L2sq(db); return false; } +#else + bool operator() (Point const& a, Point const& b) + { + // not remove this check or std::sort could generate + // a segmentation fault because it needs a strict '<' + // but due to round errors a == b doesn't mean dxy == dyx + if (a == b) return false; + Point da = a - o; + Point db = b - o; + if (da == -db) return false; + double dxy = da[X] * db[Y]; + double dyx = da[Y] * db[X]; + if (dxy > dyx) return true; + else if (dxy < dyx) return false; + return L2sq(da) < L2sq(db); + } +#endif }; void ConvexHull::find_pivot() { // Find pivot P; unsigned pivot = 0; - for(unsigned i = 1; i < boundary.size(); i++) + for (unsigned i = 1; i < boundary.size(); i++) if(boundary[i] <= boundary[pivot]) pivot = i; - + std::swap(boundary[0], boundary[pivot]); } @@ -111,12 +134,14 @@ ConvexHull::angle_sort() { std::sort(boundary.begin()+1, boundary.end(), angle_cmp(boundary[0])); } + void ConvexHull::graham_scan() { + if (boundary.size() < 4) return; unsigned stac = 2; - for(int i = 2; i < boundary.size(); i++) { - double o = SignedTriangleArea(boundary[stac-2], - boundary[stac-1], + for(unsigned int i = 2; i < boundary.size(); i++) { + double o = SignedTriangleArea(boundary[stac-2], + boundary[stac-1], boundary[i]); if(o == 0) { // colinear - dangerous... stac--; @@ -124,8 +149,8 @@ ConvexHull::graham_scan() { } else { // remove concavity while(o >= 0 && stac > 2) { stac--; - o = SignedTriangleArea(boundary[stac-2], - boundary[stac-1], + o = SignedTriangleArea(boundary[stac-2], + boundary[stac-1], boundary[i]); } } @@ -136,6 +161,8 @@ ConvexHull::graham_scan() { void ConvexHull::graham() { + if(is_degenerate()) // nothing to do + return; find_pivot(); angle_sort(); graham_scan(); @@ -143,7 +170,7 @@ ConvexHull::graham() { //Mathematically incorrect mod, but more useful. int mod(int i, int l) { - return i >= 0 ? + return i >= 0 ? i % l : (i % l) + l; } //OPT: usages can often be replaced by conditions @@ -152,6 +179,13 @@ 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; +} + +/*** 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; } @@ -166,7 +200,21 @@ ConvexHull::find_left(Point p) { } return -1; } -//OPT: do a spread iteration - quasi-random with no repeats and full coverage. + + +/*** ConvexHull::find_positive + * 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 +ConvexHull::find_strict_left(Point p) { + int l = boundary.size(); //Who knows if C++ is smart enough to optimize this? + for(int i = 0; i < l; i++) { + if(is_strict_left(p, i)) return i; + } + return -1; +} + +//OPT: do a spread iteration - quasi-random with no repeats and full coverage. /*** ConvexHull::contains_point * In order to test whether a point is inside a convex hull we can travel once around the outside making @@ -176,6 +224,14 @@ ConvexHull::contains_point(Point p) { return find_left(p) == -1; } +/*** ConvexHull::strict_contains_point + * In order to test whether a point is strictly inside (not on the boundary) a convex hull we can travel once around the outside making + * sure that each triangle made from an edge and the point has positive area. */ +bool +ConvexHull::strict_contains_point(Point p) { + return find_strict_left(p) == -1; +} + /*** ConvexHull::add_point * to add a point we need to find whether the new point extends the boundary, and if so, what it * obscures. Tarjan? Jarvis?*/ @@ -192,9 +248,9 @@ ConvexHull::merge(Point p) { bool pushed = false; - bool pre = is_left(p, -1); + bool pre = is_strict_left(p, -1); for(int i = 0; i < l; i++) { - bool cur = is_left(p, i); + bool cur = is_strict_left(p, i); if(pre) { if(cur) { if(!pushed) { @@ -211,7 +267,7 @@ ConvexHull::merge(Point p) { out.push_back(boundary[i]); pre = cur; } - + boundary = out; } //OPT: quickly find an obscured point and find the bounds by extending from there. then push all points not within the bounds in order. @@ -245,12 +301,12 @@ ConvexHull::is_clockwise() const { bool ConvexHull::top_point_first() const { std::vector::const_iterator pivot = boundary.begin(); - for(std::vector::const_iterator it(boundary.begin()+1), + for(std::vector::const_iterator it(boundary.begin()+1), e(boundary.end()); it != e; it++) { if((*it)[1] < (*pivot)[1]) pivot = it; - else if(((*it)[1] == (*pivot)[1]) && + else if(((*it)[1] == (*pivot)[1]) && ((*it)[0] < (*pivot)[0])) pivot = it; } @@ -264,7 +320,8 @@ proposed algorithm: We must be very careful about rounding here. */ bool ConvexHull::no_colinear_points() const { - + // XXX: implement me! + THROW_NOTIMPLEMENTED(); } bool @@ -292,8 +349,8 @@ bridges(ConvexHull a, ConvexHull b) { map abridges; map bbridges; - for(int ia = 0; ia < a.boundary.size(); ia++) { - for(int ib = 0; ib < b.boundary.size(); ib++) { + for(unsigned ia = 0; ia < a.boundary.size(); ia++) { + for(unsigned ib = 0; ib < b.boundary.size(); ib++) { Point d = b[ib] - a[ia]; Geom::Coord e = cross(d, a[ia - 1] - a[ia]), f = cross(d, a[ia + 1] - a[ia]); Geom::Coord g = cross(d, b[ib - 1] - a[ia]), h = cross(d, b[ib + 1] - a[ia]); @@ -301,7 +358,7 @@ bridges(ConvexHull a, ConvexHull b) { else if(e < 0 && f < 0 && g < 0 && h < 0) bbridges[ib] = ia; } } - + return make_pair(abridges, bbridges); } @@ -321,7 +378,7 @@ std::vector bridge_points(ConvexHull a, ConvexHull b) { unsigned find_bottom_right(ConvexHull const &a) { unsigned it = 1; - while(it < a.boundary.size() && + while(it < a.boundary.size() && a.boundary[it][Y] > a.boundary[it-1][Y]) it++; return it-1; @@ -335,10 +392,10 @@ unsigned find_bottom_right(ConvexHull const &a) { */ ConvexHull sweepline_intersection(ConvexHull const &a, ConvexHull const &b) { ConvexHull ret; - - int al = 0; - int bl = 0; - + + unsigned al = 0; + unsigned bl = 0; + while(al+1 < a.boundary.size() && (a.boundary[al+1][Y] > b.boundary[bl][Y])) { al++; @@ -348,8 +405,9 @@ ConvexHull sweepline_intersection(ConvexHull const &a, ConvexHull const &b) { bl++; } // al and bl now point to the top of the first pair of edges that overlap in y value - double sweep_y = std::min(a.boundary[al][Y], - b.boundary[bl][Y]); + //double sweep_y = std::min(a.boundary[al][Y], + // b.boundary[bl][Y]); + return ret; } /*** ConvexHull intersection(ConvexHull a, ConvexHull b); @@ -357,12 +415,13 @@ ConvexHull sweepline_intersection(ConvexHull const &a, ConvexHull const &b) { * (Proof: take any two points both in a and in b. Any point between them is in a by convexity, * and in b by convexity, thus in both. Need to prove still finite bounds.) */ -ConvexHull intersection(ConvexHull a, ConvexHull b) { +ConvexHull intersection(ConvexHull /*a*/, ConvexHull /*b*/) { ConvexHull ret; + /* int ai = 0, bi = 0; int aj = a.boundary.size() - 1; int bj = b.boundary.size() - 1; - + */ /*while (true) { if(a[ai] }*/ @@ -382,21 +441,21 @@ ConvexHull merge(ConvexHull a, ConvexHull b) { ab[-1] = 0; bb[-1] = 0; - int i = -1; + int i = -1; // XXX: i is int but refers to vector indices if(a.boundary[0][1] > b.boundary[0][1]) goto start_b; while(true) { for(; ab.count(i) == 0; i++) { ret.boundary.push_back(a[i]); - if(i >= a.boundary.size()) return ret; + if(i >= (int)a.boundary.size()) return ret; } if(ab[i] == 0 && i != -1) break; i = ab[i]; start_b: - + for(; bb.count(i) == 0; i++) { ret.boundary.push_back(b[i]); - if(i >= b.boundary.size()) return ret; + if(i >= (int)b.boundary.size()) return ret; } if(bb[i] == 0 && i != -1) break; i = bb[i]; @@ -406,21 +465,21 @@ ConvexHull merge(ConvexHull a, ConvexHull b) { ConvexHull graham_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(), + 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.angle_sort(); result.graham_scan(); - + return result; } //TODO: reinstate @@ -431,6 +490,66 @@ ConvexHull graham_merge(ConvexHull a, ConvexHull b) { } }*/ +double ConvexHull::centroid_and_area(Geom::Point& centroid) const { + const unsigned n = boundary.size(); + if (n < 2) + return 0; + if(n < 3) { + centroid = (boundary[0] + boundary[1])/2; + return 0; + } + 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]); + atmp += ai; + centroid_tmp += (boundary[j] + boundary[i])*ai; // first moment. + } + if (atmp != 0) { + centroid = centroid_tmp / (3 * atmp); + } + return atmp / 2; +} + +// TODO: This can be made lg(n) using golden section/fibonacci search three starting points, say 0, +// n/2, n-1 construct a new point, say (n/2 + n)/2 throw away the furthest boundary point iterate +// until interval is a single value +Point const * ConvexHull::furthest(Point direction) const { + Point const * p = &boundary[0]; + double d = dot(*p, direction); + for(unsigned i = 1; i < boundary.size(); i++) { + double dd = dot(boundary[i], direction); + if(d < dd) { + p = &boundary[i]; + d = dd; + } + } + return p; +} + + +// returns (a, (b,c)), three points which define the narrowest diameter of the hull as the pair of +// lines going through b,c, and through a, parallel to b,c TODO: This can be made linear time by +// moving point tc incrementally from the previous value (it can only move in one direction). It +// is currently n*O(furthest) +double ConvexHull::narrowest_diameter(Point &a, Point &b, Point &c) { + Point tb = boundary.back(); + double d = INFINITY; + for(unsigned i = 0; i < boundary.size(); i++) { + Point tc = boundary[i]; + Point n = -rot90(tb-tc); + Point ta = *furthest(n); + double td = dot(n, ta-tb)/dot(n,n); + if(td < d) { + a = ta; + b = tb; + c = tc; + d = td; + } + tb = tc; + } + return d; +} }; @@ -438,12 +557,9 @@ ConvexHull graham_merge(ConvexHull a, ConvexHull b) { Local Variables: mode:c++ c-file-style:"stroustrup" - c-file-offsets:((innamespace . 0)(substatement-open . 0)) + c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) indent-tabs-mode:nil - c-brace-offset:0 fill-column:99 End: - vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4 : */ - - +// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :