From b1d534560e8862902a170869ae1a781d30b79534 Mon Sep 17 00:00:00 2001 From: Ning Xu Date: Thu, 19 Jun 2014 12:59:49 -0400 Subject: [PATCH] change folders --- ...eprocessed_rotational_sweep_visibility_2.h | 106 +++ .../CGAL/Rotational_sweep_visibility_2.h | 866 ++++++++++++++++++ .../CGAL/Simple_polygon_visibility_2.h | 715 +++++++++++++++ .../CGAL/Triangular_expansion_visibility_2.h | 579 ++++++++++++ .../CGAL/Visibility_2/visibility_utils.h | 371 ++++++++ 5 files changed, 2637 insertions(+) create mode 100644 Visibility_2/include/CGAL/Preprocessed_rotational_sweep_visibility_2.h create mode 100644 Visibility_2/include/CGAL/Rotational_sweep_visibility_2.h create mode 100644 Visibility_2/include/CGAL/Simple_polygon_visibility_2.h create mode 100644 Visibility_2/include/CGAL/Triangular_expansion_visibility_2.h create mode 100644 Visibility_2/include/CGAL/Visibility_2/visibility_utils.h diff --git a/Visibility_2/include/CGAL/Preprocessed_rotational_sweep_visibility_2.h b/Visibility_2/include/CGAL/Preprocessed_rotational_sweep_visibility_2.h new file mode 100644 index 00000000000..730375cd1e6 --- /dev/null +++ b/Visibility_2/include/CGAL/Preprocessed_rotational_sweep_visibility_2.h @@ -0,0 +1,106 @@ +// Copyright (c) 2013 Technical University Braunschweig (Germany). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Author(s): Kan Huang +// + +#ifndef CGAL_PREPROCESSED_VISIBILITY_2_H +#define CGAL_PREPROCESSED_VISIBILITY_2_H + +#include +#include +#include +#include + +namespace CGAL { + +template +class Preprocessed_visibility_2 { + +public: + typedef typename Arrangement_2::Geometry_traits_2 Geometry_traits_2; + // Currently only consider with same type for both + typedef Arrangement_2 Input_Arrangement_2; + typedef Arrangement_2 Output_Arrangement_2; + + typedef typename Arrangement_2::Halfedge_const_handle Halfedge_const_handle; + typedef typename Arrangement_2::Ccb_halfedge_const_circulator Ccb_halfedge_const_circulator; + typedef typename Arrangement_2::Face_const_handle Face_const_handle; + typedef typename Arrangement_2::Kernel Kernel; + typedef typename CGAL::Arr_linear_traits_2 Linear_traits_2; + + typedef typename Geometry_traits_2::Point_2 Point_2; + typedef typename Geometry_traits_2::Ray_2 Ray_2; + typedef typename Geometry_traits_2::Segment_2 Segment_2; + typedef typename Geometry_traits_2::Line_2 Line_2; + typedef typename Geometry_traits_2::Vector_2 Vector_2; + typedef typename Geometry_traits_2::FT Number_type; + + typedef typename CGAL::Arrangement_2 Line_Arrangement_2; + + Preprocessed_visibility_2() : p_arr(NULL) {}; + + /*! Constructor given an arrangement and the Regularization tag. */ + Preprocessed_visibility_2(Input_Arrangement_2& arr/*, Regularization_tag r_t*/): p_arr(&arr) {}; + + bool is_attached() { + return (p_arr != NULL); + } + + void attach(Input_Arrangement_2& arr) { + p_arr = &arr; + } + + void detach() { + p_arr = NULL; + } + + Input_Arrangement_2 arr() { + return *p_arr; + } + + void compute_visibility(const Point_2& q, + const Face_const_handle face, + Output_Arrangement_2& out_arr + ) { + + } + + void compute_visibility(const Point_2& q, + const Halfedge_const_handle he, + Output_Arrangement_2& out_arr + ) { + +} + +private: + Input_Arrangement_2* arr; + Line_Arrangement_2 line_arr; + void preprocess() { + + } + + Line_2 dual_line(const Point_2& p) { + return Line_2(p.x(), -1, -p.y()); + } + +}; + +} // namespace CGAL + +#endif diff --git a/Visibility_2/include/CGAL/Rotational_sweep_visibility_2.h b/Visibility_2/include/CGAL/Rotational_sweep_visibility_2.h new file mode 100644 index 00000000000..fd378f8f69f --- /dev/null +++ b/Visibility_2/include/CGAL/Rotational_sweep_visibility_2.h @@ -0,0 +1,866 @@ +// Copyright (c) 2013 Technical University Braunschweig (Germany). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Author(s): Kan Huang +// + +#ifndef CGAL_ROTATIONAL_SWEEP_VISIBILITY_2_H +#define CGAL_ROTATIONAL_SWEEP_VISIBILITY_2_H + +#include +#include +#include +#include + + +namespace CGAL { + +template +class Rotational_sweep_visibility_2 { +public: + typedef Arrangement_2_ Arrangement_2; + typedef typename Arrangement_2::Traits_2 Traits_2; + typedef typename Arrangement_2::Geometry_traits_2 Geometry_traits_2; + typedef typename Arrangement_2::Vertex_const_handle Vertex_const_handle; + typedef typename Arrangement_2::Vertex_handle Vertex_handle; + typedef typename Arrangement_2::Halfedge_const_handle Halfedge_const_handle; + typedef typename Arrangement_2::Halfedge_handle Halfedge_handle; + typedef typename Arrangement_2::Ccb_halfedge_const_circulator + Ccb_halfedge_const_circulator; + typedef typename Arrangement_2::Face_const_handle Face_const_handle; + typedef typename Arrangement_2::Face_handle Face_handle; + + typedef typename Geometry_traits_2::Kernel K; + typedef typename Geometry_traits_2::Point_2 Point_2; + typedef typename Geometry_traits_2::Ray_2 Ray_2; + typedef typename Geometry_traits_2::Segment_2 Segment_2; + typedef typename Geometry_traits_2::Line_2 Line_2; + typedef typename Geometry_traits_2::Vector_2 Vector_2; + typedef typename Geometry_traits_2::Direction_2 Direction_2; + typedef typename Geometry_traits_2::FT Number_type; + typedef typename Geometry_traits_2::Object_2 Object_2; + + typedef RegularizationTag Regularization_tag; + typedef CGAL::Tag_true Supports_general_polygon_tag; + typedef CGAL::Tag_true Supports_simple_polygon_tag; + +private: + typedef std::vector Points; + typedef Vertex_const_handle VH; + typedef std::vector VHs; + typedef Halfedge_const_handle EH; + typedef std::vector EHs; + + class Less_edge: public std::binary_function { + const Geometry_traits_2* geom_traits; + public: + Less_edge() {} + Less_edge(const Geometry_traits_2* traits):geom_traits(traits) {} + bool operator() (const EH e1, const EH e2) const { + if (e1 == e2) + return false; + else { + return &(*e1)<&(*e2); +// if (e1->source() == e2->source()) +// return Visibility_2::compare_xy_2(geom_traits, e1->target()->point(), e2->target()->point()) == SMALLER; +// else +// return Visibility_2::compare_xy_2(geom_traits, e1->source()->point(), e2->source()->point()) == SMALLER; + } + } + }; + + class Less_vertex: public std::binary_function { + const Geometry_traits_2* geom_traits; + public: + Less_vertex() {} + Less_vertex(const Geometry_traits_2* traits):geom_traits(traits) {} + bool operator() (const VH v1, const VH v2) const { + if (v1 == v2) + return false; + else + // I know this is dirty but it speeds up by 25%. Michael + return &(*v1)<&(*v2); +// return Visibility_2::compare_xy_2(geom_traits, v1->point(), v2->point()) == SMALLER; + } + }; + + class Closer_edge: public std::binary_function { + const Geometry_traits_2* geom_traits; + Point_2 q; + public: + Closer_edge() {} + Closer_edge(const Geometry_traits_2* traits, const Point_2& q):geom_traits(traits), q(q) {} + + int vtype(const Point_2& c, const Point_2& p) const { + switch(Visibility_2::orientation_2(geom_traits, q, c, p)) { + case COLLINEAR: + if (Visibility_2::less_distance_to_point_2(geom_traits, q, c, p)) + return 0; + else + return 3; + case RIGHT_TURN: + return 1; + case LEFT_TURN: + return 2; + } + } + bool operator() (const EH& e1, const EH& e2) const { + if (e1 == e2) + return false; + const Point_2& s1=e1->source()->point(), + t1=e1->target()->point(), + s2=e2->source()->point(), + t2=e2->target()->point(); + if (e1->source() == e2->source()) { + + int vt1 = vtype(s1, t1), + vt2 = vtype(s1, t2); + if (vt1 != vt2) + return vt1 > vt2; + else + return (Visibility_2::orientation_2(geom_traits, s1, t2, t1)== + Visibility_2::orientation_2(geom_traits, s1, t2, q)); + } + + if (e1->target() == e2->source()) { +// const Point_2& p1 = s1, +// p2 = t2, +// c = s2; + int vt1 = vtype(t1, s1), + vt2 = vtype(t1, t2); + if (vt1 != vt2) + return vt1 > vt2; + else + return (Visibility_2::orientation_2(geom_traits, s2, t2, s1)== + Visibility_2::orientation_2(geom_traits, s2, t2, q)); + } + + if (e1->source() == e2->target()) { +// const Point_2& p1 = t1, +// p2 = s2, +// c = s1; + int vt1 = vtype(s1, t1), + vt2 = vtype(s1, s2); + if (vt1 != vt2) + return vt1 > vt2; + else return (Visibility_2::orientation_2(geom_traits, s1, s2, t1)== + Visibility_2::orientation_2(geom_traits, s1, s2, q)); + } + + if (e1->target() == e2->target()) { +// const Point_2& p1 = s1, +// p2 = s2, +// c = t1; + int vt1 = vtype(t1, s1), + vt2 = vtype(t1, s2); + if (vt1 != vt2) + return vt1 > vt2; + else return (Visibility_2::orientation_2(geom_traits, t1, s2, s1)== + Visibility_2::orientation_2(geom_traits, t1, s2, q)); + } + + Orientation e1q = Visibility_2::orientation_2(geom_traits, s1, t1, q); + switch (e1q) + { + case COLLINEAR: + if (Visibility_2::collinear(geom_traits, q, s2, t2)) { + //q is collinear with e1 and e2. + return (Visibility_2::less_distance_to_point_2(geom_traits, q, s1, s2) + || Visibility_2::less_distance_to_point_2(geom_traits, q, t1, t2)); + } + else { + //q is collinear with e1 not with e2. + if (Visibility_2::collinear(geom_traits, s2, t2, s1)) + return (Visibility_2::orientation_2(geom_traits, s2, t2, q) + == Visibility_2::orientation_2(geom_traits, s2, t2, t1)); + else + return (Visibility_2::orientation_2(geom_traits, s2, t2, q) + == Visibility_2::orientation_2(geom_traits, s2, t2, s1)); + } + case RIGHT_TURN: + switch (Visibility_2::orientation_2(geom_traits, s1, t1, s2)) { + case COLLINEAR: + return Visibility_2::orientation_2(geom_traits, s1, t1, t2)!=e1q; + case RIGHT_TURN: + if (Visibility_2::orientation_2(geom_traits, s1, t1, t2) == LEFT_TURN) + return Visibility_2::orientation_2(geom_traits, s2, t2, q) + == Visibility_2::orientation_2(geom_traits, s2, t2, s1); + else + return false; + case LEFT_TURN: + if (Visibility_2::orientation_2(geom_traits, s1, t1, t2) == RIGHT_TURN) + return Visibility_2::orientation_2(geom_traits, s2, t2, q) + == Visibility_2::orientation_2(geom_traits, s2, t2, s1); + else + return true; + } + case LEFT_TURN: + switch (Visibility_2::orientation_2(geom_traits, s1, t1, s2)) { + case COLLINEAR: + return Visibility_2::orientation_2(geom_traits, s1, t1, t2)!=e1q; + case LEFT_TURN: + if (Visibility_2::orientation_2(geom_traits, s1, t1, t2) == RIGHT_TURN) + return Visibility_2::orientation_2(geom_traits, s2, t2, q) + == Visibility_2::orientation_2(geom_traits, s2, t2, s1); + else + return false; + case RIGHT_TURN: + if (Visibility_2::orientation_2(geom_traits, s1, t1, t2) == LEFT_TURN) + return Visibility_2::orientation_2(geom_traits, s2, t2, q) + == Visibility_2::orientation_2(geom_traits, s2, t2, s1); + else + return true; + } + } + } + + }; + +// Using hash_map or edx causes a seg fault, did not have the time to see why. Michael +// class Hash_edge: public std::unary_function::result_type> { +// public: +// typename boost::hash::result_type +// operator() (const EH e1) const { +// return boost::hash()(&(e1->curve())); +// } +// }; + + + const Geometry_traits_2 *geom_traits; + const Arrangement_2 *p_arr; + Point_2 q; //query point + Points polygon; //visibility polygon + std::map incident_edges; //the edges that are + std::map edx; //index of active edges in the heap + // boost::unordered_map edx; //index of active edges in the heap + std::set active_edges; //a set of edges that intersect the current vision ray. + VHs vs; //angular sorted vertices + EHs bad_edges; //edges that pass the query point + VH cone_end1; //an end of visibility cone + VH cone_end2; //another end of visibility cone + int cone_end1_idx; //index of cone_end1->point() in visibility polygon + int cone_end2_idx; //index of cone_end2->point() in visibility polygon + + bool is_vertex_query; + bool is_edge_query; + bool is_face_query; + bool is_big_cone; //whether the angle of visibility_cone is greater than pi. + +public: + Rotational_sweep_visibility_2(): p_arr(NULL), geom_traits(NULL) {} + Rotational_sweep_visibility_2(const Arrangement_2& arr): p_arr(&arr) { + geom_traits = p_arr->geometry_traits(); + } + + const std::string name(){return std::string("R_visibility_2");} + + template + typename VARR::Face_handle + compute_visibility(const Point_2& q, const Halfedge_const_handle e, VARR& arr_out) { + arr_out.clear(); + bad_edges.clear(); + this->q = q; + + if (Visibility_2::compare_xy_2(geom_traits, q, e->target()->point())==EQUAL) { + is_vertex_query = true; + is_edge_query = false; + is_face_query = false; + cone_end1 = e->source(); + cone_end2 = e->next()->target(); + is_big_cone = CGAL::right_turn(cone_end1->point(), q, cone_end2->point()); + + typename Arrangement_2::Halfedge_around_vertex_const_circulator first, curr; + first = curr = e->target()->incident_halfedges(); + do { + if (curr->face() == e->face()) + bad_edges.push_back(curr); + else if (curr->twin()->face() == e->face()) + bad_edges.push_back(curr->twin()); + } while (++curr != first); + } + else { + is_vertex_query = false; + is_edge_query = true; + is_face_query = false; + cone_end1 = e->source(); + cone_end2 = e->target(); + bad_edges.push_back(e); + is_big_cone = false; + } + visibility_region_impl(e->face(), q); + + //decide which inside of the visibility butterfly is needed. + int small_idx, big_idx; + if ( cone_end1_idx < cone_end2_idx ) { + small_idx = cone_end1_idx; + big_idx = cone_end2_idx; + } + else { + small_idx = cone_end2_idx; + big_idx = cone_end1_idx; + } + int next_idx = small_idx + 1; + bool is_between; + //indicate whether the shape between small_idx and big_idx is the visibility region required. + if (CGAL::right_turn(cone_end1->point(), q, cone_end2->point())) { + is_between = false; + while (next_idx != big_idx) { + if (CGAL::left_turn(cone_end1->point(), q, polygon[next_idx]) || CGAL::left_turn(q, cone_end2->point(), polygon[next_idx])) { + is_between = true; + break; + } + next_idx++; + } + } + else { + is_between = true; + while (next_idx != big_idx) { + if (CGAL::right_turn(cone_end1->point(), q, polygon[next_idx]) || CGAL::right_turn(q, cone_end2->point(), polygon[next_idx])) { + is_between = false; + break; + } + next_idx++; + } + } + + typename Points::iterator first = polygon.begin() + small_idx; + typename Points::iterator last = polygon.begin() + big_idx; + if (is_between) { + Points polygon_out(first, last+1); + if (is_vertex_query) + polygon_out.push_back(q); + Visibility_2::report_while_handling_needles + (geom_traits, q, polygon_out, arr_out); + } + else { + Points polygon_out(polygon.begin(), first+1); + if (is_vertex_query) polygon_out.push_back(q); + for (int i = big_idx; i != polygon.size(); i++) { + polygon_out.push_back(polygon[i]); + } + Visibility_2::report_while_handling_needles + (geom_traits, q, polygon_out, arr_out); + } + + conditional_regularize(arr_out, Regularization_tag()); + + if (arr_out.faces_begin()->is_unbounded()) + return ++arr_out.faces_begin(); + else + return arr_out.faces_begin(); + } + + template + typename VARR::Face_handle + compute_visibility(const Point_2& q, const Face_const_handle f, VARR& arr_out) { + arr_out.clear(); + this->q = q; + is_vertex_query = false; + is_edge_query = false; + is_face_query = true; + + visibility_region_impl(f, q); + Visibility_2::report_while_handling_needles(geom_traits, q, polygon, arr_out); + conditional_regularize(arr_out, Regularization_tag()); + if (arr_out.faces_begin()->is_unbounded()) + return ++arr_out.faces_begin(); + else + return arr_out.faces_begin(); + } + +bool is_attached() { + return (p_arr != NULL); +} + +void attach(const Arrangement_2& arr) { + p_arr = &arr; + geom_traits = p_arr->geometry_traits(); +} + +void detach() { + p_arr = NULL; + geom_traits = NULL; +} + +const Arrangement_2& arr() { + return *p_arr; +} + +private: + //get the neighbor of v along edge e + VH get_neighbor(const EH e, const VH v) { + if (e->source() == v) + return e->target(); + else + return e->source(); + } + //check whether ray(q->dp) intersects segment(p1, p2) + bool do_intersect_ray(const Point_2& q, + const Point_2& dp, + const Point_2& p1, + const Point_2& p2) { + return (CGAL::orientation(q, dp, p1) != CGAL::orientation(q, dp, p2) && CGAL::orientation(q, p1, dp) == CGAL::orientation(q, p1, p2)); + } + + //arrange vertices that on a same vision ray in a 'funnel' order + void funnel(int i, int j) { + VHs right, left; + //whether the edges incident to a vertex block the left side and right side of current vision ray. + bool block_left(false), block_right(false); + VH former = vs[i], nb; + for (int l=i; lpoint(), nb->point())) + left_v = true; + else + right_v = CGAL::right_turn(q, vs[l]->point(), nb->point()); + } + if (has_predecessor) { + //if the current vertex connects to the vertex before by an edge, + //the vertex before can help it to block. + block_left = block_left || left_v; + block_right = block_right || right_v; + } + else { + block_left = left_v; + block_right = right_v; + } + if (block_left && block_right) { + //when both sides are blocked, there is no need to change the vertex after. + right.push_back(vs[l]); + break; + } + else { + if (block_left) + left.push_back(vs[l]); + else + right.push_back(vs[l]); + } + former = vs[l]; + } + for (int l=0; l!=right.size(); l++) + vs[i+l] = right[l]; + for (int l=0; l!=left.size(); l++) + vs[i+l+right.size()] = left[left.size()-1-l]; + } + + + + void visibility_region_impl(const Face_const_handle f, const Point_2& q) { + vs.clear(); + polygon.clear(); + active_edges = std::set(Closer_edge(geom_traits, q)); + incident_edges = std::map(Less_vertex(geom_traits)); + edx = std::map(Less_edge(geom_traits)); + + EHs relevant_edges; //all edges that can affect the visibility of query point. + Arrangement_2 bbox; + if (is_face_query) + input_face(f); + else + input_face(f, relevant_edges, bbox); + //the following code is the initiation of vision ray. the direction of the initial ray is between the direction + //from q to last vertex in vs and positive x-axis. By choosing this direction, we make + //sure that all plane is swept and there is not needle at the beginning of sweeping. + Vector_2 dir; + if (Direction_2(-1, 0) < Direction_2(Vector_2(q, vs.back()->point()))) + dir = Vector_2(1, 0) + Vector_2(q, vs.back()->point()); + else + dir = Vector_2(0, -1); + Point_2 dp = q + dir; + + //initiation of active_edges. for face queries, all edges on the boundary can affect visibility. + //for non-face queries, only relevant_edges has to be considered. + if (is_face_query) { + Ccb_halfedge_const_circulator curr = f->outer_ccb(); + Ccb_halfedge_const_circulator circ = curr; + do { + if (do_intersect_ray(q, dp, curr->target()->point(), curr->source()->point())) { + active_edges.insert(curr); + } + } while (++curr != circ); + + typename Arrangement_2::Hole_const_iterator hi; + for (hi = f->holes_begin(); hi != f->holes_end(); ++hi) { + Ccb_halfedge_const_circulator curr = circ = *hi; + do { + if (do_intersect_ray(q, dp, curr->target()->point(), curr->source()->point())) + active_edges.insert(curr); + } while (++curr != circ); + } + } + else { + for (int i=0; i!=relevant_edges.size(); i++) + if (do_intersect_ray(q, dp, relevant_edges[i]->source()->point(), relevant_edges[i]->target()->point())) + active_edges.insert(relevant_edges[i]); + } + + //angular sweep begins +// std::cout<(ctemp_eh); + temp_eh = insert_ehs.front(); + } + else { + for (int j=0; j!=remove_cnt; j++) + active_edges.erase(remove_ehs[j]); + for (int j=0; j!=insert_cnt; j++) + active_edges.insert(insert_ehs[j]); + } + + if (closest_e != *active_edges.begin()) { + //when the closest edge changed + if (is_face_query) { + if (remove_cnt > 0 && insert_cnt > 0) { + //some edges are added and some are deleted, which means the vertex swept is part of visibility polygon. + update_visibility(vh->point()); + } + if (remove_cnt == 0 && insert_cnt > 0) { + //only add some edges, means the view ray is blocked by new edges. + //therefore first add the intersection of view ray and former closet edge, then add the vertice swept. + update_visibility(ray_seg_intersection(q, + vh->point(), + closest_e->target()->point(), + closest_e->source()->point())); + update_visibility(vh->point()); + } + if (remove_cnt > 0 && insert_cnt == 0) { + //only delete some edges, means some block is moved and the view ray can reach the segments after the block. + update_visibility(vh->point()); + update_visibility(ray_seg_intersection(q, + vh->point(), + (*active_edges.begin())->target()->point(), + (*active_edges.begin())->source()->point())); + } + } + else { + //extra work here for edge/vertex query is the index of cone_end1 and cone_end2 will be recorded. + if (remove_cnt > 0 && insert_cnt > 0) { + //some edges are added and some are deleted, which means the vertice swept is part of visibility polygon. + if (update_visibility(vh->point())) { + if (vh == cone_end1) + cone_end1_idx = polygon.size()-1; + else if (vh == cone_end2) + cone_end2_idx = polygon.size()-1; + } + } + if (remove_cnt == 0 && insert_cnt > 0) { + //only add some edges, means the view ray is blocked by new edges. + //therefore first add the intersection of view ray and former closet edge, then add the vertice swept. + update_visibility(ray_seg_intersection(q, + vh->point(), + closest_e->target()->point(), + closest_e->source()->point())); + if (update_visibility(vh->point())) { + if (vh == cone_end1) + cone_end1_idx = polygon.size()-1; + else if (vh == cone_end2) + cone_end2_idx = polygon.size()-1; + } + } + if (remove_cnt > 0 && insert_cnt == 0) { + //only delete some edges, means some block is removed and the vision ray can reach the segments after the block. + if (update_visibility(vh->point())) { + if (vh == cone_end1) + cone_end1_idx = polygon.size()-1; + else if (vh == cone_end2) + cone_end2_idx = polygon.size()-1; + } + update_visibility(ray_seg_intersection(q, + vh->point(), + (*active_edges.begin())->target()->point(), + (*active_edges.begin())->source()->point())); + } + } + } + } + } + + void print_edge(const EH e) { + std::cout<source()->point()<<"->"<target()->point()<dp) and segment(s, t) + //if they are collinear then return the endpoint which is closer to q. + + Point_2 ray_seg_intersection( + const Point_2& q, const Point_2& dp, // the ray + const Point_2& s, const Point_2& t) // the segment + { + if (CGAL::collinear(q, dp, s)) { + if (CGAL::collinear(q, dp, t)) { + if (CGAL::compare_distance_to_point(q, s, t)==CGAL::SMALLER) + return s; + else + return t; + } + else + return s; + } + Ray_2 ray(q,dp); + Segment_2 seg(s,t); + CGAL::Object result = CGAL::intersection(ray, seg); + return *(CGAL::object_cast(&result)); + } + + //check if p has been discovered before, if not update the visibility polygon + bool update_visibility(const Point_2& p){ + if (polygon.empty()) { + polygon.push_back(p); + return true; + } + else if (Visibility_2::compare_xy_2(geom_traits, polygon.back(), p) != EQUAL) { + polygon.push_back(p); + return true; + } + return false; + } + + //functor to decide which vertex is swept earlier by the rotational sweeping ray + class Is_swept_earlier:public std::binary_function { + const Point_2& q; + const Geometry_traits_2* geom_traits; + public: + Is_swept_earlier(const Point_2& q, const Geometry_traits_2* traits):q(q), geom_traits(traits) {} + bool operator() (const VH v1, const VH v2) const { + const Point_2& p1 = v1->point(); + const Point_2& p2 = v2->point(); + int qua1 = quadrant(q, p1); + int qua2 = quadrant(q, p2); + if (qua1 < qua2) + return true; + if (qua1 > qua2) + return false; + if (collinear(q, p1, p2)) + return (CGAL::compare_distance_to_point(q, p1, p2) == CGAL::SMALLER); + else + return CGAL::right_turn(p1, q, p2); + } + + //return the quadrant of p with respect to o. + int quadrant(const Point_2& o, const Point_2& p) const { + typename Geometry_traits_2::Compare_x_2 compare_x = geom_traits->compare_x_2_object(); + typename Geometry_traits_2::Compare_y_2 compare_y = geom_traits->compare_y_2_object(); + + Comparison_result dx = compare_x(p, o); + Comparison_result dy = compare_y(p, o); + if (dx==LARGER && dy!=SMALLER) + return 1; + if (dx!=LARGER && dy==LARGER) + return 2; + if (dx==SMALLER && dy!=LARGER) + return 3; + if (dx!=SMALLER && dy==SMALLER) + return 4; + return 0; + } + }; + + //when the query point is in face, every edge is good. + void input_neighbor_f( const Halfedge_const_handle e) { + VH v = e->target(); + if (!incident_edges.count(v)) + vs.push_back(v); + incident_edges[v].push_back(e); + incident_edges[v].push_back(e->next()); + } + + //check if p is in the visibility cone + bool is_in_cone(const Point_2& p) const{ + if (is_big_cone) + return (!CGAL::right_turn(cone_end1->point(), q, p)) || (!CGAL::left_turn(cone_end2->point(), q, p)); + else + return (!CGAL::right_turn(cone_end1->point(), q, p)) && (!CGAL::left_turn(cone_end2->point(), q, p)); + } + + //for vertex and edge query: the visibility is limited in a cone. + void input_edge(const Halfedge_const_handle e, + EHs& good_edges) { + for (int i=0; itarget(); + VH v2 = e->source(); + //an edge will affect visibility only if it has an endpoint in the visibility cone or it crosses the boundary of the cone. + if (is_in_cone(v1->point()) || is_in_cone(v2->point()) || do_intersect_ray(q, cone_end1->point(), v1->point(), v2->point())) { + good_edges.push_back(e); + if (!incident_edges.count(v1)) + vs.push_back(v1); + incident_edges[v1].push_back(e); + if (!incident_edges.count(v2)) + vs.push_back(v2); + incident_edges[v2].push_back(e); + } + } + + //for face query: traverse the face to get all edges + //and sort vertices in counter-clockwise order. + void input_face (Face_const_handle fh) + { + Ccb_halfedge_const_circulator curr = fh->outer_ccb(); + Ccb_halfedge_const_circulator circ = curr; + do { + assert(curr->face() == fh); + input_neighbor_f(curr); + } while (++curr != circ); + + typename Arrangement_2::Hole_const_iterator hi; + for (hi = fh->holes_begin(); hi != fh->holes_end(); ++hi) { + Ccb_halfedge_const_circulator curr = *hi, circ = *hi; + do { + assert(curr->face() == fh); + input_neighbor_f(curr); + } while (++curr != circ); + } + + std::sort(vs.begin(), vs.end(), Is_swept_earlier(q, geom_traits)); + + for (int i=0; i!=vs.size(); i++) { + int j = i+1; + while (j != vs.size()) { + if (!CGAL::collinear(q, vs[i]->point(), vs[j]->point())) + break; + j++; + } + if (j-i>1) + funnel(i, j); + i = j-1; + } + } + //for vertex or edge query: traverse the face to get all edges + //and sort vertices in counter-clockwise order. + void input_face (Face_const_handle fh, + EHs& good_edges, + Arrangement_2& bbox) + { + Ccb_halfedge_const_circulator curr = fh->outer_ccb(); + Ccb_halfedge_const_circulator circ = curr; + do { + assert(curr->face() == fh); + input_edge(curr, good_edges); + } while (++curr != circ); + + typename Arrangement_2::Hole_const_iterator hi; + for (hi = fh->holes_begin(); hi != fh->holes_end(); ++hi) { + Ccb_halfedge_const_circulator curr = circ = *hi; + do { + assert(curr->face() == fh); + input_edge(curr, good_edges); + } while (++curr != circ); + } + + //create a box that cover all vertices such that during the sweeping, the vision ray will always intersect at least an edge. + //this box doesn't intersect any relevant_edge. + Points points; + for (int i=0; ipoint()); + } + points.push_back(q); + //first get the bounding box of all relevant points. + typename Geometry_traits_2::Iso_rectangle_2 bb = bounding_box(points.begin(), points.end()); + + Number_type xmin, xmax, ymin, ymax; + typename Geometry_traits_2::Compute_x_2 compute_x = geom_traits->compute_x_2_object(); + typename Geometry_traits_2::Compute_y_2 compute_y = geom_traits->compute_y_2_object(); + + //make the box a little bigger than bb so that it won't intersect any relevant_edge. + xmin = compute_x(bb.min())-1; + ymin = compute_y(bb.min())-1; + xmax = compute_x(bb.max())+1; + ymax = compute_y(bb.max())+1; + Point_2 box[4] = {Point_2(xmin, ymin), Point_2(xmax, ymin), + Point_2(xmax, ymax), Point_2(xmin, ymax)}; + Halfedge_handle e1 = bbox.insert_in_face_interior(Segment_2(box[0], box[1]), bbox.unbounded_face()); + Halfedge_handle e2 = bbox.insert_from_left_vertex(Segment_2(box[1], box[2]), e1->target()); + Halfedge_handle e3 = bbox.insert_from_right_vertex(Segment_2(box[2], box[3]), e2->target()); + bbox.insert_at_vertices(Segment_2(box[0], box[3]), e1->source(), e3->target()); + + circ = curr = e1->face()->outer_ccb(); + do { + VH v = curr->target(); + vs.push_back(v); + incident_edges[v].push_back(curr); + incident_edges[v].push_back(curr->next()); + good_edges.push_back(curr); + } while(++curr != circ); + + std::sort(vs.begin(), vs.end(), Is_swept_earlier(q, geom_traits)); + + for (int i=0; i!=vs.size(); i++) { + int j = i+1; + while (j != vs.size()) { + if (!CGAL::collinear(q, vs[i]->point(), vs[j]->point())) + break; + j++; + } + if (j-i>1) + funnel(i, j); + i = j-1; + } + } + + template + void conditional_regularize(VARR& arr_out, CGAL::Tag_true) { + regularize_output(arr_out); + } + + template + void conditional_regularize(VARR& arr_out, CGAL::Tag_false) { + //do nothing + } + + template + void regularize_output(VARR& arr_out) { + typename VARR::Edge_iterator e_itr; + for (e_itr = arr_out.edges_begin(); + e_itr != arr_out.edges_end(); + e_itr++) { + if (e_itr->face() == e_itr->twin()->face()) + arr_out.remove_edge(e_itr); + } + } +}; +} // end namespace CGAL + + + +#endif + + diff --git a/Visibility_2/include/CGAL/Simple_polygon_visibility_2.h b/Visibility_2/include/CGAL/Simple_polygon_visibility_2.h new file mode 100644 index 00000000000..0934d03ba1c --- /dev/null +++ b/Visibility_2/include/CGAL/Simple_polygon_visibility_2.h @@ -0,0 +1,715 @@ +// Copyright (c) 2013 Technical University Braunschweig (Germany). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$` +// $Id$ +// +// +// Author(s): Francisc Bungiu +// Michael Hemmer +// Ning Xu + +#ifndef CGAL_SIMPLE_POLYGON_VISIBILITY_2_H +#define CGAL_SIMPLE_POLYGON_VISIBILITY_2_H + +#include +#include +#include +#include +#include +#include +#include + +namespace CGAL { + +template +class Simple_polygon_visibility_2 { + +public: + // Currently only consider with same type for both + typedef Arrangement_2_ Arrangement_2; + typedef typename Arrangement_2::Traits_2 Traits_2; + typedef typename Arrangement_2::Geometry_traits_2 Geometry_traits_2; + typedef typename Geometry_traits_2::Kernel K; + + typedef typename Arrangement_2::Halfedge_const_handle + Halfedge_const_handle; + typedef typename Arrangement_2::Halfedge_handle Halfedge_handle; + typedef typename Arrangement_2::Ccb_halfedge_const_circulator + Ccb_halfedge_const_circulator; + typedef typename Arrangement_2::Face_const_handle Face_const_handle; + typedef typename Arrangement_2::Face_handle Face_handle; + typedef typename Arrangement_2::Halfedge_around_vertex_const_circulator + Halfedge_around_vertex_const_circulator; + + typedef typename Geometry_traits_2::Point_2 Point_2; + typedef typename Geometry_traits_2::Ray_2 Ray_2; + typedef typename Geometry_traits_2::Segment_2 Segment_2; + typedef typename Geometry_traits_2::Line_2 Line_2; + typedef typename Geometry_traits_2::Vector_2 Vector_2; + typedef typename Geometry_traits_2::Direction_2 Direction_2; + typedef typename Geometry_traits_2::FT Number_type; + typedef typename Geometry_traits_2::Object_2 Object_2; + + typedef RegularizationTag Regularization_tag; + typedef CGAL::Tag_false Supports_general_polygon_tag; + typedef CGAL::Tag_true Supports_simple_polygon_tag; + + Simple_polygon_visibility_2() : p_arr(NULL), geom_traits(NULL) {}; + + /*! Constructor given an arrangement and the Regularization tag. */ + Simple_polygon_visibility_2(const Arrangement_2& arr): + p_arr(&arr) { + geom_traits = p_arr->geometry_traits(); + query_pt_is_vertex = false; + query_pt_is_on_halfedge = false; + } + + + const std::string name(){return std::string("S_visibility_2");} + + /*! Method to check if the visibility object is attached or not to + an arrangement*/ + bool is_attached() { + return (p_arr != NULL); + } + + /*! Attaches the visibility object to the 'arr' arrangement */ + void attach(const Arrangement_2& arr) { + p_arr = &arr; + geom_traits = p_arr->geometry_traits(); + query_pt_is_vertex = false; + query_pt_is_on_halfedge = false; + } + + /*! Detaches the visibility object from the arrangement it is + attached to*/ + void detach() { + p_arr = NULL; + geom_traits = NULL; + vertices.clear(); + query_pt_is_vertex = false; + query_pt_is_on_halfedge = false; + p_cdt = boost::shared_ptr(); + } + + /*! Getter method for the input arrangement*/ + const Arrangement_2& arr() { + return *p_arr; + } + + /*! Computes the visibility object from the query point 'q' in the face + 'face' and constructs the output in 'out_arr'*/ + template + typename VARR::Face_handle + compute_visibility(const Point_2& q, + const Face_const_handle face, + VARR& out_arr) { + + assert(query_pt_is_vertex == false); + assert(query_pt_is_on_halfedge == false); + + // Now retrieve the circulator to first visible vertex from triangulation + Ccb_halfedge_const_circulator circ = find_visible_start(face, q); + Ccb_halfedge_const_circulator curr = circ; + Halfedge_const_handle he; + + do { + he = curr; + vertices.push_back(he->source()->point()); + } while(++curr != circ); + + vertices.push_back(vertices[0]); + + visibility_region_impl(q); + + typename std::vector points; + while (!s.empty()) { + Point_2 curr_point = s.top(); + points.push_back(curr_point); + s.pop(); + } + + std::reverse(points.begin(), points.end()); + + CGAL::Visibility_2::report_while_handling_needles + (geom_traits, + q, + points, + out_arr); + + CGAL_precondition(out_arr.number_of_isolated_vertices() == 0); + CGAL_precondition(s.size() == 0); + conditional_regularize(out_arr, Regularization_tag()); + vertices.clear(); + if (out_arr.faces_begin()->is_unbounded()) { + return ++out_arr.faces_begin(); + } + else { + return out_arr.faces_begin(); + } + } + + /*! Computes the visibility region of the query point 'q' located on the + halfedge 'he' and constructs the output in 'out_arr'*/ + template + typename VARR::Face_handle + compute_visibility( + const Point_2& q, + const Halfedge_const_handle he, + VARR& out_arr ) + { + query_pt_is_vertex = false; + query_pt_is_on_halfedge = false; + bool query_on_target = false; + + if (q != he->source()->point()) { + if (q != he->target()->point()) { + vertices.push_back(he->target()->point()); + query_pt_is_on_halfedge = true; + } + else { + query_pt_is_vertex = true; + query_on_target = true; + } + } else { + vertices.push_back( he->target()->point() ); + query_pt_is_vertex = true; + } + + Ccb_halfedge_const_circulator circ = he; + circ++; + Ccb_halfedge_const_circulator curr = circ; + + do { + Halfedge_const_handle he_handle = curr; + Point_2 curr_vertex = he_handle->target()->point(); + vertices.push_back(curr_vertex); + } while (++curr != circ); + + if ( query_on_target ) { + vertices.push_back( vertices[0] ); + } + + visibility_region_impl(q); + + typename std::vector points; + if (!s.empty()) { + Point_2 prev_pt = s.top(); + if (prev_pt != q) { + points.push_back(prev_pt); + } + else if (query_pt_is_vertex) { + points.push_back(prev_pt); + } + if (!s.empty()) { + s.pop(); + } + while(!s.empty()) { + Point_2 curr_pt = s.top(); + if (curr_pt != q) { + points.push_back(curr_pt); + } + else if (query_pt_is_vertex) { + points.push_back(curr_pt); + } + s.pop(); + } + } + + std::reverse(points.begin(), points.end()); + + CGAL::Visibility_2::report_while_handling_needles + (geom_traits, + q, + points, + out_arr); + CGAL_precondition(out_arr.number_of_isolated_vertices() == 0); + CGAL_precondition(s.size() == 0); + conditional_regularize(out_arr, Regularization_tag()); + vertices.clear(); + if (out_arr.faces_begin()->is_unbounded()) { + return ++out_arr.faces_begin(); + } + else { + return out_arr.faces_begin(); + } + } + +private: + typedef CGAL::Triangulation_vertex_base_2 Vb; + typedef CGAL::Constrained_triangulation_face_base_2 Fb; + typedef CGAL::Triangulation_data_structure_2 TDS; + typedef CGAL::No_intersection_tag Itag; + typedef CGAL::Constrained_triangulation_2 CDT; + +private: + const Arrangement_2 *p_arr; + /*! Boost pointer to the constrained Delaunay triangulation object*/ + boost::shared_ptr p_cdt; + /*! Mapping of the vertices of the input to the corresponding circulator + needed for finding the first visible vertex in case of face queries*/ + std::map + point_itr_map; + const Geometry_traits_2 *geom_traits; + /*! Stack of visibile points; manipulated when going through the sequence + of input vertices; contains the vertices of the visibility region after + the run of the algorithm*/ + std::stack s; + /*! Sequence of input vertices*/ + std::vector vertices; + /*! State of visibility region algorithm*/ + enum {LEFT, RIGHT, SCANA, SCANB, SCANC, SCAND, FINISH} upcase; + bool query_pt_is_vertex; + bool query_pt_is_on_halfedge; + + /*! Regularize output if flag is set to true*/ + template + void conditional_regularize(VARR& out_arr, CGAL::Tag_true) { + regularize_output(out_arr); + } + /*! No need to regularize output if flag is set to false*/ + template + void conditional_regularize(VARR& out_arr, CGAL::Tag_false) { + //do nothing + } + + /*! Regularizes the output - removes edges that have the same face on both + sides */ + template + void regularize_output(VARR& out_arr) { + typename VARR::Edge_iterator e_itr; + for (e_itr = out_arr.edges_begin() ; + e_itr != out_arr.edges_end() ; e_itr++) { + + typename VARR::Halfedge_handle he = e_itr; + typename VARR::Halfedge_handle he_twin = he->twin(); + if (he->face() == he_twin->face()) { + out_arr.remove_edge(he); + } + } + } + + /*! Initialized the constrained Delaunay triangulation using the edges of + the outer boundary of 'face' */ + void init_cdt(const Face_const_handle &face) { + + std::vector > constraints; + typename Arrangement_2::Ccb_halfedge_const_circulator circ = + face->outer_ccb(); + typename Arrangement_2::Ccb_halfedge_const_circulator curr = circ; + typename Arrangement_2::Halfedge_const_handle he; + + do { + he = curr; + Point_2 source = he->source()->point(); + Point_2 target = he->target()->point(); + point_itr_map.insert(std::make_pair(source, curr)); + constraints.push_back(std::make_pair(source,target)); + } while(++curr != circ); + + p_cdt = boost::shared_ptr(new CDT(constraints.begin(),constraints.end())); + } + + + /*! Finds a visible vertex from the query point 'q' in 'face' + to start the algorithm from*/ + Ccb_halfedge_const_circulator find_visible_start(Face_const_handle face, const Point_2 &q) { + init_cdt(face); + typename CDT::Face_handle fh = p_cdt->locate(q); + Point_2 start_point = fh->vertex(0)->point(); + + // Now retrieve the circulator to first visible vertex from triangulation + Ccb_halfedge_const_circulator circ = point_itr_map[start_point]; + Halfedge_const_handle he_curr = circ; + + Halfedge_around_vertex_const_circulator incident_circ = he_curr->source()->incident_halfedges(); + Halfedge_around_vertex_const_circulator incident_curr = incident_circ; + + do { + Ccb_halfedge_const_circulator curr_inc = incident_curr; + Halfedge_const_handle he_curr_inc = curr_inc; + + if (he_curr_inc->face() == face) { + Ccb_halfedge_const_circulator incident_next = incident_curr; + incident_next++; + Halfedge_const_handle he_next_inc = incident_next; + + if (CGAL::Visibility_2::orientation_2(geom_traits, + he_curr_inc->source()->point(), + he_curr_inc->target()->point(), + q) == CGAL::LEFT_TURN + || CGAL::Visibility_2::orientation_2(geom_traits, + he_next_inc->source()->point(), + he_next_inc->target()->point(), + q) == CGAL::LEFT_TURN) { + Ccb_halfedge_const_circulator result_circ = incident_next; + Halfedge_const_handle he_print = result_circ; + return result_circ; + } + } + } while (++incident_curr != incident_circ); + } + + /*! Main method of the algorithm - initializes the stack and variables + and calles the corresponding methods acc. to the algorithm's state; + 'q' - query point; + 'i' - current vertex' index + 'w' - endpoint of ray shot from query point */ + void visibility_region_impl(const Point_2& q) { + int i = 0; + Point_2 w; + CGAL::Orientation orient = CGAL::Visibility_2::orientation_2(geom_traits, + q, + vertices[0], + vertices[1]); + if ( orient != CGAL::RIGHT_TURN ) { + upcase = LEFT; + i = 1; + w = vertices[1]; + s.push(vertices[0]); + s.push(vertices[1]); + } + else { + upcase = SCANA; + i = 1; + w = vertices[1]; + s.push(vertices[0]); + } + + Ray_2 ray_origin( q, vertices[0] ); + do { + switch(upcase) { + case LEFT: + left(i, w, q); + break; + case RIGHT: + right(i, w, q); + break; + case SCANA: + scana(i, w, q); + break; + case SCANB: + scanb(i, w, q); + break; + case SCANC: + scanc(i, w, q); + break; + case SCAND: + scand(i, w, q); + break; + } + if ( upcase == LEFT ) { + Point_2 s_t = s.top(); + s.pop(); + if ( ( CGAL::Visibility_2::orientation_2 + ( geom_traits, q, vertices[0], s.top() ) == CGAL::RIGHT_TURN ) && + ( CGAL::Visibility_2::orientation_2 + ( geom_traits, q, vertices[0],s_t ) == CGAL::LEFT_TURN ) ) { + Segment_2 seg( s.top(), s_t ); + if ( CGAL::Visibility_2::do_intersect_2 + + ( geom_traits, seg, ray_origin ) ) { + Object_2 result = CGAL::Visibility_2::intersect_2 + + ( geom_traits, seg, ray_origin ); + const Point_2 * ipoint = CGAL::object_cast(&result); + assert( ipoint != NULL ); + s_t = *ipoint; + upcase = SCANB; + } + } + s.push( s_t ); + } + } while(upcase != FINISH); + } + + /*! Method that handles the left turns in the vertex algorithm */ + void left(int& i, Point_2& w, const Point_2& query_pt) { + if (i >= vertices.size() - 1) { + upcase = FINISH; + } + else { + Point_2 s_t = s.top(); + s.pop(); + Point_2 s_t_prev = s.top(); + s.push( s_t ); + CGAL::Orientation orient1 = CGAL::Visibility_2::orientation_2 + + ( geom_traits, + query_pt, + vertices[i], + vertices[i+1] ); + if ( orient1 != CGAL::RIGHT_TURN ) { + // Case L2 + upcase = LEFT; + s.push( vertices[i+1] ); + w = vertices[i+1]; + i++; + } else { + CGAL::Orientation orient2 = CGAL::Visibility_2::orientation_2 + + ( geom_traits, + s_t_prev, + vertices[i], + vertices[i+1] ); + if ( orient2 == CGAL::RIGHT_TURN ) { + // Case L3 + upcase = SCANA; + w = vertices[i+1]; + i++; + } else { + // Case L4 + upcase = RIGHT; + w = vertices[i]; + i++; + } + } + } + } + + /*! Scans the stack such that all vertices that were pushed before to the + stack and are now not visible anymore. */ + void right(int& i, Point_2& w, const Point_2& query_pt) { + Point_2 s_j; + Point_2 s_j_prev; + Point_2 u; + int mode = 0; + CGAL::Orientation orient1, orient2; + + s_j_prev = s.top(); + orient2 = CGAL::Visibility_2::orientation_2 + ( geom_traits, query_pt, s_j_prev, vertices[i] ); + while ( s.size() > 1 ) { + s_j = s_j_prev; + orient1 = orient2; + s.pop(); + s_j_prev = s.top(); + + orient2 = CGAL::Visibility_2::orientation_2 + ( geom_traits, query_pt, s_j_prev, vertices[i] ); + if ( orient1 != CGAL::LEFT_TURN && orient2 != CGAL::RIGHT_TURN ) { + mode = 1; + break; + } + + Segment_2 seg2( vertices[i-1], vertices[i] ); + Segment_2 seg( s_j_prev, s_j ); + if ( ( vertices[i-1] != s_j ) + && ( CGAL::Visibility_2::do_intersect_2 + + (geom_traits, seg, seg2) ) ) { + Object_2 result = CGAL::Visibility_2::intersect_2 + ( geom_traits, seg, seg2 ); + const Point_2 * ipoint = CGAL::object_cast(&result); + assert( ipoint != NULL ); + u = *ipoint; + mode = 2; + break; + } + } + + assert( mode != 0 ); + if ( mode == 1 ) { + orient1 = CGAL::Visibility_2::orientation_2 + ( geom_traits, query_pt, vertices[i], vertices[i+1] ); + orient2 = CGAL::Visibility_2::orientation_2 + ( geom_traits, vertices[i-1], vertices[i], vertices[i+1] ); + if ( orient1 == CGAL::RIGHT_TURN ) { + // Case R1 + // Since the next action is RIGHT, we do not compute the intersection + // of (s_j,s_j_prev) and the ray (query_pt, vertices[i]), + // thus, (s_j,s_j_prev) is not shortcutted, but it is harmless + upcase = RIGHT; + s.push( s_j ); + w = vertices[i]; + i++; + } else if ( orient2 == CGAL::RIGHT_TURN ) { + // Case R2 + Ray_2 ray( query_pt, vertices[i] ); + Segment_2 seg( s_j_prev, s_j ); + Object_2 result = CGAL::Visibility_2::intersect_2 + + ( geom_traits, seg, ray ); + const Point_2 * ipoint = CGAL::object_cast(&result); + assert( ipoint != NULL ); + u = *ipoint; + if ( s.top() != u ) { + s.push( u ); + } + upcase = LEFT; + s.push( vertices[i] ); + s.push( vertices[i+1] ); + w = vertices[i+1]; + i++; + } else { + // Case R3 + Ray_2 ray( query_pt, vertices[i] ); + Segment_2 seg( s_j_prev, s_j ); + Object_2 result = CGAL::Visibility_2::intersect_2 + + ( geom_traits, seg, ray ); + const Point_2 * ipoint = CGAL::object_cast(&result); + assert( ipoint != NULL ); + u = *ipoint; + if ( s.top() != u ) { + s.push( u ); + } + upcase = SCANC; + w = vertices[i]; + i++; + } + } else if ( mode == 2 ) { + // Case R4 + upcase = SCAND; + w = u; + } + } + + /*! Scans the vertices starting from index 'i' for the first visible vertex + out of the back hidden window */ + void scana(int& i, Point_2& w, const Point_2& query_pt) { + // Scan v_i, v_i+1, ..., v_n for the first edge to intersect (z, s_t) + Point_2 u; + int k = scan_edges( i, query_pt, s.top(), u, true ); + + CGAL::Orientation orient1 = CGAL::Visibility_2::orientation_2 + + ( geom_traits, + query_pt, + vertices[k], + vertices[k+1] ); + if ( orient1 == CGAL::RIGHT_TURN ) { + bool fwd = CGAL::Visibility_2::collinear_are_ordered_along_line_2 + + ( geom_traits, query_pt, s.top(), u ); + if ( !fwd ) { + // Case A1 + upcase = RIGHT; + i = k+1; + w = u; + } else { + // Case A2 + upcase = SCAND; + i = k+1; + w = u; + } + } else { + // Case A3 + upcase = LEFT; + i = k+1; + s.push( u ); + if ( u != vertices[k+1] ) { + s.push( vertices[k+1] ); + } + w = vertices[k+1]; + } + } + + /*! Find the first edge interecting the segment (v_0, s_t) */ + void scanb(int& i, Point_2& w, const Point_2& query_pt) { + if ( i == vertices.size() - 1 ) { + upcase = FINISH; + return; + } + Point_2 u; + int k = scan_edges( i, s.top(), vertices[0], u, false ); + if ( (k+1 == vertices.size()-1) && (vertices[0] == u) ) { + // Case B1 + upcase = FINISH; + s.push( vertices[0] ); + } else { + // Case B2 + upcase = RIGHT; + i = k+1; + w = u; + } + } + + /*! Finds the exit from a general front hidden window by finding the first + vertex to the right of the ray defined by the query_point and w*/ + void scanc(int& i, Point_2& w, const Point_2& query_pt) { + Point_2 u; + int k = scan_edges( i, s.top(), w, u, false ); + upcase = RIGHT; + i = k+1; + w = u; + } + + /*! find the first edge intersecting the given window (s_t, w) */ + void scand(int& i, Point_2& w, const Point_2& query_pt) { + Point_2 u; + int k = scan_edges( i, s.top(), w, u, false ); + upcase = LEFT; + i = k+1; + s.push( u ); + if ( u != vertices[k+1] ) { + s.push( vertices[k+1] ); + } + w = vertices[k+1]; + } + + /*! Scan edges v_i,v_{i+1},...,v_n, until find an edge intersecting given ray + or given segment. is_ray = true -> ray, false -> segment. + The intersection point is returned by u */ + int scan_edges( int i, const Point_2& ray_begin, const Point_2& ray_end, Point_2& u, bool is_ray ) { + CGAL::Orientation old_orient = CGAL::RIGHT_TURN; + Ray_2 ray( ray_begin, ray_end ); + Segment_2 s2( ray_begin, ray_end ); + int k; + Object_2 result; + for ( k = i; k+1 < vertices.size(); k++ ) { + CGAL::Orientation curr_orient = CGAL::Visibility_2::orientation_2 + ( geom_traits, + ray_begin, + ray_end, + vertices[k+1] ); + if ( curr_orient != old_orient ) { + // Orientation switch, an intersection may occur + Segment_2 seg( vertices[k], vertices[k+1] ); + if ( is_ray ) { + if (CGAL::Visibility_2::do_intersect_2 + + (geom_traits, seg, ray) ) { + result = CGAL::Visibility_2::intersect_2 + < Geometry_traits_2, Segment_2, Ray_2 > + ( geom_traits, seg, ray ); + break; + } + } else { + if (CGAL::Visibility_2::do_intersect_2 + + (geom_traits, seg, s2) ) { + result = CGAL::Visibility_2::intersect_2 + < Geometry_traits_2, Segment_2, Segment_2 > + ( geom_traits, seg, s2 ); + break; + } + } + } + old_orient = curr_orient; + } + assert( k+1( &result ); + if ( ipoint ) { + u = *ipoint; + } else { + u = vertices[k+1]; + } + return k; + } +}; + +} // namespace CGAL +#endif diff --git a/Visibility_2/include/CGAL/Triangular_expansion_visibility_2.h b/Visibility_2/include/CGAL/Triangular_expansion_visibility_2.h new file mode 100644 index 00000000000..b031b3f036a --- /dev/null +++ b/Visibility_2/include/CGAL/Triangular_expansion_visibility_2.h @@ -0,0 +1,579 @@ +// Copyright (c) 2013 Technical University Braunschweig (Germany). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Author(s): Michael Hemmer +// + +#ifndef CGAL_TRIANGULAR_EXPANSION_VISIBILITY_2__H +#define CGAL_TRIANGULAR_EXPANSION_VISIBILITY_2__H + +#include +#include +#include + +namespace CGAL { + +template +class Triangular_expansion_visibility_2 { + typedef typename Arrangement_2_::Geometry_traits_2 Geometry_traits_2; + typedef typename Geometry_traits_2::Kernel K; +public: + // Currently only consider with same type for both + typedef Arrangement_2_ Arrangement_2; + typedef typename Arrangement_2::Traits_2 Traits_2; + typedef typename Arrangement_2::Halfedge_const_handle Halfedge_const_handle; + typedef typename Arrangement_2::Halfedge_handle Halfedge_handle; + typedef typename Arrangement_2::Ccb_halfedge_const_circulator + Ccb_halfedge_const_circulator; + typedef typename Arrangement_2::Face_const_handle Face_const_handle; + typedef typename Arrangement_2::Face_handle Face_handle; + typedef typename Arrangement_2::Vertex_const_handle Vertex_const_handle; + typedef typename Arrangement_2::Vertex_handle Vertex_handle; + + typedef typename K::Point_2 Point_2; + typedef typename Geometry_traits_2::Ray_2 Ray_2; + typedef typename Geometry_traits_2::Segment_2 Segment_2; + typedef typename Geometry_traits_2::Line_2 Line_2; + typedef typename Geometry_traits_2::Vector_2 Vector_2; + typedef typename Geometry_traits_2::Direction_2 Direction_2; + typedef typename Geometry_traits_2::FT Number_type; + typedef typename Geometry_traits_2::Object_2 Object_2; + + // TODO + typedef RegularizationTag Regularization_tag; + + typedef CGAL::Tag_true Supports_general_polygon_tag; + typedef CGAL::Tag_true Supports_simple_polygon_tag; + +private: + typedef CGAL::Triangulation_vertex_base_2 Vb; + typedef CGAL::Constrained_triangulation_face_base_2 Fb; + typedef CGAL::Triangulation_data_structure_2 TDS; + typedef CGAL::No_intersection_tag Itag; + typedef CGAL::Constrained_Delaunay_triangulation_2 CDT; + +public: + const std::string name(){return std::string("T_visibility_2");} +private: + const Arrangement_2* p_arr; + boost::shared_ptr p_cdt; + std::vector needles; + +public: + Triangular_expansion_visibility_2() : p_arr(NULL){} + + /*! Constructor given an arrangement and the Regularization tag. */ + Triangular_expansion_visibility_2 (const Arrangement_2& arr) + : p_arr(&arr){ + init_cdt(); + } + + bool is_attached() { + //std::cout << "is_attached" << std::endl; + return (p_arr != NULL); + } + + void attach(const Arrangement_2& arr) { + // todo observe changes in arr; + if(p_arr != &arr){ + p_arr = &arr; + init_cdt(); + } + //std::cout << "attach done" << std::endl; + } + + void detach() { + //std::cout << "detach" << std::endl; + p_arr = NULL; + p_cdt = boost::shared_ptr(); + } + + const Arrangement_2& arr() { + return *p_arr; + } + + typename CDT::Edge get_edge(typename CDT::Face_handle fh, int i){ + return std::make_pair(fh,i); + } + + Point_2 ray_seg_intersection( + const Point_2& q, const Point_2& b, // the ray + const Point_2& s, const Point_2& t) // the segment + { + Ray_2 ray(q,b); + Segment_2 seg(s,t); + assert(typename K::Do_intersect_2()(ray,seg)); + CGAL::Object obj = typename K::Intersect_2()(ray,seg); + Point_2 result = object_cast(obj); + return result; + } + + void collect_needle( + const Point_2& q, + const typename CDT::Vertex_handle vh, + const typename CDT::Face_handle fh, + int index){ + + + // the expanded edge should not be constrained + assert(!p_cdt->is_constrained(get_edge(fh,index))); + assert(!p_cdt->is_infinite(fh)); + // go into the new face + const typename CDT::Face_handle nfh(fh->neighbor(index)); + assert(!p_cdt->is_infinite(nfh)); + + // get indices of neighbors + int nindex = nfh->index(fh); // index of new vertex and old face + int rindex = p_cdt->ccw(nindex); // index of face behind right edge + int lindex = p_cdt-> cw(nindex); // index of face behind left edge + + // get vertices seen from entering edge + const typename CDT::Vertex_handle nvh(nfh->vertex(nindex)); + const typename CDT::Vertex_handle rvh(nfh->vertex(p_cdt->cw (nindex))); + const typename CDT::Vertex_handle lvh(nfh->vertex(p_cdt->ccw(nindex))); + assert(!p_cdt->is_infinite(nvh)); + assert(!p_cdt->is_infinite(lvh)); + assert(!p_cdt->is_infinite(rvh)); + + // get edges seen from entering edge + typename CDT::Edge re = get_edge(nfh,p_cdt->ccw(nindex)); + typename CDT::Edge le = get_edge(nfh,p_cdt-> cw(nindex)); + + // do orientation computation once for new vertex + typename K::Orientation_2 orientation = + p_cdt->geom_traits().orientation_2_object(); + CGAL::Orientation orient = orientation(q,vh->point(),nvh->point()); + + + //std::cout << "\n collect_needle" <point() "<< vh->point() << std::endl; + //std::cout << "lvh->point() "<< lvh->point() << std::endl ; + //std::cout << "nvh->point() "<< nvh->point() << std::endl ; + //std::cout << "rvh->point() "<< rvh->point() << std::endl<< std::endl; + + + + switch ( orient ) { + case CGAL::COUNTERCLOCKWISE: + // looking on to the right edge + if(p_cdt->is_constrained(re)){ + if(vh!=rvh){ + Point_2 p = ray_seg_intersection(q,vh->point(),nvh->point(),rvh->point()); + //std::cout << vh->point() <<" -1- "<< p <point(),p)); + } + }else{ + collect_needle(q,vh,nfh,rindex); + } + break; + case CGAL::CLOCKWISE: + // looking on to the left edge + if(p_cdt->is_constrained(le)){ + if(vh!=lvh){ + Point_2 p = ray_seg_intersection(q,vh->point(),nvh->point(),lvh->point()); + //std::cout << vh->point() <<" -2- "<< p <point(),p)); + } + }else{ + collect_needle(q,vh,nfh,lindex); + } + break; + default: + assert(orient == CGAL::COLLINEAR); + // looking on nvh, so it must be reported + // if it wasn't already (triangles rotate around vh) + if(vh != nvh){ + //std::cout << vh->point() <<" -3- "<< nvh->point() <point(),nvh->point())); + } + // but we may also contiue looking along the vertex + if(!p_cdt->is_constrained(re)){ + collect_needle(q,nvh,nfh,rindex); + } + if(!p_cdt->is_constrained(le)){ + collect_needle(q,nvh,nfh,lindex); + } + break; + } + } + template + OIT expand_edge( + const Point_2& q, + const Point_2& left, + const Point_2& right, + typename CDT::Face_handle fh, + int index, + OIT oit){ + + // the expanded edge should not be constrained + assert(!p_cdt->is_constrained(get_edge(fh,index))); + assert(!p_cdt->is_infinite(fh)); + + // go into the new face + const typename CDT::Face_handle nfh(fh->neighbor(index)); + assert(!p_cdt->is_infinite(nfh)); + + // get indices of neighbors + int nindex = nfh->index(fh); // index of new vertex and old face + int rindex = p_cdt->ccw(nindex); // index of face behind right edge + int lindex = p_cdt-> cw(nindex); // index of face behind left edge + + // get vertices seen from entering edge + const typename CDT::Vertex_handle nvh(nfh->vertex(nindex)); + const typename CDT::Vertex_handle rvh(nfh->vertex(p_cdt->cw (nindex))); + const typename CDT::Vertex_handle lvh(nfh->vertex(p_cdt->ccw(nindex))); + assert(!p_cdt->is_infinite(nvh)); + assert(!p_cdt->is_infinite(lvh)); + assert(!p_cdt->is_infinite(rvh)); + + // get edges seen from entering edge + typename CDT::Edge re = get_edge(nfh,p_cdt->ccw(nindex)); + typename CDT::Edge le = get_edge(nfh,p_cdt-> cw(nindex)); + + // do orientation computation once for new vertex + typename K::Orientation_2 orientation = + p_cdt->geom_traits().orientation_2_object(); + CGAL::Orientation ro = orientation(q,right,nvh->point()); + CGAL::Orientation lo = orientation(q,left ,nvh->point()); + + assert(typename K::Orientation_2()(q,left ,lvh->point()) != CGAL::CLOCKWISE); + assert(typename K::Orientation_2()(q,right,rvh->point()) != CGAL::COUNTERCLOCKWISE); + + //std::cout << (ro == CGAL::COUNTERCLOCKWISE) << " " << (lo == CGAL::CLOCKWISE) << std::endl; + + + + //right edge is seen if new vertex is counter clockwise of right boarder + if(ro == CGAL::COUNTERCLOCKWISE){ + if(p_cdt->is_constrained(re)){ + // the edge is constrained + // report intersection with right boarder ray + // if it is not already the right vertex (already reported) + if(right != rvh->point()){ + *oit++ = ray_seg_intersection(q,right,nvh->point(),rvh->point()); + } + + // then report intersection with left boarder if it exists + if(lo == CGAL::COUNTERCLOCKWISE){ + *oit++ = ray_seg_intersection(q,left,nvh->point(),rvh->point()); + } + }else{ + // the edge is not a constrained + if(lo == CGAL::COUNTERCLOCKWISE){ + // no split needed and return + //std::cout<< "h1"<< std::endl; + oit = expand_edge(q,left,right,nfh,rindex,oit); + //std::cout<< "h1 done"<< std::endl; + return oit; + }else{ + // spliting at new vertex + //std::cout<< "h2"<< std::endl; + *oit++ = expand_edge(q,nvh->point(),right,nfh,rindex,oit); + //std::cout<< "h2 done"<< std::endl; + } + } + } + + + //std::cout << "q "<< q << std::endl ; + //std::cout << "lvh->point() "<< lvh->point() << std::endl; + //std::cout << "left "<< left << std::endl ; + //std::cout << "nvh->point() "<< nvh->point() << std::endl ; + //std::cout << "right "<< right << std::endl ; + //std::cout << "rvh->point() "<< rvh->point() << std::endl<< std::endl; + + + // determin whether new vertex needs to be reported + if(ro != CGAL::CLOCKWISE && lo != CGAL::COUNTERCLOCKWISE){ + *oit++ = nvh->point(); + } + if(!Regularization_tag::value){ + assert(!(ro == CGAL::COLLINEAR && lo == CGAL::COLLINEAR)); + // we have to check whether a needle starts here. + if(p_cdt->is_constrained(le) && !p_cdt->is_constrained(re) && ro == CGAL::COLLINEAR) + collect_needle(q,nvh,nfh,rindex); + if(p_cdt->is_constrained(re) && !p_cdt->is_constrained(le) && lo == CGAL::COLLINEAR) + collect_needle(q,nvh,nfh,lindex); + } + + //left edge is seen if new vertex is clockwise of left boarder + if(lo == CGAL::CLOCKWISE){ + if(p_cdt->is_constrained(le)){ + // the edge is constrained + // report interesection with right boarder if exists + if(ro == CGAL::CLOCKWISE){ + *oit++ = ray_seg_intersection(q,right,nvh->point(),lvh->point()); + } + // then report intersection with left boarder ray + // if it is not already the left vertex (already reported) + if(left != lvh->point()){ + *oit++ = ray_seg_intersection(q,left,nvh->point(),lvh->point()); + } + return oit; + }else{ + // the edge is not a constrained + if(ro == CGAL::CLOCKWISE){ + // no split needed and return + //std::cout<< "h3"<< std::endl; + oit = expand_edge(q,left,right,nfh,lindex,oit); + //std::cout<< "h3 done"<< std::endl; + return oit; + }else{ + // spliting at new vertex + //std::cout<< "h4"<< std::endl; + oit = expand_edge(q,left,nvh->point(),nfh,lindex,oit); + //std::cout<< "h4 done"<< std::endl; + return oit; + } + } + } + + return oit; + } + + template + typename VARR::Face_handle + compute_visibility(const Point_2& q, + const Face_const_handle face, + VARR& out_arr + ){ + //std::cout << "query in face interior" << std::endl; + + out_arr.clear(); + needles.clear(); + assert(!face->is_unbounded()); + + + std::vector raw_output; + typename CDT::Face_handle fh = p_cdt->locate(q); + + raw_output.push_back(fh->vertex(1)->point()); + if(!p_cdt->is_constrained(get_edge(fh,0))){ + //std::cout<< "edge 0 is not constrained" << std::endl; + expand_edge( + q, + fh->vertex(2)->point(), + fh->vertex(1)->point(), + fh,0,std::back_inserter(raw_output)); + } + + raw_output.push_back(fh->vertex(2)->point()); + if(!p_cdt->is_constrained(get_edge(fh,1))){ + //std::cout << "edge 1 is not constrained" << std::endl; + expand_edge( + q, + fh->vertex(0)->point(), + fh->vertex(2)->point(), + fh,1,std::back_inserter(raw_output)); + } + + raw_output.push_back(fh->vertex(0)->point()); + if(!p_cdt->is_constrained(get_edge(fh,2))){ + //std::cout << "edge 2 is not constrained" << std::endl; + expand_edge( + q, + fh->vertex(1)->point(), + fh->vertex(0)->point(), + fh,2,std::back_inserter(raw_output)); + } + + + return output(raw_output,out_arr); + } + + template + typename VARR::Face_handle + compute_visibility(const Point_2& q, + const Halfedge_const_handle he, + VARR& out_arr) { + //std::cout << "visibility_region he" << std::endl; + + assert(!he->face()->is_unbounded()); + out_arr.clear(); + needles.clear(); + + std::vector raw_output; + typename CDT::Locate_type location; + int index; + typename CDT::Face_handle fh = p_cdt->locate(q,location,index); + assert(location == CDT::EDGE || location == CDT::VERTEX); + // the following code tries to figure out which triangle one should start in. + + + if(location == CDT::EDGE){ + //std::cout << "query on edge" << std::endl; + // this is the easy part, there are only two possible faces + // index indicates the edge = vertex on the other side of the edge + // the next vertex in cw order should be the target of given edge + if(fh->vertex(p_cdt->cw(index))->point() != he->target()->point()){ + //std::cout << "need to swap face" << std::endl; + // take face on the other side if this is not the case + typename CDT::Face_handle nfh = fh->neighbor(index); + index = nfh->index(fh); + fh = nfh; + } + assert(fh->vertex(p_cdt->cw(index))->point() == he->target()->point()); + assert(!p_cdt->is_infinite(fh->vertex(index))); + + + // output the edge the query lies on + raw_output.push_back(he->source()->point()); + raw_output.push_back(he->target()->point()); + + if(!p_cdt->is_constrained(get_edge(fh,p_cdt->ccw(index)))){ + expand_edge( + q, + fh->vertex(index)->point(), //left + he->target()->point() , //right + fh, + p_cdt->ccw(index), + std::back_inserter(raw_output)); + } + raw_output.push_back(fh->vertex(index)->point()); + + if(!p_cdt->is_constrained(get_edge(fh,p_cdt->cw(index)))){ + expand_edge( + q, + he->source()->point() , //left + fh->vertex(index)->point(), //right + fh, + p_cdt->cw(index), + std::back_inserter(raw_output)); + } + } + + if(location == CDT::VERTEX){ + //std::cout << "query on vertex" << std::endl; + + //bool query_point_on_vertex_is_not_working_yet = false; + //assert(query_point_on_vertex_is_not_working_yet); + + assert(q == he->target()->point()); + assert(fh->vertex(index)->point() == he->target()->point()); + + // push points that are seen anyway + // raw_output.push_back(he->source()->point()); inserted last + raw_output.push_back(he->target()->point()); + raw_output.push_back(he->next()->target()->point()); + + // now start in the triangle that contains he->next() + while( + p_cdt->is_infinite(fh->vertex(p_cdt->ccw(index))) || + he->next()->target()->point() != fh->vertex(p_cdt->ccw(index))->point()){ + typename CDT::Face_handle nfh = fh->neighbor(p_cdt->ccw(index)); + int nindex = nfh->index(fh); + index = p_cdt->ccw(nindex); + fh = nfh; + assert(he->target()->point() == fh->vertex(index)->point()); + } + + + + + assert(he->next()->source()->point() == fh->vertex(index)->point()); + assert(he->next()->target()->point() == fh->vertex(p_cdt->ccw(index))->point()); + assert(!p_cdt->is_infinite(fh)); + assert(p_cdt->is_constrained(get_edge(fh,p_cdt->cw(index)))); + + while(he->source()->point() != fh->vertex(p_cdt->ccw(index))->point()){ + + if(!p_cdt->is_constrained(get_edge(fh,index))){ + expand_edge( + q, + fh->vertex(p_cdt-> cw(index))->point(), //left + fh->vertex(p_cdt->ccw(index))->point(), //right + fh, + index, + std::back_inserter(raw_output)); + } + // push left end point of edge into output + raw_output.push_back(fh->vertex(p_cdt-> cw(index))->point()); + + // take the next triangle around q in ccw order + typename CDT::Face_handle nfh = fh->neighbor(p_cdt->ccw(index)); + int nindex = nfh->index(fh); + index = p_cdt->ccw(nindex); + fh = nfh; + assert(fh->vertex(index)->point() == he->target()->point()); + } + } + return output(raw_output,out_arr); + } + + template + typename VARR::Face_handle + output(std::vector& raw_output, VARR& out_arr){ + + if(needles.size()>0){ + std::vector segments(needles.begin(),needles.end()); + for(unsigned int i = 0; i target(); + }else{ + v_last = out_arr.insert_from_right_vertex(Segment_2(raw_output[i], raw_output[i+1]), v_last)->target(); + } + } + out_arr.insert_at_vertices(Segment_2(raw_output.front(),raw_output.back()),v_last,v_first); + } + + assert(out_arr.number_of_faces()== 2); + + if(out_arr.faces_begin()->is_unbounded()) + return ++out_arr.faces_begin(); + else + return out_arr.faces_begin(); + } + + void init_cdt(){ + //std::cout<< "==============" < > constraints; + for(typename Arrangement_2::Edge_const_iterator eit = p_arr->edges_begin(); + eit != p_arr->edges_end(); eit++){ + Point_2 source = eit->source()->point(); + Point_2 target = eit->target()->point(); + //std::cout << source << " -- " << target << std::endl; + constraints.push_back(std::make_pair(source,target)); + } + //std::cout << "init_cdt new CDT" << std::endl; + p_cdt = boost::shared_ptr(new CDT(constraints.begin(),constraints.end())); + //std::cout << "init_cdt done" << std::endl; + //std::cout << std::endl; + } +}; + +} // namespace CGAL + +#endif //CGAL_TRIANGULAR_EXPANSION_VISIBILITY_2__H diff --git a/Visibility_2/include/CGAL/Visibility_2/visibility_utils.h b/Visibility_2/include/CGAL/Visibility_2/visibility_utils.h new file mode 100644 index 00000000000..90f092522f5 --- /dev/null +++ b/Visibility_2/include/CGAL/Visibility_2/visibility_utils.h @@ -0,0 +1,371 @@ +// Copyright (c) 2013 Technical University Braunschweig (Germany). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Author(s): Francisc Bungiu +// Michael Hemmer + +#ifndef CGAL_VISIBILITY_UTILS_H +#define CGAL_VISIBILITY_UTILS_H + +#include +#include +#include + +namespace CGAL { +namespace Visibility_2 { + + template +int count_edges_in_face(typename Arrangement_2::Face_const_handle fch) { + typedef typename Arrangement_2::Ccb_halfedge_const_circulator + Ccb_halfedge_const_circulator; + + Ccb_halfedge_const_circulator circ = fch->outer_ccb(); + Ccb_halfedge_const_circulator curr = circ; + + int edge_cnt(0); + do { + edge_cnt++; + } while (++curr != circ); + return edge_cnt; +} + +template +void print_edge(Edge_const_iterator eit) { + std::cout << "[" << eit->curve() << "]" << std::endl; +} +template +void print_simple_face(Face_const_handle fh) { + Ccb_halfedge_const_circulator cir = fh->outer_ccb(); + Ccb_halfedge_const_circulator curr = cir; + do { + std::cout << "[" << curr->curve() << "]" << std::endl; + } while (++ curr != cir); +} + +template +void print_arrangement(const Arrangement_2& arr) { + typedef typename Arrangement_2::Edge_const_iterator Edge_const_iterator; + Edge_const_iterator eit; + std::cout << arr.number_of_edges() << " edges:" << std::endl; + for (eit = arr.edges_begin(); eit != arr.edges_end(); ++eit) + print_edge(eit); +} + +template +void print_arrangement_by_face(const Arrangement_2& arr) { + typedef typename Arrangement_2::Face_const_iterator Face_const_iterator; + typedef typename Arrangement_2::Ccb_halfedge_const_circulator + Ccb_halfedge_const_circulator; + Face_const_iterator fit; + for (fit = arr.faces_begin() ; fit != arr.faces_end() ; fit++) { + if (!fit->is_unbounded()) { + std::cout << "FACE\n"; + print_simple_face(fit); + std::cout << "END FACE\n"; + } + } +} + +template +Orientation orientation_2(const Geometry_traits_2 *geom_traits, + const typename Geometry_traits_2::Point_2& p, + const typename Geometry_traits_2::Point_2& q, + const typename Geometry_traits_2::Point_2& r) { + + typename Geometry_traits_2::Orientation_2 orient = + geom_traits->orientation_2_object(); + return orient(p, q, r); +} + +template +bool less_distance_to_point_2(const Geometry_traits_2 *geom_traits, + const typename Geometry_traits_2::Point_2& p, + const typename Geometry_traits_2::Point_2& q, + const typename Geometry_traits_2::Point_2& r) { + + typename Geometry_traits_2::Less_distance_to_point_2 less_dist = + geom_traits->less_distance_to_point_2_object(); + return less_dist(p, q, r); +} + +template +bool collinear(const Geometry_traits_2 *geom_traits, + const typename Geometry_traits_2::Point_2& p, + const typename Geometry_traits_2::Point_2& q, + const typename Geometry_traits_2::Point_2& r) { + + typename Geometry_traits_2::Collinear_2 collinear_fnct = + geom_traits->collinear_2_object(); + return collinear_fnct(p, q, r); +} + +template +typename Geometry_traits_2::Object_2 intersect_2(const Geometry_traits_2 *geom_traits, + const _Curve_first& s1, + const _Curve_second& s2) { + + typedef typename Geometry_traits_2::Kernel Kernel; + const Kernel *kernel = static_cast (geom_traits); + typename Kernel::Intersect_2 intersect_fnct = kernel->intersect_2_object(); + return intersect_fnct(s1, s2); +} + +template +CGAL::Comparison_result compare_xy_2(const Geometry_traits_2 *geom_traits, + const typename Geometry_traits_2::Point_2 &p, + const typename Geometry_traits_2::Point_2 &q) { + + typename Geometry_traits_2::Compare_xy_2 cmp = + geom_traits->compare_xy_2_object(); + return cmp(p, q); +} + +template +typename Geometry_traits_2::FT compute_squared_distance_2( + const Geometry_traits_2 *geom_traits, + const Type1& p, + const Type2& seg) { + + typename Geometry_traits_2::Compute_squared_distance_2 compute_dist = + geom_traits->compute_squared_distance_2_object(); + return compute_dist(p, seg); +} + +template +bool do_intersect_2(const Geometry_traits_2 *geom_traits, + const Type1& c1, + const Type2& c2) { + + typename Geometry_traits_2::Do_intersect_2 intersect = + geom_traits->do_intersect_2_object(); + return intersect(c1, c2); +} + +template +bool collinear_are_ordered_along_line_2(const Geometry_traits_2 *geom_traits, + const typename Geometry_traits_2::Point_2 &p, + const typename Geometry_traits_2::Point_2 &q, + const typename Geometry_traits_2::Point_2 &r) { + + typename Geometry_traits_2::Collinear_are_ordered_along_line_2 coll = + geom_traits->collinear_are_ordered_along_line_2_object(); + return coll(p, q, r); +} + +template +typename Geometry_traits_2::Point_2 construct_projected_point_2( + const Geometry_traits_2 *geom_traits, + const Type1& s, + const Type2& p) { + + typedef typename Geometry_traits_2::Point_2 Point_2; + typedef typename Geometry_traits_2::FT Number_type; + typename Geometry_traits_2::Construct_projected_point_2 construct_proj = + geom_traits->construct_projected_point_2_object(); + Point_2 proj_pt = construct_proj(s.supporting_line(), p); + if (s.has_on(proj_pt)) { + return proj_pt; + } + else { + Number_type d_to_src = compute_squared_distance_2 + (geom_traits, proj_pt, s.source()); + Number_type d_to_trg = compute_squared_distance_2 + (geom_traits, proj_pt, s.target()); + if (d_to_src < d_to_trg) { + return s.source(); + } + else { + return s.target(); + } + } +} + +//construct an arrangement of visibility region from a vector of circular ordered vertices with respect to the query point +template +void report_while_handling_needles( + const typename Visibility_2::Arrangement_2::Geometry_traits_2 *geom_traits, + const typename Visibility_2::Arrangement_2::Point_2& q, + std::vector& points, + Visibility_arrangement_2& arr_out) { + + typedef typename Visibility_2::Arrangement_2 Arrangement_2; + typedef typename Arrangement_2::Point_2 Point_2; + typedef typename Arrangement_2::Geometry_traits_2 Geometry_traits_2; + typedef typename Arrangement_2::Halfedge_handle Halfedge_handle; + typedef typename Arrangement_2::Vertex_handle Vertex_handle; + typedef typename Geometry_traits_2::Segment_2 Segment_2; + typedef typename Geometry_traits_2::Direction_2 Direction_2; + + typename std::vector::size_type i = 0; + + if (points.front() == points.back()) { + points.pop_back(); + } + + while (collinear(geom_traits, + q, + points[i], + points.back())) { + + points.push_back(points[i]); + i++; + } + + points.push_back(points[i]); + + typename Visibility_arrangement_2::Halfedge_handle he_handle; + typename Visibility_arrangement_2::Vertex_handle v_trg; //the handle of vertex where the next segment is inserted + typename Visibility_arrangement_2::Vertex_handle v_fst; //the handle of vertex inserted first + typename Visibility_arrangement_2::Vertex_handle v_needle_end; //the handle of vertex of the end of a needle + + v_trg = v_fst = arr_out.insert_in_face_interior(points[i], arr_out.unbounded_face()); + //find a point that is right after a needle + while (i+1 < points.size()) { + if ( collinear(geom_traits, + points[i], + points[i+1], + q)) { + typename Visibility_arrangement_2::Vertex_handle v_needle_begin = v_trg; + std::vector forward_needle; //vertices of the needle that are not between q and v_needle_begin; their direction is leaving q; + std::vector backward_needle; //vertices of the needle that are not between q and v_needle_begin; their direction is towards q; + std::vector part_in_q_side; //vertices of the needle that are between q and v_needle_begin + part_in_q_side.push_back(points[i]); + forward_needle.push_back((points[i])); + + bool same_side_of_q = (compare_xy_2(geom_traits, points[i], q)==compare_xy_2(geom_traits, points[i], points[i+1])); + if (same_side_of_q) + part_in_q_side.push_back(points[i+1]); + else + forward_needle.push_back(points[i+1]); + i++; + while (i+1< points.size() && orientation_2(geom_traits, + points[i], + points[i+1], + q ) == CGAL::COLLINEAR) { + if (same_side_of_q) { + part_in_q_side.push_back(points[i+1]); + } + else { + if (compare_xy_2(geom_traits, part_in_q_side.front(), q) + == compare_xy_2(geom_traits, part_in_q_side.front(), points[i+1])) { + same_side_of_q = true; + part_in_q_side.push_back(points[i+1]); + } + else { + if (less_distance_to_point_2(geom_traits, q, points[i], points[i+1])) + forward_needle.push_back(points[i+1]); + else + backward_needle.push_back(points[i+1]); + } + } + i++; + } + //obtain the end point of a needle + Point_2 end_of_needle; + if (same_side_of_q) + end_of_needle = part_in_q_side.back(); + else { + if (backward_needle.empty()) { + end_of_needle = forward_needle.back(); + } + else { + end_of_needle = backward_needle.back(); + } + } + + std::reverse(backward_needle.begin(), backward_needle.end()); + std::vector merged_needle; + + // merge the forward_needle and backward_needle + unsigned int itr_fst = 0, itr_snd = 0; + while (itr_fst < forward_needle.size() && + itr_snd < backward_needle.size()) { + + if (less_distance_to_point_2(geom_traits, + q, + forward_needle[itr_fst], + backward_needle[itr_snd])) { + + merged_needle.push_back(forward_needle[itr_fst]); + itr_fst++; + } + else { + merged_needle.push_back(backward_needle[itr_snd]); + itr_snd++; + } + } + while (itr_fst < forward_needle.size()) { + merged_needle.push_back(forward_needle[itr_fst]); + itr_fst++; + } + while (itr_snd < backward_needle.size()) { + merged_needle.push_back(backward_needle[itr_snd]); + itr_snd++; + } + + for (unsigned int p = 0 ; p+1 < merged_needle.size() ; p++) { + if (CGAL::Visibility_2::compare_xy_2(geom_traits, merged_needle[p], merged_needle[p+1]) == CGAL::SMALLER) { + he_handle = arr_out.insert_from_left_vertex(Segment_2(merged_needle[p], merged_needle[p+1]), v_trg); + } + else { + he_handle = arr_out.insert_from_right_vertex(Segment_2(merged_needle[p], merged_needle[p+1]), v_trg); + } + v_trg = he_handle->target(); + if (merged_needle[p+1] == end_of_needle) { + v_needle_end = v_trg; + } + } + if (same_side_of_q) { + //insert the part of needle between q and v_needle_begin + v_trg = v_needle_begin; + for (unsigned int p = 0 ; p+1 < part_in_q_side.size() ; p++) { + if (CGAL::Visibility_2::compare_xy_2(geom_traits, part_in_q_side[p], part_in_q_side[p+1]) == CGAL::SMALLER) { + he_handle = arr_out.insert_from_left_vertex(Segment_2(part_in_q_side[p], part_in_q_side[p+1]), v_trg); + } + else { + he_handle = arr_out.insert_from_right_vertex(Segment_2(part_in_q_side[p], part_in_q_side[p+1]), v_trg); + } + v_trg = he_handle->target(); + } + } + else + v_trg = v_needle_end; + } + else { + if (CGAL::Visibility_2::compare_xy_2(geom_traits, v_trg->point(), points[i+1]) == CGAL::SMALLER) { + he_handle = arr_out.insert_from_left_vertex(Segment_2(points[i], points[i+1]), v_trg); + } + else { + he_handle = arr_out.insert_from_right_vertex(Segment_2(points[i], points[i+1]), v_trg); + } + v_trg = he_handle->target(); + i++; + } + if (i+2 == points.size()) { + //close the boundary + v_trg = he_handle->target(); + arr_out.insert_at_vertices(Segment_2(points[points.size()-2], points[points.size()-1]), v_trg, v_fst); + break; + } + } +} + +} // end namespace Visibility_2 +} // end namespace CGAL + +#endif