1 // Copyright (c) 2013 Technical University Braunschweig (Germany).
2 // All rights reserved.
3 //
4 // This file is part of CGAL (www.cgal.org).
5 //
6 // $URL: https://github.com/CGAL/cgal/blob/v5.3/Visibility_2/include/CGAL/Rotational_sweep_visibility_2.h $
7 // $Id: Rotational_sweep_visibility_2.h 0779373 2020-03-26T13:31:46+01:00 Sébastien Loriot
8 // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
9 //
10 //
11 // Author(s): Kan Huang <huangkandiy@gmail.com>
12 //
13
14 #ifndef CGAL_ROTATIONAL_SWEEP_VISIBILITY_2_H
15 #define CGAL_ROTATIONAL_SWEEP_VISIBILITY_2_H
16
17 #include <CGAL/license/Visibility_2.h>
18
19
20 #include <CGAL/Visibility_2/visibility_utils.h>
21 #include <CGAL/Arrangement_2.h>
22 #include <CGAL/bounding_box.h>
23 #include <CGAL/assertions.h>
24 #include <CGAL/Kernel/global_functions_2.h>
25 #include <boost/unordered_map.hpp>
26 #include <iterator>
27
28
29 namespace CGAL {
30
31 template<class Arrangement_2_ , class RegularizationCategory = CGAL::Tag_true >
32 class Rotational_sweep_visibility_2 {
33 public:
34 typedef Arrangement_2_ Arrangement_2;
35 typedef typename Arrangement_2::Traits_2 Traits_2;
36 typedef typename Arrangement_2::Geometry_traits_2 Geometry_traits_2;
37 typedef typename Arrangement_2::Vertex_const_handle Vertex_const_handle;
38 typedef typename Arrangement_2::Vertex_handle Vertex_handle;
39 typedef typename Arrangement_2::Halfedge_const_handle Halfedge_const_handle;
40 typedef typename Arrangement_2::Halfedge_handle Halfedge_handle;
41 typedef typename Arrangement_2::
42 Ccb_halfedge_const_circulator Ccb_halfedge_const_circulator;
43 typedef typename Arrangement_2::Face_const_handle Face_const_handle;
44 typedef typename Arrangement_2::Face_handle Face_handle;
45
46 typedef typename Geometry_traits_2::Kernel K;
47 typedef typename Geometry_traits_2::Point_2 Point_2;
48 typedef typename Geometry_traits_2::Ray_2 Ray_2;
49 typedef typename Geometry_traits_2::Segment_2 Segment_2;
50 typedef typename Geometry_traits_2::Line_2 Line_2;
51 typedef typename Geometry_traits_2::Vector_2 Vector_2;
52 typedef typename Geometry_traits_2::Direction_2 Direction_2;
53 typedef typename Geometry_traits_2::FT Number_type;
54 typedef typename Geometry_traits_2::Object_2 Object_2;
55
56 typedef RegularizationCategory Regularization_category;
57 typedef CGAL::Tag_true Supports_general_polygon_category;
58 typedef CGAL::Tag_true Supports_simple_polygon_category;
59
60 private:
61 typedef std::vector<Point_2> Points;
62 typedef Vertex_const_handle VH;
63 typedef std::vector<VH> VHs;
64 typedef Halfedge_const_handle EH;
65 typedef std::vector<EH> EHs;
66
67 class Less_edge: public CGAL::cpp98::binary_function<EH, EH, bool> {
68 const Geometry_traits_2* geom_traits;
69 public:
Less_edge()70 Less_edge() {}
Less_edge(const Geometry_traits_2 * traits)71 Less_edge(const Geometry_traits_2* traits):geom_traits(traits) {}
operator()72 bool operator() (const EH e1, const EH e2) const {
73 if (e1 == e2)
74 return false;
75 else {
76 return &(*e1)<&(*e2);
77 }
78 // if (e1->source() == e2->source())
79 // return Visibility_2::compare_xy_2(geom_traits,
80 // e1->target()->point(), e2->target()->point()) == SMALLER;
81 // else
82 // return Visibility_2::compare_xy_2(geom_traits,
83 // e1->source()->point(), e2->source()->point()) == SMALLER;
84
85 }
86 };
87
88 class Less_vertex: public CGAL::cpp98::binary_function<VH, VH, bool> {
89 const Geometry_traits_2* geom_traits;
90 public:
Less_vertex()91 Less_vertex() {}
Less_vertex(const Geometry_traits_2 * traits)92 Less_vertex(const Geometry_traits_2* traits):geom_traits(traits) {}
operator()93 bool operator() (const VH v1, const VH v2) const {
94 if (v1 == v2)
95 return false;
96 else
97 // I know this is dirty but it speeds up by 25%. Michael
98 return &(*v1)<&(*v2);
99 // return Visibility_2::
100 // compare_xy_2(geom_traits, v1->point(), v2->point()) == SMALLER;
101 }
102 };
103
104 class Closer_edge: public CGAL::cpp98::binary_function<EH, EH, bool> {
105 const Geometry_traits_2* geom_traits;
106 Point_2 q;
107 public:
Closer_edge()108 Closer_edge() {}
Closer_edge(const Geometry_traits_2 * traits,const Point_2 & q)109 Closer_edge(const Geometry_traits_2* traits, const Point_2& q) :
110 geom_traits(traits), q(q) {}
111
vtype(const Point_2 & c,const Point_2 & p)112 int vtype(const Point_2& c, const Point_2& p) const {
113 switch(Visibility_2::orientation_2(geom_traits, q, c, p)) {
114 case COLLINEAR:
115 if (Visibility_2::less_distance_to_point_2(geom_traits, q, c, p))
116 return 0;
117 else
118 return 3;
119 case RIGHT_TURN:
120 return 1;
121 case LEFT_TURN:
122 return 2;
123 default: CGAL_assume(false);
124 }
125 return -1;
126 }
127
operator()128 bool operator() (const EH& e1, const EH& e2) const {
129 if (e1 == e2)
130 return false;
131 const Point_2& s1=e1->source()->point(),
132 t1=e1->target()->point(),
133 s2=e2->source()->point(),
134 t2=e2->target()->point();
135 if (e1->source() == e2->source()) {
136
137 int vt1 = vtype(s1, t1),
138 vt2 = vtype(s1, t2);
139 if (vt1 != vt2)
140 return vt1 > vt2;
141 else
142 return (Visibility_2::orientation_2(geom_traits, s1, t2, t1)==
143 Visibility_2::orientation_2(geom_traits, s1, t2, q));
144 }
145
146 if (e1->target() == e2->source()) {
147 // const Point_2& p1 = s1,
148 // p2 = t2,
149 // c = s2;
150 int vt1 = vtype(t1, s1),
151 vt2 = vtype(t1, t2);
152 if (vt1 != vt2)
153 return vt1 > vt2;
154 else
155 return (Visibility_2::orientation_2(geom_traits, s2, t2, s1)==
156 Visibility_2::orientation_2(geom_traits, s2, t2, q));
157 }
158
159 if (e1->source() == e2->target()) {
160 // const Point_2& p1 = t1,
161 // p2 = s2,
162 // c = s1;
163 int vt1 = vtype(s1, t1),
164 vt2 = vtype(s1, s2);
165 if (vt1 != vt2)
166 return vt1 > vt2;
167 else return (Visibility_2::orientation_2(geom_traits, s1, s2, t1)==
168 Visibility_2::orientation_2(geom_traits, s1, s2, q));
169 }
170
171 if (e1->target() == e2->target()) {
172 // const Point_2& p1 = s1,
173 // p2 = s2,
174 // c = t1;
175 int vt1 = vtype(t1, s1),
176 vt2 = vtype(t1, s2);
177 if (vt1 != vt2)
178 return vt1 > vt2;
179 else return (Visibility_2::orientation_2(geom_traits, t1, s2, s1)==
180 Visibility_2::orientation_2(geom_traits, t1, s2, q));
181 }
182
183 Orientation e1q = Visibility_2::orientation_2(geom_traits, s1, t1, q);
184 switch (e1q)
185 {
186 case COLLINEAR:
187 if (Visibility_2::collinear(geom_traits, q, s2, t2)) {
188 //q is collinear with e1 and e2.
189 return (Visibility_2::less_distance_to_point_2(geom_traits, q, s1, s2)
190 || Visibility_2::less_distance_to_point_2(geom_traits, q, t1, t2));
191 }
192 else {
193 //q is collinear with e1 not with e2.
194 if (Visibility_2::collinear(geom_traits, s2, t2, s1))
195 return (Visibility_2::orientation_2(geom_traits, s2, t2, q)
196 == Visibility_2::orientation_2(geom_traits, s2, t2, t1));
197 else
198 return (Visibility_2::orientation_2(geom_traits, s2, t2, q)
199 == Visibility_2::orientation_2(geom_traits, s2, t2, s1));
200 }
201 break;
202 case RIGHT_TURN:
203 switch (Visibility_2::orientation_2(geom_traits, s1, t1, s2)) {
204 case COLLINEAR:
205 return Visibility_2::orientation_2(geom_traits, s1, t1, t2)!=e1q;
206 case RIGHT_TURN:
207 if (Visibility_2::orientation_2(geom_traits, s1, t1, t2) == LEFT_TURN)
208 return Visibility_2::orientation_2(geom_traits, s2, t2, q)
209 == Visibility_2::orientation_2(geom_traits, s2, t2, s1);
210 else
211 return false;
212 case LEFT_TURN:
213 if(Visibility_2::orientation_2(geom_traits, s1, t1, t2) == RIGHT_TURN)
214 return Visibility_2::orientation_2(geom_traits, s2, t2, q)
215 == Visibility_2::orientation_2(geom_traits, s2, t2, s1);
216 else
217 return true;
218 default: CGAL_assume(false);
219 }
220 break;
221 case LEFT_TURN:
222 switch (Visibility_2::orientation_2(geom_traits, s1, t1, s2)) {
223 case COLLINEAR:
224 return Visibility_2::orientation_2(geom_traits, s1, t1, t2)!=e1q;
225 case LEFT_TURN:
226 if(Visibility_2::orientation_2(geom_traits, s1, t1, t2) == RIGHT_TURN)
227 return Visibility_2::orientation_2(geom_traits, s2, t2, q)
228 == Visibility_2::orientation_2(geom_traits, s2, t2, s1);
229 else
230 return false;
231 case RIGHT_TURN:
232 if (Visibility_2::orientation_2(geom_traits, s1, t1, t2) == LEFT_TURN)
233 return Visibility_2::orientation_2(geom_traits, s2, t2, q)
234 == Visibility_2::orientation_2(geom_traits, s2, t2, s1);
235 else
236 return true;
237 default: CGAL_assume(false);
238 }
239 }
240
241 CGAL_assume(false);
242 return false;
243 }
244
245 };
246
247 const Arrangement_2 *p_arr;
248 const Geometry_traits_2 *geom_traits;
249
250
251 mutable Point_2 q; //query point
252 mutable Points polygon; //visibility polygon
253
254 mutable std::map<VH, EHs, Less_vertex> incident_edges;
255
256 mutable std::map<EH, int, Less_edge> edx; //index of active edges in
257 //the heap
258
259 mutable std::set<EH, Closer_edge> active_edges; //a set of edges that
260 //intersect the current
261 //vision ray.
262
263 mutable VHs vs; //angular sorted vertices
264 mutable EHs bad_edges; //edges that pass the query point
265 mutable VH cone_end1; //an end of visibility cone
266 mutable VH cone_end2; //another end of visibility cone
267
268 mutable typename Points::size_type cone_end1_idx;
269 //index of cone_end1->point() in
270 //visibility polygon
271
272 mutable typename Points::size_type cone_end2_idx;
273 //index of cone_end2->point() in
274 //visibility polygon
275
276 mutable bool is_vertex_query;
277 mutable bool is_edge_query;
278 mutable bool is_face_query;
279 mutable bool is_big_cone; //whether the angle of
280 //visibility_cone is greater than pi.
281
282 public:
Rotational_sweep_visibility_2()283 Rotational_sweep_visibility_2(): p_arr(nullptr), geom_traits(nullptr) {}
Rotational_sweep_visibility_2(const Arrangement_2 & arr)284 Rotational_sweep_visibility_2(const Arrangement_2& arr): p_arr(&arr) {
285 geom_traits = p_arr->geometry_traits();
286 }
287
name()288 const std::string name() const { return std::string("R_visibility_2"); }
289
290 template <typename VARR>
291 typename VARR::Face_handle
compute_visibility(const Point_2 & q,const Halfedge_const_handle e,VARR & arr_out)292 compute_visibility(
293 const Point_2& q, const Halfedge_const_handle e, VARR& arr_out) const
294 {
295 arr_out.clear();
296 bad_edges.clear();
297 this->q = q;
298
299 if (Visibility_2::compare_xy_2(geom_traits, q, e->target()->point())==EQUAL)
300 {
301 is_vertex_query = true;
302 is_edge_query = false;
303 is_face_query = false;
304 cone_end1 = e->source();
305 cone_end2 = e->next()->target();
306 is_big_cone = CGAL::right_turn(cone_end1->point(), q, cone_end2->point());
307
308 typename Arrangement_2::
309 Halfedge_around_vertex_const_circulator first, curr;
310 first = curr = e->target()->incident_halfedges();
311 do {
312 if (curr->face() == e->face())
313 bad_edges.push_back(curr);
314 else if (curr->twin()->face() == e->face())
315 bad_edges.push_back(curr->twin());
316 } while (++curr != first);
317 }
318 else {
319 is_vertex_query = false;
320 is_edge_query = true;
321 is_face_query = false;
322 cone_end1 = e->source();
323 cone_end2 = e->target();
324 bad_edges.push_back(e);
325 is_big_cone = false;
326 }
327 visibility_region_impl(e->face(), q);
328
329 //decide which inside of the visibility butterfly is needed.
330 typename Points::size_type small_idx, big_idx;
331 if ( cone_end1_idx < cone_end2_idx ) {
332 small_idx = cone_end1_idx;
333 big_idx = cone_end2_idx;
334 }
335 else {
336 small_idx = cone_end2_idx;
337 big_idx = cone_end1_idx;
338 }
339 typename Points::size_type next_idx = small_idx + 1;
340 bool is_between;
341 //indicate whether the shape between small_idx and big_idx is the visibility
342 //region required.
343 if (CGAL::right_turn(cone_end1->point(), q, cone_end2->point())) {
344 is_between = false;
345 while (next_idx != big_idx) {
346 if (CGAL::left_turn(cone_end1->point(), q, polygon[next_idx]) ||
347 CGAL::left_turn(q, cone_end2->point(), polygon[next_idx]))
348 {
349 is_between = true;
350 break;
351 }
352 next_idx++;
353 }
354 }
355 else {
356 is_between = true;
357 while (next_idx != big_idx) {
358 if (CGAL::right_turn(cone_end1->point(), q, polygon[next_idx]) ||
359 CGAL::right_turn(q, cone_end2->point(), polygon[next_idx]))
360 {
361 is_between = false;
362 break;
363 }
364 next_idx++;
365 }
366 }
367
368 typename Points::iterator first = polygon.begin();
369 std::advance(first, small_idx);
370 typename Points::iterator last = polygon.begin();
371 std::advance(last, big_idx);
372
373 if (is_between) {
374 Points polygon_out(first, last+1);
375 if (is_vertex_query)
376 polygon_out.push_back(q);
377 Visibility_2::report_while_handling_needles<Rotational_sweep_visibility_2>
378 (geom_traits, q, polygon_out, arr_out);
379 }
380 else {
381 Points polygon_out(polygon.begin(), first+1);
382 if (is_vertex_query) polygon_out.push_back(q);
383 for (typename Points::size_type i = big_idx; i != polygon.size(); i++) {
384 polygon_out.push_back(polygon[i]);
385 }
386 Visibility_2::report_while_handling_needles<Rotational_sweep_visibility_2>
387 (geom_traits, q, polygon_out, arr_out);
388 }
389
390 Visibility_2::conditional_regularize(arr_out, Regularization_category());
391
392 if (arr_out.faces_begin()->is_unbounded())
393 return ++arr_out.faces_begin();
394 else
395 return arr_out.faces_begin();
396 }
397
398 template <typename VARR>
399 typename VARR::Face_handle
compute_visibility(const Point_2 & q,const Face_const_handle f,VARR & arr_out)400 compute_visibility(
401 const Point_2& q, const Face_const_handle f, VARR& arr_out) const
402 {
403 arr_out.clear();
404 this->q = q;
405 is_vertex_query = false;
406 is_edge_query = false;
407 is_face_query = true;
408
409 visibility_region_impl(f, q);
410
411 Visibility_2::report_while_handling_needles<Rotational_sweep_visibility_2>
412 (geom_traits, q, polygon, arr_out);
413
414 Visibility_2::conditional_regularize(arr_out, Regularization_category());
415
416 if (arr_out.faces_begin()->is_unbounded())
417 return ++arr_out.faces_begin();
418 else
419 return arr_out.faces_begin();
420 }
421
is_attached()422 bool is_attached() const {
423 return (p_arr != nullptr);
424 }
425
attach(const Arrangement_2 & arr)426 void attach(const Arrangement_2& arr) {
427 p_arr = &arr;
428 geom_traits = p_arr->geometry_traits();
429 }
430
detach()431 void detach() {
432 p_arr = nullptr;
433 geom_traits = nullptr;
434 }
435
arrangement_2()436 const Arrangement_2& arrangement_2() const {
437 return *p_arr;
438 }
439
440 private:
441 //get the neighbor of v along edge e
get_neighbor(const EH e,const VH v)442 VH get_neighbor(const EH e, const VH v) const {
443 if (e->source() == v)
444 return e->target();
445 else
446 return e->source();
447 }
448
449 //check whether ray(q->dp) intersects segment(p1, p2)
do_intersect_ray(const Point_2 & q,const Point_2 & dp,const Point_2 & p1,const Point_2 & p2)450 bool do_intersect_ray(const Point_2& q,
451 const Point_2& dp,
452 const Point_2& p1,
453 const Point_2& p2) const
454 {
455 return (CGAL::orientation(q, dp, p1) != CGAL::orientation(q, dp, p2) &&
456 CGAL::orientation(q, p1, dp) == CGAL::orientation(q, p1, p2));
457 }
458
459 //arrange vertices that on a same vision ray in a 'funnel' order
funnel(typename VHs::size_type i,typename VHs::size_type j)460 void funnel(typename VHs::size_type i, typename VHs::size_type j) const {
461 VHs right, left;
462 //whether the edges incident to a vertex block the left side and right side
463 //of current vision ray.
464 bool block_left(false), block_right(false);
465 VH former = vs[i], nb;
466 for (typename VHs::size_type l=i; l<j; l++) {
467 bool left_v(false), right_v(false), has_predecessor(false);
468 EHs& edges = incident_edges[vs[l]];
469 for (typename EHs::size_type k=0; k<edges.size(); k++) {
470 nb = get_neighbor(edges[k], vs[l]);
471 if ( nb == former ) {
472 has_predecessor = true;
473 break;
474 }
475 if (CGAL::left_turn(q, vs[l]->point(), nb->point()))
476 left_v = true;
477 else
478 right_v = CGAL::right_turn(q, vs[l]->point(), nb->point());
479 }
480 if (has_predecessor) {
481 //if the current vertex connects to the vertex before by an edge,
482 //the vertex before can help it to block.
483 block_left = block_left || left_v;
484 block_right = block_right || right_v;
485 }
486 else {
487 block_left = left_v;
488 block_right = right_v;
489 }
490 if (block_left && block_right) {
491 //when both sides are blocked,
492 //there is no need to change the vertex after.
493 right.push_back(vs[l]);
494 break;
495 }
496 else {
497 if (block_left)
498 left.push_back(vs[l]);
499 else
500 right.push_back(vs[l]);
501 }
502 former = vs[l];
503 }
504 for (typename VHs::size_type l=0; l!=right.size(); l++)
505 vs[i+l] = right[l];
506 for (typename VHs::size_type l=0; l!=left.size(); l++)
507 vs[i+l+right.size()] = left[left.size()-1-l];
508 }
509
510
511
visibility_region_impl(const Face_const_handle f,const Point_2 & q)512 void visibility_region_impl(const Face_const_handle f, const Point_2& q) const
513 {
514 vs.clear();
515 polygon.clear();
516 active_edges = std::set<EH, Closer_edge>(Closer_edge(geom_traits, q));
517 incident_edges = std::map<VH, EHs, Less_vertex>(Less_vertex(geom_traits));
518 edx = std::map<EH, int, Less_edge>(Less_edge(geom_traits));
519
520 EHs relevant_edges; //edges that can affect the visibility of query point.
521 Arrangement_2 bbox;
522 if (is_face_query)
523 input_face(f);
524 else
525 input_face(f, relevant_edges, bbox);
526 //the following code is the initiation of vision ray.
527 //the direction of the initial ray is between the direction from q to last
528 //vertex in vs and positive x-axis. By choosing this direction, we make sure
529 //that all plane is swept and there is not needle at the beginning of
530 //sweeping.
531 Vector_2 dir;
532 if (Direction_2(-1, 0) < Direction_2(Vector_2(q, vs.back()->point())))
533 dir = Vector_2(1, 0) + Vector_2(q, vs.back()->point());
534 else
535 dir = Vector_2(0, -1);
536 Point_2 dp = q + dir;
537
538 //initiation of active_edges. for face queries,
539 //all edges on the boundary can affect visibility.
540 //for non-face queries, only relevant_edges has to be considered.
541 if (is_face_query) {
542 Ccb_halfedge_const_circulator curr = f->outer_ccb();
543 Ccb_halfedge_const_circulator circ = curr;
544 do {
545 if (do_intersect_ray(
546 q, dp, curr->target()->point(), curr->source()->point()))
547 active_edges.insert(curr);
548
549 } while (++curr != circ);
550
551 typename Arrangement_2::Hole_const_iterator hi;
552 for (hi = f->holes_begin(); hi != f->holes_end(); ++hi) {
553 Ccb_halfedge_const_circulator curr = circ = *hi;
554 do {
555 if (do_intersect_ray(
556 q, dp, curr->target()->point(), curr->source()->point()))
557 active_edges.insert(curr);
558 } while (++curr != circ);
559 }
560 }
561 else {
562 for (typename EHs::size_type i=0; i!=relevant_edges.size(); i++)
563 if (do_intersect_ray(q, dp, relevant_edges[i]->source()->point(),
564 relevant_edges[i]->target()->point()))
565 active_edges.insert(relevant_edges[i]);
566 }
567
568 //angular sweep begins
569 // std::cout<<active_edges.size()<<std::endl;
570 for (typename VHs::size_type i=0; i!=vs.size(); i++) {
571 VH vh = vs[i];
572 EH closest_e = *active_edges.begin();
573 EHs& edges = incident_edges[vh];
574 EHs insert_ehs, remove_ehs;
575 for (typename EHs::size_type j=0; j!=edges.size(); j++) {
576 EH& e = edges[j];
577 if (active_edges.find(e) == active_edges.end())
578 insert_ehs.push_back(e);
579 else
580 remove_ehs.push_back(e);
581 }
582 typename EHs::size_type insert_cnt = insert_ehs.size();
583 typename EHs::size_type remove_cnt = remove_ehs.size();
584 if (insert_cnt == 1 && remove_cnt == 1) {
585 const EH& ctemp_eh = *active_edges.find(remove_ehs.front());
586 EH& temp_eh = const_cast<EH&>(ctemp_eh);
587 temp_eh = insert_ehs.front();
588 }
589 else {
590 for (typename EHs::size_type j=0; j!=remove_cnt; j++)
591 active_edges.erase(remove_ehs[j]);
592 for (typename EHs::size_type j=0; j!=insert_cnt; j++)
593 active_edges.insert(insert_ehs[j]);
594 }
595
596 if (closest_e != *active_edges.begin()) {
597 //when the closest edge changed
598 if (is_face_query) {
599 if (remove_cnt > 0 && insert_cnt > 0) {
600 //some edges are added and some are deleted,
601 //which means the vertex swept is part of visibility polygon.
602 update_visibility(vh->point());
603 }
604 if (remove_cnt == 0 && insert_cnt > 0) {
605 //only add some edges, means the view ray is blocked by new edges.
606 //therefore first add the intersection of view ray and
607 //former closet edge, then add the vertice swept.
608 update_visibility(ray_seg_intersection(q,
609 vh->point(),
610 closest_e->target()->point(),
611 closest_e->source()->point())
612 );
613 update_visibility(vh->point());
614 }
615 if (remove_cnt > 0 && insert_cnt == 0) {
616 //only delete some edges, means some block is moved and the view ray
617 //can reach the segments after the block.
618 update_visibility(vh->point());
619 update_visibility(
620 ray_seg_intersection(q,
621 vh->point(),
622 (*active_edges.begin())->target()->point(),
623 (*active_edges.begin())->source()->point()
624 )
625 );
626 }
627 }
628 else {
629 //extra work here for edge/vertex query is the index of cone_end1 and
630 //cone_end2 will be recorded.
631 if (remove_cnt > 0 && insert_cnt > 0) {
632 //some edges are added and some are deleted,
633 //which means the vertice swept is part of visibility polygon.
634 if (update_visibility(vh->point())) {
635 if (vh == cone_end1)
636 cone_end1_idx = polygon.size()-1;
637 else if (vh == cone_end2)
638 cone_end2_idx = polygon.size()-1;
639 }
640 }
641 if (remove_cnt == 0 && insert_cnt > 0) {
642 //only add some edges, means the view ray is blocked by new edges.
643 //therefore first add the intersection of view ray and former closet
644 //edge, then add the vertice swept.
645 update_visibility(ray_seg_intersection(q,
646 vh->point(),
647 closest_e->target()->point(),
648 closest_e->source()->point())
649 );
650 if (update_visibility(vh->point())) {
651 if (vh == cone_end1)
652 cone_end1_idx = polygon.size()-1;
653 else if (vh == cone_end2)
654 cone_end2_idx = polygon.size()-1;
655 }
656 }
657 if (remove_cnt > 0 && insert_cnt == 0) {
658 //only delete some edges, means some block is removed and the vision
659 //ray can reach the segments after the block.
660 if (update_visibility(vh->point())) {
661 if (vh == cone_end1)
662 cone_end1_idx = polygon.size()-1;
663 else if (vh == cone_end2)
664 cone_end2_idx = polygon.size()-1;
665 }
666 update_visibility(
667 ray_seg_intersection(q,
668 vh->point(),
669 (*active_edges.begin())->target()->point(),
670 (*active_edges.begin())->source()->point())
671 );
672 }
673 }
674 }
675 }
676 }
677
print_edge(const EH e)678 void print_edge(const EH e) const {
679 std::cout << e->source()->point() <<"->"<< e->target()->point() <<std::endl;
680 }
681
682 //compute the intersection of ray(q->dp) and segment(s, t)
683 //if they are collinear then return the endpoint which is closer to q.
684
ray_seg_intersection(const Point_2 & q,const Point_2 & dp,const Point_2 & s,const Point_2 & t)685 Point_2 ray_seg_intersection(
686 const Point_2& q, const Point_2& dp, // the ray
687 const Point_2& s, const Point_2& t) // the segment
688 const
689 {
690 if (CGAL::collinear(q, dp, s)) {
691 if (CGAL::collinear(q, dp, t)) {
692 if (CGAL::compare_distance_to_point(q, s, t)==CGAL::SMALLER)
693 return s;
694 else
695 return t;
696 }
697 else
698 return s;
699 }
700 Ray_2 ray(q,dp);
701 Segment_2 seg(s,t);
702 CGAL::Object result = CGAL::intersection(ray, seg);
703 return *(CGAL::object_cast<Point_2>(&result));
704 }
705
706 //check if p has been discovered before, if not update the visibility polygon
update_visibility(const Point_2 & p)707 bool update_visibility(const Point_2& p) const {
708 if (polygon.empty()) {
709 polygon.push_back(p);
710 return true;
711 }
712 else if (Visibility_2::compare_xy_2(geom_traits, polygon.back(), p)
713 != EQUAL)
714 {
715 polygon.push_back(p);
716 return true;
717 }
718 return false;
719 }
720
721 //functor to decide which vertex is swept earlier by the rotational sweeping
722 //ray
723 class Is_swept_earlier:public CGAL::cpp98::binary_function<VH, VH, bool> {
724 const Point_2& q;
725 const Geometry_traits_2* geom_traits;
726 public:
Is_swept_earlier(const Point_2 & q,const Geometry_traits_2 * traits)727 Is_swept_earlier(const Point_2& q, const Geometry_traits_2* traits) :
728 q(q), geom_traits(traits) {}
729
operator()730 bool operator() (const VH v1, const VH v2) const {
731 const Point_2& p1 = v1->point();
732 const Point_2& p2 = v2->point();
733 int qua1 = quadrant(q, p1);
734 int qua2 = quadrant(q, p2);
735 if (qua1 < qua2)
736 return true;
737 if (qua1 > qua2)
738 return false;
739 if (collinear(q, p1, p2))
740 return (CGAL::compare_distance_to_point(q, p1, p2) == CGAL::SMALLER);
741 else
742 return CGAL::right_turn(p1, q, p2);
743 }
744
745 //return the quadrant of p with respect to o.
quadrant(const Point_2 & o,const Point_2 & p)746 int quadrant(const Point_2& o, const Point_2& p) const {
747 typename Geometry_traits_2::Compare_x_2 compare_x =
748 geom_traits->compare_x_2_object();
749
750 typename Geometry_traits_2::Compare_y_2 compare_y =
751 geom_traits->compare_y_2_object();
752
753 Comparison_result dx = compare_x(p, o);
754 Comparison_result dy = compare_y(p, o);
755 if (dx==LARGER && dy!=SMALLER)
756 return 1;
757 if (dx!=LARGER && dy==LARGER)
758 return 2;
759 if (dx==SMALLER && dy!=LARGER)
760 return 3;
761 if (dx!=SMALLER && dy==SMALLER)
762 return 4;
763 return 0;
764 }
765 };
766
767 //when the query point is in face, every edge is good.
input_neighbor_f(const Halfedge_const_handle e)768 void input_neighbor_f( const Halfedge_const_handle e) const {
769 VH v = e->target();
770 if (!incident_edges.count(v))
771 vs.push_back(v);
772 incident_edges[v].push_back(e);
773 incident_edges[v].push_back(e->next());
774 }
775
776 //check if p is in the visibility cone
is_in_cone(const Point_2 & p)777 bool is_in_cone(const Point_2& p) const{
778 if (is_big_cone)
779 return (!CGAL::right_turn(cone_end1->point(), q, p)) ||
780 (!CGAL::left_turn(cone_end2->point(), q, p));
781 else
782 return (!CGAL::right_turn(cone_end1->point(), q, p)) &&
783 (!CGAL::left_turn(cone_end2->point(), q, p));
784 }
785
786 //for vertex and edge query: the visibility is limited in a cone.
input_edge(const Halfedge_const_handle e,EHs & good_edges)787 void input_edge(const Halfedge_const_handle e,
788 EHs& good_edges) const {
789 for (typename EHs::size_type i = 0; i < bad_edges.size(); i++)
790 if (e == bad_edges[i])
791 return;
792 VH v1 = e->target();
793 VH v2 = e->source();
794 //an edge will affect visibility only if it has an endpoint in the
795 //visibility cone or it crosses the boundary of the cone.
796 if (is_in_cone(v1->point()) || is_in_cone(v2->point()) ||
797 do_intersect_ray(q, cone_end1->point(), v1->point(), v2->point()))
798 {
799 good_edges.push_back(e);
800 if (!incident_edges.count(v1))
801 vs.push_back(v1);
802 incident_edges[v1].push_back(e);
803 if (!incident_edges.count(v2))
804 vs.push_back(v2);
805 incident_edges[v2].push_back(e);
806 }
807 }
808
809 //for face query: traverse the face to get all edges
810 //and sort vertices in counter-clockwise order.
input_face(Face_const_handle fh)811 void input_face (Face_const_handle fh) const
812 {
813 Ccb_halfedge_const_circulator curr = fh->outer_ccb();
814 Ccb_halfedge_const_circulator circ = curr;
815 do {
816 CGAL_assertion(curr->face() == fh);
817 input_neighbor_f(curr);
818 } while (++curr != circ);
819
820 typename Arrangement_2::Hole_const_iterator hi;
821 for (hi = fh->holes_begin(); hi != fh->holes_end(); ++hi) {
822 Ccb_halfedge_const_circulator curr = *hi, circ = *hi;
823 do {
824 CGAL_assertion(curr->face() == fh);
825 input_neighbor_f(curr);
826 } while (++curr != circ);
827 }
828
829 std::sort(vs.begin(), vs.end(), Is_swept_earlier(q, geom_traits));
830
831 for (typename VHs::size_type i=0; i!=vs.size(); i++) {
832 typename VHs::size_type j = i+1;
833 while (j != vs.size()) {
834 if (!CGAL::collinear(q, vs[i]->point(), vs[j]->point()))
835 break;
836 j++;
837 }
838 if (j-i>1)
839 funnel(i, j);
840 i = j-1;
841 }
842 }
843 //for vertex or edge query: traverse the face to get all edges
844 //and sort vertices in counter-clockwise order.
input_face(Face_const_handle fh,EHs & good_edges,Arrangement_2 & bbox)845 void input_face (Face_const_handle fh,
846 EHs& good_edges,
847 Arrangement_2& bbox) const
848 {
849 Ccb_halfedge_const_circulator curr = fh->outer_ccb();
850 Ccb_halfedge_const_circulator circ = curr;
851 do {
852 CGAL_assertion(curr->face() == fh);
853 input_edge(curr, good_edges);
854 } while (++curr != circ);
855
856 typename Arrangement_2::Hole_const_iterator hi;
857 for (hi = fh->holes_begin(); hi != fh->holes_end(); ++hi) {
858 Ccb_halfedge_const_circulator curr = circ = *hi;
859 do {
860 CGAL_assertion(curr->face() == fh);
861 input_edge(curr, good_edges);
862 } while (++curr != circ);
863 }
864
865 //create a box that cover all vertices such that during the sweeping,
866 //the vision ray will always intersect at least an edge.
867 //this box doesn't intersect any relevant_edge.
868 Points points;
869 for (typename VHs::size_type i=0; i<vs.size(); i++) {
870 points.push_back(vs[i]->point());
871 }
872 points.push_back(q);
873 //first get the bounding box of all relevant points.
874 typename Geometry_traits_2::Iso_rectangle_2 bb =
875 bounding_box(points.begin(), points.end());
876
877 Number_type xmin, xmax, ymin, ymax;
878 typename Geometry_traits_2::Compute_x_2 compute_x =
879 geom_traits->compute_x_2_object();
880
881 typename Geometry_traits_2::Compute_y_2 compute_y =
882 geom_traits->compute_y_2_object();
883
884 //make the box a little bigger than bb so that it won't intersect any
885 //relevant_edge.
886 xmin = compute_x((bb.min)())-1;
887 ymin = compute_y((bb.min)())-1;
888 xmax = compute_x((bb.max)())+1;
889 ymax = compute_y((bb.max)())+1;
890 Point_2 box[4] = {Point_2(xmin, ymin), Point_2(xmax, ymin),
891 Point_2(xmax, ymax), Point_2(xmin, ymax)};
892
893 Halfedge_handle e1 = bbox.insert_in_face_interior(Segment_2(box[0], box[1]),
894 bbox.unbounded_face());
895
896 Halfedge_handle e2 = bbox.insert_from_left_vertex(Segment_2(box[1], box[2]),
897 e1->target());
898
899 Halfedge_handle e3 = bbox.insert_from_right_vertex(Segment_2(box[2],box[3]),
900 e2->target());
901
902 bbox.insert_at_vertices(Segment_2(box[0], box[3]),
903 e1->source(), e3->target());
904
905 circ = curr = e1->face()->outer_ccb();
906 do {
907 VH v = curr->target();
908 vs.push_back(v);
909 incident_edges[v].push_back(curr);
910 incident_edges[v].push_back(curr->next());
911 good_edges.push_back(curr);
912 } while(++curr != circ);
913
914 std::sort(vs.begin(), vs.end(), Is_swept_earlier(q, geom_traits));
915
916 for (typename VHs::size_type i=0; i!=vs.size(); i++) {
917 typename VHs::size_type j = i+1;
918 while (j != vs.size()) {
919 if (!CGAL::collinear(q, vs[i]->point(), vs[j]->point()))
920 break;
921 j++;
922 }
923 if (j-i>1)
924 funnel(i, j);
925 i = j-1;
926 }
927 }
928
929 };
930 } // end namespace CGAL
931
932
933
934 #endif
935
936
937