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 b32eeda58bb..b36d32fbec3 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 @@ -674,14 +674,12 @@ public: // Set the arc to be the full conic (and compute the orientation). Rational rat_coeffs[6]; - rat_coeffs[0] = r; rat_coeffs[1] = s; rat_coeffs[2] = t; rat_coeffs[3] = u; rat_coeffs[4] = v; rat_coeffs[5] = w; - set_full(rat_coeffs, true); } @@ -698,7 +696,7 @@ public: const Orientation& orient, const Point_2& source, const Point_2& target) const { // Make sure that the source and the taget are not the same. - const auto* alg_kernel = m_traits.m_alg_kernel; + const auto alg_kernel = m_traits.m_alg_kernel; CGAL_precondition(alg_kernel->compare_xy_2_object()(source, target) != EQUAL); // Set the arc properties (no need to compute the orientation). @@ -728,8 +726,8 @@ public: const Rational& x3 = p3.x(); const Rational& y3 = p3.y(); - const auto* nt_traits = m_traits.m_nt_traits; - const auto* alg_kernel = m_traits.m_alg_kernel; + const auto nt_traits = m_traits.m_nt_traits; + const auto alg_kernel = m_traits.m_alg_kernel; auto source = Point_2(nt_traits->convert(x1), nt_traits->convert(y1)); auto target = Point_2(nt_traits->convert(x3), nt_traits->convert(y3)); arc.set_enpoints(source, target); @@ -829,7 +827,7 @@ public: const Rational& x5 = p5.x(); const Rational& y5 = p5.y(); - const auto* nt_traits = m_traits.m_nt_traits; + const auto nt_traits = m_traits.m_nt_traits; auto source = Point_2(nt_traits->convert(x1), nt_traits->convert(y1)); auto target = Point_2(nt_traits->convert(x5), nt_traits->convert(y5)); arc.set_endpoints(source, target); @@ -919,7 +917,7 @@ public: rat_coeffs[4] = v; rat_coeffs[5] = w; - const auto* nt_traits = m_traits.m_nt_traits; + const auto nt_traits = m_traits.m_nt_traits; Integer base_coeffs[6]; nt_traits->convert_coefficients(rat_coeffs, rat_coeffs + 6, base_coeffs); @@ -1053,8 +1051,8 @@ public: * \pre p and q must not be the same. * \return A segment connecting p and q. */ - Curve_2 operator()(const Point_2& p, const Point_2& q) { - const auto* alg_kernel = m_traits.m_alg_kernel; + Curve_2 operator()(const Point_2& p, const Point_2& q) const { + const auto alg_kernel = m_traits.m_alg_kernel; CGAL_precondition(alg_kernel->compare_xy_2_object()(p, q) != EQUAL); Curve_2 arc; @@ -1080,6 +1078,96 @@ public: arc.set_extra_data(extra_data); return arc; } + + /*! Construct a conic arc from a given line segment. + * \param seg The line segment with rational endpoints. + */ + Curve_2 operator()(const Rat_segment_2& seg) const { + Curve_2 arc; + arc.set_orientation(COLLINEAR); + + // Set the source and target. + const auto rat_kernel = m_traits.m_rat_kernel; + Rat_point_2 source = rat_kernel->construct_vertex_2_object()(seg, 0); + Rat_point_2 target = rat_kernel->construct_vertex_2_object()(seg, 1); + const Rational& x1 = source.x(); + const Rational& y1 = source.y(); + const Rational& x2 = target.x(); + const Rational& y2 = target.y(); + + const auto nt_traits = m_traits.m_nt_traits; + arc.set_source(Point_2(nt_traits->convert(x1), nt_traits->convert(y1))); + arc.set_target(Point_2(nt_traits->convert(x2), nt_traits->convert(y2))); + + // Make sure that the source and the taget are not the same. + CGAL_precondition(rat_kernel->compare_xy_2_object()(source, target) != + EQUAL); + + // 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). + const Rational zero(0); + const Rational one(1); + Rational rat_coeffs[6]; + + rat_coeffs[0] = zero; + rat_coeffs[1] = zero; + rat_coeffs[2] = zero; + + if (rat_kernel->compare_x_2_object()(source, target) == EQUAL) { + // The supporting conic is a vertical line, of the form x = CONST. + rat_coeffs[3] = one; + rat_coeffs[4] = zero; + rat_coeffs[5] = -x1; + } + else { + // The supporting line is A*x + B*y + C = 0, where: + // + // A = y2 - y1, B = x1 - x2, C = x2*y1 - x1*y2 + // + rat_coeffs[3] = y2 - y1; + rat_coeffs[4] = x1 - x2; + rat_coeffs[5] = x2*y1 - x1*y2; + } + + // Set the arc properties (no need to compute the orientation). + m_traits.set(arc, rat_coeffs); + return arc; + } + + /*! Construct a conic arc that is a full circle. + * \param circ The circle with rational center and rational squared radius. + */ + Curve_2 operator()(const Rat_circle_2& circ) const { + Curve_2 arc; + arc.set_orientation(CLOCKWISE); + + // Get the circle properties. + const auto rat_kernel = m_traits.m_rat_kernel; + Rat_point_2 center = rat_kernel->construct_center_2_object()(circ); + Rational x0 = center.x(); + Rational y0 = center.y(); + Rational R_sqr = rat_kernel->compute_squared_radius_2_object()(circ); + + // Produce the correponding conic: if the circle center is (x0,y0) + // and its squared radius is R^2, that its equation is: + // x^2 + y^2 - 2*x0*x - 2*y0*y + (x0^2 + y0^2 - R^2) = 0 + // Note that this equation describes a curve with a negative (clockwise) + // orientation. + const Rational zero(0); + const Rational one(1); + const Rational minus_two(-2); + Rational rat_coeffs[6]; + rat_coeffs[0] = one; + rat_coeffs[1] = one; + rat_coeffs[2] = zero; + rat_coeffs[3] = minus_two*x0; + rat_coeffs[4] = minus_two*y0; + rat_coeffs[5] = x0*x0 + y0*y0 - R_sqr; + + // Set the arc to be the full conic (no need to compute the orientation). + m_traits.set_full(arc.rat_coeffs, false); + return arc; + } }; /*! Obtain a Construct_curve_2 functor object. */ @@ -1178,7 +1266,7 @@ public: * \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(Curve_2& arc, const Rational* rat_coeffs) { + void set(Curve_2& arc, const Rational* rat_coeffs) const { arc.set_flag(Curve_2::IS_VALID); // Convert the coefficients vector to an equivalent vector of integer @@ -1228,7 +1316,8 @@ public: if (arc.orientation() == COLLINEAR) { // Make sure the midpoint is on the line pair (thus making sure that // the two points are not taken from different lines). - Point_2 p_mid = m_alg_kernel->compare_xy_2_object()(source, target); + auto ctr_mid_point = m_alg_kernel->construct_midpoint_2_object(); + Point_2 p_mid = ctr_mid_point(source, target); // if (! is_on_supporting_conic(arc, p_mid)) if (CGAL::sign((m_nt_traits->convert(r)*p_mid.x() + m_nt_traits->convert(t)*p_mid.y() + @@ -1275,12 +1364,12 @@ public: 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()(source, target); + Point_2 p_mid = + m_alg_kernel->construct_midpoint_2_object()(source, 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); + bool finite_at_x = (points_at_x(arc, p_mid, ps) > 0); + bool finite_at_y = (points_at_y(arc, p_mid, ps) > 0); if (! finite_at_x && ! finite_at_y) { arc.reset_flags(); // inavlid arc @@ -1302,7 +1391,7 @@ public: * \param comp_orient Should we compute the orientation of the given curve. */ void set_full(Curve_2& arc, const Rational* rat_coeffs, - const bool& comp_orient) { + const bool& comp_orient) const { // Convert the coefficients vector to an equivalent vector of integer // coefficients. Integer int_coeffs[6]; @@ -1369,6 +1458,21 @@ public: 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 Curve_2& arc, const Point_2& p) const { + CGAL_precondition(! arc.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)); + } + /*! 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. @@ -1444,7 +1548,7 @@ public: /*! 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) { + void build_hyperbolic_arc_data(Curve_2& arc) const { // Let phi be the rotation angle of the conic from its canonic form. // We can write: // @@ -1537,6 +1641,132 @@ public: arc.sign_of_extra_data(target.x(), target.y())); arc.set_extra_data(extra_data); } + + /*! 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 Curve_2& arc, + 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 B = + m_nt_traits->convert(arc.t())*y + m_nt_traits->convert(arc.u()); + Algebraic C = + (m_nt_traits->convert(arc.s())*y + m_nt_traits->convert(arc.v()))*y + + m_nt_traits->convert(arc.w()); + + return solve_quadratic_equation(A, B, C, xs[0], xs[1]); + } + + /*! 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 Curve_2& arc, + 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 B = + m_nt_traits->convert(arc.t())*x + m_nt_traits->convert(arc.v()); + Algebraic C = + (m_nt_traits->convert(arc.r())*x + m_nt_traits->convert(arc.u()))*x + + m_nt_traits->convert(arc.w()); + + return solve_quadratic_equation(A, B, C, ys[0], ys[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). + * \return The number of disticnt solutions to the equation. + */ + int solve_quadratic_equation(const Algebraic& A, + const Algebraic& B, + const Algebraic& C, + Algebraic& x_minus, Algebraic& x_plus) const { + // Check if we actually have a linear equation. + if (CGAL::sign(A) == ZERO) { + if (CGAL::sign(B) == ZERO) return 0; + x_minus = x_plus = -C / B; + return 1; + } + + // Compute the discriminant and act according to its sign. + const Algebraic disc = B*B - 4*A*C; + Sign sign_disc = CGAL::sign(disc); + + // Check whether there are no real-valued solutions: + if (sign_disc == NEGATIVE) return 0; + else if (sign_disc == ZERO) { + // One distinct solution: + x_minus = x_plus = -B / (2*A); + return 1; + } + + // Compute the two distinct solutions: + Algebraic _2A = 2*A; + Nt_traits nt_traits; + Algebraic sqrt_disc = m_nt_traits->sqrt(disc); + + x_minus = -(B + sqrt_disc) / _2A; + x_plus = (sqrt_disc - B) / _2A; + return 2; + } + + /*! 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 Curve_2& arc, 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); + + // 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; + } + + // 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 Curve_2& arc, 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); + + // 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; + } + + // Return the number of points on the arc. + 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 de483e224ac..76a2f846916 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 @@ -61,27 +61,6 @@ public: typedef typename Alg_kernel::Point_2 Point_2; typedef _Conic_point_2 Conic_point_2; -protected: - Integer m_r; // - Integer m_s; // The coefficients of the supporting conic curve: - Integer m_t; // - Integer m_u; // - Integer m_v; // r*x^2 + s*y^2 + t*xy + u*x + v*y +w = 0 . - Integer m_w; // - - Orientation m_orient; // The orientation of the conic. - - // Bit masks for the m_info field. - enum { - IS_VALID = 0, - IS_FULL_CONIC, - LAST_INFO, - }; - - int m_info; // does the arc represent a full conic curve. - Conic_point_2 m_source; // the source of the arc (if not a full curve). - Conic_point_2 m_target; // the target of the arc (if not a full curve). - /*! \struct * 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 @@ -97,9 +76,29 @@ protected: Sign side; }; + // Bit masks for the m_info field. + enum { + IS_VALID = 0, + IS_FULL_CONIC, + LAST_INFO, + }; + +protected: + Integer m_r; // + Integer m_s; // The coefficients of the supporting conic curve: + Integer m_t; // + Integer m_u; // + Integer m_v; // r*x^2 + s*y^2 + t*xy + u*x + v*y +w = 0 . + Integer m_w; // + + Orientation m_orient; // The orientation of the conic. + int m_info; // does the arc represent a full conic curve. + Conic_point_2 m_source; // the source of the arc (if not a full curve). + Conic_point_2 m_target; // the target of the arc (if not a full curve). Extra_data* m_extra_data; // The extra data stored with the arc // (may be nullptr). +public: /// \name Flag manipulation functions. //@{ template @@ -864,7 +863,7 @@ public: /*! Obtain the extra data. */ - const Extra_data* extra_data() const { return m_extra_data(); } + const Extra_data* extra_data() const { return m_extra_data; } /*! Obtain a bounding box for the conic arc. * \return The bounding box. @@ -952,37 +951,13 @@ public: /*! Set the source point of the conic arc. * \param ps The new source point. - * \pre The arc is not a full conic curve. - * ps must lie on the supporting conic curve. */ - void set_source(const Point_2& ps) - { - CGAL_precondition(! is_full_conic()); - CGAL_precondition(is_on_supporting_conic(ps)); - CGAL_precondition(Alg_kernel().orientation_2_object() - (m_source, ps, m_target) == m_orient || - Alg_kernel().orientation_2_object() - (ps, m_source, m_target) == m_orient); - - m_source = ps; - } + void set_source(const Point_2& ps) { m_source = ps; } /*! Set the target point of the conic arc. * \param pt The new source point. - * \pre The arc is not a full conic curve. - * pt must lie on the supporting conic curve. */ - void set_target(const Point_2& pt) - { - CGAL_precondition(! is_full_conic()); - CGAL_precondition(is_on_supporting_conic(pt)); - CGAL_precondition(Alg_kernel().orientation_2_object() - (m_source, pt, m_target) == m_orient || - Alg_kernel().orientation_2_object() - (m_source, m_target, pt) == m_orient); - - m_target = pt; - } + void set_target(const Point_2& pt) { m_target = pt; } //@} @@ -1369,7 +1344,7 @@ private: } //@} -protected: +public: /// \name Auxiliary functions. //@{ @@ -1381,15 +1356,12 @@ protected: */ Sign sign_of_extra_data(const Algebraic& px, const Algebraic& py) const { CGAL_assertion(m_extra_data != nullptr); - if (m_extra_data == nullptr) return ZERO; - - Algebraic val = (m_extra_data->a*px + m_extra_data->b*py + - m_extra_data->c); - + Algebraic val = m_extra_data->a*px + m_extra_data->b*py + m_extra_data->c; return CGAL::sign(val); } +protected: /*! Check whether the given point lies on the supporting conic of the arc. * \param p The query point. * \return true if p lies on the supporting conic; (false) otherwise.