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 | |
56 | voronoi_builder() : index_(0) {} |
57 | |
58 | // Each point creates a single site event. |
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. |
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> |
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 | |
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 | typedef struct { |
161 | 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 | } event_comparison_type; |
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 | |
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> |
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> |
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> |
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 | |
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> |
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(site1: left_it->first.left_site(), |
312 | site2: left_it->first.right_site(), |
313 | site3: site_event, bisector_node: 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(site1: site_event, site2: right_it->first.left_site(), |
331 | site3: right_it->first.right_site(), bisector_node: 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(value: &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, site2: site_arc1, site3: site_event, bisector_node: 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(site1: site_event, site2: site_arc2, site3, bisector_node: 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> |
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(value: &it_first->second); |
417 | --it_first; |
418 | const site_event_type& site_l1 = it_first->first.left_site(); |
419 | activate_circle_event(site1: site_l1, site2: site1, site3, bisector_node: 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(value: &it_last->second); |
427 | const site_event_type& site_r1 = it_last->first.right_site(); |
428 | activate_circle_event(site1, site2: site3, site3: site_r1, bisector_node: it_last); |
429 | } |
430 | } |
431 | |
432 | // Insert new nodes into the beach line. Update the output. |
433 | template <typename 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. |
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 { |
496 | 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 | |