1 // Copyright (c) 1999 2 // Max-Planck-Institute Saarbruecken (Germany). 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/Point_set_2/include/CGAL/Point_set_2.h $ 7 // $Id: Point_set_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) : Matthias Baesken 12 13 #ifndef CGAL_POINT_SET_2_H 14 #define CGAL_POINT_SET_2_H 15 16 #include <CGAL/license/Point_set_2.h> 17 18 19 #include <CGAL/basic.h> 20 #include <CGAL/Unique_hash_map.h> 21 #include <CGAL/Delaunay_triangulation_2.h> 22 #include <CGAL/squared_distance_2_1.h> 23 #include <CGAL/compare_vertices.h> 24 #include <list> 25 #include <queue> 26 #include <map> 27 #include <stack> 28 #include <cmath> 29 #include <climits> 30 31 namespace CGAL { 32 33 template<class Gt, class Tds = Triangulation_data_structure_2 <Triangulation_vertex_base_2<Gt> > > 34 class Point_set_2 : public Delaunay_triangulation_2<Gt,Tds> 35 { 36 37 public: 38 typedef Gt Geom_traits; 39 40 typedef typename Geom_traits::Point_2 Point; 41 typedef typename Geom_traits::Segment_2 Segment; 42 43 typedef typename Geom_traits::Circle_2 Circle; 44 45 typedef typename Geom_traits::Orientation_2 Orientation_2; 46 typedef typename Geom_traits::Side_of_oriented_circle_2 Side_of_oriented_circle_2; 47 typedef typename Geom_traits::Construct_circle_2 Construct_circle_2; 48 typedef typename Geom_traits::Compute_squared_distance_2 Compute_squared_distance_2; 49 typedef typename Geom_traits::FT Numb_type; // field number type ... 50 51 52 typedef Delaunay_triangulation_2<Gt,Tds> Triangulation; 53 typedef typename Triangulation::size_type size_type; 54 typedef typename Triangulation::Locate_type Locate_type; 55 typedef typename Triangulation::Face_handle Face_handle; 56 typedef typename Triangulation::Vertex_handle Vertex_handle; 57 typedef typename Triangulation::Edge Edge; 58 typedef typename Triangulation::Vertex Vertex; 59 typedef typename Triangulation::Face Face; 60 typedef typename Triangulation::Edge_circulator Edge_circulator; 61 typedef typename Triangulation::Finite_edges_iterator Finite_edges_iterator; 62 typedef typename Triangulation::Vertex_iterator Vertex_iterator; 63 typedef typename Triangulation::Vertex_circulator Vertex_circulator; 64 typedef typename Triangulation::Edge_iterator Edge_iterator; 65 66 typedef typename Geom_traits::Bounded_side_2 Circleptori; 67 typedef typename Geom_traits::Compare_distance_2 Comparedist; 68 typedef typename Geom_traits::Construct_center_2 Circlecenter; 69 70 typedef Unique_hash_map<Vertex_handle, Numb_type> MAP_TYPE; 71 typedef Delaunay_triangulation_2<Gt,Tds> Base; 72 73 using Base::finite_vertices_begin; 74 using Base::finite_vertices_end; 75 using Base::number_of_vertices; 76 using Base::VERTEX; 77 using Base::insert; 78 using Base::remove; 79 using Base::locate; 80 using Base::is_infinite; 81 using Base::nearest_vertex; 82 using Base::incident_vertices; 83 84 Comparedist tr_comparedist; 85 Orientation_2 tr_orientation; 86 Side_of_oriented_circle_2 tr_so_circle; 87 Compute_squared_distance_2 tr_sqrdist; 88 Circleptori tr_circleptori; 89 Circlecenter tr_circlecenter; 90 91 //constructions... 92 Construct_circle_2 tr_createcircle_3p; 93 Point_set_2()94 Point_set_2() 95 { 96 init_vertex_marks(); 97 } 98 99 template<class InputIterator> Point_set_2(InputIterator first,InputIterator last)100 Point_set_2(InputIterator first, InputIterator last) 101 { 102 init_vertex_marks(); 103 insert(first,last); 104 } 105 ~Point_set_2()106 ~Point_set_2() {} 107 108 template<class OutputIterator> vertices(OutputIterator res)109 OutputIterator vertices(OutputIterator res) 110 // return vertex handles ... 111 { 112 Vertex_iterator vit = finite_vertices_begin(); 113 for (; vit != finite_vertices_end(); vit++) { *res= vit; res++; } 114 return res; 115 } 116 117 lookup(Point p)118 Vertex_handle lookup(Point p) const 119 { 120 if (number_of_vertices() == 0) return nullptr; 121 122 // locate ... 123 Locate_type lt; 124 int li; 125 Face_handle fh = locate(p,lt,li); 126 127 if (lt == VERTEX){ 128 Face f = *fh; 129 return f.vertex(li); 130 } 131 else return nullptr; 132 } 133 134 nearest_neighbor(Point p)135 Vertex_handle nearest_neighbor(Point p) 136 { 137 if (number_of_vertices() == 0) return nullptr; 138 return nearest_vertex(p); 139 } 140 141 nearest_neighbor(Vertex_handle v)142 Vertex_handle nearest_neighbor(Vertex_handle v) const 143 { 144 if (number_of_vertices() <= 1) return nullptr; 145 Point p = v->point(); 146 147 Vertex_circulator vc = incident_vertices(v); 148 Vertex_circulator start =vc; 149 150 Vertex_handle min_v = vc; 151 if (is_infinite(min_v)){ 152 vc++; 153 min_v = vc; 154 } 155 156 Vertex_handle act; 157 158 // go through the vertices ... 159 do { 160 act = vc; 161 162 if (! is_infinite(act)) { 163 if ( tr_comparedist(p,act->point(), min_v->point()) == SMALLER ) { 164 min_v = act; 165 } 166 } 167 168 vc++; 169 } while (vc != start); 170 171 return min_v; 172 } 173 174 template<class OutputIterator> nearest_neighbors(Point p,size_type k,OutputIterator res)175 OutputIterator nearest_neighbors(Point p, size_type k, OutputIterator res) 176 { 177 size_type n = number_of_vertices(); 178 179 if ( k <= 0 ) return res; 180 if ( n <= k ) { // return all finite vertices ... 181 return vertices(res); 182 } 183 184 // insert p, if nesessary 185 186 Vertex_handle vh = lookup(p); 187 bool old_node = true; 188 189 // we have to add a new vertex ... 190 if (vh == nullptr){ 191 vh = insert(p); 192 old_node = false; 193 k++; 194 } 195 196 std::list<Vertex_handle> res_list; 197 nearest_neighbors_list(vh, k, res_list); 198 199 if ( !old_node ) 200 { 201 res_list.pop_front(); 202 remove(vh); 203 } 204 205 typename std::list<Vertex_handle>::iterator it = res_list.begin(); 206 207 for (; it != res_list.end(); it++) { *res= *it; res++; } 208 209 return res; 210 } 211 212 template<class OutputIterator> nearest_neighbors(Vertex_handle v,size_type k,OutputIterator res)213 OutputIterator nearest_neighbors(Vertex_handle v, size_type k,OutputIterator res) 214 { 215 size_type n = number_of_vertices(); 216 217 if ( k <= 0 ) return res; 218 if ( n <= k ) { // return all (finite) vertices ... 219 return vertices(res); 220 } 221 222 std::list<Vertex_handle> res_list; 223 nearest_neighbors_list(v, k, res_list); 224 225 typename std::list<Vertex_handle>::iterator it = res_list.begin(); 226 227 for (; it != res_list.end(); it++) { *res= *it; res++; } 228 229 return res; 230 } 231 232 nearest_neighbors_list(Vertex_handle v,size_type k,std::list<Vertex_handle> & res)233 void nearest_neighbors_list(Vertex_handle v, 234 size_type k, 235 std::list<Vertex_handle>& res) 236 { 237 size_type n = number_of_vertices(); 238 239 if ( k <= 0 ) return; 240 if ( n <= k ) { vertices(std::back_inserter(res)); return; } 241 242 Point p = v->point(); 243 244 // "unmark" the vertices ... 245 init_search(); 246 247 MAP_TYPE priority_number; // here we save the priorities ... 248 internal::compare_vertices<Vertex_handle,Numb_type,MAP_TYPE> 249 comp(& priority_number); // comparison object ... 250 std::priority_queue<Vertex_handle, std::vector<Vertex_handle>, internal::compare_vertices<Vertex_handle,Numb_type,MAP_TYPE> > PQ(comp); 251 252 priority_number[v] = 0; 253 PQ.push(v); 254 255 mark_vertex(v); 256 257 while ( k > 0 ) 258 { 259 // find minimum from PQ ... 260 Vertex_handle w = PQ.top(); 261 PQ.pop(); 262 263 res.push_back(w); 264 k--; 265 266 // get the incident vertices of w ... 267 Vertex_circulator vc = incident_vertices(w); 268 Vertex_circulator start =vc; 269 Vertex_handle act; 270 271 do { 272 act = vc; 273 274 if ( (!is_marked(act)) && (! is_infinite(act)) ) 275 { 276 priority_number[act] = tr_sqrdist(p,act->point()); 277 PQ.push(act); 278 mark_vertex(act); 279 } 280 281 vc++; 282 } while (vc != start); 283 284 } 285 } 286 287 288 // for marking nodes in search procedures 289 size_type cur_mark; 290 291 Unique_hash_map<Vertex_handle, size_type> mark; 292 init_vertex_marks()293 void init_vertex_marks() 294 { 295 cur_mark = 0; 296 mark.clear(); 297 } 298 init_search()299 void init_search() 300 { 301 cur_mark++; 302 if (cur_mark == (std::numeric_limits<size_type>::max)()) init_vertex_marks(); 303 } 304 mark_vertex(Vertex_handle vh)305 void mark_vertex(Vertex_handle vh) 306 // mark vh as visited ... 307 { 308 mark[vh] = cur_mark; 309 } 310 is_marked(Vertex_handle vh)311 bool is_marked(Vertex_handle vh) 312 { 313 if (! mark.is_defined(vh)) return false; 314 315 return (mark[vh] == cur_mark); 316 } 317 search(Vertex_handle v,const Circle & C,std::list<Vertex_handle> & L)318 void search(Vertex_handle v,const Circle& C, std::list<Vertex_handle>& L) 319 { 320 std::stack<Vertex_handle> todo; 321 todo.push(v); 322 323 while (!todo.empty()) 324 { 325 Vertex_handle current = todo.top(); 326 todo.pop(); 327 328 if (is_marked(current)) 329 continue; 330 331 L.push_back(current); 332 mark_vertex(current); 333 334 // get incident vertices of v ... 335 Vertex_circulator vc = incident_vertices(current); 336 Vertex_circulator start =vc; 337 Vertex_handle act; 338 // go through the vertices ... 339 do { 340 act = vc; 341 342 if (! is_infinite(act)) { 343 if (!is_marked(act) && ! (tr_circleptori(C,act->point())==ON_UNBOUNDED_SIDE) ) 344 todo.push(act); 345 } 346 vc++; 347 } while (vc != start); 348 } 349 } 350 351 352 template<class OutputIterator> range_search(const Circle & C,OutputIterator res)353 OutputIterator range_search(const Circle& C, OutputIterator res) 354 { 355 if (number_of_vertices() == 0) return res; 356 if (number_of_vertices() == 1) 357 { 358 // get the one vertex ... 359 Vertex_iterator vit = finite_vertices_begin(); 360 Point p = vit->point(); 361 362 if (! (tr_circleptori(C, p) == ON_UNBOUNDED_SIDE)){ 363 *res= vit; res++; 364 } 365 return res; 366 } 367 368 // normal case ... 369 Point p = tr_circlecenter(C); 370 Vertex_handle v = lookup(p); 371 bool new_v = false; 372 373 if ( v == nullptr ) 374 { 375 new_v = true; 376 v = insert(p); 377 } 378 379 init_search(); 380 381 std::list<Vertex_handle> L; 382 search(v,C,L); 383 384 if (new_v) 385 { L.pop_front(); //first one was inserted in range_search ... 386 remove(v); 387 } 388 389 typename std::list<Vertex_handle>::const_iterator iter = L.begin(); 390 for(;iter != L.end() ;iter++){ *res= *iter; res++; } 391 return res; 392 } 393 394 395 template<class OutputIterator> range_search(const Point & a,const Point & b,const Point & c,OutputIterator res)396 OutputIterator range_search(const Point& a, const Point& b, const Point& c,OutputIterator res) 397 { int orient = (int)(tr_orientation(a,b,c)); 398 Circle C = tr_createcircle_3p(a,b,c); 399 std::list<Vertex_handle> L; 400 range_search(C,std::back_inserter(L)); 401 402 typename std::list<Vertex_handle>::const_iterator it = L.begin(); 403 404 for(;it != L.end();it++) 405 { Point p = (*it)->point(); 406 if ( ((int)(tr_orientation(a,b,p))) == - orient || 407 ((int)(tr_orientation(b,c,p))) == - orient || 408 ((int)(tr_orientation(c,a,p))) == - orient ) { } 409 else { *res = *it; res++; } 410 } 411 return res; 412 } 413 414 415 template<class OutputIterator> range_search(const Point & a1,const Point & b1,const Point & c1,const Point & d1,OutputIterator res)416 OutputIterator range_search(const Point& a1, const Point& b1, const Point& c1,const Point& 417 d1,OutputIterator res) 418 // a1 upper left, b1 lower left , c1 lower right 419 { 420 //Point b(c.xcoord(),a.ycoord()); 421 //Point d(a.xcoord(),c.ycoord()); 422 Point a=a1,b=b1,c=c1,d=d1; 423 424 if (tr_orientation(a,b,c) == RIGHT_TURN) 425 { Point tmp = b; 426 b = d; 427 d = tmp; 428 } 429 430 Circle C = tr_createcircle_3p(a,b,c); 431 432 std::list<Vertex_handle> L; 433 range_search(C,std::back_inserter(L)); 434 typename std::list<Vertex_handle>::const_iterator it = L.begin(); 435 436 for(;it != L.end();it++) 437 { Point p = (*it)->point(); 438 if ( tr_orientation(a,b,p) == RIGHT_TURN || tr_orientation(b,c,p) == RIGHT_TURN || 439 tr_orientation(c,d,p) == RIGHT_TURN || tr_orientation(d,a,p) == RIGHT_TURN ) { } 440 else { *res = *it; res++; } 441 } 442 return res; 443 } 444 445 }; 446 447 448 } //namespace CGAL 449 450 451 #endif 452