1 // Copyright (c) 2003,2004,2006 INRIA Sophia-Antipolis (France). 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/Apollonius_graph_2/include/CGAL/Apollonius_graph_2/Finite_edge_test_C2.h $ 7 // $Id: Finite_edge_test_C2.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) : Menelaos Karavelas <mkaravel@iacm.forth.gr> 12 13 14 15 #ifndef CGAL_APOLLONIUS_GRAPH_2_FINITE_EDGE_TEST_C2_H 16 #define CGAL_APOLLONIUS_GRAPH_2_FINITE_EDGE_TEST_C2_H 17 18 #include <CGAL/license/Apollonius_graph_2.h> 19 20 21 #include <CGAL/Apollonius_graph_2/basic.h> 22 23 #include <CGAL/Apollonius_graph_2/Predicate_constructions_C2.h> 24 25 #include <CGAL/functions_on_signs.h> 26 #include <CGAL/Apollonius_graph_2/compare_quadratic.h> 27 28 namespace CGAL { 29 30 namespace ApolloniusGraph_2 { 31 32 //-------------------------------------------------------------------- 33 34 template< class K > 35 class Orientation_wrt_symmetry_axis_2 36 { 37 public: 38 typedef typename K::Point_2 Point_2; 39 typedef Voronoi_circle_2<K> Voronoi_circle; 40 typedef typename K::FT FT; 41 typedef typename K::Orientation Orientation; 42 43 public: 44 45 Orientation operator()46 operator()(const Voronoi_circle& vc, const Point_2& p1, 47 const Point_2& p2, const Field_with_sqrt_tag&) const 48 { 49 FT a = vc.a1() + vc.a2() * CGAL::sqrt(vc.delta()); 50 FT b = vc.b1() + vc.b2() * CGAL::sqrt(vc.delta()); 51 FT det = a * (p2.y() - p1.y()) - b * (p2.x() - p1.x()); 52 return CGAL::sign(det); 53 } 54 55 Orientation operator()56 operator()(const Voronoi_circle& vc, const Point_2& p1, 57 const Point_2& p2, const Integral_domain_without_division_tag&) const 58 { 59 FT dx = p2.x() - p1.x(); 60 FT dy = p2.y() - p1.y(); 61 FT A = vc.a1() * dy - vc.b1() * dx; 62 FT B = vc.a2() * dy - vc.b2() * dx; 63 return sign_a_plus_b_x_sqrt_c(A, B, vc.delta()); 64 } 65 }; 66 67 68 //-------------------------------------------------------------------- 69 70 template< class K > 71 class Compare_Voronoi_radii_2 72 { 73 public: 74 typedef Voronoi_circle_2<K> Voronoi_circle; 75 typedef typename K::FT FT; 76 typedef typename K::Sign Sign; 77 typedef typename K::Comparison_result Comparison_result; 78 79 private: 80 sign_of_P4(const FT & u,const FT & v,const FT & du,const FT & dv,const FT & dr,const FT & Du,const FT & Dv,const FT & Dr)81 Sign sign_of_P4(const FT& u, const FT& v, 82 const FT& du, const FT& dv, const FT& dr, 83 const FT& Du, const FT& Dv, const FT& Dr) const 84 { 85 std::pair<FT,FT> factors = factors_of_P4(u, v, du, dv, dr, 86 Du, Dv, Dr); 87 Sign s1 = CGAL::sign(factors.first); 88 Sign s2 = CGAL::sign(factors.second); 89 return s1 * s2; 90 } 91 92 std::pair<FT,FT> factors_of_P4(const FT & u,const FT & v,const FT & du,const FT & dv,const FT & dr,const FT & Du,const FT & Dv,const FT & Dr)93 factors_of_P4(const FT& u, const FT& v, 94 const FT& du, const FT& dv, const FT& dr, 95 const FT& Du, const FT& Dv, const FT& Dr) const 96 { 97 FT u2 = CGAL::square(u); 98 FT v2 = CGAL::square(v); 99 100 FT du2 = CGAL::square(du); 101 FT Du2 = CGAL::square(Du); 102 103 FT dv2 = CGAL::square(dv); 104 FT Dv2 = CGAL::square(Dv); 105 106 FT dr2 = CGAL::square(dr); 107 FT Dr2 = CGAL::square(Dr); 108 109 FT u2_P_v2 = u2 + v2; 110 FT u2_M_v2 = u2 - v2; 111 FT uv = FT(2) * u * v; 112 113 FT drDr = FT(2) * dr * Dr; 114 115 FT du2_P_dv2 = du2 + dv2; 116 FT Du2_P_Dv2 = Du2 + Dv2; 117 118 FT uU_P_vV = du * Du + dv * Dv; 119 FT uU_M_vV = du * Du - dv * Dv; 120 121 FT uV_P_Uv = du * Dv + Du * dv; 122 FT uV_M_Uv = du * Dv - Du * dv; 123 124 125 FT F1 = du2_P_dv2 * Dr2 + Du2_P_Dv2 * dr2 126 - uU_P_vV * drDr - CGAL::square(uV_M_Uv); 127 128 FT F2 = CGAL::square(u2_P_v2) * (du2_P_dv2 * Dr2 + Du2_P_Dv2 * dr2); 129 F2 -= u2_P_v2 * (u2_M_v2 * uU_M_vV + uv * uV_P_Uv) * drDr; 130 F2 -= CGAL::square(u2_M_v2 * uV_P_Uv - uv * uU_M_vV); 131 132 std::pair<FT, FT> factors(F1,F2); 133 return factors; 134 } 135 136 public: 137 Comparison_result operator()138 operator()(const Voronoi_circle& vc1, const Voronoi_circle& vc2, 139 const Field_with_sqrt_tag&) const 140 { 141 FT c1 = (vc1.c1() + vc1.c2() * CGAL::sqrt(vc1.delta())) / vc1.d(); 142 FT c2 = (vc2.c1() + vc2.c2() * CGAL::sqrt(vc2.delta())) / vc2.d(); 143 144 Comparison_result r = CGAL::compare(c2, c1); 145 return r; 146 } 147 148 // this is the naive way but without divisions and square roots; the 149 // degree becomes 36 in this case. 150 /* 151 Comparison_result 152 operator()(const Voronoi_circle& vc1, const Voronoi_circle& vc2, 153 Integral_domain_without_division_tag) 154 { 155 FT A = vc1.c1() * vc2.d() - vc2.c1() * vc1.d(); 156 FT B = vc1.c2() * vc2.d(); 157 FT C = -vc2.c2() * vc1.d(); 158 FT E = vc1.delta(); 159 FT F = vc2.delta(); 160 161 Sign s = sign_a_plus_b_x_sqrt_e_plus_c_x_sqrt_f(A,B,C,E,F); 162 163 if ( s == ZERO ) { return EQUAL; } 164 return ( s == POSITIVE ) ? SMALLER : LARGER; 165 } 166 */ 167 168 Comparison_result operator()169 operator()(const Voronoi_circle& vc1, const Voronoi_circle& vc2, 170 const Integral_domain_without_division_tag&) const 171 { 172 bool is_first_root1 = vc1.is_first_root(); 173 bool is_first_root2 = vc2.is_first_root(); 174 175 CGAL_precondition( CGAL::is_positive(vc1.alpha()) ); 176 CGAL_precondition( CGAL::is_positive(vc2.alpha()) ); 177 178 Comparison_result r; 179 if ( is_first_root1 && is_first_root2 ) { 180 r = ke_compare_l1_l2(vc1.alpha(), vc1.beta(), vc1.gamma(), 181 vc2.alpha(), vc2.beta(), vc2.gamma()); 182 } else if ( is_first_root1 && !is_first_root2 ) { 183 r = ke_compare_l1_r2(vc1.alpha(), vc1.beta(), vc1.gamma(), 184 vc2.alpha(), vc2.beta(), vc2.gamma()); 185 } else if ( !is_first_root1 && is_first_root2 ) { 186 r = ke_compare_r1_l2(vc1.alpha(), vc1.beta(), vc1.gamma(), 187 vc2.alpha(), vc2.beta(), vc2.gamma()); 188 } else { 189 r = ke_compare_r1_r2(vc1.alpha(), vc1.beta(), vc1.gamma(), 190 vc2.alpha(), vc2.beta(), vc2.gamma()); 191 } 192 193 #ifdef COMPARATOR_PROFILER 194 if ( comparator_profiler::count_cases ) { 195 // count cases only for the tree r1-r2 196 if ( !is_first_root1 && !is_first_root2 ) { 197 comparator_profiler::count_case(vc1.alpha(), vc1.beta(), 198 vc1.gamma(), 199 vc2.alpha(), vc2.beta(), 200 vc2.gamma()); 201 } 202 } 203 #endif 204 205 if ( r == EQUAL ) { return EQUAL; } 206 return ( r == LARGER ) ? SMALLER : LARGER; 207 } 208 209 // this uses the DFMT trees; slightly slower but same degree (20). 210 /* 211 Comparison_result 212 operator()(const Voronoi_circle& vc1, const Voronoi_circle& vc2, 213 Integral_domain_without_division_tag) 214 { 215 bool is_first_root1 = vc1.is_first_root(); 216 bool is_first_root2 = vc2.is_first_root(); 217 218 CGAL_precondition( CGAL::is_positive(vc1.alpha()) ); 219 CGAL_precondition( CGAL::is_positive(vc2.alpha()) ); 220 221 Comparison_result r; 222 if ( is_first_root1 && is_first_root2 ) { 223 r = dfmt_compare_l1_l2(vc1.alpha(), vc1.beta(), vc1.gamma(), 224 vc2.alpha(), vc2.beta(), vc2.gamma()); 225 } else if ( is_first_root1 && !is_first_root2 ) { 226 r = dfmt_compare_l1_r2(vc1.alpha(), vc1.beta(), vc1.gamma(), 227 vc2.alpha(), vc2.beta(), vc2.gamma()); 228 } else if ( !is_first_root1 && is_first_root2 ) { 229 r = dfmt_compare_r1_l2(vc1.alpha(), vc1.beta(), vc1.gamma(), 230 vc2.alpha(), vc2.beta(), vc2.gamma()); 231 } else { 232 r = dfmt_compare_r1_r2(vc1.alpha(), vc1.beta(), vc1.gamma(), 233 vc2.alpha(), vc2.beta(), vc2.gamma()); 234 } 235 236 if ( r == EQUAL ) { return EQUAL; } 237 return ( r == LARGER ) ? SMALLER : LARGER; 238 } 239 */ 240 }; 241 242 243 //-------------------------------------------------------------------- 244 245 template< class K > 246 class Order_on_finite_bisector_2 247 { 248 public: 249 typedef Voronoi_circle_2<K> Voronoi_circle; 250 typedef typename K::Site_2 Site_2; 251 typedef typename K::FT FT; 252 typedef typename K::Comparison_result Comparison_result; 253 typedef typename K::Orientation Orientation; 254 255 typedef Compare_Voronoi_radii_2<K> Compare_Voronoi_radii; 256 257 typedef Orientation_wrt_symmetry_axis_2<K> 258 Orientation_wrt_symmetry_axis; 259 260 public: 261 template<class Method_tag> 262 Comparison_result operator()263 operator()(const Voronoi_circle& vc1, const Voronoi_circle& vc2, 264 const Site_2& p1, const Site_2& p2, 265 const Method_tag& tag) const 266 { 267 #ifdef AG2_PROFILE_PREDICATES 268 ag2_predicate_profiler::order_on_bisector_counter++; 269 #endif 270 271 Orientation o1 = 272 Orientation_wrt_symmetry_axis()(vc1, p1.point(), p2.point(), tag); 273 Orientation o2 = 274 Orientation_wrt_symmetry_axis()(vc2, p1.point(), p2.point(), tag); 275 276 Comparison_result cr; 277 if ( o1 == LEFT_TURN ) { 278 if ( o2 != LEFT_TURN ) { return SMALLER; } 279 Comparison_result r = Compare_Voronoi_radii()(vc1, vc2, tag); 280 281 if ( r == EQUAL ) { 282 cr = EQUAL; 283 } else { 284 cr = (r == LARGER ) ? SMALLER : LARGER; 285 } 286 } else if ( o1 == COLLINEAR ) { 287 if ( o2 == COLLINEAR ) { 288 cr = EQUAL; 289 } else { 290 cr = (o2 == LEFT_TURN) ? LARGER : SMALLER; 291 } 292 } else { 293 if ( o2 != RIGHT_TURN ) { 294 cr = LARGER; 295 } else { 296 Comparison_result r = 297 Compare_Voronoi_radii()(vc1, vc2, tag); 298 cr = r; 299 } 300 } 301 302 return cr; 303 } 304 }; 305 306 307 //-------------------------------------------------------------------- 308 309 template < class K > 310 class Finite_edge_interior_conflict 311 { 312 public: 313 typedef typename K::Site_2 Site_2; 314 typedef Weighted_point_inverter_2<K> Weighted_point_inverter; 315 typedef Inverted_weighted_point_2<K> Inverted_weighted_point; 316 typedef Voronoi_radius_2<K> Voronoi_radius; 317 typedef Voronoi_circle_2<K> Voronoi_circle; 318 typedef Bitangent_line_2<K> Bitangent_line; 319 typedef typename K::FT FT; 320 typedef typename K::Sign Sign; 321 typedef typename K::Bounded_side Bounded_side; 322 typedef typename K::Comparison_result Comparison_result; 323 324 typedef Bounded_side_of_CCW_circle_2<K> Bounded_side_of_CCW_circle; 325 typedef Order_on_finite_bisector_2<K> Order_on_finite_bisector; 326 327 typedef Sign_of_distance_from_bitangent_line_2<K> 328 Sign_of_distance_from_bitangent_line; 329 330 typedef Sign_of_distance_from_CCW_circle_2<K> 331 Sign_of_distance_from_CCW_circle; 332 333 public: 334 template<class Method_tag> 335 bool operator()336 operator()(const Site_2& p1, 337 const Site_2& p2, 338 const Site_2& p3, 339 const Site_2& p4, 340 const Site_2& q, bool b, const Method_tag& tag) const 341 { 342 #ifdef AG2_PROFILE_PREDICATES 343 ag2_predicate_profiler::shadow_region_type_counter++; 344 #endif 345 // 346 Weighted_point_inverter inverter(p1); 347 Inverted_weighted_point u2 = inverter(p2); 348 Inverted_weighted_point v = inverter(q); 349 // 350 Voronoi_radius vr_12q(u2, v); 351 Voronoi_radius vr_1q2 = vr_12q.get_symmetric(); 352 353 Bounded_side bs1 = Bounded_side_of_CCW_circle()(vr_12q, tag ); 354 Bounded_side bs2 = Bounded_side_of_CCW_circle()(vr_1q2, tag ); 355 356 bool is_bs1 = (bs1 == ON_UNBOUNDED_SIDE); 357 bool is_bs2 = (bs2 == ON_UNBOUNDED_SIDE); 358 359 // both the ccw and cw circles do not exist 360 if ( !is_bs1 && !is_bs2 ) { 361 return b; 362 } 363 364 // the ccw circle exists but not the cw 365 if ( is_bs1 && !is_bs2 ) { 366 return b; 367 } 368 369 // the cw circle exists but not the ccw 370 if ( !is_bs1 && is_bs2 ) { 371 return b; 372 } 373 374 375 // both circles exist 376 377 // check whether the shadow region is connected, i.e., whether it is 378 // of the form (a, b) or (-oo, a) U (b, +oo) 379 Bitangent_line bl_12(p1, p2); 380 381 Sign stc = 382 Sign_of_distance_from_bitangent_line()(bl_12, q, tag); 383 384 CGAL_assertion( stc != ZERO ); 385 bool is_shadow_region_connected = (stc == POSITIVE); 386 387 if ( is_shadow_region_connected ) { 388 if ( b ) { return true; } 389 390 Inverted_weighted_point u3 = inverter(p3); 391 Bitangent_line blinv_23(u2, u3); 392 393 Voronoi_circle vc_123(blinv_23); 394 Voronoi_circle vc_12q(vr_12q); 395 396 397 Comparison_result r = 398 Order_on_finite_bisector()(vc_123, vc_12q, p1, p2, tag); 399 400 if ( r != SMALLER ) { return false; } 401 402 Inverted_weighted_point u4 = inverter(p4); 403 Bitangent_line blinv_42(u4, u2); 404 405 Voronoi_circle vc_142(blinv_42); 406 Voronoi_circle vc_1q2(vr_1q2); 407 r = Order_on_finite_bisector()(vc_142, vc_1q2, p1, p2, tag); 408 409 return ( r == LARGER ); 410 } 411 412 // the shadow region is of the form (-oo, a) U (b, +oo) 413 if ( !b ) { return false; } 414 415 Inverted_weighted_point u3 = inverter(p3); 416 Bitangent_line blinv_23(u2, u3); 417 418 Voronoi_circle vc_123(blinv_23); 419 Voronoi_circle vc_1q2(vr_1q2); 420 Comparison_result r = 421 Order_on_finite_bisector()(vc_123, vc_1q2, p1, p2, tag); 422 423 if ( r != SMALLER ) { return true; } 424 425 Inverted_weighted_point u4 = inverter(p4); 426 Bitangent_line blinv_42(u4, u2); 427 428 Voronoi_circle vc_142(blinv_42); 429 Voronoi_circle vc_12q(vr_12q); 430 r = Order_on_finite_bisector()(vc_142, vc_12q, p1, p2, tag); 431 432 return ( r != LARGER ); 433 } 434 }; 435 436 //-------------------------------------------------------------------- 437 438 template < class K > 439 class Finite_edge_interior_conflict_degenerated 440 { 441 public: 442 typedef typename K::Site_2 Site_2; 443 typedef Weighted_point_inverter_2<K> Weighted_point_inverter; 444 typedef Inverted_weighted_point_2<K> Inverted_weighted_point; 445 typedef Voronoi_radius_2<K> Voronoi_radius; 446 typedef Voronoi_circle_2<K> Voronoi_circle; 447 typedef Bitangent_line_2<K> Bitangent_line; 448 typedef typename K::FT FT; 449 typedef typename K::Sign Sign; 450 typedef typename K::Comparison_result Comparison_result; 451 typedef typename K::Bounded_side Bounded_side; 452 453 typedef Bounded_side_of_CCW_circle_2<K> Bounded_side_of_CCW_circle; 454 typedef Order_on_finite_bisector_2<K> Order_on_finite_bisector; 455 456 typedef Sign_of_distance_from_bitangent_line_2<K> 457 Sign_of_distance_from_bitangent_line; 458 459 typedef Sign_of_distance_from_CCW_circle_2<K> 460 Sign_of_distance_from_CCW_circle; 461 public: 462 template<class Method_tag> 463 bool operator()464 operator()(const Site_2& p1, const Site_2& p2, const Site_2& p3, 465 const Site_2& q, bool b, const Method_tag& tag) const 466 { 467 #ifdef AG2_PROFILE_PREDICATES 468 ag2_predicate_profiler::shadow_region_type_counter++; 469 #endif 470 // 471 Weighted_point_inverter inverter(p1); 472 Inverted_weighted_point u2 = inverter(p2); 473 Inverted_weighted_point v = inverter(q); 474 475 Voronoi_radius vr_12q(u2, v); 476 Voronoi_radius vr_1q2 = vr_12q.get_symmetric(); 477 478 Bounded_side bs1 = Bounded_side_of_CCW_circle()(vr_12q, tag ); 479 Bounded_side bs2 = Bounded_side_of_CCW_circle()(vr_1q2, tag ); 480 481 bool is_bs1 = (bs1 == ON_UNBOUNDED_SIDE); 482 bool is_bs2 = (bs2 == ON_UNBOUNDED_SIDE); 483 484 // both the ccw and cw circles do not exist 485 if ( !is_bs1 && !is_bs2 ) { 486 return b; 487 } 488 489 // the ccw circle exists but not the cw 490 if ( is_bs1 && !is_bs2 ) { 491 return b; 492 } 493 494 // the cw circle exists but not the ccw 495 if ( !is_bs1 && is_bs2 ) { 496 return b; 497 } 498 499 // both circles exist 500 501 // check whether the shadow region is connected, i.e., whether it is 502 // of the form (a, b) or (-oo, a) U (b, +oo) 503 Bitangent_line bl_12(p1, p2); 504 505 Sign stc = Sign_of_distance_from_bitangent_line()(bl_12, q, tag); 506 507 Inverted_weighted_point u3 = inverter(p3); 508 Bitangent_line blinv_23(u2, u3); 509 510 CGAL_assertion( stc != ZERO ); 511 bool is_shadow_region_connected = (stc == POSITIVE); 512 513 if ( is_shadow_region_connected ) { 514 // the shadow region is of the form (a, b) 515 if ( b ) { return false; } 516 517 Voronoi_circle vc_123(blinv_23); 518 Voronoi_circle vc_12q(vr_12q); 519 520 Comparison_result r = 521 Order_on_finite_bisector()(vc_123, vc_12q, p1, p2, tag); 522 523 return ( r == SMALLER ); 524 } 525 526 // the shadow region is of the form (-oo, a) U (b, +oo) 527 if ( !b ) { return false; } 528 529 Voronoi_circle vc_123(blinv_23); 530 Voronoi_circle vc_1q2(vr_1q2); 531 Comparison_result r = 532 Order_on_finite_bisector()(vc_123, vc_1q2, p1, p2, tag); 533 534 return ( r != SMALLER ); 535 } 536 537 538 539 540 template<class Method_tag> 541 bool operator()542 operator()(const Site_2& p1, const Site_2& p2, 543 const Site_2& q, bool b, const Method_tag& tag) const 544 { 545 #ifdef AG2_PROFILE_PREDICATES 546 ag2_predicate_profiler::shadow_region_type_counter++; 547 #endif 548 // 549 Weighted_point_inverter inverter(p1); 550 Inverted_weighted_point u2 = inverter(p2); 551 Inverted_weighted_point v = inverter(q); 552 553 Voronoi_radius vr_12q(u2, v); 554 Voronoi_radius vr_1q2 = vr_12q.get_symmetric(); 555 556 Bounded_side bs1 = Bounded_side_of_CCW_circle()(vr_12q, tag ); 557 Bounded_side bs2 = Bounded_side_of_CCW_circle()(vr_1q2, tag ); 558 559 bool is_bs1 = (bs1 == ON_UNBOUNDED_SIDE); 560 bool is_bs2 = (bs2 == ON_UNBOUNDED_SIDE); 561 562 // both the ccw and cw circles do not exist 563 if ( !is_bs1 && !is_bs2 ) { 564 return b; 565 } 566 567 // only the ccw circle exists 568 if ( is_bs1 && !is_bs2 ) { return false; } 569 570 // only the cw circle exists 571 if ( !is_bs1 && is_bs2 ) { return false; } 572 573 // both circles exist 574 575 // check whether the shadow region is connected, i.e., whether it is 576 // of the form (a, b) or (-oo, a) U (b, +oo) 577 578 return !b; 579 } 580 }; 581 582 //-------------------------------------------------------------------- 583 584 585 template<class K, class MTag> 586 class Finite_edge_interior_conflict_2 587 { 588 public: 589 typedef K Kernel; 590 typedef MTag Method_tag; 591 592 typedef typename K::Site_2 Site_2; 593 594 private: 595 typedef Finite_edge_interior_conflict_degenerated<K> Test_degenerated; 596 typedef Finite_edge_interior_conflict<K> Test; 597 598 public: 599 typedef bool result_type; 600 struct argument_type {}; 601 602 inline operator()603 bool operator()(const Site_2& p1, const Site_2& p2, 604 const Site_2& q, bool b) const 605 { 606 return Test_degenerated()(p1, p2, q, b, Method_tag()); 607 } 608 609 inline operator()610 bool operator()(const Site_2& p1, const Site_2& p2, 611 const Site_2& p3, const Site_2& q, bool b) const 612 { 613 return Test_degenerated()(p1, p2, p3, q, b, Method_tag()); 614 } 615 616 inline operator()617 bool operator()(const Site_2& p1, const Site_2& p2, 618 const Site_2& p3, const Site_2& p4, 619 const Site_2& q, bool b) const 620 { 621 return Test()(p1, p2, p3, p4, q, b, Method_tag()); 622 } 623 }; 624 625 626 //-------------------------------------------------------------------- 627 628 } //namespace ApolloniusGraph_2 629 630 } //namespace CGAL 631 632 #endif // CGAL_APOLLONIUS_GRAPH_2_FINITE_EDGE_TEST_C2_H 633