voronoi_diagram.hpp 19 KB

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