mpadge.tngl.sh
Distances on Directed Graphs in R
1//
2// Author: Stanislaw Adaszewski, 2019
3// C++ port from https://github.com/mapbox/concaveman (js)
4//
5// Comments from js repo added by wheeled
6//
7
8#pragma once
9
10#include <memory>
11#include <stdexcept>
12#include <list>
13#include <array>
14#include <vector>
15#include <string>
16#include <iostream>
17#include <iomanip>
18#include <stdlib.h>
19#include <limits>
20#include <set>
21#include <queue>
22#include <assert.h>
23
24#include "Rcpp.h"
25
26//#define DEBUG // uncomment to dump debug info to screen
27//#define DEBUG_2 // uncomment to dump second-level debug info to screen
28
29template<typename T, typename... Args>
30std::unique_ptr<T> make_unq_ptr(Args&&... args) {
31 return std::unique_ptr<T>(new T(std::forward<Args>(args)...));
32}
33
34
35template<class T> class compare_first {
36public:
37 bool operator()(const T &a, const T &b) {
38 return (std::get<0>(a) < std::get<0>(b));
39 }
40};
41
42
43template<class T> T orient2d(
44 const std::array<T, 2> &p1,
45 const std::array<T, 2> &p2,
46 const std::array<T, 2> &p3) {
47
48 T res = (p2[1] - p1[1]) * (p3[0] - p2[0]) -
49 (p2[0] - p1[0]) * (p3[1] - p2[1]);
50
51 return res;
52}
53
54
55// check if the edges (p1,q1) and (p2,q2) intersect
56template<class T> bool intersects(
57 const std::array<T, 2> &p1,
58 const std::array<T, 2> &q1,
59 const std::array<T, 2> &p2,
60 const std::array<T, 2> &q2) {
61
62 auto res = (p1[0] != q2[0] || p1[1] != q2[1]) &&
63 (q1[0] != p2[0] || q1[1] != p2[1]) &&
64 (orient2d(p1, q1, p2) > 0) != (orient2d(p1, q1, q2) > 0) &&
65 (orient2d(p2, q2, p1) > 0) != (orient2d(p2, q2, q1) > 0);
66
67 return res;
68}
69
70
71// square distance between 2 points
72template<class T> T getSqDist(
73 const std::array<T, 2> &p1,
74 const std::array<T, 2> &p2) {
75
76 auto dx = p1[0] - p2[0];
77 auto dy = p1[1] - p2[1];
78 return dx * dx + dy * dy;
79}
80
81
82// square distance from a point to a segment
83template<class T> T sqSegDist(
84 const std::array<T, 2> &p,
85 const std::array<T, 2> &p1,
86 const std::array<T, 2> &p2) {
87
88 auto x = p1[0];
89 auto y = p1[1];
90 auto dx = p2[0] - x;
91 auto dy = p2[1] - y;
92
93 if (dx != 0 || dy != 0) {
94 auto t = ((p[0] - x) * dx + (p[1] - y) * dy) / (dx * dx + dy * dy);
95 if (t > 1) {
96 x = p2[0];
97 y = p2[1];
98 } else if (t > 0) {
99 x += dx * t;
100 y += dy * t;
101 }
102 }
103
104 dx = p[0] - x;
105 dy = p[1] - y;
106
107 return dx * dx + dy * dy;
108}
109
110
111// segment to segment distance, ported from http://geomalgorithms.com/a07-_distance.html by Dan Sunday
112template<class T> T sqSegSegDist(T x0, T y0,
113 T x1, T y1,
114 T x2, T y2,
115 T x3, T y3) {
116 auto ux = x1 - x0;
117 auto uy = y1 - y0;
118 auto vx = x3 - x2;
119 auto vy = y3 - y2;
120 auto wx = x0 - x2;
121 auto wy = y0 - y2;
122 auto a = ux * ux + uy * uy;
123 auto b = ux * vx + uy * vy;
124 auto c = vx * vx + vy * vy;
125 auto d = ux * wx + uy * wy;
126 auto e = vx * wx + vy * wy;
127 auto D = a * c - b * b;
128
129 T sc, sN, tc, tN;
130 auto sD = D;
131 auto tD = D;
132
133 if (D == 0) {
134 sN = 0;
135 sD = 1;
136 tN = e;
137 tD = c;
138 } else {
139 sN = b * e - c * d;
140 tN = a * e - b * d;
141 if (sN < 0) {
142 sN = 0;
143 tN = e;
144 tD = c;
145 } else if (sN > sD) {
146 sN = sD;
147 tN = e + b;
148 tD = c;
149 }
150 }
151
152 if (tN < 0) {
153 tN = 0;
154 if (-d < 0) sN = 0;
155 else if (-d > a) sN = sD;
156 else {
157 sN = -d;
158 sD = a;
159 }
160 } else if (tN > tD) {
161 tN = tD;
162 if (-d + b < 0) sN = 0;
163 else if (-d + b > a) sN = sD;
164 else {
165 sN = -d + b;
166 sD = a;
167 }
168 }
169
170 sc = ((sN == 0) ? 0 : sN / sD);
171 tc = ((tN == 0) ? 0 : tN / tD);
172
173 auto cx = (1 - sc) * x0 + sc * x1;
174 auto cy = (1 - sc) * y0 + sc * y1;
175 auto cx2 = (1 - tc) * x2 + tc * x3;
176 auto cy2 = (1 - tc) * y2 + tc * y3;
177 auto dx = cx2 - cx;
178 auto dy = cy2 - cy;
179
180 return dx * dx + dy * dy;
181}
182
183
184template<class T, int DIM, int MAX_CHILDREN, class DATA> class rtree {
185public:
186 typedef rtree<T, DIM, MAX_CHILDREN, DATA> type;
187 typedef const type const_type;
188 typedef type *type_ptr;
189 typedef const type *type_const_ptr;
190 typedef std::array<T, DIM * 2> bounds_type;
191 typedef DATA data_type;
192
193 rtree():
194 m_is_leaf(false), m_data() {
195 for (auto i = 0; i < DIM; i++) {
196 m_bounds[i] = std::numeric_limits<T>::max();
197 m_bounds[i + DIM] = std::numeric_limits<T>::min();
198 }
199 }
200
201 rtree(data_type data, const bounds_type &bounds):
202 m_is_leaf(true), m_data(data), m_bounds(bounds) {
203 for (auto i = 0; i < DIM; i++)
204 if (bounds[i] > bounds[i + DIM])
205 throw std::runtime_error("Bounds minima have to be less than maxima");
206 }
207
208 void insert(data_type data, const bounds_type &bounds) {
209 if (m_is_leaf)
210 throw std::runtime_error("Cannot insert into leaves");
211
212 m_bounds = updated_bounds(bounds);
213 if (m_children.size() < MAX_CHILDREN) {
214 auto r = make_unq_ptr<type>(data, bounds);
215 m_children.push_back(std::move(r));
216 return;
217 }
218
219 std::reference_wrapper<type> best_child = *m_children.begin()->get();
220 auto best_volume = volume(best_child.get().updated_bounds(bounds));
221 for (auto it = ++m_children.begin(); it != m_children.end(); it++) {
222 auto v = volume((*it)->updated_bounds(bounds));
223 if (v < best_volume) {
224 best_volume = v;
225 best_child = *it->get();
226 }
227 }
228 if (!best_child.get().is_leaf()) {
229 best_child.get().insert(data, bounds);
230#ifdef DEBUG
231 std::cout << "best_child: " << bounds[0] << " " << bounds[1] << std::endl;
232#endif
233 return;
234 }
235
236 auto leaf = make_unq_ptr<type>(best_child.get().data(),
237 best_child.get().bounds());
238 best_child.get().m_is_leaf = false;
239 best_child.get().m_data = data_type();
240 best_child.get().m_children.push_back(std::move(leaf));
241 best_child.get().insert(data, bounds);
242 }
243
244 void intersection(const bounds_type &bounds,
245 std::vector<std::reference_wrapper<const_type>> &res) const {
246 if (!intersects(bounds))
247 return;
248 if (m_is_leaf) {
249 res.push_back(*this);
250 return;
251 }
252 for (auto &ch : m_children)
253 ch->intersection(bounds, res);
254 }
255
256 std::vector<std::reference_wrapper<const_type>> intersection(const bounds_type& bounds) const {
257 std::vector<std::reference_wrapper<const_type>> res;
258 intersection(bounds, res);
259 return res;
260 }
261
262 bool intersects(const bounds_type &bounds) const {
263 for (auto i = 0; i < DIM; i++) {
264 if (m_bounds[i] > bounds[i + DIM])
265 return false;
266 if (m_bounds[i + DIM] < bounds[i])
267 return false;
268 }
269 return true;
270 }
271
272 void erase(data_type data, const bounds_type &bounds) {
273 if (m_is_leaf)
274 throw std::runtime_error("Cannot erase from leaves");
275
276 if (!intersects(bounds))
277 return;
278
279 for (auto it = m_children.begin(); it != m_children.end(); ) {
280 if (!(*it)->m_is_leaf) {
281 (*it)->erase(data, bounds);
282 it++;
283 } else if ((*it)->m_data == data &&
284 (*it)->m_bounds == bounds) {
285 m_children.erase(it++);
286 } else
287 it++;
288 }
289 }
290
291 void print(int level = 0) {
292 // print the entire tree
293
294 for (auto it = m_children.begin(); it != m_children.end(); ) {
295 auto bounds = (*it)->m_bounds;
296 std::string pad(level, '\t');
297 if ((*it)->m_is_leaf) {
298 printf ("%s leaf %0.6f %0.6f \n", pad.c_str(), bounds[0], bounds[1]);
299 }
300 else {
301 printf ("%s branch %0.6f %0.6f %0.6f %0.6f \n", pad.c_str(), bounds[0], bounds[1], bounds[2], bounds[3]);
302 (*it)->print(level + 1);
303 }
304 it++;
305 }
306 }
307
308 bounds_type updated_bounds(const bounds_type &child_bounds) const {
309 bounds_type res;
310 for (auto i = 0; i < DIM; i++) {
311 res[i] = std::min(child_bounds[i], m_bounds[i]);
312 res[i + DIM] = std::max(child_bounds[i + DIM], m_bounds[i + DIM]);
313 }
314 return res;
315 }
316
317 static T volume(const bounds_type &bounds) {
318 T res = 1;
319 for (auto i = 0; i < DIM; i++) {
320 auto delta = bounds[i + DIM] - bounds[i];
321 res *= delta;
322 }
323 return res;
324 }
325
326 const bounds_type& bounds() const {
327 return m_bounds;
328 }
329
330 bool is_leaf() const {
331 return m_is_leaf;
332 }
333
334 data_type data() const {
335 return m_data;
336 }
337
338 const std::list<std::unique_ptr<type>>& children() const {
339 return m_children;
340 }
341
342 static std::string bounds_to_string(const bounds_type &bounds) {
343 std::string res = "( ";
344 for (auto i = 0; i < DIM * 2; i++) {
345 if (i > 0)
346 res += ", ";
347 res += std::to_string(bounds[i]);
348 }
349 res += " )";
350 return res;
351 }
352
353 void to_string(std::string &res, int tab) const {
354 std::string pad(tab, '\t');
355
356 if (m_is_leaf) {
357 res += pad + "{ data: " + std::to_string(m_data) +
358 ", bounds: " + bounds_to_string(m_bounds) +
359 " }";
360 return;
361 }
362
363 res += pad + "{ bounds: " + bounds_to_string(m_bounds) +
364 ", children: [\n";
365 auto i = 0;
366 for (auto &ch : m_children) {
367 if (i++ > 0)
368 res += "\n";
369 ch->to_string(res, tab + 1);
370 }
371 res += "\n" + pad + "]}";
372 }
373
374 std::string to_string() const {
375 std::string res;
376 to_string(res, 0);
377 return res;
378 }
379
380private:
381 bool m_is_leaf;
382 data_type m_data;
383 std::list<std::unique_ptr<type>> m_children;
384 bounds_type m_bounds;
385};
386
387
388template<class T> struct Node {
389 typedef Node<T> type;
390 typedef type *type_ptr;
391 typedef std::array<T, 2> point_type;
392
393 Node(): p(),
394 minX(), minY(), maxX(), maxY() {
395
396 }
397
398 Node(const point_type &p): Node() {
399 this->p = p;
400 }
401
402 point_type p;
403 T minX;
404 T minY;
405 T maxX;
406 T maxY;
407};
408
409
410template <class T> class CircularList;
411
412
413template<class T> class CircularElement {
414public:
415 typedef CircularElement<T> type;
416 typedef type *ptr_type;
417
418 template<class... Args>
419 CircularElement(Args&&... args) : m_data(std::forward<Args>(args)...) {
420 }
421
422 T& data() {
423 return m_data;
424 }
425
426 template<class... Args> CircularElement<T>* insert(Args&&... args) {
427 auto elem = new CircularElement<T>(std::forward<Args>(args)...);
428 elem->m_prev = this;
429 elem->m_next = m_next;
430 m_next->m_prev = elem;
431 m_next = elem;
432 return elem;
433 }
434
435 CircularElement<T>* prev() {
436 return m_prev;
437 }
438
439 CircularElement<T>* next() {
440 return m_next;
441 }
442
443private:
444 T m_data;
445 CircularElement<T> *m_prev;
446 CircularElement<T> *m_next;
447
448 friend class CircularList<T>;
449};
450
451
452template<class T> class CircularList {
453public:
454 typedef CircularElement<T> element_type;
455
456 CircularList(): m_last(nullptr) {
457
458 }
459
460 ~CircularList() {
461#ifdef DEBUG
462 std::cout << "~CircularList()" << std::endl;
463#endif
464 auto node = m_last;
465 while (true) {
466#ifdef DEBUG
467// std::cout << (i++) << std::endl;
468#endif
469 auto tmp = node;
470 node = node->m_next;
471 delete tmp;
472 if (node == m_last)
473 break;
474 }
475 }
476
477 template<class... Args> CircularElement<T>* insert(element_type *prev, Args&&... args) {
478 auto elem = new CircularElement<T>(std::forward<Args>(args)...);
479
480 if (prev == nullptr && m_last != nullptr)
481 throw std::runtime_error("Once the list is non-empty you must specify where to insert");
482
483 if (prev == nullptr) {
484 elem->m_prev = elem->m_next = elem;
485 } else {
486 elem->m_prev = prev;
487 elem->m_next = prev->m_next;
488 prev->m_next->m_prev = elem;
489 prev->m_next = elem;
490 }
491
492 m_last = elem;
493
494 return elem;
495 }
496
497
498private:
499 element_type *m_last;
500};
501
502
503// update the bounding box of a node's edge
504template<class T> void updateBBox(typename CircularElement<T>::ptr_type elem) {
505 auto &node(elem->data());
506 auto p1 = node.p;
507 auto p2 = elem->next()->data().p;
508 node.minX = std::min(p1[0], p2[0]);
509 node.minY = std::min(p1[1], p2[1]);
510 node.maxX = std::max(p1[0], p2[0]);
511 node.maxY = std::max(p1[1], p2[1]);
512}
513
514
515#ifdef DEBUG_2
516template<class T> void snapshot(
517 const std::array<T, 2> &a,
518 const std::array<T, 2> &b,
519 const std::array<T, 2> &c,
520 const std::array<T, 2> &d,
521 const double sqLen,
522 const double maxSqLen,
523 const std::array<T, 2> &trigger,
524 const bool use_trigger) {
525
526 if ( !use_trigger || trigger == b ) {
527 if ( !use_trigger )
528 printf ("Snapshot untriggered\n");
529 else
530 printf ("Snapshot trigger: %0.6f %0.6f \n", trigger[0], trigger[1]);
531 printf ("... segment a, b: %0.6f %0.6f, %0.6f %0.6f \n", a[0], a[1], b[0], b[1]);
532 printf ("... segment c, d: %0.6f %0.6f, %0.6f %0.6f \n", c[0], c[1], d[0], d[1]);
533 printf ("... sqDist a-b, b-c, c-d: %e, %e, %e", getSqDist(a, b), getSqDist(b, c), getSqDist(c, d));
534 printf ("... sqLen, maxSqLen: %e, %e", sqLen, maxSqLen);
535 }
536}
537#endif
538
539
540template<class T, int MAX_CHILDREN> std::vector<std::array<T, 2>> concaveman(
541 const std::vector<std::array<T, 2>> &points,
542 // start with a convex hull of the points
543 const std::vector<int> &hull,
544 // a relative measure of concavity; higher value means simpler hull
545 T concavity=2,
546 // when a segment goes below this length threshold, it won't be drilled down further
547 T lengthThreshold=0
548 ) {
549
550 typedef Node<T> node_type;
551 typedef std::array<T, 2> point_type;
552 typedef CircularElement<node_type> circ_elem_type;
553 typedef CircularList<node_type> circ_list_type;
554 typedef circ_elem_type *circ_elem_ptr_type;
555
556#ifdef DEBUG
557 std::cout << "concaveman()" << std::endl;
558#endif
559
560 // exit if hull includes all points already
561 if (hull.size() == points.size()) {
562 std::vector<point_type> res;
563 for (auto &i : hull) res.push_back(points[i]);
564 return res;
565 }
566
567 // index the points with an R-tree
568 rtree<T, 2, MAX_CHILDREN, point_type> tree;
569
570 for (auto &p : points)
571 tree.insert(p, { p[0], p[1], p[0], p[1] });
572
573 circ_list_type circList;
574 circ_elem_ptr_type last = nullptr;
575
576 std::list<circ_elem_ptr_type> queue;
577
578 // turn the convex hull into a linked list and populate the initial edge queue with the nodes
579 for (auto &idx : hull) {
580 auto &p = points[idx];
581 tree.erase(p, { p[0], p[1], p[0], p[1] });
582 last = circList.insert(last, p);
583 queue.push_back(last);
584 }
585
586#ifdef DEBUG_2
587 tree.print(0);
588#endif
589
590 // loops through the hull? why?
591#ifdef DEBUG
592 std::cout << "Starting hull: ";
593#endif
594 for (auto elem = last->next(); ; elem=elem->next()) {
595#ifdef DEBUG
596 std::cout << elem->data().p[0] << " " << elem->data().p[1] << std::endl;
597#endif
598 if (elem == last)
599 break;
600 }
601
602 // index the segments with an R-tree (for intersection checks)
603 rtree<T, 2, MAX_CHILDREN, circ_elem_ptr_type> segTree;
604 for (auto &elem : queue) {
605 auto &node(elem->data());
606 updateBBox<node_type>(elem);
607 segTree.insert(elem, { node.minX,
608 node.minY, node.maxX, node.maxY });
609 }
610
611 auto sqConcavity = concavity * concavity;
612 auto sqLenThreshold = lengthThreshold * lengthThreshold;
613
614 // process edges one by one
615 while (!queue.empty()) {
616 auto elem = *queue.begin();
617 queue.pop_front();
618
619 auto a = elem->prev()->data().p;
620 auto b = elem->data().p;
621 auto c = elem->next()->data().p;
622 auto d = elem->next()->next()->data().p;
623
624 // skip the edge if it's already short enough
625 auto sqLen = getSqDist(b, c);
626 if (sqLen < sqLenThreshold)
627 continue;
628
629 auto maxSqLen = sqLen / sqConcavity;
630
631#ifdef DEBUG_2
632 // dump key parameters either on every pass or when a certain point is 'b'
633 point_type trigger = { 151.1373474787800, -33.7733192376544 };
634 snapshot(a, b, c, d, sqLen, maxSqLen, trigger, true);
635#endif
636
637 // find the best connection point for the current edge to flex inward to
638 bool ok;
639 auto p = findCandidate(tree, a, b, c, d, maxSqLen, segTree, ok);
640
641 // if we found a connection and it satisfies our concavity measure
642 if (ok && std::min(getSqDist(p, b), getSqDist(p, c)) <= maxSqLen) {
643
644#ifdef DEBUG
645 printf ("Modifying hull, p: %0.6f %0.6f \n" ,p[0], p[1]);
646#endif
647
648 // connect the edge endpoints through this point and add 2 new edges to the queue
649 queue.push_back(elem);
650 queue.push_back(elem->insert(p));
651
652 // update point and segment indexes
653 auto &node = elem->data();
654 auto &next = elem->next()->data();
655
656 tree.erase(p, { p[0], p[1], p[0], p[1] });
657 segTree.erase(elem, { node.minX, node.minY, node.maxX, node.maxY });
658
659 updateBBox<node_type>(elem);
660 updateBBox<node_type>(elem->next());
661
662 segTree.insert(elem, { node.minX, node.minY, node.maxX, node.maxY });
663 segTree.insert(elem->next(), { next.minX, next.minY, next.maxX, next.maxY });
664 }
665#ifdef DEBUG
666 else
667 printf ("No point found along segment: %0.6f %0.6f, %0.6f %0.6f \n", b[0], b[1], c[0], c[1]);
668#endif
669 }
670
671 // convert the resulting hull linked list to an array of points
672 std::vector<point_type> concave;
673 for (auto elem = last->next(); ; elem = elem->next()) {
674 concave.push_back(elem->data().p);
675 if (elem == last)
676 break;
677 }
678
679 return concave;
680}
681
682
683template<class T, int MAX_CHILDREN> std::array<T, 2> findCandidate(
684 const rtree<T, 2, MAX_CHILDREN, std::array<T, 2>> &tree,
685 const std::array<T, 2> &a,
686 const std::array<T, 2> &b,
687 const std::array<T, 2> &c,
688 const std::array<T, 2> &d,
689 T maxDist,
690 const rtree<T, 2, MAX_CHILDREN, typename CircularElement<Node<T>>::ptr_type> &segTree,
691 bool &ok) {
692
693 typedef std::array<T, 2> point_type;
694 typedef CircularElement<Node<T>> circ_elem_type;
695 typedef rtree<T, 2, MAX_CHILDREN, std::array<T, 2>> tree_type;
696 typedef const tree_type const_tree_type;
697 typedef std::reference_wrapper<const_tree_type> tree_ref_type;
698 typedef std::tuple<T, tree_ref_type> tuple_type;
699
700#ifdef DEBUG
701 std::cout << "findCandidate(), maxDist: " << maxDist << std::endl;
702#endif
703
704 ok = false;
705
706 std::priority_queue<tuple_type, std::vector<tuple_type>, compare_first<tuple_type>> queue;
707 std::reference_wrapper<const_tree_type> node = tree;
708
709 // search through the point R-tree with a depth-first search using a priority queue
710 // in the order of distance to the edge (b, c)
711 while (true) {
712 for (auto &child : node.get().children()) {
713
714 auto bounds = child->bounds();
715 point_type pt = { bounds[0], bounds[1] };
716
717 auto dist = child->is_leaf() ? sqSegDist(pt, b, c) : sqSegBoxDist(b, c, *child);
718 if (dist > maxDist)
719 continue; // skip the node if it's farther than we ever need
720
721 queue.push(tuple_type(-dist, *child));
722 }
723
724 while (!queue.empty() && std::get<1>(queue.top()).get().is_leaf()) {
725 auto item = queue.top();
726 queue.pop();
727
728 auto bounds = std::get<1>(item).get().bounds();
729 point_type p = { bounds[0], bounds[1] };
730
731 // skip all points that are as close to adjacent edges (a,b) and (c,d),
732 // and points that would introduce self-intersections when connected
733 auto d0 = sqSegDist(p, a, b);
734 auto d1 = sqSegDist(p, c, d);
735
736#ifdef DEBUG_2
737 printf (" p: %0.6f %0.6f sqSegDist: %e, %e, %e \n", bounds[0], bounds[1], d0, std::get<0>(item), d1);
738#endif
739
740 if (-std::get<0>(item) < d0 && -std::get<0>(item) < d1 &&
741 noIntersections(b, p, segTree) &&
742 noIntersections(c, p, segTree)) {
743
744 ok = true;
745 return std::get<1>(item).get().data();
746 }
747
748#ifdef DEBUG_2
749 else {
750 bool cond1 = -std::get<0>(item) < d0;
751 bool cond2 = -std::get<0>(item) < d1;
752 bool cond3 = noIntersections(b, p, segTree);
753 bool cond4 = noIntersections(c, p, segTree);
754 std::cout << "Not OK: " << cond1 << " " << cond2 << " " << cond3 << " " << cond4 << std::endl;
755 }
756#endif
757 }
758
759 if (queue.empty())
760 break;
761
762 node = std::get<1>(queue.top());
763 queue.pop();
764 }
765
766 return point_type();
767}
768
769
770// square distance from a segment bounding box to the given one
771template<class T, int MAX_CHILDREN, class USER_DATA> T sqSegBoxDist(
772 const std::array<T, 2> &a,
773 const std::array<T, 2> &b,
774 const rtree<T, 2, MAX_CHILDREN, USER_DATA> &bbox) {
775
776 if (inside(a, bbox) || inside(b, bbox))
777 return 0;
778
779 auto &bounds = bbox.bounds();
780 auto minX = bounds[0];
781 auto minY = bounds[1];
782 auto maxX = bounds[2];
783 auto maxY = bounds[3];
784
785 auto d1 = sqSegSegDist(a[0], a[1], b[0], b[1], minX, minY, maxX, minY);
786 if (d1 == 0) return 0;
787
788 auto d2 = sqSegSegDist(a[0], a[1], b[0], b[1], minX, minY, minX, maxY);
789 if (d2 == 0) return 0;
790
791 auto d3 = sqSegSegDist(a[0], a[1], b[0], b[1], maxX, minY, maxX, maxY);
792 if (d3 == 0) return 0;
793
794 auto d4 = sqSegSegDist(a[0], a[1], b[0], b[1], minX, maxY, maxX, maxY);
795 if (d4 == 0) return 0;
796
797 return std::min(std::min(d1, d2), std::min(d3, d4));
798}
799
800
801template<class T, int MAX_CHILDREN, class USER_DATA> bool inside(
802 const std::array<T, 2> &a,
803 const rtree<T, 2, MAX_CHILDREN, USER_DATA> &bbox) {
804
805 auto &bounds = bbox.bounds();
806
807 auto minX = bounds[0];
808 auto minY = bounds[1];
809 auto maxX = bounds[2];
810 auto maxY = bounds[3];
811
812 auto res = (a[0] >= minX) &&
813 (a[0] <= maxX) &&
814 (a[1] >= minY) &&
815 (a[1] <= maxY);
816 return res;
817}
818
819
820// check if the edge (a,b) doesn't intersect any other edges
821template<class T, int MAX_CHILDREN> bool noIntersections(
822 const std::array<T, 2> &a,
823 const std::array<T, 2> &b,
824 const rtree<T, 2, MAX_CHILDREN, typename CircularElement<Node<T>>::ptr_type> &segTree) {
825
826 auto minX = std::min(a[0], b[0]);
827 auto minY = std::min(a[1], b[1]);
828 auto maxX = std::max(a[0], b[0]);
829 auto maxY = std::max(a[1], b[1]);
830
831 auto isect = segTree.intersection({ minX, minY, maxX, maxY });
832
833 for (decltype(segTree) &ch : isect) {
834 auto elem = ch.data();
835
836 if (intersects(elem->data().p, elem->next()->data().p, a, b))
837 return false;
838 }
839
840 return true;
841}