| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | /** @file | ||
| 2 | * @brief Convex hull of a set of points | ||
| 3 | *//* | ||
| 4 | * Authors: | ||
| 5 | * Nathan Hurst <njh@mail.csse.monash.edu.au> | ||
| 6 | * Michael G. Sloan <mgsloan@gmail.com> | ||
| 7 | * Krzysztof Kosiński <tweenk.pl@gmail.com> | ||
| 8 | * Copyright 2006-2015 Authors | ||
| 9 | * | ||
| 10 | * This library is free software; you can redistribute it and/or | ||
| 11 | * modify it either under the terms of the GNU Lesser General Public | ||
| 12 | * License version 2.1 as published by the Free Software Foundation | ||
| 13 | * (the "LGPL") or, at your option, under the terms of the Mozilla | ||
| 14 | * Public License Version 1.1 (the "MPL"). If you do not alter this | ||
| 15 | * notice, a recipient may use your version of this file under either | ||
| 16 | * the MPL or the LGPL. | ||
| 17 | * | ||
| 18 | * You should have received a copy of the LGPL along with this library | ||
| 19 | * in the file COPYING-LGPL-2.1; if not, write to the Free Software | ||
| 20 | * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA | ||
| 21 | * You should have received a copy of the MPL along with this library | ||
| 22 | * in the file COPYING-MPL-1.1 | ||
| 23 | * | ||
| 24 | * The contents of this file are subject to the Mozilla Public License | ||
| 25 | * Version 1.1 (the "License"); you may not use this file except in | ||
| 26 | * compliance with the License. You may obtain a copy of the License at | ||
| 27 | * http://www.mozilla.org/MPL/ | ||
| 28 | * | ||
| 29 | * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY | ||
| 30 | * OF ANY KIND, either express or implied. See the LGPL or the MPL for | ||
| 31 | * the specific language governing rights and limitations. | ||
| 32 | * | ||
| 33 | */ | ||
| 34 | |||
| 35 | #include <2geom/convex-hull.h> | ||
| 36 | #include <2geom/exception.h> | ||
| 37 | #include <algorithm> | ||
| 38 | #include <map> | ||
| 39 | #include <iostream> | ||
| 40 | #include <cassert> | ||
| 41 | #include <boost/array.hpp> | ||
| 42 | |||
| 43 | /** Todo: | ||
| 44 | + modify graham scan to work top to bottom, rather than around angles | ||
| 45 | + intersection | ||
| 46 | + minimum distance between convex hulls | ||
| 47 | + maximum distance between convex hulls | ||
| 48 | + hausdorf metric? | ||
| 49 | + check all degenerate cases carefully | ||
| 50 | + check all algorithms meet all invariants | ||
| 51 | + generalise rotating caliper algorithm (iterator/circulator?) | ||
| 52 | */ | ||
| 53 | |||
| 54 | using std::vector; | ||
| 55 | using std::map; | ||
| 56 | using std::pair; | ||
| 57 | using std::make_pair; | ||
| 58 | using std::swap; | ||
| 59 | |||
| 60 | namespace Geom { | ||
| 61 | |||
| 62 | ✗ | ConvexHull::ConvexHull(Point const &a, Point const &b) | |
| 63 | ✗ | : _boundary(2) | |
| 64 | ✗ | , _lower(0) | |
| 65 | { | ||
| 66 | ✗ | _boundary[0] = a; | |
| 67 | ✗ | _boundary[1] = b; | |
| 68 | ✗ | std::sort(_boundary.begin(), _boundary.end(), Point::LexLess<X>()); | |
| 69 | ✗ | _construct(); | |
| 70 | ✗ | } | |
| 71 | |||
| 72 | ✗ | ConvexHull::ConvexHull(Point const &a, Point const &b, Point const &c) | |
| 73 | ✗ | : _boundary(3) | |
| 74 | ✗ | , _lower(0) | |
| 75 | { | ||
| 76 | ✗ | _boundary[0] = a; | |
| 77 | ✗ | _boundary[1] = b; | |
| 78 | ✗ | _boundary[2] = c; | |
| 79 | ✗ | std::sort(_boundary.begin(), _boundary.end(), Point::LexLess<X>()); | |
| 80 | ✗ | _construct(); | |
| 81 | ✗ | } | |
| 82 | |||
| 83 | ✗ | ConvexHull::ConvexHull(Point const &a, Point const &b, Point const &c, Point const &d) | |
| 84 | ✗ | : _boundary(4) | |
| 85 | ✗ | , _lower(0) | |
| 86 | { | ||
| 87 | ✗ | _boundary[0] = a; | |
| 88 | ✗ | _boundary[1] = b; | |
| 89 | ✗ | _boundary[2] = c; | |
| 90 | ✗ | _boundary[3] = d; | |
| 91 | ✗ | std::sort(_boundary.begin(), _boundary.end(), Point::LexLess<X>()); | |
| 92 | ✗ | _construct(); | |
| 93 | ✗ | } | |
| 94 | |||
| 95 | 130 | ConvexHull::ConvexHull(std::vector<Point> const &pts) | |
| 96 | 130 | : _lower(0) | |
| 97 | { | ||
| 98 | //if (pts.size() > 16) { // arbitrary threshold | ||
| 99 | // _prune(pts.begin(), pts.end(), _boundary); | ||
| 100 | //} else { | ||
| 101 |
1/2✓ Branch 1 taken 130 times.
✗ Branch 2 not taken.
|
130 | _boundary = pts; |
| 102 |
1/2✓ Branch 3 taken 130 times.
✗ Branch 4 not taken.
|
130 | std::sort(_boundary.begin(), _boundary.end(), Point::LexLess<X>()); |
| 103 | //} | ||
| 104 |
1/2✓ Branch 1 taken 130 times.
✗ Branch 2 not taken.
|
130 | _construct(); |
| 105 | 130 | } | |
| 106 | |||
| 107 | 6181 | static bool is_clockwise_turn(Point const &a, Point const &b, Point const &c) | |
| 108 | { | ||
| 109 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 6181 times.
|
6181 | if (b == c) return false; |
| 110 | 6181 | return cross(b-a, c-a) > 0; | |
| 111 | } | ||
| 112 | |||
| 113 | 913 | void ConvexHull::_construct() | |
| 114 | { | ||
| 115 | // _boundary must already be sorted in LexLess<X> order | ||
| 116 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 913 times.
|
913 | if (_boundary.empty()) { |
| 117 | ✗ | _lower = 0; | |
| 118 | ✗ | return; | |
| 119 | } | ||
| 120 |
7/8✓ Branch 1 taken 908 times.
✓ Branch 2 taken 5 times.
✓ Branch 4 taken 5 times.
✓ Branch 5 taken 903 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 5 times.
✓ Branch 11 taken 5 times.
✓ Branch 12 taken 908 times.
|
913 | if (_boundary.size() == 1 || (_boundary.size() == 2 && _boundary[0] == _boundary[1])) { |
| 121 | 5 | _boundary.resize(1); | |
| 122 | 5 | _lower = 1; | |
| 123 | 5 | return; | |
| 124 | } | ||
| 125 |
2/2✓ Branch 1 taken 5 times.
✓ Branch 2 taken 903 times.
|
908 | if (_boundary.size() == 2) { |
| 126 | 5 | _lower = 2; | |
| 127 | 5 | return; | |
| 128 | } | ||
| 129 | |||
| 130 | 903 | std::size_t k = 2; | |
| 131 |
2/2✓ Branch 1 taken 2834 times.
✓ Branch 2 taken 903 times.
|
3737 | for (std::size_t i = 2; i < _boundary.size(); ++i) { |
| 132 |
6/6✓ Branch 0 taken 3628 times.
✓ Branch 1 taken 1178 times.
✓ Branch 6 taken 1972 times.
✓ Branch 7 taken 1656 times.
✓ Branch 8 taken 1972 times.
✓ Branch 9 taken 2834 times.
|
4806 | while (k >= 2 && !is_clockwise_turn(_boundary[k-2], _boundary[k-1], _boundary[i])) { |
| 133 | 1972 | --k; | |
| 134 | } | ||
| 135 | 2834 | std::swap(_boundary[k++], _boundary[i]); | |
| 136 | } | ||
| 137 | |||
| 138 | 903 | _lower = k; | |
| 139 |
1/2✓ Branch 5 taken 903 times.
✗ Branch 6 not taken.
|
903 | std::sort(_boundary.begin() + k, _boundary.end(), Point::LexGreater<X>()); |
| 140 | 903 | _boundary.push_back(_boundary.front()); | |
| 141 |
2/2✓ Branch 1 taken 2875 times.
✓ Branch 2 taken 903 times.
|
3778 | for (std::size_t i = _lower; i < _boundary.size(); ++i) { |
| 142 |
6/6✓ Branch 0 taken 2553 times.
✓ Branch 1 taken 1368 times.
✓ Branch 6 taken 1046 times.
✓ Branch 7 taken 1507 times.
✓ Branch 8 taken 1046 times.
✓ Branch 9 taken 2875 times.
|
3921 | while (k > _lower && !is_clockwise_turn(_boundary[k-2], _boundary[k-1], _boundary[i])) { |
| 143 | 1046 | --k; | |
| 144 | } | ||
| 145 | 2875 | std::swap(_boundary[k++], _boundary[i]); | |
| 146 | } | ||
| 147 | |||
| 148 | 903 | _boundary.resize(k-1); | |
| 149 | } | ||
| 150 | |||
| 151 | 9 | double ConvexHull::area() const | |
| 152 | { | ||
| 153 |
2/2✓ Branch 1 taken 3 times.
✓ Branch 2 taken 6 times.
|
9 | if (size() <= 2) return 0; |
| 154 | |||
| 155 | 6 | double a = 0; | |
| 156 |
2/2✓ Branch 1 taken 29 times.
✓ Branch 2 taken 6 times.
|
35 | for (std::size_t i = 0; i < size()-1; ++i) { |
| 157 | 29 | a += cross(_boundary[i], _boundary[i+1]); | |
| 158 | } | ||
| 159 | 6 | a += cross(_boundary.back(), _boundary.front()); | |
| 160 | 6 | return fabs(a * 0.5); | |
| 161 | } | ||
| 162 | |||
| 163 | 7 | OptRect ConvexHull::bounds() const | |
| 164 | { | ||
| 165 |
2/2✓ Branch 1 taken 1 times.
✓ Branch 2 taken 6 times.
|
7 | if (empty()) return {}; |
| 166 |
4/8✓ Branch 2 taken 6 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 6 times.
✗ Branch 7 not taken.
✓ Branch 10 taken 6 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 6 times.
✗ Branch 14 not taken.
|
12 | return Rect(left(), top(), right(), bottom()); |
| 167 | } | ||
| 168 | |||
| 169 | 102 | std::pair<Rotate, OptRect> ConvexHull::minAreaRotation() const | |
| 170 | { | ||
| 171 |
2/2✓ Branch 1 taken 1 times.
✓ Branch 2 taken 101 times.
|
102 | if (empty()) { |
| 172 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | return {}; |
| 173 |
2/2✓ Branch 1 taken 1 times.
✓ Branch 2 taken 100 times.
|
101 | } else if (_boundary.size() == 1) { |
| 174 |
3/6✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
|
1 | return { {}, Rect(front(), front()) }; |
| 175 | } | ||
| 176 | |||
| 177 | // Move the point i cyclically along until it maximises distance in the direction n. | ||
| 178 | 100 | auto advance = [this] (size_t i, Point const &n) { | |
| 179 | 1629 | auto ih = dot(_boundary[i], n); | |
| 180 | while (true) { | ||
| 181 | 3318 | auto j = (i + 1) % size(); | |
| 182 | 3318 | auto jh = dot(_boundary[j], n); | |
| 183 |
2/2✓ Branch 0 taken 1629 times.
✓ Branch 1 taken 1689 times.
|
3318 | if (ih >= jh) { |
| 184 | 1629 | return i; | |
| 185 | } | ||
| 186 | 1689 | i = j; | |
| 187 | 1689 | ih = jh; | |
| 188 | 1689 | } | |
| 189 | 100 | }; | |
| 190 | |||
| 191 | 100 | auto min_area = std::numeric_limits<Coord>::max(); | |
| 192 | 100 | Point min_v; | |
| 193 | 100 | Rect min_rect; | |
| 194 | |||
| 195 | // Run rotating callipers. | ||
| 196 | size_t j, k, l; | ||
| 197 |
2/2✓ Branch 1 taken 543 times.
✓ Branch 2 taken 100 times.
|
643 | for (size_t i = 0; i < size(); i++) { |
| 198 | // Get the current segment. | ||
| 199 | 543 | auto &p1 = _boundary[i]; | |
| 200 | 543 | auto &p2 = _boundary[(i + 1) % size()]; | |
| 201 | 543 | auto const v = p2 - p1; | |
| 202 | 543 | auto const n = v.cw(); | |
| 203 | |||
| 204 |
2/2✓ Branch 0 taken 100 times.
✓ Branch 1 taken 443 times.
|
543 | if (i == 0) { |
| 205 | // Initialise the points. | ||
| 206 | 100 | j = advance(0, v); | |
| 207 | 100 | k = advance(j, n); | |
| 208 | 100 | l = advance(k, -v); | |
| 209 | } else { | ||
| 210 | // Advance the points. | ||
| 211 | 443 | j = advance(j, v); | |
| 212 | 443 | k = advance(k, n); | |
| 213 | 443 | l = advance(l, -v); | |
| 214 | } | ||
| 215 | |||
| 216 | // Compute the dimensions of the unconstrained rectangle, times v.length(). | ||
| 217 | 543 | auto const w = dot(_boundary[j] - _boundary[l], v); | |
| 218 | 543 | auto const h = dot(_boundary[k] - _boundary[i], n); | |
| 219 | 543 | auto const area = w * h / v.lengthSq(); | |
| 220 | |||
| 221 | // Track the minimum. | ||
| 222 |
2/2✓ Branch 0 taken 216 times.
✓ Branch 1 taken 327 times.
|
543 | if (area < min_area) { |
| 223 | 216 | min_area = area; | |
| 224 | 216 | min_v = v; | |
| 225 |
1/2✓ Branch 5 taken 216 times.
✗ Branch 6 not taken.
|
216 | min_rect = Rect::from_xywh(dot(_boundary[l], v), dot(_boundary[i], n), w, h); |
| 226 | } | ||
| 227 | } | ||
| 228 | |||
| 229 |
4/8✓ Branch 3 taken 100 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 100 times.
✗ Branch 7 not taken.
✓ Branch 15 taken 100 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 100 times.
✗ Branch 19 not taken.
|
200 | return { Rotate(min_v).inverse(), min_rect * Scale(1.0 / min_v.length()) }; |
| 230 | } | ||
| 231 | |||
| 232 | 18 | Point ConvexHull::topPoint() const | |
| 233 | { | ||
| 234 | 18 | Point ret; | |
| 235 | 18 | ret[Y] = std::numeric_limits<Coord>::infinity(); | |
| 236 | |||
| 237 |
3/4✓ Branch 2 taken 18 times.
✗ Branch 3 not taken.
✓ Branch 12 taken 48 times.
✓ Branch 13 taken 9 times.
|
57 | for (auto i : upperHull()) { |
| 238 |
2/2✓ Branch 2 taken 39 times.
✓ Branch 3 taken 9 times.
|
48 | if (ret[Y] >= i.y()) { |
| 239 | 39 | ret = i; | |
| 240 | } else { | ||
| 241 | 9 | break; | |
| 242 | } | ||
| 243 | } | ||
| 244 | |||
| 245 | 18 | return ret; | |
| 246 | } | ||
| 247 | |||
| 248 | 18 | Point ConvexHull::bottomPoint() const | |
| 249 | { | ||
| 250 | 18 | Point ret; | |
| 251 | 18 | ret[Y] = -std::numeric_limits<Coord>::infinity(); | |
| 252 | |||
| 253 |
4/6✓ Branch 2 taken 18 times.
✗ Branch 3 not taken.
✓ Branch 11 taken 39 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 51 times.
✓ Branch 15 taken 6 times.
|
96 | for (auto j : lowerHull()) { |
| 254 |
2/2✓ Branch 2 taken 39 times.
✓ Branch 3 taken 12 times.
|
51 | if (ret[Y] <= j.y()) { |
| 255 | 39 | ret = j; | |
| 256 | } else { | ||
| 257 | 12 | break; | |
| 258 | } | ||
| 259 | } | ||
| 260 | |||
| 261 | 18 | return ret; | |
| 262 | } | ||
| 263 | |||
| 264 | template <typename Iter, typename Lex> | ||
| 265 | 256 | bool below_x_monotonic_polyline(Point const &p, Iter first, Iter last, Lex lex) | |
| 266 | { | ||
| 267 | 256 | typename Lex::Secondary above; | |
| 268 |
1/2✓ Branch 2 taken 128 times.
✗ Branch 3 not taken.
|
256 | Iter f = std::lower_bound(first, last, p, lex); |
| 269 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 128 times.
|
256 | if (f == last) return false; |
| 270 |
2/2✓ Branch 1 taken 7 times.
✓ Branch 2 taken 121 times.
|
256 | if (f == first) { |
| 271 |
1/2✓ Branch 2 taken 7 times.
✗ Branch 3 not taken.
|
14 | if (p == *f) return true; |
| 272 | ✗ | return false; | |
| 273 | } | ||
| 274 | |||
| 275 | 242 | Point a = *(f-1), b = *f; | |
| 276 |
2/2✓ Branch 2 taken 6 times.
✓ Branch 3 taken 115 times.
|
242 | if (a[X] == b[X]) { |
| 277 |
5/10✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 6 times.
✓ Branch 12 taken 6 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 6 times.
✗ Branch 16 not taken.
✗ Branch 18 not taken.
✓ Branch 19 taken 6 times.
|
12 | if (above(p[Y], a[Y]) || above(b[Y], p[Y])) return false; |
| 278 | } else { | ||
| 279 | // TODO: maybe there is a more numerically stable method | ||
| 280 | 230 | Coord y = lerp((p[X] - a[X]) / (b[X] - a[X]), a[Y], b[Y]); | |
| 281 |
2/2✓ Branch 4 taken 1 times.
✓ Branch 5 taken 114 times.
|
230 | if (above(p[Y], y)) return false; |
| 282 | } | ||
| 283 | 240 | return true; | |
| 284 | } | ||
| 285 | |||
| 286 | 66 | bool ConvexHull::contains(Point const &p) const | |
| 287 | { | ||
| 288 |
2/2✓ Branch 1 taken 1 times.
✓ Branch 2 taken 65 times.
|
66 | if (_boundary.empty()) return false; |
| 289 |
2/2✓ Branch 1 taken 1 times.
✓ Branch 2 taken 64 times.
|
65 | if (_boundary.size() == 1) { |
| 290 |
1/2✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
|
1 | if (_boundary[0] == p) return true; |
| 291 | ✗ | return false; | |
| 292 | } | ||
| 293 | |||
| 294 | // 1. verify that the point is in the relevant X range | ||
| 295 |
3/6✓ Branch 3 taken 64 times.
✗ Branch 4 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 64 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 64 times.
|
64 | if (p[X] < _boundary[0][X] || p[X] > _boundary[_lower-1][X]) return false; |
| 296 | |||
| 297 | // 2. check whether it is below the upper hull | ||
| 298 |
2/4✓ Branch 3 taken 64 times.
✗ Branch 4 not taken.
✓ Branch 10 taken 64 times.
✗ Branch 11 not taken.
|
64 | UpperIterator ub = upperHull().begin(), ue = upperHull().end(); |
| 299 |
2/4✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 64 times.
|
64 | if (!below_x_monotonic_polyline(p, ub, ue, Point::LexLess<X>())) return false; |
| 300 | |||
| 301 | // 3. check whether it is above the lower hull | ||
| 302 |
2/4✓ Branch 3 taken 64 times.
✗ Branch 4 not taken.
✓ Branch 10 taken 64 times.
✗ Branch 11 not taken.
|
64 | LowerIterator lb = lowerHull().begin(), le = lowerHull().end(); |
| 303 |
3/4✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 63 times.
|
64 | if (!below_x_monotonic_polyline(p, lb, le, Point::LexGreater<X>())) return false; |
| 304 | |||
| 305 | 63 | return true; | |
| 306 | } | ||
| 307 | |||
| 308 | ✗ | bool ConvexHull::contains(Rect const &r) const | |
| 309 | { | ||
| 310 | ✗ | for (unsigned i = 0; i < 4; ++i) { | |
| 311 | ✗ | if (!contains(r.corner(i))) return false; | |
| 312 | } | ||
| 313 | ✗ | return true; | |
| 314 | } | ||
| 315 | |||
| 316 | ✗ | bool ConvexHull::contains(ConvexHull const &ch) const | |
| 317 | { | ||
| 318 | // TODO: requires interiorContains. | ||
| 319 | // We have to check all points of ch, and each point takes logarithmic time. | ||
| 320 | // If there are more points in ch that here, it is faster to make the check | ||
| 321 | // the other way around. | ||
| 322 | /*if (ch.size() > size()) { | ||
| 323 | for (iterator i = begin(); i != end(); ++i) { | ||
| 324 | if (ch.interiorContains(*i)) return false; | ||
| 325 | } | ||
| 326 | return true; | ||
| 327 | }*/ | ||
| 328 | |||
| 329 | ✗ | for (auto i : ch) { | |
| 330 | ✗ | if (!contains(i)) return false; | |
| 331 | } | ||
| 332 | ✗ | return true; | |
| 333 | } | ||
| 334 | |||
| 335 | ✗ | void ConvexHull::swap(ConvexHull &other) | |
| 336 | { | ||
| 337 | ✗ | _boundary.swap(other._boundary); | |
| 338 | ✗ | std::swap(_lower, other._lower); | |
| 339 | ✗ | } | |
| 340 | |||
| 341 | 783 | void ConvexHull::swap(std::vector<Point> &pts) | |
| 342 | { | ||
| 343 | 783 | _boundary.swap(pts); | |
| 344 | 783 | _lower = 0; | |
| 345 | 783 | std::sort(_boundary.begin(), _boundary.end(), Point::LexLess<X>()); | |
| 346 | 783 | _construct(); | |
| 347 | 783 | } | |
| 348 | |||
| 349 | #if 0 | ||
| 350 | /*** SignedTriangleArea | ||
| 351 | * returns the area of the triangle defined by p0, p1, p2. A clockwise triangle has positive area. | ||
| 352 | */ | ||
| 353 | double | ||
| 354 | SignedTriangleArea(Point p0, Point p1, Point p2) { | ||
| 355 | return cross((p1 - p0), (p2 - p0)); | ||
| 356 | } | ||
| 357 | |||
| 358 | class angle_cmp{ | ||
| 359 | public: | ||
| 360 | Point o; | ||
| 361 | angle_cmp(Point o) : o(o) {} | ||
| 362 | |||
| 363 | #if 0 | ||
| 364 | bool | ||
| 365 | operator()(Point a, Point b) { | ||
| 366 | // not remove this check or std::sort could crash | ||
| 367 | if (a == b) return false; | ||
| 368 | Point da = a - o; | ||
| 369 | Point db = b - o; | ||
| 370 | if (da == -db) return false; | ||
| 371 | |||
| 372 | #if 1 | ||
| 373 | double aa = da[0]; | ||
| 374 | double ab = db[0]; | ||
| 375 | if((da[1] == 0) && (db[1] == 0)) | ||
| 376 | return da[0] < db[0]; | ||
| 377 | if(da[1] == 0) | ||
| 378 | return true; // infinite tangent | ||
| 379 | if(db[1] == 0) | ||
| 380 | return false; // infinite tangent | ||
| 381 | aa = da[0] / da[1]; | ||
| 382 | ab = db[0] / db[1]; | ||
| 383 | if(aa > ab) | ||
| 384 | return true; | ||
| 385 | #else | ||
| 386 | //assert((ata > atb) == (aa < ab)); | ||
| 387 | double aa = atan2(da); | ||
| 388 | double ab = atan2(db); | ||
| 389 | if(aa < ab) | ||
| 390 | return true; | ||
| 391 | #endif | ||
| 392 | if(aa == ab) | ||
| 393 | return L2sq(da) < L2sq(db); | ||
| 394 | return false; | ||
| 395 | } | ||
| 396 | #else | ||
| 397 | bool operator() (Point const& a, Point const& b) | ||
| 398 | { | ||
| 399 | // not remove this check or std::sort could generate | ||
| 400 | // a segmentation fault because it needs a strict '<' | ||
| 401 | // but due to round errors a == b doesn't mean dxy == dyx | ||
| 402 | if (a == b) return false; | ||
| 403 | Point da = a - o; | ||
| 404 | Point db = b - o; | ||
| 405 | if (da == -db) return false; | ||
| 406 | double dxy = da[X] * db[Y]; | ||
| 407 | double dyx = da[Y] * db[X]; | ||
| 408 | if (dxy > dyx) return true; | ||
| 409 | else if (dxy < dyx) return false; | ||
| 410 | return L2sq(da) < L2sq(db); | ||
| 411 | } | ||
| 412 | #endif | ||
| 413 | }; | ||
| 414 | |||
| 415 | //Mathematically incorrect mod, but more useful. | ||
| 416 | int mod(int i, int l) { | ||
| 417 | return i >= 0 ? | ||
| 418 | i % l : (i % l) + l; | ||
| 419 | } | ||
| 420 | //OPT: usages can often be replaced by conditions | ||
| 421 | |||
| 422 | /*** ConvexHull::add_point | ||
| 423 | * to add a point we need to find whether the new point extends the boundary, and if so, what it | ||
| 424 | * obscures. Tarjan? Jarvis?*/ | ||
| 425 | void | ||
| 426 | ConvexHull::merge(Point p) { | ||
| 427 | std::vector<Point> out; | ||
| 428 | |||
| 429 | int len = boundary.size(); | ||
| 430 | |||
| 431 | if(len < 2) { | ||
| 432 | if(boundary.empty() || boundary[0] != p) | ||
| 433 | boundary.push_back(p); | ||
| 434 | return; | ||
| 435 | } | ||
| 436 | |||
| 437 | bool pushed = false; | ||
| 438 | |||
| 439 | bool pre = is_left(p, -1); | ||
| 440 | for(int i = 0; i < len; i++) { | ||
| 441 | bool cur = is_left(p, i); | ||
| 442 | if(pre) { | ||
| 443 | if(cur) { | ||
| 444 | if(!pushed) { | ||
| 445 | out.push_back(p); | ||
| 446 | pushed = true; | ||
| 447 | } | ||
| 448 | continue; | ||
| 449 | } | ||
| 450 | else if(!pushed) { | ||
| 451 | out.push_back(p); | ||
| 452 | pushed = true; | ||
| 453 | } | ||
| 454 | } | ||
| 455 | out.push_back(boundary[i]); | ||
| 456 | pre = cur; | ||
| 457 | } | ||
| 458 | |||
| 459 | boundary = out; | ||
| 460 | } | ||
| 461 | //OPT: quickly find an obscured point and find the bounds by extending from there. then push all points not within the bounds in order. | ||
| 462 | //OPT: use binary searches to find the actual starts/ends, use known rights as boundaries. may require cooperation of find_left algo. | ||
| 463 | |||
| 464 | /*** ConvexHull::is_clockwise | ||
| 465 | * We require that successive pairs of edges always turn right. | ||
| 466 | * We return false on collinear points | ||
| 467 | * proposed algorithm: walk successive edges and require triangle area is positive. | ||
| 468 | */ | ||
| 469 | bool | ||
| 470 | ConvexHull::is_clockwise() const { | ||
| 471 | if(is_degenerate()) | ||
| 472 | return true; | ||
| 473 | Point first = boundary[0]; | ||
| 474 | Point second = boundary[1]; | ||
| 475 | for(std::vector<Point>::const_iterator it(boundary.begin()+2), e(boundary.end()); | ||
| 476 | it != e;) { | ||
| 477 | if(SignedTriangleArea(first, second, *it) > 0) | ||
| 478 | return false; | ||
| 479 | first = second; | ||
| 480 | second = *it; | ||
| 481 | ++it; | ||
| 482 | } | ||
| 483 | return true; | ||
| 484 | } | ||
| 485 | |||
| 486 | /*** ConvexHull::top_point_first | ||
| 487 | * 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. | ||
| 488 | * proposed algorithm: track lexicographic minimum while walking the list. | ||
| 489 | */ | ||
| 490 | bool | ||
| 491 | ConvexHull::top_point_first() const { | ||
| 492 | if(size() <= 1) return true; | ||
| 493 | std::vector<Point>::const_iterator pivot = boundary.begin(); | ||
| 494 | for(std::vector<Point>::const_iterator it(boundary.begin()+1), | ||
| 495 | e(boundary.end()); | ||
| 496 | it != e; it++) { | ||
| 497 | if((*it)[1] < (*pivot)[1]) | ||
| 498 | pivot = it; | ||
| 499 | else if(((*it)[1] == (*pivot)[1]) && | ||
| 500 | ((*it)[0] < (*pivot)[0])) | ||
| 501 | pivot = it; | ||
| 502 | } | ||
| 503 | return pivot == boundary.begin(); | ||
| 504 | } | ||
| 505 | //OPT: since the Y values are orderly there should be something like a binary search to do this. | ||
| 506 | |||
| 507 | bool | ||
| 508 | ConvexHull::meets_invariants() const { | ||
| 509 | return is_clockwise() && top_point_first(); | ||
| 510 | } | ||
| 511 | |||
| 512 | /*** ConvexHull::is_degenerate | ||
| 513 | * We allow three degenerate cases: empty, 1 point and 2 points. In many cases these should be handled explicitly. | ||
| 514 | */ | ||
| 515 | bool | ||
| 516 | ConvexHull::is_degenerate() const { | ||
| 517 | return boundary.size() < 3; | ||
| 518 | } | ||
| 519 | |||
| 520 | |||
| 521 | int sgn(double x) { | ||
| 522 | if(x == 0) return 0; | ||
| 523 | return (x<0)?-1:1; | ||
| 524 | } | ||
| 525 | |||
| 526 | bool same_side(Point L[2], Point xs[4]) { | ||
| 527 | int side = 0; | ||
| 528 | for(int i = 0; i < 4; i++) { | ||
| 529 | int sn = sgn(SignedTriangleArea(L[0], L[1], xs[i])); | ||
| 530 | if(sn && !side) | ||
| 531 | side = sn; | ||
| 532 | else if(sn != side) return false; | ||
| 533 | } | ||
| 534 | return true; | ||
| 535 | } | ||
| 536 | |||
| 537 | /** find bridging pairs between two convex hulls. | ||
| 538 | * this code is based on Hormoz Pirzadeh's masters thesis. There is room for optimisation: | ||
| 539 | * 1. reduce recomputation | ||
| 540 | * 2. use more efficient angle code | ||
| 541 | * 3. write as iterator | ||
| 542 | */ | ||
| 543 | std::vector<pair<int, int> > bridges(ConvexHull a, ConvexHull b) { | ||
| 544 | vector<pair<int, int> > ret; | ||
| 545 | |||
| 546 | // 1. find maximal points on a and b | ||
| 547 | int ai = 0, bi = 0; | ||
| 548 | // 2. find first copodal pair | ||
| 549 | double ap_angle = atan2(a[ai+1] - a[ai]); | ||
| 550 | double bp_angle = atan2(b[bi+1] - b[bi]); | ||
| 551 | Point L[2] = {a[ai], b[bi]}; | ||
| 552 | while(ai < int(a.size()) || bi < int(b.size())) { | ||
| 553 | if(ap_angle == bp_angle) { | ||
| 554 | // In the case of parallel support lines, we must consider all four pairs of copodal points | ||
| 555 | { | ||
| 556 | assert(0); // untested | ||
| 557 | Point xs[4] = {a[ai-1], a[ai+1], b[bi-1], b[bi+1]}; | ||
| 558 | if(same_side(L, xs)) ret.push_back(make_pair(ai, bi)); | ||
| 559 | xs[2] = b[bi]; | ||
| 560 | xs[3] = b[bi+2]; | ||
| 561 | if(same_side(L, xs)) ret.push_back(make_pair(ai, bi)); | ||
| 562 | xs[0] = a[ai]; | ||
| 563 | xs[1] = a[ai+2]; | ||
| 564 | if(same_side(L, xs)) ret.push_back(make_pair(ai, bi)); | ||
| 565 | xs[2] = b[bi-1]; | ||
| 566 | xs[3] = b[bi+1]; | ||
| 567 | if(same_side(L, xs)) ret.push_back(make_pair(ai, bi)); | ||
| 568 | } | ||
| 569 | ai++; | ||
| 570 | ap_angle += angle_between(a[ai] - a[ai-1], a[ai+1] - a[ai]); | ||
| 571 | L[0] = a[ai]; | ||
| 572 | bi++; | ||
| 573 | bp_angle += angle_between(b[bi] - b[bi-1], b[bi+1] - b[bi]); | ||
| 574 | L[1] = b[bi]; | ||
| 575 | std::cout << "parallel\n"; | ||
| 576 | } else if(ap_angle < bp_angle) { | ||
| 577 | ai++; | ||
| 578 | ap_angle += angle_between(a[ai] - a[ai-1], a[ai+1] - a[ai]); | ||
| 579 | L[0] = a[ai]; | ||
| 580 | Point xs[4] = {a[ai-1], a[ai+1], b[bi-1], b[bi+1]}; | ||
| 581 | if(same_side(L, xs)) ret.push_back(make_pair(ai, bi)); | ||
| 582 | } else { | ||
| 583 | bi++; | ||
| 584 | bp_angle += angle_between(b[bi] - b[bi-1], b[bi+1] - b[bi]); | ||
| 585 | L[1] = b[bi]; | ||
| 586 | Point xs[4] = {a[ai-1], a[ai+1], b[bi-1], b[bi+1]}; | ||
| 587 | if(same_side(L, xs)) ret.push_back(make_pair(ai, bi)); | ||
| 588 | } | ||
| 589 | } | ||
| 590 | return ret; | ||
| 591 | } | ||
| 592 | |||
| 593 | unsigned find_bottom_right(ConvexHull const &a) { | ||
| 594 | unsigned it = 1; | ||
| 595 | while(it < a.boundary.size() && | ||
| 596 | a.boundary[it][Y] > a.boundary[it-1][Y]) | ||
| 597 | it++; | ||
| 598 | return it-1; | ||
| 599 | } | ||
| 600 | |||
| 601 | /*** ConvexHull sweepline_intersection(ConvexHull a, ConvexHull b); | ||
| 602 | * find the intersection between two convex hulls. The intersection is also a convex hull. | ||
| 603 | * (Proof: take any two points both in a and in b. Any point between them is in a by convexity, | ||
| 604 | * and in b by convexity, thus in both. Need to prove still finite bounds.) | ||
| 605 | * This algorithm works by sweeping a line down both convex hulls in parallel, working out the left and right edges of the new hull. | ||
| 606 | */ | ||
| 607 | ConvexHull sweepline_intersection(ConvexHull const &a, ConvexHull const &b) { | ||
| 608 | ConvexHull ret; | ||
| 609 | |||
| 610 | unsigned al = 0; | ||
| 611 | unsigned bl = 0; | ||
| 612 | |||
| 613 | while(al+1 < a.boundary.size() && | ||
| 614 | (a.boundary[al+1][Y] > b.boundary[bl][Y])) { | ||
| 615 | al++; | ||
| 616 | } | ||
| 617 | while(bl+1 < b.boundary.size() && | ||
| 618 | (b.boundary[bl+1][Y] > a.boundary[al][Y])) { | ||
| 619 | bl++; | ||
| 620 | } | ||
| 621 | // al and bl now point to the top of the first pair of edges that overlap in y value | ||
| 622 | //double sweep_y = std::min(a.boundary[al][Y], | ||
| 623 | // b.boundary[bl][Y]); | ||
| 624 | return ret; | ||
| 625 | } | ||
| 626 | |||
| 627 | /*** ConvexHull intersection(ConvexHull a, ConvexHull b); | ||
| 628 | * find the intersection between two convex hulls. The intersection is also a convex hull. | ||
| 629 | * (Proof: take any two points both in a and in b. Any point between them is in a by convexity, | ||
| 630 | * and in b by convexity, thus in both. Need to prove still finite bounds.) | ||
| 631 | */ | ||
| 632 | ConvexHull intersection(ConvexHull /*a*/, ConvexHull /*b*/) { | ||
| 633 | ConvexHull ret; | ||
| 634 | /* | ||
| 635 | int ai = 0, bi = 0; | ||
| 636 | int aj = a.boundary.size() - 1; | ||
| 637 | int bj = b.boundary.size() - 1; | ||
| 638 | */ | ||
| 639 | /*while (true) { | ||
| 640 | if(a[ai] | ||
| 641 | }*/ | ||
| 642 | return ret; | ||
| 643 | } | ||
| 644 | |||
| 645 | template <typename T> | ||
| 646 | T idx_to_pair(pair<T, T> p, int idx) { | ||
| 647 | return idx?p.second:p.first; | ||
| 648 | } | ||
| 649 | |||
| 650 | /*** ConvexHull merge(ConvexHull a, ConvexHull b); | ||
| 651 | * find the smallest convex hull that surrounds a and b. | ||
| 652 | */ | ||
| 653 | ConvexHull merge(ConvexHull a, ConvexHull b) { | ||
| 654 | ConvexHull ret; | ||
| 655 | |||
| 656 | std::cout << "---\n"; | ||
| 657 | std::vector<pair<int, int> > bpair = bridges(a, b); | ||
| 658 | |||
| 659 | // Given our list of bridges {(pb1, qb1), ..., (pbk, qbk)} | ||
| 660 | // we start with the highest point in p0, q0, say it is p0. | ||
| 661 | // then the merged hull is p0, ..., pb1, qb1, ..., qb2, pb2, ... | ||
| 662 | // 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. | ||
| 663 | |||
| 664 | unsigned state = (a[0][Y] < b[0][Y])?0:1; | ||
| 665 | ret.boundary.reserve(a.size() + b.size()); | ||
| 666 | ConvexHull chs[2] = {a, b}; | ||
| 667 | unsigned idx = 0; | ||
| 668 | |||
| 669 | for(unsigned k = 0; k < bpair.size(); k++) { | ||
| 670 | unsigned limit = idx_to_pair(bpair[k], state); | ||
| 671 | std::cout << bpair[k].first << " , " << bpair[k].second << "; " | ||
| 672 | << idx << ", " << limit << ", s: " | ||
| 673 | << state | ||
| 674 | << " \n"; | ||
| 675 | while(idx <= limit) { | ||
| 676 | ret.boundary.push_back(chs[state][idx++]); | ||
| 677 | } | ||
| 678 | state = 1-state; | ||
| 679 | idx = idx_to_pair(bpair[k], state); | ||
| 680 | } | ||
| 681 | while(idx < chs[state].size()) { | ||
| 682 | ret.boundary.push_back(chs[state][idx++]); | ||
| 683 | } | ||
| 684 | return ret; | ||
| 685 | } | ||
| 686 | |||
| 687 | ConvexHull graham_merge(ConvexHull a, ConvexHull b) { | ||
| 688 | ConvexHull result; | ||
| 689 | |||
| 690 | // we can avoid the find pivot step because of top_point_first | ||
| 691 | if(b.boundary[0] <= a.boundary[0]) | ||
| 692 | swap(a, b); | ||
| 693 | |||
| 694 | result.boundary = a.boundary; | ||
| 695 | result.boundary.insert(result.boundary.end(), | ||
| 696 | b.boundary.begin(), b.boundary.end()); | ||
| 697 | |||
| 698 | /** if we modified graham scan to work top to bottom as proposed in lect754.pdf we could replace the | ||
| 699 | angle sort with a simple merge sort type algorithm. furthermore, we could do the graham scan | ||
| 700 | online, avoiding a bunch of memory copies. That would probably be linear. -- njh*/ | ||
| 701 | result.angle_sort(); | ||
| 702 | result.graham_scan(); | ||
| 703 | |||
| 704 | return result; | ||
| 705 | } | ||
| 706 | |||
| 707 | ConvexHull andrew_merge(ConvexHull a, ConvexHull b) { | ||
| 708 | ConvexHull result; | ||
| 709 | |||
| 710 | // we can avoid the find pivot step because of top_point_first | ||
| 711 | if(b.boundary[0] <= a.boundary[0]) | ||
| 712 | swap(a, b); | ||
| 713 | |||
| 714 | result.boundary = a.boundary; | ||
| 715 | result.boundary.insert(result.boundary.end(), | ||
| 716 | b.boundary.begin(), b.boundary.end()); | ||
| 717 | |||
| 718 | /** if we modified graham scan to work top to bottom as proposed in lect754.pdf we could replace the | ||
| 719 | angle sort with a simple merge sort type algorithm. furthermore, we could do the graham scan | ||
| 720 | online, avoiding a bunch of memory copies. That would probably be linear. -- njh*/ | ||
| 721 | result.andrew_scan(); | ||
| 722 | |||
| 723 | return result; | ||
| 724 | } | ||
| 725 | |||
| 726 | //TODO: reinstate | ||
| 727 | /*ConvexCover::ConvexCover(Path const &sp) : path(&sp) { | ||
| 728 | cc.reserve(sp.size()); | ||
| 729 | for(Geom::Path::const_iterator it(sp.begin()), end(sp.end()); it != end; ++it) { | ||
| 730 | cc.push_back(ConvexHull((*it).begin(), (*it).end())); | ||
| 731 | } | ||
| 732 | }*/ | ||
| 733 | |||
| 734 | double ConvexHull::centroid_and_area(Geom::Point& centroid) const { | ||
| 735 | const unsigned n = boundary.size(); | ||
| 736 | if (n < 2) | ||
| 737 | return 0; | ||
| 738 | if(n < 3) { | ||
| 739 | centroid = (boundary[0] + boundary[1])/2; | ||
| 740 | return 0; | ||
| 741 | } | ||
| 742 | Geom::Point centroid_tmp(0,0); | ||
| 743 | double atmp = 0; | ||
| 744 | for (unsigned i = n-1, j = 0; j < n; i = j, j++) { | ||
| 745 | const double ai = cross(boundary[j], boundary[i]); | ||
| 746 | atmp += ai; | ||
| 747 | centroid_tmp += (boundary[j] + boundary[i])*ai; // first moment. | ||
| 748 | } | ||
| 749 | if (atmp != 0) { | ||
| 750 | centroid = centroid_tmp / (3 * atmp); | ||
| 751 | } | ||
| 752 | return atmp / 2; | ||
| 753 | } | ||
| 754 | |||
| 755 | // TODO: This can be made lg(n) using golden section/fibonacci search three starting points, say 0, | ||
| 756 | // n/2, n-1 construct a new point, say (n/2 + n)/2 throw away the furthest boundary point iterate | ||
| 757 | // until interval is a single value | ||
| 758 | Point const * ConvexHull::furthest(Point direction) const { | ||
| 759 | Point const * p = &boundary[0]; | ||
| 760 | double d = dot(*p, direction); | ||
| 761 | for(unsigned i = 1; i < boundary.size(); i++) { | ||
| 762 | double dd = dot(boundary[i], direction); | ||
| 763 | if(d < dd) { | ||
| 764 | p = &boundary[i]; | ||
| 765 | d = dd; | ||
| 766 | } | ||
| 767 | } | ||
| 768 | return p; | ||
| 769 | } | ||
| 770 | |||
| 771 | |||
| 772 | // returns (a, (b,c)), three points which define the narrowest diameter of the hull as the pair of | ||
| 773 | // lines going through b,c, and through a, parallel to b,c TODO: This can be made linear time by | ||
| 774 | // moving point tc incrementally from the previous value (it can only move in one direction). It | ||
| 775 | // is currently n*O(furthest) | ||
| 776 | double ConvexHull::narrowest_diameter(Point &a, Point &b, Point &c) { | ||
| 777 | Point tb = boundary.back(); | ||
| 778 | double d = std::numeric_limits<double>::max(); | ||
| 779 | for(unsigned i = 0; i < boundary.size(); i++) { | ||
| 780 | Point tc = boundary[i]; | ||
| 781 | Point n = -rot90(tb-tc); | ||
| 782 | Point ta = *furthest(n); | ||
| 783 | double td = dot(n, ta-tb)/dot(n,n); | ||
| 784 | if(td < d) { | ||
| 785 | a = ta; | ||
| 786 | b = tb; | ||
| 787 | c = tc; | ||
| 788 | d = td; | ||
| 789 | } | ||
| 790 | tb = tc; | ||
| 791 | } | ||
| 792 | return d; | ||
| 793 | } | ||
| 794 | #endif | ||
| 795 | |||
| 796 | }; | ||
| 797 | |||
| 798 | /* | ||
| 799 | Local Variables: | ||
| 800 | mode:c++ | ||
| 801 | c-file-style:"stroustrup" | ||
| 802 | c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) | ||
| 803 | indent-tabs-mode:nil | ||
| 804 | fill-column:99 | ||
| 805 | End: | ||
| 806 | */ | ||
| 807 | // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 : | ||
| 808 |