1 #pragma once
2
3 #ifdef DEBUG
4 #include <iostream>
5 // Example debug print for backtrace - only works on IOS
6 #include <execinfo.h>
7 #include <stdio.h>
8 //
9 // void* callstack[128];
10 // int i, frames = backtrace(callstack, 128);
11 // char** strs = backtrace_symbols(callstack, frames);
12 // for (i = 0; i < frames; ++i) {
13 // printf("%s\n", strs[i]);
14 // }
15 // free(strs);
16 #endif
17
18 #include <algorithm>
19
20 #include <mapbox/geometry/wagyu/active_bound_list.hpp>
21 #include <mapbox/geometry/wagyu/config.hpp>
22 #include <mapbox/geometry/wagyu/edge.hpp>
23 #include <mapbox/geometry/wagyu/ring.hpp>
24 #include <mapbox/geometry/wagyu/util.hpp>
25
26 namespace mapbox {
27 namespace geometry {
28 namespace wagyu {
29
30 template <typename T>
set_hole_state(bound<T> & bnd,active_bound_list<T> const & active_bounds,ring_manager<T> & rings)31 void set_hole_state(bound<T>& bnd, active_bound_list<T> const& active_bounds, ring_manager<T>& rings) {
32 auto bnd_itr = std::find(active_bounds.rbegin(), active_bounds.rend(), &bnd);
33 ++bnd_itr;
34 bound_ptr<T> bndTmp = nullptr;
35 // Find first non line ring to the left of current bound.
36 while (bnd_itr != active_bounds.rend()) {
37 if (*bnd_itr == nullptr) {
38 ++bnd_itr;
39 continue;
40 }
41 if ((*bnd_itr)->ring) {
42 if (!bndTmp) {
43 bndTmp = (*bnd_itr);
44 } else if (bndTmp->ring == (*bnd_itr)->ring) {
45 bndTmp = nullptr;
46 }
47 }
48 ++bnd_itr;
49 }
50 if (!bndTmp) {
51 bnd.ring->parent = nullptr;
52 rings.children.push_back(bnd.ring);
53 } else {
54 bnd.ring->parent = bndTmp->ring;
55 bndTmp->ring->children.push_back(bnd.ring);
56 }
57 }
58
59 template <typename T>
update_current_hp_itr(T scanline_y,ring_manager<T> & rings)60 void update_current_hp_itr(T scanline_y, ring_manager<T>& rings) {
61 while (rings.current_hp_itr->y > scanline_y) {
62 ++rings.current_hp_itr;
63 }
64 }
65
66 template <typename T>
67 struct hot_pixel_sorter {
operator ()mapbox::geometry::wagyu::hot_pixel_sorter68 inline bool operator()(mapbox::geometry::point<T> const& pt1, mapbox::geometry::point<T> const& pt2) {
69 if (pt1.y == pt2.y) {
70 return pt1.x < pt2.x;
71 } else {
72 return pt1.y > pt2.y;
73 }
74 }
75 };
76
77 template <typename T>
round_towards_min(double val)78 T round_towards_min(double val) {
79 // 0.5 rounds to 0
80 // 0.0 rounds to 0
81 // -0.5 rounds to -1
82 double half = std::floor(val) + 0.5;
83 if (values_are_equal(val, half)) {
84 return static_cast<T>(std::floor(val));
85 } else {
86 return static_cast<T>(std::llround(val));
87 }
88 }
89
90 template <typename T>
round_towards_max(double val)91 T round_towards_max(double val) {
92 // 0.5 rounds to 1
93 // 0.0 rounds to 0
94 // -0.5 rounds to 0
95 double half = std::floor(val) + 0.5;
96 if (values_are_equal(val, half)) {
97 return static_cast<T>(std::ceil(val));
98 } else {
99 return static_cast<T>(std::llround(val));
100 }
101 }
102
103 template <typename T>
get_edge_min_x(edge<T> const & edge,const T current_y)104 inline T get_edge_min_x(edge<T> const& edge, const T current_y) {
105 if (is_horizontal(edge)) {
106 if (edge.bot.x < edge.top.x) {
107 return edge.bot.x;
108 } else {
109 return edge.top.x;
110 }
111 } else if (edge.dx > 0.0) {
112 if (current_y == edge.top.y) {
113 return edge.top.x;
114 } else {
115 double lower_range_y = static_cast<double>(current_y - edge.bot.y) - 0.5;
116 double return_val = static_cast<double>(edge.bot.x) + edge.dx * lower_range_y;
117 T value = round_towards_min<T>(return_val);
118 return value;
119 }
120 } else {
121 if (current_y == edge.bot.y) {
122 return edge.bot.x;
123 } else {
124 double return_val =
125 static_cast<double>(edge.bot.x) + edge.dx * (static_cast<double>(current_y - edge.bot.y) + 0.5);
126 T value = round_towards_min<T>(return_val);
127 return value;
128 }
129 }
130 }
131
132 template <typename T>
get_edge_max_x(edge<T> const & edge,const T current_y)133 inline T get_edge_max_x(edge<T> const& edge, const T current_y) {
134 if (is_horizontal(edge)) {
135 if (edge.bot.x > edge.top.x) {
136 return edge.bot.x;
137 } else {
138 return edge.top.x;
139 }
140 } else if (edge.dx < 0.0) {
141 if (current_y == edge.top.y) {
142 return edge.top.x;
143 } else {
144 double lower_range_y = static_cast<double>(current_y - edge.bot.y) - 0.5;
145 double return_val = static_cast<double>(edge.bot.x) + edge.dx * lower_range_y;
146 T value = round_towards_max<T>(return_val);
147 return value;
148 }
149 } else {
150 if (current_y == edge.bot.y) {
151 return edge.bot.x;
152 } else {
153 double return_val =
154 static_cast<double>(edge.bot.x) + edge.dx * (static_cast<double>(current_y - edge.bot.y) + 0.5);
155 T value = round_towards_max<T>(return_val);
156 return value;
157 }
158 }
159 }
160
161 template <typename T>
hot_pixel_set_left_to_right(T y,T start_x,T end_x,bound<T> & bnd,ring_manager<T> & rings,hot_pixel_itr<T> & itr,hot_pixel_itr<T> & end,bool add_end_point)162 void hot_pixel_set_left_to_right(T y,
163 T start_x,
164 T end_x,
165 bound<T>& bnd,
166 ring_manager<T>& rings,
167 hot_pixel_itr<T>& itr,
168 hot_pixel_itr<T>& end,
169 bool add_end_point) {
170 T x_min = get_edge_min_x(*(bnd.current_edge), y);
171 x_min = std::max(x_min, start_x);
172 T x_max = get_edge_max_x(*(bnd.current_edge), y);
173 x_max = std::min(x_max, end_x);
174 for (; itr != end; ++itr) {
175 if (itr->x < x_min) {
176 continue;
177 }
178 if (itr->x > x_max) {
179 break;
180 }
181 if (!add_end_point && itr->x == end_x) {
182 continue;
183 }
184 point_ptr<T> op = bnd.ring->points;
185 bool to_front = (bnd.side == edge_left);
186 if (to_front && (*itr == *op)) {
187 continue;
188 } else if (!to_front && (*itr == *op->prev)) {
189 continue;
190 }
191 point_ptr<T> new_point = create_new_point(bnd.ring, *itr, op, rings);
192 if (to_front) {
193 bnd.ring->points = new_point;
194 }
195 }
196 }
197
198 template <typename T>
hot_pixel_set_right_to_left(T y,T start_x,T end_x,bound<T> & bnd,ring_manager<T> & rings,hot_pixel_rev_itr<T> & itr,hot_pixel_rev_itr<T> & end,bool add_end_point)199 void hot_pixel_set_right_to_left(T y,
200 T start_x,
201 T end_x,
202 bound<T>& bnd,
203 ring_manager<T>& rings,
204 hot_pixel_rev_itr<T>& itr,
205 hot_pixel_rev_itr<T>& end,
206 bool add_end_point) {
207 T x_min = get_edge_min_x(*(bnd.current_edge), y);
208 x_min = std::max(x_min, end_x);
209 T x_max = get_edge_max_x(*(bnd.current_edge), y);
210 x_max = std::min(x_max, start_x);
211 for (; itr != end; ++itr) {
212 if (itr->x > x_max) {
213 continue;
214 }
215 if (itr->x < x_min) {
216 break;
217 }
218 if (!add_end_point && itr->x == end_x) {
219 continue;
220 }
221 point_ptr<T> op = bnd.ring->points;
222 bool to_front = (bnd.side == edge_left);
223 if (to_front && (*itr == *op)) {
224 continue;
225 } else if (!to_front && (*itr == *op->prev)) {
226 continue;
227 }
228 point_ptr<T> new_point = create_new_point(bnd.ring, *itr, op, rings);
229 if (to_front) {
230 bnd.ring->points = new_point;
231 }
232 }
233 }
234
235 template <typename T>
sort_hot_pixels(ring_manager<T> & rings)236 void sort_hot_pixels(ring_manager<T>& rings) {
237 std::sort(rings.hot_pixels.begin(), rings.hot_pixels.end(), hot_pixel_sorter<T>());
238 auto last = std::unique(rings.hot_pixels.begin(), rings.hot_pixels.end());
239 rings.hot_pixels.erase(last, rings.hot_pixels.end());
240 }
241
242 template <typename T>
insert_hot_pixels_in_path(bound<T> & bnd,mapbox::geometry::point<T> const & end_pt,ring_manager<T> & rings,bool add_end_point)243 void insert_hot_pixels_in_path(bound<T>& bnd,
244 mapbox::geometry::point<T> const& end_pt,
245 ring_manager<T>& rings,
246 bool add_end_point) {
247 if (end_pt == bnd.last_point) {
248 return;
249 }
250
251 T start_y = bnd.last_point.y;
252 T start_x = bnd.last_point.x;
253 T end_y = end_pt.y;
254 T end_x = end_pt.x;
255
256 auto itr = rings.current_hp_itr;
257 while (itr->y <= start_y && itr != rings.hot_pixels.begin()) {
258 --itr;
259 }
260 if (start_x > end_x) {
261 for (; itr != rings.hot_pixels.end();) {
262 if (itr->y > start_y) {
263 ++itr;
264 continue;
265 }
266 if (itr->y < end_y) {
267 break;
268 }
269 T y = itr->y;
270 auto last_itr = hot_pixel_rev_itr<T>(itr);
271 while (itr != rings.hot_pixels.end() && itr->y == y) {
272 ++itr;
273 }
274 auto first_itr = hot_pixel_rev_itr<T>(itr);
275 bool add_end_point_itr = (y != end_pt.y || add_end_point);
276 hot_pixel_set_right_to_left(y, start_x, end_x, bnd, rings, first_itr, last_itr, add_end_point_itr);
277 }
278 } else {
279 for (; itr != rings.hot_pixels.end();) {
280 if (itr->y > start_y) {
281 ++itr;
282 continue;
283 }
284 if (itr->y < end_y) {
285 break;
286 }
287 T y = itr->y;
288 auto first_itr = itr;
289 while (itr != rings.hot_pixels.end() && itr->y == y) {
290 ++itr;
291 }
292 auto last_itr = itr;
293 bool add_end_point_itr = (y != end_pt.y || add_end_point);
294 hot_pixel_set_left_to_right(y, start_x, end_x, bnd, rings, first_itr, last_itr, add_end_point_itr);
295 }
296 }
297 bnd.last_point = end_pt;
298 }
299
300 template <typename T>
add_to_hot_pixels(mapbox::geometry::point<T> const & pt,ring_manager<T> & rings)301 void add_to_hot_pixels(mapbox::geometry::point<T> const& pt, ring_manager<T>& rings) {
302 rings.hot_pixels.push_back(pt);
303 }
304
305 template <typename T>
add_first_point(bound<T> & bnd,active_bound_list<T> & active_bounds,mapbox::geometry::point<T> const & pt,ring_manager<T> & rings)306 void add_first_point(bound<T>& bnd,
307 active_bound_list<T>& active_bounds,
308 mapbox::geometry::point<T> const& pt,
309 ring_manager<T>& rings) {
310
311 ring_ptr<T> r = create_new_ring(rings);
312 bnd.ring = r;
313 r->points = create_new_point(r, pt, rings);
314 set_hole_state(bnd, active_bounds, rings);
315 bnd.last_point = pt;
316 }
317
318 template <typename T>
add_point_to_ring(bound<T> & bnd,mapbox::geometry::point<T> const & pt,ring_manager<T> & rings)319 void add_point_to_ring(bound<T>& bnd, mapbox::geometry::point<T> const& pt, ring_manager<T>& rings) {
320 assert(bnd.ring);
321 // Handle hot pixels
322 insert_hot_pixels_in_path(bnd, pt, rings, false);
323
324 // bnd.ring->points is the 'Left-most' point & bnd.ring->points->prev is the
325 // 'Right-most'
326 point_ptr<T> op = bnd.ring->points;
327 bool to_front = (bnd.side == edge_left);
328 if (to_front && (pt == *op)) {
329 return;
330 } else if (!to_front && (pt == *op->prev)) {
331 return;
332 }
333 point_ptr<T> new_point = create_new_point(bnd.ring, pt, bnd.ring->points, rings);
334 if (to_front) {
335 bnd.ring->points = new_point;
336 }
337 }
338
339 template <typename T>
add_point(bound<T> & bnd,active_bound_list<T> & active_bounds,mapbox::geometry::point<T> const & pt,ring_manager<T> & rings)340 void add_point(bound<T>& bnd,
341 active_bound_list<T>& active_bounds,
342 mapbox::geometry::point<T> const& pt,
343 ring_manager<T>& rings) {
344 if (bnd.ring == nullptr) {
345 add_first_point(bnd, active_bounds, pt, rings);
346 } else {
347 add_point_to_ring(bnd, pt, rings);
348 }
349 }
350
351 template <typename T>
add_local_minimum_point(bound<T> & b1,bound<T> & b2,active_bound_list<T> & active_bounds,mapbox::geometry::point<T> const & pt,ring_manager<T> & rings)352 void add_local_minimum_point(bound<T>& b1,
353 bound<T>& b2,
354 active_bound_list<T>& active_bounds,
355 mapbox::geometry::point<T> const& pt,
356 ring_manager<T>& rings) {
357 if (is_horizontal(*b2.current_edge) || (b1.current_edge->dx > b2.current_edge->dx)) {
358 add_point(b1, active_bounds, pt, rings);
359 b2.last_point = pt;
360 b2.ring = b1.ring;
361 b1.side = edge_left;
362 b2.side = edge_right;
363 } else {
364 add_point(b2, active_bounds, pt, rings);
365 b1.last_point = pt;
366 b1.ring = b2.ring;
367 b1.side = edge_right;
368 b2.side = edge_left;
369 }
370 }
371
372 template <typename T>
get_dx(point<T> const & pt1,point<T> const & pt2)373 inline double get_dx(point<T> const& pt1, point<T> const& pt2) {
374 if (pt1.y == pt2.y) {
375 return std::numeric_limits<double>::infinity();
376 } else {
377 return static_cast<double>(pt2.x - pt1.x) / static_cast<double>(pt2.y - pt1.y);
378 }
379 }
380
381 template <typename T>
first_is_bottom_point(const_point_ptr<T> btmPt1,const_point_ptr<T> btmPt2)382 bool first_is_bottom_point(const_point_ptr<T> btmPt1, const_point_ptr<T> btmPt2) {
383 point_ptr<T> p = btmPt1->prev;
384 while ((*p == *btmPt1) && (p != btmPt1)) {
385 p = p->prev;
386 }
387 double dx1p = std::fabs(get_dx(*btmPt1, *p));
388
389 p = btmPt1->next;
390 while ((*p == *btmPt1) && (p != btmPt1)) {
391 p = p->next;
392 }
393 double dx1n = std::fabs(get_dx(*btmPt1, *p));
394
395 p = btmPt2->prev;
396 while ((*p == *btmPt2) && (p != btmPt2)) {
397 p = p->prev;
398 }
399 double dx2p = std::fabs(get_dx(*btmPt2, *p));
400
401 p = btmPt2->next;
402 while ((*p == *btmPt2) && (p != btmPt2)) {
403 p = p->next;
404 }
405 double dx2n = std::fabs(get_dx(*btmPt2, *p));
406
407 if (values_are_equal(std::max(dx1p, dx1n), std::max(dx2p, dx2n)) &&
408 values_are_equal(std::min(dx1p, dx1n), std::min(dx2p, dx2n))) {
409 std::size_t s = 0;
410 mapbox::geometry::box<T> bbox({ 0, 0 }, { 0, 0 });
411 return area_from_point(btmPt1, s, bbox) > 0.0; // if otherwise identical use orientation
412 } else {
413 return (greater_than_or_equal(dx1p, dx2p) && greater_than_or_equal(dx1p, dx2n)) ||
414 (greater_than_or_equal(dx1n, dx2p) && greater_than_or_equal(dx1n, dx2n));
415 }
416 }
417
418 template <typename T>
get_bottom_point(point_ptr<T> pp)419 point_ptr<T> get_bottom_point(point_ptr<T> pp) {
420 point_ptr<T> dups = nullptr;
421 point_ptr<T> p = pp->next;
422 while (p != pp) {
423 if (p->y > pp->y) {
424 pp = p;
425 dups = nullptr;
426 } else if (p->y == pp->y && p->x <= pp->x) {
427 if (p->x < pp->x) {
428 dups = nullptr;
429 pp = p;
430 } else {
431 if (p->next != pp && p->prev != pp) {
432 dups = p;
433 }
434 }
435 }
436 p = p->next;
437 }
438 if (dups) {
439 // there appears to be at least 2 vertices at bottom_point so ...
440 while (dups != p) {
441 if (!first_is_bottom_point(p, dups)) {
442 pp = dups;
443 }
444 dups = dups->next;
445 while (*dups != *pp) {
446 dups = dups->next;
447 }
448 }
449 }
450 return pp;
451 }
452
453 template <typename T>
get_lower_most_ring(ring_ptr<T> ring1,ring_ptr<T> ring2)454 ring_ptr<T> get_lower_most_ring(ring_ptr<T> ring1, ring_ptr<T> ring2) {
455 // work out which polygon fragment has the correct hole state ...
456 if (!ring1->bottom_point) {
457 ring1->bottom_point = get_bottom_point(ring1->points);
458 }
459 if (!ring2->bottom_point) {
460 ring2->bottom_point = get_bottom_point(ring2->points);
461 }
462 point_ptr<T> pt1 = ring1->bottom_point;
463 point_ptr<T> pt2 = ring2->bottom_point;
464 if (pt1->y > pt2->y) {
465 return ring1;
466 } else if (pt1->y < pt2->y) {
467 return ring2;
468 } else if (pt1->x < pt2->x) {
469 return ring1;
470 } else if (pt1->x > pt2->x) {
471 return ring2;
472 } else if (pt1->next == pt1) {
473 return ring2;
474 } else if (pt2->next == pt2) {
475 return ring1;
476 } else if (first_is_bottom_point(pt1, pt2)) {
477 return ring1;
478 } else {
479 return ring2;
480 }
481 }
482
483 template <typename T>
ring1_child_below_ring2(ring_ptr<T> ring1,ring_ptr<T> ring2)484 bool ring1_child_below_ring2(ring_ptr<T> ring1, ring_ptr<T> ring2) {
485 do {
486 ring1 = ring1->parent;
487 if (ring1 == ring2) {
488 return true;
489 }
490 } while (ring1);
491 return false;
492 }
493
494 template <typename T>
update_points_ring(ring_ptr<T> ring)495 void update_points_ring(ring_ptr<T> ring) {
496 point_ptr<T> op = ring->points;
497 do {
498 op->ring = ring;
499 op = op->prev;
500 } while (op != ring->points);
501 }
502
503 template <typename T>
append_ring(bound<T> & b1,bound<T> & b2,active_bound_list<T> & active_bounds,ring_manager<T> & manager)504 void append_ring(bound<T>& b1, bound<T>& b2, active_bound_list<T>& active_bounds, ring_manager<T>& manager) {
505 // get the start and ends of both output polygons ...
506 ring_ptr<T> outRec1 = b1.ring;
507 ring_ptr<T> outRec2 = b2.ring;
508
509 ring_ptr<T> keep_ring;
510 bound_ptr<T> keep_bound;
511 ring_ptr<T> remove_ring;
512 bound_ptr<T> remove_bound;
513 if (ring1_child_below_ring2(outRec1, outRec2)) {
514 keep_ring = outRec2;
515 keep_bound = &b2;
516 remove_ring = outRec1;
517 remove_bound = &b1;
518 } else if (ring1_child_below_ring2(outRec2, outRec1)) {
519 keep_ring = outRec1;
520 keep_bound = &b1;
521 remove_ring = outRec2;
522 remove_bound = &b2;
523 } else if (outRec1 == get_lower_most_ring(outRec1, outRec2)) {
524 keep_ring = outRec1;
525 keep_bound = &b1;
526 remove_ring = outRec2;
527 remove_bound = &b2;
528 } else {
529 keep_ring = outRec2;
530 keep_bound = &b2;
531 remove_ring = outRec1;
532 remove_bound = &b1;
533 }
534
535 // get the start and ends of both output polygons and
536 // join b2 poly onto b1 poly and delete pointers to b2 ...
537
538 point_ptr<T> p1_lft = keep_ring->points;
539 point_ptr<T> p1_rt = p1_lft->prev;
540 point_ptr<T> p2_lft = remove_ring->points;
541 point_ptr<T> p2_rt = p2_lft->prev;
542
543 // join b2 poly onto b1 poly and delete pointers to b2 ...
544 if (keep_bound->side == edge_left) {
545 if (remove_bound->side == edge_left) {
546 // z y x a b c
547 reverse_ring(p2_lft);
548 p2_lft->next = p1_lft;
549 p1_lft->prev = p2_lft;
550 p1_rt->next = p2_rt;
551 p2_rt->prev = p1_rt;
552 keep_ring->points = p2_rt;
553 } else {
554 // x y z a b c
555 p2_rt->next = p1_lft;
556 p1_lft->prev = p2_rt;
557 p2_lft->prev = p1_rt;
558 p1_rt->next = p2_lft;
559 keep_ring->points = p2_lft;
560 }
561 } else {
562 if (remove_bound->side == edge_right) {
563 // a b c z y x
564 reverse_ring(p2_lft);
565 p1_rt->next = p2_rt;
566 p2_rt->prev = p1_rt;
567 p2_lft->next = p1_lft;
568 p1_lft->prev = p2_lft;
569 } else {
570 // a b c x y z
571 p1_rt->next = p2_lft;
572 p2_lft->prev = p1_rt;
573 p1_lft->prev = p2_rt;
574 p2_rt->next = p1_lft;
575 }
576 }
577
578 keep_ring->bottom_point = nullptr;
579 bool keep_is_hole = ring_is_hole(keep_ring);
580 bool remove_is_hole = ring_is_hole(remove_ring);
581
582 remove_ring->points = nullptr;
583 remove_ring->bottom_point = nullptr;
584 if (keep_is_hole != remove_is_hole) {
585 ring1_replaces_ring2(keep_ring->parent, remove_ring, manager);
586 } else {
587 ring1_replaces_ring2(keep_ring, remove_ring, manager);
588 }
589
590 update_points_ring(keep_ring);
591
592 // nb: safe because we only get here via AddLocalMaxPoly
593 keep_bound->ring = nullptr;
594 remove_bound->ring = nullptr;
595
596 for (auto& b : active_bounds) {
597 if (b == nullptr) {
598 continue;
599 }
600 if (b->ring == remove_ring) {
601 b->ring = keep_ring;
602 b->side = keep_bound->side;
603 break; // Not sure why there is a break here but was transfered logic from angus
604 }
605 }
606 }
607
608 template <typename T>
add_local_maximum_point(bound<T> & b1,bound<T> & b2,mapbox::geometry::point<T> const & pt,ring_manager<T> & rings,active_bound_list<T> & active_bounds)609 void add_local_maximum_point(bound<T>& b1,
610 bound<T>& b2,
611 mapbox::geometry::point<T> const& pt,
612 ring_manager<T>& rings,
613 active_bound_list<T>& active_bounds) {
614 insert_hot_pixels_in_path(b2, pt, rings, false);
615 add_point(b1, active_bounds, pt, rings);
616 if (b1.ring == b2.ring) {
617 b1.ring = nullptr;
618 b2.ring = nullptr;
619 // I am not certain that order is important here?
620 } else if (b1.ring->ring_index < b2.ring->ring_index) {
621 append_ring(b1, b2, active_bounds, rings);
622 } else {
623 append_ring(b2, b1, active_bounds, rings);
624 }
625 }
626
627 enum point_in_polygon_result : std::int8_t {
628 point_on_polygon = -1,
629 point_inside_polygon = 0,
630 point_outside_polygon = 1
631 };
632
633 template <typename T>
point_in_polygon(point<T> const & pt,point_ptr<T> op)634 point_in_polygon_result point_in_polygon(point<T> const& pt, point_ptr<T> op) {
635 // returns 0 if false, +1 if true, -1 if pt ON polygon boundary
636 point_in_polygon_result result = point_outside_polygon;
637 point_ptr<T> startOp = op;
638 do {
639 if (op->next->y == pt.y) {
640 if ((op->next->x == pt.x) || (op->y == pt.y && ((op->next->x > pt.x) == (op->x < pt.x)))) {
641 return point_on_polygon;
642 }
643 }
644 if ((op->y < pt.y) != (op->next->y < pt.y)) {
645 if (op->x >= pt.x) {
646 if (op->next->x > pt.x) {
647 // Switch between point outside polygon and point inside
648 // polygon
649 if (result == point_outside_polygon) {
650 result = point_inside_polygon;
651 } else {
652 result = point_outside_polygon;
653 }
654 } else {
655 double d = static_cast<double>(op->x - pt.x) * static_cast<double>(op->next->y - pt.y) -
656 static_cast<double>(op->next->x - pt.x) * static_cast<double>(op->y - pt.y);
657 if (value_is_zero(d)) {
658 return point_on_polygon;
659 }
660 if ((d > 0) == (op->next->y > op->y)) {
661 // Switch between point outside polygon and point inside
662 // polygon
663 if (result == point_outside_polygon) {
664 result = point_inside_polygon;
665 } else {
666 result = point_outside_polygon;
667 }
668 }
669 }
670 } else {
671 if (op->next->x > pt.x) {
672 double d = static_cast<double>(op->x - pt.x) * static_cast<double>(op->next->y - pt.y) -
673 static_cast<double>(op->next->x - pt.x) * static_cast<double>(op->y - pt.y);
674 if (value_is_zero(d)) {
675 return point_on_polygon;
676 }
677 if ((d > 0) == (op->next->y > op->y)) {
678 // Switch between point outside polygon and point inside
679 // polygon
680 if (result == point_outside_polygon) {
681 result = point_inside_polygon;
682 } else {
683 result = point_outside_polygon;
684 }
685 }
686 }
687 }
688 }
689 op = op->next;
690 } while (startOp != op);
691 return result;
692 }
693
694 template <typename T>
point_in_polygon(mapbox::geometry::point<double> const & pt,point_ptr<T> op)695 point_in_polygon_result point_in_polygon(mapbox::geometry::point<double> const& pt, point_ptr<T> op) {
696 // returns 0 if false, +1 if true, -1 if pt ON polygon boundary
697 point_in_polygon_result result = point_outside_polygon;
698 point_ptr<T> startOp = op;
699 do {
700 double op_x = static_cast<double>(op->x);
701 double op_y = static_cast<double>(op->y);
702 double op_next_x = static_cast<double>(op->next->x);
703 double op_next_y = static_cast<double>(op->next->y);
704 if (values_are_equal(op_next_y, pt.y)) {
705 if (values_are_equal(op_next_x, pt.x) ||
706 (values_are_equal(op_y, pt.y) && ((op_next_x > pt.x) == (op_x < pt.x)))) {
707 return point_on_polygon;
708 }
709 }
710 if ((op_y < pt.y) != (op_next_y < pt.y)) {
711 if (greater_than_or_equal(op_x, pt.x)) {
712 if (op_next_x > pt.x) {
713 // Switch between point outside polygon and point inside
714 // polygon
715 if (result == point_outside_polygon) {
716 result = point_inside_polygon;
717 } else {
718 result = point_outside_polygon;
719 }
720 } else {
721 double d = (op_x - pt.x) * (op_next_y - pt.y) - (op_next_x - pt.x) * (op_y - pt.y);
722 if (value_is_zero(d)) {
723 return point_on_polygon;
724 }
725 if ((d > 0.0) == (op_next_y > op_y)) {
726 // Switch between point outside polygon and point inside
727 // polygon
728 if (result == point_outside_polygon) {
729 result = point_inside_polygon;
730 } else {
731 result = point_outside_polygon;
732 }
733 }
734 }
735 } else {
736 if (op_next_x > pt.x) {
737 double d = (op_x - pt.x) * (op_next_y - pt.y) - (op_next_x - pt.x) * (op_y - pt.y);
738 if (value_is_zero(d)) {
739 return point_on_polygon;
740 }
741 if ((d > 0.0) == (op_next_y > op_y)) {
742 // Switch between point outside polygon and point inside
743 // polygon
744 if (result == point_outside_polygon) {
745 result = point_inside_polygon;
746 } else {
747 result = point_outside_polygon;
748 }
749 }
750 }
751 }
752 }
753 op = op->next;
754 } while (startOp != op);
755 return result;
756 }
757
758 template <typename T>
is_convex(point_ptr<T> edge)759 bool is_convex(point_ptr<T> edge) {
760 point_ptr<T> prev = edge->prev;
761 point_ptr<T> next = edge->next;
762 T v1x = edge->x - prev->x;
763 T v1y = edge->y - prev->y;
764 T v2x = next->x - edge->x;
765 T v2y = next->y - edge->y;
766 T cross = v1x * v2y - v2x * v1y;
767 if (cross < 0 && edge->ring->area() > 0) {
768 return true;
769 } else if (cross > 0 && edge->ring->area() < 0) {
770 return true;
771 } else {
772 return false;
773 }
774 }
775
776 template <typename T>
centroid_of_points(point_ptr<T> edge)777 mapbox::geometry::point<double> centroid_of_points(point_ptr<T> edge) {
778 point_ptr<T> prev = edge->prev;
779 point_ptr<T> next = edge->next;
780 return { static_cast<double>(prev->x + edge->x + next->x) / 3.0,
781 static_cast<double>(prev->y + edge->y + next->y) / 3.0 };
782 }
783
784 template <typename T>
inside_or_outside_special(point_ptr<T> first_pt,point_ptr<T> other_poly)785 point_in_polygon_result inside_or_outside_special(point_ptr<T> first_pt, point_ptr<T> other_poly) {
786
787 // We are going to loop through all the points
788 // of the original triangle. The goal is to find a convex edge
789 // that with its next and previous forms a triangle with its centroid
790 // that is within the first ring. Then we will check the other polygon
791 // to see if it is within this polygon.
792 point_ptr<T> itr = first_pt;
793 do {
794 if (is_convex(itr)) {
795 auto pt = centroid_of_points(itr);
796 if (point_inside_polygon == point_in_polygon(pt, first_pt)) {
797 return point_in_polygon(pt, other_poly);
798 }
799 }
800 itr = itr->next;
801 } while (itr != first_pt);
802
803 throw std::runtime_error("Could not find a point within the polygon to test");
804 }
805
806 template <typename T>
box2_contains_box1(mapbox::geometry::box<T> const & box1,mapbox::geometry::box<T> const & box2)807 bool box2_contains_box1(mapbox::geometry::box<T> const& box1, mapbox::geometry::box<T> const& box2) {
808 return (box2.max.x >= box1.max.x && box2.max.y >= box1.max.y && box2.min.x <= box1.min.x &&
809 box2.min.y <= box1.min.y);
810 }
811
812 template <typename T>
poly2_contains_poly1(ring_ptr<T> ring1,ring_ptr<T> ring2)813 bool poly2_contains_poly1(ring_ptr<T> ring1, ring_ptr<T> ring2) {
814 if (!box2_contains_box1(ring1->bbox, ring2->bbox)) {
815 return false;
816 }
817 if (std::fabs(ring2->area()) < std::fabs(ring1->area())) {
818 return false;
819 }
820 point_ptr<T> outpt1 = ring1->points->next;
821 point_ptr<T> outpt2 = ring2->points->next;
822 point_ptr<T> op = outpt1;
823 do {
824 // nb: PointInPolygon returns 0 if false, +1 if true, -1 if pt on polygon
825 point_in_polygon_result res = point_in_polygon(*op, outpt2);
826 if (res != point_on_polygon) {
827 return res == point_inside_polygon;
828 }
829 op = op->next;
830 } while (op != outpt1);
831 point_in_polygon_result res = inside_or_outside_special(outpt1, outpt2);
832 return res == point_inside_polygon;
833 }
834 } // namespace wagyu
835 } // namespace geometry
836 } // namespace mapbox
837