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