![[tokkee]](http://tokkee.org/images/avatar.png){kind=avatar}
index 50c43e6d144de46e67904a21b7cf27156533dca5..e8ea2280d9d77ed808ef50e2e5f9a5efe8346d19 100644 (file)
*
*/
-#include "convex-cover.h"
+#include <2geom/convex-cover.h>
+#include <2geom/exception.h>
#include <algorithm>
#include <map>
+
/** Todo:
+ modify graham scan to work top to bottom, rather than around angles
+ intersection
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];
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]);
}
std::sort(boundary.begin()+1, boundary.end(), angle_cmp(boundary[0]));
}
+
void
ConvexHull::graham_scan() {
+ 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],
+ double o = SignedTriangleArea(boundary[stac-2],
+ boundary[stac-1],
boundary[i]);
if(o == 0) { // colinear - dangerous...
stac--;
} 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]);
}
}
void
ConvexHull::graham() {
+ if(is_degenerate()) // nothing to do
+ return;
find_pivot();
angle_sort();
graham_scan();
//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
* 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;
}
}
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
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?*/
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) {
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.
bool
ConvexHull::top_point_first() const {
std::vector<Point>::const_iterator pivot = boundary.begin();
- for(std::vector<Point>::const_iterator it(boundary.begin()+1),
+ for(std::vector<Point>::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;
}
*/
bool
ConvexHull::no_colinear_points() const {
-
+ // XXX: implement me!
+ THROW_NOTIMPLEMENTED();
}
bool
map<int, int> abridges;
map<int, int> 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]);
else if(e < 0 && f < 0 && g < 0 && h < 0) bbridges[ib] = ia;
}
}
-
+
return make_pair(abridges, bbridges);
}
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;
*/
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++;
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;
}
* (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;
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];
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
}
}*/
+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;
+}
};