diff --git a/Packages/Arrangement/include/CGAL/Arr_conic_traits_2.h b/Packages/Arrangement/include/CGAL/Arr_conic_traits_2.h index 7769076cdef..48217199b8b 100644 --- a/Packages/Arrangement/include/CGAL/Arr_conic_traits_2.h +++ b/Packages/Arrangement/include/CGAL/Arr_conic_traits_2.h @@ -1,31 +1,10 @@ -// ====================================================================== -// -// Copyright (c) 2001 The CGAL Consortium -// -// This software and related documentation is part of an INTERNAL release -// of the Computational Geometry Algorithms Library (CGAL). It is not -// intended for general use. -// -// ---------------------------------------------------------------------- -// -// release : $CGAL_Revision: CGAL-2.4-I-5 $ -// release_date : $CGAL_Date: 2001/08/31 $ -// -// file : include/CGAL/Arr_conic_traits_2.h -// package : Arrangement (2.19) -// maintainer : Eyal Flato -// author(s) : Ron Wein -// coordinator : Tel-Aviv University (Dan Halperin ) -// -// ====================================================================== #ifndef CGAL_ARR_CONIC_TRAITS_H #define CGAL_ARR_CONIC_TRAITS_H #include -#include - -#include #include +#include "CGAL/Arrangement_2/Conic_arc_2.h" +#include CGAL_BEGIN_NAMESPACE @@ -34,37 +13,39 @@ CGAL_BEGIN_NAMESPACE // template -class Arr_conic_traits_2 +class Arr_conic_traits { public: - typedef Lazy_intersection_tag Intersection_category; + typedef Lazy_intersection_tag Intersection_category; - typedef _NT NT; + typedef _NT NT; - // The difference between Curve and X_curve is semantical only, + // The difference between Curve_2 and X_curve_2 is semantical only, // NOT syntactical. - typedef Conic_arc_2 Curve_2; - typedef Curve_2 X_curve_2; + typedef Conic_arc_2 Curve_2; + typedef Curve_2 X_curve_2; // Using typename to please compiler (e.g., CC with IRIX64 on mips) - typedef typename Curve_2::R R; - typedef typename Curve_2::Point_2 Point_2; - typedef typename Curve_2::Conic_2 Conic_2; + typedef typename Curve_2::Kernel R; + typedef typename Curve_2::Point_2 Point_2; + typedef typename Curve_2::Circle_2 Circle_2; + typedef typename Curve_2::Segment_2 Segment_2; + + // For backward compatibility: + typedef Curve_2 Curve; + typedef X_curve_2 X_curve; + typedef Point_2 Point; + typedef Circle_2 Circle; + typedef Segment_2 Segment; enum Curve_point_status { - UNDER_CURVE = -1, + UNDER_CURVE = -1, + ABOVE_CURVE = 1, + ON_CURVE = 2, CURVE_NOT_IN_RANGE = 0, - ABOVE_CURVE = 1, - ON_CURVE = 2 }; - // Obsolete, for backward compatibility - typedef Point_2 Point; - typedef X_curve_2 X_curve; - typedef Curve_2 Curve; - typedef Conic_2 Conic; - #ifdef CGAL_CONIC_ARC_USE_CACHING private: @@ -75,7 +56,7 @@ class Arr_conic_traits_2 public: // Constructor. - Arr_conic_traits_2() + Arr_conic_traits() {} ////////// Planar Map methods: ////////// @@ -120,8 +101,8 @@ class Arr_conic_traits_2 // Decide wether curve1 is above, below or equal to curve2 at the // x co-ordinate of the given point. Comparison_result curve_compare_at_x (const X_curve_2& curve1, - const X_curve_2& curve2, - const Point_2& p) const + const X_curve_2& curve2, + const Point_2& p) const { CGAL_precondition(is_x_monotone(curve1)); CGAL_precondition(is_x_monotone(curve2)); @@ -135,24 +116,24 @@ class Arr_conic_traits_2 // curve1 is a vertical segment. if (compare_x (curve1.source(), p) != EQUAL) { - // The curve's x co-ordinate is different than p's. - return (EQUAL); + // The curve's x co-ordinate is different than p's. + return (EQUAL); } else { - // Both the source and target of the curve has the same x co-ordinate - // as p. Make sure that ps1[0] has a smaller y value than ps1[1]. - n1 = 2; - if (compare_y (curve1.source(), curve1.target()) == SMALLER) - { - ps1[0] = curve1.source(); - ps1[1] = curve1.target(); - } - else - { - ps1[0] = curve1.target(); - ps1[1] = curve1.source(); - } + // Both the source and target of the curve has the same x co-ordinate + // as p. Make sure that ps1[0] has a smaller y value than ps1[1]. + n1 = 2; + if (compare_y (curve1.source(), curve1.target()) == SMALLER) + { + ps1[0] = curve1.source(); + ps1[1] = curve1.target(); + } + else + { + ps1[0] = curve1.target(); + ps1[1] = curve1.source(); + } } } else @@ -161,7 +142,7 @@ class Arr_conic_traits_2 n1 = curve1.get_points_at_x (p, ps1); if (n1 == 0) // p is not in the x-range of curve1. - return (EQUAL); + return (EQUAL); CGAL_assertion(n1 == 1); } @@ -175,24 +156,24 @@ class Arr_conic_traits_2 // curve2 is a vertical segment. if (compare_x (curve2.source(), p) != EQUAL) { - // The curve's x co-ordinate is different than p's. - return (EQUAL); + // The curve's x co-ordinate is different than p's. + return (EQUAL); } else { - // Both the source and target of the curve has the same x co-ordinate - // as p. Make sure that ps2[0] has a smaller y value than ps2[1]. - n2 = 2; - if (compare_y (curve2.source(), curve2.target()) == SMALLER) - { - ps2[0] = curve2.source(); - ps2[1] = curve2.target(); - } - else - { - ps2[0] = curve2.target(); - ps2[1] = curve2.source(); - } + // Both the source and target of the curve has the same x co-ordinate + // as p. Make sure that ps2[0] has a smaller y value than ps2[1]. + n2 = 2; + if (compare_y (curve2.source(), curve2.target()) == SMALLER) + { + ps2[0] = curve2.source(); + ps2[1] = curve2.target(); + } + else + { + ps2[0] = curve2.target(); + ps2[1] = curve2.source(); + } } } else @@ -201,7 +182,7 @@ class Arr_conic_traits_2 n2 = curve2.get_points_at_x (p, ps2); if (n2 == 0) // p is not in the x-range of curve2. - return (EQUAL); + return (EQUAL); CGAL_assertion(n2 == 1); } @@ -211,27 +192,27 @@ class Arr_conic_traits_2 { // Check if the vertical segment curve1 contains ps2[0] or ps2[1]. if (compare_y (ps1[0], ps2[0]) != LARGER && - compare_y (ps1[1], ps2[0]) != SMALLER) + compare_y (ps1[1], ps2[0]) != SMALLER) { - return (EQUAL); + return (EQUAL); } if (n2 == 2) { - if (compare_y (ps1[0], ps2[1]) != LARGER && - compare_y (ps1[1], ps2[1]) != SMALLER) - { - return (EQUAL); - } + if (compare_y (ps1[0], ps2[1]) != LARGER && + compare_y (ps1[1], ps2[1]) != SMALLER) + { + return (EQUAL); + } } } else if (n2 == 2) { // Check if the vertical segment curve2 contains ps1[0]. if (compare_y (ps2[0], ps1[0]) != LARGER && - compare_y (ps2[1], ps1[0]) != SMALLER) + compare_y (ps2[1], ps1[0]) != SMALLER) { - return (EQUAL); + return (EQUAL); } } @@ -244,8 +225,8 @@ class Arr_conic_traits_2 // Decide wether curve1 is above, below or equal to curve2 immediately to // the left of the x co-ordinate of the given point. Comparison_result curve_compare_at_x_left (const X_curve_2& curve1, - const X_curve_2& curve2, - const Point_2& p) const + const X_curve_2& curve2, + const Point_2& p) const { CGAL_precondition(is_x_monotone(curve1)); CGAL_precondition(is_x_monotone(curve2)); @@ -315,119 +296,15 @@ class Arr_conic_traits_2 if (result != EQUAL) return (result); - // In case one arc is facing upwards and another facing downwards, it is - // very clear that the one facing upward is above the one facing downwards. - if (curve1.has_same_base_conic(curve2)) - { - Comparison_result fac1 = curve1.facing(); - Comparison_result fac2 = curve2.facing(); - - if (fac1 == LARGER && fac2 == SMALLER) - return (LARGER); - else if (fac1 == SMALLER && fac2 == LARGER) - return (SMALLER); - } - - // Otherwise, the two curves do intersect at p_int = ps1[0] = ps2[0]: - // make a decision based on their partial derivatives. - const Point_2& p_int = ps1[0]; - const NT _zero = 0; - - NT slope1_numer, slope1_denom; - NT slope2_numer, slope2_denom; - - curve1.derive_by_x_at (p_int, 1, slope1_numer, slope1_denom); - curve2.derive_by_x_at (p_int, 1, slope2_numer, slope2_denom); - - const bool is_vertical_slope1 = (slope1_denom == _zero); - const bool is_vertical_slope2 = (slope2_denom == _zero); - - if (!is_vertical_slope1 && !is_vertical_slope2) - { - // The two curves have derivatives at p_int: use it to determine which - // one is above the other (the one with a smaller slope in above). - Comparison_result slope_res = /*CGAL::compare*/ - eps_compare(TO_APNT(slope2_numer*slope1_denom), - TO_APNT(slope1_numer*slope2_denom)); - - if (slope_res != EQUAL) - return (slope_res); - - // Use the second order derivative. - curve1.derive_by_x_at (p_int, 2, slope1_numer, slope1_denom); - curve2.derive_by_x_at (p_int, 2, slope2_numer, slope2_denom); - - slope_res = /*CGAL::compare*/ - eps_compare (TO_APNT(slope2_numer*slope1_denom), - TO_APNT(slope1_numer*slope2_denom)); - - CGAL_assertion(slope_res != EQUAL); - - return ((slope_res == LARGER) ? SMALLER : LARGER); - } - else if (!is_vertical_slope2) - { - // The first curve has a vertical slope at p_int: check whether it is - // facing upwards or downwards and decide accordingly. - Comparison_result fac1 = curve1.facing(); - - CGAL_assertion(fac1 != EQUAL); - return (fac1); - } - else if (!is_vertical_slope1) - { - // The secong curve has a vertical slope at p_int: check whether it is - // facing upwards or downwards and decide accordingly. - Comparison_result fac2 = curve2.facing(); - - CGAL_assertion(fac2 != EQUAL); - return ((fac2 == LARGER) ? SMALLER : LARGER); - } - else - { - // The two curves have vertical slopes at p_int: - // First check whether one is facing up and one down. In this case the - // comparison result is trivial. - Comparison_result fac1 = curve1.facing(); - Comparison_result fac2 = curve2.facing(); - - if (fac1 == LARGER && fac2 == SMALLER) - return (LARGER); - else if (fac1 == SMALLER && fac2 == LARGER) - return (SMALLER); - - // Compute the second order derivative by y and act according to it. - curve1.derive_by_y_at (p_int, 2, slope1_numer, slope1_denom); - curve2.derive_by_y_at (p_int, 2, slope2_numer, slope2_denom); - - Comparison_result slope_res = /*CGAL::compare*/ - eps_compare (TO_APNT(slope2_numer*slope1_denom), - TO_APNT(slope1_numer*slope2_denom)); - - CGAL_assertion(slope_res != EQUAL); - - if (fac1 == LARGER && fac2 == LARGER) - { - // Both are facing up. - return ((slope_res == LARGER) ? SMALLER : LARGER); - } - else - { - // Both are facing down. - return (slope_res); - } - } - - // We should never reach here: - CGAL_assertion(false); - return (EQUAL); - } + // The two curves intersect at ps1[0] = ps2[0] - proceed from here. + return (_curve_compare_at_intersection_left (curve1, curve2, ps1[0])); + } // Decide wether curve1 is above, below or equal to curve2 immediately to // the right of the x co-ordinate of the given point. Comparison_result curve_compare_at_x_right (const X_curve_2& curve1, - const X_curve_2& curve2, - const Point_2& p) const + const X_curve_2& curve2, + const Point_2& p) const { CGAL_precondition(is_x_monotone(curve1)); CGAL_precondition(is_x_monotone(curve2)); @@ -497,117 +374,13 @@ class Arr_conic_traits_2 if (result != EQUAL) return (result); - // In case one arc is facing upwards and another facing downwards, it is - // very clear that the one facing upward is above the one facing downwards. - if (curve1.has_same_base_conic(curve2)) - { - Comparison_result fac1 = curve1.facing(); - Comparison_result fac2 = curve2.facing(); - - if (fac1 == LARGER && fac2 == SMALLER) - return (LARGER); - else if (fac1 == SMALLER && fac2 == LARGER) - return (SMALLER); - } - - // Otherwise, the two curves do intersect at p_int = ps1[0] = ps2[0]: - // make a decision based on their partial derivatives. - const Point_2& p_int = ps1[0]; - const NT _zero = 0; - - NT slope1_numer, slope1_denom; - NT slope2_numer, slope2_denom; - - curve1.derive_by_x_at (p_int, 1, slope1_numer, slope1_denom); - curve2.derive_by_x_at (p_int, 1, slope2_numer, slope2_denom); - - const bool is_vertical_slope1 = (slope1_denom == _zero); - const bool is_vertical_slope2 = (slope2_denom == _zero); - - if (!is_vertical_slope1 && !is_vertical_slope2) - { - // The two curves have derivatives at p_int: use it to determine which - // one is above the other (the one with a larger slope is below). - Comparison_result slope_res = /*CGAL::compare*/ - eps_compare (TO_APNT(slope1_numer*slope2_denom), - TO_APNT(slope2_numer*slope1_denom)); - - if (slope_res != EQUAL) - return (slope_res); - - // Use the second order derivative. - curve1.derive_by_x_at (p_int, 2, slope1_numer, slope1_denom); - curve2.derive_by_x_at (p_int, 2, slope2_numer, slope2_denom); - - slope_res = /*CGAL::compare*/ - eps_compare (TO_APNT(slope1_numer*slope2_denom), - TO_APNT(slope2_numer*slope1_denom)); - - CGAL_assertion(slope_res != EQUAL); - - return (slope_res); - } - else if (!is_vertical_slope2) - { - // The first curve has a vertical slope at p_int: check whether it is - // facing upwards or downwards and decide accordingly. - Comparison_result fac1 = curve1.facing(); - - CGAL_assertion(fac1 != EQUAL); - return (fac1); - } - else if (!is_vertical_slope1) - { - // The secong curve has a vertical slope at p_int: check whether it is - // facing upwards or downwards and decide accordingly. - Comparison_result fac2 = curve2.facing(); - - CGAL_assertion(fac2 != EQUAL); - return ((fac2 == LARGER) ? SMALLER : LARGER); - } - else - { - // The two curves have vertical slopes at p_int: - // First check whether one is facing up and one down. In this case the - // comparison result is trivial. - Comparison_result fac1 = curve1.facing(); - Comparison_result fac2 = curve2.facing(); - - if (fac1 == LARGER && fac2 == SMALLER) - return (LARGER); - else if (fac1 == SMALLER && fac2 == LARGER) - return (SMALLER); - - // Compute the second order derivative by y and act according to it. - curve1.derive_by_y_at (p_int, 2, slope1_numer, slope1_denom); - curve2.derive_by_y_at (p_int, 2, slope2_numer, slope2_denom); - - Comparison_result slope_res = /*CGAL::compare*/ - eps_compare (TO_APNT(slope1_numer*slope2_denom), - TO_APNT(slope2_numer*slope1_denom)); - - CGAL_assertion(slope_res != EQUAL); - - if (fac1 == LARGER && fac2 == LARGER) - { - // Both are facing up. - return ((slope_res == LARGER) ? SMALLER : LARGER); - } - else - { - // Both are facing down. - return (slope_res); - } - } - - // We should never reach here: - CGAL_assertion(false); - return (EQUAL); - } + // The two curves intersect at ps1[0] = ps2[0] - proceed from here: + return (_curve_compare_at_intersection_right (curve1, curve2, ps1[0])); + } // Check whether the given point is above, under or on the given curve. - Curve_point_status curve_get_point_status (const X_curve_2 & curve, - const Point_2 & p) const + Curve_point_status curve_get_point_status (const X_curve_2& curve, + const Point_2& p) const { CGAL_precondition(is_x_monotone(curve)); @@ -615,23 +388,23 @@ class Arr_conic_traits_2 if (curve.is_vertical_segment()) { if (compare_x (curve.source(), p) != EQUAL) - return (CURVE_NOT_IN_RANGE); + return (CURVE_NOT_IN_RANGE); // In case p has the same x c-ordinate of the vertical segment, compare // it to the segment endpoints to determine its position. if (compare_y (curve.source(), p) == SMALLER && - compare_y (curve.target(), p) == SMALLER) + compare_y (curve.target(), p) == SMALLER) { - return (ABOVE_CURVE); + return (ABOVE_CURVE); } else if (compare_y (curve.source(), p) == LARGER && - compare_y (curve.target(), p) == LARGER) + compare_y (curve.target(), p) == LARGER) { - return (UNDER_CURVE); + return (UNDER_CURVE); } else { - return (ON_CURVE); + return (ON_CURVE); } } @@ -673,8 +446,8 @@ class Arr_conic_traits_2 // clockwise direction from p from c1 to c2. // Notice that all three curves share the same end-point p. bool curve_is_between_cw (const X_curve_2& curve, - const X_curve_2& c1, const X_curve_2& c2, - const Point_2& p) const + const X_curve_2& c1, const X_curve_2& c2, + const Point_2& p) const { CGAL_precondition(is_x_monotone(curve)); CGAL_precondition(is_x_monotone(c1)); @@ -689,7 +462,7 @@ class Arr_conic_traits_2 const Point_2* p1_P = NULL; const Point_2* p2_P = NULL; const Point_2* pv_P = NULL; - int p1_flip = 1, p2_flip = 1, pv_flip = 1; + int p1_flip = 1, p2_flip = 1, pv_flip = 1; if (curve.source().equals(p)) pv_P = &(curve.target()); @@ -717,24 +490,24 @@ class Arr_conic_traits_2 // Make sure no two curves overlap. if (c1.has_same_base_conic (c2) && - (p1_flip * c1.conic().orientation() == - p2_flip * c2.conic().orientation()) && - (c1.contains_point(*p2_P) || - c2.contains_point(*p1_P))) + (p1_flip * c1.conic().orientation() == + p2_flip * c2.conic().orientation()) && + (c1.contains_point(*p2_P) || + c2.contains_point(*p1_P))) return (false); if (c1.has_same_base_conic(curve) && - (p1_flip * c1.conic().orientation() == - pv_flip * curve.conic().orientation()) && - (c1.contains_point(*pv_P) || - curve.contains_point(*p1_P))) + (p1_flip * c1.conic().orientation() == + pv_flip * curve.conic().orientation()) && + (c1.contains_point(*pv_P) || + curve.contains_point(*p1_P))) return (false); if (curve.has_same_base_conic(c2) && - (pv_flip * curve.conic().orientation() == - p2_flip * c2.conic().orientation()) && - (curve.contains_point(*p2_P) || - c2.contains_point(*pv_P))) + (pv_flip * curve.conic().orientation() == + p2_flip * c2.conic().orientation()) && + (curve.contains_point(*p2_P) || + c2.contains_point(*pv_P))) return (false); // Decide whether each arc is defined to the left or to the right of p. @@ -758,82 +531,107 @@ class Arr_conic_traits_2 // c1 is a vertical segment: if (c2_vertical != 0) { - // Both c1 and c2 are vertical segments: - if ((c1_vertical == 1) && (c2_vertical == -1)) - return (!curve_is_left); - else if ((c1_vertical == -1) && (c2_vertical == 1)) - return (curve_is_left); - else - return (false); + if (curve_vertical != 0) + return (false); + + // Both c1 and c2 are vertical segments: + if ((c1_vertical == 1) && (c2_vertical == -1)) + return (!curve_is_left); + else if ((c1_vertical == -1) && (c2_vertical == 1)) + return (curve_is_left); + else + return (false); + } + + if (curve_vertical != 0) + { + // Both c1 and curve are vertical segments: + if ((c1_vertical == 1) && (curve_vertical == -1)) + return (c2_is_left); + else if ((c1_vertical == -1) && (curve_vertical == 1)) + return (!c2_is_left); + else + return (false); } if (c1_vertical == 1) { - // c1 is a vertical segment going up: - if (c2_is_left) - { - return ((!curve_is_left) || - (curve_is_left && - curve_compare_at_x_left (c2, curve, p) == LARGER)); - } - else - { - return (!curve_is_left && - curve_compare_at_x_right (c2, curve, p) == SMALLER); - } + // c1 is a vertical segment going up: + if (c2_is_left) + { + return ((!curve_is_left) || + (curve_is_left && + _curve_compare_at_intersection_left (c2, curve, p) == LARGER)); + } + else + { + return (!curve_is_left && + _curve_compare_at_intersection_right (c2, curve, p) == SMALLER); + } } else { - // c1 is a vertical segment going down: - if (c2_is_left) - { - return (curve_is_left && - curve_compare_at_x_left (c2, curve, p) == LARGER); - } - else - { - return ((curve_is_left) || - (curve_vertical == 1) || - (!curve_is_left && - curve_compare_at_x_right (c2, curve, p) == SMALLER)); - } + // c1 is a vertical segment going down: + if (c2_is_left) + { + return (curve_is_left && + _curve_compare_at_intersection_left (c2, curve, p) == LARGER); + } + else + { + return ((curve_is_left) || + (curve_vertical == 1) || + (!curve_is_left && + _curve_compare_at_intersection_right (c2, curve, p) == SMALLER)); + } } } if (c2_vertical != 0) { + if (curve_vertical != 0) + { + // Both c2 and curve are vertical segments: + if ((c2_vertical == 1) && (curve_vertical == -1)) + return (!c1_is_left); + else if ((c2_vertical == -1) && (curve_vertical == 1)) + return (c1_is_left); + else + return (false); + } + // Only c2 is a vertical segment: if (c2_vertical == 1) { - // c2 is a vertical segment going up: - if (c1_is_left) - { - return (curve_is_left && - curve_compare_at_x_left (c1, curve, p) == SMALLER); - } - else - { - return ((curve_is_left) || - (curve_vertical == -1) || - (!curve_is_left && - curve_compare_at_x_right (c1, curve, p) == LARGER)); - } + // c2 is a vertical segment going up: + if (c1_is_left) + { + return (curve_is_left && + _curve_compare_at_intersection_left (c1, curve, p) == SMALLER); + } + else + { + return ((curve_is_left) || + (curve_vertical == -1) || + (!curve_is_left && + _curve_compare_at_intersection_right (c1, curve, p) == LARGER)); + } } else { - // c2 is a vertical segment going down: - if (c1_is_left) - { - return ((!curve_is_left) || - (curve_vertical == 1) || - (curve_is_left && - curve_compare_at_x_left (c1, curve, p) == SMALLER)); - } - else - { - return ((!curve_is_left) && - curve_compare_at_x_right (c1, curve, p) == LARGER); - } + // c2 is a vertical segment going down: + if (c1_is_left) + { + return ((!curve_is_left) || + (curve_vertical == 1) || + (curve_is_left && + _curve_compare_at_intersection_left (c1, curve, p) == SMALLER)); + } + else + { + return ((!curve_is_left) && + _curve_compare_at_intersection_right (c1, curve, p) == LARGER); + } } } @@ -842,23 +640,23 @@ class Arr_conic_traits_2 // Only curve is a vertical segment: if (curve_vertical == 1) { - // curve is a vertical segment going up: - if (c1_is_left && !c2_is_left) - return (true); - else if (c1_is_left && c2_is_left) - return (curve_compare_at_x_left(c1, c2, p) == LARGER); - else if (!c1_is_left && !c2_is_left) - return (curve_compare_at_x_right(c1, c2, p) == SMALLER); + // curve is a vertical segment going up: + if (c1_is_left && !c2_is_left) + return (true); + else if (c1_is_left && c2_is_left) + return (_curve_compare_at_intersection_left(c1, c2, p) == LARGER); + else if (!c1_is_left && !c2_is_left) + return (_curve_compare_at_intersection_right(c1, c2, p) == SMALLER); } else { - // curve is a vertical segment going down: - if (!c1_is_left && c2_is_left) - return (true); - else if (c1_is_left && c2_is_left) - return (curve_compare_at_x_left(c1, c2, p) == LARGER); - else if (!c1_is_left && !c2_is_left) - return (curve_compare_at_x_right(c1, c2, p) == SMALLER); + // curve is a vertical segment going down: + if (!c1_is_left && c2_is_left) + return (true); + else if (c1_is_left && c2_is_left) + return (_curve_compare_at_intersection_left(c1, c2, p) == LARGER); + else if (!c1_is_left && !c2_is_left) + return (_curve_compare_at_intersection_right(c1, c2, p) == SMALLER); } return (false); @@ -868,51 +666,61 @@ class Arr_conic_traits_2 // Check the following 4 cases: if (c1_is_left && c2_is_left) { + CGAL_assertion(!c1.has_same_base_conic(c2) || + c1.facing() != c2.facing()); + // Case 1: Both c1 and c2 are defined to the left of p. - if (curve_compare_at_x_left (c1, c2, p) == LARGER) + if (_curve_compare_at_intersection_left (c1, c2, p) == LARGER) { - // c1 is above c2: - return (!(curve_compare_at_x_left (c2, curve, p) == SMALLER && - curve_compare_at_x_left (c1, curve, p) == LARGER)); + // c1 is above c2: + return (!curve_is_left || + !(_curve_compare_at_intersection_left (c2, curve, p) == SMALLER && + _curve_compare_at_intersection_left (c1, curve, p) == LARGER)); } else { - // c2 is above c1: - return (curve_compare_at_x_left (c1, curve, p) == SMALLER && - curve_compare_at_x_left (c2, curve, p) == LARGER); + // c2 is above c1: + return (curve_is_left && + _curve_compare_at_intersection_left (c1, curve, p) == SMALLER && + _curve_compare_at_intersection_left (c2, curve, p) == LARGER); } } else if (!c1_is_left && !c2_is_left) { + CGAL_assertion(!c1.has_same_base_conic(c2) || + c1.facing() != c2.facing()); + // Case 2: Both c1 and c2 are defined to the right of p. - if (curve_compare_at_x_right (c1, c2, p) == LARGER) + if (_curve_compare_at_intersection_right (c1, c2, p) == LARGER) { - // c1 is above c2: - return (curve_compare_at_x_right (c2, curve, p) == SMALLER && - curve_compare_at_x_right (c1, curve, p) == LARGER); + // c1 is above c2: + return (!curve_is_left && + _curve_compare_at_intersection_right (c2, curve, p) == SMALLER && + _curve_compare_at_intersection_right (c1, curve, p) == LARGER); } else { - // c2 is above c1: - return (!(curve_compare_at_x_right (c1, curve, p) == SMALLER && - curve_compare_at_x_right (c2, curve, p) == LARGER)); + // c2 is above c1: + return (curve_is_left || + !(_curve_compare_at_intersection_right (c1, curve, p) == SMALLER && + _curve_compare_at_intersection_right (c2, curve, p) == LARGER)); } } else if (c1_is_left && !c2_is_left) { // Case 3: c1 is defined to the left and c2 is to the right of p. if (curve_is_left) - return (curve_compare_at_x_left(c1, curve, p) == SMALLER); + return (_curve_compare_at_intersection_left(c1, curve, p) == SMALLER); else - return (curve_compare_at_x_right(c2, curve, p) == SMALLER); + return (_curve_compare_at_intersection_right(c2, curve, p) == SMALLER); } else if (!c1_is_left && c2_is_left) { // Case 4: c1 is defined to the right and c2 is to the left of p. if (curve_is_left) - return (curve_compare_at_x_left (c2, curve, p) == LARGER); + return (_curve_compare_at_intersection_left (c2, curve, p) == LARGER); else - return (curve_compare_at_x_right(c1, curve, p) == LARGER); + return (_curve_compare_at_intersection_right(c1, curve, p) == LARGER); } // We should never reach here. @@ -920,6 +728,12 @@ class Arr_conic_traits_2 return false; } + // Cehck whether the two points are identical. + bool point_is_same (const Point_2& p1, const Point_2& p2) const + { + return (p1.equals(p2)); + } + // Check whether the two curves are identical. bool curve_is_same (const X_curve_2& curve1, const X_curve_2& curve2) const { @@ -931,15 +745,15 @@ class Arr_conic_traits_2 { // Same orientation: return (curve1.has_same_base_conic(curve2) && - curve1.source().equals(curve2.source()) && - curve1.target().equals(curve2.target())); + curve1.source().equals(curve2.source()) && + curve1.target().equals(curve2.target())); } else { // Check the flip case: return (curve1.has_same_base_conic(curve2) && - curve1.source().equals(curve2.target()) && - curve1.target().equals(curve2.source())); + curve1.source().equals(curve2.target()) && + curve1.target().equals(curve2.source())); } } @@ -954,6 +768,21 @@ class Arr_conic_traits_2 return (curve.target()); } + // Return a point to the left or to the right of p. + Point_2 point_to_left (const Point_2& p) const + { + NT x = CGAL::to_double(p.x()) - 1; + return (Point_2(x , p.y(), + Point_2::User_defined)); + } + + Point_2 point_to_right (const Point_2& p) const + { + NT x = CGAL::to_double(p.x()) + 1; + return (Point_2(x , p.y(), + Point_2::User_defined)); + } + // Reflect a point in y. Point_2 point_reflect_in_y (const Point_2& p) const { @@ -980,6 +809,13 @@ class Arr_conic_traits_2 ////////// Arrangement methods: ////////// + // Change the orientation of the curve (swap the source and the target). + X_curve_2 curve_flip (const X_curve_2& curve) const + { + // Flip the arc. + return (curve.flip()); + } + // Check whether the curve is x-monotone. bool is_x_monotone (const Curve_2& curve) const { @@ -988,7 +824,7 @@ class Arr_conic_traits_2 // Cut the curve to several x-monotone sub-curves. void make_x_monotone (const Curve_2& curve, - std::list& x_curves) const + std::list& x_curves) const { CGAL_precondition(!is_x_monotone(curve)); @@ -1022,57 +858,57 @@ class Arr_conic_traits_2 if (n == 1) { - // Split the arc into two x-monotone sub-curves: one going from the - // arc source to ps[0], and the other from ps[0] to the target. - _curve_split (curve, - sub_curve1, sub_curve2, + // Split the arc into two x-monotone sub-curves: one going from the + // arc source to ps[0], and the other from ps[0] to the target. + _curve_split (curve, + sub_curve1, sub_curve2, ps[0]); - x_curves.push_back(sub_curve1); - x_curves.push_back(sub_curve2); + x_curves.push_back(sub_curve1); + x_curves.push_back(sub_curve2); } else if (n == 2) { - // Split the arc into three x-monotone sub-curves: one going from the - // arc source to ps[0], one from ps[0] to ps[1], and the last one - // from ps[1] to the target. - // Notice that ps[0] and ps[1] might switch places. - X_curve_2 temp; + // Split the arc into three x-monotone sub-curves: one going from the + // arc source to ps[0], one from ps[0] to ps[1], and the last one + // from ps[1] to the target. + // Notice that ps[0] and ps[1] might switch places. + X_curve_2 temp; - _curve_split (curve, - sub_curve1, sub_curve2, + _curve_split (curve, + sub_curve1, sub_curve2, ps[0]); - if (sub_curve2.contains_point(ps[1])) - { - temp = sub_curve2; - _curve_split (temp, - sub_curve2, sub_curve3, - ps[1]); - } - else if (sub_curve1.contains_point(ps[1])) - { - // Actually we switch between ps[0] and ps[1]. - temp = sub_curve1; - sub_curve3 = sub_curve2; - _curve_split (temp, - sub_curve1, sub_curve2, - ps[1]); - } - else - { - // We should never reach here: - CGAL_assertion(false); - } + if (sub_curve2.contains_point(ps[1])) + { + temp = sub_curve2; + _curve_split (temp, + sub_curve2, sub_curve3, + ps[1]); + } + else if (sub_curve1.contains_point(ps[1])) + { + // Actually we switch between ps[0] and ps[1]. + temp = sub_curve1; + sub_curve3 = sub_curve2; + _curve_split (temp, + sub_curve1, sub_curve2, + ps[1]); + } + else + { + // We should never reach here: + CGAL_assertion(false); + } - x_curves.push_back(sub_curve1); - x_curves.push_back(sub_curve2); - x_curves.push_back(sub_curve3); + x_curves.push_back(sub_curve1); + x_curves.push_back(sub_curve2); + x_curves.push_back(sub_curve3); } else { - // We should never reach here: - CGAL_assertion(false); + // We should never reach here: + CGAL_assertion(false); } } @@ -1081,7 +917,7 @@ class Arr_conic_traits_2 // Split the given curve into two sub-curves at the given point. void curve_split (const X_curve_2& curve, - X_curve_2& sub_curve1, X_curve_2& sub_curve2, + X_curve_2& sub_curve1, X_curve_2& sub_curve2, const Point_2& p) const { CGAL_precondition(is_x_monotone(curve)); @@ -1093,8 +929,8 @@ class Arr_conic_traits_2 // Split the curve. _curve_split (curve, - sub_curve1, sub_curve2, - p); + sub_curve1, sub_curve2, + p); return; } @@ -1118,7 +954,7 @@ class Arr_conic_traits_2 // Check if at least one end-point of the overlapping curve is to the // right of p. return (ovlp_arcs[0].source().compare_lex_xy(p) == LARGER || - ovlp_arcs[0].target().compare_lex_xy(p) == LARGER); + ovlp_arcs[0].target().compare_lex_xy(p) == LARGER); } // In case there are no overlaps and the base conics are the same, @@ -1127,17 +963,17 @@ class Arr_conic_traits_2 if (curve1.has_same_base_conic(curve2)) { if ((curve1.source().equals(curve2.source()) || - curve1.source().equals(curve2.target())) && - curve1.source().compare_lex_xy(p) == LARGER) + curve1.source().equals(curve2.target())) && + curve1.source().compare_lex_xy(p) == LARGER) { - return (true); + return (true); } if ((curve1.target().equals(curve2.source()) || - curve1.target().equals(curve2.target())) && - curve1.target().compare_lex_xy(p) == LARGER) + curve1.target().equals(curve2.target())) && + curve1.target().compare_lex_xy(p) == LARGER) { - return (true); + return (true); } // No overlaps at all: the two curves do not intersect. @@ -1157,7 +993,7 @@ class Arr_conic_traits_2 for (int i = 0; i < n; i++) { if (ps[i].compare_lex_xy(p) == LARGER) - return (true); + return (true); } return (false); @@ -1185,7 +1021,7 @@ class Arr_conic_traits_2 // Check if at least one end-point of the overlapping curve is to the // left of p. return (ovlp_arcs[0].source().compare_lex_xy(p) == SMALLER || - ovlp_arcs[0].target().compare_lex_xy(p) == SMALLER); + ovlp_arcs[0].target().compare_lex_xy(p) == SMALLER); } // In case there are no overlaps and the base conics are the same, @@ -1194,17 +1030,17 @@ class Arr_conic_traits_2 if (curve1.has_same_base_conic(curve2)) { if ((curve1.source().equals(curve2.source()) || - curve1.source().equals(curve2.target())) && - curve1.source().compare_lex_xy(p) == SMALLER) + curve1.source().equals(curve2.target())) && + curve1.source().compare_lex_xy(p) == SMALLER) { - return (true); + return (true); } if ((curve1.target().equals(curve2.source()) || - curve1.target().equals(curve2.target())) && - curve1.target().compare_lex_xy(p) == SMALLER) + curve1.target().equals(curve2.target())) && + curve1.target().compare_lex_xy(p) == SMALLER) { - return (true); + return (true); } // No overlaps at all: the two curves do not intersect. @@ -1224,7 +1060,7 @@ class Arr_conic_traits_2 for (int i = 0; i < n; i++) { if (ps[i].compare_lex_xy(p) == SMALLER) - return (true); + return (true); } return (false); @@ -1236,8 +1072,8 @@ class Arr_conic_traits_2 // In case of an overlap, p1 and p2 are the source and destination of the // overlapping curve. Otherwise p1=p2 is the calculated intersection point. bool nearest_intersection_to_right (const X_curve_2& curve1, - const X_curve_2& curve2, - const Point_2& p, + const X_curve_2& curve2, + const Point_2& p, Point_2& p1, Point_2& p2) const { @@ -1257,33 +1093,33 @@ class Arr_conic_traits_2 Point_2 ovlp_target = ovlp_arcs[0].target(); if (ovlp_source.compare_lex_xy(p) == LARGER && - ovlp_target.compare_lex_xy(p) == LARGER) + ovlp_target.compare_lex_xy(p) == LARGER) { - // The entire overlapping arc is to the right of p: - p1 = ovlp_source; - p2 = ovlp_target; - return (true); + // The entire overlapping arc is to the right of p: + p1 = ovlp_source; + p2 = ovlp_target; + return (true); } else if (ovlp_source.compare_lex_xy(p) != LARGER && - ovlp_target.compare_lex_xy(p) == LARGER) + ovlp_target.compare_lex_xy(p) == LARGER) { - // The source is to the left of p, and the traget is to its right. - p1 = p; - p2 = ovlp_target; - return (true); + // The source is to the left of p, and the traget is to its right. + p1 = p; + p2 = ovlp_target; + return (true); } else if (ovlp_source.compare_lex_xy(p) == LARGER && - ovlp_target.compare_lex_xy(p) != LARGER) + ovlp_target.compare_lex_xy(p) != LARGER) { - // The source is to the right of p, and the traget is to its left. - p1 = ovlp_source; - p2 = p; - return (true); + // The source is to the right of p, and the traget is to its left. + p1 = ovlp_source; + p2 = p; + return (true); } else { - // The entire overlapping arc is to the left of p: - return (false); + // The entire overlapping arc is to the left of p: + return (false); } } @@ -1295,33 +1131,33 @@ class Arr_conic_traits_2 const Point_2 *nearest_end_P = NULL; if ((curve1.source().equals(curve2.source()) || - curve1.source().equals(curve2.target())) && - curve1.source().compare_lex_xy(p) == LARGER) + curve1.source().equals(curve2.target())) && + curve1.source().compare_lex_xy(p) == LARGER) { - nearest_end_P = &(curve1.source()); + nearest_end_P = &(curve1.source()); } if ((curve1.target().equals(curve2.source()) || - curve1.target().equals(curve2.target())) && - curve1.target().compare_lex_xy(p) == LARGER) + curve1.target().equals(curve2.target())) && + curve1.target().compare_lex_xy(p) == LARGER) { - if (nearest_end_P == NULL || - nearest_end_P->compare_lex_xy (curve1.target()) == LARGER) - { - nearest_end_P = &(curve1.target()); - } + if (nearest_end_P == NULL || + nearest_end_P->compare_lex_xy (curve1.target()) == LARGER) + { + nearest_end_P = &(curve1.target()); + } } if (nearest_end_P != NULL) { - // A common end point was found: - p1 = p2 = *nearest_end_P; - return (true); + // A common end point was found: + p1 = p2 = *nearest_end_P; + return (true); } else { - // No intersection: - return (false); + // No intersection: + return (false); } } @@ -1341,12 +1177,12 @@ class Arr_conic_traits_2 // Check if the current point is to the right of p. if (ps[i].compare_lex_xy(p) == LARGER) { - // Compare with the nearest point so far. - if (nearest_inter_P == NULL || - nearest_inter_P->compare_lex_xy (ps[i]) == LARGER) - { - nearest_inter_P = &(ps[i]); - } + // Compare with the nearest point so far. + if (nearest_inter_P == NULL || + nearest_inter_P->compare_lex_xy (ps[i]) == LARGER) + { + nearest_inter_P = &(ps[i]); + } } } @@ -1368,10 +1204,10 @@ class Arr_conic_traits_2 // In case of an overlap, p1 and p2 are the source and destination of the // overlapping curve. Otherwise p1=p2 is the calculated intersection point. bool nearest_intersection_to_left (const X_curve_2& curve1, - const X_curve_2& curve2, - const Point_2& p, - Point_2& p1, - Point_2& p2) const + const X_curve_2& curve2, + const Point_2& p, + Point_2& p1, + Point_2& p2) const { CGAL_precondition(is_x_monotone(curve1)); CGAL_precondition(is_x_monotone(curve2)); @@ -1389,33 +1225,33 @@ class Arr_conic_traits_2 Point_2 ovlp_target = ovlp_arcs[0].target(); if (ovlp_source.compare_lex_xy(p) == SMALLER && - ovlp_target.compare_lex_xy(p) == SMALLER) + ovlp_target.compare_lex_xy(p) == SMALLER) { - // The entire overlapping arc is to the left of p: - p1 = ovlp_source; - p2 = ovlp_target; - return (true); + // The entire overlapping arc is to the left of p: + p1 = ovlp_source; + p2 = ovlp_target; + return (true); } else if (ovlp_source.compare_lex_xy(p) != SMALLER && - ovlp_target.compare_lex_xy(p) == SMALLER) + ovlp_target.compare_lex_xy(p) == SMALLER) { - // The source is to the right of p, and the traget is to its left. - p1 = p; - p2 = ovlp_target; - return (true); + // The source is to the right of p, and the traget is to its left. + p1 = p; + p2 = ovlp_target; + return (true); } else if (ovlp_source.compare_lex_xy(p) == SMALLER && - ovlp_target.compare_lex_xy(p) != SMALLER) + ovlp_target.compare_lex_xy(p) != SMALLER) { - // The source is to the left of p, and the traget is to its right. - p1 = ovlp_source; - p2 = p; - return (true); + // The source is to the left of p, and the traget is to its right. + p1 = ovlp_source; + p2 = p; + return (true); } else { - // The entire overlapping arc is to the right of p: - return (false); + // The entire overlapping arc is to the right of p: + return (false); } } @@ -1427,33 +1263,33 @@ class Arr_conic_traits_2 const Point_2 *nearest_end_P = NULL; if ((curve1.source().equals(curve2.source()) || - curve1.source().equals(curve2.target())) && - curve1.source().compare_lex_xy(p) == SMALLER) + curve1.source().equals(curve2.target())) && + curve1.source().compare_lex_xy(p) == SMALLER) { - nearest_end_P = &(curve1.source()); + nearest_end_P = &(curve1.source()); } if ((curve1.target().equals(curve2.source()) || - curve1.target().equals(curve2.target())) && - curve1.target().compare_lex_xy(p) == SMALLER) + curve1.target().equals(curve2.target())) && + curve1.target().compare_lex_xy(p) == SMALLER) { - if (nearest_end_P == NULL || - nearest_end_P->compare_lex_xy (curve1.target()) == SMALLER) - { - nearest_end_P = &(curve1.target()); - } + if (nearest_end_P == NULL || + nearest_end_P->compare_lex_xy (curve1.target()) == SMALLER) + { + nearest_end_P = &(curve1.target()); + } } if (nearest_end_P != NULL) { - // A common end point was found: - p1 = p2 = *nearest_end_P; - return (true); + // A common end point was found: + p1 = p2 = *nearest_end_P; + return (true); } else { - // No intersection: - return (false); + // No intersection: + return (false); } } @@ -1473,12 +1309,12 @@ class Arr_conic_traits_2 // Check if the current point is to the right of p. if (ps[i].compare_lex_xy(p) == SMALLER) { - // Compare with the nearest point so far. - if (nearest_inter_P == NULL || - nearest_inter_P->compare_lex_xy (ps[i]) == SMALLER) - { - nearest_inter_P = &(ps[i]); - } + // Compare with the nearest point so far. + if (nearest_inter_P == NULL || + nearest_inter_P->compare_lex_xy (ps[i]) == SMALLER) + { + nearest_inter_P = &(ps[i]); + } } } @@ -1512,16 +1348,287 @@ class Arr_conic_traits_2 // Split the given curve into two sub-curves at the given point. // Since this is a private function, there are no preconditions. void _curve_split (const X_curve_2& curve, - X_curve_2& sub_curve1, X_curve_2& sub_curve2, - const Point_2& p) const + X_curve_2& sub_curve1, X_curve_2& sub_curve2, + const Point_2& p) const { + CGAL_precondition(! p.equals(curve.source())); + CGAL_precondition(! p.equals(curve.target())); + // Split the curve to source->p and p->target. - sub_curve1 = X_curve_2(curve, curve.source(), p, false); - sub_curve2 = X_curve_2(curve, p, curve.target(), false); + sub_curve1 = X_curve_2 (curve, curve.source(), p, false); + sub_curve2 = X_curve_2 (curve, p, curve.target(), false); + + CGAL_assertion(sub_curve1.is_x_monotone()); + CGAL_assertion(sub_curve2.is_x_monotone()); return; } + // Decide wether curve1 is above, below or equal to curve2 immediately to + // the left of p_int, which is assumed to be an intersection of both curves. + // Furthermore, the two curves are assumed to be defined to p_int's left. + Comparison_result _curve_compare_at_intersection_left + (const X_curve_2& curve1, + const X_curve_2& curve2, + const Point_2& p_int) const + { + // In case one arc is facing upwards and another facing downwards, it is + // very clear that the one facing upward is above the one facing downwards. + if (curve1.has_same_base_conic(curve2)) + { + Comparison_result fac1 = curve1.facing(); + Comparison_result fac2 = curve2.facing(); + + if (fac1 == LARGER && fac2 == SMALLER) + return (LARGER); + else if (fac1 == SMALLER && fac2 == LARGER) + return (SMALLER); + } + + // Otherwise, the two curves do intersect at p_int = ps1[0] = ps2[0]: + // make a decision based on their partial derivatives. + //const Point_2& p_int = ps1[0]; + const NT _zero = 0; + + NT slope1_numer, slope1_denom; + NT slope2_numer, slope2_denom; + + curve1.derive_by_x_at (p_int, 1, slope1_numer, slope1_denom); + curve2.derive_by_x_at (p_int, 1, slope2_numer, slope2_denom); + + /*const*/ bool is_vertical_slope1 = (slope1_denom == _zero); + /*const*/ bool is_vertical_slope2 = (slope2_denom == _zero); + + // RWRW - CHECK THIS !!! + if (! is_vertical_slope1 && + p_int.is_approximate() && + eps_compare(TO_APNT(slope1_denom), 0) == EQUAL) + { + is_vertical_slope1 = true; + } + if (! is_vertical_slope2 && + p_int.is_approximate() && + eps_compare(TO_APNT(slope2_denom), 0) == EQUAL) + { + is_vertical_slope2 = true; + } + // TO HERE ... + + if (!is_vertical_slope1 && !is_vertical_slope2) + { + // The two curves have derivatives at p_int: use it to determine which + // one is above the other (the one with a smaller slope in above). + + // RWRW - CHECK THIS !!! + Comparison_result slope_res = + (!p_int.is_approximate()) ? + CGAL::compare(slope2_numer*slope1_denom, slope1_numer*slope2_denom) : + eps_compare(TO_APNT(slope2_numer*slope1_denom), + TO_APNT(slope1_numer*slope2_denom)); + + if (slope_res != EQUAL) + return (slope_res); + + // Use the second order derivative. + curve1.derive_by_x_at (p_int, 2, slope1_numer, slope1_denom); + curve2.derive_by_x_at (p_int, 2, slope2_numer, slope2_denom); + + slope_res = CGAL::compare (slope2_numer*slope1_denom, + slope1_numer*slope2_denom); + + CGAL_assertion(slope_res != EQUAL); + + return ((slope_res == LARGER) ? SMALLER : LARGER); + } + else if (!is_vertical_slope2) + { + // The first curve has a vertical slope at p_int: check whether it is + // facing upwards or downwards and decide accordingly. + Comparison_result fac1 = curve1.facing(); + + CGAL_assertion(fac1 != EQUAL); + return (fac1); + } + else if (!is_vertical_slope1) + { + // The second curve has a vertical slope at p_int: check whether it is + // facing upwards or downwards and decide accordingly. + Comparison_result fac2 = curve2.facing(); + + CGAL_assertion(fac2 != EQUAL); + return ((fac2 == LARGER) ? SMALLER : LARGER); + } + else + { + // The two curves have vertical slopes at p_int: + // First check whether one is facing up and one down. In this case the + // comparison result is trivial. + Comparison_result fac1 = curve1.facing(); + Comparison_result fac2 = curve2.facing(); + + if (fac1 == LARGER && fac2 == SMALLER) + return (LARGER); + else if (fac1 == SMALLER && fac2 == LARGER) + return (SMALLER); + + // Compute the second order derivative by y and act according to it. + curve1.derive_by_y_at (p_int, 2, slope1_numer, slope1_denom); + curve2.derive_by_y_at (p_int, 2, slope2_numer, slope2_denom); + + Comparison_result slope_res = CGAL::compare(slope2_numer*slope1_denom, + slope1_numer*slope2_denom); + + CGAL_assertion(slope_res != EQUAL); + + if (fac1 == LARGER && fac2 == LARGER) + { + // Both are facing up. + return ((slope_res == LARGER) ? SMALLER : LARGER); + } + else + { + // Both are facing down. + return (slope_res); + } + } + + // We should never reach here: + CGAL_assertion(false); + return (EQUAL); + } + + // Decide wether curve1 is above, below or equal to curve2 immediately to + // the right of p_int, which is assumed to be an intersection of both curves. + // Furthermore, the two curves are assumed to be defined to p_int's right. + Comparison_result _curve_compare_at_intersection_right + (const X_curve_2& curve1, + const X_curve_2& curve2, + const Point_2& p_int) const + { + // In case one arc is facing upwards and another facing downwards, it is + // very clear that the one facing upward is above the one facing downwards. + if (curve1.has_same_base_conic(curve2)) + { + Comparison_result fac1 = curve1.facing(); + Comparison_result fac2 = curve2.facing(); + + if (fac1 == LARGER && fac2 == SMALLER) + return (LARGER); + else if (fac1 == SMALLER && fac2 == LARGER) + return (SMALLER); + } + + // Otherwise, the two curves do intersect at p_int = ps1[0] = ps2[0]: + // make a decision based on their partial derivatives. + //const Point_2& p_int = ps1[0]; + const NT _zero = 0; + + NT slope1_numer, slope1_denom; + NT slope2_numer, slope2_denom; + + curve1.derive_by_x_at (p_int, 1, slope1_numer, slope1_denom); + curve2.derive_by_x_at (p_int, 1, slope2_numer, slope2_denom); + + /*const*/ bool is_vertical_slope1 = (slope1_denom == _zero); + /*const*/ bool is_vertical_slope2 = (slope2_denom == _zero); + + // RWRW - CHECK THIS !!! + if (! is_vertical_slope1 && + p_int.is_approximate() && + eps_compare(TO_APNT(slope1_denom), 0) == EQUAL) + { + is_vertical_slope1 = true; + } + if (! is_vertical_slope2 && + p_int.is_approximate() && + eps_compare(TO_APNT(slope2_denom), 0) == EQUAL) + { + is_vertical_slope2 = true; + } + // TO HERE ... + + if (!is_vertical_slope1 && !is_vertical_slope2) + { + // The two curves have derivatives at p_int: use it to determine which + // one is above the other (the one with a larger slope is below). + // RWRW - CHECK THIS !!! + Comparison_result slope_res = + (! p_int.is_approximate()) ? + CGAL::compare (slope1_numer*slope2_denom, slope2_numer*slope1_denom) : + eps_compare (TO_APNT(slope1_numer*slope2_denom), + TO_APNT(slope2_numer*slope1_denom)); + + if (slope_res != EQUAL) + return (slope_res); + + // Use the second order derivative. + curve1.derive_by_x_at (p_int, 2, slope1_numer, slope1_denom); + curve2.derive_by_x_at (p_int, 2, slope2_numer, slope2_denom); + + slope_res = CGAL::compare (slope1_numer*slope2_denom, + slope2_numer*slope1_denom); + + CGAL_assertion(slope_res != EQUAL); + + return (slope_res); + } + else if (!is_vertical_slope2) + { + // The first curve has a vertical slope at p_int: check whether it is + // facing upwards or downwards and decide accordingly. + Comparison_result fac1 = curve1.facing(); + + CGAL_assertion(fac1 != EQUAL); + return (fac1); + } + else if (!is_vertical_slope1) + { + // The second curve has a vertical slope at p_int: check whether it is + // facing upwards or downwards and decide accordingly. + Comparison_result fac2 = curve2.facing(); + + CGAL_assertion(fac2 != EQUAL); + return ((fac2 == LARGER) ? SMALLER : LARGER); + } + else + { + // The two curves have vertical slopes at p_int: + // First check whether one is facing up and one down. In this case the + // comparison result is trivial. + Comparison_result fac1 = curve1.facing(); + Comparison_result fac2 = curve2.facing(); + + if (fac1 == LARGER && fac2 == SMALLER) + return (LARGER); + else if (fac1 == SMALLER && fac2 == LARGER) + return (SMALLER); + + // Compute the second order derivative by y and act according to it. + curve1.derive_by_y_at (p_int, 2, slope1_numer, slope1_denom); + curve2.derive_by_y_at (p_int, 2, slope2_numer, slope2_denom); + + Comparison_result slope_res = CGAL::compare (slope1_numer*slope2_denom, + slope2_numer*slope1_denom); + + CGAL_assertion(slope_res != EQUAL); + + if (fac1 == LARGER && fac2 == LARGER) + { + // Both are facing up. + return ((slope_res == LARGER) ? SMALLER : LARGER); + } + else + { + // Both are facing down. + return (slope_res); + } + } + + // We should never reach here: + CGAL_assertion(false); + return (EQUAL); + } + }; CGAL_END_NAMESPACE diff --git a/Packages/Arrangement/include/CGAL/Arrangement_2/Conic_arc_2.h b/Packages/Arrangement/include/CGAL/Arrangement_2/Conic_arc_2.h new file mode 100644 index 00000000000..b3fa1797b0b --- /dev/null +++ b/Packages/Arrangement/include/CGAL/Arrangement_2/Conic_arc_2.h @@ -0,0 +1,2853 @@ +#ifndef CGAL_CONIC_ARC_2_H +#define CGAL_CONIC_ARC_2_H + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#define CGAL_CONIC_ARC_USE_BOUNDING_BOX +#define CGAL_CONIC_ARC_USE_CACHING +#define CGAL_CONIC_ARC_USE_FILTER + +CGAL_BEGIN_NAMESPACE + +enum +{ + REFLECT_IN_X = 1, + REFLECT_IN_Y = 2, + REFLECTION_FACTOR = 4 +}; + +// ---------------------------------------------------------------------------- +// Representation of a conic curve. +// + +static int _conics_count = 0; +template class Arr_conic_traits; + +template +class Conic_arc_2 +{ + friend class Arr_conic_traits; + + public: + + typedef Cartesian Kernel; + typedef Point_2_ex Point_2; + typedef Conic_2 Conic_2; + typedef Circle_2 Circle_2; + typedef Segment_2 Segment_2; + + private: + + enum + { + DEGREE_0 = 0, + DEGREE_1 = 1, + DEGREE_2 = 2, + DEGREE_MASK = 1 + 2, + FULL_CONIC = 4, + X_MONOTONE = 8, + X_MON_UNDEFINED = 8 + 16, + FACING_UP = 32, + FACING_DOWN = 64, + FACING_MASK = 32 + 64, + IS_CIRCLE = 128, + IS_HYPERBOLA = 256 + }; + + Conic_2 _conic; // The conic that contains the arc. + int _conic_id; // The id of the conic. + Point_2 _source; // The source of the arc. + Point_2 _target; // The target of the arc. + int _info; // A bit array with extra information: + // Bit 0 & 1 - The degree of the conic + // (either 1 or 2). + // Bit 2 - Whether the arc is a full conic. + // Bit 3 & 4 - Whether the arc is x-monotone + // (00, 01 or 11 - undefined). + // Bit 5 & 6 - Indicate whether the arc is + // facing upwards or downwards (for + // x-monotone curves of degree 2). + // Bit 7 - Is the underlying curve a circle. + // Bit 8 - Is the underlying curve a hyperbola. +#ifdef CGAL_CONIC_ARC_USE_BOUNDING_BOX + Bbox_2 _bbox; // A bounding box for the arc. +#endif + + // For arcs whose base is a hyperbola we store the axis (a*x + b*y + c = 0) + // which separates the two bracnes of the hyperbola. We also store the side + // (-1 or 1) that the arc occupies. + struct Hyperbolic_arc_data + { + NT a; + NT b; + NT c; + int side; + }; + + // In case of a circle, is is convinient to store its center and radius. + struct Circular_arc_data + { + NT x0; + NT y0; + NT r; + }; + +#ifdef CGAL_CONIC_ARC_USE_CACHING + struct Intersections + { + int id1; + int id2; + int n_points; + Point_2 ps[4]; + }; +#endif + + union + { + Hyperbolic_arc_data *hyper_P; + Circular_arc_data *circ_P; + } _data; + + // Produce a unique id for a new conic. + int _get_new_conic_id () + { + _conics_count++; + return (REFLECTION_FACTOR * _conics_count); + } + + // Private constructor. + Conic_arc_2 (const Conic_arc_2& arc, + const Point_2& source, const Point_2& target, + const bool& is_full) : + _conic(arc._conic), + _conic_id(arc._conic_id), + _source(source), + _target(target) + { + CGAL_precondition(is_full || ! source.equals(target) ); + + _info = (arc._info & DEGREE_MASK) | + (is_full ? FULL_CONIC : 0) | + X_MON_UNDEFINED; + + // Check whether the conic is x-monotone. + if (is_x_monotone()) + { + _info = (_info & ~X_MON_UNDEFINED) | X_MONOTONE; + + // Check whether the facing information is set for the orginating arc. + // If it is, just copy it - otherwise calculate it if the degree is 2. + Comparison_result facing_res = arc.facing(); + + if (facing_res != EQUAL) + _info = _info | (facing_res == LARGER ? FACING_UP : FACING_DOWN); + else if ((_info & DEGREE_MASK) == DEGREE_2) + _set_facing(); + } + else + { + _info = (_info & ~X_MON_UNDEFINED); + } + + // Copy the hyperbolic or circular data, if necessary. + if ((arc._info & IS_HYPERBOLA) != 0) + { + _info = _info | IS_HYPERBOLA; + _data.hyper_P = new Hyperbolic_arc_data (*(arc._data.hyper_P)); + } + else if ((arc._info & IS_CIRCLE) != 0) + { + _info = _info | IS_CIRCLE; + _data.circ_P = new Circular_arc_data (*(arc._data.circ_P)); + } + else + { + _data.hyper_P = NULL; + } + +#ifdef CGAL_CONIC_ARC_USE_BOUNDING_BOX + _bbox = bounding_box(); // Compute the bounding box. +#endif + } + + public: + + // Default constructor. + Conic_arc_2 () : + _conic_id(0), + _info(X_MON_UNDEFINED) + { + _data.hyper_P = NULL; + } + + // Copy constructor. + Conic_arc_2 (const Conic_arc_2& arc) : + _conic(arc._conic), + _conic_id(arc._conic_id), + _source(arc._source), + _target(arc._target), + _info(arc._info) + { + // Copy the hyperbolic or circular data, if necessary. + if ((arc._info & IS_HYPERBOLA) != 0) + { + _data.hyper_P = new Hyperbolic_arc_data (*(arc._data.hyper_P)); + } + else if ((arc._info & IS_CIRCLE) != 0) + { + _data.circ_P = new Circular_arc_data (*(arc._data.circ_P)); + } + else + { + _data.hyper_P = NULL; + } + +#ifdef CGAL_CONIC_ARC_USE_BOUNDING_BOX + _bbox = arc._bbox; // Copy the bounding box. +#endif + } + + // Construct a conic arc which lies on the conic: + // C: r*x^2 + s*y^2 + t*xy + u*x + v*y + w = 0 + // The source and the target must be on the conic boundary and must + // not be the same. + Conic_arc_2 (const NT& r, const NT& s, const NT& t, + const NT& u, const NT& v, const NT& w, + const Point_2& source, const Point_2& target) : + _conic_id(0), + _source(source), + _target(target), + _info(X_MON_UNDEFINED) + { + // Create a conic from the coefficients. + Conic_2 conic (r, s, t, u, v, w); + + // Make sure the conic contains the two end-points on its boundary. + CGAL_precondition(conic.has_on_boundary(source)); + CGAL_precondition(conic.has_on_boundary(target)); + + // Make sure that the source and the taget are not the same. + CGAL_precondition(source != target); + + // Set the arc properties. + _set (conic, source, target); + } + + // Construct a conic arc which is the full conic: + // C: r*x^2 + s*y^2 + t*xy + u*x + v*y + w = 0 + // The conic C must be an ellipse (so 4rs - t^2 > 0). + Conic_arc_2 (const NT& r, const NT& s, const NT& t, + const NT& u, const NT& v, const NT& w) + { + // Make sure the given curve is an ellipse. + CGAL_precondition(4*r*s - t*t > 0); + + // Create a conic from the coefficients. + Conic_2 conic (r, s, t, u, v, w); + + // Set the arc to be the full conic. + _set_full (conic); + } + + // Construct a conic arc from the given line segment. + Conic_arc_2 (const Segment_2& segment) + { + // Use the source and target of the given segment. + Point_2 source = Point_2(segment.source().x(), segment.source().y()); + Point_2 target = Point_2(segment.target().x(), segment.target().y()); + Conic_2 conic; + static const NT _zero = 0; + static const NT _one = 1; + + // Make sure that the source and the taget are not the same. + CGAL_precondition(source != target); + + // The supporting conic is r=s=t=0, and u*x + v*y + w = 0 should hold + // for both the source (x1,y1) and the target (x2, y2). + if (source.x() == target.x()) + { + // The supporting conic is a vertical line, of the form x = CONST. + conic.set (_zero, _zero, _zero, // r = s = t = 0 + _one, // u = 1 + _zero, // v = 0 + -source.x()); // w = -CONST + } + else + { + // The supporting line is A*x + B*y + C = 0, where: + // + // A = y2 - y1, B = x1 - x2, C = x2*y1 - x1*y2 + // + const NT A = (target.y() - source.y()); + const NT B = (source.x() - target.x()); + const NT C = (target.x()*source.y() - source.x()*target.y()); + + // Now we can set: + conic.set (_zero, _zero, _zero, // r = s = t = 0 + A, // u = A + B, // v = B + C); // w = C + } + + // Set the arc properties. + _set (conic, source, target); + } + + // Set a circular arc that lies on the given circle with the given + // end-point. + // Note that the orientation of the input circle is preserved. + Conic_arc_2 (const Circle_2& circle, + const Point_2& source, const Point_2& target) + { + // Make sure the circle contains the two endpoints on its boundary. + CGAL_precondition(circle.has_on_boundary(source)); + CGAL_precondition(circle.has_on_boundary(target)); + + // Make sure that the source and the taget are not the same. + CGAL_precondition(source != target); + + // Produce the correponding conic: if the circle centre is (x0,y0) + // and it radius is r, that its equation is: + // x^2 + y^2 - 2*x0*x - 2*y0*y + (x0^2 + y0^2 - r^2) = 0 + // Since this equation describes a curve with a negative orientation, + // we multiply it by -1 if necessary to preserve the original orientation + // of the input circle. + static const NT _zero = 0; + static const NT _one = 1; + static const NT _minus_one = -1; + static const NT _two = 2; + static const NT _minus_two = -2; + const NT x0 = circle.center().x(); + const NT y0 = circle.center().y(); + const NT r_squared = circle.squared_radius(); + Conic_2 conic; + + if (circle.orientation() == CGAL::COUNTERCLOCKWISE) + { + conic.set (_minus_one, _minus_one, // r = s = -1 + _zero, // t = 0 + _two*x0, + _two*y0, + r_squared - x0*x0 - y0*y0); + } + else + { + conic.set (_one, _one, // r = s = 1 + _zero, // t = 0 + _minus_two*x0, + _minus_two*y0, + x0*x0 + y0*y0 - r_squared); + } + + // Prepare the auxiliary data structure. + Circular_arc_data *circ_data_P = new Circular_arc_data; + + circ_data_P->x0 = x0; + circ_data_P->y0 = y0; + circ_data_P->r = CGAL::sqrt(circle.squared_radius()); + + // Set the arc properties. + _set (conic, source, target, circ_data_P); + } + + // Construct an arc from the full circle. + Conic_arc_2 (const Circle_2& circle) + { + // Produce the correponding conic: if the circle centre is (x0,y0) + // and it radius is r, that its equation is: + // x^2 + y^2 - 2*x0*x - 2*y0*y + (x0^2 + y0^2 - r^2) = 0 + static const NT _zero = 0; + static const NT _one = 1; + static const NT _minus_two = -2; + const NT x0 = circle.center().x(); + const NT y0 = circle.center().y(); + const NT r_squared = circle.squared_radius(); + Conic_2 conic (_one, _one, // r = s = 1 + _zero, // t = 0 + _minus_two*x0, + _minus_two*y0, + x0*x0 + y0*y0 - r_squared); + + // Prepare the auxiliary data structure. + Circular_arc_data *circ_data_P = new Circular_arc_data; + + circ_data_P->x0 = x0; + circ_data_P->y0 = y0; + circ_data_P->r = CGAL::sqrt(circle.squared_radius()); + + // Set the arc to be the full conic. + _set_full (conic, circ_data_P); + } + + // Construct a conic arc which lies on the conic: + // C: r*x^2 + s*y^2 + t*xy + u*x + v*y + w = 0 + // The source and the target are specified by the intersection of the + // conic with: + // C_1: r_1*x^2 + s_1*y^2 + t_1*xy + u_1*x + v_1*y + w_1 = 0 + // C_2: r_2*x^2 + s_2*y^2 + t_2*xy + u_2*x + v_2*y + w_2 = 0 + // The user must also specify the source and the target with approximated + // coordinates. The actual intersection points that best fits the source + // (or the target) will be selected. + Conic_arc_2 (const NT& r, const NT& s, const NT& t, + const NT& u, const NT& v, const NT& w, + const Point_2& app_source, + const NT& r_1, const NT& s_1, const NT& t_1, + const NT& u_1, const NT& v_1, const NT& w_1, + const Point_2& app_target, + const NT& r_2, const NT& s_2, const NT& t_2, + const NT& u_2, const NT& v_2, const NT& w_2) : + _conic_id(0), + _info(X_MON_UNDEFINED) + { + // Create a conic from the given coefficients. + Conic_2 conic (r, s, t, u, v, w); + + _conic = conic; + if (r == 0 && s == 0 && t == 0) + _info = _info | DEGREE_1; + else + _info = _info | DEGREE_2; + + // Create fictitious arcs for the source and target computation. + Conic_2 conic_s (r_1, s_1, t_1, u_1, v_1, w_1); + Conic_2 conic_t (r_2, s_2, t_2, u_2, v_2, w_2); + Conic_arc_2 arc_s; + Conic_arc_2 arc_t; + + arc_s._conic = conic_s; + if (r_1 == 0 && s_1 == 0 && t_1 == 0) + arc_s._info = arc_s._info | DEGREE_1; + else + arc_s._info = arc_s._info | DEGREE_2; + + arc_t._conic = conic_t; + if (r_2 == 0 && s_2 == 0 && t_2 == 0) + arc_t._info = arc_t._info | DEGREE_1; + else + arc_t._info = arc_t._info | DEGREE_2; + + // Compute the source and the target. + Point_2 source; + Point_2 target; + Conic_arc_2 *arc_P; + const Point_2 *app_P; + Point_2 *end_P; + int i, j; + Point_2 ipts[4]; // The intersection points. + int n_points = 0; // Their number. + NT xs[4]; + int x_mults[4]; + int n_xs; // Total number of x co-ordinates. + int n_approx_xs; // Number of approximate x co-ordinates. + NT ys[4]; + int y_mults[4]; + int n_ys; // Total number of y co-ordinates. + int n_approx_ys; // Number of approximate y co-ordinates. + int x_deg[2]; // The generating polynomial for the x values. + std::vector x_coeffs[2]; + int y_deg[2]; // The generating polynomial for the y values. + std::vector y_coeffs[2]; + APNT dist, best_dist = 0; + int index; + + for (i = 0; i < 2; i++) + { + if (i == 0) + { + arc_P = &arc_s; + app_P = &app_source; + end_P = &source; + } + else + { + arc_P = &arc_t; + app_P = &app_target; + end_P = ⌖ + } + + n_xs = _x_coordinates_of_intersections_with (*arc_P, + xs, x_mults, + n_approx_xs, + x_deg[i], + x_coeffs[i]); + + n_ys = _y_coordinates_of_intersections_with (*arc_P, + ys, y_mults, + n_approx_ys, + y_deg[i], + y_coeffs[i]); + + n_points = _pair_intersection_points (*arc_P, + n_xs, + xs, x_mults, + n_approx_xs, + n_ys, + ys, y_mults, + n_approx_ys, + ipts); + + // Match the best point. + index = -1; + for (j = 0; j < n_points; j++) + { + dist = TO_APNT((ipts[j].x() - app_P->x())*(ipts[j].x() - app_P->x()) + + (ipts[j].y() - app_P->y())*(ipts[j].y() - app_P->y())); + + if (index == -1 || dist < best_dist) + { + index = j; + best_dist = dist; + } + } + + CGAL_assertion(index != -1); + *end_P = ipts[index]; + } + + // Set the arc. + _set (conic, source, target); + + // Make sure the end-point carry all information about their background + // polynomials. + _source.attach_polynomials (x_deg[0], x_coeffs[0], + y_deg[0], y_coeffs[0]); + + _target.attach_polynomials (x_deg[1], x_coeffs[1], + y_deg[1], y_coeffs[1]); + } + + // Destructor. + virtual ~Conic_arc_2 () + { + if ((_info & IS_HYPERBOLA) != 0) + delete _data.hyper_P; + else if ((_info & IS_CIRCLE) != 0) + delete _data.circ_P; + _data.hyper_P = NULL; + } + + // Assignment operator. + const Conic_arc_2& operator= (const Conic_arc_2& arc) + { + if (this == &arc) + return (*this); + + // Free any existing data. + if ((_info & IS_HYPERBOLA) != 0) + delete _data.hyper_P; + else if ((_info & IS_CIRCLE) != 0) + delete _data.circ_P; + _data.hyper_P = NULL; + + // Copy the arc's attributes. + _conic = arc._conic; + _conic_id = arc._conic_id; + _source = arc._source; + _target = arc._target; + _info = arc._info; + + // Duplicate the data for hyperbolic or circular arcs. + if ((arc._info & IS_HYPERBOLA) != 0) + { + _data.hyper_P = new Hyperbolic_arc_data (*(arc._data.hyper_P)); + } + else if ((arc._info & IS_CIRCLE) != 0) + { + _data.circ_P = new Circular_arc_data (*(arc._data.circ_P)); + } + +#ifdef CGAL_CONIC_ARC_USE_BOUNDING_BOX + _bbox = arc._bbox; // Copy the bounding box. +#endif + + return (*this); + } + + // Get the arc's base conic. + const Conic_2& conic () const + { + return (_conic); + } + + // Get the arc's source. + const Point_2& source () const + { + return (_source); + } + + // Get the arc's target. + const Point_2& target () const + { + return (_target); + } + + // Check whether the two arcs have the same base conic. + bool has_same_base_conic (const Conic_arc_2& arc) const + { + if (_conic_id == arc._conic_id) + return (true); + else + return (_conic == arc._conic); + } + + // Check whether the arc is a full conic (i.e. a non-degenerate ellipse). + bool is_full_conic () const + { + return ((_info & FULL_CONIC) != 0); + } + + // Check whether the curve is a sgement. + bool is_segment () const + { + return ((_info & DEGREE_MASK) == 1); + } + + // Check whether the curve is a vertical segment. + bool is_vertical_segment () const + { + // A vertical segment is contained in the degenerate conic: u*x + w = 0. + static const NT _zero = 0; + + return ((_info & DEGREE_MASK) == 1 && _conic.v() == _zero); + } + + // Check whether the curve is a horizontal segment. + bool is_horizontal_segment () const + { + // A vertical segment is contained in the degenerate conic: v*y + w = 0. + static const NT _zero = 0; + + return ((_info & DEGREE_MASK) == 1 && _conic.u() == _zero); + } + + // Get a bounding box for the conic arc. + Bbox_2 bounding_box () const + { + // Use the source and target to initialize the exterme points. + bool source_left = + CGAL::to_double(_source.x()) < CGAL::to_double(_target.x()); + double x_min = source_left ? + CGAL::to_double(_source.x()) : CGAL::to_double(_target.x()); + double x_max = source_left ? + CGAL::to_double(_target.x()) : CGAL::to_double(_source.x()); + bool source_down = + CGAL::to_double(_source.y()) < CGAL::to_double(_target.y()); + double y_min = source_down ? + CGAL::to_double(_source.y()) : CGAL::to_double(_target.y()); + double y_max = source_down ? + CGAL::to_double(_target.y()) : CGAL::to_double(_source.y()); + + // Go over the vertical tangency points and try to update the x-points. + Point_2 tps[2]; + int n_tps; + int i; + + n_tps = vertical_tangency_points (tps); + for (i = 0; i < n_tps; i++) + { + if (CGAL::to_double(tps[i].x()) < x_min) + x_min = CGAL::to_double(tps[i].x()); + if (CGAL::to_double(tps[i].x()) > x_max) + x_max = CGAL::to_double(tps[i].x()); + } + + // Go over the horizontal tangency points and try to update the y-points. + n_tps = horizontal_tangency_points (tps); + for (i = 0; i < n_tps; i++) + { + if (CGAL::to_double(tps[i].y()) < y_min) + y_min = CGAL::to_double(tps[i].y()); + if (CGAL::to_double(tps[i].y()) > y_max) + y_max = CGAL::to_double(tps[i].y()); + } + + // Return the resulting bounding box. + return (Bbox_2 (x_min, y_min, x_max, y_max)); + } + + // Check whether the given point is on the arc. + bool contains_point (const Point_2& p) const + { + // Check whether the conic contains the point (x,y). + if (p.is_generating_conic_id(_conic_id) || + (p.is_approximate() && _conic_has_approx_point_on_boundary (p)) || + _conic.has_on_boundary(p)) + { + // If the point is on the conic boundary, it is contained in the arc + // either if the arc is a full conic, or if it is between the two + // endpoints of the arc. + return (is_full_conic() || _is_between_endpoints(p)); + } + else + { + // If the point is not on the conic boundary, it cannot be on the arc. + return (false); + } + } + + // Calculate the vertical tangency points of the arc (ps should be allocated + // at the size of 2). + // The function return the number of vertical tangency points. + int vertical_tangency_points (Point_2* vpts) const + { + // No vertical tangency points for segments or for x-monotone curves: + if ((_info & DEGREE_MASK) < 2 || + (_info & X_MON_UNDEFINED) == X_MONOTONE) + return (0); + + // Calculate the vertical tangency points of the conic. + Point_2 ps[2]; + int n; + + n = _conic_vertical_tangency_points (ps); + + // Return only the points that are contained in the arc interior. + int m = 0; + + for (int i = 0; i < n; i++) + { + if (is_full_conic() || _is_strictly_between_endpoints(ps[i])) + { + vpts[m] = ps[i]; + m++; + } + } + + // Return the number of vertical tangency points found. + return (m); + } + + // Calculate the horizontal tangency points of the arc (ps should be + // allocated to the size of 2). + // The function return the number of vertical tangency points. + int horizontal_tangency_points (Point_2* hpts) const + { + // No horizontal tangency points for segments: + if ((_info & DEGREE_MASK) < 2) + return (0); + + // Calculate the horizontal tangency points of the conic. + Point_2 ps[2]; + int n; + + n = _conic_horizontal_tangency_points (ps); + + // Return only the points that are contained in the arc interior. + int m = 0; + + for (int i = 0; i < n; i++) + { + if (is_full_conic() || _is_strictly_between_endpoints(ps[i])) + { + hpts[m] = ps[i]; + m++; + } + } + + // Return the number of horizontal tangency points found. + return (m); + } + + // Check whether the arc is x-monotone. + bool is_x_monotone() const + { + // If the answer is pre-calculated (and stored in the _info field), just + // return it: + int is_x_mon = _info & X_MON_UNDEFINED; + + if (is_x_mon == 0) + return (false); + else if (is_x_mon == X_MONOTONE) + return (true); + + // Check the number of vertical tangency points. + Point_2 vpts[2]; + + return (vertical_tangency_points(vpts) == 0); + } + + // Find all points on the arc with a given x-coordinate: ps should be + // allocated to the size of 2. + // The function return the number of points found. + int get_points_at_x (const Point_2& p, + Point_2 *ps) const + { + // Make sure the conic is not a vertical segment. + CGAL_precondition(!is_vertical_segment()); + + // Get the y co-ordinates of the points on the conic. + NT ys[2]; + int n; + + n = _conic_get_y_coordinates (p.x(), ys); + + // Find all the points that are contained in the arc. + int m = 0; + + for (int i = 0; i < n; i++) + { + ps[m] = Point_2 (p.x(), ys[i], + p.is_approximate() ? Point_2::Ray_shooting_approx : + Point_2::Ray_shooting_exact, + _conic_id); + + if (is_full_conic() || _is_between_endpoints(ps[m])) + m++; + } + + // Return the number of points on the arc. + return (m); + } + + // Find all points on the arc with a given y-coordinate: ps should be + // allocated to the size of 2. + // The function return the number of points found. + int get_points_at_y (const Point_2& p, + Point_2 *ps) const + { + // Make sure the conic is not a horizontal segment. + CGAL_precondition(!is_horizontal_segment()); + + // Get the x co-ordinates of the points on the conic. + NT xs[2]; + int n; + + n = _conic_get_x_coordinates (p.y(), xs); + + // Find all the points that are contained in the arc. + int m = 0; + + for (int i = 0; i < n; i++) + { + ps[m] = Point_2 (xs[i], p.y(), + p.is_approximate() ? Point_2::Ray_shooting_approx : + Point_2::Ray_shooting_exact, + _conic_id); + + if (is_full_conic() || _is_between_endpoints(ps[m])) + m++; + } + + // Return the number of points on the arc. + return (m); + } + + // Return a flipped conic arc. + Conic_arc_2 flip () const + { + // Create an identical conic with an opposite orientation. + Conic_2 opp_conic (-_conic.r(), -_conic.s(), -_conic.t(), + -_conic.u(), -_conic.v(), -_conic.w()); + + + // Create the reflected curve (exchange the source and the target): + Conic_arc_2 opp_arc; + + opp_arc._conic = opp_conic; + opp_arc._conic_id = _conic_id; + opp_arc._source = _target; + opp_arc._target = _source; + opp_arc._info = _info; // These properties do not change. + + if ((_info & IS_HYPERBOLA) != 0) + { + opp_arc._data.hyper_P = new Hyperbolic_arc_data (*_data.hyper_P); + } + else if ((_info & IS_CIRCLE) != 0) + { + opp_arc._data.circ_P = new Circular_arc_data (*_data.circ_P); + } + else + { + opp_arc._data.hyper_P = NULL; + } + +#ifdef CGAL_CONIC_ARC_USE_BOUNDING_BOX + opp_arc._bbox = _bbox; // The bounding box is the same. +#endif + + return (opp_arc); + } + + // Reflect the curve in the y axis. + Conic_arc_2 reflect_in_y () const + { + // Reflect the base conic in y: + Conic_2 ref_conic ( _conic.r(), + _conic.s(), + - _conic.t(), + - _conic.u(), + _conic.v(), + _conic.w()); + + // Create the reflected curve: + Conic_arc_2 ref_arc; + + ref_arc._conic = ref_conic; + ref_arc._conic_id = _conic_id ^ REFLECT_IN_Y; + ref_arc._source = _source.reflect_in_y(); + ref_arc._target = _target.reflect_in_y(); + ref_arc._info = _info; // These properties do not change. + + if ((_info & IS_HYPERBOLA) != 0) + { + ref_arc._data.hyper_P = new Hyperbolic_arc_data (*_data.hyper_P); + ref_arc._data.hyper_P->a = - _data.hyper_P->a; + } + else if ((_info & IS_CIRCLE) != 0) + { + ref_arc._data.circ_P = new Circular_arc_data (*_data.circ_P); + ref_arc._data.circ_P->x0 = - _data.circ_P->x0; + } + else + { + ref_arc._data.hyper_P = NULL; + } + +#ifdef CGAL_CONIC_ARC_USE_BOUNDING_BOX + ref_arc._bbox = ref_arc.bounding_box(); // Compute the bounding box. +#endif + + return (ref_arc); + } + + // Reflect the curve in the x and y axes. + Conic_arc_2 reflect_in_x_and_y () const + { + // Reflect the base conic in x and y: + Conic_2 ref_conic ( _conic.r(), + _conic.s(), + _conic.t(), + - _conic.u(), + - _conic.v(), + _conic.w()); + + // Create the reflected curve: + Conic_arc_2 ref_arc; + + ref_arc._conic = ref_conic; + ref_arc._conic_id = _conic_id ^ (REFLECT_IN_X | REFLECT_IN_Y); + ref_arc._source = _source.reflect_in_x_and_y(); + ref_arc._target = _target.reflect_in_x_and_y(); + ref_arc._info = _info; // These properties do not change. + + if ((_info & IS_HYPERBOLA) != 0) + { + ref_arc._data.hyper_P = new Hyperbolic_arc_data (*_data.hyper_P); + ref_arc._data.hyper_P->a = - _data.hyper_P->a; + ref_arc._data.hyper_P->b = - _data.hyper_P->b; + } + else if ((_info & IS_CIRCLE) != 0) + { + ref_arc._data.circ_P = new Circular_arc_data (*_data.circ_P); + ref_arc._data.circ_P->x0 = - _data.circ_P->x0; + ref_arc._data.circ_P->y0 = - _data.circ_P->y0; + } + else + { + ref_arc._data.hyper_P = NULL; + } + + +#ifdef CGAL_CONIC_ARC_USE_BOUNDING_BOX + ref_arc._bbox = ref_arc.bounding_box(); // Compute the bounding box. +#endif + + return (ref_arc); + } + + // Get the i'th order derivative by x of the conic at the point p=(x,y). + // Note that i should be either 1 (first order) or 2 (second order). + void derive_by_x_at (const Point_2& p, const int& i, + NT& slope_numer, NT& slope_denom) const + { + // Make sure i is either 1 or 2. + CGAL_precondition(i == 1 || i == 2); + + // Make sure p is contained in the arc. + CGAL_precondition(contains_point(p)); + + // The derivative by x of the conic + // C: {r*x^2 + s*y^2 + t*xy + u*x + v*y + w = 0} + // at the point p=(x,y) is given by: + // + // 2r*x + t*y + u alpha + // y' = - ---------------- = - ------- + // 2s*y + t*x + v beta + // + const NT _two = 2; + const NT sl_numer = _two*_conic.r()*p.x() + _conic.t()*p.y() + _conic.u(); + const NT sl_denom = _two*_conic.s()*p.y() + _conic.t()*p.x() + _conic.v(); + + if (i == 1) + { + if (sl_denom > 0) + { + slope_numer = -sl_numer; + slope_denom = sl_denom; + } + else + { + slope_numer = sl_numer; + slope_denom = -sl_denom; + } + + return; + } + + // The second derivative is given by: + // + // s*alpha^2 - t*alpha*beta + r*beta^2 + // y'' = -2 ------------------------------------- + // beta^3 + // + const NT sl2_numer = _conic.s() * sl_numer * sl_numer - + _conic.t() * sl_numer * sl_denom + + _conic.r() * sl_denom * sl_denom; + const NT sl2_denom = sl_denom * sl_denom * sl_denom; + + if (sl_denom > 0) // so sl2_denom > 0 as well ... + { + slope_numer = -_two *sl2_numer; + slope_denom = sl2_denom; + } + else + { + slope_numer = _two *sl2_numer; + slope_denom = -sl2_denom; + } + + return; + } + + // Get the i'th order derivative by y of the conic at the point p=(x,y). + // Note that i should be either 1 (first order) or 2 (second order). + void derive_by_y_at (const Point_2& p, const int& i, + NT& slope_numer, NT& slope_denom) const + { + // Make sure i is either 1 or 2. + CGAL_precondition(i == 1 || i == 2); + + // Make sure p is contained in the arc. + CGAL_precondition(contains_point(p)); + + // The derivative by y of the conic + // C: {r*x^2 + s*y^2 + t*xy + u*x + v*y + w = 0} + // at the point p=(x,y) is given by: + // + // 2s*y + t*x + v alpha + // x' = - ---------------- = ------- + // 2r*x + t*y + u beta + // + const NT _two = 2; + const NT sl_numer = _two*_conic.s()*p.y() + _conic.t()*p.x() + _conic.v(); + const NT sl_denom = _two*_conic.r()*p.x() + _conic.t()*p.y() + _conic.u(); + + if (i == 1) + { + if (sl_denom > 0) + { + slope_numer = -sl_numer; + slope_denom = sl_denom; + } + else + { + slope_numer = sl_numer; + slope_denom = -sl_denom; + } + + return; + } + + // The second derivative is given by: + // + // r*alpha^2 - t*alpha*beta + s*beta^2 + // x'' = -2 ------------------------------------- + // beta^3 + // + const NT sl2_numer = _conic.r() * sl_numer * sl_numer - + _conic.t() * sl_numer * sl_denom + + _conic.s() * sl_denom * sl_denom; + const NT sl2_denom = sl_denom * sl_denom * sl_denom; + + if (sl_denom > 0) // so sl2_denom > 0 as well ... + { + slope_numer = -_two *sl2_numer; + slope_denom = sl2_denom; + } + else + { + slope_numer = _two *sl2_numer; + slope_denom = -sl2_denom; + } + + return; + } + + // Calculate the intersection points between the arc and the given arc. + // ps must be allocated at the size of 4. + // The function returns the number of the actual intersection points. + int intersections_with (const Conic_arc_2& arc, + Point_2* ps +#ifdef CGAL_CONIC_ARC_USE_CACHING + ,std::list *inter_list_P = NULL +#endif + ) const + { +#ifdef CGAL_CONIC_ARC_USE_BOUNDING_BOX + // Perform quick rejections, if possible. + if (! do_overlap(_bbox, arc._bbox)) + return (false); +#endif + + // First make sure that (this->degree) is >= than (arc.degree). + if ((arc._info & DEGREE_MASK) == DEGREE_2 && + (_info & DEGREE_MASK) == DEGREE_1) + { + return (arc.intersections_with (*this, ps +#ifdef CGAL_CONIC_ARC_USE_CACHING + ,inter_list_P +#endif + )); + } + + static const NT _zero = 0; + + // The two conics must not be the same. + CGAL_precondition(_conic != arc._conic); + + // Deal with vertical segments. + if (arc.is_vertical_segment()) + { + if (is_vertical_segment()) + { + // Two vertical segments intersect only if they overlap. + return (0); + } + + // Find all points on our arc that have the same x co-ordinate as + // the other vertical segment. + int n_ys; + Point_2 xps[2]; + int j; + int n = 0; + + n_ys = get_points_at_x (arc._source, xps); + + for (j = 0; j < n_ys; j++) + { + // Store this point only if it is contained on the other arc. + if (arc.contains_point(xps[j])) + { + // Return an exact point: + ps[n] = Point_2 (xps[j].x(), xps[j].y(), + Point_2::Intersection_exact, + _conic_id, arc._conic_id); + n++; + } + } + + return (n); + } + else if (is_vertical_segment()) + { + // Find all points on the other arc that have the same x co-ordinate as + // our vertical segment. + int n_ys; + Point_2 xps[2]; + int j; + int n = 0; + + n_ys = arc.get_points_at_x (_source, xps); + + for (j = 0; j < n_ys; j++) + { + // Store this point only if it is contained on the other arc. + if (contains_point(xps[j])) + { + // Return an exact point: + ps[n] = Point_2 (xps[j].x(), xps[j].y(), + Point_2::Intersection_exact, + _conic_id, arc._conic_id); + n++; + } + } + + return (n); + } + + // Find all intersection points between the two base conic curves. + Point_2 ipts[4]; // The intersection points. + int n_points = 0; // Their number. + bool calc_points = true; + // REF_TRICK + int reflect_this = (_conic_id & (REFLECT_IN_X | REFLECT_IN_Y)); + int k; + + CGAL_assertion (reflect_this == + (arc._conic_id & (REFLECT_IN_X | REFLECT_IN_Y))); + +#ifdef CGAL_CONIC_ARC_USE_CACHING + Intersections inter; + + if (inter_list_P != NULL && + (_info & DEGREE_MASK) != DEGREE_1) + { + // REF-TRICK + int id1 = _conic_id / REFLECTION_FACTOR; + int id2 = arc._conic_id / REFLECTION_FACTOR; + //int id1 = _conic_id; + //int id2 = arc._conic_id; + + inter.id1 = id1 < id2 ? id1 : id2; + inter.id2 = id1 > id2 ? id1 : id2; + + typename std::list::iterator iter; + for (iter = inter_list_P->begin(); iter != inter_list_P->end(); iter++) + { + if ((*iter).id1 == inter.id1 && (*iter).id2 == inter.id2) + { + n_points = (*iter).n_points; + + for (k = 0; k < n_points; k++) + { + // REF-TRICK + if (reflect_this == (REFLECT_IN_X | REFLECT_IN_Y)) + ipts[k] = (*iter).ps[k].reflect_in_x_and_y(); + else if (reflect_this == (REFLECT_IN_Y)) + ipts[k] = (*iter).ps[k].reflect_in_y(); + else + ipts[k] = (*iter).ps[k]; + //ipts[k] = (*iter).ps[k]; + } + calc_points = false; + } + } + } +#endif // (of ifdef CGAL_CONIC_ARC_USE_CACHING) + + if (calc_points) + { + // Find all potential x co-ordinates and y co-ordinates of the + // intersection points. + NT xs[4]; + int x_mults[4]; + int n_xs; // Total number of x co-ordinates. + int n_approx_xs; // Number of approximate x co-ordinates. + NT ys[4]; + int y_mults[4]; + int n_ys; // Total number of y co-ordinates. + int n_approx_ys; // Number of approximate y co-ordinates. + int x_deg; // The generating polynomial for the x values. + std::vector x_coeffs; + int y_deg; // The generating polynomial for the y values. + std::vector y_coeffs; + + if (_conic.s() == _zero && arc._conic.s() != _zero) + { + n_xs = arc._x_coordinates_of_intersections_with (*this, + xs, x_mults, + n_approx_xs, + x_deg, + x_coeffs); + } + else + { + n_xs = _x_coordinates_of_intersections_with (arc, + xs, x_mults, + n_approx_xs, + x_deg, + x_coeffs); + } + + if (_conic.r() == _zero && arc._conic.r() != _zero) + { + n_ys = arc._y_coordinates_of_intersections_with (*this, + ys, y_mults, + n_approx_ys, + y_deg, y_coeffs); + } + else + { + n_ys = _y_coordinates_of_intersections_with (arc, + ys, y_mults, + n_approx_ys, + y_deg, y_coeffs); + } + + // Perform the pairing process od the x and y coordinates. + n_points = _pair_intersection_points (arc, + n_xs, + xs, x_mults, + n_approx_xs, + n_ys, + ys, y_mults, + n_approx_ys, + ipts); + + // Attach the generating polynomials for the intersection points. + for (k = 0; k < n_points; k++) + { + ipts[k].attach_polynomials (x_deg, x_coeffs, + y_deg, y_coeffs); + } + +#ifdef CGAL_CONIC_ARC_USE_CACHING + if (inter_list_P != NULL && + (_info & DEGREE_MASK) != DEGREE_1) + { + inter.n_points = n_points; + + for (k = 0; k < n_points; k++) + { + // REF-TRICK + if (reflect_this == (REFLECT_IN_X | REFLECT_IN_Y)) + inter.ps[k] = ipts[k].reflect_in_x_and_y(); + else if (reflect_this == (REFLECT_IN_Y)) + inter.ps[k] = ipts[k].reflect_in_y(); + else + inter.ps[k] = ipts[k]; + //inter.ps[k] = ipts[k]; + } + + inter_list_P->push_front(inter); + } + +#endif // (of ifdef CGAL_CONIC_ARC_USE_CACHING) + } + + // Go over all intersection points between the two base conics and return + // only those located on both arcs. + int n = 0; + int i; + + for (i = 0; i < n_points; i++) + { + // Check for an exact point. + if (contains_point(ipts[i]) && + arc.contains_point(ipts[i])) + { + ps[n] = ipts[i]; + n++; + } + } + + return (n); + } + + // Check whether the two arcs overlap. + // The function computes the number of overlapping arcs (2 at most), and + // return their number (0 means there is no overlap). + int overlaps (const Conic_arc_2& arc, + Conic_arc_2* ovlp_arcs) const + { + // Two arcs can overlap only if their base conics are identical. + if (_conic != arc._conic) + return (0); + + // If the two arcs are completely equal, return one of them as the + // overlapping arc. + int orient1 = _conic.orientation(); + int orient2 = arc._conic.orientation(); + bool same_or = (orient1 == orient2); + bool identical = false; + + if (orient1 == 0) + { + // That mean both arcs are really segments, so they are identical + // if their endpoints are the same. + if ((_source.equals(arc._source) && _target.equals(arc._target)) || + (_source.equals(arc._target) && _target.equals(arc._source))) + identical = true; + } + else + { + // If those are really curves of degree 2, than the points curves + // are identical only if their source and target are the same and the + // orientation is the same, or vice-versa if the orientation is opposite. + if ((same_or && + _source.equals(arc._source) && _target.equals(arc._target)) || + (!same_or && + _source.equals(arc._target) && _target.equals(arc._source))) + identical = true; + } + + if (identical) + { + ovlp_arcs[0] = arc; + return (1); + } + + // In case one of the arcs is a full conic, return the whole other conic. + if (arc.is_full_conic()) + { + ovlp_arcs[0] = *this; + return (1); + } + else if (is_full_conic()) + { + ovlp_arcs[0] = arc; + return (1); + } + + // In case the other arc has an opposite orientation, switch its source + // and target (notice that in case of segments, when the orientation is 0, + // we make sure the two segments have the same direction). + const Point_2 *arc_sourceP; + const Point_2 *arc_targetP; + + if (orient1 == 0) + orient1 = (_source.compare_lex_xy(_target) + == LARGER) ? 1 : -1; + if (orient2 == 0) + orient2 = (arc._source.compare_lex_xy(arc._target) + == LARGER) ? 1 : -1; + + // Check the overlap cases: + if (orient1 == orient2) + { + arc_sourceP = &(arc._source); + arc_targetP = &(arc._target); + } + else + { + arc_sourceP = &(arc._target); + arc_targetP = &(arc._source); + } + + if (_is_strictly_between_endpoints(*arc_sourceP)) + { + if (_is_strictly_between_endpoints(*arc_targetP)) + { + // Check the next special case (when there are 2 overlapping arcs): + if (arc._is_strictly_between_endpoints(_source) && + arc._is_strictly_between_endpoints(_target)) + { + ovlp_arcs[0] = Conic_arc_2(*this,_source, *arc_targetP, false); + ovlp_arcs[1] = Conic_arc_2(*this, *arc_sourceP, _target, false); + //ovlp_arcs[0] = Conic_arc_2(_conic, _source, *arc_targetP); + //ovlp_arcs[1] = Conic_arc_2(_conic, *arc_sourceP, _target); + return (2); + } + + // Case 1 - *this: +-----------> + // arc: +=====> + ovlp_arcs[0] = Conic_arc_2(*this, *arc_sourceP,*arc_targetP, false); + //ovlp_arcs[0] = Conic_arc_2(_conic, *arc_sourceP,*arc_targetP); + return (1); + } + else + { + // Case 2 - *this: +-----------> + // arc: +=====> + ovlp_arcs[0] = Conic_arc_2(*this, *arc_sourceP, _target, false); + //ovlp_arcs[0] = Conic_arc_2(_conic, *arc_sourceP, _target); + return (1); + } + } + else if (_is_strictly_between_endpoints(*arc_targetP)) + { + // Case 3 - *this: +-----------> + // arc: +=====> + ovlp_arcs[0] = Conic_arc_2(*this, _source, *arc_targetP, false); + //ovlp_arcs[0] = Conic_arc_2(_conic, _source, *arc_targetP); + return (1); + } + else if (arc._is_between_endpoints(_source) && + arc._is_between_endpoints(_target) && + (arc._is_strictly_between_endpoints(_source) || + arc._is_strictly_between_endpoints(_target))) + { + // Case 4 - *this: +-----------> + // arc: +================> + ovlp_arcs[0] = *this; + return (1); + } + + // If we reached here, there are no overlaps: + return (0); + } + + // Check whether the arc is facing up (LARGER is then returned), + // facing down (SMALLER is then returned). + // At any other case the function returns EQUAL. + Comparison_result facing () const + { + if ((_info & FACING_MASK) == 0) + return (EQUAL); + else if ((_info & FACING_UP) != 0) + return (LARGER); + else + return (SMALLER); + } + + private: + + // Set the properties of a conic arc (for the usage of the constructors). + // The source and the target are assumed be on the conic boundary. + void _set (const Conic_2& conic, + const Point_2& source, const Point_2& target, + Circular_arc_data *circ_data_P = NULL) + { + // Set the data members. + _conic = conic; + _conic_id = 0; + _source = source; + _target = target; + _info = X_MON_UNDEFINED; + + // Find the degree and make sure the conic is not invalid. + static const NT _zero = 0; + int deg; + + if (_conic.r() != _zero || _conic.s() != _zero || _conic.t() != _zero) + { + // In case one of the coefficients of x^2,y^2 or xy is not zero, the + // degree is 2. + deg = 2; + } + else if (_conic.u() != _zero || _conic.v() != _zero) + { + // In case of a line - the degree is 1. + deg = 1; + } + else + { + // Empty conic! + deg = 0; + } + + CGAL_precondition(deg > 0); + + _info = _info | deg; + + // In case the base conic is a hyperbola, build the hyperbolic data. + if (deg == 2 && _conic.is_hyperbola()) + { + _info = _info | IS_HYPERBOLA; + _build_hyperbolic_arc_data (); + } + // In case the base conic is a circle, set the circular data. + else if (deg == 2 && circ_data_P != NULL) + { + _info = _info | IS_CIRCLE; + _data.circ_P = circ_data_P; + } + else + { + _data.hyper_P = NULL; + } + + // In case of a non-degenerate parabola or a hyperbola, make sure + // the arc is not infinite. + if (deg == 2 && ! _conic.is_ellipse()) + { + CGAL_precondition_code( + const NT _two = 2; + const Point_2 p_mid ((source.x() + target.x()) / _two, + (source.y() + target.y()) / _two); + Point_2 ps[2]; + + bool finite_at_x = (this->get_points_at_x(p_mid, ps) > 0); + bool finite_at_y = (this->get_points_at_y(p_mid, ps) > 0); + ); + CGAL_precondition(finite_at_x && finite_at_y); + } + + // If we reached here, the conic arc is legal: Get a new id for the conic. + _conic_id = _get_new_conic_id(); + + _source = Point_2 (_source.x(), _source.y(), + Point_2::User_defined, + _conic_id); + + _target = Point_2 (_target.x(), _target.y(), + Point_2::User_defined, + _conic_id); + + // Check whether the conic is x-monotone. + if (is_x_monotone()) + { + _info = (_info & ~X_MON_UNDEFINED) | X_MONOTONE; + + // In case the conic is od degree 2, determine where is it facing. + if ((_info & DEGREE_MASK) == DEGREE_2) + _set_facing(); + } + else + { + _info = (_info & ~X_MON_UNDEFINED); + } + +#ifdef CGAL_CONIC_ARC_USE_BOUNDING_BOX + _bbox = bounding_box(); // Compute the bounding box. +#endif + + return; + } + + // Set an arc which is basically a full conic (an ellipse). + void _set_full (const Conic_2& conic, + Circular_arc_data *circ_data_P = NULL) + { + // Set the data members. + _conic = conic; + _conic_id = 0; + _info = 0; + + static const NT _zero = 0; + + // Make sure the conic is a non-degenerate ellipse. + CGAL_precondition(_conic.is_ellipse()); + + // Set the information: a full conic, which is obvoiusly not x-monotone. + _info = DEGREE_2 | FULL_CONIC; + + // Check if the conic is really a circle. + if (circ_data_P != NULL) + { + _info = _info | IS_CIRCLE; + _data.circ_P = circ_data_P; + } + else + { + _data.hyper_P = NULL; + } + + // Assign one of the vertical tangency points as both the source and + // the target of the conic arc. + Point_2 vpts[2]; + int n_vpts; + + if (circ_data_P != NULL) + { + vpts[0] = Point_2 (circ_data_P->x0 + circ_data_P->r, + circ_data_P->y0); + } + else + { + n_vpts = _conic_vertical_tangency_points (vpts); + + CGAL_assertion(n_vpts > 0); + CGAL_assertion(_conic.has_on_boundary(vpts[0])); + } + + // If we reached here, the conic arc is legal: Get a new id for the conic. + _conic_id = _get_new_conic_id(); + + _source = Point_2 (vpts[0].x(), vpts[0].y(), + Point_2::User_defined, + _conic_id); + + _target = Point_2 (vpts[0].x(), vpts[0].y(), + Point_2::User_defined, + _conic_id); + +#ifdef CGAL_CONIC_ARC_USE_BOUNDING_BOX + _bbox = bounding_box(); // Compute the bounding box. +#endif + + return; + } + + // Build the data for hyperbolic arc, contaning the characterization of the + // hyperbolic branch the arc is placed on. + void _build_hyperbolic_arc_data () + { + // Let phi be the rotation angle of the conic from its canonic form. + // We can write: + // + // t + // sin(2*phi) = ----------------------- + // sqrt((r - s)^2 + t^2) + // + // r - s + // cos(2*phi) = ----------------------- + // sqrt((r - s)^2 + t^2) + // + const NT r = _conic.orientation() * _conic.r(); + const NT s = _conic.orientation() * _conic.s(); + const NT t = _conic.orientation() * _conic.t(); + const NT cos_2phi = (r - s) / CGAL::sqrt((r-s)*(r-s) + t*t); + static const NT _zero(0); + static const NT _half(0.5); + static const NT _one(1); + static const NT _two(2); + NT sin_phi; + NT cos_phi; + + // Calculate sin(phi) and cos(phi) according to the half-edge formulae: + // + // sin(phi)^2 = 0.5 * (1 - cos(2*phi)) + // cos(phi)^2 = 0.5 * (1 + cos(2*phi)) + if (t == _zero) + { + // sin(2*phi) == 0, so phi = 0 or phi = PI/2 + if (cos_2phi > _zero) + { + // phi = 0. + sin_phi = _zero; + cos_phi = _one; + } + else + { + // phi = PI. + sin_phi = _zero; + cos_phi = -_one; + } + } + else if (t > _zero) + { + // sin(2*phi) > 0 so 0 < phi < PI/2. + sin_phi = CGAL::sqrt((_one + cos_2phi) / _two); + cos_phi = CGAL::sqrt((_one - cos_2phi) / _two); + } + else + { + // sin(2*phi) < 0 so PI/2 < phi < PI. + sin_phi = CGAL::sqrt((_one + cos_2phi) / _two); + cos_phi = -CGAL::sqrt((_one - cos_2phi) / _two); + } + + // Calculate the centre (x0, y0) of the conic, given by the formulae: + // + // t*v - 2*s*u t*u - 2*r*v + // x0 = ------------- , y0 = ------------- + // 4*r*s - t^2 4*r*s - t^2 + // + // The denominator (4*r*s - t^2) must be negative for hyperbolas. + const NT u = _conic.orientation() * _conic.u(); + const NT v = _conic.orientation() * _conic.v(); + const NT det = NT(4)*r*s - t*t; + NT x0, y0; + + CGAL_assertion (det < _zero); + + x0 = (t*v - _two*s*u) / det; + y0 = (t*u - _two*r*v) / det; + + // The axis separating the two branches of the hyperbola is now given by: + // + // cos(phi)*x + sin(phi)*y - (cos(phi)*x0 + sin(phi)*y0) = 0 + // + _data.hyper_P = new Hyperbolic_arc_data; + + _data.hyper_P->a = cos_phi; + _data.hyper_P->b = sin_phi; + _data.hyper_P->c = - (cos_phi*x0 + sin_phi*y0); + + // Make sure that the two endpoints are located on the same branch + // of the hyperbola. + _data.hyper_P->side = _hyperbolic_arc_side(_source); + + CGAL_assertion (_data.hyper_P->side = _hyperbolic_arc_side(_target)); + + return; + } + + // Find on which branch of the hyperbola is the given point located. + // The point is assumed to be on the hyperbola. + int _hyperbolic_arc_side (const Point_2& p) const + { + if (_data.hyper_P == NULL) + return (0); + + NT val; + + val = _data.hyper_P->a*p.x() + _data.hyper_P->b*p.y() + _data.hyper_P->c; + return ((val > 0) ? 1 : -1); + } + + // Check whether the given point is between the source and the target. + // The point is assumed to be on the conic's boundary. + bool _is_between_endpoints (const Point_2& p) const + { + if (p.equals(_source) || p.equals(_target)) + return (true); + else + return (_is_strictly_between_endpoints(p)); + } + + // Check whether the given point is strictly between the source and the + // target (but not any of them). + // The point is assumed to be on the conic's boundary. + bool _is_strictly_between_endpoints (const Point_2& p) const + { + // In case this is a full conic, any point on its boundary is between + // its end points. + if (is_full_conic()) + return (true); + + // In case of a hyperbolic arc, make sure the point is located on the + // same branch as the arc. + if ((_info & IS_HYPERBOLA) != 0) + { + if (_hyperbolic_arc_side(p) != _data.hyper_P->side) + return (false); + } + + // Act according to the conic degree. + if ((_info & DEGREE_MASK) == DEGREE_1) + { + if (is_vertical_segment()) + { + // In case of a vertical segment - just check whether the y co-ordinate + // of p is between those of the source's and of the target's. + Comparison_result r1 = compare_y (p, _source); + Comparison_result r2 = compare_y (p, _target); + + return ((r1 == SMALLER && r2 == LARGER) || + (r1 == LARGER && r2 == SMALLER)); + } + else + { + // Otherwise, since the segment is x-monotone, just check whether the + // x co-ordinate of p is between those of the source's and of the + // target's. + Comparison_result r1 = compare_x (p, _source); + Comparison_result r2 = compare_x (p, _target); + + return ((r1 == SMALLER && r2 == LARGER) || + (r1 == LARGER && r2 == SMALLER)); + } + } + else + { + // In case of a conic of degree 2, make a decision based on the conic's + // orientation and whether (source,p,target) is a right or a left turn. + static Kernel ker; + typename Kernel::Orientation_2 orient_f = ker.orientation_2_object(); + + if (_conic.orientation() == 1) + return (orient_f(_source, p, _target) == LEFTTURN); + else + return (orient_f(_source, p, _target) == RIGHTTURN); + } + } + + // Find the y-coordinates of the conic at a given x-coordinate. + int _conic_get_y_coordinates (const NT& x, + NT *ys) const + { + int mults[2]; + + // Solve the quadratic equation for a given x and find the y values: + // s*y^2 + (t*x + v)*y + (r*x^2 + u*x + w) = 0 + return (solve_quadratic_eq (_conic.s(), + x*_conic.t() + _conic.v(), + x*(x*_conic.r() + _conic.u()) + _conic.w(), + ys, mults)); + } + + // Find the x-coordinates of the conic at a given y-coordinate. + int _conic_get_x_coordinates (const NT& y, + NT *xs) const + { + int mults[2]; + + // Solve the quadratic equation for a given y and find the x values: + // r*x^2 + (t*y + u)*x + (s*y^2 + v*y + w) = 0 + return (solve_quadratic_eq (_conic.r(), + y*_conic.t() + _conic.u(), + y*(y*_conic.s() + _conic.v()) + _conic.w(), + xs, mults)); + } + + // Find the vertical tangency points of the conic. + int _conic_vertical_tangency_points (Point_2* ps) const + { + // In case the base conic is of degree 1 (and not 2), the arc has no + // vertical tangency points. + static const NT _zero = 0; + static const NT _two = 2; + static const NT _four = 4; + + if ((_info & DEGREE_MASK) == DEGREE_1 || _conic.s() == _zero) + return (0); + + // Special treatment for circles, where the vertical tangency points + // are simply (x0-r,y0) and (x0+r,y0). + if ((_info & IS_CIRCLE) != 0) + { + ps[0] = Point_2 (_data.circ_P->x0 - _data.circ_P->r, + _data.circ_P->y0, + Point_2::Tangency, + _conic_id); + ps[1] = Point_2 (_data.circ_P->x0 + _data.circ_P->r, + _data.circ_P->y0, + Point_2::Tangency, + _conic_id); + + return (2); + } + + // We are interested in the x co-ordinates were the quadratic equation: + // s*y^2 + (t*x + v)*y + (r*x^2 + u*x + w) = 0 + // has a single solution (obviously if s = 0, there are no such points). + // We therefore demand that the discriminant of this equation is zero: + // (t*x + v)^2 - 4*s*(r*x^2 + u*x + w) = 0 + const NT r = _conic.r(); + const NT s = _conic.s(); + const NT t = _conic.t(); + const NT u = _conic.u(); + const NT v = _conic.v(); + const NT w = _conic.w(); + int n; + + NT xs[2]; + int n_xs; + int x_mults[2]; + NT ys[2]; + int n_ys; + int y_mults[2]; + + n_xs = solve_quadratic_eq (t*t - _four*r*s, + _two*t*v - _four*s*u, + v*v - _four*s*w, + xs, x_mults); + + // Store the generating polynomial for the x-coordinates. + int x_deg = 2; + std::vector x_coeffs(3); + int y_deg; + std::vector y_coeffs; + + x_coeffs[2] = t*t - _four*r*s; + x_coeffs[1] = _two*t*v - _four*s*u; + x_coeffs[0] = v*v - _four*s*w; + + // Find the y-coordinates of the vertical tangency points. + if (t == _zero) + { + n_ys = 1; + ys[0] = -v / (_two*s); + y_mults[0] = 2; + + y_deg = 0; + } + else + { + n_ys = solve_quadratic_eq (_four*r*s*s - s*t*t, + _four*r*s*v - _two*s*t*u, + r*v*v - t*u*v + t*t*w, + ys, y_mults); + + // Store the generating polynomial for the y-coordinates. + y_deg = 2; + y_coeffs.resize(3); + y_coeffs[2] = _four*r*s*s - s*t*t; + y_coeffs[1] = _four*r*s*v - _two*s*t*u; + y_coeffs[0] = r*v*v - t*u*v + t*t*w; + } + + // Pair the x and y coordinates and obtain the vertical tangency points. + n = 0; + for (int i = 0; i < n_xs; i++) + { + if (n_ys == 1) + { + ps[n] = Point_2 (xs[i], ys[0], + Point_2::Tangency, + _conic_id); + n++; + } + else + { + for (int j = 0; j < n_ys; j++) + { + if (ys[j] == -(t*xs[i] + v) / (_two*s)) + { + ps[n] = Point_2 (xs[i], ys[j], + Point_2::Tangency, + _conic_id); + n++; + break; + } + } + } + } + + // Attach the generating polynomials. + for (int i = 0; i < n; i++) + { + ps[i].attach_polynomials (x_deg, x_coeffs, + y_deg, y_coeffs); + } + + return (n); + } + + // Find the horizontal tangency points of the conic. + int _conic_horizontal_tangency_points (Point_2* ps) const + { + // In case the base conic is of degree 1 (and not 2), the arc has no + // vertical tangency points. + static const NT _zero = 0; + static const NT _two = 2; + static const NT _four = 4; + + if ((_info & DEGREE_MASK) == DEGREE_1 || _conic.r() == _zero) + return (0); + + // Special treatment for circles, where the horizontal tangency points + // are simply (x0,y0-r) and (x0,y0+r). + if ((_info & IS_CIRCLE) != 0) + { + ps[0] = Point_2 (_data.circ_P->x0, + _data.circ_P->y0 - _data.circ_P->r, + Point_2::Tangency, + _conic_id); + ps[1] = Point_2 (_data.circ_P->x0, + _data.circ_P->y0 + _data.circ_P->r, + Point_2::Tangency, + _conic_id); + + return (2); + } + + // We are interested in the y co-ordinates were the quadratic equation: + // r*x^2 + (t*y + u)*x + (s*y^2 + v*y + w) = 0 + // has a single solution (obviously if r = 0, there are no such points). + // We therefore demand that the discriminant of this equation is zero: + // (t*y + u)^2 - 4*r*(s*y^2 + v*y + w) = 0 + const NT r = _conic.r(); + const NT s = _conic.s(); + const NT t = _conic.t(); + const NT u = _conic.u(); + const NT v = _conic.v(); + const NT w = _conic.w(); + int n; + NT ys[2]; + int y_mults[2]; + NT x; + + n = solve_quadratic_eq (t*t - _four*r*s, + _two*t*u - _four*r*v, + u*u - _four*r*w, + ys, y_mults); + + for (int i = 0; i < n; i++) + { + // Having computed y, x is the simgle solution to the quadratic equation + // above, and since its discriminant is 0, x is simply given by: + x = -(t*ys[i] + u) / (_two*r); + + ps[i] = Point_2 (x, ys[i], + Point_2::Tangency, + _conic_id); + } + + return (n); + } + + // Check whether the base conic contains the given approximate point on its + // boundary. + bool _conic_has_approx_point_on_boundary (const Point_2& p) const + { + APNT r = TO_APNT(_conic.r()); + APNT s = TO_APNT(_conic.s()); + APNT t = TO_APNT(_conic.t()); + APNT u = TO_APNT(_conic.u()); + APNT v = TO_APNT(_conic.v()); + APNT w = TO_APNT(_conic.w()); + APNT x = TO_APNT(p.x()); + APNT y = TO_APNT(p.y()); + + APNT value = (r*x + u)*x + (s*y + t*x + v)*y + w; + + // RWRW: Correct this !!!! + return (APNT_ABS(value) < 0.001); + //return (eps_compare(value, 0) == EQUAL); + } + + // Set the facing information for the arc. + void _set_facing () + { + // Check whether the arc (which is x-monotone of degree 2) lies above or + // below the segement that contects its two end-points (x1,y1) and (x2,y2). + // To do that, we find the y co-ordinate of a point on the arc whose x + // co-ordinate is (x1+x2)/2 and compare it to (y1+y2)/2. + static const NT _two = 2; + const NT x_mid = (_source.x() + _target.x()) / _two; + const NT y_mid = (_source.y() + _target.y()) / _two; + Point_2 p_mid (x_mid, y_mid); + Point_2 ps[2]; + int n_ps; + + n_ps = get_points_at_x (p_mid, ps); + CGAL_assertion (n_ps == 1); + + Comparison_result res = CGAL::compare (ps[0].y(), y_mid); + + if (res == LARGER) + { + // The arc is above the connecting segment, so it is facing upwards. + _info = _info | FACING_UP; + } + else if (res == SMALLER) + { + // The arc is below the connecting segment, so it is facing downwards. + _info = _info | FACING_DOWN; + } + + CGAL_assertion(res != EQUAL); + return; + } + + // Calculate all x co-ordinates of intersection points between the two + // base curves of (*this) and the given arc. + int _x_coordinates_of_intersections_with (const Conic_arc_2& arc, + NT* xs, int* x_mults, + int& n_approx, + int& x_deg, + std::vector& x_coeffs) const + { + static const NT _zero = 0; + int n_roots; // The number of distinct x values. + + if ((_info & DEGREE_MASK) == DEGREE_1 && + (arc._info & DEGREE_MASK) == DEGREE_1) + { + // The two conics are: u*x + v*y + w = 0 + // and: u'*x + v'*y + w' = 0 + // There's a single solution for x, which is: + const NT denom = _conic.v()*arc._conic.w() - _conic.w()*arc._conic.v(); + const NT numer = _conic.u()*arc._conic.v() - _conic.v()*arc._conic.u(); + + if (numer == _zero) + { + n_roots = 0; + } + else + { + xs[0] = denom / numer; + x_mults[0] = 1; + n_roots = 1; + n_approx = 0; + } + + // Do not set any coefficients' information in this case. + x_deg = 0; + } + else if ((arc._info & DEGREE_MASK) == DEGREE_1) + { + // The two conics are: r*x^2 + s*y^2 + t*xy + u*x + v*y + w = 0 + // and: a*x + b*y + c = 0 + // There are therefore 2 possible x-values, the solutions for: + const NT a = arc._conic.u(); + const NT b = arc._conic.v(); + const NT c = arc._conic.w(); + const NT r = _conic.r(); + const NT s = _conic.s(); + const NT t = _conic.t(); + const NT u = _conic.u(); + const NT v = _conic.v(); + const NT w = _conic.w(); + + // Set the generating polynomial. + x_deg = 2; + x_coeffs.resize(3); + x_coeffs[2] = a*a*s + b*b*r - a*b*t; + x_coeffs[1] = 2*a*c*s + b*b*u - a*b*v - b*c*t; + x_coeffs[0] = c*c*s + b*b*w - b*c*v; + + // Find the roots. + n_roots = solve_quadratic_eq (x_coeffs[2], + x_coeffs[1], + x_coeffs[0], + xs, x_mults); + n_approx = 0; + } + else if ((_info & IS_CIRCLE) != 0 && (arc._info & IS_CIRCLE) != 0) + { + // Special treatment for two circles. + // The two curves are: r*x^2 + r*y^2 + u*x + v*y + w = 0 + // and: r'*x^2 + r'*y^2 + u'*x + v'*y + w' = 0 + // + // Thus, r'*C1-r*C2 is a line whose equation is: a*x + b*y + = 0, where: + const NT r = _conic.r(); + const NT s = _conic.s(); + const NT t = _conic.t(); + const NT u = _conic.u(); + const NT v = _conic.v(); + const NT w = _conic.w(); + const NT a = arc._conic.r()*u - r*arc._conic.u(); + const NT b = arc._conic.r()*v - r*arc._conic.v(); + const NT c = arc._conic.r()*w - r*arc._conic.w(); + + if (b == _zero) + { + // The line a*x + c = 0 connects both intersection points of the two + // circles, so the both have an x-coordinate of -c/a. + if (a == _zero) + { + n_roots = 0; + n_approx = 0; + } + else + { + xs[0] = -c / a; + x_mults[0] = 2; + n_roots = 1; + n_approx = 0; + } + // Do not set any coefficients' information in this case. + x_deg = 0; + } + else + { + // The intersection points of the two circles are the same as the + // intersection points of one of the circles with a*x + b*y + c = 0. + // Set the generating polynomial. + x_deg = 2; + x_coeffs.resize(3); + x_coeffs[2] = a*a*s + b*b*r - a*b*t; + x_coeffs[1] = 2*a*c*s + b*b*u - a*b*v - b*c*t; + x_coeffs[0] = c*c*s + b*b*w - b*c*v; + + // Find the roots. + n_roots = solve_quadratic_eq (x_coeffs[2], + x_coeffs[1], + x_coeffs[0], + xs, x_mults); + n_approx = 0; + } + } + else if (_conic.s() == _zero && _conic.t() == _zero && + arc._conic.s() == _zero && arc._conic.t() == _zero) + { + // Special treatment for canonic parabolas whose axes are parallel + // to the y axis. + // The two curves are: r*x^2 + u*x + v*y + w = 0 + // and: r'*x^2 + u'*x + v'*y + w' = 0 + // There are therefore 2 possible x-values, the solutions for: + + // Set the generating polynomial. + x_deg = 2; + x_coeffs.resize(3); + x_coeffs[2] = _conic.r()*arc._conic.v() - _conic.v()*arc._conic.r(); + x_coeffs[1] = _conic.u()*arc._conic.v() - _conic.v()*arc._conic.u(); + x_coeffs[0] = _conic.w()*arc._conic.v() - _conic.v()*arc._conic.w(); + + // Solve the equation. + n_roots = solve_quadratic_eq (x_coeffs[2], + x_coeffs[1], + x_coeffs[0], + xs, x_mults); + n_approx = 0; + } + else + { + // If the two conics are: r*x^2 + s*y^2 + t*xy + u*x + v*y + w = 0 + // and: r'*x^2 + s'*y^2 + t'*xy + u'*x + v'*y + w' = 0 + // Let us define: + // A = s'*r - s*r' if (s' != 0), otherwise A = r' + // B = s'*u - s*u' if (s' != 0), otherwise B = u' + // C = s'*w - s*w' if (s' != 0), otherwise C = w' + // D = s'*t - s*t' if (s' != 0), otherwise A = t' + // E = s'*v - s*v' if (s' != 0), otherwise A = v' + const NT r = _conic.r(); + const NT s = _conic.s(); + const NT t = _conic.t(); + const NT u = _conic.u(); + const NT v = _conic.v(); + const NT w = _conic.w(); + NT A, B, C, D, E; + + if (arc._conic.s() != _zero) + { + const NT s_tag = arc._conic.s(); + + A = s_tag*r - s*arc._conic.r(); + B = s_tag*u - s*arc._conic.u(); + C = s_tag*w - s*arc._conic.w(); + D = s_tag*t - s*arc._conic.t(); + E = s_tag*v - s*arc._conic.v(); + } + else + { + A = arc._conic.r(); + B = arc._conic.u(); + C = arc._conic.w(); + D = arc._conic.t(); + E = arc._conic.v(); + } + + if (D == _zero && E == _zero) + { + // In this case: A*x^2 + B*x + C = 0, so: + x_deg = 2; + x_coeffs.resize(3); + + x_coeffs[2] = A; + x_coeffs[1] = B; + x_coeffs[0] = C; + + // Solve the equation. + n_roots = solve_quadratic_eq (x_coeffs[2], + x_coeffs[1], + x_coeffs[0], + xs, x_mults); + n_approx = 0; + } + else + { + // Now we have: + // + // A*x^2 + B*x + C + // y = - ----------------- + // D*x + E + // + // Applying this two our conic's equation yields a quartic equation: + // + // c[4]*x^4 + c[3]*x^3 + c[2]*x^2 + c[1]*x + c[0] = 0 + const NT _two = 2; + + // Set the generating polynomial. + x_deg = 4; + x_coeffs.resize(5); + + if (t == _zero && arc._conic.t() == _zero) + { + x_coeffs[4] = s*A*A; + x_coeffs[3] = _two*s*A*B; + x_coeffs[2] = r*E*E + _two*s*A*C + s*B*B - v*A*E; + x_coeffs[1] = u*E*E + _two*s*B*C - v*B*E; + x_coeffs[0] = w*E*E + s*C*C - v*C*E; + } + else + { + const NT F = t*E + v*D; + + x_coeffs[4] = r*D*D + s*A*A - t*A*D; + x_coeffs[3] = _two*r*D*E + u*D*D + _two*s*A*B - t*B*D - F*A; + x_coeffs[2] = r*E*E + _two*u*D*E + w*D*D + + _two*s*A*C + s*B*B - t*C*D - F*B - v*A*E; + x_coeffs[1] = u*E*E + _two*w*D*E + _two*s*B*C - F*C - v*B*E; + x_coeffs[0] = w*E*E + s*C*C - v*C*E; + } + + // Try to factor out intersection points that occur at vertical + // tangency points of one of the curves. + Point_2 vpts[4]; + int n_vpts; + NT red_coeffs[5], temp_coeffs[5]; + int red_degree = x_deg, temp_degree; + int k, l; + bool reduced = false; + bool first_time; + + n_vpts = this->_conic_vertical_tangency_points (vpts); + n_vpts += arc._conic_vertical_tangency_points (vpts + n_vpts); + + for (l = 0; l <= x_deg; l++) + red_coeffs[l] = x_coeffs[l]; + + n_roots = 0; + for (k = 0; k < n_vpts; k++) + { + first_time = true; + while (factor_root (red_coeffs, red_degree, + vpts[k].x(), + temp_coeffs, temp_degree)) + { + if (first_time) + { + xs[n_roots] = vpts[k].x(); + x_mults[n_roots] = 1; + n_roots++; + first_time = false; + } + else + { + x_mults[n_roots-1]++; + } + + red_coeffs[red_degree] = _zero; + for (l = 0; l <= temp_degree; l++) + red_coeffs[l] = temp_coeffs[l]; + red_degree = temp_degree; + reduced = true; + } + } + + if (reduced) + { + n_roots += solve_quartic_eq (red_coeffs[4], + red_coeffs[3], + red_coeffs[2], + red_coeffs[1], + red_coeffs[0], + xs + n_roots, x_mults + n_roots, + n_approx); + } + else + { + // Solve the quartic equation. + n_roots = solve_quartic_eq (x_coeffs[4], + x_coeffs[3], + x_coeffs[2], + x_coeffs[1], + x_coeffs[0], + xs, x_mults, + n_approx); + } + } + } + + return (n_roots); + } + + // Calculate all y co-ordinates of intersection points between the two + // base curves of (*this) and the given arc. + int _y_coordinates_of_intersections_with (const Conic_arc_2& arc, + NT* ys, int* y_mults, + int& n_approx, + int& y_deg, + std::vector& y_coeffs) const + { + static const NT _zero = 0; + int n_roots; // The number of distinct y values. + + if ((_info & DEGREE_MASK) == DEGREE_1 && + (arc._info & DEGREE_MASK) == DEGREE_1) + { + // The two conics are: u*x + v*y + w = 0 + // and: u'*x + v'*y + w' = 0 + // There's a single solution for y, which is: + const NT denom = _conic.u()*arc._conic.w() - _conic.w()*arc._conic.u(); + const NT numer = _conic.v()*arc._conic.u() - _conic.u()*arc._conic.v(); + + if (numer == _zero) + { + n_roots = 0; + } + else + { + ys[0] = denom / numer; + y_mults[0] = 1; + n_roots = 1; + n_approx = 0; + } + + // Do not store coefficients information in this case. + y_deg = 0; + } + else if ((arc._info & DEGREE_MASK) == DEGREE_1) + { + // The two conics are: r*x^2 + s*y^2 + t*xy + u*x + v*y + w = 0 + // and: a*x + b*y + c = 0 + // There are therefore 2 possible y-values, the solutions for: + const NT a = arc._conic.u(); + const NT b = arc._conic.v(); + const NT c = arc._conic.w(); + const NT r = _conic.r(); + const NT s = _conic.s(); + const NT t = _conic.t(); + const NT u = _conic.u(); + const NT v = _conic.v(); + const NT w = _conic.w(); + + // Store the generating polynomial information in this case. + y_deg = 2; + y_coeffs.resize(3); + y_coeffs[2] = b*b*r + a*a*s - a*b*t; + y_coeffs[1] = 2*b*c*r + a*a*v - a*b*u - a*c*t; + y_coeffs[0] = c*c*r + a*a*w - a*c*u; + + // Solve the equation. + n_roots = solve_quadratic_eq (y_coeffs[2], + y_coeffs[1], + y_coeffs[0], + ys, y_mults); + n_approx = 0; + } + else if ((_info & IS_CIRCLE) != 0 && (arc._info & IS_CIRCLE) != 0) + { + // Special treatment for two circles. + // The two curves are: r*x^2 + r*y^2 + u*x + v*y + w = 0 + // and: r'*x^2 + r'*y^2 + u'*x + v'*y + w' = 0 + // + // Thus, r'*C1-r*C2 is a line whose equation is: a*x + b*y + = 0, where: + const NT r = _conic.r(); + const NT s = _conic.s(); + const NT t = _conic.t(); + const NT u = _conic.u(); + const NT v = _conic.v(); + const NT w = _conic.w(); + const NT a = arc._conic.r()*u - r*arc._conic.u(); + const NT b = arc._conic.r()*v - r*arc._conic.v(); + const NT c = arc._conic.r()*w - r*arc._conic.w(); + + if (a == _zero) + { + // The line b*y + c = 0 connects both intersection points of the two + // circles, so the both have a y-coordinate of -c/b. + if (b == _zero) + { + n_roots = 0; + n_approx = 0; + } + else + { + ys[0] = -c / b; + y_mults[0] = 2; + n_roots = 1; + n_approx = 0; + } + + // Do not set any coefficients' information in this case. + y_deg = 0; + } + else + { + // The intersection points of the two circles are the same as the + // intersection points of one of the circles with a*x + b*y + c = 0. + // Set the generating polynomial. + y_deg = 2; + y_coeffs.resize(3); + y_coeffs[2] = b*b*r + a*a*s - a*b*t; + y_coeffs[1] = 2*b*c*r + a*a*v - a*b*u - a*c*t; + y_coeffs[0] = c*c*r + a*a*w - a*c*u; + + // Find the roots. + n_roots = solve_quadratic_eq (y_coeffs[2], + y_coeffs[1], + y_coeffs[0], + ys, y_mults); + n_approx = 0; + } + } + else if (_conic.r() == _zero && _conic.t() == _zero && + arc._conic.r() == _zero && arc._conic.t() == _zero) + { + // Special treatment for canonic parabolas whose axes are parallel + // to the x axis. + // The two curves are: s*y^2 + u*x + v*y + w = 0 + // and: s'*y^2 + u'*x + v'*y + w' = 0 + // There are therefore 2 possible y-values, the solutions for: + + // Store the generating polynomial information in this case. + y_deg = 2; + y_coeffs.resize(3); + y_coeffs[2] = _conic.s()*arc._conic.u() - _conic.u()*arc._conic.s(); + y_coeffs[1] = _conic.v()*arc._conic.u() - _conic.u()*arc._conic.v(); + y_coeffs[0] = _conic.w()*arc._conic.u() - _conic.u()*arc._conic.w(); + + // Solve the equation: + n_roots = solve_quadratic_eq (y_coeffs[2], + y_coeffs[1], + y_coeffs[0], + ys, y_mults); + n_approx = 0; + } + else + { + // If the two conics are: r*x^2 + s*y^2 + t*xy + u*x + v*y + w = 0 + // and: r'*x^2 + s'*y^2 + t'*xy + u'*x + v'*y + w' = 0 + // Let us define: + // A = r'*s - r*s' if (r' != 0), otherwise A = s' + // B = r'*v - r*v' if (r' != 0), otherwise B = v' + // C = r'*w - r*w' if (r' != 0), otherwise C = w' + // D = r'*t - r*t' if (r' != 0), otherwise A = t' + // E = r'*u - r*u' if (r' != 0), otherwise A = u' + const NT r = _conic.r(); + const NT s = _conic.s(); + const NT t = _conic.t(); + const NT u = _conic.u(); + const NT v = _conic.v(); + const NT w = _conic.w(); + NT A, B, C, D, E; + + if (arc._conic.r() != _zero) + { + const NT r_tag = arc._conic.r(); + + A = r_tag*s - r*arc._conic.s(); + B = r_tag*v - r*arc._conic.v(); + C = r_tag*w - r*arc._conic.w(); + D = r_tag*t - r*arc._conic.t(); + E = r_tag*u - r*arc._conic.u(); + } + else + { + A = arc._conic.s(); + B = arc._conic.v(); + C = arc._conic.w(); + D = arc._conic.t(); + E = arc._conic.u(); + } + + if (D == _zero && E == _zero) + { + // In this case: A*y^2 + B*y + C = 0, so: + y_deg = 2; + y_coeffs.resize(3); + + y_coeffs[2] = A; + y_coeffs[1] = B; + y_coeffs[0] = C; + + // Solve the equation. + n_roots = solve_quadratic_eq (y_coeffs[2], + y_coeffs[1], + y_coeffs[0], + ys, y_mults); + n_approx = 0; + } + else + { + // Now we have: + // + // A*y^2 + B*y + C + // x = - ----------------- + // D*y + E + // + // Applying this two our conic's equation yields a quartic equation: + // + // c[4]*y^4 + c[3]*y^3 + c[2]*y^2 + c[1]*y + c[0] = 0 + const NT _two = 2; + + // Store the generating polynomial information in this case. + y_deg = 4; + y_coeffs.resize(5); + + if (t == _zero && arc._conic.t() == _zero) + { + y_coeffs[4] = r*A*A; + y_coeffs[3] = _two*r*A*B; + y_coeffs[2] = s*E*E + _two*r*A*C + r*B*B - u*A*E; + y_coeffs[1] = v*E*E + _two*r*B*C - u*B*E; + y_coeffs[0] = w*E*E + r*C*C - u*C*E; + } + else + { + const NT F = t*E + u*D; + + y_coeffs[4] = s*D*D + r*A*A - t*A*D; + y_coeffs[3] = _two*s*D*E + v*D*D + _two*r*A*B - t*B*D - F*A; + y_coeffs[2] = s*E*E + _two*v*D*E + w*D*D + + _two*r*A*C + r*B*B - t*C*D - F*B - u*A*E; + y_coeffs[1] = v*E*E + _two*w*D*E + _two*r*B*C - F*C - u*B*E; + y_coeffs[0] = w*E*E + r*C*C - u*C*E; + } + + // Try to factor out intersection points that occur at vertical + // tangency points of one of the curves. + Point_2 vpts[4]; + int n_vpts; + NT red_coeffs[5], temp_coeffs[5]; + int red_degree = y_deg, temp_degree; + int k, l; + bool reduced = false; + bool first_time; + + n_vpts = this->_conic_vertical_tangency_points (vpts); + n_vpts += arc._conic_vertical_tangency_points (vpts + n_vpts); + + for (l = 0; l <= y_deg; l++) + red_coeffs[l] = y_coeffs[l]; + + n_roots = 0; + for (k = 0; k < n_vpts; k++) + { + first_time = true; + while (factor_root (red_coeffs, red_degree, + vpts[k].y(), + temp_coeffs, temp_degree)) + { + if (first_time) + { + ys[n_roots] = vpts[k].y(); + y_mults[n_roots] = 1; + n_roots++; + first_time = false; + } + else + { + y_mults[n_roots-1]++; + } + + red_coeffs[red_degree] = _zero; + for (l = 0; l <= temp_degree; l++) + red_coeffs[l] = temp_coeffs[l]; + red_degree = temp_degree; + reduced = true; + } + } + + if (reduced) + { + n_roots += solve_quartic_eq (red_coeffs[4], + red_coeffs[3], + red_coeffs[2], + red_coeffs[1], + red_coeffs[0], + ys + n_roots, y_mults + n_roots, + n_approx); + } + else + { + // Solve the equation: + n_roots = solve_quartic_eq (y_coeffs[4], + y_coeffs[3], + y_coeffs[2], + y_coeffs[1], + y_coeffs[0], + ys, y_mults, + n_approx); + } + } + } + + return (n_roots); + } + + // Pair the x coordinates and the y coordinates of the intersection point + // of (*this) and arc and return a vector of intersection points. + int _pair_intersection_points (const Conic_arc_2& arc, + const int& n_xs, + const NT* xs, const int* x_mults, + const int& n_approx_xs, + const int& n_ys, + const NT* ys, const int* y_mults, + const int& n_approx_ys, + Point_2* ipts) const + { + // Calculate the minimal number of intersection points. + int x_total = 0, y_total = 0; + int min_points; + int i, j; + + for (i = 0; i < n_xs; i++) + x_total += x_mults[i]; + + for (j = 0; j < n_ys; j++) + y_total += y_mults[j]; + + min_points = (x_total < y_total) ? x_total : y_total; + + // Go over all x coordinates. + const APNT r1 = TO_APNT(_conic.r()); + const APNT s1 = TO_APNT(_conic.s()); + const APNT t1 = TO_APNT(_conic.t()); + const APNT u1 = TO_APNT(_conic.u()); + const APNT v1 = TO_APNT(_conic.v()); + const APNT w1 = TO_APNT(_conic.w()); + const APNT r2 = TO_APNT(arc._conic.r()); + const APNT s2 = TO_APNT(arc._conic.s()); + const APNT t2 = TO_APNT(arc._conic.t()); + const APNT u2 = TO_APNT(arc._conic.u()); + const APNT v2 = TO_APNT(arc._conic.v()); + const APNT w2 = TO_APNT(arc._conic.w()); + APNT x, y; + int n_ipts = 0; + APNT results[16], res; + int x_indices[16], y_indices[16], temp; + bool is_approx; + int k, l; + + k = 0; + for (i = 0; i < n_xs; i++) + { + // For each y, calculate: c[j] = |C1(x[i],y[j])| + |C2(x[i],y[j])|. + x = TO_APNT(xs[i]); + for (j = 0; j < n_ys; j++) + { + y = TO_APNT(ys[j]); + results[k] = APNT_ABS(r1*x*x + s1*y*y + t1*x*y + u1*x + v1*y + w1) + + APNT_ABS(r2*x*x + s2*y*y + t2*x*y + u2*x + v2*y + w2); + x_indices[k] = i; + y_indices[k] = j; + k++; + } + } + + // Perform bubble-sort. + for (k = 0; k < n_xs*n_ys - 1; k++) + { + for (l = k + 1; l < n_xs*n_ys; l++) + { + if (results[k] > results[l]) + { + res = results[k]; + results[k] = results[l]; + results[l] = res; + + temp = x_indices[k]; + x_indices[k] = x_indices[l]; + x_indices[l] = temp; + temp = y_indices[k]; + y_indices[k] = y_indices[l]; + y_indices[l] = temp; + } + } + } + + // Use the x,y pairs for which the c[j] values are the lowest. + n_ipts = 0; + k = 0; + while (k < n_xs*n_ys && n_ipts < 4) + { + if (n_ipts >= min_points && + eps_compare(results[k]*results[k], 0) != EQUAL) + { + break; + } + + is_approx = (x_indices[k] >= (n_xs - n_approx_xs)) || + (y_indices[k] >= (n_ys - n_approx_ys)); + + ipts[n_ipts] = Point_2 (xs[x_indices[k]], ys[y_indices[k]], + (is_approx ? + Point_2::Intersection_approx : + Point_2::Intersection_exact), + _conic_id, arc._conic_id); + k++; + n_ipts++; + } + + return (n_ipts); + } +}; + +#ifndef NO_OSTREAM_INSERT_CONIC_ARC_2 +template +std::ostream& operator<< (std::ostream& os, const Conic_arc_2& arc) +{ + typename Conic_arc_2::Conic_2 conic = arc.conic(); + typename Conic_arc_2::Point_2 source = arc.source(), + target = arc.target(); + + os << "{" << CGAL::to_double(conic.r()) << "*x^2 + " + << CGAL::to_double(conic.s()) << "*y^2 + " + << CGAL::to_double(conic.t()) << "*xy + " + << CGAL::to_double(conic.u()) << "*x + " + << CGAL::to_double(conic.v()) << "*y + " + << CGAL::to_double(conic.w()) << "} :" + << "(" << CGAL::to_double(source.x()) << "," + << CGAL::to_double(source.y()) << ") -> " + << "(" << CGAL::to_double(target.x()) << "," + << CGAL::to_double(target.y()) << ")"; + + return (os); +} +#endif // NO_OSTREAM_INSERT_CONIC_ARC_2 + +CGAL_END_NAMESPACE + +#endif diff --git a/Packages/Arrangement/include/CGAL/Arrangement_2/Conic_arc_2_eq.h b/Packages/Arrangement/include/CGAL/Arrangement_2/Conic_arc_2_eq.h new file mode 100644 index 00000000000..586a5023885 --- /dev/null +++ b/Packages/Arrangement/include/CGAL/Arrangement_2/Conic_arc_2_eq.h @@ -0,0 +1,996 @@ +#ifndef CGAL_CONIC_ARC_2_EQ_H +#define CGAL_CONIC_ARC_2_EQ_H + +#include +#include + +CGAL_BEGIN_NAMESPACE + +typedef double APNT; +#define APNT_EPS 0.0000001 +#define COMP_EPS 0.0001 +#define TO_APNT(x) CGAL::to_double(x) +#define APNT_ABS(x) fabs(x) + +// ---------------------------------------------------------------------------- +// Compare two approximate values. +// +template +Comparison_result eps_compare (const APNT& val1, const APNT& val2, + const APNT& _eps = APNT_EPS) +{ + APNT denom = APNT_ABS(val1 - val2); + APNT numer = APNT_ABS(val1) + APNT_ABS(val2); + + if (numer < CGAL::sqrt(_eps)) + { + // The two numbers are very close to zero: + if (denom < _eps) + return (EQUAL); + else + return ((val1 < val2) ? SMALLER : LARGER); + } + else + { + // | x - y | + // x = y <==> ----------- < _eps + // |x| + |y| + if (denom / numer < _eps) + return (EQUAL); + else + return ((val1 < val2) ? SMALLER : LARGER); + } +} + +/*****************************************************************************/ +// Some functions for quick-solving equations when working with doubles. + +//-------------------------------------------------------------------- +// Solve the cubic equation: a*x^3 + b*x^2 + c*x + d = 0 (a != 0), +// with double precision floats and using Tartaglia's method. +// The equation is assumed not to have multiple roots. +// The function returns the number of real solutions. +// The roots area must be at least of size 3. +// +template +int _Tartaglia_cubic_eq (const NT& a, + const NT& b, + const NT& c, + const NT& d, + NT *roots) +{ + // Some constants. + static const NT _Zero = 0; + static const NT _One = 1; + static const NT _Two = 2; + static const NT _Three = 3; + static const NT _Nine = 9; + static const NT _Twenty_seven = 27; + static const NT _Fifty_four = 54; + static const NT _One_third = _One / _Three; + + // Normalize the equation to obtain: x^3 + a2*x^2 + a1*x^1 + a0 = 0. + const NT a2 = b / a; + const NT a1 = c / a; + const NT a0 = d / a; + + // Caclculate the "cubic discriminant" D. + const NT Q = (_Three*a1 - a2*a2) / + _Nine; + const NT R = (_Nine*a1*a2 - _Twenty_seven*a0 - _Two*a2*a2*a2) / + _Fifty_four; + const NT D = Q*Q*Q + R*R; + NT S, T; + + if (D > _Zero) + { + // One real-valued root (the other two are complex conjugates): + NT temp = R + CGAL::sqrt(D); + + if (temp > _Zero) + S = ::pow(temp, _One_third); + else + S = - ::pow(-temp, _One_third); + + temp = R - CGAL::sqrt(D); + if (temp > _Zero) + T = ::pow(temp, _One_third); + else + T = - ::pow(-temp, _One_third); + + roots[0] = S + T - a2/_Three; + return (1); + } + else + { + // We have three real-valued roots: + const NT phi = ::atan2 (CGAL::sqrt(-D), R) / _Three; + const NT sin_phi = ::sin(phi); + const NT cos_phi = ::cos(phi); + + roots[0] = _Two*CGAL::sqrt(-Q)*cos_phi - a2/_Three; + roots[1] = -CGAL::sqrt(-Q)*cos_phi - + CGAL::sqrt(-_Three*Q)*sin_phi - a2/_Three; + roots[2] = -CGAL::sqrt(-Q)*cos_phi + + CGAL::sqrt(-_Three*Q)*sin_phi - a2/_Three; + + return (3); + } +} + +//-------------------------------------------------------------------- +// Calculate the square root of a complex number. +// +template +void _complex_sqrt (const NT& real, const NT& imag, + NT& sqrt_re, NT& sqrt_im, + const NT& _eps) +{ + static const NT _Zero = 0; + static const NT _Two = 2; + + // Check the special case where imag=0: + if (eps_compare(imag, _Zero, _eps) == EQUAL) + { + if (real >= _Zero) + { + sqrt_re = CGAL::sqrt(real); + sqrt_im = _Zero; + } + else + { + sqrt_re = _Zero; + sqrt_im = CGAL::sqrt(-real); + } + + return; + } + + // We seek x,y such that: (x^2 - y^2 = 0) and (2xy = imag). + // So, we get the equation: x^4 - real*x^2 - (imag/2)^2 = 0. + // This equation has one positive root (and one irrelevant negative root). + const NT sqrt_disc = CGAL::sqrt(real*real + imag*imag); + const NT x_2 = (real + sqrt_disc) / _Two; + + // Calculate on of the roots (the other root is -sqrt_re - i*sqrt_im). + sqrt_re = CGAL::sqrt(x_2); + sqrt_im = imag / (_Two*sqrt_re); + + return; +} + +//-------------------------------------------------------------------- +// Solve the quartic equation: a*x^4 + b*x^3 + c*x^2 + d*x + e = 0 (a != 0), +// with double precision floats and using Ferrari's method. +// The equation is assumed not to have multiple roots. +// The function returns the number of real solutions. +// The roots area must be at least of size 4. +// +template +int _Ferrari_quartic_eq (const NT& a, + const NT& b, + const NT& c, + const NT& d, + const NT& e, + NT *roots, + const NT& _eps) +{ + // Some constants. + static const NT _Zero = 0; + static const NT _One = 1; + static const NT _Two = 2; + static const NT _Three = 3; + static const NT _Four = 4; + static const NT _Eight = 8; + static const NT _One_quarter = _One/_Four; + static const NT _One_half = _One/_Two; + static const NT _Three_quarters = _Three/_Four; + + // Normalize the equation to obtain: x^4 + a3*x^3 + a2*x^2 + a1*x^1 + a0 = 0. + const NT a3 = b / a; + const NT a2 = c / a; + const NT a1 = d / a; + const NT a0 = e / a; + + // Find a real root of the resolvent cubic equation: + // y^3 - a2*y^2 + (a1*a3 - 4*a0)*y + (2*a2*a0 - a1^2 - a0*a3^2) = 0. + NT cubic_roots[3]; + NT y1; + + _Tartaglia_cubic_eq (_One, + -a2, + a1*a3 - _Four*a0, + _Four*a2*a0 - a1*a1 - a3*a3*a0, + cubic_roots); + y1 = cubic_roots[0]; + + // Calculate R^2 and act accordingly. + const NT R_square = _One_quarter*a3*a3 - a2 + y1; + int n_roots = 0; + + if (eps_compare(R_square, _Zero, _eps) == EQUAL) + { + // R is _Zero: In this case, D and E must be real-valued. + const NT disc_square = y1*y1 - _Four*a0; + + if (disc_square < _Zero) + { + // No real roots at all. + return (0); + } + + const NT disc = CGAL::sqrt(disc_square); + const NT D_square = _Three_quarters*a3*a3 + _Two*(disc - a2); + const NT E_square = _Three_quarters*a3*a3 - _Two*(disc + a2); + + if (D_square > _Zero) + { + // Add two real-valued roots. + const NT D = CGAL::sqrt(D_square); + + roots[n_roots] = -_One_quarter*a3 + _One_half*D; + n_roots++; + roots[n_roots] = -_One_quarter*a3 - _One_half*D; + n_roots++; + } + + if (E_square > _Zero) + { + // Add two real-valued roots. + const NT E = CGAL::sqrt(E_square); + + roots[n_roots] = -_One_quarter*a3 + _One_half*E; + n_roots++; + roots[n_roots] = -_One_quarter*a3 - _One_half*E; + n_roots++; + } + } + else if (R_square > _Zero) + { + // R is real-valued: In this case, D and E must also be real-valued. + const NT R = CGAL::sqrt(R_square); + const NT D_square = _Three_quarters*a3*a3 - R_square - _Two*a2 + + (_Four*a3*a2 - _Eight*a1 - a3*a3*a3) / (_Four*R); + const NT E_square = _Three_quarters*a3*a3 - R_square - _Two*a2 - + (_Four*a3*a2 - _Eight*a1 - a3*a3*a3) / (_Four*R); + + if (D_square > _Zero) + { + // Add two real-valued roots. + const NT D = CGAL::sqrt(D_square); + + roots[n_roots] = -_One_quarter*a3 + _One_half*(R + D); + n_roots++; + roots[n_roots] = -_One_quarter*a3 + _One_half*(R - D); + n_roots++; + } + + if (E_square > _Zero) + { + // Add two real-valued roots. + const NT E = CGAL::sqrt(E_square); + + roots[n_roots] = -_One_quarter*a3 + _One_half*(E - R); + n_roots++; + roots[n_roots] = -_One_quarter*a3 - _One_half*(E + R); + n_roots++; + } + } + else + { + // R is imaginary: Calclualte D and E as complex numbers. + const NT R_im = CGAL::sqrt(-R_square); + const NT R_inv_im = -_One / R_im; // The inverse of R. + const NT D_square_re = _Three_quarters*a3*a3 - R_square - _Two*a2; + const NT D_square_im = + (a3*a2 - _Two*a1 - _One_quarter*a3*a3*a3) * R_inv_im; + NT D_re, D_im; + const NT E_square_re = D_square_re; + const NT E_square_im = -D_square_im; + NT E_re, E_im; + + _complex_sqrt (D_square_re, D_square_im, + D_re, D_im, _eps); + _complex_sqrt (E_square_re, E_square_im, + E_re, E_im, _eps); + + // Use only the combinations that yield real-valued numbers. + if (eps_compare(R_im, -D_im, _eps) == EQUAL) // (R_im == -D_im) + { + roots[n_roots] = -_One_quarter*a3 + _One_half*D_re; + n_roots++; + } + if (eps_compare(R_im, D_im, _eps) == EQUAL) // (R_im == D_im) + { + roots[n_roots] = -_One_quarter*a3 - _One_half*D_re; + n_roots++; + } + + if (eps_compare(R_im, E_im, _eps) == EQUAL) // (R_im == E_im) + { + roots[n_roots] = -_One_quarter*a3 + _One_half*E_re; + n_roots++; + } + if (eps_compare(R_im, -E_im, _eps) == EQUAL) // (R_im == -E_im) + { + roots[n_roots] = -_One_quarter*a3 - _One_half*D_re; + n_roots++; + } + } + + return (n_roots); +} + +/*****************************************************************************/ + +// ---------------------------------------------------------------------------- +// Solve the quadratic equation: a*x^2 + b*x + c = 0. +// The function returns the number of distinct solutions. +// The roots area must be at least of size 2. +// +template +int solve_quadratic_eq (const NT& a, const NT& b, const NT& c, + NT* roots, int* mults) +{ + static const NT _zero = 0; + static const NT _two = 2; + static const NT _four = 4; + + if (a == _zero) + { + // Solve a linear equation. + if (b != _zero) + { + roots[0] = -c/b; + mults[0] = 1; + return (1); + } + else + return (0); + } + + // Act according to the discriminant. + const NT disc = b*b - _four*a*c; + + if (disc < _zero) + { + // No real roots. + return (0); + } + else if (disc == _zero) + { + // One real root with mutliplicity 2. + roots[0] = roots[1] = -b/(_two*a); + mults[0] = 2; + return (1); + } + else + { + // Two real roots. + const NT sqrt_disc = CGAL::sqrt(disc); + + roots[0] = (sqrt_disc - b)/(_two*a); + mults[0] = 1; + roots[1] = -(sqrt_disc + b)/(_two*a); + mults[1] = 1; + return (2); + } +} + +// ---------------------------------------------------------------------------- +// Solve the cubic equation: _a*x^3 + _b*x^2 + _c*x + _d = 0. +// Notice this is an auxiliary function (used only by solve_quartic_eq()). +// The function returns the number of distinct solutions. +// The roots area must be at least of size 3. +// +template +static int _solve_cubic_eq (const NT& _a, const NT& _b, + const NT& _c, const NT& _d, + NT* roots, int* mults, + int& n_approx) +{ + n_approx = 0; + + // In case we have an equation of a lower degree: + static const NT _zero = 0; + + if (_a == _zero) + { + // We have to solve _b*x^2 + _c*x + _d = 0. + return (solve_quadratic_eq (_b, _c, _d, + roots, mults)); + } + + // Normalize the equation, so that the leading coefficient is 1 and we have: + // x^3 + b*x^2 + c*x + d = 0 + // + const NT b = _b/_a; + const NT c = _c/_a; + const NT d = _d/_a; + + // Check whether the equation has multiple roots. + // If we write: + // p(x) = x^3 + b*x^2 + c*x + d = 0 + // + // Then: + // p'(x) = 3*x^2 + 2*b*x + c + // + // We know that there are multiple roots iff GCD(p(x), p'(x)) != 0. + // In order to check the GCD, let us denote: + // p(x) mod p'(x) = A*x + B + // + // Then: + static const NT _two = 2; + static const NT _three = 3; + static const NT _nine = 9; + + const NT A = _two*(c - b*b/_three)/_three; + const NT B = d - b*c/_nine; + + if (A == _zero) + { + // In case A,B == 0, then p'(x) divides p(x). + // This means the equation has one solution with multiplicity of 3. + // We can obtain this root using the fact that -b is the sum of p(x)'s + // roots. + if (B == _zero) + { + roots[0] = -b / _three; + mults[0] = 3; + return (1); + } + } + else + { + // Check whether A*x + B divides p'(x). + const NT x0 = -B/A; + + if (c + x0*(_two*b + _three*x0) == _zero) + { + // Since GCD(p(x), p'(x)) = (x - x0), then x0 is a root of p(x) with + // multiplicity 2. The other root is obtained from b. + roots[0] = x0; + mults[0] = 2; + roots[1] = -(b + 2*x0); + mults[1] = 1; + return (2); + } + } + + // If we reached here, we can be sure that there are no multiple roots. + // We use Ferrari's method to approximate the solutions. + APNT app_roots[4]; + + n_approx = _Tartaglia_cubic_eq (1, + TO_APNT(b), + TO_APNT(c), + TO_APNT(d), + app_roots); + + // Copy the approximations. + for (int i = 0; i < n_approx; i++) + { + roots[i] = NT(app_roots[i]); + mults[i] = 1; + } + + return (n_approx); +} + +// ---------------------------------------------------------------------------- +// Solve a quartic equation: _a*x^4 + _b*x^3 + _c*x^2 + _d*x + _e = 0. +// The function returns the number of distinct solutions. +// The roots area must be at least of size 4. +// +template +int solve_quartic_eq (const NT& _a, const NT& _b, const NT& _c, + const NT& _d, const NT& _e, + NT* roots, int* mults, + int& n_approx) +{ + // Right now we have no approximate solutions. + n_approx = 0; + + // First check whether we have roots that are 0. + static const NT _zero = 0; + + if (_e == _zero) + { + if (_d == _zero) + { + if (_c == _zero) + { + if (_b == _zero) + { + // Unless the equation is trivial, we have the root 0, + // with multiplicity of 4. + if (_a == _zero) + return (0); + + roots[0] = _zero; + mults[0] = 4; + return (1); + } + else + { + // Add the solution 0 (with multiplicity 3), and add the solution + // to _a*x + _b = 0. + if (_a == _zero) + return (0); + + roots[0] = _zero; + mults[0] = 3; + roots[1] = -(_b/_a); + mults[1] = 1; + return (2); + } + } + else + { + // Add the solution 0 (with multiplicity 2), + // and solve _a*x^2 + _b*x + _c*x = 0. + roots[0] = _zero; + mults[0] = 2; + + return (1 + solve_quadratic_eq (_a, _b, _c, + roots + 1, mults + 1)); + } + } + else + { + // Add the solution 0, and solve _a*x^3 + _b*x^2 + _c*x + _d = 0. + roots[0] = _zero; + mults[0] = 1; + + return (1 + _solve_cubic_eq (_a, _b, _c, _d, + roots + 1, mults + 1, + n_approx)); + } + } + + // In case we have an equation of a lower degree: + if (_a == _zero) + { + if (_b == _zero) + { + // We have to solve _c*x^2 + _d*x + _e = 0. + return (solve_quadratic_eq (_c, _d, _e, + roots, mults)); + } + else + { + // We have to solve _b*x^3 + _c*x^2 + _d*x + _e = 0. + return (_solve_cubic_eq (_b, _c, _d, _e, + roots, mults, + n_approx)); + } + } + + // In case we have an equation of the form: + // _a*x^4 + _c*x^2 + _e = 0 + // + // Then by substituting y = x^2, we obtain a quadratic equation: + // _a*y^2 + _c*y + _e = 0 + // + if (_b == _zero && _d == _zero) + { + // Solve the equation for y = x^2. + NT sq_roots[2]; + int sq_mults[2]; + int n_sq; + int n = 0; + + n_sq = solve_quadratic_eq (_a, _c, _e, + sq_roots, sq_mults); + + // Convert to roots of the original equation: only positive roots for y + // are relevant of course. + for (int j = 0; j < n_sq; j++) + { + if (sq_roots[j] > _zero) + { + roots[n] = CGAL::sqrt(sq_roots[j]); + mults[n] = sq_mults[j]; + roots[n+1] = -roots[n]; + mults[n+1] = sq_mults[j]; + n += 2; + } + } + + return (n); + } + + // Normalize the equation, so that the leading coefficient is 1 and we have: + // x^4 + b*x^3 + c*x^2 + d*x + e = 0 + // + const NT b = _b/_a; + const NT c = _c/_a; + const NT d = _d/_a; + const NT e = _e/_a; + + // Check whether the equation has multiple roots. + // If we write: + // p(x) = x^4 + b*x^3 + c*x^2 + d*x + e = 0 + // + // Then: + // p'(x) = 4*x^3 + 3*b*x^2 + 2*c*x + d + // + // We know that there are multiple roots iff GCD(p(x), p'(x)) != 0. + // In order to check the GCD, let us denote: + // p(x) mod p'(x) = alpha*x^2 + beta*x + gamma + // + // Then: + static const NT _two = 2; + static const NT _three = 3; + static const NT _four = 4; + static const NT _eight = 8; + static const NT _sixteen = 16; + + const NT alpha = c/_two - _three*b*b/_sixteen; + const NT beta = _three*d/_four - b*c/_eight; + const NT gamma = e - b*d/_sixteen; + + if (alpha == _zero) + { + if (beta == _zero) + { + if (gamma == _zero) + { + // p'(x) divides p(x). This can happen only if there is a single + // root with multiplicity 4. + // Since b = -(sum of p(x)'s roots) then this root is -b/4. + roots[0] = -b/_four; + mults[0] = 4; + return (1); + } + } + else + { + // In case alpha is 0, check whether (beta*x + gamma) divides p'(x) - + // this happens iff -gamma/beta is a root of p'(x). + if (beta == _zero) + return (0); + + const NT x0 = -gamma/beta; + + if (d + x0*(_two*c + x0*(_three*b + x0*_four)) == _zero) + { + // We have two ditinct solutions: x0 with multiplicity of 3, and + // p'(x)/(x - x0)^2 with multiplicity of 1. + // Since b = -(sum of p(x)'s roots), the other root is -(b + 3*x0). + roots[0] = x0; + mults[0] = 3; + roots[1] = -(b + 3*x0); + mults[1] = 1; + return (2); + } + } + } + else + { + // If alpha is not 0, let us denote: + // p'(x) mod (alpha*x^2 + beta*x + gamma) = A*x + B + // + // And so: + const NT A = _two*c - (_four*gamma + _three*b*beta)/alpha + + _four*beta*beta/(alpha*alpha); + const NT B = d - _three*b*gamma/alpha + _four*beta*gamma/(alpha*alpha); + + // In case A,B == 0, then GCD(p(x), p'(x)) = (alpha*x^2 + beta*x + gamma). + // This means the equation has two distinct solution, each one of them + // has a multiplicity of 2. + if (A == _zero) + { + if (B == _zero) + { + // There are two cases: + // - If (alpha*x^2 + beta*x + gamma) has two distinct roots, that + // these roots are roots of p(x), each with multiplicity 2. + // - If there is one root x0, than this root is a root of p(x) with + // multiplicity 3. The other root can be obtained using the fact + // that b = -(sum of p(x)'s roots). + int m = solve_quadratic_eq (alpha, beta, gamma, + roots, mults); + + if (m == 1) + { + // Add the second root. + mults[0] = 3; + roots[1] = -(b + 3*roots[0]); + mults[1] = 1; + m++; + } + else if (m == 2) + { + mults[0] = 2; + mults[1] = 2; + } + + return (m); + } + } + else + { + // Check whether (A*x + B) divides (alpha*x^2 + beta*x + gamma). + // This happend iff -B/A is a root of (alpha*x^2 + beta*x + gamma). + const NT x0 = -B / A; + + if (gamma + x0*(beta + x0*alpha) == _zero) + { + // We have one solution with multiplicity of 2. + roots[0] = x0; + mults[0] = 2; + + // The other two solutions are obtained from solving: + // p(x) / (x - x0)^2 = x^2 + (b + 2*x0)*x + (c + 2*b*x0 + 3*x0^2) + return (1 + + solve_quadratic_eq (NT(1), b + 2*x0, c + x0*(2*b + 3*x0), + roots + 1, mults + 1)); + } + } + } + + // Compute the number of real roots we should compute, using Sturm sequences. + Polynom p; + p.set (4, _a); + p.set (3, _b); + p.set (2, _c); + p.set (1, _d); + p.set (0, _e); + + Sturm_seq sts (p); + int n_expected = sts.sign_changes_at_infinity(-1) - + sts.sign_changes_at_infinity(1); + + // If we reached here, we can be sure that there are no multiple roots. + // We use Ferrari's method to approximate the solutions. + APNT app_roots[4]; + APNT _eps = APNT_EPS; + + do + { + n_approx = _Ferrari_quartic_eq (1, + TO_APNT(b), + TO_APNT(c), + TO_APNT(d), + TO_APNT(e), + app_roots, + _eps); + + _eps *= 10; + } while (n_approx != n_expected && _eps < 1); + + CGAL_assertion(n_approx == n_expected); + + // Copy the approximations. + for (int i = 0; i < n_approx; i++) + { + roots[i] = NT(app_roots[i]); + mults[i] = 1; + } + + return (n_approx); +} + +// ---------------------------------------------------------------------------- +// Find the GCD of the two given polynomials. +// +template +void _mod_polynomials (const NT* a_coeffs, const int& a_deg, + const NT* b_coeffs, const int& b_deg, + NT* r_coeffs, int& r_deg) +{ + const static NT _zero = 0; + NT work[5]; + int deg = a_deg; + NT factor; + int k; + + // Copy the polynomial a(x) to the work area. + for (k = 0; k <= a_deg; k++) + work[k] = a_coeffs[k]; + + // Perform long division of a(x) by b(x). + while (deg >= b_deg) + { + factor = work[deg] / b_coeffs[b_deg]; + + for (k = b_deg - 1; k >= 0; k--) + work[deg - b_deg + k] -= factor*b_coeffs[k]; + + do + { + deg--; + } while (deg >= 0 && work[deg] == _zero); + } + + // Copy back the residue. + for (k = 0; k <= deg; k++) + r_coeffs[k] = work[k]; + r_deg = deg; + + return; +} + +template +void gcd_polynomials (const std::vector& a_coeffs, const int& a_deg, + const std::vector& b_coeffs, const int& b_deg, + std::vector& gcd_coeffs, int& gcd_deg) +{ + // Copy the polynomials to the work areas (make sure that A has a higher + // degree than B). + NT a[5], b[5], r[5]; + NT *aP = &(a[0]), *bP=&(b[0]), *rP = &(r[0]); + NT *swapP; + int adeg, bdeg, rdeg; + int k; + + if (a_deg >= b_deg) + { + for (k = 0; k <= a_deg; k++) + aP[k] = a_coeffs[k]; + adeg = a_deg; + + for (k = 0; k <= b_deg; k++) + bP[k] = b_coeffs[k]; + bdeg = b_deg; + } + else + { + for (k = 0; k <= b_deg; k++) + aP[k] = b_coeffs[k]; + adeg = b_deg; + + for (k = 0; k <= a_deg; k++) + bP[k] = a_coeffs[k]; + bdeg = a_deg; + } + + // Perform Euclid's algorithm. + while (true) + { + _mod_polynomials (aP, adeg, bP, bdeg, + rP, rdeg); + + if (rdeg <= 0) + break; + + swapP = aP; + aP = bP; + adeg = bdeg; + bP = rP; + bdeg = rdeg; + rP = swapP; + } + + if (rdeg < 0) + { + // b(x) divides a(x), so gcd(a,b) = b. + gcd_coeffs.resize(bdeg + 1); + for (k = 0; k <= bdeg; k++) + gcd_coeffs[k] = bP[k]; + gcd_deg = bdeg; + } + else + { + // The two polynomials have gcd(a,b) = 1. + gcd_deg = 0; + } + + return; +} + +// ---------------------------------------------------------------------------- +// Check whether alpha is a root of the given polynomial p(x). +// If so, compute: q(x) = p(x)/(x - alpha). +// +template +bool factor_root (const NT* p_coeffs, const int& p_deg, + const NT& alpha, + NT* q_coeffs, int& q_deg) +{ + const static NT _zero = 0; + NT coeff = p_coeffs[p_deg]; + int k; + + for (k = p_deg; k > 0; k--) + { + q_coeffs[k-1] = coeff; + coeff = p_coeffs[k-1] + alpha*coeff; + } + + if (coeff == _zero) + { + q_deg = p_deg - 1; + return (true); + } + else + { + q_deg = 0; + return (false); + } +} + +// ---------------------------------------------------------------------------- +// Check which one of the polynomial roots is greater than the other. +// +template +Comparison_result compare_roots (const std::vector& p1_coeffs, + const int& p1_deg, + const NT& alpha1, const NT& delta1, + const std::vector& p2_coeffs, + const int& p2_deg, + const NT& alpha2, const NT& delta2) +{ + NT l1 = alpha1 - delta1; + NT r1 = alpha1 + delta1; + NT l2 = alpha2 - delta2; + NT r2 = alpha2 + delta2; + + if (r1 < l2) + return (SMALLER); + else if (l1 > r2) + return (LARGER); + + // Generate Sturm sequences. + Polynom p1; + Polynom p2; + APNT max_coeff = 0; + int i; + + for (i = p1_deg; i >= 0; i--) + { + p1.set (i, p1_coeffs[i]); + if (TO_APNT(p1_coeffs[i]) > max_coeff) + max_coeff = TO_APNT(p1_coeffs[i]); + } + + for (i = p2_deg; i >= 0; i--) + { + p2.set (i, p2_coeffs[i]); + if (TO_APNT(p2_coeffs[i]) > max_coeff) + max_coeff = TO_APNT(p2_coeffs[i]); + } + + Sturm_seq sts1 (p1); + Sturm_seq sts2 (p2); + + // Determine the number of iterations. + const int max_deg = (p1_deg > p2_deg) ? p1_deg : p2_deg; + const int max_iters = max_deg * + static_cast(log(max_coeff) / log(2) + 1); + NT m1, m2; + const NT _two = NT(2); + + for (i = 0; i < max_iters; i++) + { + // Bisect [l1, r1] and select either [l1, m1] or [m1, r1]. + m1 = to_bigfloat((r1 + l1) / _two); + + if (sts1.sign_changes_at_x(l1) - sts1.sign_changes_at_x(m1) == 1) + r1 = m1; + else + l1 = m1; + + // Bisect [l2, r2] and select either [l2, m2] or [m2, r2]. + m2 = to_bigfloat((r2 + l2) / _two); + + if (sts2.sign_changes_at_x(l2) - sts2.sign_changes_at_x(m2) == 1) + r2 = m2; + else + l2 = m2; + + // If [l1, r1] and [l2, r2] do not overlap, return the comparison result. + if (r1 < l2) + return (SMALLER); + else if (l1 > r2) + return (LARGER); + } + + // According to Loos, if we reached here, the two roots are equal. + return (EQUAL); +} + +CGAL_END_NAMESPACE + +#endif diff --git a/Packages/Arrangement/include/CGAL/Arrangement_2/Conic_arc_2_point.h b/Packages/Arrangement/include/CGAL/Arrangement_2/Conic_arc_2_point.h new file mode 100644 index 00000000000..49572c86ab6 --- /dev/null +++ b/Packages/Arrangement/include/CGAL/Arrangement_2/Conic_arc_2_point.h @@ -0,0 +1,468 @@ +CGAL_BEGIN_NAMESPACE + +template +class Point_2_ex : public Point_2 +{ + public: + + typedef typename R::FT NT; + + enum Type + { + User_defined, + Tangency, + Intersection_exact, + Intersection_approx, + Ray_shooting_exact, + Ray_shooting_approx + }; + + private: + + Type _type; + int conic_id1; + int conic_id2; + APNT x_app; + APNT x_err; + int x_deg; + std::vector x_coeffs; + APNT y_app; + APNT y_err; + int y_deg; + std::vector y_coeffs; + + public: + + // Constructors. + Point_2_ex () : + Point_2(), + _type(User_defined), + conic_id1(0), + conic_id2(0), + x_app(0), + x_err(0), + x_deg(0), + y_app(0), + y_err(0), + y_deg(0) + {} + + Point_2_ex (const NT& hx, const NT& hy, const NT& hz, + const Type& type) : + Point_2(hx,hy,hz), + _type(type), + conic_id1(0), + conic_id2(0), + x_deg(0), + y_deg(0) + { + x_app = TO_APNT(x()); + x_err = 2 * x().get_double_error(); + y_app = TO_APNT(y()); + y_err = 2 * y().get_double_error(); + } + + Point_2_ex (const NT& hx, const NT& hy, + const Type& type, + const int& id1 = 0, const int& id2 = 0) : + Point_2(hx,hy), + _type(type), + conic_id1(id1), + conic_id2(id2), + x_deg(0), + y_deg(0) + { + x_app = TO_APNT(x()); + x_err = 2 * x().get_double_error(); + y_app = TO_APNT(y()); + y_err = 2 * y().get_double_error(); + } + + Point_2_ex (const NT& hx, const NT& hy) : + Point_2(hx,hy), + _type(User_defined), + conic_id1(0), + conic_id2(0), + x_deg(0), + y_deg(0) + { + x_app = TO_APNT(x()); + x_err = 2 * x().get_double_error(); + y_app = TO_APNT(y()); + y_err = 2 * y().get_double_error(); + } + + // Attach the generating polynomials information. + void attach_polynomials (const int& _x_deg, std::vector _x_coeffs, + const int& _y_deg, std::vector _y_coeffs) + { + // Copy the polynomials. + x_deg = _x_deg; + x_coeffs = _x_coeffs; + y_deg = _y_deg; + y_coeffs = _y_coeffs; + + // Adjust the degrees. + if (x_deg > 0) + { + while (x_deg > 0 && x_coeffs[x_deg] == 0) + x_deg--; + CGAL_assertion(x_deg > 0); + } + + if (y_deg > 0) + { + while (y_deg > 0 && y_coeffs[y_deg] == 0) + y_deg--; + CGAL_assertion(y_deg > 0); + } + + // Adjust the error bounds for roots of polynomials of degree > 2. + double slope; + NT val; + int k; + + if (x_deg > 2) + { + // Calculate p'(x). + slope = x_deg*TO_APNT(x_coeffs[x_deg]); + for (k = x_deg - 1; k >= 1 ; k--) + slope = slope*x_app + k*TO_APNT(x_coeffs[k]); + + slope = APNT_ABS(slope); + if (slope > 1) + slope = 1; + + if (slope > APNT_EPS) + { + // Calculate p(x_app). + val = x_coeffs[x_deg]; + for (k = x_deg - 1; k >= 0; k--) + val = val*x() + x_coeffs[k]; + + // Now since |x - x_app|*p'(x) ~= p(x_app) we get: + x_err = TO_APNT(x_deg*val / NT(slope)); + } + } + + if (y_deg > 2) + { + // Calculate p'(y). + slope = x_deg*TO_APNT(y_coeffs[y_deg]); + for (k = y_deg - 1; k >= 1 ; k--) + slope = slope*y_app + k*TO_APNT(y_coeffs[k]); + + slope = APNT_ABS(slope); + if (slope > 1) + slope = 1; + + if (slope > APNT_EPS) + { + // Calculate p(y_app). + val = y_coeffs[y_deg]; + for (k = y_deg - 1; k >= 0; k--) + val = val*y() + y_coeffs[k]; + + // Now since |y - y_app|*p'(y) ~= p(y_app) we get: + y_err = TO_APNT(y_deg*val / NT(slope)); + } + } + + return; + } + + // Get the approximated co-ordinates. + const APNT& x_approx () const + { + return (x_app); + } + + const APNT& y_approx () const + { + return (y_app); + } + + // Get the approximation error of the co-ordinates. + const APNT& x_approx_err () const + { + return (x_err); + } + + const APNT& y_approx_err () const + { + return (y_err); + } + + // Check if the point is approximate. + bool is_approximate () const + { + return (_type == Intersection_approx || + _type == Ray_shooting_approx); + } + + // Check if the given conic generates the points. + bool is_generating_conic_id (const int& id) const + { + return (id == conic_id1 || id == conic_id2); + } + + // Compare the co-ordinates of two given points. + Comparison_result compare_x (const Point_2_ex& p) const + { +#ifdef CGAL_CONIC_ARC_USE_FILTER + if (APNT_ABS(x_app - p.x_app) >= (x_err + p.x_err)) + { + // Using the approximated values is safe in this case! + if (x_app > p.x_app) + return (LARGER); + else if (x_app < p.x_app) + return (SMALLER); + } +#endif + + if ((is_approximate() || p.is_approximate()) && + eps_compare(TO_APNT(x()), TO_APNT(p.x())) == EQUAL) + { + if (x_deg > 0 && (p.x_deg == 0 || !p.is_approximate())) + { + if (_is_root (p.x(), x_coeffs, x_deg)) + return (EQUAL); + } + else if ((x_deg == 0 || ! is_approximate()) && p.x_deg > 0) + { + if (_is_root (x(), p.x_coeffs, p.x_deg)) + return (EQUAL); + } + else if (x_deg > 0 && p.x_deg > 0) + { + if ((conic_id1 == p.conic_id1 && conic_id2 == p.conic_id2) || + (conic_id1 == p.conic_id2 && conic_id2 == p.conic_id1)) + return (EQUAL); + + std::vector gcd_coeffs; + int deg_gcd; + + gcd_polynomials (x_coeffs, x_deg, + p.x_coeffs, p.x_deg, + gcd_coeffs, deg_gcd); + + if (deg_gcd == x_deg) + return (EQUAL); + + if (deg_gcd > 0 && _is_approx_root (x_app, x_err, + gcd_coeffs, deg_gcd)) + { + return (EQUAL); + } + + // The two roots are very close: perform exhaustive comparison. + return (compare_roots(x_coeffs, x_deg, + NT(x_app), NT(x_err), + p.x_coeffs, p.x_deg, + NT(p.x_app), NT(p.x_err))); + } + } + + return (CGAL::compare (x(), p.x())); + } + + Comparison_result compare_y (const Point_2_ex& p) const + { +#ifdef CGAL_CONIC_ARC_USE_FILTER + if (APNT_ABS(y_app - p.y_app) >= (y_err + p.y_err)) + { + // Using the approximated values is safe in this case! + if (y_app > p.y_app) + return (LARGER); + else if (y_app < p.y_app) + return (SMALLER); + } +#endif + + if ((is_approximate() || p.is_approximate()) && + eps_compare(TO_APNT(y()), TO_APNT(p.y())) == EQUAL) + { + if (y_deg > 0 && (p.y_deg == 0 || ! p.is_approximate())) + { + if (_is_root (p.y(), y_coeffs, y_deg)) + return (EQUAL); + } + else if ((y_deg == 0 || ! is_approximate()) && p.y_deg > 0) + { + if (_is_root (y(), p.y_coeffs, p.y_deg)) + return (EQUAL); + } + else if (y_deg > 0 && p.y_deg > 0) + { + if ((conic_id1 == p.conic_id1 && conic_id2 == p.conic_id2) || + (conic_id1 == p.conic_id2 && conic_id2 == p.conic_id1)) + return (EQUAL); + + std::vector gcd_coeffs; + int deg_gcd; + + gcd_polynomials (y_coeffs, y_deg, + p.y_coeffs, p.y_deg, + gcd_coeffs, deg_gcd); + + if (deg_gcd == y_deg) + return (EQUAL); + + if (deg_gcd > 0 && _is_approx_root (y_app, y_err, + gcd_coeffs, deg_gcd)) + { + return (EQUAL); + } + + // The two roots are very close: perform exhaustive comparison. + return (compare_roots(y_coeffs, y_deg, + NT(y_app), NT(y_err), + p.y_coeffs, p.y_deg, + NT(p.y_app), NT(p.y_err))); + } + } + + return (CGAL::compare (y(), p.y())); + } + + // Equality operators. + bool equals (const Point_2_ex& p) const + { + return (compare_x(p) == EQUAL && compare_y(p) == EQUAL); + } + + bool operator== (const Point_2_ex& p) const + { + return (compare_x(p) == EQUAL && compare_y(p) == EQUAL); + } + + bool operator!= (const Point_2_ex& p) const + { + return (compare_x(p) != EQUAL || compare_y(p) != EQUAL); + } + + // Copmare two points lexicographically. + Comparison_result compare_lex_xy (const Point_2_ex& p) const + { + Comparison_result res = this->compare_x (p); + + if (res != EQUAL) + return (res); + + return (this->compare_y (p)); + } + + Comparison_result compare_lex_yx (const Point_2_ex& p) const + { + Comparison_result res = this->compare_y (p); + + if (res != EQUAL) + return (res); + + return (this->compare_x (p)); + } + + // Reflect a point. + Point_2_ex reflect_in_y () const + { + Point_2_ex ref_point (-hx(), hy(), hw(), _type); + int i; + + if (conic_id1 != 0) + ref_point.conic_id1 = conic_id1 ^ REFLECT_IN_Y; + if (conic_id2 != 0) + ref_point.conic_id2 = conic_id2 ^ REFLECT_IN_Y; + + ref_point.x_app = -x_app; + ref_point.x_err = x_err; + ref_point.x_deg = x_deg; + if (x_deg > 0) + { + ref_point.x_coeffs.resize (x_deg + 1); + for (i = 0; i <= x_deg; i++) + ref_point.x_coeffs[i] = + (i % 2 == 0) ? x_coeffs[i] : -x_coeffs[i]; + } + + ref_point.y_app = y_app; + ref_point.y_err = y_err; + ref_point.y_deg = y_deg; + ref_point.y_coeffs = y_coeffs; + + return (ref_point); + } + + Point_2_ex reflect_in_x_and_y () const + { + Point_2_ex ref_point (-hx(), -hy(), hw(), _type); + int i; + + if (conic_id1 != 0) + ref_point.conic_id1 = conic_id1 ^ (REFLECT_IN_X | REFLECT_IN_Y); + if (conic_id2 != 0) + ref_point.conic_id2 = conic_id2 ^ (REFLECT_IN_X | REFLECT_IN_Y); + + ref_point.x_app = -x_app; + ref_point.x_err = x_err; + ref_point.x_deg = x_deg; + if (x_deg > 0) + { + ref_point.x_coeffs.resize(x_deg + 1); + for (i = 0; i <= x_deg; i++) + ref_point.x_coeffs[i] = + (i % 2 == 0) ? x_coeffs[i] : -x_coeffs[i]; + } + + ref_point.y_app = -y_app; + ref_point.y_err = y_err; + ref_point.y_deg = y_deg; + if (y_deg > 0) + { + ref_point.y_coeffs.resize(y_deg + 1); + for (i = 0; i <= y_deg; i++) + ref_point.y_coeffs[i] = + (i % 2 == 0) ? y_coeffs[i] : -y_coeffs[i]; + } + + return (ref_point); + } + +private: + + bool _is_root (const NT& z, + const std::vector& p, const int& deg) const + { + NT val = p[deg]; + int k; + + for (k = deg - 1; k >= 0 ; k--) + val = val*z + p[k]; + + return (CGAL::compare (val, NT(0)) == EQUAL); + } + + bool _is_approx_root (const APNT& z_app, const APNT& z_err, + const std::vector& p, const int& deg) const + { + APNT val; + APNT slope; + int k; + + // Compute p(z). + val = TO_APNT(p[deg]); + for (k = deg - 1; k >= 0 ; k--) + val = val*z_app + TO_APNT(p[k]); + + // Compute p'(z). + slope = deg*TO_APNT(p[deg]); + for (k = deg - 1; k >= 1 ; k--) + slope = slope*z_app + k*TO_APNT(p[k]); + slope = APNT_ABS(slope); + + // z approximates a root of p if p(z) < err(z)*p'(z). + return (APNT_ABS(val) < z_err*slope); + } + +}; + +CGAL_END_NAMESPACE diff --git a/Packages/Arrangement/include/CGAL/Arrangement_2/Polynom.h b/Packages/Arrangement/include/CGAL/Arrangement_2/Polynom.h new file mode 100644 index 00000000000..a2bcb8c1216 --- /dev/null +++ b/Packages/Arrangement/include/CGAL/Arrangement_2/Polynom.h @@ -0,0 +1,431 @@ +#ifndef __Polynom__ +#define __Polynom__ + +#include +#include + +template +class Polynom +{ + private: + + typename std::vector coeffs; // The vector of coefficients. + int degree; // The degree of the polynomial. + + public: + + // Constructors. + Polynom () : + coeffs(), + degree(-1) + {} + + // Destructor. + virtual ~Polynom () + {} + + // Get the degree of the polynomial (-1 means the polynomial equals 0). + int deg () const + { + return (degree); + } + + // Get the i'th coefficient of the polynomial. + const NT& get (const int& i) const + { + return (coeffs[i]); + } + + // Set the i'th coefficient of the polynomial. + void set (const int& i, const NT& c) + { + // Resize the coefficients vector, if needed. + if (i > degree) + { + coeffs.resize(i+1); + for (int j = degree + 1; j < i; j++) + coeffs[j] = NT(0); + degree = i; + } + + // Set the coefficient. + coeffs[i] = c; + + // If needed, update the degree. + _reduce_degree(); + + return; + } + + // Evaluate p(x0) (where p(x) is our polynomial). + NT eval_at (const NT& x0) const + { + // Check the case of a zero polynomial. + if (degree < 0) + return (NT(0)); + + // Use Horner's method: + // p(x) = a[n]*x^n + a[n-1]*x^(n-1) + ... + a[1]*x + a[0] = + // = (( ... (a[n]*x + a[n-1])*x + ... )*x + a[1])*x + a[0] + NT result = coeffs[degree]; + + for (int i = degree - 1; i >= 0; i--) + result = result*x0 + coeffs[i]; + + return (result); + } + + // Add two polynomials. + Polynom operator+ (const Polynom& p) const + { + // Find the polynomial with the higher degree. + const Polynom *hdeg_P, *ldeg_P; + + if (degree > p.degree) + { + hdeg_P = this; + ldeg_P = &p; + } + else + { + hdeg_P = &p; + ldeg_P = this; + } + + // Start with a copy of the polynomial of the higher degree. + Polynom res (*hdeg_P); + + // Add the coefficients of the other polynomial. + for (int j = 0; j <= ldeg_P->degree; j++) + { + res.coeffs[j] += ldeg_P->coeffs[j]; + } + + // If needed, reduce the degree of the resulting polynomial. + res._reduce_degree(); + + // Return the result. + return (res); + } + + // Subtract two polynomials. + Polynom operator- (const Polynom& p) const + { + // Find the polynomial with the higher degree. + const Polynom *hdeg_P, *ldeg_P; + int sign; + + if (degree > p.degree) + { + hdeg_P = this; + ldeg_P = &p; + sign = 1; + } + else + { + hdeg_P = &p; + ldeg_P = this; + sign = -1; + } + + // Start with a copy of the polynomial of the higher degree. + Polynom res (hdeg_P->degree); + int j; + + for (j = 0; j <= hdeg_P->degree; j++) + { + if (sign == 1) + res.coeffs[j] = hdeg_P->coeffs[j]; + else + res.coeffs[j] = -(hdeg_P->coeffs[j]); + } + + // Subtract the coefficients of the other polynomial. + for (int j = 0; j <= ldeg_P->degree; j++) + { + if (sign == 1) + res.coeffs[j] -= ldeg_P->coeffs[j]; + else + res.coeffs[j] += ldeg_P->coeffs[j]; + } + + // If needed, reduce the degree of the resulting polynomial. + res._reduce_degree(); + + // Return the result. + return (res); + } + + // Multiply two polynomials. + Polynom operator* (const Polynom& p) const + { + // In case one of the polynomial is zero, return a zero polynomial. + if (degree < 0 || p.degree < 0) + { + Polynom zero_poly; + return (zero_poly); + } + + // Allocate a polynomial whose degrees is the sum of the two degrees + // and initialize it with zeroes. + int sum_deg = degree + p.degree; + Polynom res (sum_deg); + int i, j; + const NT _zero (0); + + for (i = 0; i <= sum_deg; i++) + res.coeffs[i] = _zero; + + // Perform the convolution. + for (i = 0; i <= degree; i++) + for (j = 0; j <= p.degree; j++) + res.coeffs[i + j] += (coeffs[i] * p.coeffs[j]); + + // Return the result (no need for degree reduction). + return (res); + } + + // Divide two polynomials. + Polynom operator/ (const Polynom& p) const + { + Polynom q, r; + + // Find q(x), r(x) such that: + // (*this)(x) = q(x)*p(x) + r(x) + _divide (p, q, r); + + // Return the quontient polynomial. + return (q); + } + + // Return the residue from the division of the two polynomials. + Polynom operator% (const Polynom& p) const + { + Polynom q, r; + + // Find q(x), r(x) such that: + // (*this)(x) = q(x)*p(x) + r(x) + _divide (p, q, r); + + // Return the residue polynomial. + return (r); + } + + // Unary minus. + Polynom operator- () const + { + // Negate all coefficients. + Polynom res (degree); + + for (int i = 0; i <= degree; i++) + res.coeffs[i] = -(coeffs[i]); + + return (res); + } + + // Return the derivative of the polynomial. + Polynom derive () const + { + // In case of a constant (or a zero) polynomial, the derivative is zero. + if (degree <= 0) + { + Polynom zero_poly; + return (zero_poly); + } + + // If out polynomial is: + // p(x) = a[n]*x^n + a[n-1]*x^(n-1) + ... + a[1]*x + a[0] + // Then its derivative is: + // p(x) = n*a[n]*x^(n-1) + (n-1)*a[n-1]*x^(n-2) + ... + 2*a[2]*x + a[1] + Polynom res (degree - 1); + + for (int i = 1; i <= degree; i++) + res.coeffs[i - 1] = NT(i)*coeffs[i]; + + // Return the result (no need for degree reduction). + return (res); + } + + protected: + + // Protected constructor: set the initial size of the coefficients vector. + Polynom (const int& _degree) : + coeffs(_degree + 1), + degree(_degree) + {} + + // Reduce the degree of the polynomial (omit zero coefficients). + void _reduce_degree () + { + int old_deg = degree; + + // Omit any leading zero coefficients. + while (degree >= 0 && coeffs[degree] == NT(0)) + degree--; + + // If the degree has been reduced, resize the coefficients vector. + if (degree != old_deg) + coeffs.resize(degree + 1); + + return; + } + + // Calculate the quontient and residue polynomials of the division, + // i.e. find q(x) and r(x) such that: + // (*this)(x) = q(x)*p(x) + r(x) + void _divide (const Polynom& p, + Polynom& q, Polynom& r) const + { + // In case p(x) has a higher degree, then q(x) is zero and r(x) is (*this). + if (p.degree > degree) + { + Polynom zero_poly; + q = zero_poly; + r = *this; + return; + } + + // Otherwise, q(x)'s degree is the difference between the two degrees, and + // the final degree of r(x) should be less than p(x)'s degree. + const NT _zero (0); + int i; + + q.degree = degree - p.degree; + q.coeffs.resize (q.degree + 1); + + for (i = 0; i <= q.degree; i++) + q.coeffs[i] = _zero; + + // Perform "long division". + Polynom work1(*this), work2; + Polynom *work1_P = &work1, *work2_P = &work2, *swap_P; + int deg_diff; + NT factor; + int r_deg = degree; + + while (r_deg >= p.degree) + { + deg_diff = r_deg - p.degree; + factor = work1_P->coeffs[r_deg] / p.coeffs[p.degree]; + + q.coeffs[deg_diff] = factor; + + work2_P->coeffs.resize(r_deg); + for (i = p.degree - 1; i >= 0; i--) + work2_P->coeffs[i + deg_diff] = work1_P->coeffs[i + deg_diff] - + (factor * p.coeffs[i]); + for (i = 0; i < deg_diff; i++) + work2_P->coeffs[i] = work1_P->coeffs[i]; + + do + { + r_deg--; + } while (r_deg >= 0 && work2_P->coeffs[r_deg] == _zero); + + swap_P = work1_P; + work1_P = work2_P; + work2_P = swap_P; + } + + // Assign the residue polynomial. + if (r_deg < 0) + { + Polynom zero_residue; + r = zero_residue; + } + else + { + r.degree = r_deg; + r.coeffs.resize (r.degree + 1); + + for (i = 0; i <= r.degree; i++) + r.coeffs[i] = work1_P->coeffs[i]; + } + + return; + } + +}; + +// Print the polynomial. +template +std::ostream& operator<< (std::ostream& os, const Polynom& p) +{ + const NT _zero(0); + int degree = p.deg(); + + if (degree < 0) + { + os << '(' << _zero << ')'; + return (os); + } + + os << '(' << p.get(degree) << ')'; + if (degree == 1) + os << "*x"; + else if (degree > 1) + os << "*x^" << degree; + + for (int i = degree - 1; i >= 0; i--) + { + const NT& c = p.get(i); + + if (c == 0) + continue; + os << " + (" << c << ')'; + + if (i == 1) + os << "*x"; + else if (i > 1) + os << "*x^" << i; + } + + return (os); +} + +// Calculate the GCD of the two polynomials. +template +Polynom polynom_gcd (const Polynom& f, const Polynom& g) +{ + // Let p0(x) be the polynomial with the higher degree, and p1(x) + // be the polynomial with a lower degree. + Polynom p0, p1; + + if (f.deg() > g.deg()) + { + p0 = f; + p1 = g; + } + else + { + p0 = g; + p1 = f; + } + + // Use Euclid's algorithm. + Polynom r; + + do + { + r = p0 % p1; + + p0 = p1; + p1 = r; + } while (r.deg() > 0); + + if (r.deg() < 0) + { + // f(x), g(x) are not disjoint, and p0(x) is their GCD. + return (p0); + } + else + { + // f(x), g(x) are disjoint, therefore their GCD is 1. + Polynom unit_gcd; + unit_gcd.set (0, NT(1)); + + return (unit_gcd); + } +} + +#endif diff --git a/Packages/Arrangement/include/CGAL/Arrangement_2/Sturm_seq.h b/Packages/Arrangement/include/CGAL/Arrangement_2/Sturm_seq.h new file mode 100644 index 00000000000..ba458fe0c55 --- /dev/null +++ b/Packages/Arrangement/include/CGAL/Arrangement_2/Sturm_seq.h @@ -0,0 +1,160 @@ +#ifndef __STURM_SEQ__ +#define __STURM_SEQ__ + +#include +#include +#include + +template +class Sturm_seq +{ + private: + + std::list > fs; + Sturm_seq *next_P; + + // Copy constructor and assignment operator - not supported. + Sturm_seq (const Sturm_seq& ); + void operator= (const Sturm_seq& ); + + public: + + // Construct the Strum sequence for the given polynomial. + Sturm_seq (const Polynom& p) : + fs(), + next_P(NULL) + { + // If the polynomial is contant, just push it to the list and finish. + if (p.deg() <= 0) + { + fs.push_back(p); + return; + } + + // We shall start with f[0](x) = p(x), f[1](x) = p'(x). + Polynom f0 (p); + Polynom f1 = p.derive(); + + if (f1.deg() > 0) + { + // Calculate the greatest common divisor of p(x) and p'(x). + Polynom gcd = polynom_gcd (f0, f1); + + if (gcd.deg() >= 0) + { + // Because two polynomials p(x) and p'(x) are not disjoint, then p(x) + // has at least one root with multiplicity >= 1: + // First, construct a sequence for the GCD. + next_P = new Sturm_seq (gcd); + + // We know that q(x) = p(x)/GCD(p(x),p'(x)) has only simple roots, + // so we start the sequence with f[0](x) = q(x), f[1](x) = q'(x). + f0 = p / gcd; + f1 = f0.derive(); + } + } + + fs.push_back(f0); + fs.push_back(f1); + + // Complete the rest of the sequence. + const Polynom *fi_P = &f0; + const Polynom *fi1_P = &f1; + Polynom fi2; + + while (fi1_P->deg() > 0) + { + // The next polynomial f[i+2] is defined by -(f[i](x) % f[i+1](x)). + fs.push_back (-((*fi_P) % (*fi1_P))); + + fi_P = fi1_P; + fi1_P = &(fs.back()); + } + } + + // Destructor. + ~Sturm_seq () + { + if (next_P != NULL) + delete next_P; + } + + // Get the number of sign-changes in the Sturm sequence for a given x. + int sign_changes_at_x (const NT& x) const + { + // Count the sign changes in the sequence: + // f[0](x), f[1](x), ... , f[s](x). + typename std::list >::const_iterator iter; + NT y; + const NT _zero (0); + int prev_sign = 0; + int curr_sign; + int result = 0; + + for (iter = fs.begin(); iter != fs.end(); ++iter) + { + // Determine the sign of f[i](x). + y = (*iter).eval_at(x); + + if (y == _zero) + curr_sign = 0; + else + curr_sign = (y > _zero) ? 1 : -1; + + + // Count if there was a sign change. + if (prev_sign != 0 && curr_sign != prev_sign) + result++; + + prev_sign = curr_sign; + } + + // Sum up the result with the number of sign changes in the next sequence. + if (next_P != NULL) + return (result + next_P->sign_changes_at_x (x)); + else + return (result); + } + + // Get the number of sign-changes at Infinity (if inf > 0) or at -Infinity + // (if inf < 0). + int sign_changes_at_infinity (const int& inf) const + { + // Count the sign changes in the sequence: + // f[0](-Inf), f[1](-Inf), ... , f[s](-Inf). + typename std::list >::const_iterator iter; + int deg; + NT leading_coeff; + const NT _zero(0); + int prev_sign = 0; + int curr_sign; + int result = 0; + + for (iter = fs.begin(); iter != fs.end(); ++iter) + { + // Determine the sign of f[i](-Infinity): This is the sign of the leading + // coefficient if the polynomial has an even dergee, otherwise it has an + // opposite sign. + deg = (*iter).deg(); + leading_coeff = (*iter).get(deg); + + curr_sign = (leading_coeff < _zero) ? -1 : 1; + if (inf < 0 && deg % 2 == 1) + curr_sign = -curr_sign; + + // Count if there was a sign change. + if (prev_sign != 0 && curr_sign != prev_sign) + result++; + + prev_sign = curr_sign; + } + + // Sum up the result with the number of sign changes in the next sequence. + if (next_P != NULL) + return (result + next_P->sign_changes_at_infinity(inf)); + else + return (result); + } +}; + +#endif