voronoi_builder.hpp 19 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521
  1. // Boost.Polygon library voronoi_builder.hpp header file
  2. // Copyright Andrii Sydorchuk 2010-2012.
  3. // Distributed under the Boost Software License, Version 1.0.
  4. // (See accompanying file LICENSE_1_0.txt or copy at
  5. // http://www.boost.org/LICENSE_1_0.txt)
  6. // See http://www.boost.org for updates, documentation, and revision history.
  7. #ifndef BOOST_POLYGON_VORONOI_BUILDER
  8. #define BOOST_POLYGON_VORONOI_BUILDER
  9. #include <algorithm>
  10. #include <map>
  11. #include <queue>
  12. #include <utility>
  13. #include <vector>
  14. #include "detail/voronoi_ctypes.hpp"
  15. #include "detail/voronoi_predicates.hpp"
  16. #include "detail/voronoi_structures.hpp"
  17. #include "voronoi_geometry_type.hpp"
  18. namespace boost {
  19. namespace polygon {
  20. // GENERAL INFO:
  21. // The sweepline algorithm implementation to compute Voronoi diagram of
  22. // points and non-intersecting segments (excluding endpoints).
  23. // Complexity - O(N*logN), memory usage - O(N), where N is the total number
  24. // of input geometries.
  25. //
  26. // CONTRACT:
  27. // 1) Input geometries should have integral (e.g. int32, int64) coordinate type.
  28. // 2) Input geometries should not intersect except their endpoints.
  29. //
  30. // IMPLEMENTATION DETAILS:
  31. // Each input point creates one input site. Each input segment creates three
  32. // input sites: two for its endpoints and one for the segment itself (this is
  33. // made to simplify output construction). All the site objects are constructed
  34. // and sorted at the algorithm initialization step. Priority queue is used to
  35. // dynamically hold circle events. At each step of the algorithm execution the
  36. // leftmost event is retrieved by comparing the current site event and the
  37. // topmost element from the circle event queue. STL map (red-black tree)
  38. // container was chosen to hold state of the beach line. The keys of the map
  39. // correspond to the neighboring sites that form a bisector and values map to
  40. // the corresponding Voronoi edges in the output data structure.
  41. template <typename T,
  42. typename CTT = detail::voronoi_ctype_traits<T>,
  43. typename VP = detail::voronoi_predicates<CTT> >
  44. class voronoi_builder {
  45. public:
  46. typedef typename CTT::int_type int_type;
  47. typedef typename CTT::fpt_type fpt_type;
  48. voronoi_builder() : index_(0) {}
  49. // Each point creates a single site event.
  50. std::size_t insert_point(const int_type& x, const int_type& y) {
  51. site_events_.push_back(site_event_type(x, y));
  52. site_events_.back().initial_index(index_);
  53. site_events_.back().source_category(SOURCE_CATEGORY_SINGLE_POINT);
  54. return index_++;
  55. }
  56. // Each segment creates three site events that correspond to:
  57. // 1) the start point of the segment;
  58. // 2) the end point of the segment;
  59. // 3) the segment itself defined by its start point.
  60. std::size_t insert_segment(
  61. const int_type& x1, const int_type& y1,
  62. const int_type& x2, const int_type& y2) {
  63. // Set up start point site.
  64. point_type p1(x1, y1);
  65. site_events_.push_back(site_event_type(p1));
  66. site_events_.back().initial_index(index_);
  67. site_events_.back().source_category(SOURCE_CATEGORY_SEGMENT_START_POINT);
  68. // Set up end point site.
  69. point_type p2(x2, y2);
  70. site_events_.push_back(site_event_type(p2));
  71. site_events_.back().initial_index(index_);
  72. site_events_.back().source_category(SOURCE_CATEGORY_SEGMENT_END_POINT);
  73. // Set up segment site.
  74. if (point_comparison_(p1, p2)) {
  75. site_events_.push_back(site_event_type(p1, p2));
  76. site_events_.back().source_category(SOURCE_CATEGORY_INITIAL_SEGMENT);
  77. } else {
  78. site_events_.push_back(site_event_type(p2, p1));
  79. site_events_.back().source_category(SOURCE_CATEGORY_REVERSE_SEGMENT);
  80. }
  81. site_events_.back().initial_index(index_);
  82. return index_++;
  83. }
  84. // Run sweepline algorithm and fill output data structure.
  85. template <typename OUTPUT>
  86. void construct(OUTPUT* output) {
  87. // Init structures.
  88. output->_reserve(site_events_.size());
  89. init_sites_queue();
  90. init_beach_line(output);
  91. // The algorithm stops when there are no events to process.
  92. event_comparison_predicate event_comparison;
  93. while (!circle_events_.empty() ||
  94. !(site_event_iterator_ == site_events_.end())) {
  95. if (circle_events_.empty()) {
  96. process_site_event(output);
  97. } else if (site_event_iterator_ == site_events_.end()) {
  98. process_circle_event(output);
  99. } else {
  100. if (event_comparison(*site_event_iterator_,
  101. circle_events_.top().first)) {
  102. process_site_event(output);
  103. } else {
  104. process_circle_event(output);
  105. }
  106. }
  107. while (!circle_events_.empty() &&
  108. !circle_events_.top().first.is_active()) {
  109. circle_events_.pop();
  110. }
  111. }
  112. beach_line_.clear();
  113. // Finish construction.
  114. output->_build();
  115. }
  116. void clear() {
  117. index_ = 0;
  118. site_events_.clear();
  119. }
  120. private:
  121. typedef detail::point_2d<int_type> point_type;
  122. typedef detail::site_event<int_type> site_event_type;
  123. typedef typename std::vector<site_event_type>::const_iterator
  124. site_event_iterator_type;
  125. typedef detail::circle_event<fpt_type> circle_event_type;
  126. typedef typename VP::template point_comparison_predicate<point_type>
  127. point_comparison_predicate;
  128. typedef typename VP::
  129. template event_comparison_predicate<site_event_type, circle_event_type>
  130. event_comparison_predicate;
  131. typedef typename VP::
  132. template circle_formation_predicate<site_event_type, circle_event_type>
  133. circle_formation_predicate_type;
  134. typedef void edge_type;
  135. typedef detail::beach_line_node_key<site_event_type> key_type;
  136. typedef detail::beach_line_node_data<edge_type, circle_event_type>
  137. value_type;
  138. typedef typename VP::template node_comparison_predicate<key_type>
  139. node_comparer_type;
  140. typedef std::map< key_type, value_type, node_comparer_type > beach_line_type;
  141. typedef typename beach_line_type::iterator beach_line_iterator;
  142. typedef std::pair<circle_event_type, beach_line_iterator> event_type;
  143. struct event_comparison_type {
  144. bool operator()(const event_type& lhs, const event_type& rhs) const {
  145. return predicate(rhs.first, lhs.first);
  146. }
  147. event_comparison_predicate predicate;
  148. };
  149. typedef detail::ordered_queue<event_type, event_comparison_type>
  150. circle_event_queue_type;
  151. typedef std::pair<point_type, beach_line_iterator> end_point_type;
  152. void init_sites_queue() {
  153. // Sort site events.
  154. std::sort(site_events_.begin(), site_events_.end(),
  155. event_comparison_predicate());
  156. // Remove duplicates.
  157. site_events_.erase(std::unique(
  158. site_events_.begin(), site_events_.end()), site_events_.end());
  159. // Index sites.
  160. for (std::size_t cur = 0; cur < site_events_.size(); ++cur) {
  161. site_events_[cur].sorted_index(cur);
  162. }
  163. // Init site iterator.
  164. site_event_iterator_ = site_events_.begin();
  165. }
  166. template <typename OUTPUT>
  167. void init_beach_line(OUTPUT* output) {
  168. if (site_events_.empty())
  169. return;
  170. if (site_events_.size() == 1) {
  171. // Handle single site event case.
  172. output->_process_single_site(site_events_[0]);
  173. ++site_event_iterator_;
  174. } else {
  175. int skip = 0;
  176. while (site_event_iterator_ != site_events_.end() &&
  177. VP::is_vertical(site_event_iterator_->point0(),
  178. site_events_.begin()->point0()) &&
  179. VP::is_vertical(*site_event_iterator_)) {
  180. ++site_event_iterator_;
  181. ++skip;
  182. }
  183. if (skip == 1) {
  184. // Init beach line with the first two sites.
  185. init_beach_line_default(output);
  186. } else {
  187. // Init beach line with collinear vertical sites.
  188. init_beach_line_collinear_sites(output);
  189. }
  190. }
  191. }
  192. // Init beach line with the two first sites.
  193. // The first site is always a point.
  194. template <typename OUTPUT>
  195. void init_beach_line_default(OUTPUT* output) {
  196. // Get the first and the second site event.
  197. site_event_iterator_type it_first = site_events_.begin();
  198. site_event_iterator_type it_second = site_events_.begin();
  199. ++it_second;
  200. insert_new_arc(
  201. *it_first, *it_first, *it_second, beach_line_.end(), output);
  202. // The second site was already processed. Move the iterator.
  203. ++site_event_iterator_;
  204. }
  205. // Init beach line with collinear sites.
  206. template <typename OUTPUT>
  207. void init_beach_line_collinear_sites(OUTPUT* output) {
  208. site_event_iterator_type it_first = site_events_.begin();
  209. site_event_iterator_type it_second = site_events_.begin();
  210. ++it_second;
  211. while (it_second != site_event_iterator_) {
  212. // Create a new beach line node.
  213. key_type new_node(*it_first, *it_second);
  214. // Update the output.
  215. edge_type* edge = output->_insert_new_edge(*it_first, *it_second).first;
  216. // Insert a new bisector into the beach line.
  217. beach_line_.insert(beach_line_.end(),
  218. std::pair<key_type, value_type>(new_node, value_type(edge)));
  219. // Update iterators.
  220. ++it_first;
  221. ++it_second;
  222. }
  223. }
  224. void deactivate_circle_event(value_type* value) {
  225. if (value->circle_event()) {
  226. value->circle_event()->deactivate();
  227. value->circle_event(NULL);
  228. }
  229. }
  230. template <typename OUTPUT>
  231. void process_site_event(OUTPUT* output) {
  232. // Get next site event to process.
  233. site_event_type site_event = *site_event_iterator_;
  234. // Move site iterator.
  235. site_event_iterator_type last = site_event_iterator_ + 1;
  236. // If a new site is an end point of some segment,
  237. // remove temporary nodes from the beach line data structure.
  238. if (!site_event.is_segment()) {
  239. while (!end_points_.empty() &&
  240. end_points_.top().first == site_event.point0()) {
  241. beach_line_iterator b_it = end_points_.top().second;
  242. end_points_.pop();
  243. beach_line_.erase(b_it);
  244. }
  245. } else {
  246. while (last != site_events_.end() &&
  247. last->is_segment() && last->point0() == site_event.point0())
  248. ++last;
  249. }
  250. // Find the node in the binary search tree with left arc
  251. // lying above the new site point.
  252. key_type new_key(*site_event_iterator_);
  253. beach_line_iterator right_it = beach_line_.lower_bound(new_key);
  254. for (; site_event_iterator_ != last; ++site_event_iterator_) {
  255. site_event = *site_event_iterator_;
  256. beach_line_iterator left_it = right_it;
  257. // Do further processing depending on the above node position.
  258. // For any two neighboring nodes the second site of the first node
  259. // is the same as the first site of the second node.
  260. if (right_it == beach_line_.end()) {
  261. // The above arc corresponds to the second arc of the last node.
  262. // Move the iterator to the last node.
  263. --left_it;
  264. // Get the second site of the last node
  265. const site_event_type& site_arc = left_it->first.right_site();
  266. // Insert new nodes into the beach line. Update the output.
  267. right_it = insert_new_arc(
  268. site_arc, site_arc, site_event, right_it, output);
  269. // Add a candidate circle to the circle event queue.
  270. // There could be only one new circle event formed by
  271. // a new bisector and the one on the left.
  272. activate_circle_event(left_it->first.left_site(),
  273. left_it->first.right_site(),
  274. site_event, right_it);
  275. } else if (right_it == beach_line_.begin()) {
  276. // The above arc corresponds to the first site of the first node.
  277. const site_event_type& site_arc = right_it->first.left_site();
  278. // Insert new nodes into the beach line. Update the output.
  279. left_it = insert_new_arc(
  280. site_arc, site_arc, site_event, right_it, output);
  281. // If the site event is a segment, update its direction.
  282. if (site_event.is_segment()) {
  283. site_event.inverse();
  284. }
  285. // Add a candidate circle to the circle event queue.
  286. // There could be only one new circle event formed by
  287. // a new bisector and the one on the right.
  288. activate_circle_event(site_event, right_it->first.left_site(),
  289. right_it->first.right_site(), right_it);
  290. right_it = left_it;
  291. } else {
  292. // The above arc corresponds neither to the first,
  293. // nor to the last site in the beach line.
  294. const site_event_type& site_arc2 = right_it->first.left_site();
  295. const site_event_type& site3 = right_it->first.right_site();
  296. // Remove the candidate circle from the event queue.
  297. deactivate_circle_event(&right_it->second);
  298. --left_it;
  299. const site_event_type& site_arc1 = left_it->first.right_site();
  300. const site_event_type& site1 = left_it->first.left_site();
  301. // Insert new nodes into the beach line. Update the output.
  302. beach_line_iterator new_node_it =
  303. insert_new_arc(site_arc1, site_arc2, site_event, right_it, output);
  304. // Add candidate circles to the circle event queue.
  305. // There could be up to two circle events formed by
  306. // a new bisector and the one on the left or right.
  307. activate_circle_event(site1, site_arc1, site_event, new_node_it);
  308. // If the site event is a segment, update its direction.
  309. if (site_event.is_segment()) {
  310. site_event.inverse();
  311. }
  312. activate_circle_event(site_event, site_arc2, site3, right_it);
  313. right_it = new_node_it;
  314. }
  315. }
  316. }
  317. // In general case circle event is made of the three consecutive sites
  318. // that form two bisectors in the beach line data structure.
  319. // Let circle event sites be A, B, C, two bisectors that define
  320. // circle event are (A, B), (B, C). During circle event processing
  321. // we remove (A, B), (B, C) and insert (A, C). As beach line comparison
  322. // works correctly only if one of the nodes is a new one we remove
  323. // (B, C) bisector and change (A, B) bisector to the (A, C). That's
  324. // why we use const_cast there and take all the responsibility that
  325. // map data structure keeps correct ordering.
  326. template <typename OUTPUT>
  327. void process_circle_event(OUTPUT* output) {
  328. // Get the topmost circle event.
  329. const event_type& e = circle_events_.top();
  330. const circle_event_type& circle_event = e.first;
  331. beach_line_iterator it_first = e.second;
  332. beach_line_iterator it_last = it_first;
  333. // Get the C site.
  334. site_event_type site3 = it_first->first.right_site();
  335. // Get the half-edge corresponding to the second bisector - (B, C).
  336. edge_type* bisector2 = it_first->second.edge();
  337. // Get the half-edge corresponding to the first bisector - (A, B).
  338. --it_first;
  339. edge_type* bisector1 = it_first->second.edge();
  340. // Get the A site.
  341. site_event_type site1 = it_first->first.left_site();
  342. if (!site1.is_segment() && site3.is_segment() &&
  343. site3.point1() == site1.point0()) {
  344. site3.inverse();
  345. }
  346. // Change the (A, B) bisector node to the (A, C) bisector node.
  347. const_cast<key_type&>(it_first->first).right_site(site3);
  348. // Insert the new bisector into the beach line.
  349. it_first->second.edge(output->_insert_new_edge(
  350. site1, site3, circle_event, bisector1, bisector2).first);
  351. // Remove the (B, C) bisector node from the beach line.
  352. beach_line_.erase(it_last);
  353. it_last = it_first;
  354. // Pop the topmost circle event from the event queue.
  355. circle_events_.pop();
  356. // Check new triplets formed by the neighboring arcs
  357. // to the left for potential circle events.
  358. if (it_first != beach_line_.begin()) {
  359. deactivate_circle_event(&it_first->second);
  360. --it_first;
  361. const site_event_type& site_l1 = it_first->first.left_site();
  362. activate_circle_event(site_l1, site1, site3, it_last);
  363. }
  364. // Check the new triplet formed by the neighboring arcs
  365. // to the right for potential circle events.
  366. ++it_last;
  367. if (it_last != beach_line_.end()) {
  368. deactivate_circle_event(&it_last->second);
  369. const site_event_type& site_r1 = it_last->first.right_site();
  370. activate_circle_event(site1, site3, site_r1, it_last);
  371. }
  372. }
  373. // Insert new nodes into the beach line. Update the output.
  374. template <typename OUTPUT>
  375. beach_line_iterator insert_new_arc(
  376. const site_event_type& site_arc1, const site_event_type &site_arc2,
  377. const site_event_type& site_event, beach_line_iterator position,
  378. OUTPUT* output) {
  379. // Create two new bisectors with opposite directions.
  380. key_type new_left_node(site_arc1, site_event);
  381. key_type new_right_node(site_event, site_arc2);
  382. // Set correct orientation for the first site of the second node.
  383. if (site_event.is_segment()) {
  384. new_right_node.left_site().inverse();
  385. }
  386. // Update the output.
  387. std::pair<edge_type*, edge_type*> edges =
  388. output->_insert_new_edge(site_arc2, site_event);
  389. position = beach_line_.insert(position,
  390. typename beach_line_type::value_type(
  391. new_right_node, value_type(edges.second)));
  392. if (site_event.is_segment()) {
  393. // Update the beach line with temporary bisector, that will
  394. // disappear after processing site event corresponding to the
  395. // second endpoint of the segment site.
  396. key_type new_node(site_event, site_event);
  397. new_node.right_site().inverse();
  398. position = beach_line_.insert(position,
  399. typename beach_line_type::value_type(new_node, value_type(NULL)));
  400. // Update the data structure that holds temporary bisectors.
  401. end_points_.push(std::make_pair(site_event.point1(), position));
  402. }
  403. position = beach_line_.insert(position,
  404. typename beach_line_type::value_type(
  405. new_left_node, value_type(edges.first)));
  406. return position;
  407. }
  408. // Add a new circle event to the event queue.
  409. // bisector_node corresponds to the (site2, site3) bisector.
  410. void activate_circle_event(const site_event_type& site1,
  411. const site_event_type& site2,
  412. const site_event_type& site3,
  413. beach_line_iterator bisector_node) {
  414. circle_event_type c_event;
  415. // Check if the three input sites create a circle event.
  416. if (circle_formation_predicate_(site1, site2, site3, c_event)) {
  417. // Add the new circle event to the circle events queue.
  418. // Update bisector's circle event iterator to point to the
  419. // new circle event in the circle event queue.
  420. event_type& e = circle_events_.push(
  421. std::pair<circle_event_type, beach_line_iterator>(
  422. c_event, bisector_node));
  423. bisector_node->second.circle_event(&e.first);
  424. }
  425. }
  426. private:
  427. point_comparison_predicate point_comparison_;
  428. struct end_point_comparison {
  429. bool operator() (const end_point_type& end1,
  430. const end_point_type& end2) const {
  431. return point_comparison(end2.first, end1.first);
  432. }
  433. point_comparison_predicate point_comparison;
  434. };
  435. std::vector<site_event_type> site_events_;
  436. site_event_iterator_type site_event_iterator_;
  437. std::priority_queue< end_point_type, std::vector<end_point_type>,
  438. end_point_comparison > end_points_;
  439. circle_event_queue_type circle_events_;
  440. beach_line_type beach_line_;
  441. circle_formation_predicate_type circle_formation_predicate_;
  442. std::size_t index_;
  443. // Disallow copy constructor and operator=
  444. voronoi_builder(const voronoi_builder&);
  445. void operator=(const voronoi_builder&);
  446. };
  447. typedef voronoi_builder<detail::int32> default_voronoi_builder;
  448. } // polygon
  449. } // boost
  450. #endif // BOOST_POLYGON_VORONOI_BUILDER