From 5094c75a19a81544de34da47b2ecb4d1c23b91eb Mon Sep 17 00:00:00 2001 From: Ron Wein Date: Tue, 27 Jun 2006 08:22:20 +0000 Subject: [PATCH] Updated the Bezier traits class. --- .../include/CGAL/Arr_Bezier_curve_traits_2.h | 30 ++-- .../CGAL/Arr_traits_2/Bezier_curve_2.h | 148 ++++++++++++++++-- .../CGAL/Arr_traits_2/Bezier_point_2.h | 59 ++++++- .../CGAL/Arr_traits_2/Bezier_x_monotone_2.h | 32 ++-- 4 files changed, 217 insertions(+), 52 deletions(-) diff --git a/Arrangement_2/include/CGAL/Arr_Bezier_curve_traits_2.h b/Arrangement_2/include/CGAL/Arr_Bezier_curve_traits_2.h index 98fbe5d628f..a5c2adcc6a7 100644 --- a/Arrangement_2/include/CGAL/Arr_Bezier_curve_traits_2.h +++ b/Arrangement_2/include/CGAL/Arr_Bezier_curve_traits_2.h @@ -267,22 +267,15 @@ public: } // Compare the slopes of the two arcs. - Comparison_result res; - unsigned int mult; - - res = cv1.compare_slopes (cv2, p, mult); + Comparison_result res = cv1.compare_slopes (cv2, p); + CGAL_assertion (res != EQUAL); + // TODO: Compare slopes may return EQUAL in case of tangency. + // We should use other methods. - // The comparison result is to the right of p. In case the multiplicity - // of the intersection point p is odd, reverse this result. - if (mult % 2 == 1) - { - if (res == SMALLER) - res = LARGER; - else if (res == LARGER) - res = SMALLER; - } - - return (res); + // The comparison result is to the right of p, so if the slopes are + // not equals this means that p is a simple intersection point and + // we therefore have to reverse the result to the left of p. + return ((res == LARGER) ? SMALLER : LARGER); } }; @@ -327,9 +320,12 @@ public: // Compare the slopes of the two arcs to determine thir relative // position immediately to the right of p. - unsigned int mult; + Comparison_result res = cv1.compare_slopes (cv2, p); + CGAL_assertion (res != EQUAL); + // TODO: Compare slopes may return EQUAL in case of tangency. + // We should use other methods. - return (cv1.compare_slopes (cv2, p, mult)); + return (res); } }; diff --git a/Arrangement_2/include/CGAL/Arr_traits_2/Bezier_curve_2.h b/Arrangement_2/include/CGAL/Arr_traits_2/Bezier_curve_2.h index 1a411dd532f..837562ef68f 100644 --- a/Arrangement_2/include/CGAL/Arr_traits_2/Bezier_curve_2.h +++ b/Arrangement_2/include/CGAL/Arr_traits_2/Bezier_curve_2.h @@ -84,6 +84,10 @@ private: Polynomial _polyY; // The polynomial for y. Integer _normY; // Normalizing factor for y. + typedef std::vector Control_point_vec; + + Control_point_vec _ctrl_pts; // The control points. + public: /*! Default constructor. */ @@ -110,6 +114,8 @@ public: CGAL_precondition_msg (n > 0, "There must be at least 2 control points."); + _ctrl_pts.resize (n + 1); + for (j = 0; j <= n; j++) coeffsX[j] = coeffsY[j] = rat_zero; @@ -129,6 +135,9 @@ public: for (k = 0; pts_begin != pts_end; ++pts_begin, k++) { + // Store the current control point and obtain its coordinates. + _ctrl_pts[k] = *pts_begin; + px = pts_begin->x(); py = pts_begin->y(); @@ -180,7 +189,7 @@ public: _polyY, _normY); delete[] coeffsY; - // TODO: Normalize the polynomials by GCD computation (?) + // TODO: Should we normalize the polynomials by GCD computation (?) } private: @@ -229,6 +238,8 @@ private: Nt_traits> Bcv_rep; typedef Handle_for Bcv_handle; + typedef typename Bcv_rep::Control_point_vec Control_pt_vec; + public: typedef typename Bcv_rep::Rat_point_2 Rat_point_2; @@ -239,6 +250,8 @@ public: typedef typename Bcv_rep::Algebraic Algebraic; typedef typename Bcv_rep::Polynomial Polynomial; + typedef typename Control_pt_vec::const_iterator Control_point_iterator; + public: /*! @@ -266,6 +279,8 @@ public: Bcv_handle (Bcv_rep (pts_begin, pts_end)) {} + // TODO: Assignment operator??? (Also for Bezier_curve_2) + /*! * Get the polynomial for the x-coordinates of the curve. */ @@ -298,6 +313,42 @@ public: return (_rep()._normY); } + /*! + * Get the number of control points inducing the Bezier curve. + */ + unsigned int number_of_control_points () const + { + return (_rep()._ctrl_pts.size()); + } + + /*! + * Get the i'th control point. + * \pre i must be between 0 and n - 1, where n is the number of control + * points. + */ + const Rat_point_2& control_point (unsigned int i) const + { + CGAL_precondition (i < number_of_control_points()); + + return ((_rep()._ctrl_pts)[i]); + } + + /*! + * Get an interator for the first control point. + */ + Control_point_iterator control_points_begin () const + { + return (_rep()._ctrl_pts.begin()); + } + + /*! + * Get a past-the-end interator for control points. + */ + Control_point_iterator control_points_end () const + { + return (_rep()._ctrl_pts.end()); + } + /*! * Check if both curve handles refer to the same object. */ @@ -331,25 +382,68 @@ public: } /*! - * Compute a point of the Bezier curve given a t-value. + * Compute a point of the Bezier curve given a rational t-value. + * \param t The given t-value. + * \pre t must be between 0 and 1. + */ + Rat_point_2 operator() (const Rational& t) const + { + CGAL_precondition (CGAL::sign (t) != NEGATIVE && + CGAL::compare (t, Rational(1)) != LARGER); + + // Compute the x and y coordinates. + const Rational x = nt_traits.evaluate_at (_rep()._polyX, t) / + Rational (_rep()._normX, 1); + const Rational y = nt_traits.evaluate_at (_rep()._polyY, t) / + Rational (_rep()._normY, 1); + + return (Rat_point_2 (x, y)); + } + + /*! + * Compute a point of the Bezier curve given an algebraic t-value. * \param t The given t-value. * \pre t must be between 0 and 1. */ Alg_point_2 operator() (const Algebraic& t) const { - CGAL_precondition (CGAL::compare (t, Algebraic(0)) != SMALLER && - CGAL::compare (t, Algebraic(1)) != LARGER); + // Check for extermal t values (either 0 or 1). + Nt_traits nt_traits; + const CGAL::Sign sign_t = CGAL::sign (t); - // Compute the x and y coordinates. - Nt_traits nt_traits; - const Algebraic x = nt_traits.evaluate_at (_rep()._polyX, t) / - nt_traits.convert (_rep()._normX); - const Algebraic y = nt_traits.evaluate_at (_rep()._polyY, t) / - nt_traits.convert (_rep()._normY); + CGAL_precondition (sign_t != NEGATIVE); - return Alg_point_2 (x, y); + if (sign_t == ZERO) + { + // Is t is 0, simply return the first control point. + const Rat_point_2& p_0 = _rep()._ctrl_pts[0]; + + return (Alg_point_2 (nt_traits.convert (p_0.x()), + nt_traits.convert (p_0.y()))); + } + + Comparison_result res = CGAL::compare (t, Algebraic(1)); + + CGAL_precondition (res != LARGER); + + if (res == EQUAL) + { + // Is t is 0, simply return the first control point. + const Rat_point_2& p_n = _rep()._ctrl_pts[_rep()._ctrl_pts.size() - 1]; + + return (Alg_point_2 (nt_traits.convert (p_n.x()), + nt_traits.convert (p_n.y()))); + } + + // The t-value is between 0 and 1: Compute the x and y coordinates. + const Algebraic x = nt_traits.evaluate_at (_rep()._polyX, t) / + nt_traits.convert (_rep()._normX); + const Algebraic y = nt_traits.evaluate_at (_rep()._polyY, t) / + nt_traits.convert (_rep()._normY); + + return (Alg_point_2 (x, y)); } - + /*! * Compute the points with vertical tangents on the curve. The function * actually returns t-values such that the tangent at (*this)(t) is vertical. @@ -366,6 +460,7 @@ public: const Algebraic one = Algebraic(1); nt_traits.compute_polynomial_roots (polyX_der, std::back_inserter(t_vals)); + // TODO: Compute polynomial roots in interval [0,1] // Take only t-values strictly between 0 and 1. Note that we use the // fact that the list of roots we obtain is sorted in ascending order. @@ -402,6 +497,7 @@ public: const Algebraic one = Algebraic(1); nt_traits.compute_polynomial_roots (polyY_der, std::back_inserter(t_vals)); + // TODO: compute roots in [0, 1]. // Take only t-values strictly between 0 and 1. Note that we use the // fact that the list of roots we obtain is sorted in ascending order. @@ -481,6 +577,7 @@ public: const Algebraic one = Algebraic(1); nt_traits.compute_polynomial_roots (res, std::back_inserter(s_vals)); + // TODO: Compute roots only in the interval [0,1]. typename std::list::iterator s_iter = s_vals.begin(); @@ -551,11 +648,17 @@ public: // Sample the approximated curve. const int n = (n_samples >= 2) ? n_samples : 2; const double delta_t = (t_end - t_start) / (n - 1); + double x, y; double t; for (k = 0; k < n; k++) { t = t_start + k * delta_t; + x = _evaluate_at (coeffsX, degX, t) / normX; + y = _evaluate_at (coeffsY, degY, t) / normY; + + *oi = std::make_pair (x, y); + ++oi; } return (oi); @@ -662,7 +765,6 @@ private: } } - // Zero the whole i'th column of the following rows. for (k = i + 1; k < dim; k++) { @@ -737,6 +839,26 @@ private: return (os); } + + /*! + * Evaluate a polynomial with double-precision coefficient at a given point. + * \param coeffs The coefficients. + * \param deg The degree of the polynomial. + * \param t The value to evaluate at. + * \return The value of the polynomial at t. + */ + double _evaluate_at (const std::vector& coeffs, + int deg, const double& t) const + { + // Use Horner's rule to evaluate the polynomial at t. + double val = coeffs[deg]; + int k; + + for (k = deg - 1; k >= 0; k--) + val = val*t + coeffs[k]; + + return (val); + } }; /*! diff --git a/Arrangement_2/include/CGAL/Arr_traits_2/Bezier_point_2.h b/Arrangement_2/include/CGAL/Arr_traits_2/Bezier_point_2.h index f561a806ebc..38a890928bf 100644 --- a/Arrangement_2/include/CGAL/Arr_traits_2/Bezier_point_2.h +++ b/Arrangement_2/include/CGAL/Arr_traits_2/Bezier_point_2.h @@ -53,7 +53,9 @@ public: typedef Alg_kernel_ Alg_kernel; typedef Nt_traits_ Nt_traits; + typedef typename Rat_kernel::Point_2 Rat_point_2; typedef typename Alg_kernel::Point_2 Alg_point_2; + typedef typename Nt_traits::Rational Rational; typedef typename Nt_traits::Algebraic Algebraic; private: @@ -85,6 +87,16 @@ public: _y (y) {} + /*! + * Constructor from a point with rational coordinates. + */ + _Bezier_point_2_rep (const Rat_point_2& p) + { + Nt_traits nt_traits; + _x = nt_traits.convert (p.x()); + _y = nt_traits.convert (p.y()); + } + /*! * Constructor from a point with algebraic coordinates. */ @@ -94,14 +106,28 @@ public: {} /*! - * Constructor given an originating curve and a t0 value. + * Constructor given an originating curve and a rational t0 value. + * \pre t0 must be between 0 and 1. + */ + _Bezier_point_2_rep (const Curve_2& B, const Rational& t0) + { + // Set the point coordinates. + Nt_traits nt_traits; + const Rat_point_2 p = B(t0); + + _x = nt_traits.convert (p.x()); + _y = nt_traits.convert (p.y()); + + // Create the originator pair . + _origs.push_back (Originator (B, nt_traits.convert (t0))); + } + + /*! + * Constructor given an originating curve and an algebraic t0 value. * \pre t0 must be between 0 and 1. */ _Bezier_point_2_rep (const Curve_2& B, const Algebraic& t0) { - CGAL_precondition (CGAL::sign (t0) != NEGATIVE && - CGAL::compare (t0, Algebraic(1)) != LARGER); - // Set the point coordinates. const Alg_point_2 p = B(t0); @@ -138,10 +164,12 @@ private: public: + typedef typename Bpt_rep::Rat_point_2 Rat_point_2; typedef typename Bpt_rep::Alg_point_2 Alg_point_2; + typedef typename Bpt_rep::Rational Rational; typedef typename Bpt_rep::Algebraic Algebraic; typedef typename Bpt_rep::Curve_2 Curve_2; - typedef typename Bpt_rep::Orig_iter Originator_const_iterator; + typedef typename Bpt_rep::Orig_iter Originator_iterator; /*! * Default constructor. @@ -164,6 +192,13 @@ public: Bpt_handle (Bpt_rep (x, y)) {} + /*! + * Constructor from a point with rational coordinates. + */ + _Bezier_point_2 (const Rat_point_2& p) : + Bpt_handle (Bpt_rep (p)) + {} + /*! * Constructor from a point with algebraic coordinates. */ @@ -172,7 +207,15 @@ public: {} /*! - * Constructor given an originating curve and a t0 value. + * Constructor given an originating curve and a rational t0 value. + * \pre t0 must be between 0 and 1. + */ + _Bezier_point_2 (const Curve_2& B, const Rational& t0) : + Bpt_handle (Bpt_rep (B, t0)) + {} + + /*! + * Constructor given an originating curve and an algebraic t0 value. * \pre t0 must be between 0 and 1. */ _Bezier_point_2 (const Curve_2& B, const Algebraic& t0) : @@ -259,12 +302,12 @@ public: /*! * Get the range of originators (pair of ). */ - Originator_const_iterator originators_begin () const + Originator_iterator originators_begin () const { return (_rep()._origs.begin()); } - Originator_const_iterator originators_end () const + Originator_iterator originators_end () const { return (_rep()._origs.end()); } diff --git a/Arrangement_2/include/CGAL/Arr_traits_2/Bezier_x_monotone_2.h b/Arrangement_2/include/CGAL/Arr_traits_2/Bezier_x_monotone_2.h index 052d9e9e481..1d11dea9b35 100644 --- a/Arrangement_2/include/CGAL/Arr_traits_2/Bezier_x_monotone_2.h +++ b/Arrangement_2/include/CGAL/Arr_traits_2/Bezier_x_monotone_2.h @@ -194,7 +194,10 @@ public: const Comparison_result res2 = CGAL::compare (p.x(), _pt.x()); ); CGAL_precondition (res1 == EQUAL || res2 == EQUAL || res1 != res2); - // RWRW - A hueristic solution. Find a robust one! + // TODO: First of all, check intersections with the originating curves + // of p. This way we find if p is on our curve. + + // TODO: A hueristic solution. Find a robust one! double t_low, t_high; const double init_eps = 0.00000001; bool dir_match; @@ -271,7 +274,6 @@ public: * their given intersection point. * \param cv The other subcurve. * \param p The intersection point. - * \param mult Output: The mutiplicity of the intersection point. * \pre p lies of both subcurves. * \pre Neither of the subcurves is a vertical segment. * \return SMALLER if (*this) slope is less than cv's; @@ -279,8 +281,7 @@ public: * LARGER if (*this) slope is greater than cv's. */ Comparison_result compare_slopes (const Self& cv, - const Point_2& p, - unsigned int& mult) const + const Point_2& p) const { if (_has_same_support (cv)) return (EQUAL); @@ -313,12 +314,7 @@ public: nt_traits.convert (cv._curve.y_norm())); // Compare the slopes. - Comparison_result res = CGAL::compare (slope1, slope2); - - CGAL_assertion (res != EQUAL); // RWRW: what should we do here? - - mult = 1; - return (res); + return (CGAL::compare (slope1, slope2)); } /*! @@ -345,7 +341,7 @@ public: OutputIterator intersect (const Self& cv, OutputIterator oi) const { - // RWRW: Handle overlapping curves ... + // TODO: Handle overlapping curves ... // Let B_1 be the supporting curve of (*this) and B_2 be the supporting // curve of cv. We compute s-values and t-values such that B_1(s) and @@ -369,10 +365,11 @@ public: // If the point lies in the range of (*this), find a t-value for cv // that matches this point. + // TODO: Be more carful here ... Point_2 p (_curve, *s_it); const double px = CGAL::to_double (p.x()); const double py = CGAL::to_double (p.y()); - unsigned int mult = 1; + unsigned int mult = 0; double sqr_dist; double best_sqr_dist = 0; @@ -459,7 +456,9 @@ public: */ bool can_merge_with (const Self& cv) const { - return (_has_same_support (cv) && + // Note that we only allow merging subcurves of the same originating + // Bezier curve (overlapping curves will not do in this case). + return (_curve.is_same (cv._curve) && (right().equals (cv.left()) || left().equals (cv.right()))); return (false); @@ -473,7 +472,7 @@ public: */ Self merge (const Self& cv) const { - CGAL_precondition (_has_same_support (cv)); + CGAL_precondition (_curve.is_same (cv._curve)); Self res = cv; @@ -517,6 +516,8 @@ public: */ Self flip () const { + // TODO: Is this "legal"? Should we touch the Bezier curve instead + // so that _trg > _src in all cases? Self cv = *this; cv._src = this->_trg; @@ -535,6 +536,9 @@ private: */ bool _has_same_support (const Self& cv) const { + // TODO: Sample (m + n + 1) rational points on one curve. + // If they are all on the other curve, the two are equal. + return (_curve.equals (cv._curve)); }