diff options
Diffstat (limited to '')
| -rw-r--r-- | src/2geom/convex-hull.cpp | 746 |
1 files changed, 746 insertions, 0 deletions
diff --git a/src/2geom/convex-hull.cpp b/src/2geom/convex-hull.cpp new file mode 100644 index 0000000..f801fcc --- /dev/null +++ b/src/2geom/convex-hull.cpp @@ -0,0 +1,746 @@ +/** @file + * @brief Convex hull of a set of points + *//* + * Authors: + * Nathan Hurst <njh@mail.csse.monash.edu.au> + * Michael G. Sloan <mgsloan@gmail.com> + * Krzysztof Kosiński <tweenk.pl@gmail.com> + * Copyright 2006-2015 Authors + * + * This library is free software; you can redistribute it and/or + * modify it either under the terms of the GNU Lesser General Public + * License version 2.1 as published by the Free Software Foundation + * (the "LGPL") or, at your option, under the terms of the Mozilla + * Public License Version 1.1 (the "MPL"). If you do not alter this + * notice, a recipient may use your version of this file under either + * the MPL or the LGPL. + * + * You should have received a copy of the LGPL along with this library + * in the file COPYING-LGPL-2.1; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * You should have received a copy of the MPL along with this library + * in the file COPYING-MPL-1.1 + * + * The contents of this file are subject to the Mozilla Public License + * Version 1.1 (the "License"); you may not use this file except in + * compliance with the License. You may obtain a copy of the License at + * http://www.mozilla.org/MPL/ + * + * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY + * OF ANY KIND, either express or implied. See the LGPL or the MPL for + * the specific language governing rights and limitations. + * + */ + +#include <2geom/convex-hull.h> +#include <2geom/exception.h> +#include <algorithm> +#include <map> +#include <iostream> +#include <cassert> +#include <boost/array.hpp> + +/** Todo: + + modify graham scan to work top to bottom, rather than around angles + + intersection + + minimum distance between convex hulls + + maximum distance between convex hulls + + hausdorf metric? + + check all degenerate cases carefully + + check all algorithms meet all invariants + + generalise rotating caliper algorithm (iterator/circulator?) +*/ + +using std::vector; +using std::map; +using std::pair; +using std::make_pair; +using std::swap; + +namespace Geom { + +ConvexHull::ConvexHull(Point const &a, Point const &b) + : _boundary(2) + , _lower(0) +{ + _boundary[0] = a; + _boundary[1] = b; + std::sort(_boundary.begin(), _boundary.end(), Point::LexLess<X>()); + _construct(); +} + +ConvexHull::ConvexHull(Point const &a, Point const &b, Point const &c) + : _boundary(3) + , _lower(0) +{ + _boundary[0] = a; + _boundary[1] = b; + _boundary[2] = c; + std::sort(_boundary.begin(), _boundary.end(), Point::LexLess<X>()); + _construct(); +} + +ConvexHull::ConvexHull(Point const &a, Point const &b, Point const &c, Point const &d) + : _boundary(4) + , _lower(0) +{ + _boundary[0] = a; + _boundary[1] = b; + _boundary[2] = c; + _boundary[3] = d; + std::sort(_boundary.begin(), _boundary.end(), Point::LexLess<X>()); + _construct(); +} + +ConvexHull::ConvexHull(std::vector<Point> const &pts) + : _lower(0) +{ + //if (pts.size() > 16) { // arbitrary threshold + // _prune(pts.begin(), pts.end(), _boundary); + //} else { + _boundary = pts; + std::sort(_boundary.begin(), _boundary.end(), Point::LexLess<X>()); + //} + _construct(); +} + +bool ConvexHull::_is_clockwise_turn(Point const &a, Point const &b, Point const &c) +{ + if (b == c) return false; + return cross(b-a, c-a) > 0; +} + +void ConvexHull::_construct() +{ + // _boundary must already be sorted in LexLess<X> order + if (_boundary.empty()) { + _lower = 0; + return; + } + if (_boundary.size() == 1 || (_boundary.size() == 2 && _boundary[0] == _boundary[1])) { + _boundary.resize(1); + _lower = 1; + return; + } + if (_boundary.size() == 2) { + _lower = 2; + return; + } + + std::size_t k = 2; + for (std::size_t i = 2; i < _boundary.size(); ++i) { + while (k >= 2 && !_is_clockwise_turn(_boundary[k-2], _boundary[k-1], _boundary[i])) { + --k; + } + std::swap(_boundary[k++], _boundary[i]); + } + + _lower = k; + std::sort(_boundary.begin() + k, _boundary.end(), Point::LexGreater<X>()); + _boundary.push_back(_boundary.front()); + for (std::size_t i = _lower; i < _boundary.size(); ++i) { + while (k > _lower && !_is_clockwise_turn(_boundary[k-2], _boundary[k-1], _boundary[i])) { + --k; + } + std::swap(_boundary[k++], _boundary[i]); + } + + _boundary.resize(k-1); +} + +double ConvexHull::area() const +{ + if (size() <= 2) return 0; + + double a = 0; + for (std::size_t i = 0; i < size()-1; ++i) { + a += cross(_boundary[i], _boundary[i+1]); + } + a += cross(_boundary.back(), _boundary.front()); + return fabs(a * 0.5); +} + +OptRect ConvexHull::bounds() const +{ + OptRect ret; + if (empty()) return ret; + ret = Rect(left(), top(), right(), bottom()); + return ret; +} + +Point ConvexHull::topPoint() const +{ + Point ret; + ret[Y] = std::numeric_limits<Coord>::infinity(); + + for (auto i : upperHull()) { + if (ret[Y] >= i.y()) { + ret = i; + } else { + break; + } + } + + return ret; +} + +Point ConvexHull::bottomPoint() const +{ + Point ret; + ret[Y] = -std::numeric_limits<Coord>::infinity(); + + for (auto j : lowerHull()) { + if (ret[Y] <= j.y()) { + ret = j; + } else { + break; + } + } + + return ret; +} + +template <typename Iter, typename Lex> +bool below_x_monotonic_polyline(Point const &p, Iter first, Iter last, Lex lex) +{ + typename Lex::Secondary above; + Iter f = std::lower_bound(first, last, p, lex); + if (f == last) return false; + if (f == first) { + if (p == *f) return true; + return false; + } + + Point a = *(f-1), b = *f; + if (a[X] == b[X]) { + if (above(p[Y], a[Y]) || above(b[Y], p[Y])) return false; + } else { + // TODO: maybe there is a more numerically stable method + Coord y = lerp((p[X] - a[X]) / (b[X] - a[X]), a[Y], b[Y]); + if (above(p[Y], y)) return false; + } + return true; +} + +bool ConvexHull::contains(Point const &p) const +{ + if (_boundary.empty()) return false; + if (_boundary.size() == 1) { + if (_boundary[0] == p) return true; + return false; + } + + // 1. verify that the point is in the relevant X range + if (p[X] < _boundary[0][X] || p[X] > _boundary[_lower-1][X]) return false; + + // 2. check whether it is below the upper hull + UpperIterator ub = upperHull().begin(), ue = upperHull().end(); + if (!below_x_monotonic_polyline(p, ub, ue, Point::LexLess<X>())) return false; + + // 3. check whether it is above the lower hull + LowerIterator lb = lowerHull().begin(), le = lowerHull().end(); + if (!below_x_monotonic_polyline(p, lb, le, Point::LexGreater<X>())) return false; + + return true; +} + +bool ConvexHull::contains(Rect const &r) const +{ + for (unsigned i = 0; i < 4; ++i) { + if (!contains(r.corner(i))) return false; + } + return true; +} + +bool ConvexHull::contains(ConvexHull const &ch) const +{ + // TODO: requires interiorContains. + // We have to check all points of ch, and each point takes logarithmic time. + // If there are more points in ch that here, it is faster to make the check + // the other way around. + /*if (ch.size() > size()) { + for (iterator i = begin(); i != end(); ++i) { + if (ch.interiorContains(*i)) return false; + } + return true; + }*/ + + for (auto i : ch) { + if (!contains(i)) return false; + } + return true; +} + +void ConvexHull::swap(ConvexHull &other) +{ + _boundary.swap(other._boundary); + std::swap(_lower, other._lower); +} + +void ConvexHull::swap(std::vector<Point> &pts) +{ + _boundary.swap(pts); + _lower = 0; + std::sort(_boundary.begin(), _boundary.end(), Point::LexLess<X>()); + _construct(); +} + +#if 0 +/*** SignedTriangleArea + * returns the area of the triangle defined by p0, p1, p2. A clockwise triangle has positive area. + */ +double +SignedTriangleArea(Point p0, Point p1, Point p2) { + return cross((p1 - p0), (p2 - p0)); +} + +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]; + if((da[1] == 0) && (db[1] == 0)) + return da[0] < db[0]; + if(da[1] == 0) + return true; // infinite tangent + if(db[1] == 0) + return false; // infinite tangent + aa = da[0] / da[1]; + ab = db[0] / db[1]; + if(aa > ab) + return true; +#else + //assert((ata > atb) == (aa < ab)); + double aa = atan2(da); + double ab = atan2(db); + if(aa < ab) + return true; +#endif + if(aa == ab) + 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 +}; + +//Mathematically incorrect mod, but more useful. +int mod(int i, int l) { + return i >= 0 ? + i % l : (i % l) + l; +} +//OPT: usages can often be replaced by conditions + +/*** 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?*/ +void +ConvexHull::merge(Point p) { + std::vector<Point> out; + + int len = boundary.size(); + + if(len < 2) { + if(boundary.empty() || boundary[0] != p) + boundary.push_back(p); + return; + } + + bool pushed = false; + + bool pre = is_left(p, -1); + for(int i = 0; i < len; i++) { + bool cur = is_left(p, i); + if(pre) { + if(cur) { + if(!pushed) { + out.push_back(p); + pushed = true; + } + continue; + } + else if(!pushed) { + out.push_back(p); + pushed = true; + } + } + 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. + //OPT: use binary searches to find the actual starts/ends, use known rights as boundaries. may require cooperation of find_left algo. + +/*** 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 +ConvexHull::is_clockwise() const { + if(is_degenerate()) + return true; + Point first = boundary[0]; + Point second = boundary[1]; + for(std::vector<Point>::const_iterator it(boundary.begin()+2), e(boundary.end()); + it != e;) { + if(SignedTriangleArea(first, second, *it) > 0) + return false; + first = second; + second = *it; + ++it; + } + return true; +} + +/*** ConvexHull::top_point_first + * We require that the first point in the convex hull has the least y coord, and that off all such points on the hull, it has the least x coord. + * proposed algorithm: track lexicographic minimum while walking the list. + */ +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()); + it != e; it++) { + if((*it)[1] < (*pivot)[1]) + pivot = it; + else if(((*it)[1] == (*pivot)[1]) && + ((*it)[0] < (*pivot)[0])) + pivot = it; + } + return pivot == boundary.begin(); +} +//OPT: since the Y values are orderly there should be something like a binary search to do this. + +bool +ConvexHull::meets_invariants() const { + return is_clockwise() && top_point_first(); +} + +/*** ConvexHull::is_degenerate + * We allow three degenerate cases: empty, 1 point and 2 points. In many cases these should be handled explicitly. + */ +bool +ConvexHull::is_degenerate() const { + return boundary.size() < 3; +} + + +int sgn(double x) { + if(x == 0) return 0; + return (x<0)?-1:1; +} + +bool same_side(Point L[2], Point xs[4]) { + int side = 0; + for(int i = 0; i < 4; i++) { + int sn = sgn(SignedTriangleArea(L[0], L[1], xs[i])); + if(sn && !side) + side = sn; + else if(sn != side) return false; + } + return true; +} + +/** find bridging pairs between two convex hulls. + * this code is based on Hormoz Pirzadeh's masters thesis. There is room for optimisation: + * 1. reduce recomputation + * 2. use more efficient angle code + * 3. write as iterator + */ +std::vector<pair<int, int> > bridges(ConvexHull a, ConvexHull b) { + vector<pair<int, int> > ret; + + // 1. find maximal points on a and b + int ai = 0, bi = 0; + // 2. find first copodal pair + double ap_angle = atan2(a[ai+1] - a[ai]); + double bp_angle = atan2(b[bi+1] - b[bi]); + Point L[2] = {a[ai], b[bi]}; + while(ai < int(a.size()) || bi < int(b.size())) { + if(ap_angle == bp_angle) { + // In the case of parallel support lines, we must consider all four pairs of copodal points + { + assert(0); // untested + Point xs[4] = {a[ai-1], a[ai+1], b[bi-1], b[bi+1]}; + if(same_side(L, xs)) ret.push_back(make_pair(ai, bi)); + xs[2] = b[bi]; + xs[3] = b[bi+2]; + if(same_side(L, xs)) ret.push_back(make_pair(ai, bi)); + xs[0] = a[ai]; + xs[1] = a[ai+2]; + if(same_side(L, xs)) ret.push_back(make_pair(ai, bi)); + xs[2] = b[bi-1]; + xs[3] = b[bi+1]; + if(same_side(L, xs)) ret.push_back(make_pair(ai, bi)); + } + ai++; + ap_angle += angle_between(a[ai] - a[ai-1], a[ai+1] - a[ai]); + L[0] = a[ai]; + bi++; + bp_angle += angle_between(b[bi] - b[bi-1], b[bi+1] - b[bi]); + L[1] = b[bi]; + 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]); + L[0] = a[ai]; + Point xs[4] = {a[ai-1], a[ai+1], b[bi-1], b[bi+1]}; + if(same_side(L, xs)) ret.push_back(make_pair(ai, bi)); + } else { + bi++; + bp_angle += angle_between(b[bi] - b[bi-1], b[bi+1] - b[bi]); + L[1] = b[bi]; + Point xs[4] = {a[ai-1], a[ai+1], b[bi-1], b[bi+1]}; + if(same_side(L, xs)) ret.push_back(make_pair(ai, bi)); + } + } + return ret; +} + +unsigned find_bottom_right(ConvexHull const &a) { + unsigned it = 1; + while(it < a.boundary.size() && + a.boundary[it][Y] > a.boundary[it-1][Y]) + it++; + return it-1; +} + +/*** ConvexHull sweepline_intersection(ConvexHull a, ConvexHull b); + * find the intersection between two convex hulls. The intersection is also a convex hull. + * (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.) + * This algorithm works by sweeping a line down both convex hulls in parallel, working out the left and right edges of the new hull. + */ +ConvexHull sweepline_intersection(ConvexHull const &a, ConvexHull const &b) { + ConvexHull ret; + + unsigned al = 0; + unsigned bl = 0; + + while(al+1 < a.boundary.size() && + (a.boundary[al+1][Y] > b.boundary[bl][Y])) { + al++; + } + while(bl+1 < b.boundary.size() && + (b.boundary[bl+1][Y] > a.boundary[al][Y])) { + 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]); + return ret; +} + +/*** ConvexHull intersection(ConvexHull a, ConvexHull b); + * find the intersection between two convex hulls. The intersection is also a convex hull. + * (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 ret; + /* + int ai = 0, bi = 0; + int aj = a.boundary.size() - 1; + int bj = b.boundary.size() - 1; + */ + /*while (true) { + if(a[ai] + }*/ + return ret; +} + +template <typename T> +T idx_to_pair(pair<T, T> p, int idx) { + return idx?p.second:p.first; +} + +/*** ConvexHull merge(ConvexHull a, ConvexHull b); + * find the smallest convex hull that surrounds a and b. + */ +ConvexHull merge(ConvexHull a, ConvexHull b) { + ConvexHull ret; + + std::cout << "---\n"; + std::vector<pair<int, int> > bpair = bridges(a, b); + + // Given our list of bridges {(pb1, qb1), ..., (pbk, qbk)} + // we start with the highest point in p0, q0, say it is p0. + // then the merged hull is p0, ..., pb1, qb1, ..., qb2, pb2, ... + // In other words, either of the two polygons vertices are added in order until the vertex coincides with a bridge point, at which point we swap. + + unsigned state = (a[0][Y] < b[0][Y])?0:1; + ret.boundary.reserve(a.size() + b.size()); + ConvexHull chs[2] = {a, b}; + unsigned idx = 0; + + for(unsigned k = 0; k < bpair.size(); k++) { + unsigned limit = idx_to_pair(bpair[k], state); + std::cout << bpair[k].first << " , " << bpair[k].second << "; " + << idx << ", " << limit << ", s: " + << state + << " \n"; + while(idx <= limit) { + ret.boundary.push_back(chs[state][idx++]); + } + state = 1-state; + idx = idx_to_pair(bpair[k], state); + } + while(idx < chs[state].size()) { + ret.boundary.push_back(chs[state][idx++]); + } + return ret; +} + +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]) + 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.angle_sort(); + result.graham_scan(); + + 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]) + 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()); + for(Geom::Path::const_iterator it(sp.begin()), end(sp.end()); it != end; ++it) { + cc.push_back(ConvexHull((*it).begin(), (*it).end())); + } +}*/ + +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 = std::numeric_limits<double>::max(); + 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; +} +#endif + +}; + +/* + Local Variables: + mode:c++ + c-file-style:"stroustrup" + c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) + indent-tabs-mode:nil + fill-column:99 + End: +*/ +// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 : |