diff --git a/Arrangement_on_surface_2/include/CGAL/Arr_conic_traits_2.h b/Arrangement_on_surface_2/include/CGAL/Arr_conic_traits_2.h index 54a3db3fe01..e794c25c97f 100644 --- a/Arrangement_on_surface_2/include/CGAL/Arr_conic_traits_2.h +++ b/Arrangement_on_surface_2/include/CGAL/Arr_conic_traits_2.h @@ -64,6 +64,7 @@ public: typedef typename Rat_kernel::Circle_2 Rat_circle_2; typedef typename Alg_kernel::FT Algebraic; + typedef typename Alg_kernel::Point_2 Alg_point_2; typedef typename Nt_traits::Integer Integer; @@ -88,19 +89,35 @@ public: private: // Type definition for the intersection points mapping. - typedef typename X_monotone_curve_2::Conic_id Conic_id; - typedef typename X_monotone_curve_2::Intersection_point Intersection_point; - typedef typename X_monotone_curve_2::Intersection_map Intersection_map; + typedef typename Point_2::Conic_id Conic_id; + typedef std::pair Conic_pair; - typedef std::shared_ptr Shared_rat_kernel; - typedef std::shared_ptr Shared_alg_kernel; - typedef std::shared_ptr Shared_nt_Traits; + /*! \struct Less functor for Conic_pair. + */ + struct Less_conic_pair { + bool operator()(const Conic_pair& cp1, const Conic_pair& cp2) const { + // Compare the pairs of IDs lexicographically. + return ((cp1.first < cp2.first) || + ((cp1.first == cp2.first) && (cp1.second < cp2.second))); + } + }; + + typedef std::pair Intersection_point; + typedef std::list Intersection_list; + typedef std::map + Intersection_map; + typedef typename Intersection_map::iterator Intersection_map_iterator; + + + typedef std::shared_ptr Shared_rat_kernel; + typedef std::shared_ptr Shared_alg_kernel; + typedef std::shared_ptr Shared_nt_Traits; const Shared_rat_kernel m_rat_kernel; const Shared_alg_kernel m_alg_kernel; const Shared_nt_Traits m_nt_traits; - mutable Intersection_map inter_map; // Mapping conic pairs to their + mutable Intersection_map m_inter_map; // Mapping conic pairs to their // intersection points. public: @@ -253,36 +270,35 @@ public: * LARGER if y(p) > cv(x(p)), i.e. the point is above the curve; * EQUAL if p lies on the curve. */ - Comparison_result operator()(const Point_2& p, - const X_monotone_curve_2& cv) const - { + Comparison_result operator()(const Point_2& p, const X_monotone_curve_2& xcv) + const { auto cmp_y = m_traits.m_alg_kernel->compare_y_2_object(); - if (cv.is_vertical()) { + if (xcv.is_vertical()) { // A special treatment for vertical segments: // In case p has the same x c-ordinate of the vertical segment, compare // it to the segment endpoints to determine its position. - Comparison_result res1 = cmp_y(p, cv.left()); - Comparison_result res2 = cmp_y(p, cv.right()); + Comparison_result res1 = cmp_y(p, xcv.left()); + Comparison_result res2 = cmp_y(p, xcv.right()); return (res1 == res2) ? res1 : EQUAL; } // Check whether the point is exactly on the curve. - if (cv.contains_point(p)) return EQUAL; + if (m_traits.contains_point(xcv, p)) return EQUAL; // Obtain a point q on the x-monotone arc with the same x coordinate as p. Point_2 q; auto cmp_x = m_traits.m_alg_kernel->compare_x_2_object(); - Comparison_result x_res_left = cmp_x(p, cv.left()); - if (x_res_left == EQUAL) q = cv.left(); + Comparison_result x_res_left = cmp_x(p, xcv.left()); + if (x_res_left == EQUAL) q = xcv.left(); else { CGAL_precondition(x_res_left != SMALLER); - auto x_res_right = cmp_x(p, cv.right()); - if (x_res_right == EQUAL) q = cv.right(); + auto x_res_right = cmp_x(p, xcv.right()); + if (x_res_right == EQUAL) q = xcv.right(); else { CGAL_precondition(x_res_right != LARGER); - q = point_at_x(cv, p); + q = point_at_x(xcv, p); } } @@ -299,9 +315,10 @@ public: */ Point_2 point_at_x(const X_monotone_curve_2& xcv, const Point_2& p) const { // Make sure that p is in the x-range of the arc. - CGAL_precondition(! xcv.test_flag(X_monotone_curve_2::IS_VERTICAL_SEGMENT)); + CGAL_precondition(! xcv.is_vertical()); - CGAL_precondition_code(auto cmp_x = m_traits.m_alg_kernel->compare_x_2_object()); + CGAL_precondition_code(auto cmp_x = + m_traits.m_alg_kernel->compare_x_2_object()); CGAL_precondition((cmp_x(p, xcv.left()) != SMALLER) && (cmp_x(p, xcv.right()) != LARGER)); @@ -320,7 +337,7 @@ public: // conic curve. Algebraic y; - if (xcv.degree_mask() == xcv.flag_mask(X_monotone_curve_2::DEGREE_1)) { + if (xcv.degree_mask() == X_monotone_curve_2::degree_1_mask()) { // In case of a linear curve, the y-coordinate is a simple linear // expression of x(p) (note that v is not 0 as the arc is not vertical): // y = -(u*x(p) + w) / v @@ -335,7 +352,7 @@ public: y = -(extra_data->a * p.x() + extra_data->c) / extra_data->b; } else { - CGAL_assertion(xcv.degree_mask() == xcv.flag_mask(X_monotone_curve_2::DEGREE_2)); + CGAL_assertion(xcv.degree_mask() == X_monotone_curve_2::degree_2_mask()); // In this case the y-coordinate is one of solutions to the quadratic // equation: @@ -406,18 +423,15 @@ public: { // Make sure that p lies on both curves, and that both are defined to its // left (so their left endpoint is lexicographically smaller than p). - CGAL_precondition(xcv1.contains_point(p) && - xcv2.contains_point(p)); + CGAL_precondition(m_traits.contains_point(xcv1, p) && + m_traits.contains_point(xcv2, p)); CGAL_precondition_code(const auto ker = m_traits.m_alg_kernel); CGAL_precondition(ker->compare_xy_2_object()(p, xcv1.left()) == LARGER && ker->compare_xy_2_object()(p, xcv2.left()) == LARGER); // If one of the curves is vertical, it is below the other one. - if (xcv1.is_vertical()) { - // Check whether both are vertical: - return (xcv2.is_vertical()) ? EQUAL : SMALLER; - } + if (xcv1.is_vertical()) return (xcv2.is_vertical()) ? EQUAL : SMALLER; else if (xcv2.is_vertical()) return LARGER; // Compare the two curves immediately to the left of p: @@ -433,11 +447,8 @@ public: */ Comparison_result compare_to_left(const X_monotone_curve_2& xcv1, const X_monotone_curve_2& xcv2, - const Point_2& p) - const - { - CGAL_precondition(! xcv1.test_flag(X_monotone_curve_2::IS_VERTICAL_SEGMENT) && - ! xcv2.test_flag(X_monotone_curve_2::IS_VERTICAL_SEGMENT)); + const Point_2& p) const { + CGAL_precondition(! xcv1.is_vertical() && ! xcv2.is_vertical()); // In case one arc is facing upwards and another facing downwards, it is // clear that the one facing upward is above the one facing downwards. @@ -587,16 +598,15 @@ public: { // Make sure that p lies on both curves, and that both are defined to its // left (so their left endpoint is lexicographically smaller than p). - CGAL_precondition(xcv1.contains_point(p) && xcv2.contains_point(p)); + CGAL_precondition(m_traits.contains_point(xcv1, p) && + m_traits.contains_point(xcv2, p)); CGAL_precondition_code(const auto ker = m_traits.m_alg_kernel); CGAL_precondition_code(auto cmp_xy = ker->compare_xy_2_object()); CGAL_precondition(cmp_xy(p, xcv1.right()) == SMALLER && cmp_xy(p, xcv2.right()) == SMALLER); // If one of the curves is vertical, it is above the other one. - if (xcv1.is_vertical()) - // Check whether both are vertical: - return (xcv2.is_vertical()) ? EQUAL : LARGER; + if (xcv1.is_vertical()) return (xcv2.is_vertical()) ? EQUAL : LARGER; else if (xcv2.is_vertical()) return SMALLER; // Compare the two curves immediately to the right of p: @@ -613,8 +623,7 @@ public: Comparison_result compare_to_right(const X_monotone_curve_2& xcv1, const X_monotone_curve_2& xcv2, const Point_2& p) const { - CGAL_precondition(! xcv1.test_flag(X_monotone_curve_2::IS_VERTICAL_SEGMENT) && - ! xcv2.test_flag(X_monotone_curve_2::IS_VERTICAL_SEGMENT)); + CGAL_precondition(! xcv1.is_vertical() && ! xcv2.is_vertical()); // In case one arc is facing upwards and another facing downwards, it is // clear that the one facing upward is above the one facing downwards. @@ -752,7 +761,7 @@ public: const X_monotone_curve_2& xcv2) const { if (&xcv1 == &xcv2) return true; - return xcv1.equals(xcv2); + return equals(xcv1, xcv2); } /*! Check whether the two points are the same. @@ -764,6 +773,38 @@ public: if (&p1 == &p2) return (true); return(m_traits.m_alg_kernel->compare_xy_2_object()(p1, p2) == EQUAL); } + + private: + /*! Check whether the two arcs are equal (have the same graph). + * \param arc The compared arc. + * \return (true) if the two arcs have the same graph; (false) otherwise. + */ + bool equals(const X_monotone_curve_2& xcv1, + const X_monotone_curve_2& xcv2) const { + // The two arc must have the same supporting conic curves. + if (! m_traits.has_same_supporting_conic(xcv1, xcv2)) return false; + + auto eq = m_traits.m_alg_kernel->equal_2_object(); + + // Check that the arc endpoints are the same. + if (xcv1.orientation() == COLLINEAR) { + CGAL_assertion(xcv2.orientation() == COLLINEAR); + return((eq(xcv1.source(), xcv2.source()) && + eq(xcv1.target(), xcv2.target())) || + (eq(xcv1.source(), xcv2.target()) && + eq(xcv1.target(), xcv2.source()))); + } + + if (xcv1.orientation() == xcv2.m_orient) { + // Same orientation - the source and target points must be the same. + return (eq(xcv1.source(), xcv2.source()) && + eq(xcv1.target(), xcv2.target())); + } + + // Reverse orientation - the source and target points must be swapped. + return (eq(xcv1.source(), xcv2.target()) && + eq(xcv1.target(), xcv2.source())); + } }; /*! Obtain an Equal_2 functor object. */ @@ -812,8 +853,8 @@ public: Conic_id conic_id(index); // Find the points of vertical tangency to cv and act accordingly. - typename Curve_2::Point_2 vtan_ps[2]; - auto n_vtan_ps = cv.vertical_tangency_points(vtan_ps); + Alg_point_2 vtan_ps[2]; + auto n_vtan_ps = m_traits.vertical_tangency_points(cv, vtan_ps); if (n_vtan_ps == 0) { // In case the given curve is already x-monotone: *oi++ = Make_x_monotone_result(ctr_xcv(cv, conic_id)); @@ -828,19 +869,19 @@ public: // In case the curve is a full conic, split it into two x-monotone // arcs, one going from ps[0] to ps[1], and the other from ps[1] to // ps[0]. - *oi++ = Make_x_monotone_result(ctr_xcv(cv, vtan_ps[0], - vtan_ps[1], conic_id)); - *oi++ = Make_x_monotone_result(ctr_xcv(cv, vtan_ps[1], - vtan_ps[0], conic_id)); + *oi++ = Make_x_monotone_result(ctr_xcv(cv, vtan_ps[0], vtan_ps[1], + conic_id)); + *oi++ = Make_x_monotone_result(ctr_xcv(cv, vtan_ps[1], vtan_ps[0], + conic_id)); } else { if (n_vtan_ps == 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. - *oi++ = Make_x_monotone_result(ctr_xcv(cv, cv.source(), - vtan_ps[0], conic_id)); - *oi++ = Make_x_monotone_result(ctr_xcv(cv, vtan_ps[0], - cv.target(), conic_id)); + *oi++ = Make_x_monotone_result(ctr_xcv(cv, cv.source(), vtan_ps[0], + conic_id)); + *oi++ = Make_x_monotone_result(ctr_xcv(cv, vtan_ps[0], cv.target(), + conic_id)); } else { CGAL_assertion(n_vtan_ps == 2); @@ -928,7 +969,7 @@ public: X_monotone_curve_2& xcv1, X_monotone_curve_2& xcv2) const { // Make sure that p lies on the interior of the arc. CGAL_precondition_code(auto eq = m_traits.m_alg_kernel->equal_2_object()); - CGAL_precondition(xcv.contains_point(p) && + CGAL_precondition(m_traits.contains_point(xcv, p) && ! eq(p, xcv.source()) && ! eq(p, xcv.target())); // Make copies of the current arc. @@ -966,13 +1007,25 @@ public: Split_2 split_2_object() const { return Split_2(*this); } class Intersect_2 { - private: + protected: + using Traits = Arr_conic_traits_2; + Intersection_map& m_inter_map; // The map of intersection points. - public: - /*! Constructor. */ - Intersect_2(Intersection_map& map) : m_inter_map(map) {} + /*! The traits (in case it has state) */ + const Traits& m_traits; + /*! Constructor. + * \param traits the traits. + */ + Intersect_2(Intersection_map& map, const Traits& traits) : + m_inter_map(map), + m_traits(traits) + {} + + friend class Arr_conic_traits_2; + + public: /*! Find the intersections of the two given curves and insert them to the * given output iterator. As two segments may itersect only once, only a * single will be contained in the iterator. @@ -982,14 +1035,351 @@ public: * \return The past-the-end iterator. */ template - OutputIterator operator()(const X_monotone_curve_2& cv1, - const X_monotone_curve_2& cv2, + OutputIterator operator()(const X_monotone_curve_2& xcv1, + const X_monotone_curve_2& xcv2, OutputIterator oi) const - { return cv1.intersect(cv2, m_inter_map, oi); } + { return intersect(xcv1, xcv2, m_inter_map, oi); } + + private: + /*! Compute the overlap with a given arc, which is supposed to have the same + * supporting conic curve as this arc. + * \param arc The given arc. + * \param overlap Output: The overlapping arc (if any). + * \return Whether we found an overlap. + */ + bool compute_overlap(const X_monotone_curve_2& xcv1, + const X_monotone_curve_2& xcv2, + X_monotone_curve_2& overlap) const { + // Check if the two arcs are identical. + if (m_traits.equal_2_object()(xcv1,xcv2)) { + overlap = xcv2; + return true; + } + + if (m_traits.is_strictly_between_endpoints(xcv1, xcv2.left())) { + if (m_traits.is_strictly_between_endpoints(xcv1, xcv2.right())) { + // Case 1 - *this: +-----------> + // arc: +=====> + overlap = xcv2; + return true; + } + else { + // Case 2 - *this: +-----------> + // arc: +=====> + overlap = xcv1; + + if (overlap.is_directed_right()) overlap.m_source = xcv2.left(); + else overlap.m_target = xcv2.left(); + + return true; + } + } + else if (m_traits.is_strictly_between_endpoints(xcv1, xcv2.right())) { + // Case 3 - *this: +-----------> + // arc: +=====> + overlap = xcv1; + if (overlap.is_directed_right()) overlap.m_target = xcv2.right(); + else overlap.m_source = xcv2.right(); + return true; + } + else if (m_traits.is_between_endpoints(xcv2, xcv1.source()) && + m_traits.is_between_endpoints(xcv2, xcv1.target()) && + (m_traits.is_strictly_between_endpoints(xcv2, xcv1.source()) || + m_traits.is_strictly_between_endpoints(xcv2, xcv1.target()))) + { + // Case 4 - *this: +-----------> + // arc: +================> + overlap = xcv1; + return true; + } + + // If we reached here, there are no overlaps: + return false; + } + + /*! Intersect the supporing conic curves of this arc and the given arc. + * \param arc The arc to intersect with. + * \param inter_list The list of intersection points. + */ + void intersect_supporting_conics(const X_monotone_curve_2& xcv1, + const X_monotone_curve_2& xcv2, + Intersection_list& inter_list) const { + if (xcv1.is_special_segment() && ! xcv2.is_special_segment()) { + // If one of the arcs is a special segment, make sure it is (arc). + intersect_supporting_conics(xcv2, xcv1, inter_list); + return; + } + + const int deg1 = (xcv1.degree_mask() == X_monotone_curve_2::degree_1_mask()) ? 1 : 2; + const int deg2 = (xcv2.degree_mask() == X_monotone_curve_2::degree_1_mask()) ? 1 : 2; + Nt_traits nt_traits; + Algebraic xs[4]; + int n_xs = 0; + Algebraic ys[4]; + int n_ys = 0; + + if (xcv2.is_special_segment()) { + // The second arc is a special segment (a*x + b*y + c = 0). + if (xcv1.is_special_segment()) { + // Both arc are sepcial segment, so they have at most one intersection + // point. + const auto* extra_data1 = xcv1.extra_data(); + const auto* extra_data2 = xcv2.extra_data(); + Algebraic denom = + extra_data1->a * extra_data2->b - extra_data1->b * extra_data2->a; + + if (CGAL::sign (denom) != CGAL::ZERO) { + xs[0] = (extra_data1->b * extra_data2->c - + extra_data1->c * extra_data2->b) / denom; + n_xs = 1; + + ys[0] = (extra_data1->c * extra_data2->a - + extra_data1->a * extra_data2->c) / denom; + n_ys = 1; + } + } + else { + const auto* extra_data2 = xcv2.extra_data(); + + // Compute the x-coordinates of the intersection points. + n_xs = compute_resultant_roots(nt_traits, + xcv1.alg_r(), xcv1.alg_s(), + xcv1.alg_t(), xcv1.alg_u(), + xcv1.alg_v(), xcv1.alg_w(), + deg1, + extra_data2->a, + extra_data2->b, + extra_data2->c, + xs); + CGAL_assertion(n_xs <= 2); + + // Compute the y-coordinates of the intersection points. + n_ys = compute_resultant_roots(nt_traits, + xcv1.alg_s(), xcv1.alg_r(), + xcv1.alg_t(), xcv1.alg_v(), + xcv1.alg_u(), xcv1.alg_w(), + deg1, + extra_data2->b, + extra_data2->a, + extra_data2->c, + ys); + CGAL_assertion(n_ys <= 2); + } + } + else { + // Compute the x-coordinates of the intersection points. + n_xs = compute_resultant_roots(nt_traits, + xcv1.r(), xcv1.s(), xcv1.t(), + xcv1.u(), xcv1.v(), xcv1.w(), + deg1, + xcv2.r(), xcv2.s(), xcv2.t(), + xcv2.u(), xcv2.v(), xcv2.w(), + deg2, + xs); + CGAL_assertion(n_xs <= 4); + + // Compute the y-coordinates of the intersection points. + n_ys = compute_resultant_roots(nt_traits, + xcv1.s(), xcv1.r(), xcv1.t(), + xcv1.v(), xcv1.u(), xcv1.w(), + deg1, + xcv2.s(), xcv2.r(), xcv2.t(), + xcv2.v(), xcv2.u(), xcv2.w(), + deg2, + ys); + CGAL_assertion(n_ys <= 4); + } + + // Pair the coordinates of the intersection points. As the vectors of + // x and y-coordinates are sorted in ascending order, we output the + // intersection points in lexicographically ascending order. + unsigned int mult; + int i, j; + + if (xcv2.is_special_segment()) { + if ((n_xs == 0) || (n_ys == 0)) return; + + if ((n_xs == 1) && (n_ys == 1)) { + // Single intersection. + Point_2 ip(xs[0], ys[0]); + ip.set_generating_conic(xcv1.id()); + ip.set_generating_conic(xcv2.id()); + + // In case the other curve is of degree 2, this is a tangency point. + mult = ((deg1 == 1) || xcv1.is_special_segment()) ? 1 : 2; + inter_list.push_back(Intersection_point(ip, mult)); + } + else if ((n_xs == 1) && (n_ys == 2)) { + Point_2 ip1(xs[0], ys[0]); + ip1.set_generating_conic(xcv1.id()); + ip1.set_generating_conic(xcv2.id()); + + inter_list.push_back(Intersection_point(ip1, 1)); + + Point_2 ip2(xs[0], ys[1]); + ip2.set_generating_conic(xcv1.id()); + ip2.set_generating_conic(xcv2.id()); + + inter_list.push_back(Intersection_point(ip2, 1)); + } + else if ((n_xs == 2) && (n_ys == 1)) { + Point_2 ip1(xs[0], ys[0]); + ip1.set_generating_conic(xcv1.id()); + ip1.set_generating_conic(xcv2.id()); + + inter_list.push_back(Intersection_point(ip1, 1)); + + Point_2 ip2(xs[1], ys[0]); + ip2.set_generating_conic(xcv1.id()); + ip2.set_generating_conic(xcv2.id()); + + inter_list.push_back(Intersection_point(ip2, 1)); + } + else { + CGAL_assertion((n_xs == 2) && (n_ys == 2)); + + // The x-coordinates and the y-coordinates are given in ascending + // order. If the slope of the segment is positive, we pair the + // coordinates as is - otherwise, we swap the pairs. + int ind_first_y(0), ind_second_y(1); + + const auto* extra_data2 = xcv2.extra_data(); + if (CGAL::sign(extra_data2->b) == CGAL::sign(extra_data2->a)) { + ind_first_y = 1; + ind_second_y = 0; + } + + Point_2 ip1(xs[0], ys[ind_first_y]); + ip1.set_generating_conic(xcv1.id()); + ip1.set_generating_conic(xcv2.id()); + + inter_list.push_back(Intersection_point(ip1, 1)); + + Point_2 ip2(xs[1], ys[ind_second_y]); + ip2.set_generating_conic(xcv1.id()); + ip2.set_generating_conic(xcv2.id()); + + inter_list.push_back(Intersection_point(ip2, 1)); + } + + return; + } + + for (i = 0; i < n_xs; ++i) { + for (j = 0; j < n_ys; ++j) { + if (xcv1.is_on_supporting_conic(xs[i], ys[j]) && + xcv2.is_on_supporting_conic(xs[i], ys[j])) + { + // Create the intersection point and set its generating conics. + Point_2 ip(xs[i], ys[j]); + + ip.set_generating_conic(xcv1.id()); + ip.set_generating_conic(xcv2.id()); + + // Compute the multiplicity of the intersection point. + if (deg1 == 1 && deg2 == 1) mult = 1; + else mult = xcv1.multiplicity_of_intersection_point(xcv2, ip); + + // Insert the intersection point to the output list. + inter_list.push_back(Intersection_point(ip, mult)); + } + } + } + } + + /*! Compute the intersections with the given arc. + * \param arc The given intersecting arc. + * \param inter_map Maps conic pairs to lists of their intersection points. + * \param oi The output iterator. + * \return The past-the-end iterator. + */ + template + OutputIterator intersect(const X_monotone_curve_2& xcv1, + const X_monotone_curve_2& xcv2, + Intersection_map& inter_map, + OutputIterator oi) const { + typedef boost::variant + Intersection_result; + + if (m_traits.has_same_supporting_conic(xcv1, xcv2)) { + // Check for overlaps between the two arcs. + X_monotone_curve_2 overlap; + + if (compute_overlap(xcv1, xcv2, overlap)) { + // There can be just a single overlap between two x-monotone arcs: + *oi++ = Intersection_result(overlap); + return oi; + } + + // In case there is not overlap and the supporting conics are the same, + // there cannot be any intersection points, unless the two arcs share + // an end point. + // Note that in this case we do not define the multiplicity of the + // intersection points we report. + auto alg_kernel = m_traits.m_alg_kernel; + auto eq = alg_kernel->equal_2_object(); + if (eq(xcv1.left(), xcv2.left())) { + Intersection_point ip(xcv1.left(), 0); + *oi++ = Intersection_result(ip); + } + + if (eq(xcv1.right(), xcv2.right())) { + Intersection_point ip(xcv1.right(), 0); + *oi++ = Intersection_result(ip); + } + + return oi; + } + + // Search for the pair of supporting conics in the map (the first conic + // ID in the pair should be smaller than the second one, to guarantee + // uniqueness). + Conic_pair conic_pair; + Intersection_map_iterator map_iter; + Intersection_list inter_list; + bool invalid_ids = false; + + if (xcv1.id().is_valid() && xcv2.id().is_valid()) { + if (xcv1.id() < xcv2.id()) conic_pair = Conic_pair(xcv1.id(), xcv2.id()); + else conic_pair = Conic_pair(xcv2.id(), xcv1.id()); + map_iter = inter_map.find(conic_pair); + } + else { + // In case one of the IDs is invalid, we do not look in the map neither + // we cache the results. + map_iter = inter_map.end(); + invalid_ids = true; + } + + if (map_iter == inter_map.end()) { + // In case the intersection points between the supporting conics have + // not been computed before, compute them now and store them in the map. + intersect_supporting_conics(xcv1, xcv2, inter_list); + + if (! invalid_ids) inter_map[conic_pair] = inter_list; + } + else { + // Obtain the precomputed intersection points from the map. + inter_list = (*map_iter).second; + } + + // Go over the list of intersection points and report those that lie on + // both x-monotone arcs. + for (auto iter = inter_list.begin(); iter != inter_list.end(); ++iter) { + if (m_traits.is_between_endpoints(xcv1, (*iter).first) && + m_traits.is_between_endpoints(xcv2, (*iter).first)) + { + *oi++ = Intersection_result(*iter); + } + } + + return oi; + } }; /*! Obtain an Intersect_2 functor object. */ - Intersect_2 intersect_2_object() const { return (Intersect_2(inter_map)); } + Intersect_2 intersect_2_object() const + { return Intersect_2(m_inter_map, *this); } class Are_mergeable_2 { protected: @@ -1432,8 +1822,8 @@ public: * \pre cv is x-monotone. */ X_monotone_curve_2 operator()(const Curve_2& cv) const { + CGAL_precondition(cv.is_valid() && is_x_monotone(cv)); X_monotone_curve_2 xcv(cv); - CGAL_precondition(xcv.is_valid() && xcv.is_x_monotone()); m_traits.set_x_monotone(xcv); return xcv; } @@ -1442,8 +1832,7 @@ public: * \param xcv The given (Curve_2) curve. * \param id The ID of the base curve. */ - X_monotone_curve_2 operator()(const Curve_2& cv, const Conic_id& id) const - { + X_monotone_curve_2 operator()(const Curve_2& cv, const Conic_id& id) const { X_monotone_curve_2 xcv(cv, id); CGAL_precondition(xcv.is_valid() && id.is_valid()); m_traits.set_x_monotone(xcv); @@ -1554,6 +1943,26 @@ public: xcv.set_flag(X_monotone_curve_2::IS_SPECIAL_SEGMENT); return xcv; } + + private: + /*! Determine whether the arc is x-monotone. + */ + bool is_x_monotone(const Curve_2& cv) const { + // Collect vertical tangency points. + Alg_point_2 vtan_ps[2]; + auto res = m_traits.vertical_tangency_points(cv, vtan_ps); + return (res == 0); + } + + /*! Determine whether the arc is y-monotone. + */ + bool is_y_monotone(const Curve_2& cv) const { + // Collect horizontal tangency points. + Alg_point_2 htan_ps[2]; + auto res = m_traits.horizontal_tangency_points(cv, htan_ps); + return (res == 0); + } + }; /*! Obtain a Construct_x_monotone_curve_2 functor object. */ @@ -1723,8 +2132,7 @@ public: Curve_2 arc; // Make sure that no three points are collinear. - Rat_kernel ker; - auto orient_f = m_rat_kernel.orientation_2_object(); + auto orient_f = m_traits.m_rat_kernel->orientation_2_object(); const bool point_collinear = (orient_f(p1, p2, p3) == COLLINEAR || orient_f(p1, p2, p4) == COLLINEAR || @@ -1795,9 +2203,9 @@ public: Point_2 mp4 = Point_2(nt_traits->convert(p4.x()), nt_traits->convert(p4.y())); - if (! is_strictly_between_endpoints(arc, mp2) || - ! is_strictly_between_endpoints(arc, mp3) || - ! is_strictly_between_endpoints(arc, mp4)) + if (! m_traits.is_strictly_between_endpoints(arc, mp2) || + ! m_traits.is_strictly_between_endpoints(arc, mp3) || + ! m_traits.is_strictly_between_endpoints(arc, mp4)) { arc.reset_flags(); // inavlid arc return arc; @@ -2101,8 +2509,8 @@ public: Curve_2 operator()(const Rat_circle_2& circ, const Orientation& orient, const Point_2& source, const Point_2& target) const { // Make sure that the source and the taget are not the same. - CGAL_precondition(Alg_kernel().compare_xy_2_object()(source, target) != - EQUAL); + CGAL_precondition_code(auto cmp_xy = m_traits.m_alg_kernel->compare_xy_2_object()); + CGAL_precondition(cmp_xy(source, target) != EQUAL); CGAL_precondition(orient != COLLINEAR); Curve_2 arc; @@ -2150,9 +2558,8 @@ public: // Set the arc properties (no need to compute the orientation). m_traits.set(arc, rat_coeffs); - return arc; - } + } }; /*! Obtain a Construct_curve_2 functor object. */ @@ -2222,8 +2629,7 @@ public: * \pre both points must be interior and must lie on \c cv */ X_monotone_curve_2 operator()(const X_monotone_curve_2& xcv, - const Point_2& src, - const Point_2& tgt)const + const Point_2& src, const Point_2& tgt) const { // make functor objects CGAL_precondition_code(Compare_y_at_x_2 compare_y_at_x_2 = @@ -2255,7 +2661,8 @@ public: const Point_2& ps, const Point_2& pt) const { // Make sure that both ps and pt lie on the arc. - CGAL_precondition(xcv.contains_point(ps) && xcv.contains_point(pt)); + CGAL_precondition(m_traits.contains_point(xcv, ps) && + m_traits.contains_point(xcv, pt)); // Make sure that the endpoints conform with the direction of the arc. X_monotone_curve_2 res_xcv = xcv; @@ -2498,13 +2905,13 @@ public: * \return true if the point is between the two endpoints, * (false) if it is not. */ - bool is_between_endpoints(const Curve_2& arc, const Point_2& p) const { - CGAL_precondition(! arc.is_full_conic()); + bool is_between_endpoints(const Curve_2& cv, const Point_2& p) const { + CGAL_precondition(! cv.is_full_conic()); // Check if p is one of the endpoints. auto eq = m_alg_kernel->equal_2_object(); - if (eq(p, arc.source()) || eq(p, arc.target())) return true; - else return (is_strictly_between_endpoints(arc, p)); + if (eq(p, cv.source()) || eq(p, cv.target())) return true; + else return is_strictly_between_endpoints(cv, p); } /*! Check whether the given point is strictly between the source and the @@ -2514,33 +2921,33 @@ public: * \return true if the point is strictly between the two endpoints, * (false) if it is not. */ - bool is_strictly_between_endpoints(const Curve_2& arc, const Point_2& p) const + bool is_strictly_between_endpoints(const Curve_2& cv, const Point_2& p) const { // In case this is a full conic, any point on its boundary is between // its end points. - if (arc.is_full_conic()) return true; + if (cv.is_full_conic()) return true; // Check if we have extra data available. - const auto* extra_data = arc.extra_data(); + const auto* extra_data = cv.extra_data(); if (extra_data != nullptr) { if (extra_data->side != ZERO) { // In case of a hyperbolic arc, make sure the point is located on the // same branch as the arc. - if (arc.sign_of_extra_data(p.x(), p.y()) != extra_data->side) + if (cv.sign_of_extra_data(p.x(), p.y()) != extra_data->side) return false; } else { // In case we have a segment of a line pair, make sure that p really // satisfies the equation of the line. - if (arc.sign_of_extra_data(p.x(), p.y()) != ZERO) return false; + if (cv.sign_of_extra_data(p.x(), p.y()) != ZERO) return false; } } // Act according to the conic degree. auto orient_f = m_alg_kernel->orientation_2_object(); - const auto& source = arc.source(); - const auto& target = arc.target(); - if (arc.orientation() == COLLINEAR) { + const auto& source = cv.source(); + const auto& target = cv.target(); + if (cv.orientation() == COLLINEAR) { Comparison_result res1; Comparison_result res2; @@ -2569,20 +2976,19 @@ public: // source and the target. return (orient_f(source, p, target) == COLLINEAR); } - 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. - if (arc.orientation() == COUNTERCLOCKWISE) - return (orient_f(source, p, target) == LEFT_TURN); - else - return (orient_f(source, p, target) == RIGHT_TURN); - } + + // 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. + if (cv.orientation() == COUNTERCLOCKWISE) + return (orient_f(source, p, target) == LEFT_TURN); + else + return (orient_f(source, p, target) == RIGHT_TURN); } /*! Build the data for hyperbolic arc, contaning the characterization of the * hyperbolic branch the arc is placed on. */ - void build_hyperbolic_arc_data(Curve_2& arc) const { + void build_hyperbolic_arc_data(Curve_2& cv) const { // Let phi be the rotation angle of the conic from its canonic form. // We can write: // @@ -2594,10 +3000,10 @@ public: // cos(2*phi) = ----------------------- // sqrt((r - s)^2 + t^2) // - const int or_fact = (arc.orientation() == CLOCKWISE) ? -1 : 1; - const Algebraic r = m_nt_traits->convert(or_fact * arc.r()); - const Algebraic s = m_nt_traits->convert(or_fact * arc.s()); - const Algebraic t = m_nt_traits->convert(or_fact * arc.t()); + const int or_fact = (cv.orientation() == CLOCKWISE) ? -1 : 1; + const Algebraic r = m_nt_traits->convert(or_fact * cv.r()); + const Algebraic s = m_nt_traits->convert(or_fact * cv.s()); + const Algebraic t = m_nt_traits->convert(or_fact * cv.t()); const Algebraic cos_2phi = (r - s) / m_nt_traits->sqrt((r-s)*(r-s) + t*t); const Algebraic zero = 0; const Algebraic one = 1; @@ -2642,8 +3048,8 @@ public: // 4*r*s - t^2 4*r*s - t^2 // // The denominator (4*r*s - t^2) must be negative for hyperbolas. - const Algebraic u = m_nt_traits->convert(or_fact * arc.u()); - const Algebraic v = m_nt_traits->convert(or_fact * arc.v()); + const Algebraic u = m_nt_traits->convert(or_fact * cv.u()); + const Algebraic v = m_nt_traits->convert(or_fact * cv.v()); const Algebraic det = 4*r*s - t*t; Algebraic x0, y0; @@ -2666,14 +3072,14 @@ public: // Make sure that the two endpoints are located on the same branch // of the hyperbola. - const auto& source = arc.source(); - const auto& target = arc.target(); - arc.set_extra_data(extra_data); - extra_data->side = arc.sign_of_extra_data(source.x(), source.y()); + const auto& source = cv.source(); + const auto& target = cv.target(); + cv.set_extra_data(extra_data); + extra_data->side = cv.sign_of_extra_data(source.x(), source.y()); CGAL_assertion(extra_data->side != ZERO); CGAL_assertion(extra_data->side == - arc.sign_of_extra_data(target.x(), target.y())); + cv.sign_of_extra_data(target.x(), target.y())); } /*! Find the x coordinates of the underlying conic at a given y coordinate. @@ -2682,16 +3088,16 @@ public: * \pre The vector xs must be allocated at the size of 2. * \return The number of x coordinates computed (either 0, 1 or 2). */ - int conic_get_x_coordinates(const Curve_2& arc, + int conic_get_x_coordinates(const Curve_2& cv, const Algebraic& y, Algebraic* xs) const { // 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 - Algebraic A = m_nt_traits->convert(arc.r()); + Algebraic A = m_nt_traits->convert(cv.r()); Algebraic B = - m_nt_traits->convert(arc.t())*y + m_nt_traits->convert(arc.u()); + m_nt_traits->convert(cv.t())*y + m_nt_traits->convert(cv.u()); Algebraic C = - (m_nt_traits->convert(arc.s())*y + m_nt_traits->convert(arc.v()))*y + - m_nt_traits->convert(arc.w()); + (m_nt_traits->convert(cv.s())*y + m_nt_traits->convert(cv.v()))*y + + m_nt_traits->convert(cv.w()); return solve_quadratic_equation(A, B, C, xs[0], xs[1]); } @@ -2702,16 +3108,16 @@ public: * \pre The vector ys must be allocated at the size of 2. * \return The number of y coordinates computed (either 0, 1 or 2). */ - int conic_get_y_coordinates(const Curve_2& arc, + int conic_get_y_coordinates(const Curve_2& cv, const Algebraic& x, Algebraic* ys) const { // 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 - Algebraic A = m_nt_traits->convert(arc.s()); + Algebraic A = m_nt_traits->convert(cv.s()); Algebraic B = - m_nt_traits->convert(arc.t())*x + m_nt_traits->convert(arc.v()); + m_nt_traits->convert(cv.t())*x + m_nt_traits->convert(cv.v()); Algebraic C = - (m_nt_traits->convert(arc.r())*x + m_nt_traits->convert(arc.u()))*x + - m_nt_traits->convert(arc.w()); + (m_nt_traits->convert(cv.r())*x + m_nt_traits->convert(cv.u()))*x + + m_nt_traits->convert(cv.w()); return solve_quadratic_equation(A, B, C, ys[0], ys[1]); } @@ -2760,17 +3166,17 @@ public: * \pre The vector ps should be allocated at the size of 2. * \return The number of points found. */ - int points_at_x(const Curve_2& arc, const Point_2& p, Point_2* ps) const { + int points_at_x(const Curve_2& cv, const Point_2& p, Point_2* ps) const { // Get the y coordinates of the points on the conic. Algebraic ys[2]; - int n = conic_get_y_coordinates(arc, p.x(), ys); + int n = conic_get_y_coordinates(cv, 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]); - if (arc.is_full_conic() || is_between_endpoints(arc, ps[m])) ++m; + if (cv.is_full_conic() || is_between_endpoints(cv, ps[m])) ++m; } // Return the number of points on the arc. @@ -2784,17 +3190,17 @@ public: * \pre The vector ps should be allocated at the size of 2. * \return The number of points found. */ - int points_at_y(const Curve_2& arc, const Point_2& p, Point_2* ps) const { + int points_at_y(const Curve_2& cv, const Point_2& p, Point_2* ps) const { // Get the y coordinates of the points on the conic. Algebraic xs[2]; - int n = conic_get_x_coordinates(arc, p.y(), xs); + int n = conic_get_x_coordinates(cv, 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()); - if (arc.is_full_conic() || is_between_endpoints(arc, ps[m])) ++m; + if (cv.is_full_conic() || is_between_endpoints(cv, ps[m])) ++m; } // Return the number of points on the arc. @@ -2917,24 +3323,19 @@ public: if ((xcv1.orientation() == COLLINEAR) && (xcv2.orientation() == COLLINEAR)) { // Construct the two supporting lines and compare them. - auto construct_line = m_alg_kernel->construct_line_2_object(); - typename Alg_kernel::Line_2 l1 = - construct_line(xcv1.source(), xcv1.target()); - typename Alg_kernel::Line_2 l2 = - construct_line(xcv2.source(), xcv2.target()); + auto ctr_line = m_alg_kernel->construct_line_2_object(); + typename Alg_kernel::Line_2 l1 = ctr_line(xcv1.source(), xcv1.target()); + typename Alg_kernel::Line_2 l2 = ctr_line(xcv2.source(), xcv2.target()); auto equal = m_alg_kernel->equal_2_object(); - if (equal(l1, l2)) return true; - // Try to compare l1 with the opposite of l2. - l2 = construct_line(xcv2.target(), xcv2.source()); - + l2 = ctr_line(xcv2.target(), xcv2.source()); return equal(l1, l2); } else if ((xcv1.orientation() == COLLINEAR) || (xcv2.orientation() == COLLINEAR)) { - // Only one arc is collinear, so the supporting curves cannot be the - // same: + // Only one arc is collinear; therefore so the supporting curves cannot + // be the same: return false; } @@ -2964,6 +3365,83 @@ public: CGAL::compare(xcv1.v() * factor2, xcv2.v() * factor1) == EQUAL && CGAL::compare(xcv1.w() * factor2, xcv2.w() * factor1) == EQUAL); } + + /*! Check whether the given point lies on the arc. + * \param p The qury point. + * \param (true) if p lies on the arc; (false) otherwise. + */ + bool contains_point(const X_monotone_curve_2& xcv, const Point_2& p) const { + // First check if p lies on the supporting conic. We first check whether + // it is one of p's generating conic curves. + bool p_on_conic(false); + if (p.is_generating_conic(xcv.id())) p_on_conic = true; + else { + // Check whether p satisfies the supporting conic equation. + p_on_conic = xcv.is_on_supporting_conic(p.x(), p.y()); + if (p_on_conic) { + // As p lies on the supporting conic of our arc, add its ID to + // the list of generating conics for p. + Point_2& p_non_const = const_cast(p); + p_non_const.set_generating_conic(xcv.id()); + } + } + + if (! p_on_conic) return false; + + // Check if p is between the endpoints of the arc. + return is_between_endpoints(xcv, p); + } + + /*! Calculate the vertical tangency points of the arc. + * \param vpts The vertical tangency points. + * \pre The vpts vector should be allocated at the size of 2. + * \return The number of vertical tangency points. + */ + int vertical_tangency_points(const Curve_2& cv, Alg_point_2* vpts) const { + // No vertical tangency points for line segments: + if (cv.orientation() == COLLINEAR) return 0; + + // Calculate the vertical tangency points of the supporting conic. + Alg_point_2 ps[2]; + auto n = cv.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) { + std::cout << ps[i] << std::endl; + if (cv.is_full_conic() || is_strictly_between_endpoints(cv, ps[i])) + vpts[m++] = ps[i]; + } + + // Return the number of vertical tangency points found. + CGAL_assertion(m <= 2); + return m; + } + + /*! Calculate the horizontal tangency points of the arc. + * \param hpts The horizontal tangency points. + * \pre The hpts vector should be allocated at the size of 2. + * \return The number of horizontal tangency points. + */ + int horizontal_tangency_points(const Curve_2& cv, Point_2* hpts) const { + // No horizontal tangency points for line segments: + if (cv.orientation() == COLLINEAR) return 0; + + // Calculate the horizontal tangency points of the conic. + Alg_point_2 ps[2]; + int n = cv.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 (cv.is_full_conic() || is_strictly_between_endpoints(cv, ps[i])) + hpts[m++] = ps[i]; + } + + // Return the number of horizontal tangency points found. + CGAL_assertion(m <= 2); + return m; + } }; #include diff --git a/Arrangement_on_surface_2/include/CGAL/Arr_geometry_traits/Conic_arc_2.h b/Arrangement_on_surface_2/include/CGAL/Arr_geometry_traits/Conic_arc_2.h index ac9a63b356e..5f17cad0d82 100644 --- a/Arrangement_on_surface_2/include/CGAL/Arr_geometry_traits/Conic_arc_2.h +++ b/Arrangement_on_surface_2/include/CGAL/Arr_geometry_traits/Conic_arc_2.h @@ -163,15 +163,15 @@ public: /// \name Get the arc properties. //@{ - /*! Check wheather the arc is valid. + /*! Determine wheather the arc is valid. */ bool is_valid() const { return test_flag(IS_VALID); } - /*! Check whether the arc represents a full conic curve. + /*! Determine whether the arc represents a full conic curve. */ bool is_full_conic() const { return test_flag(IS_FULL_CONIC); } - /*! Get the coefficients of the underlying conic. + /*! Obtain the coefficients of the underlying conic. */ const Integer& r() const { return (m_r); } const Integer& s() const { return (m_s); } @@ -180,22 +180,6 @@ public: const Integer& v() const { return (m_v); } const Integer& w() const { return (m_w); } - /*! Check whether the arc is x-monotone. - */ - bool is_x_monotone() const { - // Check if the arc contains no vertical tangency points. - Point_2 vtan_ps[2]; - return (vertical_tangency_points(vtan_ps) == 0); - } - - /*! Check whether the arc is y-monotone. - */ - bool is_y_monotone() const { - // Check if the arc contains no horizontal tangency points. - Point_2 htan_ps[2]; - return (horizontal_tangency_points(htan_ps) == 0); - } - /*! Obtain the arc's source. * \return The source point. * \pre The arc does not represent a full conic curve. @@ -305,11 +289,13 @@ public: //@} -public: +protected: + template friend class Arr_conic_traits_2; + /// \name Flag manipulation functions. //@{ template - constexpr size_t flag_mask(const T flag) const { return 0x1 << flag; } + static constexpr size_t flag_mask(const T flag) { return 0x1 << flag; } void reset_flags() { m_info = 0; } @@ -369,390 +355,6 @@ public: /*! Set the extra data field. */ void set_extra_data(Extra_data* extra_data) { m_extra_data = extra_data; } - - //@} - - /// \name Compute points on the arc. - //@{ - - /*! Calculate the vertical tangency points of the arc. - * \param vpts The vertical tangency points. - * \pre The vpts vector should be allocated at the size of 2. - * \return The number of vertical tangency points. - */ - int vertical_tangency_points(Point_2* vpts) const { - // No vertical tangency points for line segments: - if (m_orient == COLLINEAR) return 0; - - // Calculate the vertical tangency points of the supporting conic. - Point_2 ps[2]; - auto 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. - CGAL_assertion(m <= 2); - return m; - } - - /*! Calculate the horizontal tangency points of the arc. - * \param hpts The horizontal tangency points. - * \pre The hpts vector should be allocated at the size of 2. - * \return The number of horizontal tangency points. - */ - int horizontal_tangency_points(Point_2* hpts) const { - // No horizontal tangency points for line segments: - if (m_orient == COLLINEAR) return 0; - - // Calculate the horizontal tangency points of the conic. - Point_2 ps[2]; - int 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. - CGAL_assertion(m <= 2); - return m; - } - - /*! Find all points on the arc with a given x-coordinate. - * \param p A placeholder for the x-coordinate. - * \param ps The point on the arc at x(p). - * \pre The vector ps should be allocated at the size of 2. - * \return The number of points found. - */ - int points_at_x(const Point_2& p, Point_2* ps) const { - // Get the y coordinates of the points on the conic. - Algebraic ys[2]; - int 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]); - - if (is_full_conic() || is_between_endpoints(ps[m])) ++m; - } - - // Return the number of points on the arc. - CGAL_assertion(m <= 2); - return m; - } - - /*! Find all points on the arc with a given y-coordinate. - * \param p A placeholder for the y-coordinate. - * \param ps The point on the arc at x(p). - * \pre The vector ps should be allocated at the size of 2. - * \return The number of points found. - */ - int points_at_y(const Point_2& p, Point_2* ps) const { - // Get the y coordinates of the points on the conic. - Algebraic xs[2]; - int 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()); - if (is_full_conic() || is_between_endpoints(ps[m])) ++m; - } - - // Return the number of points on the arc. - CGAL_assertion(m <= 2); - return m; - } - //@} - -private: - /// \name Auxiliary construction functions. - //@{ - - /*! Set the properties of a conic arc (for the usage of the constructors). - * \param rat_coeffs A vector of size 6, storing the rational coefficients - * of x^2, y^2, xy, x, y and the free coefficient resp. - */ - void set(const Rational* rat_coeffs) { - set_flag(IS_VALID); - - // Convert the coefficients vector to an equivalent vector of integer - // coefficients. - Nt_traits nt_traits; - Integer int_coeffs[6]; - - nt_traits.convert_coefficients(rat_coeffs, rat_coeffs + 6, int_coeffs); - - // Check the orientation of conic curve, and negate the conic coefficients - // if its given orientation. - typename Rat_kernel::Conic_2 temp_conic(rat_coeffs[0], rat_coeffs[1], - rat_coeffs[2], rat_coeffs[3], - rat_coeffs[4], rat_coeffs[5]); - - if (m_orient == temp_conic.orientation()) { - m_r = int_coeffs[0]; - m_s = int_coeffs[1]; - m_t = int_coeffs[2]; - m_u = int_coeffs[3]; - m_v = int_coeffs[4]; - m_w = int_coeffs[5]; - } - else { - m_r = -int_coeffs[0]; - m_s = -int_coeffs[1]; - m_t = -int_coeffs[2]; - m_u = -int_coeffs[3]; - m_v = -int_coeffs[4]; - m_w = -int_coeffs[5]; - } - - // Make sure both endpoint lie on the supporting conic. - if (! is_on_supporting_conic(m_source) || - ! is_on_supporting_conic(m_target)) - { - reset_flags(); // inavlid arc - return; - } - - m_extra_data = nullptr; - - // Check whether we have a degree 2 curve. - if ((CGAL::sign(m_r) != ZERO) || (CGAL::sign(m_s) != ZERO) || - (CGAL::sign(m_t) != ZERO)) - { - if (m_orient == COLLINEAR) { - // We have a segment of a line pair with rational coefficients. - // Compose the equation of the underlying line - // (with algebraic coefficients). - const Algebraic x1 = m_source.x(), y1 = m_source.y(); - const Algebraic x2 = m_target.x(), y2 = m_target.y(); - - // The supporting line is A*x + B*y + C = 0, where: - // - // A = y2 - y1, B = x1 - x2, C = x2*y1 - x1*y2 - // - // We use the extra dat field to store the equation of this line. - m_extra_data = new Extra_data; - m_extra_data->a = y2 - y1; - m_extra_data->b = x1 - x2; - m_extra_data->c = x2*y1 - x1*y2; - m_extra_data->side = ZERO; - - // Make sure the midpoint is on the line pair (thus making sure that - // the two points are not taken from different lines). - Alg_kernel ker; - Point_2 p_mid = ker.construct_midpoint_2_object()(m_source, m_target); - - if (CGAL::sign((nt_traits.convert(m_r)*p_mid.x() + - nt_traits.convert(m_t)*p_mid.y() + - nt_traits.convert(m_u)) * p_mid.x() + - (nt_traits.convert(m_s)*p_mid.y() + - nt_traits.convert(m_v)) * p_mid.y() + - nt_traits.convert(m_w)) != ZERO) - { - reset_flags(); // inavlid arc - return; - } - } - else { - // The sign of (4rs - t^2) detetmines the conic type: - // - if it is possitive, the conic is an ellipse, - // - if it is negative, the conic is a hyperbola, - // - if it is zero, the conic is a parabola. - CGAL::Sign sign_conic = CGAL::sign(4*m_r*m_s - m_t*m_t); - - // Build the extra hyperbolic data if necessary - if (sign_conic == NEGATIVE) build_hyperbolic_arc_data(); - - if (sign_conic != POSITIVE) { - // In case of a non-degenerate parabola or a hyperbola, make sure - // the arc is not infinite. - Alg_kernel ker; - Point_2 p_mid = ker.construct_midpoint_2_object()(m_source, m_target); - Point_2 ps[2]; - - bool finite_at_x = (points_at_x(p_mid, ps) > 0); - bool finite_at_y = (points_at_y(p_mid, ps) > 0); - - if (! finite_at_x && ! finite_at_y) { - reset_flags(); // inavlid arc - return; - } - } - } - } - - - set_flag(IS_VALID); // mark that this arc valid - reset_flag(IS_FULL_CONIC); // mark that this arc is not a full conic - } - - /*! Set the properties of a conic arc that is really a full curve - * (that is, an ellipse). - * \param rat_coeffs A vector of size 6, storing the rational coefficients - * of x^2, y^2, xy, x, y and the free coefficient resp. - * \param comp_orient Should we compute the orientation of the given curve. - */ - void set_full(const Rational* rat_coeffs, const bool& comp_orient) { - // Convert the coefficients vector to an equivalent vector of integer - // coefficients. - Nt_traits nt_traits; - Integer int_coeffs[6]; - - nt_traits.convert_coefficients(rat_coeffs, rat_coeffs + 6, int_coeffs); - - // Check the orientation of conic curve, and negate the conic coefficients - // if its given orientation. - typename Rat_kernel::Conic_2 temp_conic(rat_coeffs[0], rat_coeffs[1], - rat_coeffs[2], rat_coeffs[3], - rat_coeffs[4], rat_coeffs[5]); - const Orientation temp_orient = temp_conic.orientation(); - - if (comp_orient) m_orient = temp_orient; - - if (m_orient == temp_orient) { - m_r = int_coeffs[0]; - m_s = int_coeffs[1]; - m_t = int_coeffs[2]; - m_u = int_coeffs[3]; - m_v = int_coeffs[4]; - m_w = int_coeffs[5]; - } - else { - m_r = -int_coeffs[0]; - m_s = -int_coeffs[1]; - m_t = -int_coeffs[2]; - m_u = -int_coeffs[3]; - m_v = -int_coeffs[4]; - m_w = -int_coeffs[5]; - } - - // Make sure the conic is a non-degenerate ellipse: - // The coefficients should satisfy (4rs - t^2) > 0. - const bool is_ellipse = (CGAL::sign(4*m_r*m_s - m_t*m_t) == POSITIVE); - CGAL_assertion(is_ellipse); - - // We do not have to store any extra data with the arc. - m_extra_data = nullptr; - - // Mark that this arc is a full conic curve. - if (is_ellipse) { - set_flag(IS_VALID); - set_flag(IS_FULL_CONIC); - } - else reset_flags(); // inavlid arc - } - - /*! 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) - // - Nt_traits nt_traits; - const int or_fact = (m_orient == CLOCKWISE) ? -1 : 1; - const Algebraic r = nt_traits.convert(or_fact * m_r); - const Algebraic s = nt_traits.convert(or_fact * m_s); - const Algebraic t = nt_traits.convert(or_fact * m_t); - const Algebraic cos_2phi = (r - s) / nt_traits.sqrt((r-s)*(r-s) + t*t); - const Algebraic zero = 0; - const Algebraic one = 1; - const Algebraic two = 2; - Algebraic sin_phi; - Algebraic cos_phi; - - // Calculate sin(phi) and cos(phi) according to the half-angle formulae: - // - // sin(phi)^2 = 0.5 * (1 - cos(2*phi)) - // cos(phi)^2 = 0.5 * (1 + cos(2*phi)) - Sign sign_t = CGAL::sign(t); - - if (sign_t == ZERO) { - // sin(2*phi) == 0, so phi = 0 or phi = PI/2 - if (CGAL::sign(cos_2phi) == POSITIVE) { - // phi = 0. - sin_phi = zero; - cos_phi = one; - } - else { - // phi = PI/2. - sin_phi = one; - cos_phi = zero; - } - } - else if (sign_t == POSITIVE) { - // sin(2*phi) > 0 so 0 < phi < PI/2. - sin_phi = nt_traits.sqrt((one + cos_2phi) / two); - cos_phi = nt_traits.sqrt((one - cos_2phi) / two); - } - else { - // sin(2*phi) < 0 so PI/2 < phi < PI. - sin_phi = nt_traits.sqrt((one + cos_2phi) / two); - cos_phi = -nt_traits.sqrt((one - cos_2phi) / two); - } - - // Calculate the center (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 Algebraic u = nt_traits.convert(or_fact * m_u); - const Algebraic v = nt_traits.convert(or_fact * m_v); - const Algebraic det = 4*r*s - t*t; - Algebraic x0, y0; - - CGAL_assertion(CGAL::sign(det) == NEGATIVE); - - 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 - // - // We store the equation of this line in the extra data structure and also - // the sign (side of half-plane) our arc occupies with respect to the line. - m_extra_data = new Extra_data; - - m_extra_data->a = cos_phi; - m_extra_data->b = sin_phi; - m_extra_data->c = - (cos_phi*x0 + sin_phi*y0); - - // Make sure that the two endpoints are located on the same branch - // of the hyperbola. - m_extra_data->side = sign_of_extra_data(m_source.x(), m_source.y()); - - CGAL_assertion(m_extra_data->side != ZERO); - CGAL_assertion(m_extra_data->side == - sign_of_extra_data(m_target.x(), m_target.y())); - } //@} public: @@ -790,90 +392,6 @@ protected: return (CGAL::sign(val) == ZERO); } - /*! Check whether the given point is between the source and the target. - * The point is assumed to be on the conic's boundary. - * \param p The query point. - * \return true if the point is between the two endpoints, - * (false) if it is not. - */ - bool is_between_endpoints(const Point_2& p) const { - CGAL_precondition(! is_full_conic()); - - // Check if p is one of the endpoints. - Alg_kernel ker; - - if (ker.equal_2_object()(p, m_source) || ker.equal_2_object()(p, m_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. - * \param p The query point. - * \return true if the point is strictly between the two endpoints, - * (false) if it is not. - */ - 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; - - // Check if we have extra data available. - if (m_extra_data != nullptr) { - if (m_extra_data->side != ZERO) { - // In case of a hyperbolic arc, make sure the point is located on the - // same branch as the arc. - if (sign_of_extra_data(p.x(), p.y()) != m_extra_data->side) - return false; - } - else { - // In case we have a segment of a line pair, make sure that p really - // satisfies the equation of the line. - if (sign_of_extra_data(p.x(), p.y()) != ZERO) return false; - } - } - - // Act according to the conic degree. - Alg_kernel ker; - - if (m_orient == COLLINEAR) { - Comparison_result res1; - Comparison_result res2; - - if (ker.compare_x_2_object()(m_source, m_target) == EQUAL) { - // In case of a vertical segment - just check whether the y coordinate - // of p is between those of the source's and of the target's. - res1 = ker.compare_y_2_object()(p, m_source); - res2 = ker.compare_y_2_object()(p, m_target); - } - else { - // Otherwise, since the segment is x-monotone, just check whether the - // x coordinate of p is between those of the source's and of the - // target's. - res1 = ker.compare_x_2_object()(p, m_source); - res2 = ker.compare_x_2_object()(p, m_target); - } - - // If p is not in the (open) x-range (or y-range) of the segment, it - // cannot be contained in the segment. - if ((res1 == EQUAL) || (res2 == EQUAL) || (res1 == res2)) return false; - - // Perform an orientation test: This is crucial for segment of line - // pairs, as we want to make sure that p lies on the same line as the - // source and the target. - return (ker.orientation_2_object()(m_source, p, m_target) == COLLINEAR); - } - 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. - if (m_orient == COUNTERCLOCKWISE) - return (ker.orientation_2_object()(m_source, p, m_target) == LEFT_TURN); - else - return (ker.orientation_2_object()(m_source, p, m_target) == RIGHT_TURN); - } - } - /*! Find the vertical tangency points of the undelying conic. * \param ps The output points of vertical tangency. * This area must be allocated at the size of 2. @@ -930,7 +448,7 @@ protected: ++n; } else { - for (int j = 0; j < n_ys; j++) { + for (int j = 0; j < n_ys; ++j) { if (CGAL::compare(nt_traits.convert(two*m_s) * ys[j], -(nt_traits.convert(m_t) * xs[i] + nt_traits.convert(m_v))) == EQUAL) @@ -953,7 +471,7 @@ protected: * \return The number of horizontal tangency points. */ int conic_horizontal_tangency_points(Point_2* ps) const { - const Integer zero = 0; + const Integer zero(0); // In case the base conic is of degree 1 (and not 2), the arc has no // vertical tangency points. @@ -968,7 +486,7 @@ protected: const Integer four = 4; int n; Algebraic ys[2]; - Algebraic *ys_end; + Algebraic* ys_end; Nt_traits nt_traits; ys_end = nt_traits.solve_quadratic_equation(m_t*m_t - four*m_r*m_s, @@ -990,42 +508,6 @@ protected: return n; } - /*! Find the y coordinates of the underlying conic at a given x coordinate. - * \param x The x coordinate. - * \param ys The output y coordinates. - * \pre The vector ys must be allocated at the size of 2. - * \return The number of y coordinates computed (either 0, 1 or 2). - */ - int conic_get_y_coordinates(const Algebraic& x, Algebraic* ys) const { - // 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 - Nt_traits nt_traits; - Algebraic A = nt_traits.convert(m_s); - Algebraic B = nt_traits.convert(m_t)*x + nt_traits.convert(m_v); - Algebraic C = (nt_traits.convert(m_r)*x + nt_traits.convert(m_u))*x + - nt_traits.convert(m_w); - - return (solve_quadratic_equation(A, B, C, ys[0], ys[1])); - } - - /*! Find the x coordinates of the underlying conic at a given y coordinate. - * \param y The y coordinate. - * \param xs The output x coordinates. - * \pre The vector xs must be allocated at the size of 2. - * \return The number of x coordinates computed (either 0, 1 or 2). - */ - int conic_get_x_coordinates(const Algebraic& y, Algebraic* xs) const { - // 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 - Nt_traits nt_traits; - Algebraic A = nt_traits.convert(m_r); - Algebraic B = nt_traits.convert(m_t)*y + nt_traits.convert(m_u); - Algebraic C = (nt_traits.convert(m_s)*y + nt_traits.convert(m_v))*y + - nt_traits.convert(m_w); - - return (solve_quadratic_equation(A, B, C, xs[0], xs[1])); - } - /*! Solve the given quadratic equation: Ax^2 + B*x + C = 0. * \param x_minus The root obtained from taking -sqrt(discriminant). * \param x_plus The root obtained from taking -sqrt(discriminant). @@ -1064,7 +546,6 @@ protected: return 2; } //@} - }; /*! Exporter for conic arcs. diff --git a/Arrangement_on_surface_2/include/CGAL/Arr_geometry_traits/Conic_x_monotone_arc_2.h b/Arrangement_on_surface_2/include/CGAL/Arr_geometry_traits/Conic_x_monotone_arc_2.h index 16eef6de002..cd993d7e7d0 100644 --- a/Arrangement_on_surface_2/include/CGAL/Arr_geometry_traits/Conic_x_monotone_arc_2.h +++ b/Arrangement_on_surface_2/include/CGAL/Arr_geometry_traits/Conic_x_monotone_arc_2.h @@ -45,14 +45,8 @@ public: // Type definition for the intersection points mapping. typedef typename Conic_point_2::Conic_id Conic_id; - typedef std::pair Conic_pair; - typedef std::pair Intersection_point; - typedef std::list Intersection_list; using Conic_arc_2::sign_of_extra_data; - using Conic_arc_2::is_between_endpoints; - using Conic_arc_2::is_strictly_between_endpoints; - using Conic_arc_2::conic_get_y_coordinates; using Conic_arc_2::IS_VALID; using Conic_arc_2::IS_FULL_CONIC; @@ -64,21 +58,6 @@ public: using Conic_arc_2::flip_flag; using Conic_arc_2::test_flag; - /*! \struct Less functor for Conic_pair. - */ - struct Less_conic_pair { - bool operator()(const Conic_pair& cp1, const Conic_pair& cp2) const { - // Compare the pairs of IDs lexicographically. - return ((cp1.first < cp2.first) || - ((cp1.first == cp2.first) && (cp1.second < cp2.second))); - } - }; - - typedef std::map - Intersection_map; - typedef typename Intersection_map::value_type Intersection_map_entry; - typedef typename Intersection_map::iterator Intersection_map_iterator; - private: typedef Conic_arc_2 Base; @@ -112,7 +91,7 @@ public: FACING_MASK = (0x1 << FACING_UP) | (0x1 << FACING_DOWN) }; - /// \name Constrcution methods. + /// \name Public constrcutors, assignment operators, and destructors. //@{ /*! Default constructor. @@ -136,14 +115,42 @@ public: m_id(arc.m_id) {} + /*! Assignment operator. + * \param arc The copied arc. + */ + const Self& operator=(const Self& arc) { + CGAL_precondition (arc.is_valid()); + + if (this == &arc) return (*this); + + // Copy the base arc. + Base::operator= (arc); + + // Set the rest of the properties. + m_alg_r = arc.m_alg_r; + m_alg_s = arc.m_alg_s; + m_alg_t = arc.m_alg_t; + m_alg_u = arc.m_alg_u; + m_alg_v = arc.m_alg_v; + m_alg_w = arc.m_alg_w; + + m_id = arc.m_id; + + return (*this); + } + //@} + +private: + template friend class Arr_conic_traits_2; + + /// \name private constrcutors to be used only by the traits class template. + //@{ + /*! Construct an x-monotone arc from a conic arc. * \param arc The given (base) arc. * \pre The given arc is x-monotone. */ - _Conic_x_monotone_arc_2(const Base& arc) : - Base(arc), - m_id() - {} + _Conic_x_monotone_arc_2(const Base& arc) : Base(arc), m_id() {} /*! Construct an x-monotone arc from a conic arc. * \param arc The given (base) arc. @@ -189,31 +196,16 @@ public: Base() {} - /*! Assignment operator. - * \param arc The copied arc. + /*! Obtain the non const reference to the curve source point. */ - const Self& operator=(const Self& arc) { - CGAL_precondition (arc.is_valid()); + Conic_point_2& source() { return (this->m_source); } - if (this == &arc) return (*this); - - // Copy the base arc. - Base::operator= (arc); - - // Set the rest of the properties. - m_alg_r = arc.m_alg_r; - m_alg_s = arc.m_alg_s; - m_alg_t = arc.m_alg_t; - m_alg_u = arc.m_alg_u; - m_alg_v = arc.m_alg_v; - m_alg_w = arc.m_alg_w; - - m_id = arc.m_id; - - return (*this); - } + /*! Obtain the non const reference to the curve source point. + */ + Conic_point_2& target() { return this->m_target; } //@} +public: /// \name Accessing the arc properties. //@{ @@ -238,13 +230,11 @@ public: * \return The source point. */ const Conic_point_2& source() const { return (this->m_source); } - Conic_point_2& source() { return (this->m_source); } /*! Obtain the arc's target. * \return The target point. */ const Conic_point_2& target() const { return this->m_target; } - Conic_point_2& target() { return this->m_target; } /*! Obtain the orientation of the arc. * \return The orientation. @@ -265,11 +255,37 @@ public: else return this->m_source; } - /*! Return true iff the conic arc is directed right iexicographically. + /*! Determine whether the conic arc is directed iexicographically right. */ - bool is_directed_right() const - { return test_flag(IS_DIRECTED_RIGHT); } + bool is_directed_right() const { return test_flag(IS_DIRECTED_RIGHT); } + /*! Determine whether the conic arc is a vertical segment. + */ + bool is_vertical() const { return test_flag(IS_VERTICAL_SEGMENT); } + + /*! Determine whether the conic arc is a facing up. + */ + bool is_upper() const { return test_flag(FACING_UP); } + + /*! Determine whether the conic arc is a facing down. + */ + bool is_lower() const { return test_flag(FACING_DOWN); } + + /*! Check whether the arc is a special segment connecting two algebraic + * endpoints (and has no undelying integer conic coefficients). + */ + bool is_special_segment() const { return test_flag(IS_SPECIAL_SEGMENT); } + + /*! Obtain the mask of the DEGREE_1 flag. + */ + static constexpr size_t degree_1_mask() { return flag_mask(DEGREE_1); } + + /*! Obtain the mask of the DEGREE_1 flag. + */ + static constexpr size_t degree_2_mask() { return flag_mask(DEGREE_2); } + + /*! Obtain the algebraic coefficients. + */ Algebraic alg_r() const { return m_alg_r; } Algebraic alg_s() const { return m_alg_s; } Algebraic alg_t() const { return m_alg_t; } @@ -289,6 +305,9 @@ public: // Setters //@{ + + /*! Set the algebraic coefficients. + */ void set_alg_coefficients(const Algebraic& alg_r, const Algebraic& alg_s, const Algebraic& alg_t, const Algebraic& alg_u, const Algebraic& alg_v, const Algebraic& alg_w) @@ -301,49 +320,12 @@ public: m_alg_w = alg_w; } - /*! Add a generating conic ID. */ + /*! Add a generating conic ID. + */ void set_generating_conic(const Conic_id& id) { this->m_source.set_generating_conic(id); this->m_target.set_generating_conic(id); } - - //@} - - /// \name Predicates. - //@{ - - /*! Check if the conic arc is a vertical segment. - */ - bool is_vertical() const - { return test_flag(IS_VERTICAL_SEGMENT); } - - /*! Check whether the given point lies on the arc. - * \param p The qury point. - * \param (true) if p lies on the arc; (false) otherwise. - */ - bool contains_point(const Conic_point_2& p) const { - // First check if p lies on the supporting conic. We first check whether - // it is one of p's generating conic curves. - bool p_on_conic(false); - - if (p.is_generating_conic(m_id)) p_on_conic = true; - else { - // Check whether p satisfies the supporting conic equation. - p_on_conic = is_on_supporting_conic(p.x(), p.y()); - - if (p_on_conic) { - // As p lies on the supporting conic of our arc, add its ID to - // the list of generating conics for p. - Conic_point_2& p_non_const = const_cast(p); - p_non_const.set_generating_conic(m_id); - } - } - - if (! p_on_conic) return false; - - // Check if p is between the endpoints of the arc. - return is_between_endpoints(p); - } //@} /// \name Constructing points on the arc. @@ -413,92 +395,6 @@ public: return oi; } - /*! Compute the intersections with the given arc. - * \param arc The given intersecting arc. - * \param inter_map Maps conic pairs to lists of their intersection points. - * \param oi The output iterator. - * \return The past-the-end iterator. - */ - template - OutputIterator intersect(const Self& arc, Intersection_map& inter_map, - OutputIterator oi) const - { - typedef boost::variant Intersection_result; - - if (has_same_supporting_conic(arc)) { - // Check for overlaps between the two arcs. - Self overlap; - - if (compute_overlap(arc, overlap)) { - // There can be just a single overlap between two x-monotone arcs: - *oi++ = Intersection_result(overlap); - return oi; - } - - // In case there is not overlap and the supporting conics are the same, - // there cannot be any intersection points, unless the two arcs share - // an end point. - // Note that in this case we do not define the multiplicity of the - // intersection points we report. - Alg_kernel ker; - - if (ker.equal_2_object()(left(), arc.left())) { - Intersection_point ip(left(), 0); - *oi++ = Intersection_result(ip); - } - - if (ker.equal_2_object()(right(), arc.right())) { - Intersection_point ip(right(), 0); - *oi++ = Intersection_result(ip); - } - - return oi; - } - - // Search for the pair of supporting conics in the map (the first conic - // ID in the pair should be smaller than the second one, to guarantee - // uniqueness). - Conic_pair conic_pair; - Intersection_map_iterator map_iter; - Intersection_list inter_list; - bool invalid_ids = false; - - if (m_id.is_valid() && arc.m_id.is_valid()) { - if (m_id < arc.m_id) conic_pair = Conic_pair(m_id, arc.m_id); - else conic_pair = Conic_pair(arc.m_id, m_id); - map_iter = inter_map.find(conic_pair); - } - else { - // In case one of the IDs is invalid, we do not look in the map neither - // we cache the results. - map_iter = inter_map.end(); - invalid_ids = true; - } - - if (map_iter == inter_map.end()) { - // In case the intersection points between the supporting conics have - // not been computed before, compute them now and store them in the map. - intersect_supporting_conics(arc, inter_list); - - if (! invalid_ids) inter_map[conic_pair] = inter_list; - } - else { - // Obtain the precomputed intersection points from the map. - inter_list = (*map_iter).second; - } - - // Go over the list of intersection points and report those that lie on - // both x-monotone arcs. - for (auto iter = inter_list.begin(); iter != inter_list.end(); ++iter) { - if (is_between_endpoints((*iter).first) && - arc.is_between_endpoints((*iter).first)) - { - *oi++ = Intersection_result(*iter); - } - } - - return oi; - } //@} /// \name Constructing x-monotone arcs. @@ -524,161 +420,12 @@ public: return arc; } - - /*! Check whether the two arcs are equal (have the same graph). - * \param arc The compared arc. - * \return (true) if the two arcs have the same graph; (false) otherwise. - */ - bool equals(const Self& arc) const { - // The two arc must have the same supporting conic curves. - if (! has_same_supporting_conic (arc)) return false; - - // Check that the arc endpoints are the same. - Alg_kernel ker; - - if (this->m_orient == COLLINEAR) { - CGAL_assertion(arc.m_orient == COLLINEAR); - return((ker.equal_2_object()(this->m_source, arc.m_source) && - ker.equal_2_object()(this->m_target, arc.m_target)) || - (ker.equal_2_object()(this->m_source, arc.m_target) && - ker.equal_2_object()(this->m_target, arc.m_source))); - } - - if (this->m_orient == arc.m_orient) { - // Same orientation - the source and target points must be the same. - return (ker.equal_2_object()(this->m_source, arc.m_source) && - ker.equal_2_object()(this->m_target, arc.m_target)); - } - else { - // Reverse orientation - the source and target points must be swapped. - return (ker.equal_2_object()(this->m_source, arc.m_target) && - ker.equal_2_object()(this->m_target, arc.m_source)); - } - } - - bool is_upper() const { return test_flag(FACING_UP); } - - bool is_lower() const { return test_flag(FACING_DOWN); } - - /*! Check whether the arc is a special segment connecting two algebraic - * endpoints (and has no undelying integer conic coefficients). - */ - bool is_special_segment() const { return test_flag(IS_SPECIAL_SEGMENT); } - //@} private: /// \name Auxiliary (private) functions. //@{ - /*! Set the properties of the x-monotone conic arc (for the usage of the - * constructors). - */ - void set_x_monotone() { - // Convert the coefficients of the supporting conic to algebraic numbers. - Nt_traits nt_traits; - - m_alg_r = nt_traits.convert(this->m_r); - m_alg_s = nt_traits.convert(this->m_s); - m_alg_t = nt_traits.convert(this->m_t); - m_alg_u = nt_traits.convert(this->m_u); - m_alg_v = nt_traits.convert(this->m_v); - m_alg_w = nt_traits.convert(this->m_w); - - // Set the generating conic ID for the source and target points. - this->m_source.set_generating_conic(m_id); - this->m_target.set_generating_conic(m_id); - - // Update the m_info bits. - set_flag(IS_VALID); - reset_flag(IS_FULL_CONIC); - - // Check if the arc is directed right (the target is lexicographically - // greater than the source point), or to the left. - Alg_kernel ker; - Comparison_result dir_res = - ker.compare_xy_2_object()(this->m_source, this->m_target); - - CGAL_assertion(dir_res != EQUAL); - - if (dir_res == SMALLER) set_flag(IS_DIRECTED_RIGHT); - - // Compute the degree of the underlying conic. - if ((CGAL::sign(this->m_r) != ZERO) || - (CGAL::sign(this->m_s) != ZERO) || - (CGAL::sign(this->m_t) != ZERO)) - { - set_flag(DEGREE_2); - - if (this->m_orient == COLLINEAR) { - set_flag(IS_SPECIAL_SEGMENT); - - // Check whether the arc is a vertical segment: - if (ker.compare_x_2_object()(this->m_source, this->m_target) == EQUAL) - set_flag(IS_VERTICAL_SEGMENT); - - return; - } - } - else { - CGAL_assertion(CGAL::sign(this->m_u) != ZERO || - CGAL::sign(this->m_v) != ZERO); - - if (CGAL::sign(this->m_v) == ZERO) { - - // The supporting curve is of the form: _u*x + _w = 0 - set_flag(IS_VERTICAL_SEGMENT); - } - - set_flag(DEGREE_1); - - return; - } - - if (this->m_orient == COLLINEAR) return; - - // Compute a midpoint between the source and the target and get the y-value - // of the arc at its x-coordiante. - Point_2 p_mid = - ker.construct_midpoint_2_object()(this->m_source, this->m_target); - Algebraic ys[2]; - CGAL_assertion_code(int n_ys = ) - conic_get_y_coordinates(p_mid.x(), ys); - - CGAL_assertion(n_ys != 0); - - // Check which solution lies on the x-monotone arc. - Point_2 p_arc_mid(p_mid.x(), ys[0]); - - if (is_strictly_between_endpoints(p_arc_mid)) { - // Mark that we should use the -sqrt(disc) root for points on this - // x-monotone arc. - reset_flag(PLUS_SQRT_DISC_ROOT); - } - else { - CGAL_assertion (n_ys == 2); - p_arc_mid = Point_2 (p_mid.x(), ys[1]); - - CGAL_assertion(is_strictly_between_endpoints(p_arc_mid)); - - // Mark that we should use the +sqrt(disc) root for points on this - // x-monotone arc. - set_flag(PLUS_SQRT_DISC_ROOT); - } - - // Check whether the conic is facing up or facing down: - // 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 coordinate of a point on the arc whose x - // coordinate is (x1+x2)/2 and compare it to (y1+y2)/2. - Comparison_result res = ker.compare_y_2_object() (p_arc_mid, p_mid); - - // If the arc is above the connecting segment, so it is facing upwards. - if (res == LARGER) set_flag(FACING_UP); - // If the arc is below the connecting segment, so it is facing downwards. - else if (res == SMALLER) set_flag(FACING_DOWN); - } - /*! Check whether the given point lies on the supporting conic of the arc. * \param px The x-coordinate of query point. * \param py The y-coordinate of query point. @@ -702,66 +449,6 @@ private: return (_sign == ZERO); } - /*! Check whether the two arcs have the same supporting conic. - * \param arc The compared arc. - * \return (true) if the two supporting conics are the same. - */ - bool has_same_supporting_conic(const Self& arc) const { - // Check if the two arcs originate from the same conic: - if (m_id == arc.m_id && m_id.is_valid() && arc.m_id.is_valid()) return true; - - // In case both arcs are collinear, check if they have the same - // supporting lines. - if ((this->m_orient == COLLINEAR) && (arc.m_orient == COLLINEAR)) { - // Construct the two supporting lines and compare them. - Alg_kernel ker; - auto construct_line = ker.construct_line_2_object(); - typename Alg_kernel::Line_2 l1 = - construct_line(this->m_source, this->m_target); - typename Alg_kernel::Line_2 l2 = - construct_line(arc.m_source, arc.m_target); - auto equal = ker.equal_2_object(); - - if (equal(l1, l2)) return true; - - // Try to compare l1 with the opposite of l2. - l2 = construct_line(arc.m_target, arc.m_source); - - return equal(l1, l2); - } - else if ((this->m_orient == COLLINEAR) || (arc.m_orient == COLLINEAR)) { - // Only one arc is collinear, so the supporting curves cannot be the - // same: - return false; - } - - // Check whether the coefficients of the two supporting conics are equal - // up to a constant factor. - Integer factor1 = 1; - Integer factor2 = 1; - - if (CGAL::sign(this->m_r) != ZERO) factor1 = this->m_r; - else if (CGAL::sign(this->m_s) != ZERO) factor1 = this->m_s; - else if (CGAL::sign(this->m_t) != ZERO) factor1 = this->m_t; - else if (CGAL::sign(this->m_u) != ZERO) factor1 = this->m_u; - else if (CGAL::sign(this->m_v) != ZERO) factor1 = this->m_v; - else if (CGAL::sign(this->m_w) != ZERO) factor1 = this->m_w; - - if (CGAL::sign(arc.m_r) != ZERO) factor2 = arc.m_r; - else if (CGAL::sign(arc.m_s) != ZERO) factor2 = arc.m_s; - else if (CGAL::sign(arc.m_t) != ZERO) factor2 = arc.m_t; - else if (CGAL::sign(arc.m_u) != ZERO) factor2 = arc.m_u; - else if (CGAL::sign(arc.m_v) != ZERO) factor2 = arc.m_v; - else if (CGAL::sign(arc.m_w) != ZERO) factor2 = arc.m_w; - - return (CGAL::compare(this->m_r * factor2, arc.m_r * factor1) == EQUAL && - CGAL::compare(this->m_s * factor2, arc.m_s * factor1) == EQUAL && - CGAL::compare(this->m_t * factor2, arc.m_t * factor1) == EQUAL && - CGAL::compare(this->m_u * factor2, arc.m_u * factor1) == EQUAL && - CGAL::compare(this->m_v * factor2, arc.m_v * factor1) == EQUAL && - CGAL::compare(this->m_w * factor2, arc.m_w * factor1) == EQUAL); - } - public: /*! Obtain the i'th order derivative by x of the conic at the point p=(x,y). * \param p The point where we derive. @@ -886,7 +573,7 @@ public: // so their first-order derivative by x is simply -b/a. The higher-order // derivatives are all 0. if (i == 1) { - if (CGAL::sign (this->m_extra_data->a) != NEGATIVE) { + if (CGAL::sign(this->m_extra_data->a) != NEGATIVE) { slope_numer = - this->m_extra_data->b; slope_denom = this->m_extra_data->a; } @@ -982,264 +669,14 @@ public: } private: - /*! Compute the overlap with a given arc, which is supposed to have the same - * supporting conic curve as this arc. - * \param arc The given arc. - * \param overlap Output: The overlapping arc (if any). - * \return Whether we found an overlap. - */ - bool compute_overlap(const Self& arc, Self& overlap) const { - // Check if the two arcs are identical. - if (equals(arc)) { - overlap = arc; - return true; - } - - if (is_strictly_between_endpoints(arc.left())) { - if (is_strictly_between_endpoints(arc.right())) { - // Case 1 - *this: +-----------> - // arc: +=====> - overlap = arc; - return true; - } - else { - // Case 2 - *this: +-----------> - // arc: +=====> - overlap = *this; - - if (overlap.test_flag(IS_DIRECTED_RIGHT)) - overlap.m_source = arc.left(); - else - overlap.m_target = arc.left(); - - return true; - } - } - else if (is_strictly_between_endpoints(arc.right())) { - // Case 3 - *this: +-----------> - // arc: +=====> - overlap = *this; - - if (overlap.test_flag(IS_DIRECTED_RIGHT)) - overlap.m_target = arc.right(); - else - overlap.m_source = arc.right(); - - return true; - } - else if (arc.is_between_endpoints(this->m_source) && - arc.is_between_endpoints(this->m_target) && - (arc.is_strictly_between_endpoints(this->m_source) || - arc.is_strictly_between_endpoints(this->m_target))) - { - // Case 4 - *this: +-----------> - // arc: +================> - overlap = *this; - return true; - } - - // If we reached here, there are no overlaps: - return false; - } - - /*! Intersect the supporing conic curves of this arc and the given arc. - * \param arc The arc to intersect with. - * \param inter_list The list of intersection points. - */ - void intersect_supporting_conics(const Self& arc, - Intersection_list& inter_list) const { - if (is_special_segment() && ! arc.is_special_segment()) { - // If one of the arcs is a special segment, make sure it is (arc). - arc.intersect_supporting_conics(*this, inter_list); - return; - } - - const int deg1 = ((this->m_info & DEGREE_MASK) == flag_mask(DEGREE_1)) ? 1 : 2; - const int deg2 = ((arc.m_info & DEGREE_MASK) == flag_mask(DEGREE_1)) ? 1 : 2; - Nt_traits nt_traits; - Algebraic xs[4]; - int n_xs = 0; - Algebraic ys[4]; - int n_ys = 0; - - if (arc.is_special_segment()) { - // The second arc is a special segment (a*x + b*y + c = 0). - if (is_special_segment()) { - // Both arc are sepcial segment, so they have at most one intersection - // point. - Algebraic denom = this->m_extra_data->a * arc.m_extra_data->b - - this->m_extra_data->b * arc.m_extra_data->a; - - if (CGAL::sign (denom) != CGAL::ZERO) { - xs[0] = (this->m_extra_data->b * arc.m_extra_data->c - - this->m_extra_data->c * arc.m_extra_data->b) / denom; - n_xs = 1; - - ys[0] = (this->m_extra_data->c * arc.m_extra_data->a - - this->m_extra_data->a * arc.m_extra_data->c) / denom; - n_ys = 1; - } - } - else { - // Compute the x-coordinates of the intersection points. - n_xs = compute_resultant_roots(nt_traits, - m_alg_r, m_alg_s, m_alg_t, - m_alg_u, m_alg_v, m_alg_w, - deg1, - arc.m_extra_data->a, - arc.m_extra_data->b, - arc.m_extra_data->c, - xs); - CGAL_assertion (n_xs <= 2); - - // Compute the y-coordinates of the intersection points. - n_ys = compute_resultant_roots(nt_traits, - m_alg_s, m_alg_r, m_alg_t, - m_alg_v, m_alg_u, m_alg_w, - deg1, - arc.m_extra_data->b, - arc.m_extra_data->a, - arc.m_extra_data->c, - ys); - CGAL_assertion(n_ys <= 2); - } - } - else { - // Compute the x-coordinates of the intersection points. - n_xs = compute_resultant_roots(nt_traits, - this->m_r, this->m_s, this->m_t, - this->m_u, this->m_v, this->m_w, - deg1, - arc.m_r, arc.m_s, arc.m_t, - arc.m_u, arc.m_v, arc.m_w, - deg2, - xs); - CGAL_assertion(n_xs <= 4); - - // Compute the y-coordinates of the intersection points. - n_ys = compute_resultant_roots(nt_traits, - this->m_s, this->m_r, this->m_t, - this->m_v, this->m_u, this->m_w, - deg1, - arc.m_s, arc.m_r, arc.m_t, - arc.m_v, arc.m_u, arc.m_w, - deg2, - ys); - CGAL_assertion(n_ys <= 4); - } - - // Pair the coordinates of the intersection points. As the vectors of - // x and y-coordinates are sorted in ascending order, we output the - // intersection points in lexicographically ascending order. - unsigned int mult; - int i, j; - - if (arc.is_special_segment()) { - if ((n_xs == 0) || (n_ys == 0)) return; - - if ((n_xs == 1) && (n_ys == 1)) { - // Single intersection. - Conic_point_2 ip (xs[0], ys[0]); - - ip.set_generating_conic(m_id); - ip.set_generating_conic (arc.m_id); - - // In case the other curve is of degree 2, this is a tangency point. - mult = ((deg1 == 1) || is_special_segment()) ? 1 : 2; - inter_list.push_back(Intersection_point(ip, mult)); - } - else if ((n_xs == 1) && (n_ys == 2)) { - Conic_point_2 ip1(xs[0], ys[0]); - - ip1.set_generating_conic(m_id); - ip1.set_generating_conic(arc.m_id); - - inter_list.push_back(Intersection_point(ip1, 1)); - - Conic_point_2 ip2 (xs[0], ys[1]); - - ip2.set_generating_conic(m_id); - ip2.set_generating_conic(arc.m_id); - - inter_list.push_back(Intersection_point(ip2, 1)); - } - else if ((n_xs == 2) && (n_ys == 1)) { - Conic_point_2 ip1 (xs[0], ys[0]); - - ip1.set_generating_conic(m_id); - ip1.set_generating_conic (arc.m_id); - - inter_list.push_back(Intersection_point(ip1, 1)); - - Conic_point_2 ip2 (xs[1], ys[0]); - - ip2.set_generating_conic(m_id); - ip2.set_generating_conic(arc.m_id); - - inter_list.push_back(Intersection_point(ip2, 1)); - - } - else { - CGAL_assertion((n_xs == 2) && (n_ys == 2)); - - // The x-coordinates and the y-coordinates are given in ascending - // order. If the slope of the segment is positive, we pair the - // coordinates as is - otherwise, we swap the pairs. - int ind_first_y = 0, ind_second_y = 1; - - if (CGAL::sign(arc.m_extra_data->b) == CGAL::sign(arc.m_extra_data->a)) { - ind_first_y = 1; - ind_second_y = 0; - } - - Conic_point_2 ip1(xs[0], ys[ind_first_y]); - - ip1.set_generating_conic(m_id); - ip1.set_generating_conic(arc.m_id); - - inter_list.push_back(Intersection_point(ip1, 1)); - - Conic_point_2 ip2(xs[1], ys[ind_second_y]); - - ip2.set_generating_conic(m_id); - ip2.set_generating_conic(arc.m_id); - - inter_list.push_back(Intersection_point(ip2, 1)); - } - - return; - } - - for (i = 0; i < n_xs; ++i) { - for (j = 0; j < n_ys; ++j) { - if (is_on_supporting_conic (xs[i], ys[j]) && - arc.is_on_supporting_conic(xs[i], ys[j])) - { - // Create the intersection point and set its generating conics. - Conic_point_2 ip(xs[i], ys[j]); - - ip.set_generating_conic(m_id); - ip.set_generating_conic(arc.m_id); - - // Compute the multiplicity of the intersection point. - if (deg1 == 1 && deg2 == 1) mult = 1; - else mult = multiplicity_of_intersection_point(arc, ip); - - // Insert the intersection point to the output list. - inter_list.push_back(Intersection_point(ip, mult)); - } - } - } - } - /*! Compute the multiplicity of an intersection point. * \param arc The arc to intersect with. * \param p The intersection point. * \return The multiplicity of the intersection point. */ - unsigned int multiplicity_of_intersection_point(const Self& arc, + unsigned int multiplicity_of_intersection_point(const Self& xcv, const Point_2& p) const { - CGAL_assertion(! is_special_segment() || ! arc.is_special_segment()); + CGAL_assertion(! is_special_segment() || ! xcv.is_special_segment()); // Compare the slopes of the two arcs at p, using their first-order // partial derivatives. @@ -1247,7 +684,7 @@ private: Algebraic slope2_numer, slope2_denom; derive_by_x_at(p, 1, slope1_numer, slope1_denom); - arc.derive_by_x_at(p, 1, slope2_numer, slope2_denom); + xcv.derive_by_x_at(p, 1, slope2_numer, slope2_denom); if (CGAL::compare(slope1_numer*slope2_denom, slope2_numer*slope1_denom) != EQUAL) { @@ -1260,13 +697,13 @@ private: // The curves do not have a vertical slope at p. // Compare their second-order derivative by x: derive_by_x_at(p, 2, slope1_numer, slope1_denom); - arc.derive_by_x_at(p, 2, slope2_numer, slope2_denom); + xcv.derive_by_x_at(p, 2, slope2_numer, slope2_denom); } else { // Both curves have a vertical slope at p. // Compare their second-order derivative by y: derive_by_y_at(p, 2, slope1_numer, slope1_denom); - arc.derive_by_y_at(p, 2, slope2_numer, slope2_denom); + xcv.derive_by_y_at(p, 2, slope2_numer, slope2_denom); } if (CGAL::compare(slope1_numer*slope2_denom, @@ -1287,28 +724,28 @@ private: */ template std::ostream& operator<<(std::ostream& os, - const _Conic_x_monotone_arc_2& arc) + const _Conic_x_monotone_arc_2& xcv) { // Output the supporting conic curve. - os << "{" << CGAL::to_double(arc.r()) << "*x^2 + " - << CGAL::to_double(arc.s()) << "*y^2 + " - << CGAL::to_double(arc.t()) << "*xy + " - << CGAL::to_double(arc.u()) << "*x + " - << CGAL::to_double(arc.v()) << "*y + " - << CGAL::to_double(arc.w()) << "}"; + os << "{" << CGAL::to_double(xcv.r()) << "*x^2 + " + << CGAL::to_double(xcv.s()) << "*y^2 + " + << CGAL::to_double(xcv.t()) << "*xy + " + << CGAL::to_double(xcv.u()) << "*x + " + << CGAL::to_double(xcv.v()) << "*y + " + << CGAL::to_double(xcv.w()) << "}"; // Output the endpoints. - os << " : (" << CGAL::to_double(arc.source().x()) << "," - << CGAL::to_double(arc.source().y()) << ") "; + os << " : (" << CGAL::to_double(xcv.source().x()) << "," + << CGAL::to_double(xcv.source().y()) << ") "; - if (arc.orientation() == CLOCKWISE) os << "--cw-->"; - else if (arc.orientation() == COUNTERCLOCKWISE) os << "--ccw-->"; + if (xcv.orientation() == CLOCKWISE) os << "--cw-->"; + else if (xcv.orientation() == COUNTERCLOCKWISE) os << "--ccw-->"; else os << "--l-->"; - os << " (" << CGAL::to_double(arc.target().x()) << "," - << CGAL::to_double(arc.target().y()) << ")"; + os << " (" << CGAL::to_double(xcv.target().x()) << "," + << CGAL::to_double(xcv.target().y()) << ")"; - return (os); + return os; } } //namespace CGAL