1 // Boost.Polygon library voronoi_diagram.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_DIAGRAM 11 #define BOOST_POLYGON_VORONOI_DIAGRAM 12 13 #include <vector> 14 #include <utility> 15 16 #include "detail/voronoi_ctypes.hpp" 17 #include "detail/voronoi_structures.hpp" 18 19 #include "voronoi_geometry_type.hpp" 20 21 namespace boost { 22 namespace polygon { 23 24 // Forward declarations. 25 template <typename T> 26 class voronoi_edge; 27 28 // Represents Voronoi cell. 29 // Data members: 30 // 1) index of the source within the initial input set 31 // 2) pointer to the incident edge 32 // 3) mutable color member 33 // Cell may contain point or segment site inside. 34 template <typename T> 35 class voronoi_cell { 36 public: 37 typedef T coordinate_type; 38 typedef std::size_t color_type; 39 typedef voronoi_edge<coordinate_type> voronoi_edge_type; 40 typedef std::size_t source_index_type; 41 typedef SourceCategory source_category_type; 42 voronoi_cell(source_index_type source_index,source_category_type source_category)43 voronoi_cell(source_index_type source_index, 44 source_category_type source_category) : 45 source_index_(source_index), 46 incident_edge_(NULL), 47 color_(source_category) {} 48 49 // Returns true if the cell contains point site, false else. contains_point() const50 bool contains_point() const { 51 source_category_type source_category = this->source_category(); 52 return belongs(source_category, GEOMETRY_CATEGORY_POINT); 53 } 54 55 // Returns true if the cell contains segment site, false else. contains_segment() const56 bool contains_segment() const { 57 source_category_type source_category = this->source_category(); 58 return belongs(source_category, GEOMETRY_CATEGORY_SEGMENT); 59 } 60 source_index() const61 source_index_type source_index() const { 62 return source_index_; 63 } 64 source_category() const65 source_category_type source_category() const { 66 return static_cast<source_category_type>(color_ & SOURCE_CATEGORY_BITMASK); 67 } 68 69 // Degenerate cells don't have any incident edges. is_degenerate() const70 bool is_degenerate() const { return incident_edge_ == NULL; } 71 incident_edge()72 voronoi_edge_type* incident_edge() { return incident_edge_; } incident_edge() const73 const voronoi_edge_type* incident_edge() const { return incident_edge_; } incident_edge(voronoi_edge_type * e)74 void incident_edge(voronoi_edge_type* e) { incident_edge_ = e; } 75 color() const76 color_type color() const { return color_ >> BITS_SHIFT; } color(color_type color) const77 void color(color_type color) const { 78 color_ &= BITS_MASK; 79 color_ |= color << BITS_SHIFT; 80 } 81 82 private: 83 // 5 color bits are reserved. 84 enum Bits { 85 BITS_SHIFT = 0x5, 86 BITS_MASK = 0x1F 87 }; 88 89 source_index_type source_index_; 90 voronoi_edge_type* incident_edge_; 91 mutable color_type color_; 92 }; 93 94 // Represents Voronoi vertex. 95 // Data members: 96 // 1) vertex coordinates 97 // 2) pointer to the incident edge 98 // 3) mutable color member 99 template <typename T> 100 class voronoi_vertex { 101 public: 102 typedef T coordinate_type; 103 typedef std::size_t color_type; 104 typedef voronoi_edge<coordinate_type> voronoi_edge_type; 105 voronoi_vertex(const coordinate_type & x,const coordinate_type & y)106 voronoi_vertex(const coordinate_type& x, const coordinate_type& y) : 107 x_(x), 108 y_(y), 109 incident_edge_(NULL), 110 color_(0) {} 111 x() const112 const coordinate_type& x() const { return x_; } y() const113 const coordinate_type& y() const { return y_; } 114 is_degenerate() const115 bool is_degenerate() const { return incident_edge_ == NULL; } 116 incident_edge()117 voronoi_edge_type* incident_edge() { return incident_edge_; } incident_edge() const118 const voronoi_edge_type* incident_edge() const { return incident_edge_; } incident_edge(voronoi_edge_type * e)119 void incident_edge(voronoi_edge_type* e) { incident_edge_ = e; } 120 color() const121 color_type color() const { return color_ >> BITS_SHIFT; } color(color_type color) const122 void color(color_type color) const { 123 color_ &= BITS_MASK; 124 color_ |= color << BITS_SHIFT; 125 } 126 127 private: 128 // 5 color bits are reserved. 129 enum Bits { 130 BITS_SHIFT = 0x5, 131 BITS_MASK = 0x1F 132 }; 133 134 coordinate_type x_; 135 coordinate_type y_; 136 voronoi_edge_type* incident_edge_; 137 mutable color_type color_; 138 }; 139 140 // Half-edge data structure. Represents Voronoi edge. 141 // Data members: 142 // 1) pointer to the corresponding cell 143 // 2) pointer to the vertex that is the starting 144 // point of the half-edge 145 // 3) pointer to the twin edge 146 // 4) pointer to the CCW next edge 147 // 5) pointer to the CCW prev edge 148 // 6) mutable color member 149 template <typename T> 150 class voronoi_edge { 151 public: 152 typedef T coordinate_type; 153 typedef voronoi_cell<coordinate_type> voronoi_cell_type; 154 typedef voronoi_vertex<coordinate_type> voronoi_vertex_type; 155 typedef voronoi_edge<coordinate_type> voronoi_edge_type; 156 typedef std::size_t color_type; 157 voronoi_edge(bool is_linear,bool is_primary)158 voronoi_edge(bool is_linear, bool is_primary) : 159 cell_(NULL), 160 vertex_(NULL), 161 twin_(NULL), 162 next_(NULL), 163 prev_(NULL), 164 color_(0) { 165 if (is_linear) 166 color_ |= BIT_IS_LINEAR; 167 if (is_primary) 168 color_ |= BIT_IS_PRIMARY; 169 } 170 cell()171 voronoi_cell_type* cell() { return cell_; } cell() const172 const voronoi_cell_type* cell() const { return cell_; } cell(voronoi_cell_type * c)173 void cell(voronoi_cell_type* c) { cell_ = c; } 174 vertex0()175 voronoi_vertex_type* vertex0() { return vertex_; } vertex0() const176 const voronoi_vertex_type* vertex0() const { return vertex_; } vertex0(voronoi_vertex_type * v)177 void vertex0(voronoi_vertex_type* v) { vertex_ = v; } 178 vertex1()179 voronoi_vertex_type* vertex1() { return twin_->vertex0(); } vertex1() const180 const voronoi_vertex_type* vertex1() const { return twin_->vertex0(); } 181 twin()182 voronoi_edge_type* twin() { return twin_; } twin() const183 const voronoi_edge_type* twin() const { return twin_; } twin(voronoi_edge_type * e)184 void twin(voronoi_edge_type* e) { twin_ = e; } 185 next()186 voronoi_edge_type* next() { return next_; } next() const187 const voronoi_edge_type* next() const { return next_; } next(voronoi_edge_type * e)188 void next(voronoi_edge_type* e) { next_ = e; } 189 prev()190 voronoi_edge_type* prev() { return prev_; } prev() const191 const voronoi_edge_type* prev() const { return prev_; } prev(voronoi_edge_type * e)192 void prev(voronoi_edge_type* e) { prev_ = e; } 193 194 // Returns a pointer to the rotation next edge 195 // over the starting point of the half-edge. rot_next()196 voronoi_edge_type* rot_next() { return prev_->twin(); } rot_next() const197 const voronoi_edge_type* rot_next() const { return prev_->twin(); } 198 199 // Returns a pointer to the rotation prev edge 200 // over the starting point of the half-edge. rot_prev()201 voronoi_edge_type* rot_prev() { return twin_->next(); } rot_prev() const202 const voronoi_edge_type* rot_prev() const { return twin_->next(); } 203 204 // Returns true if the edge is finite (segment, parabolic arc). 205 // Returns false if the edge is infinite (ray, line). is_finite() const206 bool is_finite() const { return vertex0() && vertex1(); } 207 208 // Returns true if the edge is infinite (ray, line). 209 // Returns false if the edge is finite (segment, parabolic arc). is_infinite() const210 bool is_infinite() const { return !vertex0() || !vertex1(); } 211 212 // Returns true if the edge is linear (segment, ray, line). 213 // Returns false if the edge is curved (parabolic arc). is_linear() const214 bool is_linear() const { 215 return (color_ & BIT_IS_LINEAR) ? true : false; 216 } 217 218 // Returns true if the edge is curved (parabolic arc). 219 // Returns false if the edge is linear (segment, ray, line). is_curved() const220 bool is_curved() const { 221 return (color_ & BIT_IS_LINEAR) ? false : true; 222 } 223 224 // Returns false if edge goes through the endpoint of the segment. 225 // Returns true else. is_primary() const226 bool is_primary() const { 227 return (color_ & BIT_IS_PRIMARY) ? true : false; 228 } 229 230 // Returns true if edge goes through the endpoint of the segment. 231 // Returns false else. is_secondary() const232 bool is_secondary() const { 233 return (color_ & BIT_IS_PRIMARY) ? false : true; 234 } 235 color() const236 color_type color() const { return color_ >> BITS_SHIFT; } color(color_type color) const237 void color(color_type color) const { 238 color_ &= BITS_MASK; 239 color_ |= color << BITS_SHIFT; 240 } 241 242 private: 243 // 5 color bits are reserved. 244 enum Bits { 245 BIT_IS_LINEAR = 0x1, // linear is opposite to curved 246 BIT_IS_PRIMARY = 0x2, // primary is opposite to secondary 247 248 BITS_SHIFT = 0x5, 249 BITS_MASK = 0x1F 250 }; 251 252 voronoi_cell_type* cell_; 253 voronoi_vertex_type* vertex_; 254 voronoi_edge_type* twin_; 255 voronoi_edge_type* next_; 256 voronoi_edge_type* prev_; 257 mutable color_type color_; 258 }; 259 260 template <typename T> 261 struct voronoi_diagram_traits { 262 typedef T coordinate_type; 263 typedef voronoi_cell<coordinate_type> cell_type; 264 typedef voronoi_vertex<coordinate_type> vertex_type; 265 typedef voronoi_edge<coordinate_type> edge_type; 266 class vertex_equality_predicate_type { 267 public: 268 enum { ULPS = 128 }; operator ()(const vertex_type & v1,const vertex_type & v2) const269 bool operator()(const vertex_type& v1, const vertex_type& v2) const { 270 return (ulp_cmp(v1.x(), v2.x(), ULPS) == 271 detail::ulp_comparison<T>::EQUAL) && 272 (ulp_cmp(v1.y(), v2.y(), ULPS) == 273 detail::ulp_comparison<T>::EQUAL); 274 } 275 private: 276 typename detail::ulp_comparison<T> ulp_cmp; 277 }; 278 }; 279 280 // Voronoi output data structure. 281 // CCW ordering is used on the faces perimeter and around the vertices. 282 template <typename T, typename TRAITS = voronoi_diagram_traits<T> > 283 class voronoi_diagram { 284 public: 285 typedef typename TRAITS::coordinate_type coordinate_type; 286 typedef typename TRAITS::cell_type cell_type; 287 typedef typename TRAITS::vertex_type vertex_type; 288 typedef typename TRAITS::edge_type edge_type; 289 290 typedef std::vector<cell_type> cell_container_type; 291 typedef typename cell_container_type::const_iterator const_cell_iterator; 292 293 typedef std::vector<vertex_type> vertex_container_type; 294 typedef typename vertex_container_type::const_iterator const_vertex_iterator; 295 296 typedef std::vector<edge_type> edge_container_type; 297 typedef typename edge_container_type::const_iterator const_edge_iterator; 298 voronoi_diagram()299 voronoi_diagram() {} 300 clear()301 void clear() { 302 cells_.clear(); 303 vertices_.clear(); 304 edges_.clear(); 305 } 306 cells() const307 const cell_container_type& cells() const { 308 return cells_; 309 } 310 vertices() const311 const vertex_container_type& vertices() const { 312 return vertices_; 313 } 314 edges() const315 const edge_container_type& edges() const { 316 return edges_; 317 } 318 num_cells() const319 std::size_t num_cells() const { 320 return cells_.size(); 321 } 322 num_edges() const323 std::size_t num_edges() const { 324 return edges_.size(); 325 } 326 num_vertices() const327 std::size_t num_vertices() const { 328 return vertices_.size(); 329 } 330 _reserve(std::size_t num_sites)331 void _reserve(std::size_t num_sites) { 332 cells_.reserve(num_sites); 333 vertices_.reserve(num_sites << 1); 334 edges_.reserve((num_sites << 2) + (num_sites << 1)); 335 } 336 337 template <typename CT> _process_single_site(const detail::site_event<CT> & site)338 void _process_single_site(const detail::site_event<CT>& site) { 339 cells_.push_back(cell_type(site.initial_index(), site.source_category())); 340 } 341 342 // Insert a new half-edge into the output data structure. 343 // Takes as input left and right sites that form a new bisector. 344 // Returns a pair of pointers to a new half-edges. 345 template <typename CT> _insert_new_edge(const detail::site_event<CT> & site1,const detail::site_event<CT> & site2)346 std::pair<void*, void*> _insert_new_edge( 347 const detail::site_event<CT>& site1, 348 const detail::site_event<CT>& site2) { 349 // Get sites' indexes. 350 std::size_t site_index1 = site1.sorted_index(); 351 std::size_t site_index2 = site2.sorted_index(); 352 353 bool is_linear = is_linear_edge(site1, site2); 354 bool is_primary = is_primary_edge(site1, site2); 355 356 // Create a new half-edge that belongs to the first site. 357 edges_.push_back(edge_type(is_linear, is_primary)); 358 edge_type& edge1 = edges_.back(); 359 360 // Create a new half-edge that belongs to the second site. 361 edges_.push_back(edge_type(is_linear, is_primary)); 362 edge_type& edge2 = edges_.back(); 363 364 // Add the initial cell during the first edge insertion. 365 if (cells_.empty()) { 366 cells_.push_back(cell_type( 367 site1.initial_index(), site1.source_category())); 368 } 369 370 // The second site represents a new site during site event 371 // processing. Add a new cell to the cell records. 372 cells_.push_back(cell_type( 373 site2.initial_index(), site2.source_category())); 374 375 // Set up pointers to cells. 376 edge1.cell(&cells_[site_index1]); 377 edge2.cell(&cells_[site_index2]); 378 379 // Set up twin pointers. 380 edge1.twin(&edge2); 381 edge2.twin(&edge1); 382 383 // Return a pointer to the new half-edge. 384 return std::make_pair(&edge1, &edge2); 385 } 386 387 // Insert a new half-edge into the output data structure with the 388 // start at the point where two previously added half-edges intersect. 389 // Takes as input two sites that create a new bisector, circle event 390 // that corresponds to the intersection point of the two old half-edges, 391 // pointers to those half-edges. Half-edges' direction goes out of the 392 // new Voronoi vertex point. Returns a pair of pointers to a new half-edges. 393 template <typename CT1, typename CT2> _insert_new_edge(const detail::site_event<CT1> & site1,const detail::site_event<CT1> & site3,const detail::circle_event<CT2> & circle,void * data12,void * data23)394 std::pair<void*, void*> _insert_new_edge( 395 const detail::site_event<CT1>& site1, 396 const detail::site_event<CT1>& site3, 397 const detail::circle_event<CT2>& circle, 398 void* data12, void* data23) { 399 edge_type* edge12 = static_cast<edge_type*>(data12); 400 edge_type* edge23 = static_cast<edge_type*>(data23); 401 402 // Add a new Voronoi vertex. 403 vertices_.push_back(vertex_type(circle.x(), circle.y())); 404 vertex_type& new_vertex = vertices_.back(); 405 406 // Update vertex pointers of the old edges. 407 edge12->vertex0(&new_vertex); 408 edge23->vertex0(&new_vertex); 409 410 bool is_linear = is_linear_edge(site1, site3); 411 bool is_primary = is_primary_edge(site1, site3); 412 413 // Add a new half-edge. 414 edges_.push_back(edge_type(is_linear, is_primary)); 415 edge_type& new_edge1 = edges_.back(); 416 new_edge1.cell(&cells_[site1.sorted_index()]); 417 418 // Add a new half-edge. 419 edges_.push_back(edge_type(is_linear, is_primary)); 420 edge_type& new_edge2 = edges_.back(); 421 new_edge2.cell(&cells_[site3.sorted_index()]); 422 423 // Update twin pointers. 424 new_edge1.twin(&new_edge2); 425 new_edge2.twin(&new_edge1); 426 427 // Update vertex pointer. 428 new_edge2.vertex0(&new_vertex); 429 430 // Update Voronoi prev/next pointers. 431 edge12->prev(&new_edge1); 432 new_edge1.next(edge12); 433 edge12->twin()->next(edge23); 434 edge23->prev(edge12->twin()); 435 edge23->twin()->next(&new_edge2); 436 new_edge2.prev(edge23->twin()); 437 438 // Return a pointer to the new half-edge. 439 return std::make_pair(&new_edge1, &new_edge2); 440 } 441 _build()442 void _build() { 443 // Remove degenerate edges. 444 edge_iterator last_edge = edges_.begin(); 445 for (edge_iterator it = edges_.begin(); it != edges_.end(); it += 2) { 446 const vertex_type* v1 = it->vertex0(); 447 const vertex_type* v2 = it->vertex1(); 448 if (v1 && v2 && vertex_equality_predicate_(*v1, *v2)) { 449 remove_edge(&(*it)); 450 } else { 451 if (it != last_edge) { 452 edge_type* e1 = &(*last_edge = *it); 453 edge_type* e2 = &(*(last_edge + 1) = *(it + 1)); 454 e1->twin(e2); 455 e2->twin(e1); 456 if (e1->prev()) { 457 e1->prev()->next(e1); 458 e2->next()->prev(e2); 459 } 460 if (e2->prev()) { 461 e1->next()->prev(e1); 462 e2->prev()->next(e2); 463 } 464 } 465 last_edge += 2; 466 } 467 } 468 edges_.erase(last_edge, edges_.end()); 469 470 // Set up incident edge pointers for cells and vertices. 471 for (edge_iterator it = edges_.begin(); it != edges_.end(); ++it) { 472 it->cell()->incident_edge(&(*it)); 473 if (it->vertex0()) { 474 it->vertex0()->incident_edge(&(*it)); 475 } 476 } 477 478 // Remove degenerate vertices. 479 vertex_iterator last_vertex = vertices_.begin(); 480 for (vertex_iterator it = vertices_.begin(); it != vertices_.end(); ++it) { 481 if (it->incident_edge()) { 482 if (it != last_vertex) { 483 *last_vertex = *it; 484 vertex_type* v = &(*last_vertex); 485 edge_type* e = v->incident_edge(); 486 do { 487 e->vertex0(v); 488 e = e->rot_next(); 489 } while (e != v->incident_edge()); 490 } 491 ++last_vertex; 492 } 493 } 494 vertices_.erase(last_vertex, vertices_.end()); 495 496 // Set up next/prev pointers for infinite edges. 497 if (vertices_.empty()) { 498 if (!edges_.empty()) { 499 // Update prev/next pointers for the line edges. 500 edge_iterator edge_it = edges_.begin(); 501 edge_type* edge1 = &(*edge_it); 502 edge1->next(edge1); 503 edge1->prev(edge1); 504 ++edge_it; 505 edge1 = &(*edge_it); 506 ++edge_it; 507 508 while (edge_it != edges_.end()) { 509 edge_type* edge2 = &(*edge_it); 510 ++edge_it; 511 512 edge1->next(edge2); 513 edge1->prev(edge2); 514 edge2->next(edge1); 515 edge2->prev(edge1); 516 517 edge1 = &(*edge_it); 518 ++edge_it; 519 } 520 521 edge1->next(edge1); 522 edge1->prev(edge1); 523 } 524 } else { 525 // Update prev/next pointers for the ray edges. 526 for (cell_iterator cell_it = cells_.begin(); 527 cell_it != cells_.end(); ++cell_it) { 528 if (cell_it->is_degenerate()) 529 continue; 530 // Move to the previous edge while 531 // it is possible in the CW direction. 532 edge_type* left_edge = cell_it->incident_edge(); 533 while (left_edge->prev() != NULL) { 534 left_edge = left_edge->prev(); 535 // Terminate if this is not a boundary cell. 536 if (left_edge == cell_it->incident_edge()) 537 break; 538 } 539 540 if (left_edge->prev() != NULL) 541 continue; 542 543 edge_type* right_edge = cell_it->incident_edge(); 544 while (right_edge->next() != NULL) 545 right_edge = right_edge->next(); 546 left_edge->prev(right_edge); 547 right_edge->next(left_edge); 548 } 549 } 550 } 551 552 private: 553 typedef typename cell_container_type::iterator cell_iterator; 554 typedef typename vertex_container_type::iterator vertex_iterator; 555 typedef typename edge_container_type::iterator edge_iterator; 556 typedef typename TRAITS::vertex_equality_predicate_type 557 vertex_equality_predicate_type; 558 559 template <typename SEvent> is_primary_edge(const SEvent & site1,const SEvent & site2) const560 bool is_primary_edge(const SEvent& site1, const SEvent& site2) const { 561 bool flag1 = site1.is_segment(); 562 bool flag2 = site2.is_segment(); 563 if (flag1 && !flag2) { 564 return (site1.point0() != site2.point0()) && 565 (site1.point1() != site2.point0()); 566 } 567 if (!flag1 && flag2) { 568 return (site2.point0() != site1.point0()) && 569 (site2.point1() != site1.point0()); 570 } 571 return true; 572 } 573 574 template <typename SEvent> is_linear_edge(const SEvent & site1,const SEvent & site2) const575 bool is_linear_edge(const SEvent& site1, const SEvent& site2) const { 576 if (!is_primary_edge(site1, site2)) { 577 return true; 578 } 579 return !(site1.is_segment() ^ site2.is_segment()); 580 } 581 582 // Remove degenerate edge. remove_edge(edge_type * edge)583 void remove_edge(edge_type* edge) { 584 // Update the endpoints of the incident edges to the second vertex. 585 vertex_type* vertex = edge->vertex0(); 586 edge_type* updated_edge = edge->twin()->rot_next(); 587 while (updated_edge != edge->twin()) { 588 updated_edge->vertex0(vertex); 589 updated_edge = updated_edge->rot_next(); 590 } 591 592 edge_type* edge1 = edge; 593 edge_type* edge2 = edge->twin(); 594 595 edge_type* edge1_rot_prev = edge1->rot_prev(); 596 edge_type* edge1_rot_next = edge1->rot_next(); 597 598 edge_type* edge2_rot_prev = edge2->rot_prev(); 599 edge_type* edge2_rot_next = edge2->rot_next(); 600 601 // Update prev/next pointers for the incident edges. 602 edge1_rot_next->twin()->next(edge2_rot_prev); 603 edge2_rot_prev->prev(edge1_rot_next->twin()); 604 edge1_rot_prev->prev(edge2_rot_next->twin()); 605 edge2_rot_next->twin()->next(edge1_rot_prev); 606 } 607 608 cell_container_type cells_; 609 vertex_container_type vertices_; 610 edge_container_type edges_; 611 vertex_equality_predicate_type vertex_equality_predicate_; 612 613 // Disallow copy constructor and operator= 614 voronoi_diagram(const voronoi_diagram&); 615 void operator=(const voronoi_diagram&); 616 }; 617 } // polygon 618 } // boost 619 620 #endif // BOOST_POLYGON_VORONOI_DIAGRAM 621