1 // Copyright (c) 2006,2007,2008,2009,2010,2011 Tel-Aviv University (Israel). 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/Arrangement_on_surface_2/include/CGAL/Arr_geometry_traits/Conic_x_monotone_arc_2.h $ 7 // $Id: Conic_x_monotone_arc_2.h 6fcbee1 2020-04-21T17:12:21+03:00 Efi Fogel 8 // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial 9 // 10 // Author(s): Ron Wein <wein@post.tau.ac.il> 11 12 #ifndef CGAL_CONIC_X_MONOTONE_ARC_2_H 13 #define CGAL_CONIC_X_MONOTONE_ARC_2_H 14 15 #include <CGAL/license/Arrangement_on_surface_2.h> 16 17 /*! \file 18 * Header file for the _Conic_x_monotone_arc_2<Conic_arc_2> class. 19 */ 20 21 #include <map> 22 #include <ostream> 23 24 #include <boost/variant.hpp> 25 26 #include <CGAL/Arr_geometry_traits/Conic_intersections_2.h> 27 28 namespace CGAL { 29 30 /*! Representation of an x-monotone conic arc. 31 * The class is templated by a representation of a general bounded conic arc. 32 */ 33 34 template <typename Conic_arc_> 35 class _Conic_x_monotone_arc_2 : private Conic_arc_ { 36 public: 37 38 typedef Conic_arc_ Conic_arc_2; 39 typedef _Conic_x_monotone_arc_2<Conic_arc_2> Self; 40 41 typedef typename Conic_arc_2::Alg_kernel Alg_kernel; 42 typedef typename Conic_arc_2::Algebraic Algebraic; 43 44 typedef typename Conic_arc_2::Point_2 Point_2; 45 typedef typename Conic_arc_2::Conic_point_2 Conic_point_2; 46 47 // Type definition for the intersection points mapping. 48 typedef typename Conic_point_2::Conic_id Conic_id; 49 typedef std::pair<Conic_id, Conic_id> Conic_pair; 50 typedef std::pair<Conic_point_2, unsigned int> Intersection_point; 51 typedef std::list<Intersection_point> Intersection_list; 52 53 using Conic_arc_2::_sign_of_extra_data; 54 using Conic_arc_2::_is_between_endpoints; 55 using Conic_arc_2::_is_strictly_between_endpoints; 56 using Conic_arc_2::_conic_get_y_coordinates; 57 /*! 58 * \struct Less functor for Conic_pair. 59 */ 60 struct Less_conic_pair { operatorLess_conic_pair61 bool operator()(const Conic_pair& cp1, const Conic_pair& cp2) const 62 { 63 // Compare the pairs of IDs lexicographically. 64 return ((cp1.first < cp2.first) || 65 ((cp1.first == cp2.first) && (cp1.second < cp2.second))); 66 } 67 }; 68 69 typedef std::map<Conic_pair, 70 Intersection_list, 71 Less_conic_pair> Intersection_map; 72 typedef typename Intersection_map::value_type Intersection_map_entry; 73 typedef typename Intersection_map::iterator Intersection_map_iterator; 74 75 protected: 76 77 typedef Conic_arc_2 Base; 78 79 typedef typename Conic_arc_2::Integer Integer; 80 typedef typename Conic_arc_2::Nt_traits Nt_traits; 81 typedef typename Conic_arc_2::Rat_kernel Rat_kernel; 82 83 // Bit masks for the _info field (the two least significant bits are already 84 // used by the base class). 85 enum { 86 IS_VERTICAL_SEGMENT = 4, 87 IS_DIRECTED_RIGHT = 8, 88 DEGREE_1 = 16, 89 DEGREE_2 = 32, 90 DEGREE_MASK = 16 + 32, 91 PLUS_SQRT_DISC_ROOT = 64, 92 FACING_UP = 128, 93 FACING_DOWN = 256, 94 FACING_MASK = 128 + 256, 95 IS_SPECIAL_SEGMENT = 512 96 }; 97 98 Algebraic alg_r; // The coefficients of the supporting conic curve: 99 Algebraic alg_s; // 100 Algebraic alg_t; // r*x^2 + s*y^2 + t*xy + u*x + v*y +w = 0 , 101 Algebraic alg_u; // 102 Algebraic alg_v; // converted to algebraic numbers. 103 Algebraic alg_w; // 104 105 Conic_id _id; // The ID number of the supporting conic curve. 106 107 public: 108 109 /// \name Constrcution methods. 110 //@{ 111 112 /*! 113 * Default constructor. 114 */ _Conic_x_monotone_arc_2()115 _Conic_x_monotone_arc_2 () : 116 Base(), 117 _id() 118 {} 119 120 /*! 121 * Copy constructor. 122 * \param arc The copied arc. 123 */ _Conic_x_monotone_arc_2(const Self & arc)124 _Conic_x_monotone_arc_2(const Self& arc) : 125 Base(arc), 126 alg_r(arc.alg_r), 127 alg_s(arc.alg_s), 128 alg_t(arc.alg_t), 129 alg_u(arc.alg_u), 130 alg_v(arc.alg_v), 131 alg_w(arc.alg_w), 132 _id(arc._id) 133 {} 134 135 /*! Construct an x-monotone arc from a conic arc. 136 * \param arc The given (base) arc. 137 * \pre The given arc is x-monotone. 138 */ _Conic_x_monotone_arc_2(const Base & arc)139 _Conic_x_monotone_arc_2(const Base& arc) : 140 Base(arc), 141 _id() 142 { 143 CGAL_precondition(arc.is_valid() && arc.is_x_monotone()); 144 _set (); 145 } 146 147 /*! Construct an x-monotone arc from a conic arc. 148 * \param arc The given (base) arc. 149 * \param id The ID of the base arc. 150 */ _Conic_x_monotone_arc_2(const Base & arc,const Conic_id & id)151 _Conic_x_monotone_arc_2(const Base& arc, const Conic_id& id) : 152 Base(arc), 153 _id(id) 154 { 155 CGAL_precondition(arc.is_valid() && id.is_valid()); 156 _set (); 157 } 158 159 /*! Construct an x-monotone sub-arc from a conic arc. 160 * \param arc The given (base) arc. 161 * \param source The source point. 162 * \param target The target point. 163 * \param id The ID of the base arc. 164 */ _Conic_x_monotone_arc_2(const Base & arc,const Point_2 & source,const Point_2 & target,const Conic_id & id)165 _Conic_x_monotone_arc_2(const Base& arc, 166 const Point_2& source, const Point_2& target, 167 const Conic_id& id) : 168 Base(arc), 169 _id(id) 170 { 171 CGAL_precondition(arc.is_valid() && id.is_valid()); 172 173 // Set the two endpoints. 174 this->_source = source; 175 this->_target = target; 176 177 _set(); 178 } 179 180 /*! 181 * Construct a special segment connecting to given endpoints (for the usage 182 * of the landmarks point-location strategy). 183 * \param source The source point. 184 * \param target The target point. 185 */ _Conic_x_monotone_arc_2(const Point_2 & source,const Point_2 & target)186 _Conic_x_monotone_arc_2(const Point_2& source, const Point_2& target) : 187 Base() 188 { 189 // Set the basic properties and clear the _info bits. 190 this->_source = source; 191 this->_target = target; 192 this->_orient = COLLINEAR; 193 this->_info = 0; 194 195 // Check if the arc is directed right (the target is lexicographically 196 // greater than the source point), or to the left. 197 Alg_kernel ker; 198 Comparison_result dir_res = 199 ker.compare_xy_2_object()(this->_source, this->_target); 200 201 CGAL_precondition (dir_res != EQUAL); 202 // Invalid arc: 203 if (dir_res == EQUAL) return; 204 205 this->_info = (Conic_arc_2::IS_VALID | DEGREE_1); 206 if (dir_res == SMALLER) 207 this->_info = (this->_info | IS_DIRECTED_RIGHT); 208 209 // Compose the equation of the underlying line. 210 const Algebraic x1 = source.x(), y1 = source.y(); 211 const Algebraic x2 = target.x(), y2 = target.y(); 212 213 // The supporting line is A*x + B*y + C = 0, where: 214 // 215 // A = y2 - y1, B = x1 - x2, C = x2*y1 - x1*y2 216 // 217 // We use the extra data field to store the equation of this line. 218 this->_extra_data_P = new typename Base::Extra_data; 219 this->_extra_data_P->a = y2 - y1; 220 this->_extra_data_P->b = x1 - x2; 221 this->_extra_data_P->c = x2*y1 - x1*y2; 222 this->_extra_data_P->side = ZERO; 223 224 // Check if the segment is vertical. 225 if (CGAL::sign (this->_extra_data_P->b) == ZERO) 226 this->_info = (this->_info | IS_VERTICAL_SEGMENT); 227 228 // Mark that this is a special segment. 229 this->_info = (this->_info | IS_SPECIAL_SEGMENT); 230 231 return; 232 } 233 234 /*! 235 * Construct a special segment of a given line connecting to given 236 * endpoints. 237 * \param a, b, c The coefficients of the supporting line (ax + by + c = 0). 238 * \param source The source point. 239 * \param target The target point. 240 */ _Conic_x_monotone_arc_2(const Algebraic & a,const Algebraic & b,const Algebraic & c,const Point_2 & source,const Point_2 & target)241 _Conic_x_monotone_arc_2(const Algebraic& a, 242 const Algebraic& b, 243 const Algebraic& c, 244 const Point_2& source, const Point_2& target) : 245 Base() 246 { 247 // Make sure the two endpoints lie on the supporting line. 248 CGAL_precondition(CGAL::sign(a * source.x() + 249 b * source.y() + c) == CGAL::ZERO); 250 251 CGAL_precondition(CGAL::sign(a * target.x() + 252 b * target.y() + c) == CGAL::ZERO); 253 254 // Set the basic properties and clear the _info bits. 255 this->_source = source; 256 this->_target = target; 257 this->_orient = COLLINEAR; 258 this->_info = 0; 259 260 // Check if the arc is directed right (the target is lexicographically 261 // greater than the source point), or to the left. 262 Alg_kernel ker; 263 Comparison_result res = 264 ker.compare_x_2_object()(this->_source, this->_target); 265 266 this->_info = (Conic_arc_2::IS_VALID | DEGREE_1); 267 if (res == EQUAL) { 268 // Mark that the segment is vertical. 269 this->_info = (this->_info | IS_VERTICAL_SEGMENT); 270 271 // Compare the endpoints lexicographically. 272 res = ker.compare_y_2_object()(this->_source, this->_target); 273 274 CGAL_precondition (res != EQUAL); 275 if (res == EQUAL) { 276 // Invalid arc: 277 this->_info = 0; 278 return; 279 } 280 } 281 282 if (res == SMALLER) 283 this->_info = (this->_info | IS_DIRECTED_RIGHT); 284 285 // Store the coefficients of the line. 286 this->_extra_data_P = new typename Base::Extra_data; 287 this->_extra_data_P->a = a; 288 this->_extra_data_P->b = b; 289 this->_extra_data_P->c = c; 290 this->_extra_data_P->side = ZERO; 291 292 // Mark that this is a special segment. 293 this->_info = (this->_info | IS_SPECIAL_SEGMENT); 294 295 return; 296 } 297 298 /*! 299 * Assignment operator. 300 * \param arc The copied arc. 301 */ 302 const Self& operator=(const Self& arc) 303 { 304 CGAL_precondition (arc.is_valid()); 305 306 if (this == &arc) return (*this); 307 308 // Copy the base arc. 309 Base::operator= (arc); 310 311 // Set the rest of the properties. 312 alg_r = arc.alg_r; 313 alg_s = arc.alg_s; 314 alg_t = arc.alg_t; 315 alg_u = arc.alg_u; 316 alg_v = arc.alg_v; 317 alg_w = arc.alg_w; 318 319 _id = arc._id; 320 321 return (*this); 322 } 323 //@} 324 325 /// \name Accessing the arc properties. 326 //@{ 327 328 /*! 329 * Get the coefficients of the underlying conic. 330 */ r()331 const Integer& r() const { return (this->_r); } s()332 const Integer& s() const { return (this->_s); } t()333 const Integer& t() const { return (this->_t); } u()334 const Integer& u() const { return (this->_u); } v()335 const Integer& v() const { return (this->_v); } w()336 const Integer& w() const { return (this->_w); } 337 338 /*! 339 * Get the arc's source. 340 * \return The source point. 341 */ source()342 const Conic_point_2& source() const { return (this->_source); } 343 344 /*! Get the arc's target. 345 * \return The target point. 346 */ target()347 const Conic_point_2& target() const { return (this->_target); } 348 349 /*! Get the orientation of the arc. 350 * \return The orientation. 351 */ orientation()352 Orientation orientation() const { return (this->_orient); } 353 354 /*! Get the left endpoint of the arc. 355 */ left()356 const Conic_point_2& left () const 357 { 358 if ((this->_info & IS_DIRECTED_RIGHT) != 0) return (this->_source); 359 else return (this->_target); 360 } 361 362 /*! 363 * Get the right endpoint of the arc. 364 */ right()365 const Conic_point_2& right () const 366 { 367 if ((this->_info & IS_DIRECTED_RIGHT) != 0) return (this->_target); 368 else return (this->_source); 369 } 370 371 /*! 372 * Return true iff the conic arc is directed right iexicographically. 373 */ is_directed_right()374 bool is_directed_right() const 375 { return ((this->_info & IS_DIRECTED_RIGHT) != 0); } 376 377 /*! 378 * Get a bounding box for the conic arc. 379 * \return The bounding box. 380 */ bbox()381 Bbox_2 bbox() const { return (Base::bbox()); } 382 //@} 383 384 /// \name Predicates. 385 //@{ 386 387 /*! 388 * Check if the conic arc is a vertical segment. 389 */ is_vertical()390 bool is_vertical () const 391 { 392 return ((this->_info & IS_VERTICAL_SEGMENT) != 0); 393 } 394 395 /*! 396 * Check whether the given point lies on the arc. 397 * \param p The qury point. 398 * \param (true) if p lies on the arc; (false) otherwise. 399 */ contains_point(const Conic_point_2 & p)400 bool contains_point(const Conic_point_2& p) const 401 { 402 // First check if p lies on the supporting conic. We first check whether 403 // it is one of p's generating conic curves. 404 bool p_on_conic = false; 405 406 if (p.is_generating_conic(_id)) { 407 p_on_conic = true; 408 } 409 else { 410 // Check whether p satisfies the supporting conic equation. 411 p_on_conic = _is_on_supporting_conic(p.x(), p.y()); 412 413 if (p_on_conic) { 414 // As p lies on the supporting conic of our arc, add its ID to 415 // the list of generating conics for p. 416 Conic_point_2& p_non_const = const_cast<Conic_point_2&> (p); 417 p_non_const.set_generating_conic (_id); 418 } 419 } 420 421 if (! p_on_conic) 422 return (false); 423 424 // Check if p is between the endpoints of the arc. 425 return (_is_between_endpoints (p)); 426 } 427 //@} 428 429 /// \name Constructing points on the arc. 430 //@{ 431 432 /*! 433 * Compute a point on the arc with the same x-coordiante as the given point. 434 * \param p The given point. 435 * \pre The arc is not vertical and p is in the x-range of the arc. 436 * \return A point on the arc with the same x-coordiante as p. 437 */ point_at_x(const Point_2 & p)438 Point_2 point_at_x (const Point_2& p) const 439 { 440 // Make sure that p is in the x-range of the arc. 441 CGAL_precondition ((this->_info & IS_VERTICAL_SEGMENT) == 0); 442 443 CGAL_precondition_code ( 444 Alg_kernel ker; 445 ); 446 447 CGAL_precondition(ker.compare_x_2_object() (p, left()) != SMALLER && 448 ker.compare_x_2_object() (p, right()) != LARGER); 449 450 if (_is_special_segment()) { 451 // In case of a special segment, the equation of the supported line 452 // (a*x + b*y + c) = 0 is stored with the extra data field, and we 453 // simply have: 454 Algebraic _y = -(this->_extra_data_P->a*p.x() + 455 this->_extra_data_P->c) / 456 this->_extra_data_P->b; 457 458 // Return the computed point. 459 return (Point_2 (p.x(), _y)); 460 } 461 462 // Compute the y-coordinate according to the degree of the supporting 463 // conic curve. 464 Nt_traits nt_traits; 465 Algebraic y; 466 467 if ((this->_info & DEGREE_MASK) == DEGREE_1) { 468 // In case of a linear curve, the y-coordinate is a simple linear 469 // expression of x(p) (note that v is not 0 as the arc is not vertical): 470 // y = -(u*x(p) + w) / v 471 y = -(alg_u*p.x() + alg_w) / alg_v; 472 } 473 else if (this->_orient == COLLINEAR) { 474 CGAL_assertion (this->_extra_data_P != nullptr); 475 476 // In this case the equation of the supporting line is given by the 477 // extra data structure. 478 y = -(this->_extra_data_P->a * p.x() + 479 this->_extra_data_P->c) / this->_extra_data_P->b; 480 } 481 else { 482 CGAL_assertion((this->_info & DEGREE_MASK) == DEGREE_2); 483 484 // In this case the y-coordinate is one of solutions to the quadratic 485 // equation: 486 // s*y^2 + (t*x(p) + v)*y + (r*x(p)^2 + u*x(p) + w) = 0 487 Algebraic A = alg_s; 488 Algebraic B = alg_t*p.x() + alg_v; 489 Algebraic C = (alg_r*p.x() + alg_u)*p.x() + alg_w; 490 491 if (CGAL::sign(this->_s) == ZERO) { 492 // In this case A is 0 and we have a linear equation. 493 CGAL_assertion (CGAL::sign (B) != ZERO); 494 495 y = -C / B; 496 } 497 else { 498 // Solve the quadratic equation. 499 Algebraic disc = B*B - 4*A*C; 500 501 CGAL_assertion (CGAL::sign (disc) != NEGATIVE); 502 503 // We take either the root involving -sqrt(disc) or +sqrt(disc) 504 // based on the information flags. 505 if ((this->_info & PLUS_SQRT_DISC_ROOT) != 0) { 506 y = (nt_traits.sqrt (disc) - B) / (2*A); 507 } 508 else { 509 y = -(B + nt_traits.sqrt (disc)) / (2*A); 510 } 511 } 512 } 513 514 // Return the computed point. 515 return (Point_2 (p.x(), y)); 516 } 517 518 /*! 519 * Get a polyline approximating the conic arc. 520 * \param n The maximal number of sample points. 521 * \param oi An output iterator, whose value-type is pair<double,double> 522 * (representing an approximated point). 523 * In case the arc is a line segment, there are 2 output points, 524 * otherwise the arc is approximated by the polyline defined by 525 * (p_0, p_1, ..., p_n), where p_0 and p_n are the left and right 526 * endpoints of the arc, respectively. 527 */ 528 template <class OutputIterator> polyline_approximation(size_t n,OutputIterator oi)529 OutputIterator polyline_approximation (size_t n, 530 OutputIterator oi) const 531 { 532 CGAL_precondition (n != 0); 533 534 const double x_left = CGAL::to_double (left().x()); 535 const double y_left = CGAL::to_double (left().y()); 536 const double x_right = CGAL::to_double (right().x()); 537 const double y_right = CGAL::to_double (right().y()); 538 539 if (this->_orient == COLLINEAR) { 540 // In case of a line segment, return the two endpoints. 541 *oi++ = std::pair<double, double> (x_left, y_left); 542 *oi++ = std::pair<double, double> (x_right, y_right); 543 return oi; 544 } 545 546 // Otherwise, sample (n - 1) equally-spaced points in between. 547 const double app_r = CGAL::to_double (this->_r); 548 const double app_s = CGAL::to_double (this->_s); 549 const double app_t = CGAL::to_double (this->_t); 550 const double app_u = CGAL::to_double (this->_u); 551 const double app_v = CGAL::to_double (this->_v); 552 const double app_w = CGAL::to_double (this->_w); 553 const double x_jump = (x_right - x_left) / n; 554 double x, y; 555 const bool A_is_zero = (CGAL::sign(this->_s) == ZERO); 556 double A = app_s, B, C; 557 double disc; 558 size_t i; 559 560 *oi = std::pair<double, double>(x_left, y_left); // The left point. 561 ++oi; 562 for (i = 1; i < n; i++) { 563 x = x_left + x_jump*i; 564 565 // Solve the quadratic equation: A*x^2 + B*x + C = 0: 566 B = app_t*x + app_v; 567 C = (app_r*x + app_u)*x + app_w; 568 569 if (A_is_zero) { 570 y = -C / B; 571 } 572 else { 573 disc = B*B - 4*A*C; 574 575 if (disc < 0) disc = 0; 576 577 // We take either the root involving -sqrt(disc) or +sqrt(disc) 578 // based on the information flags. 579 if ((this->_info & PLUS_SQRT_DISC_ROOT) != 0) { 580 y = (std::sqrt(disc) - B) / (2*A); 581 } 582 else { 583 y = -(B + std::sqrt (disc)) / (2*A); 584 } 585 } 586 587 *oi++ = std::pair<double, double> (x, y); 588 } 589 *oi++ = std::pair<double, double> (x_right, y_right); // The right point. 590 591 return oi; 592 } 593 594 /*! Compare to arcs immediately to the right of their intersection point. 595 * \param arc The compared arc. 596 * \param p The reference intersection point. 597 * \return The relative position of the arcs to the right of p. 598 * \pre Both arcs we compare are not vertical segments. 599 */ compare_to_right(const Self & arc,const Conic_point_2 & p)600 Comparison_result compare_to_right(const Self& arc, 601 const Conic_point_2& p) const 602 { 603 CGAL_precondition((this->_info & IS_VERTICAL_SEGMENT) == 0 && 604 (arc._info & IS_VERTICAL_SEGMENT) == 0); 605 606 // In case one arc is facing upwards and another facing downwards, it is 607 // clear that the one facing upward is above the one facing downwards. 608 if (_has_same_supporting_conic (arc)) { 609 if ((this->_info & FACING_UP) != 0 && (arc._info & FACING_DOWN) != 0) 610 return LARGER; 611 else if ((this->_info & FACING_DOWN)!= 0 && (arc._info & FACING_UP) != 0) 612 return SMALLER; 613 614 // In this case the two arcs overlap. 615 CGAL_assertion((this->_info & FACING_MASK) == (arc._info & FACING_MASK)); 616 617 return EQUAL; 618 } 619 620 // Compare the slopes of the two arcs at p, using their first-order 621 // partial derivatives. 622 Algebraic slope1_numer, slope1_denom; 623 Algebraic slope2_numer, slope2_denom; 624 625 _derive_by_x_at (p, 1, slope1_numer, slope1_denom); 626 arc._derive_by_x_at (p, 1, slope2_numer, slope2_denom); 627 628 // Check if any of the slopes is vertical. 629 const bool is_vertical_slope1 = (CGAL::sign(slope1_denom) == ZERO); 630 const bool is_vertical_slope2 = (CGAL::sign(slope2_denom) == ZERO); 631 632 if (!is_vertical_slope1 && !is_vertical_slope2) { 633 // The two derivatives at p are well-defined: use them to determine 634 // which arc is above the other (the one with a larger slope is below). 635 Comparison_result slope_res = 636 CGAL::compare(slope1_numer*slope2_denom, slope2_numer*slope1_denom); 637 638 if (slope_res != EQUAL) return (slope_res); 639 640 // Use the second-order derivative. 641 _derive_by_x_at(p, 2, slope1_numer, slope1_denom); 642 arc._derive_by_x_at(p, 2, slope2_numer, slope2_denom); 643 644 slope_res = 645 CGAL::compare(slope1_numer*slope2_denom, slope2_numer*slope1_denom); 646 647 if (slope_res != EQUAL) return (slope_res); 648 649 // Use the third-order derivative. 650 _derive_by_x_at(p, 3, slope1_numer, slope1_denom); 651 arc._derive_by_x_at(p, 3, slope2_numer, slope2_denom); 652 653 slope_res = 654 CGAL::compare(slope1_numer*slope2_denom, slope2_numer*slope1_denom); 655 656 // \todo Handle higher-order derivatives: 657 CGAL_assertion (slope_res != EQUAL); 658 659 return (slope_res); 660 } 661 else if (!is_vertical_slope2) { 662 // The first arc has a vertical slope at p: check whether it is 663 // facing upwards or downwards and decide accordingly. 664 CGAL_assertion ((this->_info & FACING_MASK) != 0); 665 666 if ((this->_info & FACING_UP) != 0) return (LARGER); 667 return SMALLER; 668 } 669 else if (!is_vertical_slope1) { 670 // The second arc has a vertical slope at p_int: check whether it is 671 // facing upwards or downwards and decide accordingly. 672 CGAL_assertion ((arc._info & FACING_MASK) != 0); 673 674 if ((arc._info & FACING_UP) != 0) return (SMALLER); 675 return LARGER; 676 } 677 678 // The two arcs have vertical slopes at p_int: 679 // First check whether one is facing up and one down. In this case the 680 // comparison result is trivial. 681 if ((this->_info & FACING_UP) != 0 && (arc._info & FACING_DOWN) != 0) 682 return (LARGER); 683 else if ((this->_info & FACING_DOWN)!= 0 && (arc._info & FACING_UP)!= 0) 684 return SMALLER; 685 686 // Compute the second-order derivative by y and act according to it. 687 _derive_by_y_at (p, 2, slope1_numer, slope1_denom); 688 arc._derive_by_y_at (p, 2, slope2_numer, slope2_denom); 689 690 Comparison_result slope_res = 691 CGAL::compare(slope1_numer*slope2_denom, slope2_numer*slope1_denom); 692 693 // If necessary, use the third-order derivative by y. 694 if (slope_res == EQUAL) { 695 // \todo Check this! 696 _derive_by_y_at(p, 3, slope1_numer, slope1_denom); 697 arc._derive_by_y_at(p, 3, slope2_numer, slope2_denom); 698 699 slope_res = 700 CGAL::compare(slope2_numer*slope1_denom, slope1_numer*slope2_denom); 701 } 702 703 // \todo Handle higher-order derivatives: 704 CGAL_assertion(slope_res != EQUAL); 705 706 if ((this->_info & FACING_UP) != 0 && (arc._info & FACING_UP) != 0) { 707 // Both are facing up. 708 return ((slope_res == LARGER) ? SMALLER : LARGER); 709 } 710 // Both are facing down. 711 return (slope_res); 712 } 713 714 /*! 715 * Compare to arcs immediately to the leftt of their intersection point. 716 * \param arc The compared arc. 717 * \param p The reference intersection point. 718 * \return The relative position of the arcs to the left of p. 719 * \pre Both arcs we compare are not vertical segments. 720 */ compare_to_left(const Self & arc,const Conic_point_2 & p)721 Comparison_result compare_to_left(const Self& arc, 722 const Conic_point_2& p) const 723 { 724 CGAL_precondition((this->_info & IS_VERTICAL_SEGMENT) == 0 && 725 (arc._info & IS_VERTICAL_SEGMENT) == 0); 726 727 // In case one arc is facing upwards and another facing downwards, it is 728 // clear that the one facing upward is above the one facing downwards. 729 if (_has_same_supporting_conic (arc)) { 730 if ((this->_info & FACING_UP) != 0 && (arc._info & FACING_DOWN) != 0) 731 return LARGER; 732 else if ((this->_info & FACING_DOWN)!= 0 && (arc._info & FACING_UP)!= 0) 733 return SMALLER; 734 735 // In this case the two arcs overlap. 736 CGAL_assertion((this->_info & FACING_MASK) == (arc._info & FACING_MASK)); 737 738 return EQUAL; 739 } 740 741 // Compare the slopes of the two arcs at p, using their first-order 742 // partial derivatives. 743 Algebraic slope1_numer, slope1_denom; 744 Algebraic slope2_numer, slope2_denom; 745 746 _derive_by_x_at(p, 1, slope1_numer, slope1_denom); 747 arc._derive_by_x_at(p, 1, slope2_numer, slope2_denom); 748 749 // Check if any of the slopes is vertical. 750 const bool is_vertical_slope1 = (CGAL::sign (slope1_denom) == ZERO); 751 const bool is_vertical_slope2 = (CGAL::sign (slope2_denom) == ZERO); 752 753 if (!is_vertical_slope1 && !is_vertical_slope2) { 754 // The two derivatives at p are well-defined: use them to determine 755 // which arc is above the other (the one with a larger slope is below). 756 Comparison_result slope_res = CGAL::compare(slope2_numer*slope1_denom, 757 slope1_numer*slope2_denom); 758 759 if (slope_res != EQUAL) return (slope_res); 760 761 // Use the second-order derivative. 762 _derive_by_x_at (p, 2, slope1_numer, slope1_denom); 763 arc._derive_by_x_at (p, 2, slope2_numer, slope2_denom); 764 765 slope_res = CGAL::compare (slope1_numer*slope2_denom, 766 slope2_numer*slope1_denom); 767 768 if (slope_res != EQUAL) return (slope_res); 769 770 // Use the third-order derivative. 771 _derive_by_x_at(p, 3, slope1_numer, slope1_denom); 772 arc._derive_by_x_at(p, 3, slope2_numer, slope2_denom); 773 774 slope_res = CGAL::compare(slope2_numer*slope1_denom, 775 slope1_numer*slope2_denom); 776 777 // \todo Handle higher-order derivatives: 778 CGAL_assertion (slope_res != EQUAL); 779 780 return (slope_res); 781 } 782 else if (!is_vertical_slope2) { 783 // The first arc has a vertical slope at p: check whether it is 784 // facing upwards or downwards and decide accordingly. 785 CGAL_assertion ((this->_info & FACING_MASK) != 0); 786 787 if ((this->_info & FACING_UP) != 0) return (LARGER); 788 return SMALLER; 789 } 790 else if (!is_vertical_slope1) { 791 // The second arc has a vertical slope at p_int: check whether it is 792 // facing upwards or downwards and decide accordingly. 793 CGAL_assertion ((arc._info & FACING_MASK) != 0); 794 795 if ((arc._info & FACING_UP) != 0) return (SMALLER); 796 return LARGER; 797 } 798 799 // The two arcs have vertical slopes at p_int: 800 // First check whether one is facing up and one down. In this case the 801 // comparison result is trivial. 802 if ((this->_info & FACING_UP) != 0 && (arc._info & FACING_DOWN) != 0) 803 return LARGER; 804 else if ((this->_info & FACING_DOWN)!= 0 && (arc._info & FACING_UP)!= 0) 805 return SMALLER; 806 807 // Compute the second-order derivative by y and act according to it. 808 _derive_by_y_at(p, 2, slope1_numer, slope1_denom); 809 arc._derive_by_y_at(p, 2, slope2_numer, slope2_denom); 810 811 Comparison_result slope_res = 812 CGAL::compare(slope2_numer*slope1_denom, slope1_numer*slope2_denom); 813 814 // If necessary, use the third-order derivative by y. 815 if (slope_res == EQUAL) { 816 // \todo Check this! 817 _derive_by_y_at(p, 3, slope1_numer, slope1_denom); 818 arc._derive_by_y_at(p, 3, slope2_numer, slope2_denom); 819 820 slope_res = 821 CGAL::compare(slope2_numer*slope1_denom, slope1_numer*slope2_denom); 822 } 823 824 // \todo Handle higher-order derivatives: 825 CGAL_assertion(slope_res != EQUAL); 826 827 if ((this->_info & FACING_UP) != 0 && (arc._info & FACING_UP) != 0) { 828 // Both are facing up. 829 return ((slope_res == LARGER) ? SMALLER : LARGER); 830 } 831 // Both are facing down. 832 return (slope_res); 833 } 834 835 /*! 836 * Compute the intersections with the given arc. 837 * \param arc The given intersecting arc. 838 * \param inter_map Maps conic pairs to lists of their intersection points. 839 * \param oi The output iterator. 840 * \return The past-the-end iterator. 841 */ 842 template <typename OutputIterator> intersect(const Self & arc,Intersection_map & inter_map,OutputIterator oi)843 OutputIterator intersect(const Self& arc, 844 Intersection_map& inter_map, 845 OutputIterator oi) const 846 { 847 typedef boost::variant<Intersection_point, Self> Intersection_result; 848 849 if (_has_same_supporting_conic(arc)) { 850 // Check for overlaps between the two arcs. 851 Self overlap; 852 853 if (_compute_overlap(arc, overlap)) { 854 // There can be just a single overlap between two x-monotone arcs: 855 *oi++ = Intersection_result(overlap); 856 return oi; 857 } 858 859 // In case there is not overlap and the supporting conics are the same, 860 // there cannot be any intersection points, unless the two arcs share 861 // an end point. 862 // Note that in this case we do not define the multiplicity of the 863 // intersection points we report. 864 Alg_kernel ker; 865 866 if (ker.equal_2_object()(left(), arc.left())) { 867 Intersection_point ip(left(), 0); 868 *oi++ = Intersection_result(ip); 869 } 870 871 if (ker.equal_2_object()(right(), arc.right())) { 872 Intersection_point ip(right(), 0); 873 *oi++ = Intersection_result(ip); 874 } 875 876 return oi; 877 } 878 879 // Search for the pair of supporting conics in the map (the first conic 880 // ID in the pair should be smaller than the second one, to guarantee 881 // uniqueness). 882 Conic_pair conic_pair; 883 Intersection_map_iterator map_iter; 884 Intersection_list inter_list; 885 bool invalid_ids = false; 886 887 if (_id.is_valid() && arc._id.is_valid()) { 888 if (_id < arc._id) conic_pair = Conic_pair (_id, arc._id); 889 else conic_pair = Conic_pair (arc._id, _id); 890 map_iter = inter_map.find (conic_pair); 891 } 892 else { 893 // In case one of the IDs is invalid, we do not look in the map neither 894 // we cache the results. 895 map_iter = inter_map.end(); 896 invalid_ids = true; 897 } 898 899 if (map_iter == inter_map.end()) { 900 // In case the intersection points between the supporting conics have 901 // not been computed before, compute them now and store them in the map. 902 _intersect_supporting_conics(arc, inter_list); 903 904 if (! invalid_ids) inter_map[conic_pair] = inter_list; 905 } 906 else { 907 // Obtain the precomputed intersection points from the map. 908 inter_list = (*map_iter).second; 909 } 910 911 // Go over the list of intersection points and report those that lie on 912 // both x-monotone arcs. 913 for (auto iter = inter_list.begin(); iter != inter_list.end(); ++iter) { 914 if (_is_between_endpoints((*iter).first) && 915 arc._is_between_endpoints((*iter).first)) 916 { 917 *oi++ = Intersection_result(*iter); 918 } 919 } 920 921 return oi; 922 } 923 //@} 924 925 /// \name Constructing x-monotone arcs. 926 //@{ 927 928 /*! 929 * Split the arc into two at a given split point. 930 * \param p The split point. 931 * \param c1 Output: The first resulting arc, lying to the left of p. 932 * \param c2 Output: The first resulting arc, lying to the right of p. 933 * \pre p lies in the interior of the arc (not one of its endpoints). 934 */ split(const Conic_point_2 & p,Self & c1,Self & c2)935 void split(const Conic_point_2& p, Self& c1, Self& c2) const 936 { 937 // Make sure that p lies on the interior of the arc. 938 CGAL_precondition_code(Alg_kernel ker); 939 940 CGAL_precondition (this->contains_point (p) && 941 ! ker.equal_2_object() (p, this->_source) && 942 ! ker.equal_2_object() (p, this->_target)); 943 944 // Make copies of the current arc. 945 c1 = *this; 946 c2 = *this; 947 948 // Assign the endpoints of the arc. 949 if ((this->_info & IS_DIRECTED_RIGHT) != 0) 950 { 951 // The arc is directed from left to right, so p becomes c1's target 952 // and c2's source. 953 c1._target = p; 954 c2._source = p; 955 956 if (! p.is_generating_conic (_id)) { 957 c1._target.set_generating_conic (_id); 958 c2._source.set_generating_conic (_id); 959 } 960 } 961 else 962 { 963 // The arc is directed from right to left, so p becomes c2's target 964 // and c1's source. 965 c1._source = p; 966 c2._target = p; 967 968 if (! p.is_generating_conic (_id)) { 969 c1._source.set_generating_conic (_id); 970 c2._target.set_generating_conic (_id); 971 } 972 } 973 974 return; 975 } 976 977 /*! 978 * Flip the arc. 979 * \return An arc with swapped source and target and a reverse orienation. 980 */ flip()981 Self flip() const 982 { 983 // Make a copy of the current arc. 984 Self arc = *this; 985 986 // Reverse the orientation. 987 if (this->_orient == CLOCKWISE) arc._orient = COUNTERCLOCKWISE; 988 else if (this->_orient == COUNTERCLOCKWISE) arc._orient = CLOCKWISE; 989 990 // Swap the source and the target. 991 arc._source = this->_target; 992 arc._target = this->_source; 993 994 // Change the direction bit among the information flags. 995 arc._info = (this->_info ^ IS_DIRECTED_RIGHT); 996 997 return arc; 998 } 999 1000 /*! 1001 * Trim the arc given its new endpoints. 1002 * \param ps The new source point. 1003 * \param pt The new target point. 1004 * \return The new trimmed arc. 1005 * \pre Both ps and pt lies on the arc and must conform with the current 1006 * direction of the arc. 1007 */ trim(const Conic_point_2 & ps,const Conic_point_2 & pt)1008 Self trim(const Conic_point_2& ps, const Conic_point_2& pt) const 1009 { 1010 // Make sure that both ps and pt lie on the arc. 1011 CGAL_precondition(this->contains_point (ps) && 1012 this->contains_point (pt)); 1013 1014 // Make sure that the endpoints conform with the direction of the arc. 1015 Self arc = *this; 1016 Alg_kernel ker; 1017 1018 if (! ((((this->_info & IS_DIRECTED_RIGHT) != 0) && 1019 ker.compare_xy_2_object() (ps, pt) == SMALLER) || 1020 (((this->_info & IS_DIRECTED_RIGHT) == 0) && 1021 ker.compare_xy_2_object() (ps, pt) == LARGER))) 1022 { 1023 // We are allowed to change the direction only in case of a segment. 1024 CGAL_assertion (this->_orient == COLLINEAR); 1025 arc._info = (this->_info ^ IS_DIRECTED_RIGHT); 1026 } 1027 1028 // Make a copy of the current arc and assign its endpoints. 1029 if (! ker.equal_2_object() (ps, this->_source)) { 1030 arc._source = ps; 1031 1032 if (! ps.is_generating_conic (_id)) 1033 arc._source.set_generating_conic (_id); 1034 } 1035 1036 if (! ker.equal_2_object() (pt, this->_target)) 1037 { 1038 arc._target = pt; 1039 1040 if (! pt.is_generating_conic (_id)) 1041 arc._target.set_generating_conic (_id); 1042 } 1043 1044 return (arc); 1045 } 1046 1047 /*! 1048 * Check whether the two arcs are equal (have the same graph). 1049 * \param arc The compared arc. 1050 * \return (true) if the two arcs have the same graph; (false) otherwise. 1051 */ equals(const Self & arc)1052 bool equals (const Self& arc) const 1053 { 1054 // The two arc must have the same supporting conic curves. 1055 if (! _has_same_supporting_conic (arc)) 1056 return false; 1057 1058 // Check that the arc endpoints are the same. 1059 Alg_kernel ker; 1060 1061 if (this->_orient == COLLINEAR) { 1062 CGAL_assertion(arc._orient == COLLINEAR); 1063 return((ker.equal_2_object()(this->_source, arc._source) && 1064 ker.equal_2_object()(this->_target, arc._target)) || 1065 (ker.equal_2_object()(this->_source, arc._target) && 1066 ker.equal_2_object()(this->_target, arc._source))); 1067 } 1068 1069 if (this->_orient == arc._orient) { 1070 // Same orientation - the source and target points must be the same. 1071 return (ker.equal_2_object()(this->_source, arc._source) && 1072 ker.equal_2_object()(this->_target, arc._target)); 1073 } 1074 else { 1075 // Reverse orientation - the source and target points must be swapped. 1076 return (ker.equal_2_object()(this->_source, arc._target) && 1077 ker.equal_2_object()(this->_target, arc._source)); 1078 } 1079 } 1080 1081 /*! Check whether it is possible to merge the arc with the given arc. 1082 * \param arc The query arc. 1083 * \return (true) if it is possible to merge the two arcs; 1084 * (false) otherwise. 1085 */ can_merge_with(const Self & arc)1086 bool can_merge_with(const Self& arc) const 1087 { 1088 // In order to merge the two arcs, they should have the same supporting 1089 // conic. 1090 if (! _has_same_supporting_conic(arc)) return false; 1091 1092 // Check if the left endpoint of one curve is the right endpoint of the 1093 // other. 1094 Alg_kernel ker; 1095 1096 return (ker.equal_2_object() (right(), arc.left()) || 1097 ker.equal_2_object() (left(), arc.right())); 1098 } 1099 1100 /*! Merge the current arc with the given arc. 1101 * \param arc The arc to merge with. 1102 * \pre The two arcs are mergeable. 1103 */ merge(const Self & arc)1104 void merge(const Self& arc) 1105 { 1106 CGAL_precondition (this->can_merge_with (arc)); 1107 1108 // Check if we should extend the arc to the left or to the right. 1109 Alg_kernel ker; 1110 1111 if (ker.equal_2_object() (right(), arc.left())) { 1112 // Extend the arc to the right. 1113 if ((this->_info & IS_DIRECTED_RIGHT) != 0) this->_target = arc.right(); 1114 else this->_source = arc.right(); 1115 } 1116 else { 1117 CGAL_precondition (ker.equal_2_object() (left(), arc.right())); 1118 1119 // Extend the arc to the left. 1120 if ((this->_info & IS_DIRECTED_RIGHT) != 0) 1121 this->_source = arc.left(); 1122 else 1123 this->_target = arc.left(); 1124 } 1125 1126 return; 1127 } 1128 is_upper()1129 bool is_upper() const 1130 { 1131 return ((this->_info & FACING_UP) != 0); 1132 } 1133 is_lower()1134 bool is_lower() const 1135 { 1136 return ((this->_info & FACING_DOWN) != 0); 1137 } 1138 //@} 1139 1140 private: 1141 1142 /// \name Auxiliary (private) functions. 1143 //@{ 1144 1145 /*! 1146 * Set the properties of the x-monotone conic arc (for the usage of the 1147 * constructors). 1148 */ _set()1149 void _set () 1150 { 1151 // Convert the coefficients of the supporting conic to algebraic numbers. 1152 Nt_traits nt_traits; 1153 1154 alg_r = nt_traits.convert (this->_r); 1155 alg_s = nt_traits.convert (this->_s); 1156 alg_t = nt_traits.convert (this->_t); 1157 alg_u = nt_traits.convert (this->_u); 1158 alg_v = nt_traits.convert (this->_v); 1159 alg_w = nt_traits.convert (this->_w); 1160 1161 // Set the generating conic ID for the source and target points. 1162 this->_source.set_generating_conic (_id); 1163 this->_target.set_generating_conic (_id); 1164 1165 // Clear the _info bits. 1166 this->_info = Conic_arc_2::IS_VALID; 1167 1168 // Check if the arc is directed right (the target is lexicographically 1169 // greater than the source point), or to the left. 1170 Alg_kernel ker; 1171 Comparison_result dir_res = ker.compare_xy_2_object() (this->_source, 1172 this->_target); 1173 1174 CGAL_assertion (dir_res != EQUAL); 1175 1176 if (dir_res == SMALLER) 1177 this->_info = (this->_info | IS_DIRECTED_RIGHT); 1178 1179 // Compute the degree of the underlying conic. 1180 if (CGAL::sign (this->_r) != ZERO || 1181 CGAL::sign (this->_s) != ZERO || 1182 CGAL::sign (this->_t) != ZERO) 1183 { 1184 this->_info = (this->_info | DEGREE_2); 1185 1186 if (this->_orient == COLLINEAR) 1187 { 1188 this->_info = (this->_info | IS_SPECIAL_SEGMENT); 1189 1190 if (ker.compare_x_2_object() (this->_source, this->_target) == EQUAL) 1191 { 1192 // The arc is a vertical segment: 1193 this->_info = (this->_info | IS_VERTICAL_SEGMENT); 1194 } 1195 1196 return; 1197 } 1198 } 1199 else 1200 { 1201 CGAL_assertion (CGAL::sign (this->_u) != ZERO || 1202 CGAL::sign (this->_v) != ZERO); 1203 1204 if (CGAL::sign (this->_v) == ZERO) 1205 { 1206 1207 // The supporting curve is of the form: _u*x + _w = 0 1208 this->_info = (this->_info | IS_VERTICAL_SEGMENT); 1209 } 1210 1211 this->_info = (this->_info | DEGREE_1); 1212 1213 return; 1214 } 1215 1216 if (this->_orient == COLLINEAR) 1217 return; 1218 1219 // Compute a midpoint between the source and the target and get the y-value 1220 // of the arc at its x-coordiante. 1221 Point_2 p_mid = ker.construct_midpoint_2_object() (this->_source, 1222 this->_target); 1223 Algebraic ys[2]; 1224 CGAL_assertion_code(int n_ys = ) 1225 _conic_get_y_coordinates (p_mid.x(), ys); 1226 1227 CGAL_assertion (n_ys != 0); 1228 1229 // Check which solution lies on the x-monotone arc. 1230 Point_2 p_arc_mid (p_mid.x(), ys[0]); 1231 1232 if (_is_strictly_between_endpoints (p_arc_mid)) 1233 { 1234 // Mark that we should use the -sqrt(disc) root for points on this 1235 // x-monotone arc. 1236 this->_info = (this->_info & ~PLUS_SQRT_DISC_ROOT); 1237 } 1238 else 1239 { 1240 CGAL_assertion (n_ys == 2); 1241 p_arc_mid = Point_2 (p_mid.x(), ys[1]); 1242 1243 CGAL_assertion (_is_strictly_between_endpoints (p_arc_mid)); 1244 1245 // Mark that we should use the +sqrt(disc) root for points on this 1246 // x-monotone arc. 1247 this->_info = (this->_info | PLUS_SQRT_DISC_ROOT); 1248 } 1249 1250 // Check whether the conic is facing up or facing down: 1251 // Check whether the arc (which is x-monotone of degree 2) lies above or 1252 // below the segement that contects its two end-points (x1,y1) and (x2,y2). 1253 // To do that, we find the y coordinate of a point on the arc whose x 1254 // coordinate is (x1+x2)/2 and compare it to (y1+y2)/2. 1255 Comparison_result res = ker.compare_y_2_object() (p_arc_mid, p_mid); 1256 1257 if (res == LARGER) 1258 { 1259 // The arc is above the connecting segment, so it is facing upwards. 1260 this->_info = (this->_info | FACING_UP); 1261 } 1262 else if (res == SMALLER) 1263 { 1264 // The arc is below the connecting segment, so it is facing downwards. 1265 this->_info = (this->_info | FACING_DOWN); 1266 } 1267 1268 return; 1269 } 1270 1271 /*! 1272 * Check if the arc is a special segment connecting two algebraic endpoints 1273 * (and has no undelying integer conic coefficients). 1274 */ _is_special_segment()1275 bool _is_special_segment () const 1276 { 1277 return ((this->_info & IS_SPECIAL_SEGMENT) != 0); 1278 } 1279 1280 /*! 1281 * Check whether the given point lies on the supporting conic of the arc. 1282 * \param px The x-coordinate of query point. 1283 * \param py The y-coordinate of query point. 1284 * \return (true) if p lies on the supporting conic; (false) otherwise. 1285 */ _is_on_supporting_conic(const Algebraic & px,const Algebraic & py)1286 bool _is_on_supporting_conic (const Algebraic& px, 1287 const Algebraic& py) const 1288 { 1289 CGAL::Sign _sign; 1290 1291 if (! _is_special_segment()) 1292 { 1293 // Check whether p satisfies the conic equation. 1294 // The point must satisfy: r*x^2 + s*y^2 + t*xy + u*x + v*y + w = 0. 1295 _sign = CGAL::sign ((alg_r*px + alg_t*py + alg_u) * px + 1296 (alg_s*py + alg_v) * py + 1297 alg_w); 1298 } 1299 else 1300 { 1301 // Check whether p satisfies the equation of the line stored with the 1302 // extra data. 1303 _sign = _sign_of_extra_data (px, py); 1304 } 1305 1306 return (_sign == ZERO); 1307 } 1308 1309 /*! 1310 * Check whether the two arcs have the same supporting conic. 1311 * \param arc The compared arc. 1312 * \return (true) if the two supporting conics are the same. 1313 */ _has_same_supporting_conic(const Self & arc)1314 bool _has_same_supporting_conic (const Self& arc) const 1315 { 1316 // Check if the two arcs originate from the same conic: 1317 if (_id == arc._id && _id.is_valid() && arc._id.is_valid()) 1318 return (true); 1319 1320 // In case both arcs are collinear, check if they have the same 1321 // supporting lines. 1322 if (this->_orient == COLLINEAR && arc._orient == COLLINEAR) 1323 { 1324 // Construct the two supporting lines and compare them. 1325 Alg_kernel ker; 1326 typename Alg_kernel::Construct_line_2 construct_line = 1327 ker.construct_line_2_object(); 1328 typename Alg_kernel::Line_2 l1 = construct_line (this->_source, 1329 this->_target); 1330 typename Alg_kernel::Line_2 l2 = construct_line (arc._source, 1331 arc._target); 1332 typename Alg_kernel::Equal_2 equal = ker.equal_2_object(); 1333 1334 if (equal (l1, l2)) 1335 return (true); 1336 1337 // Try to compare l1 with the opposite of l2. 1338 l2 = construct_line (arc._target, arc._source); 1339 1340 return (equal (l1, l2)); 1341 } 1342 else if (this->_orient == COLLINEAR || arc._orient == COLLINEAR) 1343 { 1344 // Only one arc is collinear, so the supporting curves cannot be the 1345 // same: 1346 return (false); 1347 } 1348 1349 // Check whether the coefficients of the two supporting conics are equal 1350 // up to a constant factor. 1351 Integer factor1 = 1; 1352 Integer factor2 = 1; 1353 1354 if (CGAL::sign (this->_r) != ZERO) 1355 factor1 = this->_r; 1356 else if (CGAL::sign (this->_s) != ZERO) 1357 factor1 = this->_s; 1358 else if (CGAL::sign (this->_t) != ZERO) 1359 factor1 = this->_t; 1360 else if (CGAL::sign (this->_u) != ZERO) 1361 factor1 = this->_u; 1362 else if (CGAL::sign (this->_v) != ZERO) 1363 factor1 = this->_v; 1364 else if (CGAL::sign (this->_w) != ZERO) 1365 factor1 = this->_w; 1366 1367 if (CGAL::sign (arc._r) != ZERO) 1368 factor2 = arc._r; 1369 else if (CGAL::sign (arc._s) != ZERO) 1370 factor2 = arc._s; 1371 else if (CGAL::sign (arc._t) != ZERO) 1372 factor2 = arc._t; 1373 else if (CGAL::sign (arc._u) != ZERO) 1374 1375 factor2 = arc._u; 1376 else if (CGAL::sign (arc._v) != ZERO) 1377 factor2 = arc._v; 1378 else if (CGAL::sign (arc._w) != ZERO) 1379 factor2 = arc._w; 1380 1381 return (CGAL::compare (this->_r * factor2, arc._r * factor1) == EQUAL && 1382 CGAL::compare (this->_s * factor2, arc._s * factor1) == EQUAL && 1383 CGAL::compare (this->_t * factor2, arc._t * factor1) == EQUAL && 1384 CGAL::compare (this->_u * factor2, arc._u * factor1) == EQUAL && 1385 CGAL::compare (this->_v * factor2, arc._v * factor1) == EQUAL && 1386 CGAL::compare (this->_w * factor2, arc._w * factor1) == EQUAL); 1387 } 1388 1389 /*! 1390 * Get the i'th order derivative by x of the conic at the point p=(x,y). 1391 * \param p The point where we derive. 1392 * \param i The order of the derivatives (either 1, 2 or 3). 1393 * \param slope_numer The numerator of the slope. 1394 * \param slope_denom The denominator of the slope. 1395 * \todo Allow higher order derivatives. 1396 */ _derive_by_x_at(const Point_2 & p,const unsigned int & i,Algebraic & slope_numer,Algebraic & slope_denom)1397 void _derive_by_x_at (const Point_2& p, const unsigned int& i, 1398 Algebraic& slope_numer, Algebraic& slope_denom) const 1399 { 1400 if (_is_special_segment()) 1401 { 1402 // Special treatment for special segments, given by (a*x + b*y + c = 0), 1403 // so their first-order derivative by x is simply -a/b. The higher-order 1404 // derivatives are all 0. 1405 if (i == 1) 1406 { 1407 if (CGAL::sign (this->_extra_data_P->b) != NEGATIVE) 1408 { 1409 slope_numer = - this->_extra_data_P->a; 1410 slope_denom = this->_extra_data_P->b; 1411 } 1412 else 1413 { 1414 slope_numer = this->_extra_data_P->a; 1415 slope_denom = - this->_extra_data_P->b; 1416 } 1417 } 1418 else 1419 { 1420 slope_numer = 0; 1421 slope_denom = 1; 1422 } 1423 1424 return; 1425 } 1426 1427 // The derivative by x of the conic 1428 // C: {r*x^2 + s*y^2 + t*xy + u*x + v*y + w = 0} 1429 // at the point p=(x,y) is given by: 1430 // 1431 // 2r*x + t*y + u alpha 1432 // y' = - ---------------- = - ------- 1433 // 2s*y + t*x + v beta 1434 // 1435 const Algebraic _two = 2; 1436 const Algebraic sl_numer = _two*alg_r*p.x() + alg_t*p.y() + alg_u; 1437 const Algebraic sl_denom = _two*alg_s*p.y() + alg_t*p.x() + alg_v; 1438 1439 if (i == 1) 1440 { 1441 // Make sure that the denominator is always positive. 1442 if (CGAL::sign (sl_denom) != NEGATIVE) 1443 { 1444 slope_numer = -sl_numer; 1445 slope_denom = sl_denom; 1446 } 1447 else 1448 { 1449 slope_numer = sl_numer; 1450 slope_denom = -sl_denom; 1451 } 1452 1453 return; 1454 } 1455 1456 // The second-order derivative is given by: 1457 // 1458 // s*alpha^2 - t*alpha*beta + r*beta^2 gamma 1459 // y'' = -2 ------------------------------------- = ------- 1460 // beta^3 delta 1461 // 1462 const Algebraic sl2_numer = alg_s * sl_numer*sl_numer - 1463 alg_t * sl_numer*sl_denom + 1464 alg_r * sl_denom*sl_denom; 1465 const Algebraic sl2_denom = sl_denom*sl_denom*sl_denom; 1466 1467 if (i == 2) 1468 { 1469 // Make sure that the denominator is always positive. 1470 if (CGAL::sign (sl_denom) != NEGATIVE) 1471 { 1472 slope_numer = -_two *sl2_numer; 1473 slope_denom = sl2_denom; 1474 } 1475 else 1476 { 1477 slope_numer = _two *sl2_numer; 1478 slope_denom = -sl2_denom; 1479 } 1480 1481 return; 1482 } 1483 1484 // The third-order derivative is given by: 1485 // 1486 // (2s*alpha - t*beta) * gamma 1487 // y''' = -6 ------------------------------ 1488 // beta^2 * delta 1489 // 1490 const Algebraic sl3_numer = (_two * alg_s * sl_numer - 1491 alg_t * sl_denom) * sl2_numer; 1492 const Algebraic sl3_denom = sl_denom*sl_denom * sl2_denom; 1493 1494 if (i == 3) 1495 { 1496 // Make sure that the denominator is always positive. 1497 if (CGAL::sign (sl_denom) != NEGATIVE) 1498 { 1499 slope_numer = -6 * sl3_numer; 1500 slope_denom = sl3_denom; 1501 } 1502 else 1503 { 1504 slope_numer = 6 * sl3_numer; 1505 slope_denom = -sl2_denom; 1506 } 1507 1508 return; 1509 } 1510 1511 // \todo Handle higher-order derivatives as well. 1512 CGAL_error(); 1513 return; 1514 } 1515 1516 /*! 1517 * Get the i'th order derivative by y of the conic at the point p=(x,y). 1518 * \param p The point where we derive. 1519 * \param i The order of the derivatives (either 1, 2 or 3). 1520 * \param slope_numer The numerator of the slope. 1521 * \param slope_denom The denominator of the slope. 1522 * \todo Allow higher order derivatives. 1523 */ _derive_by_y_at(const Point_2 & p,const int & i,Algebraic & slope_numer,Algebraic & slope_denom)1524 void _derive_by_y_at (const Point_2& p, const int& i, 1525 Algebraic& slope_numer, Algebraic& slope_denom) const 1526 { 1527 if (_is_special_segment()) 1528 { 1529 // Special treatment for special segments, given by (a*x + b*y + c = 0), 1530 // so their first-order derivative by x is simply -b/a. The higher-order 1531 // derivatives are all 0. 1532 if (i == 1) 1533 { 1534 if (CGAL::sign (this->_extra_data_P->a) != NEGATIVE) 1535 { 1536 slope_numer = - this->_extra_data_P->b; 1537 slope_denom = this->_extra_data_P->a; 1538 } 1539 else 1540 { 1541 slope_numer = this->_extra_data_P->b; 1542 slope_denom = - this->_extra_data_P->a; 1543 } 1544 } 1545 else 1546 { 1547 slope_numer = 0; 1548 slope_denom = 1; 1549 } 1550 1551 return; 1552 } 1553 1554 // The derivative by y of the conic 1555 // C: {r*x^2 + s*y^2 + t*xy + u*x + v*y + w = 0} 1556 // at the point p=(x,y) is given by: 1557 // 1558 // 2s*y + t*x + v alpha 1559 // x' = - ---------------- = ------- 1560 // 2r*x + t*y + u beta 1561 // 1562 const Algebraic _two = 2; 1563 const Algebraic sl_numer = _two*alg_s*p.y() + alg_t*p.x() + alg_v; 1564 const Algebraic sl_denom = _two*alg_r*p.x() + alg_t*p.y() + alg_u; 1565 1566 if (i == 1) 1567 { 1568 // Make sure that the denominator is always positive. 1569 if (CGAL::sign (sl_denom) != NEGATIVE) 1570 { 1571 slope_numer = -sl_numer; 1572 slope_denom = sl_denom; 1573 } 1574 else 1575 { 1576 slope_numer = sl_numer; 1577 slope_denom = -sl_denom; 1578 } 1579 1580 1581 return; 1582 } 1583 1584 // The second-order derivative is given by: 1585 // 1586 // r*alpha^2 - t*alpha*beta + s*beta^2 1587 // x'' = -2 ------------------------------------- 1588 // beta^3 1589 // 1590 const Algebraic sl2_numer = alg_r * sl_numer*sl_numer - 1591 alg_t * sl_numer*sl_denom + 1592 alg_s * sl_denom*sl_denom; 1593 const Algebraic sl2_denom = sl_denom*sl_denom*sl_denom; 1594 1595 if (i == 2) 1596 { 1597 1598 // Make sure that the denominator is always positive. 1599 if (CGAL::sign (sl_denom) != NEGATIVE) 1600 { 1601 slope_numer = -_two *sl2_numer; 1602 slope_denom = sl2_denom; 1603 } 1604 else 1605 { 1606 slope_numer = _two *sl2_numer; 1607 slope_denom = -sl2_denom; 1608 } 1609 1610 return; 1611 } 1612 1613 // The third-order derivative is given by: 1614 // 1615 // (2t*alpha - t*beta) * gamma 1616 // y''' = -6 ------------------------------ 1617 // beta^2 * delta 1618 // 1619 const Algebraic sl3_numer = (_two * alg_r * sl_numer - 1620 alg_t * sl_denom) * sl2_numer; 1621 const Algebraic sl3_denom = sl_denom*sl_denom * sl2_denom; 1622 1623 if (i == 3) 1624 { 1625 // Make sure that the denominator is always positive. 1626 if (CGAL::sign (sl_denom) != NEGATIVE) 1627 { 1628 slope_numer = -6 * sl3_numer; 1629 slope_denom = sl3_denom; 1630 } 1631 else 1632 { 1633 slope_numer = 6 * sl3_numer; 1634 slope_denom = -sl2_denom; 1635 } 1636 1637 return; 1638 } 1639 1640 // \todo Handle higher-order derivatives as well. 1641 CGAL_error(); 1642 return; 1643 } 1644 1645 /*! 1646 * Compute the overlap with a given arc, which is supposed to have the same 1647 * supporting conic curve as this arc. 1648 * \param arc The given arc. 1649 * \param overlap Output: The overlapping arc (if any). 1650 * \return Whether we found an overlap. 1651 */ _compute_overlap(const Self & arc,Self & overlap)1652 bool _compute_overlap (const Self& arc, Self& overlap) const 1653 { 1654 // Check if the two arcs are identical. 1655 if (equals (arc)) 1656 { 1657 overlap = arc; 1658 return (true); 1659 } 1660 1661 if (_is_strictly_between_endpoints (arc.left())) 1662 { 1663 if (_is_strictly_between_endpoints(arc.right())) 1664 { 1665 // Case 1 - *this: +-----------> 1666 // arc: +=====> 1667 overlap = arc; 1668 return (true); 1669 } 1670 else 1671 { 1672 // Case 2 - *this: +-----------> 1673 // arc: +=====> 1674 overlap = *this; 1675 1676 if ((overlap._info & IS_DIRECTED_RIGHT) != 0) 1677 overlap._source = arc.left(); 1678 else 1679 overlap._target = arc.left(); 1680 1681 return (true); 1682 } 1683 } 1684 else if (_is_strictly_between_endpoints (arc.right())) 1685 { 1686 // Case 3 - *this: +-----------> 1687 // arc: +=====> 1688 overlap = *this; 1689 1690 if ((overlap._info & IS_DIRECTED_RIGHT) != 0) 1691 overlap._target = arc.right(); 1692 else 1693 overlap._source = arc.right(); 1694 1695 return (true); 1696 } 1697 else if (arc._is_between_endpoints (this->_source) && 1698 arc._is_between_endpoints (this->_target) && 1699 (arc._is_strictly_between_endpoints(this->_source) || 1700 arc._is_strictly_between_endpoints(this->_target))) 1701 { 1702 // Case 4 - *this: +-----------> 1703 // arc: +================> 1704 overlap = *this; 1705 return (true); 1706 } 1707 1708 // If we reached here, there are no overlaps: 1709 return (false); 1710 } 1711 1712 /*! 1713 * Intersect the supporing conic curves of this arc and the given arc. 1714 * \param arc The arc to intersect with. 1715 * \param inter_list The list of intersection points. 1716 */ _intersect_supporting_conics(const Self & arc,Intersection_list & inter_list)1717 void _intersect_supporting_conics(const Self& arc, 1718 Intersection_list& inter_list) const 1719 { 1720 if (_is_special_segment() && ! arc._is_special_segment()) { 1721 // If one of the arcs is a special segment, make sure it is (arc). 1722 arc._intersect_supporting_conics(*this, inter_list); 1723 return; 1724 } 1725 1726 const int deg1 = ((this->_info & DEGREE_MASK) == DEGREE_1) ? 1 : 2; 1727 const int deg2 = ((arc._info & DEGREE_MASK) == DEGREE_1) ? 1 : 2; 1728 Nt_traits nt_traits; 1729 Algebraic xs[4]; 1730 int n_xs = 0; 1731 Algebraic ys[4]; 1732 int n_ys = 0; 1733 1734 if (arc._is_special_segment()) { 1735 // The second arc is a special segment (a*x + b*y + c = 0). 1736 if (_is_special_segment()) { 1737 // Both arc are sepcial segment, so they have at most one intersection 1738 // point. 1739 Algebraic denom = this->_extra_data_P->a * arc._extra_data_P->b - 1740 this->_extra_data_P->b * arc._extra_data_P->a; 1741 1742 if (CGAL::sign (denom) != CGAL::ZERO) { 1743 xs[0] = (this->_extra_data_P->b * arc._extra_data_P->c - 1744 this->_extra_data_P->c * arc._extra_data_P->b) / denom; 1745 n_xs = 1; 1746 1747 ys[0] = (this->_extra_data_P->c * arc._extra_data_P->a - 1748 this->_extra_data_P->a * arc._extra_data_P->c) / denom; 1749 n_ys = 1; 1750 } 1751 } 1752 else { 1753 // Compute the x-coordinates of the intersection points. 1754 n_xs = _compute_resultant_roots (nt_traits, 1755 alg_r, alg_s, alg_t, 1756 alg_u, alg_v, alg_w, 1757 deg1, 1758 arc._extra_data_P->a, 1759 arc._extra_data_P->b, 1760 arc._extra_data_P->c, 1761 xs); 1762 CGAL_assertion (n_xs <= 2); 1763 1764 // Compute the y-coordinates of the intersection points. 1765 n_ys = _compute_resultant_roots (nt_traits, 1766 alg_s, alg_r, alg_t, 1767 alg_v, alg_u, alg_w, 1768 deg1, 1769 arc._extra_data_P->b, 1770 arc._extra_data_P->a, 1771 arc._extra_data_P->c, 1772 ys); 1773 CGAL_assertion (n_ys <= 2); 1774 } 1775 } 1776 else { 1777 // Compute the x-coordinates of the intersection points. 1778 n_xs = _compute_resultant_roots (nt_traits, 1779 this->_r, this->_s, this->_t, 1780 this->_u, this->_v, this->_w, 1781 deg1, 1782 arc._r, arc._s, arc._t, 1783 arc._u, arc._v, arc._w, 1784 deg2, 1785 xs); 1786 CGAL_assertion (n_xs <= 4); 1787 1788 // Compute the y-coordinates of the intersection points. 1789 n_ys = _compute_resultant_roots (nt_traits, 1790 this->_s, this->_r, this->_t, 1791 this->_v, this->_u, this->_w, 1792 deg1, 1793 arc._s, arc._r, arc._t, 1794 arc._v, arc._u, arc._w, 1795 deg2, 1796 ys); 1797 CGAL_assertion(n_ys <= 4); 1798 } 1799 1800 // Pair the coordinates of the intersection points. As the vectors of 1801 // x and y-coordinates are sorted in ascending order, we output the 1802 // intersection points in lexicographically ascending order. 1803 unsigned int mult; 1804 int i, j; 1805 1806 if (arc._is_special_segment()) 1807 { 1808 if (n_xs == 0 || n_ys == 0) 1809 return; 1810 1811 if (n_xs == 1 && n_ys == 1) 1812 { 1813 // Single intersection. 1814 Conic_point_2 ip (xs[0], ys[0]); 1815 1816 ip.set_generating_conic (_id); 1817 ip.set_generating_conic (arc._id); 1818 1819 // In case the other curve is of degree 2, this is a tangency point. 1820 mult = (deg1 == 1 || _is_special_segment()) ? 1 : 2; 1821 inter_list.push_back(Intersection_point (ip, mult)); 1822 } 1823 else if (n_xs == 1 && n_ys == 2) 1824 { 1825 Conic_point_2 ip1 (xs[0], ys[0]); 1826 1827 ip1.set_generating_conic (_id); 1828 ip1.set_generating_conic (arc._id); 1829 1830 inter_list.push_back(Intersection_point (ip1, 1)); 1831 1832 Conic_point_2 ip2 (xs[0], ys[1]); 1833 1834 ip2.set_generating_conic (_id); 1835 ip2.set_generating_conic (arc._id); 1836 1837 inter_list.push_back(Intersection_point (ip2, 1)); 1838 } 1839 else if (n_xs == 2 && n_ys == 1) 1840 { 1841 Conic_point_2 ip1 (xs[0], ys[0]); 1842 1843 ip1.set_generating_conic (_id); 1844 ip1.set_generating_conic (arc._id); 1845 1846 inter_list.push_back(Intersection_point (ip1, 1)); 1847 1848 Conic_point_2 ip2 (xs[1], ys[0]); 1849 1850 ip2.set_generating_conic (_id); 1851 ip2.set_generating_conic (arc._id); 1852 1853 inter_list.push_back(Intersection_point (ip2, 1)); 1854 1855 } 1856 else { 1857 CGAL_assertion (n_xs == 2 && n_ys == 2); 1858 1859 // The x-coordinates and the y-coordinates are given in ascending 1860 // order. If the slope of the segment is positive, we pair the 1861 // coordinates as is - otherwise, we swap the pairs. 1862 int ind_first_y = 0, ind_second_y = 1; 1863 1864 if (CGAL::sign (arc._extra_data_P->b) == 1865 CGAL::sign(arc._extra_data_P->a)) 1866 { 1867 ind_first_y = 1; 1868 ind_second_y = 0; 1869 } 1870 1871 Conic_point_2 ip1(xs[0], ys[ind_first_y]); 1872 1873 ip1.set_generating_conic(_id); 1874 ip1.set_generating_conic(arc._id); 1875 1876 inter_list.push_back(Intersection_point (ip1, 1)); 1877 1878 Conic_point_2 ip2(xs[1], ys[ind_second_y]); 1879 1880 ip2.set_generating_conic (_id); 1881 ip2.set_generating_conic (arc._id); 1882 1883 inter_list.push_back(Intersection_point(ip2, 1)); 1884 } 1885 1886 return; 1887 } 1888 1889 for (i = 0; i < n_xs; i++) { 1890 for (j = 0; j < n_ys; j++) { 1891 if (_is_on_supporting_conic (xs[i], ys[j]) && 1892 arc._is_on_supporting_conic (xs[i], ys[j])) 1893 { 1894 // Create the intersection point and set its generating conics. 1895 Conic_point_2 ip(xs[i], ys[j]); 1896 1897 ip.set_generating_conic (_id); 1898 ip.set_generating_conic (arc._id); 1899 1900 // Compute the multiplicity of the intersection point. 1901 if (deg1 == 1 && deg2 == 1) mult = 1; 1902 else mult = _multiplicity_of_intersection_point(arc, ip); 1903 1904 // Insert the intersection point to the output list. 1905 inter_list.push_back(Intersection_point(ip, mult)); 1906 } 1907 } 1908 } 1909 1910 return; 1911 } 1912 1913 /*! 1914 * Compute the multiplicity of an intersection point. 1915 * \param arc The arc to intersect with. 1916 * \param p The intersection point. 1917 * \return The multiplicity of the intersection point. 1918 */ _multiplicity_of_intersection_point(const Self & arc,const Point_2 & p)1919 unsigned int _multiplicity_of_intersection_point (const Self& arc, 1920 const Point_2& p) const 1921 { 1922 CGAL_assertion (! _is_special_segment() || ! arc._is_special_segment()); 1923 1924 // Compare the slopes of the two arcs at p, using their first-order 1925 // partial derivatives. 1926 Algebraic slope1_numer, slope1_denom; 1927 Algebraic slope2_numer, slope2_denom; 1928 1929 _derive_by_x_at (p, 1, slope1_numer, slope1_denom); 1930 arc._derive_by_x_at (p, 1, slope2_numer, slope2_denom); 1931 1932 if (CGAL::compare (slope1_numer*slope2_denom, 1933 slope2_numer*slope1_denom) != EQUAL) 1934 { 1935 // Different slopes at p - the mutiplicity of p is 1: 1936 return (1); 1937 } 1938 1939 if (CGAL::sign (slope1_denom) != ZERO && 1940 CGAL::sign (slope2_denom) != ZERO) 1941 { 1942 // The curves do not have a vertical slope at p. 1943 // Compare their second-order derivative by x: 1944 _derive_by_x_at (p, 2, slope1_numer, slope1_denom); 1945 arc._derive_by_x_at (p, 2, slope2_numer, slope2_denom); 1946 } 1947 else 1948 { 1949 // Both curves have a vertical slope at p. 1950 // Compare their second-order derivative by y: 1951 _derive_by_y_at (p, 2, slope1_numer, slope1_denom); 1952 arc._derive_by_y_at (p, 2, slope2_numer, slope2_denom); 1953 } 1954 1955 if (CGAL::compare (slope1_numer*slope2_denom, 1956 slope2_numer*slope1_denom) != EQUAL) 1957 { 1958 // Different curvatures at p - the mutiplicity of p is 2: 1959 return (2); 1960 } 1961 1962 // If we reached here, the multiplicity of the intersection point is 3: 1963 return (3); 1964 } 1965 //@} 1966 1967 }; 1968 1969 /*! 1970 * Exporter for x-monotone conic arcs. 1971 */ 1972 template <class Conic_arc_2> 1973 std::ostream& operator<< (std::ostream& os, 1974 const _Conic_x_monotone_arc_2<Conic_arc_2>& arc) 1975 { 1976 // Output the supporting conic curve. 1977 os << "{" << CGAL::to_double(arc.r()) << "*x^2 + " 1978 << CGAL::to_double(arc.s()) << "*y^2 + " 1979 << CGAL::to_double(arc.t()) << "*xy + " 1980 << CGAL::to_double(arc.u()) << "*x + " 1981 << CGAL::to_double(arc.v()) << "*y + " 1982 << CGAL::to_double(arc.w()) << "}"; 1983 1984 // Output the endpoints. 1985 os << " : (" << CGAL::to_double(arc.source().x()) << "," 1986 << CGAL::to_double(arc.source().y()) << ") "; 1987 1988 if (arc.orientation() == CLOCKWISE) 1989 os << "--cw-->"; 1990 else if (arc.orientation() == COUNTERCLOCKWISE) 1991 os << "--ccw-->"; 1992 else 1993 os << "--l-->"; 1994 1995 os << " (" << CGAL::to_double(arc.target().x()) << "," 1996 << CGAL::to_double(arc.target().y()) << ")"; 1997 1998 return (os); 1999 } 2000 2001 } //namespace CGAL 2002 2003 #endif 2004