Subject: CGAL users discussion list
List archive
Re: [cgal-discuss] Arrangement_2: assertion violation in Arr_traits_2/Bezier_bounding_rational_traits.h
Chronological Thread
- From: Efraim Fogel <>
- To:
- Cc: Ron Wein <>
- Subject: Re: [cgal-discuss] Arrangement_2: assertion violation in Arr_traits_2/Bezier_bounding_rational_traits.h
- Date: Sun, 16 Sep 2007 18:01:53 +0200
Hi David,
Ron Wein, has provided a fix. Attached is a patch you can use to apply the fix. Apply it from the include/CGAL directory ('patch < Bezier.diff').
It suppose to affect the following files:
Arr_traits_2/Bezier_cache.h
Arr_traits_2/Bezier_curve_2.h
Arr_traits_2/Bezier_x_monotone_2.h
Arr_traits_2/Bezier_point_2.h
Arr_traits_2/Bezier_bounding_rational_traits.h
Arr_Bezier_curve_traits_2.h
The problem in test case #1 can simply be avoided if you comment out the assertion on line 247 in Bezier_bounding_rational_traits.h
The problem in test case #2 happens when there exists a vertical tangency point that occurs exactly at the middle of the curve. As you can see from the size of the patch, it was a bit more difficult to fix. We hope that we eliminated the bugs or at least reduced their number...
Naturally, the fix will make it into the next release.
David Keller wrote:
Hello,
I'm encountering problems with the Arrangement_2 package for Bezier-curves. From the documentation it seems the failure has not to be expected (please give me a hint if I missed something). So I think it would be a good decision to either fix it or give a more detailed report for which input the code is supposed to work.
There are two scenarios for which errors occur. Both can be reproduced by calling the example application in example/Arrangement_2/Bezier_curves.cpp.
TEST CASE 1:
-------------------
Input:
3
3 0 0 1 1 13/9 1
3 13/9 1 2 1 2 -7/8
3 2 -7/8 2 -2 0 3
Output:
B = {3 0 0 1 1 13/9 1}
B = {3 13/9 1 2 1 2 -7/8}
B = {3 2 -7/8 2 -2 0 3}
CGAL error: assertion violation!
Expr: intersection_pairs.size() == 1
File: /opt/CGAL-3.3.1/include/CGAL/Arr_traits_2/Bezier_bounding_rational_traits.h
Line: 247
Explanation:
terminate called after throwing an instance of 'CGAL::Assertion_exception'
what(): CGAL ERROR: assertion violation!
Expr: intersection_pairs.size() == 1
File: /opt/CGAL-3.3.1/include/CGAL/Arr_traits_2/Bezier_bounding_rational_traits.h
Line: 247
Aborted
TEST CASE 2:
-------------------
Input:
2
3 0 0 8 2 0 4
3 4 0 0 2 4 4
Output:
B = {3 0 0 8 2 0 4}
B = {3 4 0 0 2 4 4}
CGAL error: assertion violation!
Expr: CGAL::compare(right[0].x(), right[1].x()) != EQUAL
File: /opt/CGAL-3.3.1/include/CGAL/Arr_traits_2/Bezier_bounding_rational_traits.h
Line: 378
Explanation:
terminate called after throwing an instance of 'CGAL::Assertion_exception'
what(): CGAL ERROR: assertion violation!
Expr: CGAL::compare(right[0].x(), right[1].x()) != EQUAL
File: /opt/CGAL-3.3.1/include/CGAL/Arr_traits_2/Bezier_bounding_rational_traits.h
Line: 378
Aborted
Thanks and kind regards,
David Keller
--
____ _ ____ _
/_____/_) o /__________ __ //
(____ ( ( ( (_/ (_/-(-'_(/
_/
Index: Arr_traits_2/Bezier_cache.h
===================================================================
--- Arr_traits_2/Bezier_cache.h (revision 40349)
+++ Arr_traits_2/Bezier_cache.h (working copy)
@@ -222,10 +222,19 @@
const Polynomial& polyY_2, const Integer& normY_2,
bool& do_ovlp);
+ /*!
+ * Mark two curves as overlapping.
+ * \param id1 The ID of the first curve.
+ * \param id2 The ID of the second curve.
+ */
+ void mark_as_overlapping (const Curve_id& id1,
+ const Curve_id& id2);
+
private:
/*!
- * Compute all parameter values of (X_1(s), Y_1(s)) with (X_2(t), Y_2(t)).
+ * Compute all s-parameter values such that the curve (X_1(s), Y_1(s))
+ * intersects with (X_2(t), Y_2(t)) for some t.
* \param polyX_1 The coefficients of X_1.
* \param normX_1 The normalizing factor of X_1.
* \param polyY_1 The coefficients of Y_1.
@@ -247,6 +256,18 @@
Parameter_list& s_vals) const;
/*!
+ * Compute all s-parameter values of the self intersection of (X(s), Y(s))
+ * with (X(t), Y(t)) for some t.
+ * \param polyX The coefficients of X.
+ * \param polyY The coefficients of Y.
+ * \param s_vals Output: A list of s-values of the self-intersection points.
+ * Note that the values are not bounded to [0,1].
+ */
+ void _self_intersection_params (const Polynomial& polyX,
+ const Polynomial& polyY,
+ Parameter_list& s_vals) const;
+
+ /*!
* Compute the resultant of two bivariate polynomials in x and y with
* respect to y. The bivariate polynomials are given as vectors of
* polynomials, where bp1[i] is a coefficient of y^i, which is in turn a
@@ -302,7 +323,7 @@
const Polynomial& polyY_2, const Integer& normY_2,
bool& do_ovlp)
{
- CGAL_precondition (id1 < id2);
+ CGAL_precondition (id1 <= id2);
// Construct the pair of curve IDs, and try to find it in the map.
Curve_pair curve_pair (id1, id2);
@@ -318,6 +339,57 @@
// We need to compute the intersection-parameter pairs.
Intersection_info& info = intersect_map[curve_pair];
+ // Check if we have to compute a self intersection (a special case),
+ // or a regular intersection between two curves.
+ if (id1 == id2)
+ {
+ // Compute all parameter values that lead to a self intersection.
+ Parameter_list s_vals;
+
+ _self_intersection_params (polyX_1, polyY_1,
+ s_vals);
+
+ // Match pairs of parameter values that correspond to the same point.
+ // Note that we make sure that both parameter pairs are in the range [0, 1]
+ // (if not, the self-intersection point is imaginary).
+ typename Parameter_list::iterator s_it;
+ typename Parameter_list::iterator t_it;
+ const Algebraic one (1);
+ const Algebraic& denX = nt_traits.convert (normX_1);
+ const Algebraic& denY = nt_traits.convert (normY_1);
+ Point_list pts1;
+
+ for (s_it = s_vals.begin(); s_it != s_vals.end(); ++s_it)
+ {
+ if (CGAL::sign (*s_it) == NEGATIVE)
+ continue;
+
+ if (CGAL::compare (*s_it, one) == LARGER)
+ break;
+
+ const Algebraic& x = nt_traits.evaluate_at (polyX_1, *s_it);
+ const Algebraic& y = nt_traits.evaluate_at (polyY_1, *s_it);
+
+ for (t_it = s_it; t_it != s_vals.end(); ++t_it)
+ {
+ if (CGAL::compare (*t_it, one) == LARGER)
+ break;
+
+ if (CGAL::compare (nt_traits.evaluate_at (polyX_1, *t_it),
+ x) == EQUAL &&
+ CGAL::compare (nt_traits.evaluate_at (polyY_1, *t_it),
+ y) == EQUAL)
+ {
+ info.first.push_back (Intersection_point_2 (*s_it, *t_it,
+ x / denX, y / denY));
+ }
+ }
+ }
+
+ info.second = false;
+ return (info.first);
+ }
+
// Compute s-values and t-values such that (X_1(s), Y_1(s)) and
// (X_2(t), Y_2(t)) are the intersection points.
Parameter_list s_vals;
@@ -462,13 +534,42 @@
}
}
+ info.second = false;
return (info.first);
}
// ---------------------------------------------------------------------------
-// Compute all parameter values of (X_1(s), Y_1(s)) with (X_2(t), Y_2(t)).
+// Mark two curves as overlapping.
//
template<class NtTraits>
+void _Bezier_cache<NtTraits>::mark_as_overlapping (const Curve_id& id1,
+ const Curve_id& id2)
+{
+ CGAL_precondition (id1 < id2);
+
+ // Construct the pair of curve IDs, and try to find it in the map.
+ Curve_pair curve_pair (id1, id2);
+ Intersect_map_iterator map_iter = intersect_map.find (curve_pair);
+
+ if (map_iter != intersect_map.end())
+ {
+ // Found in the map: Make sure the curves are marked as overlapping.
+ CGAL_assertion (map_iter->second.second);
+ return;
+ }
+
+ // Add a new entry and mark the curves as overlapping.
+ Intersection_info& info = intersect_map[curve_pair];
+
+ info.second = true;
+ return;
+}
+
+// ---------------------------------------------------------------------------
+// Compute all s-parameter values of the intersection of (X_1(s), Y_1(s))
+// with (X_2(t), Y_2(t)) for some t.
+//
+template<class NtTraits>
bool _Bezier_cache<NtTraits>::_intersection_params
(const Polynomial& polyX_1, const Integer& normX_1,
const Polynomial& polyY_1, const Integer& normY_1,
@@ -531,6 +632,75 @@
}
// ---------------------------------------------------------------------------
+// Compute all s-parameter values of the self intersection of (X(s), Y(s))
+// with (X(t), Y(t)) for some t.
+//
+template<class NtTraits>
+void _Bezier_cache<NtTraits>::_self_intersection_params
+ (const Polynomial& polyX,
+ const Polynomial& polyY,
+ Parameter_list& s_vals) const
+{
+ // Clear the output parameter list.
+ if (! s_vals.empty())
+ s_vals.clear();
+
+ // We are looking for the solutions of the following system of bivariate
+ // polynomials in s and t:
+ // I: X(t) - X(s) / (t - s) = 0
+ // II: Y(t) - Y(s) / (t - s) = 0
+ //
+ Integer *coeffs;
+ int i, k;
+
+ // Consruct the bivariate polynomial that corresponds to Equation I.
+ // Note that we represent a bivariate polynomial as a vector of univariate
+ // polynomials, whose i'th entry corresponds to the coefficient of t^i,
+ // which is in turn a polynomial it s.
+ const int degX = nt_traits.degree (polyX);
+ std::vector<Polynomial> coeffsX_st (degX);
+
+ coeffs = new Integer [degX];
+
+ for (i = 0; i < degX; i++)
+ {
+ for (k = i + 1; k < degX; k++)
+ coeffs[k - i - 1] = nt_traits.get_coefficient (polyX, k);
+
+ coeffsX_st[i] = nt_traits.construct_polynomial (coeffs, degX - i - 1);
+ }
+
+ delete[] coeffs;
+
+ // Consruct the bivariate polynomial that corresponds to Equation II.
+ const int degY = nt_traits.degree (polyY);
+ std::vector<Polynomial> coeffsY_st (degY);
+
+ coeffs = new Integer [degY];
+
+ for (i = 0; i < degY; i++)
+ {
+ for (k = i + 1; k < degY; k++)
+ coeffs[k - i - 1] = nt_traits.get_coefficient (polyY, k);
+
+ coeffsY_st[i] = nt_traits.construct_polynomial (coeffs, degY - i - 1);
+ }
+
+ delete[] coeffs;
+
+ // Compute the resultant of the two bivariate polynomials and obtain
+ // a polynomial in s.
+ Polynomial res = _compute_resultant (coeffsX_st, coeffsY_st);
+
+ if (nt_traits.degree (res) < 0)
+ return;
+
+ // Compute the roots of the resultant polynomial.
+ nt_traits.compute_polynomial_roots (res, std::back_inserter (s_vals));
+ return;
+}
+
+// ---------------------------------------------------------------------------
// Compute the resultant of two bivariate polynomials.
//
template<class NtTraits>
@@ -698,4 +868,3 @@
CGAL_END_NAMESPACE
#endif
-
Index: Arr_traits_2/Bezier_curve_2.h
===================================================================
--- Arr_traits_2/Bezier_curve_2.h (revision 40349)
+++ Arr_traits_2/Bezier_curve_2.h (working copy)
@@ -28,8 +28,6 @@
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Handle_for.h>
#include <CGAL/Arr_traits_2/de_Casteljau_2.h>
-#include <CGAL/Sweep_line_2_algorithms.h>
-
#include <algorithm>
#include <deque>
#include <vector>
@@ -87,11 +85,14 @@
private:
- typedef std::deque<Rat_point_2> Control_point_vec;
+ typedef std::deque<Rat_point_2> Control_point_vec;
- Control_point_vec _ctrl_pts; // The control points (we prefer deque to
- // a vector, as it enables push_front()).
- Bbox_2 _bbox; // A bounding box for the curve.
+ Control_point_vec _ctrl_pts; /*!< The control points (we prefer deque
+ to a vector, as it enables
+ push_front()). */
+ Bbox_2 _bbox; /*!< A bounding box for the curve. */
+ bool _no_self_inter; /*!< Whether the curve surely has no
+ self intersections. */
/// \name Lazily-evaluated values of the polynomials B(t) = (X(t), Y(t)).
//@{
@@ -109,16 +110,18 @@
/*! Default constructor. */
_Bezier_curve_2_rep () :
+ _no_self_inter (true),
p_polyX(NULL),
p_normX(NULL),
p_polyY(NULL),
- p_normY(NULL)
+ p_normY(NULL)
{}
/*! Copy constructor (isn't really used). */
_Bezier_curve_2_rep (const _Bezier_curve_2_rep& other) :
_ctrl_pts(other._ctrl_pts),
_bbox(other._bbox),
+ _no_self_inter(other._no_self_inter),
p_polyX(NULL),
p_normX(NULL),
p_polyY(NULL),
@@ -138,9 +141,7 @@
* Constructor from a given range of control points.
* \param pts_begin An iterator pointing to the first point in the range.
* \param pts_end A past-the-end iterator for the range.
- * The value-type of the input iterators must be Rat_kernel::Point_2.
- * \pre There must be at least 2 control points, and the polyline defined
- * by these points may not be self intersecting.
+ * \pre The value-type of the input iterator must be Rat_kernel::Point_2.
*/
template <class InputIterator>
_Bezier_curve_2_rep (InputIterator pts_begin, InputIterator pts_end) :
@@ -149,25 +150,6 @@
p_polyY(NULL),
p_normY(NULL)
{
- // Make sure that the control polyline is not self-intersecting.
- CGAL_expensive_precondition_code (
- Rat_kernel ker;
- InputIterator pt_curr = pts_begin;
- InputIterator pt_next = pt_curr;
- std::list<typename Rat_kernel::Segment_2> segs;
-
- ++pt_next;
- while (pt_next != pts_end)
- {
- segs.push_back (ker.construct_segment_2_object() (*pt_curr, *pt_next));
- pt_curr = pt_next;
- ++pt_next;
- }
- );
- CGAL_expensive_precondition_msg
- (! do_curves_intersect (segs.begin(), segs.end()),
- "The control polyline may be be self-intersecting.");
-
// Copy the control points and compute their bounding box.
const int pts_size = std::distance (pts_begin, pts_end);
double x, y;
@@ -210,6 +192,12 @@
// Construct the bounding box.
_bbox = Bbox_2 (x_min, y_min, x_max, y_max);
+
+ // Use the bounding traits to determine whether the curve surely has
+ // not self intersections.
+ Bounding_traits bound_tr;
+
+ _no_self_inter = ! bound_tr.may_have_self_intersections (_ctrl_pts);
}
/*! Destructor. */
@@ -349,8 +337,6 @@
* \param pts_begin An iterator pointing to the first point in the range.
* \param pts_end A past-the-end iterator for the range.
* \pre The value-type of the input iterator must be Rat_kernel::Point_2.
- * \pre There must be at least 2 control points, and the polyline defined
- * by these points may not be self intersecting.
*/
template <class InputIterator>
_Bezier_curve_2 (InputIterator pts_begin, InputIterator pts_end) :
@@ -559,6 +545,16 @@
return (_rep()._bbox);
}
+ /*!
+ * Check if the curve contains not self intersections.
+ * Note that there may not be any self intersections even if the
+ * function returns true (but not vice versa).
+ */
+ bool has_no_self_intersections () const
+ {
+ return (_rep()._no_self_inter);
+ }
+
private:
// Get the representation.
@@ -595,7 +591,9 @@
int k;
for (k = 1; k <= deg; k++)
+ {
coeffs[k] = nt_traits.get_coefficient (poly, k) * denom;
+ }
coeffs[0] = nt_traits.get_coefficient (poly, 0) * denom -
numer * norm;
@@ -752,8 +750,8 @@
*p_polyY, *p_normY);
delete[] coeffsY;
- CGAL_assertion (nt_traits.degree (*p_polyX) > 0);
- CGAL_assertion (nt_traits.degree (*p_polyY) > 0);
+ CGAL_assertion (nt_traits.degree (*p_polyX) >= 0);
+ CGAL_assertion (nt_traits.degree (*p_polyY) >= 0);
return;
}
@@ -779,8 +777,8 @@
for (i = 2; i <= k; i++)
k_fact *= Integer (i);
-
- return (reduced_fact / (j_fact * k_fact));
+
+ return (CGAL::div (reduced_fact, (j_fact * k_fact)));
}
// ---------------------------------------------------------------------------
@@ -826,9 +824,9 @@
Nt_traits nt_traits;
x = nt_traits.evaluate_at (_rep().x_polynomial(), t) /
- Rational (_rep().x_norm(), 1);
+ Rational (_rep().x_norm(), 1);
y = nt_traits.evaluate_at (_rep().y_polynomial(), t) /
- Rational (_rep().y_norm(), 1);
+ Rational (_rep().y_norm(), 1);
// Return the point.
return (Rat_point_2 (x, y));
@@ -870,11 +868,6 @@
nt_traits.convert (p_n.y())));
}
- // IDDO: Check if it is worth evaluating this using de Casteljau
- // (since the parameter value is algebraic).
- // IDDO: Check if there is a way to know that an Algebraic is a rational
- // number...
-
// The t-value is between 0 and 1: Compute the x and y coordinates.
const Algebraic x = nt_traits.evaluate_at (_rep().x_polynomial(), t) /
nt_traits.convert (_rep().x_norm());
@@ -922,7 +915,7 @@
for (t_iter = t_vals.begin(); t_iter != t_vals.end(); ++t_iter)
{
- const Alg_point_2& p2 = bc (*t_iter, false);
+ const Alg_point_2& p2 = bc (*t_iter);
if (CGAL::compare (y1, p2.y()) == CGAL::EQUAL)
{
Index: Arr_traits_2/Bezier_x_monotone_2.h
===================================================================
--- Arr_traits_2/Bezier_x_monotone_2.h (revision 40349)
+++ Arr_traits_2/Bezier_x_monotone_2.h (working copy)
@@ -106,6 +106,7 @@
bool operator() (const Intersection_point_2& ip1,
const Intersection_point_2& ip2) const
{
+ // Use an xy-lexicographic comparison.
return (ip1.first.compare_xy (ip2.first, *p_cache) == SMALLER);
}
};
@@ -113,8 +114,8 @@
// Type definition for the bounded intersection-point mapping.
typedef std::list<Point_2> Intersection_list;
- /*! \struct
- * An auxiliary functor for comparing pair of curve IDs.
+ /*! \struct Less_curve_pair
+ * An auxiliary functor for comparing a pair of curve IDs.
*/
struct Less_curve_pair
{
@@ -175,19 +176,23 @@
private:
// Data members:
- Curve_2 _curve; /*! The supporting Bezier curve. */
- Point_2 _ps; /*! The source point. */
- Point_2 _pt; /*! The target point. */
- bool _dir_right; /*! Is the subcurve directed right (or left). */
- bool _inc_to_right; /*! Does the parameter value increase when
- traversing the subcurve from left to
- right. */
- bool _is_vert; /*! Is the subcurve a vertical segment. */
+ Curve_2 _curve; /*!< The supporting Bezier curve. */
+ unsigned int _xid; /*!< The serial number of the basic x-monotone
+ subcurve of the Bezier curve. */
+ Point_2 _ps; /*!< The source point. */
+ Point_2 _pt; /*!< The target point. */
+ bool _dir_right; /*!< Is the subcurve directed right
+ (or left). */
+ bool _inc_to_right; /*!< Does the parameter value increase when
+ traversing the subcurve from left to
+ right. */
+ bool _is_vert; /*!< Is the subcurve a vertical segment. */
public:
/*! Default constructor. */
_Bezier_x_monotone_2 () :
+ _xid(0),
_dir_right (false),
_is_vert (false)
{}
@@ -195,12 +200,18 @@
/*!
* Constructor given two endpoints.
* \param B The supporting Bezier curve.
+ * \param xid The serial number of the x-monotone subcurve with respect to
+ * the parameter range of the Bezier curve.
+ * For example, if B is split to two x-monotone subcurves at t',
+ * the subcurve defined over [0, t'] has a serial number 1,
+ * and the other, defined over [t', 1] has a serial number 2.
* \param ps The source point.
* \param pt The target point.
* \param cache Caches the vertical tangency points and intersection points.
* \pre B should be an originator of both ps and pt.
+ * \pre xid is a non-zero serial number.
*/
- _Bezier_x_monotone_2 (const Curve_2& B,
+ _Bezier_x_monotone_2 (const Curve_2& B, unsigned int xid,
const Point_2& ps, const Point_2& pt,
Bezier_cache& cache);
@@ -213,6 +224,14 @@
}
/*!
+ * Get the x-monotone ID of the curve.
+ */
+ unsigned int xid () const
+ {
+ return (_xid);
+ }
+
+ /*!
* Get the source point.
*/
const Point_2& source () const
@@ -285,7 +304,7 @@
* \param cv The other subcurve.
* \param p The intersection point.
* \param cache Caches the vertical tangency points and intersection points.
- * \pre p is the common left endpoint of both subcurves.
+ * \pre p is incident to both subcurves, and both are defined to its right.
* \return SMALLER if (*this) lies below cv to the right of p;
* EQUAL in case of an overlap (should not happen);
* LARGER if (*this) lies above cv to the right of p.
@@ -300,7 +319,7 @@
* \param cv The other subcurve.
* \param p The intersection point.
* \param cache Caches the vertical tangency points and intersection points.
- * \pre p is the common right endpoint of both subcurves.
+ * \pre p is incident to both subcurves, and both are defined to its left.
* \return SMALLER if (*this) lies below cv to the right of p;
* EQUAL in case of an overlap (should not happen);
* LARGER if (*this) lies above cv to the right of p.
@@ -334,7 +353,7 @@
OutputIterator oi) const
{
// In case we have two x-monotone subcurves of the same Bezier curve,
- // the only intersections that may occur are common endpoints.
+ // check if they have a common left endpoint.
if (_curve.is_same (cv._curve))
{
if (left().is_same (cv.left()) || left().is_same (cv.right()))
@@ -342,16 +361,8 @@
*oi = CGAL::make_object (Intersection_point_2 (left(), 0));
++oi;
}
+ }
- if (right().is_same (cv.left()) || right().is_same (cv.right()))
- {
- *oi = CGAL::make_object (Intersection_point_2 (right(), 0));
- ++oi;
- }
-
- return (oi);
- }
-
// Compute the intersections of the two sucurves. Note that for caching
// purposes we always apply the _intersect() function on the subcurve whose
// curve ID is smaller.
@@ -359,8 +370,7 @@
Self ovlp_cv;
bool do_ovlp;
- CGAL_assertion (_curve.id() != cv._curve.id());
- if (_curve.id() < cv._curve.id())
+ if (_curve.id() <= cv._curve.id())
do_ovlp = _intersect (cv, inter_map, cache, ipts, ovlp_cv);
else
do_ovlp = cv._intersect (*this, inter_map, cache, ipts, ovlp_cv);
@@ -383,6 +393,18 @@
*oi = CGAL::make_object (*ip_it);
++oi;
}
+
+ // In case we have two x-monotone subcurves of the same Bezier curve,
+ // check if they have a common right endpoint.
+ if (_curve.is_same (cv._curve))
+ {
+ if (right().is_same (cv.left()) || right().is_same (cv.right()))
+ {
+ *oi = CGAL::make_object (Intersection_point_2 (right(), 0));
+ ++oi;
+ }
+ }
+
return (oi);
}
@@ -417,8 +439,8 @@
*/
Self flip () const
{
- // TODO: Is this "legal"? Should we touch the Bezier curve instead
- // so that _trg > _src in all cases?
+ // Note that we just swap the source and target of the original subcurve
+ // and do not touch the supporting Beizer curve.
Self cv = *this;
cv._ps = this->_pt;
@@ -448,6 +470,21 @@
bool _is_in_range (const Point_2& p, bool& is_certain) const;
/*!
+ * Given a point p that lies on the supporting Bezier curve (X(t), Y(t)),
+ * determine whether p lies within the t-range of the x-monotone subcurve.
+ * If so, the value t0 such that p = (X(t0), Y(t0)) is also computed.
+ * \param p The point, which lies on the supporting Bezier curve.
+ * \param cache Caches the vertical tangency points and intersection points.
+ * \param t0 Output: A parameter value such that p = (X(t0), Y(t0)).
+ * \param is_endpoint Output: Whether p equals on the of the endpoints.
+ * \return Whether p lies in the t-range of the subcurve.
+ */
+ bool _is_in_range (const Point_2& p,
+ Bezier_cache& cache,
+ Algebraic& t0,
+ bool& is_endpoint) const;
+
+ /*!
* Compute a y-coordinate of a point on the x-monotone subcurve with a
* given x-coordinate.
* \param x0 The given x-coodinate.
@@ -462,6 +499,7 @@
* their given intersection point.
* \param cv The other subcurve.
* \param p The intersection point.
+ * \param cache Caches the vertical tangency points and intersection points.
* \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;
@@ -469,7 +507,8 @@
* LARGER if (*this) slope is greater than cv's.
*/
Comparison_result _compare_slopes (const Self& cv,
- const Point_2& p) const;
+ const Point_2& p,
+ Bezier_cache& cache) const;
/*!
* Get the range of t-value over which the subcurve is defined.
@@ -489,11 +528,11 @@
* \param cache Caches the vertical tangency points and intersection points.
* \pre The x-coordinate of the right endpoint of *this is smaller than
* (or equal to) the x-coordinate of the right endpoint of cv.
- * \pre p is the common left endpoint of both subcurves.
+ * \pre p is a common point of both subcurves.
* \pre Neither of the subcurves is a vertical segment.
- * \return SMALLER if (*this) lies below cv to the right of p;
+ * \return SMALLER if (*this) lies below cv next to p;
* EQUAL in case of an overlap (should not happen);
- * LARGER if (*this) lies above cv to the right of p.
+ * LARGER if (*this) lies above cv next to p.
*/
Comparison_result _compare_to_side (const Self& cv,
const Point_2& p,
@@ -501,14 +540,29 @@
Bezier_cache& cache) const;
/*!
- * Approximate the intersection points between the two given curves.
- * \param B1 The first Bezier curve.
- * \param B2 The second Bezier curve.
- * \param inter_pts Output: An output list of intersection points.
+ * Clip the control polygon of the supporting Bezier curve such that it
+ * fits the current x-monotone subcurve.
+ * \param ctrl Output: The clipped control polygon.
+ * \param t_min Output: The minimal t-value of the clipped curve.
+ * \param t_max Output: The maximal t-value of the clipped curve.
+ */
+ void _clip_control_polygon
+ (typename Bounding_traits::Control_points& ctrl,
+ typename Bounding_traits::NT& t_min,
+ typename Bounding_traits::NT& t_max) const;
+
+ /*!
+ * Approximate the intersection points between the supporting Bezier curves
+ * of the given x-monotone curves.
+ * \param cv The x-monotone curve we intersect.
+ * \param inter_pts Output: An output list of intersection points between
+ * the supporting Bezier curves.
+ * In case (*this) and cv are subcurve of the same
+ * Bezier curve, only the self-intersection points
+ * between the subcurves are approximated.
* \return Whether all intersection points where successfully approximated.
*/
- bool _approximate_intersection_points (const Curve_2& B1,
- const Curve_2& B2,
+ bool _approximate_intersection_points (const Self& cv,
std::list<Point_2>& inter_pts) const;
/*!
@@ -531,9 +585,9 @@
* the x-monotone curve.
* \param p The point.
* \param force_exact Sould we force an exact result.
- * \return SMALLER if the point is below the arc;
- * LARGER if the point is above the arc;
- * EQUAL if p lies on the arc.
+ * \return SMALLER if the point is below the curve;
+ * LARGER if the point is above the curve;
+ * EQUAL if p lies on the curve.
*/
Comparison_result _exact_vertical_position (const Point_2& p,
bool force_exact) const;
@@ -551,7 +605,8 @@
Bounding_traits>& cv)
{
os << cv.supporting_curve()
- << " | " << cv.source()
+ << " [" << cv.xid()
+ << "] | " << cv.source()
<< " --> " << cv.target();
return (os);
@@ -562,31 +617,31 @@
//
template <class RatKer, class AlgKer, class NtTrt, class BndTrt>
_Bezier_x_monotone_2<RatKer, AlgKer, NtTrt, BndTrt>::_Bezier_x_monotone_2
- (const Curve_2& B,
+ (const Curve_2& B, unsigned int xid,
const Point_2& ps, const Point_2& pt,
Bezier_cache& cache) :
_curve (B),
+ _xid (xid),
_ps (ps),
_pt (pt),
_is_vert (false)
{
+ CGAL_precondition (xid > 0);
+
// Get the originators of the point that correspond to the curve B.
- Originator_iterator ps_org = ps.get_originator(B);
+ Originator_iterator ps_org = ps.get_originator (B, _xid);
CGAL_precondition (ps_org != ps.originators_end());
- Originator_iterator pt_org = pt.get_originator(B);
+ Originator_iterator pt_org = pt.get_originator (B, _xid);
CGAL_precondition (pt_org != pt.originators_end());
// Check if the subcurve is directed left or right.
const Comparison_result res = _ps.compare_x (_pt, cache);
- // Iddo: TODO - alternative check if the original curve is vertical,
- // a check that can work on intervals.
-
+
if (res == EQUAL)
{
// We have a vertical segment. Check if the source is below the target.
_is_vert = true;
- // Iddo: TODO change the check to use compare_xy (or point bbox).
_dir_right = (CGAL::compare (_ps.y(), _pt.y()) == SMALLER);
}
else
@@ -617,9 +672,9 @@
t_res = CGAL::compare (ps_org->parameter(), pt_org->parameter());
}
-
+
CGAL_precondition (t_res != EQUAL);
-
+
if (_dir_right)
_inc_to_right = (t_res == SMALLER);
else
@@ -630,14 +685,14 @@
// Get the approximate parameter range defining the curve.
//
template <class RatKer, class AlgKer, class NtTrt, class BndTrt>
- std::pair<double, double>
+std::pair<double, double>
_Bezier_x_monotone_2<RatKer, AlgKer, NtTrt, BndTrt>::parameter_range () const
{
// First try to use the approximate representation of the endpoints.
- Originator_iterator s_org = _ps.get_originator (_curve);
+ Originator_iterator s_org = _ps.get_originator (_curve, _xid);
CGAL_assertion (s_org != _ps.originators_end());
- Originator_iterator t_org = _pt.get_originator (_curve);
+ Originator_iterator t_org = _pt.get_originator (_curve, _xid);
CGAL_assertion (t_org != _pt.originators_end());
double t_src = (CGAL::to_double (s_org->point_bound().t_min) +
@@ -702,14 +757,14 @@
CGAL_precondition (res1 != res2);
// Check for the case when curve is an originator of the point.
- Originator_iterator p_org = p.get_originator(_curve);
+ Originator_iterator p_org = p.get_originator (_curve, _xid);
if (p_org != p.originators_end())
{
- Originator_iterator ps_org = _ps.get_originator(_curve);
+ Originator_iterator ps_org = _ps.get_originator (_curve, _xid);
CGAL_assertion (ps_org != _ps.originators_end());
- Originator_iterator pt_org = _pt.get_originator(_curve);
+ Originator_iterator pt_org = _pt.get_originator (_curve, _xid);
CGAL_assertion (pt_org != _pt.originators_end());
// Check if the point is in the parameter range of this subcurve.
@@ -736,15 +791,15 @@
// Call the vertical-position function that uses the bounding-boxes
// to evaluate the comparsion result.
- typename Bounding_traits::Control_point_vec cp;
+ typename Bounding_traits::Control_points cp;
std::copy (_curve.control_points_begin(), _curve.control_points_end(),
std::back_inserter(cp));
- Originator_iterator ps_org = _ps.get_originator(_curve);
+ Originator_iterator ps_org = _ps.get_originator (_curve, _xid);
CGAL_assertion (ps_org != _ps.originators_end());
- Originator_iterator pt_org = _pt.get_originator(_curve);
+ Originator_iterator pt_org = _pt.get_originator (_curve, _xid);
CGAL_assertion (pt_org != _pt.originators_end());
Comparison_result res_bound = EQUAL;
@@ -853,7 +908,7 @@
CGAL_assertion (p.originators_begin() != p.originators_end());
- // Iddo: If the point is a rational point (e.g., ray shooting)
+ // \todo If the point is a rational point (e.g., ray shooting)
// use comparison between Y(root_of(X0-X(t))) and Y0.
Originator org = *(p.originators_begin());
bool do_ovlp;
@@ -893,16 +948,16 @@
CGAL_assertion (_is_in_range (t, cache));
Point_2& pt = const_cast<Point_2&> (p);
- pt.add_originator (Originator (_curve, t));
+ pt.add_originator (Originator (_curve, _xid, t));
// The point p lies on the subcurve.
return (EQUAL);
}
}
- // We now that p is not on the subcurve. We therefore subdivide the curve
- // using exact rational arithmetic, until we reach a separation between the
- // curve and the point. (This case should be very rare.)
+ // We now know that p is not on the subcurve. We therefore subdivide the
+ // curve using exact rational arithmetic, until we reach a separation
+ // between the curve and the point. (This case should be very rare.)
// Note that we first try to work with inexact endpoint representation, and
// only if we fail we make the endpoints of the x-monotone curves exact.
if (! p.is_exact())
@@ -935,10 +990,10 @@
{
CGAL_precondition (p.compare_xy (right(), cache) != LARGER);
CGAL_precondition (p.compare_xy (cv.right(), cache) != LARGER);
-
+
if (this == &cv)
return (EQUAL);
-
+
// Make sure that p is incident to both curves (either equals the left
// endpoint or lies in the curve interior). Note that this is important to
// carry out these tests, as it assures us the eventually both curves are
@@ -979,58 +1034,60 @@
if (_curve.is_same (cv._curve))
{
- // Get the originator, and make sure that p is a vertical tangency
- // point of this originator.
- Originator_iterator org = p.get_originator(_curve);
-
- CGAL_assertion (org != p.originators_end());
- CGAL_assertion (_inc_to_right != cv._inc_to_right);
-
- if (! p.is_exact())
+ // Get the originator, and check whether p is a vertical tangency
+ // point of this originator (otherwise it is a self-intersection point,
+ // and we proceed as if it is a regular intersection point).
+ Originator_iterator org = p.get_originator(_curve, _xid);
+
+ CGAL_assertion (org != p.originators_end());
+
+ if (org->point_bound().type == Bez_point_bound::VERTICAL_TANGENCY_PT)
{
- CGAL_assertion (org->point_bound().point_type ==
- Bez_point_bound::VERTICAL_TANGENCY_PT);
+ CGAL_assertion (_inc_to_right != cv._inc_to_right);
- // Comparison based on the control polygon of the bounded vertical
- // tangency point, using the fact this polygon is y-monotone.
- const typename Bounding_traits::Control_point_vec& cp =
- org->point_bound().bounding_polyline;
-
- if (_inc_to_right)
+ if (! p.is_exact())
{
- return (CGAL::compare (cp.back().y(), cp.front().y()));
+ // Comparison based on the control polygon of the bounded vertical
+ // tangency point, using the fact this polygon is y-monotone.
+ const typename Bounding_traits::Control_points& cp =
+ org->point_bound().ctrl;
+
+ if (_inc_to_right)
+ {
+ return (CGAL::compare (cp.back().y(), cp.front().y()));
+ }
+ else
+ {
+ return (CGAL::compare (cp.front().y(), cp.back().y()));
+ }
}
- else
- {
- return (CGAL::compare (cp.front().y(), cp.back().y()));
- }
- }
- // Iddo: Handle rational points (using de Casteljau derivative)?
+ // Iddo: Handle rational points (using de Casteljau derivative)?
- // In this case we know that we have a vertical tangency at t0, so
- // X'(t0) = 0. We evaluate the sign of Y'(t0) in order to find the
- // vertical position of the two subcurves to the right of this point.
- CGAL_assertion (org->has_parameter());
+ // In this case we know that we have a vertical tangency at t0, so
+ // X'(t0) = 0. We evaluate the sign of Y'(t0) in order to find the
+ // vertical position of the two subcurves to the right of this point.
+ CGAL_assertion (org->has_parameter());
- const Algebraic& t0 = org->parameter();
- Polynomial polyY_der = nt_traits.derive (_curve.y_polynomial());
- const CGAL::Sign sign_der =
- CGAL::sign (nt_traits.evaluate_at (polyY_der, t0));
+ const Algebraic& t0 = org->parameter();
+ Polynomial polyY_der = nt_traits.derive (_curve.y_polynomial());
+ const CGAL::Sign sign_der =
+ CGAL::sign (nt_traits.evaluate_at (polyY_der, t0));
- CGAL_assertion (sign_der != CGAL::ZERO);
+ CGAL_assertion (sign_der != CGAL::ZERO);
- if (_inc_to_right)
- return ((sign_der == CGAL::POSITIVE) ? LARGER : SMALLER);
- else
- return ((sign_der == CGAL::NEGATIVE) ? LARGER : SMALLER);
+ if (_inc_to_right)
+ return ((sign_der == CGAL::POSITIVE) ? LARGER : SMALLER);
+ else
+ return ((sign_der == CGAL::NEGATIVE) ? LARGER : SMALLER);
+ }
}
-
+
// Compare the slopes of the two supporting curves at p. In the general
// case, the slopes are not equal and their comparison gives us the
- // vertical order to p's right.
- Comparison_result slope_res = _compare_slopes (cv, p);
-
+ // vertical order to p's right.
+ Comparison_result slope_res = _compare_slopes (cv, p, cache);
+
if (slope_res != EQUAL)
return (slope_res);
@@ -1069,7 +1126,7 @@
{
CGAL_precondition (p.compare_xy (left(), cache) != SMALLER);
CGAL_precondition (p.compare_xy (cv.left(), cache) != SMALLER);
-
+
if (this == &cv)
return (EQUAL);
@@ -1113,62 +1170,63 @@
if (_curve.is_same (cv._curve))
{
- // Get the originator, and make sure that p is a vertical tangency
- // point of this originator.
- Originator_iterator org = p.get_originator(_curve);
-
- CGAL_assertion (org != p.originators_end());
+ // Get the originator, and check whether p is a vertical tangency
+ // point of this originator (otherwise it is a self-intersection point,
+ // and we proceed as if it is a regular intersection point).
+ Originator_iterator org = p.get_originator (_curve, _xid);
+
+ CGAL_assertion (org != p.originators_end());
CGAL_assertion (_inc_to_right != cv._inc_to_right);
- if (! p.is_exact())
+ if (org->point_bound().type == Bez_point_bound::VERTICAL_TANGENCY_PT)
{
- CGAL_assertion (org->point_bound().point_type ==
- Bez_point_bound::VERTICAL_TANGENCY_PT);
-
- // Comparison based on the control polygon of the bounded vertical
- // tangency point, using the fact this polygon is y-monotone.
- const typename Bounding_traits::Control_point_vec& cp =
- org->point_bound().bounding_polyline;
+ if (! p.is_exact())
+ {
+ // Comparison based on the control polygon of the bounded vertical
+ // tangency point, using the fact this polygon is y-monotone.
+ const typename Bounding_traits::Control_points& cp =
+ org->point_bound().ctrl;
- if (_inc_to_right)
- {
- return (CGAL::compare(cp.front().y(), cp.back().y()));
+ if (_inc_to_right)
+ {
+ return (CGAL::compare(cp.front().y(), cp.back().y()));
+ }
+ else
+ {
+ return (CGAL::compare(cp.back().y(), cp.front().y()));
+ }
}
- else
- {
- return (CGAL::compare(cp.back().y(), cp.front().y()));
- }
- }
- // Iddo: Handle rational points (using de Casteljau derivative)?
+ // Iddo: Handle rational points (using de Casteljau derivative)?
- // In this case we know that we have a vertical tangency at t0, so
- // X'(t0) = 0. We evaluate the sign of Y'(t0) in order to find the
- // vertical position of the two subcurves to the right of this point.
- CGAL_assertion (org->has_parameter());
+ // In this case we know that we have a vertical tangency at t0, so
+ // X'(t0) = 0. We evaluate the sign of Y'(t0) in order to find the
+ // vertical position of the two subcurves to the right of this point.
+ CGAL_assertion (org->has_parameter());
- const Algebraic& t0 = org->parameter();
- Polynomial polyY_der = nt_traits.derive (_curve.y_polynomial());
- const CGAL::Sign sign_der =
- CGAL::sign (nt_traits.evaluate_at (polyY_der, t0));
+ const Algebraic& t0 = org->parameter();
+ Polynomial polyY_der = nt_traits.derive (_curve.y_polynomial());
+ const CGAL::Sign sign_der =
+ CGAL::sign (nt_traits.evaluate_at (polyY_der, t0));
- CGAL_assertion (sign_der != CGAL::ZERO);
-
- if (_inc_to_right)
- return ((sign_der == CGAL::NEGATIVE) ? LARGER : SMALLER);
- else
- return ((sign_der == CGAL::POSITIVE) ? LARGER : SMALLER);
+ CGAL_assertion (sign_der != CGAL::ZERO);
+
+ if (_inc_to_right)
+ return ((sign_der == CGAL::NEGATIVE) ? LARGER : SMALLER);
+ else
+ return ((sign_der == CGAL::POSITIVE) ? LARGER : SMALLER);
+ }
}
-
+
// Compare the slopes of the two supporting curves at p. In the general
// case, the slopes are not equal and their comparison gives us the
// vertical order to p's right; note that we swap the order of the curves
- // to obtain their position to the left.
- Comparison_result slope_res = cv._compare_slopes (*this, p);
-
+ // to obtains their position to the left.
+ Comparison_result slope_res = cv._compare_slopes (*this, p, cache);
+
if (slope_res != EQUAL)
return (slope_res);
-
+
// Compare the two subcurves by choosing some point to the left of p
// and compareing the vertical position there.
Comparison_result left_res;
@@ -1201,27 +1259,18 @@
// Check if the two subcurve have overlapping supporting curves.
if (! _curve.is_same (cv._curve))
{
- // Ron: For the time being, we do not deal with this case ...
- /*
+ // Check whether the two curves have the same support:
if (! _curve.has_same_support (cv._curve))
return (false);
- // Mark that the two curves overlap.
+ // Mark that the two curves overlap in the cache.
const Curve_id id1 = _curve.id();
const Curve_id id2 = cv._curve.id();
- Curve_pair curve_pair;
-
+
if (id1 < id2)
- curve_pair = Curve_pair (id1, id2);
+ cache.mark_as_overlapping (id1, id2);
else
- curve_pair = Curve_pair (id2, id1);
-
- if (inter_map.find (curve_pair) == inter_map.end())
- {
- Intersection_info& info = inter_map[curve_pair];
- info.second = true;
- }
- */
+ cache.mark_as_overlapping (id2, id1);
}
// Check for equality of the endpoints.
@@ -1237,7 +1286,7 @@
(const Point_2& p,
Self& c1, Self& c2) const
{
- CGAL_precondition (p.get_originator(_curve) != p.originators_end());
+ CGAL_precondition (p.get_originator (_curve, _xid) != p.originators_end());
// Duplicate the curve.
c1 = c2 = *this;
@@ -1245,7 +1294,7 @@
// Perform the split.
if (_dir_right)
{
- c1._pt = p;
+ c1._pt = p;
c2._ps = p;
}
else
@@ -1267,6 +1316,7 @@
// 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) &&
+ _xid == cv._xid &&
(right().is_same (cv.left()) || left().is_same (cv.right())));
return (false);
@@ -1281,20 +1331,17 @@
(const Self& cv) const
{
CGAL_precondition (_curve.is_same (cv._curve));
-
+ CGAL_precondition (_xid == cv._xid);
+
Self res = cv;
-
+
if (right().is_same (cv.left()))
{
// Extend the subcurve to the right.
if (_dir_right)
- {
res._pt = cv.right();
- }
else
- {
res._ps = cv.right();
- }
}
else
{
@@ -1302,13 +1349,9 @@
// Extend the subcurve to the left.
if (_dir_right)
- {
res._ps = cv.left();
- }
else
- {
res._pt = cv.left();
- }
}
return (res);
@@ -1323,25 +1366,27 @@
Bezier_cache& cache) const
{
// First try to use the approximate representation of the endpoints.
- Originator_iterator s_org = _ps.get_originator (_curve);
+ Originator_iterator s_org = _ps.get_originator (_curve, _xid);
CGAL_assertion (s_org != _ps.originators_end());
- Originator_iterator t_org = _pt.get_originator (_curve);
+ Originator_iterator t_org = _pt.get_originator (_curve, _xid);
CGAL_assertion (t_org != _pt.originators_end());
- bool p_lt_ps = (CGAL::compare (t,
- Algebraic (s_org->point_bound().t_min)) ==
- SMALLER);
- bool p_gt_ps = (CGAL::compare (t,
- Algebraic (s_org->point_bound().t_max)) ==
- LARGER);
- bool p_lt_pt = (CGAL::compare (t,
- Algebraic (t_org->point_bound().t_min)) ==
- SMALLER);
- bool p_gt_pt = (CGAL::compare (t,
- Algebraic (t_org->point_bound().t_max)) ==
- LARGER);
+ Nt_traits nt_traits;
+ bool p_lt_ps =
+ (CGAL::compare (t, nt_traits.convert (s_org->point_bound().t_min)) ==
+ SMALLER);
+ bool p_gt_ps =
+ (CGAL::compare (t, nt_traits.convert (s_org->point_bound().t_max)) ==
+ LARGER);
+ bool p_lt_pt =
+ (CGAL::compare (t, nt_traits.convert (t_org->point_bound().t_min)) ==
+ SMALLER);
+ bool p_gt_pt =
+ (CGAL::compare (t, nt_traits.convert (t_org->point_bound().t_max)) ==
+ LARGER);
+
if ((p_gt_ps && p_lt_pt) || (p_lt_ps && p_gt_pt))
{
// The point p is definately in the x-range of the subcurve, as its
@@ -1383,13 +1428,22 @@
return true;
// Compare the parameter of p with the parameters of the endpoints.
- Originator_iterator p_org = p.get_originator (_curve);
- CGAL_assertion (p_org != p.originators_end());
-
- Originator_iterator s_org = _ps.get_originator (_curve);
+ Originator_iterator p_org = p.get_originator (_curve, _xid);
+
+ if (p_org == p.originators_end())
+ {
+ CGAL_assertion (p.get_originator (_curve) != p.originators_end());
+
+ // In this case a different x-monotone curve of the supporting Bezier
+ // curve is an originator of the point, so we know that p does not
+ // lie in the range of our x-monotone subcurve.
+ return (false);
+ }
+
+ Originator_iterator s_org = _ps.get_originator (_curve, _xid);
CGAL_assertion (s_org != _ps.originators_end());
- Originator_iterator t_org = _pt.get_originator (_curve);
+ Originator_iterator t_org = _pt.get_originator (_curve, _xid);
CGAL_assertion (t_org != _pt.originators_end());
bool can_refine_p = ! p.is_exact();
@@ -1439,6 +1493,102 @@
}
// ---------------------------------------------------------------------------
+// Given a point p that lies on the supporting Bezier curve (X(t), Y(t)),
+// determine whether p lies within the t-range of the x-monotone subcurve.
+// If so, the value t0 such that p = (X(t0), Y(t0)) is also computed.
+//
+template <class RatKer, class AlgKer, class NtTrt, class BndTrt>
+bool _Bezier_x_monotone_2<RatKer, AlgKer, NtTrt, BndTrt>::_is_in_range
+ (const Point_2& p,
+ Bezier_cache& cache,
+ Algebraic& t0,
+ bool& is_endpoint) const
+{
+ // The given point p must be rational, otherwise there is no point checking
+ // whether it lies in the interior of the curve.
+ if (! p.is_rational())
+ {
+ is_endpoint = false;
+ return (false);
+ }
+
+ const Rat_point_2& rat_p = (Rat_point_2) p;
+
+ // Determine the parameter range [t_min, t_max] for our x-monotone
+ // subcurve.
+ std::pair<Algebraic, Algebraic> t_range = _t_range (cache);
+ Algebraic t_min, t_max;
+
+ if ((_dir_right && _inc_to_right) || (! _dir_right && ! _inc_to_right))
+ {
+ t_min = t_range.first;
+ t_max = t_range.second;
+ }
+ else
+ {
+ t_min = t_range.second;
+ t_max = t_range.first;
+ }
+
+ // The given point p must lie on (X(t), Y(t)) for some t-value. Obtain the
+ // parameter value t0 for that point. We start by computing all t-values
+ // such that X(t) equals the x-coordinate of p.
+ Nt_traits nt_traits;
+ std::list<Algebraic> t_vals;
+ typename std::list<Algebraic>::iterator t_iter;
+ Comparison_result res1, res2;
+ Algebraic y0;
+
+ _curve.get_t_at_x (rat_p.x(), std::back_inserter(t_vals));
+ CGAL_assertion (! t_vals.empty());
+
+ for (t_iter = t_vals.begin(); t_iter != t_vals.end(); ++t_iter)
+ {
+ // Compare the current t-value with t_min.
+ res1 = CGAL::compare (t_min, *t_iter);
+
+ if (res1 == LARGER)
+ continue;
+
+ // Make sure the y-coordinates match.
+ y0 = nt_traits.evaluate_at (_curve.y_polynomial(), *t_iter) /
+ nt_traits.convert (_curve.y_norm());
+
+ if (CGAL::compare (nt_traits.convert (rat_p.y()), y0) == EQUAL)
+ {
+ if (res1 == EQUAL)
+ {
+ t0 = t_min;
+ is_endpoint = true;
+ return (true);
+ }
+
+ // Compare the current t-value with t_max.
+ res2 = CGAL::compare (t_max, *t_iter);
+
+ if (res2 == EQUAL)
+ {
+ t0 = t_max;
+ is_endpoint = true;
+ return (true);
+ }
+
+ if (res2 == LARGER)
+ {
+ t0 = *t_iter;
+ is_endpoint = false;
+ return (true);
+ }
+ }
+ }
+
+ // In this case, we have not found a t-value in the range of our subcurve,
+ // so p does not lie on the subcurve:
+ is_endpoint = false;
+ return (false);
+}
+
+// ---------------------------------------------------------------------------
// Compute a y-coordinate of a point on the x-monotone subcurve with a
// given x-coordinate.
//
@@ -1451,9 +1601,9 @@
// Obtain the t-values for with the x-coordinates of the supporting
// curve equal x0.
std::list<Algebraic> t_vals;
-
+
_curve.get_t_at_x (x0, std::back_inserter(t_vals));
-
+
// Find a t-value that is in the range of the current curve.
Nt_traits nt_traits;
typename std::list<Algebraic>::iterator t_iter;
@@ -1465,15 +1615,15 @@
for (t_iter = t_vals.begin(); t_iter != t_vals.end(); ++t_iter)
{
res1 = CGAL::compare (*t_iter, t_src);
-
+
if (res1 == EQUAL)
{
// Return the y-coordinate of the source point:
return (_ps.y());
}
-
+
res2 = CGAL::compare (*t_iter, t_trg);
-
+
if (res2 == EQUAL)
{
// Return the y-coordinate of the source point:
@@ -1497,22 +1647,26 @@
// ---------------------------------------------------------------------------
// Compare the slopes of the subcurve with another given Bezier subcurve at
// their given intersection point.
+//
template <class RatKer, class AlgKer, class NtTrt, class BndTrt>
Comparison_result
_Bezier_x_monotone_2<RatKer, AlgKer, NtTrt, BndTrt>::_compare_slopes
(const Self& cv,
- const Point_2& p) const
+ const Point_2& p,
+ Bezier_cache& cache) const
{
// Get the originators of p.
- Originator_iterator org1 = p.get_originator (_curve);
- CGAL_assertion (org1 != p.originators_end());
+ Originator_iterator org1 = p.get_originator (_curve, _xid);
+ const bool valid_org1 = (org1 != p.originators_end());
- Originator_iterator org2 = p.get_originator (cv._curve);
- CGAL_assertion (org2 != p.originators_end());
-
+ Originator_iterator org2 = p.get_originator (cv._curve, cv._xid);
+ const bool valid_org2 = (org2 != p.originators_end());
+
+ CGAL_assertion (valid_org1 || valid_org2);
+
// If the point is only approximated, we can carry out a comparison using
// an approximate number type.
- if (! p.is_exact())
+ if (valid_org1 && valid_org2 && ! p.is_exact())
{
// If the point is inexact, we assume it is a bounded intersection
// point of two curves, and therefore the bounding angle these curves
@@ -1521,21 +1675,52 @@
const Bez_point_bound& bound2 = org2->point_bound();
Bounding_traits bound_tr;
- return (bound_tr.compare_slopes_of_bounded_intersection_point (bound1,
- bound2));
+ return (bound_tr.compare_slopes_at_intersection_point (bound1,
+ bound2));
}
-
- // Iddo: Implement slope comparison at a Rational point
- // (using de Casteljau).
-
- // The slope of (X(t), Y(t)) at t0 is Y'(t)/X'(t).
+
+ // Obtain the parameter values t1 and t2 that correspond to the point p.
+ // Note that it is possible that one of the curves is not an originator
+ // of p. This can happen if p is an endpoint of the other curve (hence
+ // it must be a ratioal point!) and lies in its interior. In this
+ // (degenerate) case we compute the parameter value and set the appropriate
+ // originator for p.
+ Nt_traits nt_traits;
+ Algebraic t1;
+ Algebraic t2;
+
+ if (valid_org1)
+ {
+ CGAL_assertion (org1->has_parameter());
+ t1 = org1->parameter();
+ }
+ else
+ {
+ bool is_endpoint1;
+ CGAL_assertion_code (bool in_range1 =)
+ _is_in_range (p, cache, t1, is_endpoint1);
+ CGAL_assertion (in_range1);
+ p.add_originator (Originator (_curve, _xid, t1));
+ }
+
+ if (valid_org2)
+ {
+ CGAL_assertion (org2->has_parameter());
+ t2 = org2->parameter();
+ }
+ else
+ {
+ bool is_endpoint2;
+ CGAL_assertion_code (bool in_range2 =)
+ cv._is_in_range (p, cache, t2, is_endpoint2);
+ CGAL_assertion (in_range2);
+ p.add_originator (Originator (cv._curve, cv._xid, t2));
+ }
+
+ // The slope of (X(t), Y(t)) at t0 is given by Y'(t0)/X'(t0).
// Compute the slope of (*this).
- // Note that we take special care of the case X'(t) = 0, when the tangent
+ // Note that we take special care of the case X'(t0) = 0, when the tangent
// is vertical and its slope is +/- oo.
- CGAL_assertion (org1->has_parameter() && org2->has_parameter());
-
- Nt_traits nt_traits;
- const Algebraic& t1 = org1->parameter();
Polynomial derivX = nt_traits.derive (_curve.x_polynomial());
Polynomial derivY = nt_traits.derive (_curve.y_polynomial());
Algebraic numer1 = nt_traits.evaluate_at (derivY, t1) *
@@ -1549,7 +1734,7 @@
{
inf_slope1 = CGAL::sign (numer1);
- // If also Y'(t) = 0, we cannot perform the comparison:
+ // If both derivatives are zero, we cannot perform the comparison:
if (inf_slope1 == CGAL::ZERO)
return (EQUAL);
}
@@ -1559,7 +1744,6 @@
}
// Compute the slope of the other subcurve.
- const Algebraic& t2 = org2->parameter();
derivX = nt_traits.derive (cv._curve.x_polynomial());
derivY = nt_traits.derive (cv._curve.y_polynomial());
Algebraic numer2 = nt_traits.evaluate_at (derivY, t2) *
@@ -1573,7 +1757,7 @@
{
inf_slope2 = CGAL::sign (numer2);
- // If also Y'(t) = 0, we cannot perform the comparison:
+ // If both derivatives are zero, we cannot perform the comparison:
if (inf_slope2 == CGAL::ZERO)
return (EQUAL);
}
@@ -1610,10 +1794,10 @@
_Bezier_x_monotone_2<RatKer, AlgKer, NtTrt, BndTrt>::_t_range
(Bezier_cache& cache) const
{
- Originator_iterator ps_org = _ps.get_originator(_curve);
+ Originator_iterator ps_org = _ps.get_originator (_curve, _xid);
CGAL_assertion(ps_org != _ps.originators_end());
- Originator_iterator pt_org = _pt.get_originator(_curve);
+ Originator_iterator pt_org = _pt.get_originator (_curve, _xid);
CGAL_assertion(pt_org != _pt.originators_end());
// Make sure that the two endpoints are exact.
@@ -1642,7 +1826,7 @@
{
// Get the intersection points of the two curves from the cache. Note that
// we make sure that the ID of this->_curve is smaller than of cv's curve ID.
- const bool no_swap_curves = (_curve.id() < cv._curve.id());
+ const bool no_swap_curves = (_curve.id() <= cv._curve.id());
bool do_ovlp;
const Intersect_list& inter_list =
(no_swap_curves ?
@@ -1661,11 +1845,8 @@
_curve.y_polynomial(), _curve.y_norm(),
do_ovlp)));
- if (do_ovlp)
- return (EQUAL);
-
// Get the parameter value for the point p.
- Originator_iterator org = p.get_originator(_curve);
+ Originator_iterator org = p.get_originator (_curve, _xid);
CGAL_assertion (org != p.originators_end());
CGAL_assertion (org->has_parameter());
@@ -1752,36 +1933,158 @@
}
// ---------------------------------------------------------------------------
+// Clip the control polygon of the supporting Bezier curve such that it fits
+// the current x-monotone subcurve.
+//
+template <class RatKer, class AlgKer, class NtTrt, class BndTrt>
+void _Bezier_x_monotone_2<RatKer, AlgKer, NtTrt,
+ BndTrt>::_clip_control_polygon
+ (typename Bounding_traits::Control_points& ctrl,
+ typename Bounding_traits::NT& t_min,
+ typename Bounding_traits::NT& t_max) const
+{
+ // Start from the control polygon of the supporting curve.
+ ctrl.clear();
+ std::copy (_curve.control_points_begin(), _curve.control_points_end(),
+ std::back_inserter (ctrl));
+
+ // The x-monotone subcurve is defined over a parameter range
+ // 0 <= t_min < t_max <= 1. Determine the endpoint with minimal t-value and
+ // the one with maximal t-value.
+ const Point_2& p_min = (_inc_to_right ? left() : right());
+ Originator_iterator org_min = p_min.get_originator (_curve, _xid);
+ const Point_2& p_max = (_inc_to_right ? right() : left());
+ Originator_iterator org_max = p_max.get_originator (_curve, _xid);
+ bool clipped_min = false;
+
+ CGAL_assertion (org_min != p_min.originators_end());
+ CGAL_assertion (org_max != p_max.originators_end());
+
+ // Check if t_min = 0. If so, there is no need to clip.
+ if (! (org_min->point_bound().type == Bez_point_bound::RATIONAL_PT &&
+ CGAL::sign (org_min->point_bound().t_min) == CGAL::ZERO))
+ {
+ // It is possible that the paramater range of the originator is too large.
+ // We therefore make sure it fits the current bounding box of the point
+ // (which we know is tight enough).
+ p_min.fit_to_bbox();
+
+ // Obtain two control polygons, the first for [0, t_min] and the other for
+ // [t_min, 1] and take the second one.
+ typename Bounding_traits::Control_points cp_a;
+ typename Bounding_traits::Control_points cp_b;
+
+ t_min = org_min->point_bound().t_max;
+ de_Casteljau_2 (ctrl.begin(), ctrl.end(),
+ t_min,
+ std::back_inserter(cp_a),
+ std::front_inserter(cp_b));
+
+ ctrl.clear();
+ std::copy (cp_b.begin(), cp_b.end(),
+ std::back_inserter (ctrl));
+ clipped_min = true;
+ }
+ else
+ {
+ t_min = 0;
+ }
+
+ // Check if t_max = 1. If so, there is no need to clip.
+ if (! (org_max->point_bound().type == Bez_point_bound::RATIONAL_PT &&
+ CGAL::compare (org_max->point_bound().t_max, 1) == CGAL::EQUAL))
+ {
+ // It is possible that the paramater range of the originator is too large.
+ // We therefore make sure it fits the current bounding box of the point
+ // (which we know is tight enough).
+ p_max.fit_to_bbox();
+
+ // Obtain two control polygons, the first for [t_min, t_max] and the other
+ // for [t_max, 1] and take the first one.
+ typename Bounding_traits::Control_points cp_a;
+ typename Bounding_traits::Control_points cp_b;
+
+ if (clipped_min)
+ {
+ t_max = (org_max->point_bound().t_min - t_min) / (1 - t_min);
+ }
+ else
+ {
+ t_max = org_max->point_bound().t_min;
+ }
+
+ de_Casteljau_2 (ctrl.begin(), ctrl.end(),
+ t_max,
+ std::back_inserter(cp_a),
+ std::front_inserter(cp_b));
+
+ ctrl.clear();
+ std::copy (cp_a.begin(), cp_a.end(),
+ std::back_inserter (ctrl));
+ t_max = org_max->point_bound().t_min;
+ }
+ else
+ {
+ t_max = 1;
+ }
+
+ return;
+}
+
+// ---------------------------------------------------------------------------
// Approximate the intersection points between the two given curves.
//
template <class RatKer, class AlgKer, class NtTrt, class BndTrt>
bool _Bezier_x_monotone_2<RatKer, AlgKer, NtTrt,
BndTrt>::_approximate_intersection_points
- (const Curve_2& B1,
- const Curve_2& B2,
+ (const Self& cv,
std::list<Point_2>& inter_pts) const
{
+ typedef typename Bounding_traits::Intersection_point Intersection_point;
+
inter_pts.clear();
- // Make local copies of the control polygons.
- typename Bounding_traits::Control_point_vec cp1, cp2;
-
- std::copy (B1.control_points_begin(), B1.control_points_end(),
- std::back_inserter (cp1));
- std::copy (B2.control_points_begin(), B2.control_points_end(),
- std::back_inserter(cp2));
+ // Get the supporting Bezier curves, and make local copies of their control
+ // polygons.
+ const Curve_2& B1 = this->_curve;
+ const Curve_2& B2 = cv._curve;
+ typename Bounding_traits::Control_points cp1;
+ typename Bounding_traits::NT t_min1 = 0, t_max1 = 1;
+ typename Bounding_traits::Control_points cp2;
+ typename Bounding_traits::NT t_min2 = 0, t_max2 = 1;
+ bool is_self_intersection = false;
+ if (! B1.is_same (B2))
+ {
+ // In case B1 and B2 are different curves, use their full control polygons
+ // in order to approximate all intersection points between the two
+ // supporting Bezier curves.
+ std::copy (B1.control_points_begin(), B1.control_points_end(),
+ std::back_inserter (cp1));
+ std::copy (B2.control_points_begin(), B2.control_points_end(),
+ std::back_inserter(cp2));
+ }
+ else
+ {
+ // In this case we need to approximate the (self-)intersection points of
+ // two subcurves of the same Bezier curve. Clip the control polygons and
+ // obtain the control polygons of the two subcurves.
+ _clip_control_polygon (cp1, t_min1, t_max1);
+ cv._clip_control_polygon (cp2, t_min2, t_max2);
+ is_self_intersection = true;
+ }
+
// Use the bounding traits to isolate the intersection points.
- Bounding_traits bound_tr;
- typename Bounding_traits::Bound_pair_lst inter_pairs;
+ Bounding_traits bound_tr;
+ std::list<Intersection_point> ipt_bounds;
- bound_tr.CrvCrvInter (cp1, cp2,
- inter_pairs);
-
+ bound_tr.compute_intersection_points (cp1, cp2,
+ std::back_inserter (ipt_bounds));
+
// Construct the approximated points.
- typename Bounding_traits::Bound_pair_lst::const_iterator iter;
+ typename std::list<Intersection_point>::const_iterator iter;
- for (iter = inter_pairs.begin(); iter != inter_pairs.end(); ++iter)
+ for (iter = ipt_bounds.begin(); iter != ipt_bounds.end(); ++iter)
{
const Bez_point_bound& bound1 = iter->bound1;
const Bez_point_bound& bound2 = iter->bound2;
@@ -1794,24 +2097,63 @@
// Create the approximated intersection point.
Point_2 pt;
- if (bound1.point_type == Bounding_traits::Bez_point_bound::RATIONAL_PT &&
- bound2.point_type == Bounding_traits::Bez_point_bound::RATIONAL_PT)
+ if (bound1.type == Bounding_traits::Bez_point_bound::RATIONAL_PT &&
+ bound2.type == Bounding_traits::Bez_point_bound::RATIONAL_PT)
{
CGAL_assertion (CGAL::compare (bound1.t_min, bound1.t_max) == EQUAL);
CGAL_assertion (CGAL::compare (bound2.t_min, bound2.t_max) == EQUAL);
- Rational t1 = bound1.t_min;
- Rational t2 = bound2.t_min;
+ Rational t1 = bound1.t_min;
+ Rational t2 = bound2.t_min;
+ Nt_traits nt_traits;
- pt = Point_2 (B1, t1);
- pt.add_originator (Originator (B2, t2));
+ if (is_self_intersection)
+ {
+ // Set the originators with the curve x-monotone IDs.
+ // Note that the parameter values we have computed relate to the
+ // parameter range [t_min1, t_max1] and [t_min2, t_max2], respectively,
+ // so we scale them back to the parameter range [0, 1] that represent
+ // the entire curve.
+ t1 = t_min1 + t1 * (t_max1 - t_min1);
+ t2 = t_min2 + t2 * (t_max2 - t_min2);
+
+ pt = Point_2 (B1, _xid, t1);
+ pt.add_originator (Originator (B2, cv._xid, t2));
+ }
+ else
+ {
+ // Set the originators referring to the entire supporting curves.
+ pt = Point_2 (B1, t1);
+ pt.add_originator (Originator (B2, nt_traits.convert (t2)));
+ }
}
else
{
- pt.add_originator (Originator (B1, bound1));
- pt.add_originator (Originator (B2, bound2));
+ if (is_self_intersection)
+ {
+ // Set the originators with the curve x-monotone IDs.
+ // Note that the parameter values we have computed relate to the
+ // parameter range [t_min1, t_max1] and [t_min2, t_max2], respectively,
+ // so we scale them back to the parameter range [0, 1] that represent
+ // the entire curve.
+ Bez_point_bound sc_bound1 = bound1;
+ sc_bound1.t_min = t_min1 + bound1.t_min * (t_max1 - t_min1);
+ sc_bound1.t_max = t_min1 + bound1.t_max * (t_max1 - t_min1);
+ pt.add_originator (Originator (B1, _xid, sc_bound1));
+
+ Bez_point_bound sc_bound2 = bound2;
+ sc_bound2.t_min = t_min2 + bound2.t_min * (t_max2 - t_min2);
+ sc_bound2.t_max = t_min2 + bound2.t_max * (t_max2 - t_min2);
+ pt.add_originator (Originator (B2, cv._xid, sc_bound2));
+ }
+ else
+ {
+ // Set the originators with the curve x-monotone IDs.
+ pt.add_originator (Originator (B1, bound1));
+ pt.add_originator (Originator (B2, bound2));
+ }
}
pt.set_bbox (bbox);
-
+
inter_pts.push_back (pt);
}
@@ -1828,12 +2170,27 @@
Intersection_map& inter_map,
Bezier_cache& cache,
std::vector<Intersection_point_2>& ipts,
- Self& /* ovlp_cv */) const
+ Self& ovlp_cv) const
{
- CGAL_precondition (_curve.id() < cv._curve.id());
+ CGAL_precondition (_curve.id() <= cv._curve.id());
ipts.clear();
+ // In case the two x-monotone curves are subcurves of the same Bezier curve,
+ // first check if this base curve is not self-intersecting. If this is the
+ // case we can avoid any attempt of computing intersection points between
+ // the two subcurves.
+ const bool self_intersect = (_curve.id() == cv._curve.id());
+
+ if (self_intersect)
+ {
+ if (_xid == cv._xid)
+ return (false);
+
+ if (_curve.has_no_self_intersections())
+ return (false);
+ }
+
// Construct the pair of curve IDs and look for it in the intersection map.
Curve_pair curve_pair (_curve.id(), cv._curve.id());
Intersection_map_iterator map_iter = inter_map.find (curve_pair);
@@ -1849,11 +2206,13 @@
else
{
// Approximate the intersection points and store them in the map.
- app_ok = _approximate_intersection_points (_curve,
- cv._curve,
+ // Note that we do not store approximated self-intersections in the map,
+ // as they realte only to the pecific x-monotone curves, and not to the
+ // entire curve.
+ app_ok = _approximate_intersection_points (cv,
inter_pts);
- if (app_ok)
+ if (app_ok && ! self_intersect)
inter_map[curve_pair] = inter_pts;
}
@@ -1866,38 +2225,56 @@
// If the approximation went OK, then we know that we just have simple
// intersection points (with multiplicity 1). We go over the points
// and report the ones lying in the parameter ranges of both curves.
+ // Note that in case of self-intersections, all points we get are in
+ // the respective parameter range of the curves.
typename std::list<Point_2>::iterator pit;
for (pit = inter_pts.begin(); pit != inter_pts.end(); ++pit)
{
// Check if the point is in the range of this curve - first using
// its parameter bounds, and if we fail we perform an exact check.
- in_range1 = _is_in_range (*pit, correct_res);
-
+ if (! self_intersect)
+ {
+ in_range1 = _is_in_range (*pit, correct_res);
+ }
+ else
+ {
+ in_range1 = true;
+ correct_res = true;
+ }
+
if (! correct_res)
{
if (! pit->is_exact())
pit->make_exact (cache);
- Originator_iterator p_org = pit->get_originator (_curve);
+ Originator_iterator p_org = pit->get_originator (_curve, _xid);
CGAL_assertion (p_org != pit->originators_end());
in_range1 = _is_in_range (p_org->parameter(), cache);
}
-
+
if (! in_range1)
continue;
// Check if the point is in the range of the other curve - first using
// its parameter bounds, and if we fail we perform an exact check.
- in_range2 = cv._is_in_range (*pit, correct_res);
-
+ if (! self_intersect)
+ {
+ in_range2 = cv._is_in_range (*pit, correct_res);
+ }
+ else
+ {
+ in_range2 = true;
+ correct_res = true;
+ }
+
if (! correct_res)
{
if (! pit->is_exact())
pit->make_exact (cache);
- Originator_iterator p_org = pit->get_originator (cv._curve);
+ Originator_iterator p_org = pit->get_originator (cv._curve, cv._xid);
CGAL_assertion (p_org != pit->originators_end());
in_range2 = cv._is_in_range (p_org->parameter(), cache);
@@ -1905,13 +2282,28 @@
if (in_range1 && in_range2)
{
+ // In case the originators of the intersection point are not marked
+ // with x-monotone identifiers, mark them now as we know in which
+ // subcurves they lie.
+ Originator_iterator p_org1 = pit->get_originator (_curve, _xid);
+ CGAL_assertion (p_org1 != pit->originators_end());
+
+ if (p_org1->xid() == 0)
+ pit->update_originator_xid (*p_org1, _xid);
+
+ Originator_iterator p_org2 = pit->get_originator (cv._curve, cv._xid);
+ CGAL_assertion (p_org2 != pit->originators_end());
+
+ if (p_org2->xid() == 0)
+ pit->update_originator_xid (*p_org2, cv._xid);
+
// The point lies within the parameter range of both curves, so we
// report it as a valid intersection point with multiplicity 1.
ipts.push_back (Intersection_point_2 (*pit, 1));
}
}
- // Since the apporximation went fine we cannot possible have an overlap:
+ // Since the apporximation went fine we cannot possibly have an overlap:
return (false);
}
@@ -1929,9 +2321,174 @@
if (do_ovlp)
{
- // Ron: TODO -- Handle overlapping curves ...
- CGAL_assertion (false);
- return (true);
+ // Check the case of co-inciding endpoints
+ if (left().equals (cv.left(), cache))
+ {
+ if (right().equals (cv.right(), cache))
+ {
+ // The two curves entirely overlap one another:
+ ovlp_cv = cv;
+ return (true);
+ }
+
+ Algebraic t_right;
+ bool is_endpoint;
+
+ if (_is_in_range (cv.right(), cache, t_right, is_endpoint))
+ {
+ CGAL_assertion (! is_endpoint);
+
+ // Case 1 - *this: s +-----------+ t
+ // cv: s'+=====+ t'
+ //
+ // Take cv as the overlapping subcurve, and add originators for its
+ // right endpoint referring to *this.
+ ovlp_cv = cv;
+ ovlp_cv.right().add_originator (Originator (_curve, _xid, t_right));
+
+ return (true);
+ }
+ else if (cv._is_in_range (right(), cache, t_right, is_endpoint))
+ {
+ CGAL_assertion (! is_endpoint);
+
+ // Case 2 - *this: s +----+ t
+ // cv: s'+==========+ t'
+ //
+ // Take this as the overlapping subcurve, and add originators for its
+ // right endpoint referring to cv.
+ ovlp_cv = *this;
+ ovlp_cv.right().add_originator (Originator (cv._curve, cv._xid,
+ t_right));
+
+ return (true);
+ }
+
+ // In this case the two curves do not overlap, but have a common left
+ // endpoint.
+ ipts.push_back (Intersection_point_2 (left(), 0));
+ return (false);
+ }
+ else if (right().equals (cv.right(), cache))
+ {
+ Algebraic t_left;
+ bool is_endpoint;
+
+ if (_is_in_range (cv.left(), cache, t_left, is_endpoint))
+ {
+ CGAL_assertion (! is_endpoint);
+
+ // Case 3 - *this: s +-----------+ t
+ // cv: s'+=====+ t'
+ //
+ // Take cv as the overlapping subcurve, and add originators for its
+ // left endpoint referring to *this.
+ ovlp_cv = cv;
+ ovlp_cv.left().add_originator (Originator (_curve, _xid, t_left));
+
+ return (true);
+ }
+ else if (cv._is_in_range (left(), cache, t_left, is_endpoint))
+ {
+ CGAL_assertion (! is_endpoint);
+
+ // Case 4 - *this: s +----+ t
+ // cv: s'+==========+ t'
+ //
+ // Take this as the overlapping subcurve, and add originators for its
+ // left endpoint referring to cv.
+ ovlp_cv = *this;
+ ovlp_cv.left().add_originator (Originator (cv._curve, cv._xid,
+ t_left));
+
+ return (true);
+ }
+
+ // In this case the two curves do not overlap, but have a common right
+ // endpoint.
+ ipts.push_back (Intersection_point_2 (right(), 0));
+ return (false);
+ }
+
+ // If we reached here, none of the endpoints coincide.
+ // Check the possible overlap scenarios.
+ Point_2 ovrp_src, ovlp_trg;
+ Algebraic t_cv_src;
+ Algebraic t_cv_trg;
+ bool is_endpoint = false;
+
+ if (_is_in_range (cv._ps, cache, t_cv_src, is_endpoint) &&
+ ! is_endpoint)
+ {
+ if (_is_in_range (cv._pt, cache, t_cv_trg, is_endpoint) &&
+ ! is_endpoint)
+ {
+ // Case 5 - *this: s +-----------+ t
+ // cv: s' +=====+ t'
+ //
+ // Take cv as the overlapping subcurve, and add originators for its
+ // endpoints referring to *this.
+ ovlp_cv = cv;
+ ovlp_cv._ps.add_originator (Originator (_curve, _xid, t_cv_src));
+ ovlp_cv._pt.add_originator (Originator (_curve, _xid, t_cv_trg));
+
+ return (true);
+ }
+ else
+ {
+ // Case 6 - *this: s +-----------+ t
+ // cv: s' +=====+ t'
+ //
+ // Use *this as a base, and replace its source point.
+ ovlp_cv = *this;
+ ovlp_cv._ps = cv._ps;
+ ovlp_cv._ps.add_originator (Originator (_curve, _xid, t_cv_src));
+
+ // Add an originator to the target point, referring to cv:
+ CGAL_assertion_code (bool pt_in_cv_range =)
+ cv._is_in_range (ovlp_cv._pt, cache, t_cv_trg, is_endpoint);
+ CGAL_assertion (pt_in_cv_range);
+ ovlp_cv._pt.add_originator (Originator (cv._curve, cv._xid, t_cv_trg));
+
+ return (true);
+ }
+ }
+ else if (_is_in_range (cv._pt, cache, t_cv_trg, is_endpoint) &&
+ ! is_endpoint)
+ {
+ // Case 7 - *this: s +-----------+ t
+ // cv: s' +=====+ t'
+ //
+ // Use *this as a base, and replace its target point.
+ ovlp_cv = *this;
+ ovlp_cv._pt = cv._pt;
+ ovlp_cv._pt.add_originator (Originator (_curve, _xid, t_cv_trg));
+
+ // Add an originator to the source point, referring to cv:
+ CGAL_assertion_code (bool ps_in_cv_range =)
+ cv._is_in_range (ovlp_cv._ps, cache, t_cv_src, is_endpoint);
+ CGAL_assertion (ps_in_cv_range);
+ ovlp_cv._ps.add_originator (Originator (cv._curve, cv._xid, t_cv_src));
+
+ return (true);
+ }
+ else if (cv._is_in_range (_ps, cache, t_cv_src, is_endpoint) &&
+ cv._is_in_range (_pt, cache, t_cv_trg, is_endpoint))
+ {
+ // Case 8 - *this: s +---------+ t
+ // cv: s' +================+ t'
+ //
+ // Take *this as the overlapping subcurve, and add originators for its
+ // endpoints referring to cv.
+ ovlp_cv = *this;
+ ovlp_cv._ps.add_originator (Originator (cv._curve, cv._xid, t_cv_src));
+ ovlp_cv._pt.add_originator (Originator (cv._curve, cv._xid, t_cv_trg));
+
+ return (true);
+ }
+
+ // If we reached here, there are no overlaps:
+ return (false);
}
// Go over the points and report the ones lying in the parameter ranges
@@ -1945,10 +2502,10 @@
{
// Construct an intersection point with unknown multiplicity.
Point_2 pt (iit->x, iit->y);
-
- pt.add_originator (Originator (_curve, iit->s));
- pt.add_originator (Originator (cv._curve, iit->t));
-
+
+ pt.add_originator (Originator (_curve, _xid, iit->s));
+ pt.add_originator (Originator (cv._curve, cv._xid, iit->t));
+
ipts.push_back (Intersection_point_2 (pt, 0));
}
}
@@ -1975,10 +2532,10 @@
// Get a rational approximation of the parameter values at the endpoints.
Nt_traits nt_traits;
- Originator_iterator ps_org = _ps.get_originator(_curve);
+ Originator_iterator ps_org = _ps.get_originator (_curve, _xid);
CGAL_assertion (ps_org != _ps.originators_end());
- Originator_iterator pt_org = _pt.get_originator(_curve);
+ Originator_iterator pt_org = _pt.get_originator (_curve, _xid);
CGAL_assertion (pt_org != _pt.originators_end());
Rational my_t_min;
Index: Arr_traits_2/Bezier_point_2.h
===================================================================
--- Arr_traits_2/Bezier_point_2.h (revision 40349)
+++ Arr_traits_2/Bezier_point_2.h (working copy)
@@ -88,33 +88,57 @@
{
private:
- Curve_2 _curve; /*! The originating curve. */
- Bez_point_bound _bpb; /*! Bounding information for the
- point: bouding control polygon,
- point type, etc. */
- Algebraic *p_t; /*! The algebraic parameter for the
- point (if available). */
+ Curve_2 _curve; /*!< The originating curve. */
+ unsigned int _xid; /*!< Serial number of the originating
+ x-monotone curve. */
+ Bez_point_bound _bpb; /*!< Bounding information for the
+ point: bouding control polygon,
+ point type, etc. */
+ Algebraic *p_t; /*!< The algebraic parameter for the
+ point (if available). */
public:
/*! Constructor, given an exact algebraic representation. */
- Originator(const Curve_2& c, const Algebraic& t) :
- _curve(c),
+ Originator (const Curve_2& c, const Algebraic& t) :
+ _curve (c),
+ _xid (0),
p_t (NULL)
{
set_parameter (t);
}
+ /*! Constructor, given an exact algebraic representation. */
+ Originator (const Curve_2& c, unsigned int xid,
+ const Algebraic& t) :
+ _curve (c),
+ _xid (xid),
+ p_t (NULL)
+ {
+ set_parameter (t);
+ }
+
/*! Constructor with bounding information and no exact representation. */
Originator (const Curve_2& c, const Bez_point_bound& bpb) :
_curve (c),
+ _xid (0),
_bpb (bpb),
p_t (NULL)
{}
+ /*! Constructor with bounding information and no exact representation. */
+ Originator (const Curve_2& c, unsigned int xid,
+ const Bez_point_bound& bpb) :
+ _curve (c),
+ _xid (xid),
+ _bpb (bpb),
+ p_t (NULL)
+ {}
+
/*! Copy constructor. */
Originator (const Originator& other) :
_curve (other._curve),
+ _xid (other._xid),
_bpb (other._bpb),
p_t (NULL)
{
@@ -136,7 +160,7 @@
// Avoid self assignments.
if (this == &other)
return (*this);
-
+
// Free memory, if necessary.
if (p_t != NULL)
delete p_t;
@@ -144,6 +168,7 @@
// Copy the data members.
_curve = other._curve;
+ _xid = other._xid;
_bpb = other._bpb;
// Deep copy of lazy instantiation
@@ -159,6 +184,12 @@
return (_curve);
}
+ /*! Get the serial number of the originating x-monotone curve. */
+ unsigned int xid () const
+ {
+ return (_xid);
+ }
+
/*! Get the bounding information. */
const Bez_point_bound& point_bound () const
{
@@ -198,8 +229,8 @@
p_t = new Algebraic (t);
- // Update the Bez_point_bound, probably by converting t to
- // an interval of doubles and setting _bpb accordingly.
+ // Update the Bez_point_bound by converting t to an interval of doubles
+ // and setting _bpb accordingly.
Nt_traits nt_traits;
const std::pair<double, double>& t_bnd = nt_traits.double_interval (t);
@@ -208,28 +239,50 @@
return;
}
+
+ /*!
+ * Set the serial number of the originating x-monotone curve.
+ * \param xid the new serial number of the originating x-monotone curve.
+ * \pre The current xid() is 0.
+ * \pre xid is possitive.
+ */
+ void set_xid (unsigned int xid)
+ {
+ CGAL_precondition (_xid == 0);
+ CGAL_precondition (xid > 0);
+
+ _xid = xid;
+ return;
+ }
};
/*! \struct Subcurve
* Auxilary structure for the vertical_position() function.
*/
- typedef typename Bounding_traits::Control_point_vec Control_point_vec;
- typedef typename Bounding_traits::NT BoundNT;
+ typedef typename Bounding_traits::Control_points Control_points;
+ typedef typename Bounding_traits::NT BoundNT;
struct Subcurve
{
- Control_point_vec cp;
- BoundNT l;
- BoundNT r;
+ Control_points ctrl; /*!< The control points. */
+ BoundNT t_min; /*!< Minimal parameter value. */
+ BoundNT t_max; /*!< Maximal parameter value. */
- /*! Constructor. */
- Subcurve (const Control_point_vec& _cp,
- const BoundNT& _l,
- const BoundNT& _r) :
- cp (_cp),
- l (_l),
- r (_r)
- {}
+ /*! Constructor given control points an a t-range. */
+ Subcurve (const Control_points& _ctrl,
+ const BoundNT& _tmin,
+ const BoundNT& _tmax) :
+ ctrl (_ctrl),
+ t_min (_tmin),
+ t_max (_tmax)
+ {}
+
+ /*! Constructor given a t-range. */
+ Subcurve (const BoundNT& _tmin,
+ const BoundNT& _tmax) :
+ t_min (_tmin),
+ t_max (_tmax)
+ {}
};
typedef std::list<Originator> Orig_list;
@@ -304,11 +357,25 @@
_Bezier_point_2_rep (const Curve_2& B, const Rational& t0);
/*!
+ * Constructor given an x-monotone curve and a rational t0 value.
+ * \pre t0 must be between 0 and 1.
+ */
+ _Bezier_point_2_rep (const Curve_2& B, unsigned int xid,
+ const Rational& 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);
+ /*!
+ * Constructor given an x-monotone curve and an algebraic t0 value.
+ * \pre t0 must be between 0 and 1.
+ */
+ _Bezier_point_2_rep (const Curve_2& B, unsigned int xid,
+ const Algebraic& t0);
+
/*! Destructor. */
~_Bezier_point_2_rep ()
{
@@ -396,7 +463,7 @@
* LARGER if the point is located above the curve;
* EQUAL if we cannot determine its precise position.
*/
- Comparison_result vertical_position (const Control_point_vec& cp,
+ Comparison_result vertical_position (const Control_points& cp,
const BoundNT& t_min,
const BoundNT& t_max);
@@ -409,6 +476,11 @@
bool _refine ();
/*!
+ * Make sure the originator parameters fit the bound box.
+ */
+ void _fit_to_bbox ();
+
+ /*!
* Compute the exact representation of the point.
* \param cache A cache for the vertical tangency points and the
* intersection points.
@@ -487,6 +559,15 @@
{}
/*!
+ * Constructor given an x-monotone curve and a rational t0 value.
+ * \pre t0 must be between 0 and 1.
+ */
+ _Bezier_point_2 (const Curve_2& B, unsigned int xid,
+ const Rational& t0) :
+ Bpt_handle (Bpt_rep (B, xid, t0))
+ {}
+
+ /*!
* Constructor given an originating curve and an algebraic t0 value.
* \pre t0 must be between 0 and 1.
*/
@@ -495,6 +576,15 @@
{}
/*!
+ * Constructor given an x-monotone curve and an algebraic t0 value.
+ * \pre t0 must be between 0 and 1.
+ */
+ _Bezier_point_2 (const Curve_2& B, unsigned int xid,
+ const Algebraic& t0) :
+ Bpt_handle (Bpt_rep (B, xid, t0))
+ {}
+
+ /*!
* Assignment operator.
*/
Self& operator= (const Self& pt)
@@ -594,6 +684,16 @@
}
/*!
+ * Make sure the originator parameters fit the bound box.
+ */
+ void fit_to_bbox () const
+ {
+ Bpt_rep& rep = const_cast<Bpt_rep&> (_rep());
+
+ return (rep._fit_to_bbox());
+ }
+
+ /*!
* Compute the exact coordinates of the point.
*/
void make_exact (Bezier_cache& cache) const
@@ -674,7 +774,7 @@
* EQUAL if we cannot determine its precise position.
*/
Comparison_result vertical_position
- (const typename Bounding_traits::Control_point_vec& cp,
+ (const typename Bounding_traits::Control_points& cp,
const typename Bounding_traits::NT& t_min,
const typename Bounding_traits::NT& t_max) const
{
@@ -683,16 +783,12 @@
}
/*!
- * Get the originator of the point that is associates with the given curve.
+ * Get the originator of the point that is associated with the given curve.
* \param B The Bezier curve.
* \return An iterator pointing to the requested originator;
- * or originators_end() if B is not an originator of the point.
+ * originators_end() if B is not an originator of the point.
*/
- // Iddo: this is a bit const-problematic since it should return
- // const_iterator, currently Originator_iterator is typedefed to a
- // const_iterator.
- // (TODO - Originator_const_iterator and Originator_iterator)
- Originator_iterator get_originator(const Curve_2& B) const
+ Originator_iterator get_originator (const Curve_2& B) const
{
// Scan the list of originators and look for B.
typename Bpt_rep::Orig_const_iter it = _rep()._origs.begin();
@@ -701,14 +797,49 @@
while (it != end)
{
if (B.is_same (it->curve()))
- break;
+ return (it);
+
++it;
}
- return it;
+ // If we reached here, we have not found an originator:
+ return (it);
}
/*!
+ * Get the originator of the point that is associated with the given
+ * x-monotone curve.
+ * \param B The Bezier curve.
+ * \param xid The serial number of the x-monotone subcurve.
+ * \return An iterator pointing to the requested originator;
+ * originators_end() if B is not an originator of the point.
+ */
+ Originator_iterator get_originator (const Curve_2& B,
+ unsigned int xid) const
+ {
+ // Scan the list of originators and look for B.
+ typename Bpt_rep::Orig_const_iter it = _rep()._origs.begin();
+ typename Bpt_rep::Orig_const_iter end = _rep()._origs.end();
+
+ while (it != end)
+ {
+ if (B.is_same (it->curve()))
+ {
+ // An x-monotone id that equals 0 means that the originator is not
+ // associated with a specific x-monotone curve. Otherwise, we require
+ // that the IDs match.
+ if (it->xid() == 0 || it->xid() == xid)
+ return (it);
+ }
+
+ ++it;
+ }
+
+ // If we reached here, we have not found an originator:
+ return (it);
+ }
+
+ /*!
* Get the range of originators.
*/
Originator_iterator originators_begin () const
@@ -733,6 +864,18 @@
}
/*!
+ * Update the xid field of the given orinigator.
+ */
+ void update_originator_xid (const Originator& o,
+ unsigned int xid) const
+ {
+ Originator& orig = const_cast<Originator&> (o);
+
+ orig.set_xid (xid);
+ return;
+ }
+
+ /*!
* Add the originators of the given point.
*/
void merge_originators (const Self& pt) const
@@ -749,7 +892,6 @@
return;
}
- // Iddo: workaround the ctr problems
/*! Set the bounding box for the point. */
void set_bbox (const Bez_point_bbox& bbox)
{
@@ -820,12 +962,57 @@
// Insert an originator of a rational point.
// Note that this constructor also takes care of the Bez_bound
// for the originator.
- _origs.push_back (Originator(B, t0));
+ Nt_traits nt_traits;
+ Originator org (B, nt_traits.convert (t0));
+ Bez_point_bound bound = org.point_bound();
+ bound.type = Bez_point_bound::RATIONAL_PT;
+ org.update_point_bound (bound);
+ _origs.push_back (org);
+
// Evaluate the point coordinates.
const Rat_point_2& p = B (t0);
+
+ p_rat_x = new Rational (p.x());
+ p_alg_x = new Algebraic (nt_traits.convert (*p_rat_x));
+ p_rat_y = new Rational (p.y());
+ p_alg_y = new Algebraic (nt_traits.convert (*p_rat_y));
+
+ // Also set the bounding box for this point, by converting x, y
+ // to two ranges of doubles.
+ const std::pair<double, double>& x_bnd =
+ nt_traits.double_interval (*p_alg_x);
+ const std::pair<double, double>& y_bnd =
+ nt_traits.double_interval (*p_alg_y);
+
+ _bbox.min_x = x_bnd.first;
+ _bbox.max_x = x_bnd.second;
+ _bbox.min_y = y_bnd.first;
+ _bbox.max_y = y_bnd.second;
+}
+
+// ---------------------------------------------------------------------------
+// Constructor given an x-monotone curve and a rational t0 value.
+//
+template <class RatKer, class AlgKer, class NtTrt, class BndTrt>
+_Bezier_point_2_rep<RatKer, AlgKer, NtTrt, BndTrt>::_Bezier_point_2_rep
+ (const Curve_2& B, unsigned int xid,
+ const Rational& t0)
+{
+ // Insert an originator of a rational point.
+ // Note that this constructor also takes care of the Bez_bound
+ // for the originator.
Nt_traits nt_traits;
+ Originator org (B, nt_traits.convert (t0));
+ Bez_point_bound bound = org.point_bound();
+ bound.type = Bez_point_bound::RATIONAL_PT;
+ org.update_point_bound (bound);
+ _origs.push_back (org);
+
+ // Evaluate the point coordinates.
+ const Rat_point_2& p = B (t0);
+
p_rat_x = new Rational (p.x());
p_alg_x = new Algebraic (nt_traits.convert (*p_rat_x));
p_rat_y = new Rational (p.y());
@@ -879,6 +1066,41 @@
}
// ---------------------------------------------------------------------------
+// Constructor given an x-monotone curve and an algebraic t0 value.
+//
+template <class RatKer, class AlgKer, class NtTrt, class BndTrt>
+_Bezier_point_2_rep<RatKer, AlgKer, NtTrt, BndTrt>::_Bezier_point_2_rep
+ (const Curve_2& B, unsigned int xid,
+ const Algebraic& t0) :
+ p_rat_x (NULL),
+ p_rat_y (NULL)
+{
+ // Create the originator pair <B(t), t0>.
+ // Note that this constructor also takes care of the Bez_bound
+ // for the originator.
+ _origs.push_back (Originator (B, xid, t0));
+
+ // Set the point coordinates.
+ const Alg_point_2 p = B (t0);
+
+ p_alg_x = new Algebraic (p.x());
+ p_alg_y = new Algebraic (p.y());
+
+ // Also set the bounding box for this point, by converting x, y
+ // to two ranges of doubles.
+ Nt_traits nt_traits;
+ const std::pair<double, double>& x_bnd =
+ nt_traits.double_interval (*p_alg_x);
+ const std::pair<double, double>& y_bnd =
+ nt_traits.double_interval (*p_alg_y);
+
+ _bbox.min_x = x_bnd.first;
+ _bbox.max_x = x_bnd.second;
+ _bbox.min_y = y_bnd.first;
+ _bbox.max_y = y_bnd.second;
+}
+
+// ---------------------------------------------------------------------------
// Compare the x-coordinate to this of the given point.
//
template <class RatKer, class AlgKer, class NtTrt, class BndTrt>
@@ -890,7 +1112,7 @@
// Handle rational points first.
if (is_rational() && pt.is_rational())
{
- return (CGAL::compare(*p_rat_x, *(pt.p_rat_x)));
+ return (CGAL::compare (*p_rat_x, *(pt.p_rat_x)));
}
// Try to handle the comparison using the x-range of the bounding boxes.
@@ -907,21 +1129,27 @@
return (LARGER);
// Check if only one of the points is exactly represented.
+ Nt_traits nt_traits;
+
if (is_exact())
{
// Compare the exact x-coordinate to pt's bounding box.
- if (CGAL::compare (*p_alg_x, Algebraic (pt._bbox.min_x)) == SMALLER)
+ if (CGAL::compare (*p_alg_x,
+ nt_traits.convert (pt._bbox.min_x)) == SMALLER)
return (SMALLER);
- if (CGAL::compare (*p_alg_x, Algebraic (pt._bbox.max_x)) == LARGER)
+ if (CGAL::compare (*p_alg_x,
+ nt_traits.convert (pt._bbox.max_x)) == LARGER)
return (LARGER);
}
if (pt.is_exact())
{
// Compare the bounding box to pt's exact x-coordinate.
- if (CGAL::compare (Algebraic (_bbox.max_x), *(pt.p_alg_x)) == SMALLER)
+ if (CGAL::compare (nt_traits.convert (_bbox.max_x),
+ *(pt.p_alg_x)) == SMALLER)
return (SMALLER);
- if (CGAL::compare (Algebraic (_bbox.min_x), *(pt.p_alg_x)) == LARGER)
+ if (CGAL::compare (nt_traits.convert (_bbox.min_x),
+ *(pt.p_alg_x)) == LARGER)
return (LARGER);
}
@@ -996,26 +1224,23 @@
template <class RatKer, class AlgKer, class NtTrt, class BndTrt>
Comparison_result
_Bezier_point_2_rep<RatKer, AlgKer, NtTrt, BndTrt>::vertical_position
- (const Control_point_vec& cp,
+ (const Control_points& cp,
const BoundNT& t_min,
const BoundNT& t_max)
{
+ // Initialize a list of subcurves of the original curve.
+ // We start from the entire curve.
Bounding_traits bound_tr;
std::list<Subcurve> subcurves;
- subcurves.push_back(Subcurve(cp,0,1));
+ subcurves.push_back (Subcurve (cp, 0, 1));
- // Iddo: maybe make here a comparison with an incrementor (e.g., loop 4
- // times first..)
bool can_refine_pt = true;
bool can_refine_cv = true;
Bez_point_bbox scv_bbox;
while (can_refine_pt || can_refine_cv)
{
- // Iddo: Implement Algebraic comparisons with scv_bbox if is_exact
- // (or use the current bbox of the exact rep as is done now).
-
// Go over the list of subcurves and consider only those lying in the
// given [t_min, t_max] bound.
typename std::list<Subcurve>::iterator iter = subcurves.begin();
@@ -1024,8 +1249,8 @@
while (iter != subcurves.end())
{
- if (CGAL::compare (iter->r, t_min) == SMALLER ||
- CGAL::compare (iter->l, t_max) == LARGER)
+ if (CGAL::compare (iter->t_max, t_min) == SMALLER ||
+ CGAL::compare (iter->t_min, t_max) == LARGER)
{
// Subcurve out of bounds - erase it and continue to next subcurve.
subcurves.erase(iter++);
@@ -1034,52 +1259,54 @@
// Construct the bounding box of the subcurve and compare it to
// the bounding box of the point.
- // Iddo: In the future it might be not a bez_point_bbox
- // but a bbox<Rational>.
- bound_tr.cp_bbox (iter->cp, scv_bbox);
+ bound_tr.construct_bbox (iter->ctrl, scv_bbox);
- if (! _bbox.Overlaps_x(scv_bbox))
+ if (! _bbox.overlaps_x (scv_bbox))
{
// Subcurve out of x bounds - erase it and continue to next subcurve.
subcurves.erase(iter++);
continue;
}
- is_fully_in_t_range = (CGAL::compare (iter->l, t_min) != SMALLER) &&
- (CGAL::compare (iter->r, t_max) != LARGER);
+ is_fully_in_t_range = (CGAL::compare (iter->t_min, t_min) != SMALLER) &&
+ (CGAL::compare (iter->t_max, t_max) != LARGER);
- if (_bbox.Overlaps(scv_bbox) || ! is_fully_in_t_range)
+ if (_bbox.overlaps (scv_bbox) || ! is_fully_in_t_range)
{
- // Iddo: This is a special case of subdividing the curve
+ // \todo This is a special case of subdividing the curve
// not as part of an Originator. Think again if the can_refine
// and de Casteljau should not be different here!!
- can_refine_cv = bound_tr.can_refine (iter->cp, iter->l, iter->r);
+ can_refine_cv = bound_tr.can_refine (iter->ctrl,
+ iter->t_min, iter->t_max);
if (! can_refine_pt && ! can_refine_cv)
+ {
// It is not possible to refine the point or the subcurve anymore:
return (EQUAL);
+ }
if (! can_refine_cv)
{
+ // We are not able to refine the curve. However, we keep it in the
+ // list as we are able to refine the point - so in the future we
+ // can compare the refined point to this curve.
++iter;
continue;
}
- // Subdivide the current subcurve and replace iter with the two
- // resulting subcurves.
- Control_point_vec left_cp, right_cp;
+ // Subdivide the current subcurve and replace it with the two
+ // resulting subcurves (note that we insert the two new subcurves
+ // before iter and remove the current one, pointed by iter).
+ const BoundNT t_mid = (iter->t_min + iter->t_max) / 2;
+ Subcurve scv_left (iter->t_min, t_mid);
+ Subcurve scv_right (t_mid, iter->t_max);
- bisect_control_polygon_2 (iter->cp.begin(), iter->cp.end(),
- std::back_inserter(left_cp),
- std::front_inserter(right_cp));
-
- //bound_tr.DeCasteljau (iter->cp, 0.5, left, right);
-
- // Insert two new subcurves before iter and remove iter.
- BoundNT t_mid = BoundNT(0.5) * (iter->l + iter->r);
+ bisect_control_polygon_2 (iter->ctrl.begin(), iter->ctrl.end(),
+ std::back_inserter(scv_left.ctrl),
+ std::front_inserter(scv_right.ctrl));
- subcurves.insert (iter, Subcurve (left_cp, iter->l, t_mid));
- subcurves.insert (iter, Subcurve (right_cp, t_mid, iter->r));
+ subcurves.insert (iter, scv_left);
+ subcurves.insert (iter, scv_right);
subcurves.erase(iter++);
was_split = true;
@@ -1090,8 +1317,8 @@
// We found a subcurve whose x-range contains our point, but whose
// bounding box is disjoint from the bounding box of the point.
// We can therefore compare the y-positions of the bounding boxes.
- CGAL_assertion (! _bbox.Overlaps (scv_bbox) &&
- _bbox.Overlaps_x (scv_bbox) &&
+ CGAL_assertion (! _bbox.overlaps (scv_bbox) &&
+ _bbox.overlaps_x (scv_bbox) &&
is_fully_in_t_range);
return (CGAL::compare (_bbox.max_y, scv_bbox.max_y));
@@ -1113,8 +1340,10 @@
can_refine_pt = _refine();
if (! can_refine_pt && ! can_refine_cv)
+ {
// It is not possible to refine the point or the subcurve anymore:
return (EQUAL);
+ }
}
// If we reached here, we cannot refine any more:
@@ -1131,29 +1360,49 @@
Bounding_traits bound_tr;
Originator& orig1 = *(_origs.begin());
- if (orig1.point_bound().point_type == Bez_point_bound::VERTICAL_TANGENCY_PT)
+ if (orig1.point_bound().type == Bez_point_bound::VERTICAL_TANGENCY_PT)
{
CGAL_assertion(_origs.size() == 1);
// Refine the vertical tangency point.
- std::pair<Bez_point_bound, Bez_point_bbox> refined_tang_pt;
+ typename Bounding_traits::Vertical_tangency_point vpt (orig1.point_bound(),
+ _bbox);
+ typename Bounding_traits::Vertical_tangency_point ref_vpt;
- bound_tr.refine_tangency_point (orig1.point_bound().bounding_polyline,
- orig1.point_bound().t_min,
- orig1.point_bound().t_max,
- refined_tang_pt);
+ bound_tr.refine_vertical_tangency_point (vpt,
+ ref_vpt);
- if (! refined_tang_pt.first.can_refine)
+ if (! ref_vpt.bound.can_refine)
+ {
// Indicate that it was not possible to refine the point.
return (false);
+ }
// Update the originator and the bounding box of the point.
- orig1.update_point_bound (refined_tang_pt.first);
- _bbox = refined_tang_pt.second;
+ orig1.update_point_bound (ref_vpt.bound);
+ _bbox = ref_vpt.bbox;
+
+ if (ref_vpt.bound.type == Bez_point_bound::RATIONAL_PT)
+ {
+ // In this case the point is exactly computed.
+ Nt_traits nt_traits;
+
+ p_rat_x = new Rational (ref_vpt.bbox.min_x);
+ p_alg_x = new Algebraic (nt_traits.convert (*p_rat_x));
+ p_rat_y = new Rational (ref_vpt.bbox.min_y);
+ p_alg_y = new Algebraic (nt_traits.convert (*p_rat_y));
+
+ orig1.set_parameter (nt_traits.convert (ref_vpt.bound.t_min));
+
+ // Make sure the point is marked as a vertical tangency point.
+ ref_vpt.bound.type = Bez_point_bound::VERTICAL_TANGENCY_PT;
+ orig1.update_point_bound (ref_vpt.bound);
+ }
+
return (true);
}
- if (orig1.point_bound().point_type == Bez_point_bound::INTERSECTION_PT)
+ if (orig1.point_bound().type == Bez_point_bound::INTERSECTION_PT)
{
CGAL_assertion(_origs.size() == 2);
@@ -1163,30 +1412,111 @@
++org_it;
Originator& orig2 = *org_it;
- typename Bounding_traits::Bound_pair intersect_pt (orig1.point_bound(),
- orig2.point_bound(),
- _bbox);
- typename Bounding_traits::Bound_pair refined_pt;
+ typename Bounding_traits::Intersection_point ipt (orig1.point_bound(),
+ orig2.point_bound(),
+ _bbox);
+ typename Bounding_traits::Intersection_point ref_ipt;
- bound_tr.refine_intersection_point (intersect_pt, refined_pt);
+ bound_tr.refine_intersection_point (ipt, ref_ipt);
- if (! refined_pt.bound1.can_refine || ! refined_pt.bound2.can_refine)
+ if (! ref_ipt.bound1.can_refine || ! ref_ipt.bound2.can_refine)
+ {
// Indicate that it was not possible to refine the point.
return (false);
+ }
// Update the originators and the bounding box of the point.
- orig1.update_point_bound (refined_pt.bound1);
- orig2.update_point_bound (refined_pt.bound2);
- _bbox = refined_pt.bbox;
+ orig1.update_point_bound (ref_ipt.bound1);
+ orig2.update_point_bound (ref_ipt.bound2);
+ _bbox = ref_ipt.bbox;
+
+ if (ref_ipt.bound1.type == Bez_point_bound::RATIONAL_PT ||
+ ref_ipt.bound2.type == Bez_point_bound::RATIONAL_PT)
+ {
+ // In this case the point is exactly computed.
+ Nt_traits nt_traits;
+
+ p_rat_x = new Rational (ref_ipt.bbox.min_x);
+ p_alg_x = new Algebraic (nt_traits.convert (*p_rat_x));
+ p_rat_y = new Rational (ref_ipt.bbox.min_y);
+ p_alg_y = new Algebraic (nt_traits.convert (*p_rat_y));
+
+ orig1.set_parameter (nt_traits.convert (ref_ipt.bound1.t_min));
+ orig2.set_parameter (nt_traits.convert (ref_ipt.bound2.t_min));
+
+ // Make sure the point is marked as an intersection point.
+ ref_ipt.bound1.type = Bez_point_bound::INTERSECTION_PT;
+ ref_ipt.bound2.type = Bez_point_bound::INTERSECTION_PT;
+ orig1.update_point_bound (ref_ipt.bound1);
+ orig2.update_point_bound (ref_ipt.bound2);
+ }
+
return (true);
}
- // We should never reach here:
- //CGAL_assertion (false);
+ // If we reached here, the point is computed in an exact manner, so there
+ // is not need to refine its approximation.
return (false);
}
// ---------------------------------------------------------------------------
+// Make sure the originator parameters fit the bound box.
+//
+template <class RatKer, class AlgKer, class NtTrt, class BndTrt>
+void _Bezier_point_2_rep<RatKer, AlgKer, NtTrt, BndTrt>::_fit_to_bbox ()
+{
+ // Get the bounding box of the point.
+ const typename Bounding_traits::NT x_min = _bbox.min_x;
+ const typename Bounding_traits::NT x_max = _bbox.max_x;
+ const typename Bounding_traits::NT y_min = _bbox.min_y;
+ const typename Bounding_traits::NT y_max = _bbox.max_y;
+
+ // Go over all originators, and make sure that the bounding box of the
+ // control polygon of each originator is not larger than the bounding box
+ // of the point.
+ Bounding_traits bound_tr;
+ Bez_point_bbox org_bbox;
+ Orig_iter org_it;
+ bool refined;
+ bool all_fit;
+
+ do
+ {
+ all_fit = true;
+ refined = false;
+
+ for (org_it = _origs.begin(); org_it != _origs.end(); ++org_it)
+ {
+ // Skip rational points.
+ if (org_it->point_bound().type == Bez_point_bound::RATIONAL_PT ||
+ org_it->point_bound().ctrl.empty())
+ {
+ continue;
+ }
+
+ // Get the bounding box of the control polygon.
+ // In case this bounding box is larger than the bounding box
+ // of the point, we refine the point (hence refine the control polygon;
+ // note however we keep the original bounding box of the point, in order
+ // to avoid an infinitive loop here ...)
+ bound_tr.construct_bbox (org_it->point_bound().ctrl,
+ org_bbox);
+
+ if (CGAL::compare (org_bbox.min_x, x_min) == SMALLER ||
+ CGAL::compare (org_bbox.max_x, x_max) == LARGER ||
+ CGAL::compare (org_bbox.min_y, y_min) == SMALLER ||
+ CGAL::compare (org_bbox.max_y, y_max) == LARGER)
+ {
+ all_fit = false;
+ refined |= _refine();
+ }
+ }
+ } while (! all_fit && refined);
+
+ return;
+}
+
+// ---------------------------------------------------------------------------
// Compute the point in an exact manner.
//
template <class RatKer, class AlgKer, class NtTrt, class BndTrt>
@@ -1196,17 +1526,14 @@
if (is_exact())
return;
- /* Ron: For informational purposes ...
- std::cout << "***** MAKE EXACT (" << _origs.size() << " ORIGINATORS) *****"
- << std::endl;
- */
+ // Check if the point is a vertical tangency point of the originator.
+ Nt_traits nt_traits;
- // Check if the point is a vertical tangency point of the originator.
if (_origs.size() == 1)
{
Orig_iter org_it = _origs.begin();
- CGAL_assertion (org_it->point_bound().point_type ==
+ CGAL_assertion (org_it->point_bound().type ==
Bez_point_bound::VERTICAL_TANGENCY_PT);
// Compute (using the cache) the vertical tangency parameters of
@@ -1218,8 +1545,8 @@
typename Bezier_cache::Vertical_tangency_iter vt_it;
// Look for a parameter within the range of the bounding interval.
- const Algebraic t_min (org_it->point_bound().t_min);
- const Algebraic t_max (org_it->point_bound().t_max);
+ const Algebraic& t_min = nt_traits.convert (org_it->point_bound().t_min);
+ const Algebraic& t_max = nt_traits.convert (org_it->point_bound().t_max);
for (vt_it = vt_list.begin(); vt_it != vt_list.end(); ++vt_it)
{
@@ -1295,7 +1622,6 @@
CGAL_assertion (! do_ovlp);
// Look for a parameter pair within the ranges of the bounding intervals.
- Nt_traits nt_traits;
const Algebraic s_min = nt_traits.convert (orig1.point_bound().t_min);
const Algebraic s_max = nt_traits.convert (orig1.point_bound().t_max);
const Algebraic t_min = nt_traits.convert (orig2.point_bound().t_min);
Index: Arr_traits_2/Bezier_bounding_rational_traits.h
===================================================================
--- Arr_traits_2/Bezier_bounding_rational_traits.h (revision 40349)
+++ Arr_traits_2/Bezier_bounding_rational_traits.h (working copy)
@@ -16,6 +16,7 @@
//
//
// Author(s) : Iddo Hanniel <>
+// Ron Wein <>
#ifndef CGAL_BEZIER_BOUNDING_RATIONAL_TRAITS_H
#define CGAL_BEZIER_BOUNDING_RATIONAL_TRAITS_H
@@ -35,12 +36,129 @@
CGAL_BEGIN_NAMESPACE
-template <class _NT>
-struct _Bez_point_bound ;
+/*! \struct _Bez_point_bound
+ * Representation of a bounding interval for a point on a Bezier curve.
+ * Basically, we store an interval [t_min, t_max] such that the point p
+ * equals B(t) for some point t_min <= t <= t_max. We also store a bounding
+ * polynomial for the point (obtain by de Casteljau's algorithm).
+ */
+template <class Kernel_>
+struct _Bez_point_bound
+{
+ /*! \enum Type
+ * The point type.
+ */
+ enum Type
+ {
+ RATIONAL_PT, /*!< A point with rational coordinates .*/
+ VERTICAL_TANGENCY_PT, /*!< A vertical tangency point .*/
+ INTERSECTION_PT, /*!< An intersection point .*/
+ UNDEFINED
+ };
-template <class _NT>
-struct _Bez_point_bbox ;
+ typedef Kernel_ Kernel;
+ typedef typename Kernel::FT NT;
+ typedef typename Kernel::Point_2 Point_2;
+ typedef std::deque<Point_2> Control_points;
+ Type type; /*!< The point type. */
+ Control_points ctrl; /*!< The control point whose convex
+ hull contains the point. */
+ NT t_min; /*!< Minimal t value. */
+ NT t_max; /*!< Maximal t value. */
+ bool can_refine; /*!< Can we refine the current
+ representation. */
+
+ /*! Default constructor. */
+ _Bez_point_bound() :
+ type (UNDEFINED),
+ ctrl(),
+ t_min(), t_max(),
+ can_refine (false)
+ {}
+
+ /*! Constructor with parameters. */
+ _Bez_point_bound (Type _type,
+ const Control_points& _ctrl,
+ const NT& _t_min, const NT& _t_max,
+ bool _can_refine) :
+ type (_type),
+ ctrl (_ctrl),
+ t_min (_t_min), t_max (_t_max),
+ can_refine (_can_refine)
+ {}
+};
+
+/*! \struct _Bez_point_bbox
+ * A bounding box for a point on a Bezier curve.
+ */
+template <class Kernel_>
+struct _Bez_point_bbox
+{
+ typedef Kernel_ Kernel;
+ typedef typename Kernel::FT NT;
+ typedef _Bez_point_bbox<Kernel> Self;
+
+ NT min_x; /*!< Lower bound for the x-coordinate. */
+ NT max_x; /*!< Upper bound for the x-coordinate. */
+ NT min_y; /*!< Lower bound for the y-coordinate. */
+ NT max_y; /*!< Upper bound for the y-coordinate. */
+
+ /*! Default constructor. */
+ _Bez_point_bbox() :
+ min_x(0), max_x(0),
+ min_y(0), max_y(0)
+ {}
+
+ /*! Constructor with parameters. */
+ _Bez_point_bbox (const NT& _min_x, const NT& _max_x,
+ const NT& _min_y, const NT& _max_y) :
+ min_x(_min_x), max_x(_max_x),
+ min_y(_min_y), max_y(_max_y)
+ {}
+
+ /*! Add two bounding boxes, and obtain a box bounding them both. */
+ Self operator+ (const Self& other) const
+ {
+ return (Self (std::min(min_x, other.min_x), std::max(max_x, other.max_x),
+ std::min(min_y, other.min_y), std::max(max_y, other.max_y)));
+ }
+
+ /*! Addition and assignment. */
+ void operator+= (const Self& other)
+ {
+ min_x = std::min(min_x, other.min_x);
+ max_x = std::max(max_x, other.max_x);
+ min_y = std::min(min_y, other.min_y);
+ max_y = std::max(max_y, other.max_y);
+ return;
+ }
+
+ /*! Check whether the two bounding boxed overlap. */
+ bool overlaps (const Self& other) const
+ {
+ if ((CGAL::compare (max_x, other.min_x) == CGAL::SMALLER) ||
+ (CGAL::compare (min_x, other.max_x) == CGAL::LARGER) ||
+ (CGAL::compare (max_y, other.min_y) == CGAL::SMALLER) ||
+ (CGAL::compare (min_y, other.max_y) == CGAL::LARGER))
+ {
+ return (false);
+ }
+ return (true);
+ }
+
+ /*! Check whether the two bounding boxed overlap in their x-range. */
+ bool overlaps_x (const Self& other) const
+ {
+ if ((CGAL::compare (max_x, other.min_x) == CGAL::SMALLER) ||
+ (CGAL::compare (min_x, other.max_x) == CGAL::LARGER))
+ {
+ return (false);
+ }
+ return (true);
+ }
+};
+
/*!\ class Bezier_bounding_rational_traits
* A traits class for performing bounding operations on Bezier curves,
* assuming the NT is rational.
@@ -49,1134 +167,1347 @@
class Bezier_bounding_rational_traits
{
public:
- typedef _Kernel Kernel; // Used in polygon and point/vector operations
- typedef typename Kernel::FT NT; // The number type used to represent control points and t-parameters.
- typedef typename Kernel::Point_2 Point_2; // The control point type.
- typedef typename Kernel::Vector_2 Vector_2;
- typedef typename Kernel::Line_2 Line_2;
- typedef std::deque<Point_2> Control_point_vec;
- typedef _Bez_point_bound<NT> Bez_point_bound;
- typedef _Bez_point_bbox<NT> Bez_point_bbox;
+ typedef _Kernel Kernel;
+ typedef typename Kernel::FT NT;
- // Auxilary structure for output of intersection points.
- struct Bound_pair
+ typedef _Bez_point_bound<Kernel> Bez_point_bound;
+ typedef _Bez_point_bbox<Kernel> Bez_point_bbox;
+ typedef typename Bez_point_bound::Control_points Control_points;
+
+ /*! \struct Vertical_tangency_point
+ * Representation of an approximated vertical tangency point.
+ */
+ struct Vertical_tangency_point
{
- Bez_point_bound bound1;
- Bez_point_bound bound2;
- Bez_point_bbox bbox;
-
- Bound_pair(const Bez_point_bound& b1, const Bez_point_bound& b2,
- const Bez_point_bbox& bb) : bound1(b1), bound2(b2), bbox(bb) {}
- Bound_pair() : bound1(), bound2(), bbox() {}
+ Bez_point_bound bound; /*!< A point-on-curve bound. */
+ Bez_point_bbox bbox; /*!< Bounding box for the point. */
+
+ /*! Default constructor. */
+ Vertical_tangency_point () :
+ bound (),
+ bbox ()
+ {}
+
+ /*! Constructor. */
+ Vertical_tangency_point (const Bez_point_bound& _bound,
+ const Bez_point_bbox& _bbox) :
+ bound (_bound),
+ bbox (_bbox)
+ {}
};
- typedef std::list<Bound_pair> Bound_pair_lst;
+ /*! \struct Vertical_tangency_point
+ * Representation of an approximated intersection point.
+ */
+ struct Intersection_point
+ {
+ Bez_point_bound bound1; /*!< A point-on-curve bound for the
+ first curve. */
+ Bez_point_bound bound2; /*!< A point-on-curve bound for the
+ second curve. */
+ Bez_point_bbox bbox; /*!< Bounding box for the point. */
- //TODO - add here typedefs for vertical tangency output (pair<..>).
+ /*! Default constructor. */
+ Intersection_point() :
+ bound1(),
+ bound2(),
+ bbox()
+ {}
- Kernel kernel;
- NT _can_refine_bound;
+ /*! Constructor. */
+ Intersection_point (const Bez_point_bound& _bound1,
+ const Bez_point_bound& _bound2,
+ const Bez_point_bbox& _bbox) :
+ bound1 (_bound1),
+ bound2 (_bound2),
+ bbox (_bbox)
+ {}
- // a ctr that enables to change can_refine_bound.
- // iddo: important! we need to make the arr traits receive this traits as a reference,
- // since it now contains state.
- // The line below currently causes a runtime bug, so for now.
- //Bezier_bounding_rational_traits(double can_refine_bound = std::pow(2,-53))
- //0.00000000000000000000001)
- Bezier_bounding_rational_traits (double can_refine_bound =
- 0.00000000000000011)
- : _can_refine_bound(can_refine_bound)
- {}
-
- // iddo: can_refine() function can be based on Euclidean or parametric space.
- // For rationals it is a simple parametric-space check.
- bool can_refine(const Control_point_vec& , const NT& l, const NT& r)
+ };
+
+protected:
+
+ typedef typename Kernel::Point_2 Point_2;
+ typedef typename Kernel::Vector_2 Vector_2;
+ typedef typename Kernel::Line_2 Line_2;
+ typedef typename Kernel::Direction_2 Direction_2;
+
+ // Data members:
+ NT m_accuracy_bound; /*!< An accuracy bound. */
+ Kernel m_kernel; /*!< A geometry kernel. */
+ unsigned int m_active_nodes; /*!< Limit for the number of recursive
+ calls for the approximated
+ intersection procedure. */
+
+ // Kernel functors:
+ typename Kernel::Construct_vector_2 f_construct_vector;
+ typename Kernel::Construct_line_2 f_construct_line;
+ typename Kernel::Construct_direction_2 f_construct_direction;
+ typename Kernel::Construct_opposite_direction_2 f_opposite_direction;
+ typename Kernel::Construct_translated_point_2 f_construct_point;
+ typename Kernel::Equal_2 f_equal;
+ typename Kernel::Compare_x_2 f_compare_x;
+ typename Kernel::Compare_y_2 f_compare_y;
+ typename Kernel::Orientation_2 f_orientation;
+ typename Kernel::Counterclockwise_in_between_2 f_ccw_in_between;
+ typename Kernel::Less_signed_distance_to_line_2 f_less_signed_distance;
+ typename Kernel::Oriented_side_2 f_oriented_side;
+ typename Kernel::Intersect_2 f_intersect;
+
+public:
+
+ /*!
+ * Constructor.
+ * \param bound Accuracy bound. We cannot determine the order of two values
+ * whose absolute difference is smaller than this bound.
+ */
+ Bezier_bounding_rational_traits (double bound = 0.00000000000000011) :
+ m_accuracy_bound (bound),
+ m_active_nodes (0)
{
- // iddo: might be made more efficient if needed (based on denom...)
- if (r-l < _can_refine_bound)
- return false;
- return true;
+ // Construct the kernel functors.
+ f_construct_vector = m_kernel.construct_vector_2_object();
+ f_construct_line = m_kernel.construct_line_2_object();
+ f_construct_direction = m_kernel.construct_direction_2_object();
+ f_opposite_direction = m_kernel.construct_opposite_direction_2_object();
+ f_construct_point = m_kernel.construct_translated_point_2_object();
+ f_equal = m_kernel.equal_2_object();
+ f_compare_x = m_kernel.compare_x_2_object();
+ f_compare_y = m_kernel.compare_y_2_object();
+ f_orientation = m_kernel.orientation_2_object();
+ f_ccw_in_between = m_kernel.counterclockwise_in_between_2_object();
+ f_less_signed_distance = m_kernel.less_signed_distance_to_line_2_object();
+ f_oriented_side = m_kernel.oriented_side_2_object();
+ f_intersect = m_kernel.intersect_2_object();
}
- // Can be made a template function <Control_point_vec, NT>
- /*
- void DeCasteljau(const Control_point_vec& cp, const NT& t0,
- Control_point_vec& lcp, Control_point_vec& rcp) const
+ /// \name Handling intersection points.
+ //@{
+
+ /*!
+ * Check whether the curve specified by the given control polygon may
+ * have self intersections.
+ * \param cp The control polygon of the curve.
+ * \return Whether the curve may be self-intersecting (but not necessarily).
+ * In case the function returns false, there can be absolutely no
+ * self-intersections.
+ */
+ bool may_have_self_intersections (const Control_points& cp)
{
- const NT t1 = NT(1) - t0;
+ // In case the control polygon is convex, the Bezier curve cannot be
+ // self-intersecting.
+ if (is_convex_2 (cp.begin(), cp.end(), m_kernel))
+ return (false);
- lcp.push_back(cp.front());
- rcp.push_front(cp.back());
+ // In case the angular span of the control polygon is less than 180
+ // degrees, it cannot be self-intersecting.
+ Vector_2 v_min;
+ Vector_2 v_max;
- Control_point_vec aux(cp.begin(), cp.end());
- unsigned int i=0;
- for (; i<cp.size()-1; ++i)
- {
- Control_point_vec aux1(cp.size()-i-1); //in every step of the subdivision we reduce by 1 the size of the auxilary vector
- unsigned int j=0;
- for (; j<cp.size()-i-1; ++j)
- {
- aux1[j] = Point_2(t1*(aux[j].x()) + t0*(aux[j+1].x()), t1*(aux[j].y()) + t0*(aux[j+1].y()));
- }
- lcp.push_back(aux1.front());
- rcp.push_front(aux1.back());
+ if (_compute_angular_span (cp, v_min, v_max))
+ return (false);
- aux = aux1;
- }
+ // Otherwise, the curve may be self-intersecting.
+ return (true);
}
- */
- ////////////////////////////////////////////////////////////////
- //Intersection points
- ///////////////////////////////////////////////////////////////
- void CrvCrvInter(const Control_point_vec& cp1, const Control_point_vec& cp2,
- Bound_pair_lst& intersection_pairs)
+ /*!
+ * Compute bounds for the intersection points between two Bezier curves.
+ * \param cp1 The control polygon of the first curve.
+ * \param cp2 The control polygon of the second curve.
+ * \param oi Output: The bounds of the intersection points.
+ * The value-type of this iterator is Intersection_point.
+ */
+ template <class OutputIterator>
+ OutputIterator compute_intersection_points (const Control_points& cp1,
+ const Control_points& cp2,
+ OutputIterator oi)
{
- //iddo: in order to support interval, we will need to convert Rational to
- // interval do we need to do it outside the call to the function or
- // here (needs conversion from Rat_Point...).
-
- std::set<NT> endpoint_intersections; // Auxilary variable needed in the recursion
- // to identify already existing points due to endpoint intersections.
+ // Call the recursive function.
+ std::set<NT> dummy;
+ std::list<Intersection_point> ipts;
- // Call to recursive function.
- _CrvCrvInterAux (cp1, 0, 1,
- cp2, 0, 1,
- true,
- endpoint_intersections,
- intersection_pairs);
+ m_active_nodes = 1;
+ _compute_intersection_points (cp1, 0, 1,
+ cp2, 0, 1,
+ true, dummy, // Check span.
+ ipts);
+
+ // Copy the computed points to the output iterator.
+ typename std::list<Intersection_point>::iterator it;
+
+ for (it = ipts.begin(); it != ipts.end(); ++it)
+ {
+ *oi = *it;
+ ++oi;
+ }
+
+ return (oi);
}
- void refine_intersection_point (const Bound_pair& intersection_point,
- Bound_pair& refined_point)
+ /*!
+ * Refine the approximation of a given intersection point.
+ * \param in_pt The intersection point.
+ * \param ref_pt Output: The refined representation of the point.
+ * \return Whether the point has been successfully refined (if not, the
+ * output point is the same as the input).
+ */
+ bool refine_intersection_point (const Intersection_point& in_pt,
+ Intersection_point& ref_pt)
{
- // Assumes we do not need to check span, intersection_point has all the
- //necessary data.
- CGAL_precondition(intersection_point.bound1.can_refine);
- CGAL_precondition(intersection_point.bound2.can_refine);
- // The following precondition makes sure the point is not rational already
- // (in which case it would not need refinement).
- //CGAL_precondition(intersection_point.bbox.min_x != intersection_point.bbox.max_x);
+ // Check if it is not possible to refine.
+ if (! in_pt.bound1.can_refine || ! in_pt.bound2.can_refine)
+ {
+ ref_pt = in_pt;
+ return (false);
+ }
- const Control_point_vec& cv1 = intersection_point.bound1.bounding_polyline;
- const NT& l1 = intersection_point.bound1.t_min;
- const NT& r1 = intersection_point.bound1.t_max;
- const Control_point_vec& cv2 = intersection_point.bound2.bounding_polyline;
- const NT& l2 = intersection_point.bound2.t_min;
- const NT& r2 = intersection_point.bound2.t_max;
+ // In case the point is already rational, there is not point in refining
+ // it.
+ if (in_pt.bound1.type == Bez_point_bound::RATIONAL_PT ||
+ in_pt.bound2.type == Bez_point_bound::RATIONAL_PT)
+ {
+ ref_pt = in_pt;
+ return (false);
+ }
- /* Subdivide the two curves and recurse. */
- // This is probably not needed as it will be handled
- // in the recursion.
- bool can_refine1 = can_refine(cv1, l1, r1);
- bool can_refine2 = can_refine(cv2, l2, r2);
- if (!can_refine1 || !can_refine2)
+ // Check whether we can refine the parametric ranges of the two
+ // originating Bezier curves.
+ const Control_points& cp1 = in_pt.bound1.ctrl;
+ const NT& t_min1 = in_pt.bound1.t_min;
+ const NT& t_max1 = in_pt.bound1.t_max;
+ bool can_refine1 = can_refine (cp1, t_min1, t_max1);
+
+ const Control_points& cp2 = in_pt.bound2.ctrl;
+ const NT& t_min2 = in_pt.bound2.t_min;
+ const NT& t_max2 = in_pt.bound2.t_max;
+ bool can_refine2 = can_refine (cp2, t_min2, t_max2);
+
+ if (! can_refine1 || ! can_refine2)
{
// Failed in the bounded approach - stop the subdivision and indicate
// that the subdivision has failed.
- refined_point = intersection_point;
- refined_point.bound1.can_refine = false;
- refined_point.bound2.can_refine = false;
+ ref_pt = in_pt;
+ ref_pt.bound1.can_refine = false;
+ ref_pt.bound2.can_refine = false;
- return;
+ return (false);
}
- Control_point_vec cv1a, cv1b;
-
- bisect_control_polygon_2 (cv1.begin(), cv1.end(),
- std::back_inserter(cv1a),
- std::front_inserter(cv1b));
- //DeCasteljau(cv1, 0.5, cv1a, cv1b);
- NT t_mid1 = NT(0.5) * (l1 + r1);
+ // Apply de Casteljau's algorithm and bisect both bounding polylines.
+ Control_points cp1a, cp1b;
+ const NT t_mid1 = (t_min1 + t_max1) / 2;
- Control_point_vec cv2a, cv2b;
+ bisect_control_polygon_2 (cp1.begin(), cp1.end(),
+ std::back_inserter(cp1a),
+ std::front_inserter(cp1b));
- bisect_control_polygon_2 (cv2.begin(), cv2.end(),
- std::back_inserter(cv2a),
- std::front_inserter(cv2b));
- //DeCasteljau(cv2, 0.5, cv2a, cv2b);
- NT t_mid2 = NT(0.5) * (l2 + r2);
+ Control_points cp2a, cp2b;
+ const NT t_mid2 = (t_min2 + t_max2) / 2;
- // Iddo: Catch situations that we cannot refine, in roder to prevent
- // having multiple points as the refined point.
- bool can_refine_pt_cv1a = can_refine (cv1a, l1, t_mid1);
- bool can_refine_pt_cv1b = can_refine (cv1b, t_mid1, r1);
- bool can_refine_pt_cv2a = can_refine (cv2a, l2, t_mid2);
- bool can_refine_pt_cv2b = can_refine (cv2b, t_mid2, r2);
+ bisect_control_polygon_2 (cp2.begin(), cp2.end(),
+ std::back_inserter(cp2a),
+ std::front_inserter(cp2b));
+
+ // Catch situations that we cannot further refine, in order to prevent
+ // having multiple points as the refined point.
+ const bool can_refine_pt_cp1a = can_refine (cp1a, t_min1, t_mid1);
+ const bool can_refine_pt_cp1b = can_refine (cp1b, t_mid1, t_max1);
+ const bool can_refine_pt_cp2a = can_refine (cp2a, t_min2, t_mid2);
+ const bool can_refine_pt_cp2b = can_refine (cp2b, t_mid2, t_max2);
- if (!can_refine_pt_cv1a || !can_refine_pt_cv1b ||
- !can_refine_pt_cv2a || !can_refine_pt_cv2b)
+ if (! can_refine_pt_cp1a || ! can_refine_pt_cp1b ||
+ ! can_refine_pt_cp2a || ! can_refine_pt_cp2b)
{
- // Construct a Bez_point_bound that is inconclusive:
- // all parameters found so far + can_refine = false, with the bounding
- // box of cv1, and add to output.
- Bez_point_bbox bbox1;
+ // Construct an inconclusive point bound, which includes all parameters
+ // found so far, with can_refine = false, and with the bounding
+ // box of cp1.
+ Bez_point_bound bound1 (Bez_point_bound::INTERSECTION_PT,
+ cp1, t_min1, t_max1,
+ false); // Cannot refine further.
+ Bez_point_bound bound2 (Bez_point_bound::INTERSECTION_PT,
+ cp2, t_min2, t_max2,
+ false); // Cannot refine further.
+ Bez_point_bbox ref_bbox;
- cp_bbox (cv1, bbox1);
- refined_point =
- Bound_pair (Bez_point_bound (cv1, l1, r1,
- Bez_point_bound::INTERSECTION_PT, false),
- Bez_point_bound (cv2, l2, r2,
- Bez_point_bound::INTERSECTION_PT, false),
- bbox1);
+ construct_bbox (cp1, ref_bbox);
+ ref_pt = Intersection_point (bound1, bound2, ref_bbox);
+ return (true);
+ }
- return;
+ // Use the bisection method to refine the intersection point.
+ // We assume that there is no need to check the span, as the input point
+ // already has all the necessary data.
+ std::list<Intersection_point> ipts;
+ std::set<NT> dummy;
+
+ m_active_nodes = 1;
+ _compute_intersection_points (cp1a, t_min1, t_mid1,
+ cp2a, t_min2, t_mid2,
+ false, dummy, // Don't check span.
+ ipts);
+ if (! ipts.empty())
+ {
+ if (ipts.size() == 1)
+ {
+ ref_pt = ipts.front();
+ return (true);
+ }
+ ipts.clear();
}
- //Recursion:
- Bound_pair_lst intersection_pairs;
- std::set<NT> endpoint_intersections;
- bool should_check_span = false; //this precondition is hard to check efficiently.
+ m_active_nodes = 1;
+ _compute_intersection_points (cp1a, t_min1, t_mid1,
+ cp2b, t_mid2, t_max2,
+ false, dummy, // Don't check span.
+ ipts);
+ if (! ipts.empty())
+ {
+ if (ipts.size() == 1)
+ {
+ ref_pt = ipts.front();
+ return (true);
+ }
+ ipts.clear();
+ }
- _CrvCrvInterAux (cv1a, l1, t_mid1, cv2a, l2, t_mid2,
- should_check_span, endpoint_intersections,
- intersection_pairs);
- _CrvCrvInterAux (cv1a, l1, t_mid1, cv2b, t_mid2, r2,
- should_check_span, endpoint_intersections,
- intersection_pairs);
- _CrvCrvInterAux (cv1b, t_mid1, r1, cv2a, l2, t_mid2,
- should_check_span, endpoint_intersections,
- intersection_pairs);
- _CrvCrvInterAux (cv1b, t_mid1, r1, cv2b, t_mid2, r2,
- should_check_span, endpoint_intersections,
- intersection_pairs);
+ m_active_nodes = 1;
+ _compute_intersection_points (cp1b, t_mid1, t_max1,
+ cp2a, t_min2, t_mid2,
+ false, dummy, // Don't check span.
+ ipts);
+ if (! ipts.empty())
+ {
+ if (ipts.size() == 1)
+ {
+ ref_pt = ipts.front();
+ return (true);
+ }
+ ipts.clear();
+ }
- refined_point = intersection_pairs.front();
+ m_active_nodes = 1;
+ _compute_intersection_points (cp1b, t_mid1, t_max1,
+ cp2b, t_mid2, t_max2,
+ false, dummy, // Don't check span.
+ ipts);
- CGAL_assertion (intersection_pairs.size() == 1);
-
- return;
+ CGAL_assertion (ipts.size() == 1);
+ ref_pt = ipts.front();
+ return (true);
}
+ //@}
- ////////////////////////////////////////////////////////////////
- // Vertical tangency points
- ///////////////////////////////////////////////////////////////
+ /// \name Handling vertcial tangency points.
+ //@{
- // Vertical tangency points (cp is a Bez curve between t0-t1).
- // Output are pairs of <bound,bbox> representing the vertical tangency points.
- void vertical_tangency_points
- (const Control_point_vec& cp,
- const NT& t0, const NT& t1,
- std::list<std::pair<Bez_point_bound, Bez_point_bbox> >& tangency_points)
+ /*!
+ * Compute bounds for the vertical tangency points of a Bezier curves.
+ * The points are returned sorted by their parameter value.
+ * \param cp The control polygon of the curve.
+ * \param oi Output: The bounds of the vertical tangency points.
+ * The value-type is Vertical_tangency_point.
+ */
+ template <class OutputIterator>
+ OutputIterator
+ compute_vertical_tangency_points (const Control_points& cp,
+ OutputIterator oi)
{
- //TODO - handle special case of degree two curves.
+ // Call the recursive function on the entire curve (0 <= t <= 1).
+ std::list<Vertical_tangency_point> vpts;
- bool can_refine_pt = can_refine(cp, t0, t1);
- if (!can_refine_pt)
- {
- // Construct a Bez_point_bound that is inconclusive:
- // all parameters found so far + can_refine = false,
- // and add to output.
- Bez_point_bbox bbox;
- cp_bbox(cp, bbox);
- tangency_points.push_back
- (std::make_pair(Bez_point_bound(cp, t0, t1,
- Bez_point_bound::VERTICAL_TANGENCY_PT,
- false),
- bbox));
- return;
- }
+ _compute_vertical_tangency_points (cp,
+ 0, 1,
+ vpts);
- // Check for x-monoticity (if exists return)
- bool is_x_monotone = _is_x_monotone(cp);
- if (is_x_monotone)
- return;
+ // Copy the computed points to the output iterator, sorted by their
+ // parameter value (note we compare the t_min values of the bounds; this
+ // is possible as different vertical tangency points are expected to
+ // have disjoint parametric ranges).
+ typename std::list<Vertical_tangency_point>::iterator vpt_it;
+ typename std::list<Vertical_tangency_point>::iterator vpt_end;
+ typename std::list<Vertical_tangency_point>::iterator vpt_min;
- // Check for y-monoticity (else recurse)
- bool is_y_monotone = _is_y_monotone(cp);
- if (!is_y_monotone)
+ while (! vpts.empty())
{
- // Recurse
- Control_point_vec left, right;
+ // Locate the vertical tangency point with minimal t-value.
+ vpt_min = vpt_it = vpts.begin();
+ vpt_end = vpts.end();
- bisect_control_polygon_2 (cp.begin(), cp.end(),
- std::back_inserter(left),
- std::front_inserter(right));
- //DeCasteljau(cp, 0.5, left, right);
- NT t_mid = NT(0.5) * (t0 + t1);
- // TODO - handle the case where t_mid is a vertical (rational) tangency
- // point, by checking whether the first vector of right is vertical.
- // Ron: This happens in the Bezier.dat input file!
- CGAL_assertion(CGAL::compare(right[0].x(), right[1].x()) != EQUAL);
+ ++vpt_it;
+ while (vpt_it != vpt_end)
+ {
+ if (CGAL::compare (vpt_it->bound.t_min,
+ vpt_min->bound.t_min) == SMALLER)
+ {
+ vpt_min = vpt_it;
+ }
- vertical_tangency_points(left, t0, t_mid, tangency_points);
- vertical_tangency_points(right, t_mid, t1, tangency_points);
- return;
+ ++vpt_it;
+ }
+
+ // Copy this point to the output iterator.
+ *oi = *vpt_min;
+ ++oi;
+
+ // Remove it from the list.
+ vpts.erase (vpt_min);
}
- // Check for convexity (since it is y-monotone and not x-monotone)
- // if convex - add to list, else recurse.
- if (is_convex_2 (cp.begin(), cp.end(), kernel))
+ return (oi);
+ }
+
+ /*!
+ * Refine the approximation of a given vertical tangency point.
+ * \param in_pt The vertical tangency point.
+ * \param ref_pt Output: The refined representation of the point.
+ * \return Whether the point has been successfully refined (if not, the
+ * output point is the same as the input).
+ */
+ bool refine_vertical_tangency_point (const Vertical_tangency_point& in_pt,
+ Vertical_tangency_point& ref_pt)
+ {
+ // Check if it is not possible to refine.
+ if (! in_pt.bound.can_refine ||
+ in_pt.bound.type == Bez_point_bound::RATIONAL_PT)
{
- Bez_point_bbox bbox;
- typename Control_point_vec::const_iterator res;
- res = bottom_vertex_2 (cp.begin(), cp.end(), kernel);
- bbox.min_y = res -> y();
- res = top_vertex_2 (cp.begin(), cp.end(), kernel);
- bbox.max_y = res -> y();
- res = left_vertex_2 (cp.begin(), cp.end(), kernel);
- bbox.min_x = res -> x();
- res = right_vertex_2 (cp.begin(), cp.end(), kernel);
- bbox.max_x = res -> x();
-
- tangency_points.push_back(std::make_pair(
- Bez_point_bound(cp, t0, t1, Bez_point_bound::VERTICAL_TANGENCY_PT, true),
- bbox));
+ ref_pt = in_pt;
+ return (false);
}
- else
- {
- // Recurse
- Control_point_vec left, right;
- bisect_control_polygon_2 (cp.begin(), cp.end(),
- std::back_inserter(left),
- std::front_inserter(right));
- //DeCasteljau(cp, 0.5, left, right);
- NT t_mid = NT(0.5) * (t0 + t1);
+ // Check whether we can refine the parametric ranges of the originating
+ // Bezier curve.
+ const Control_points& cp = in_pt.bound.ctrl;
+ const NT& t_min = in_pt.bound.t_min;
+ const NT& t_max = in_pt.bound.t_max;
- // TODO - handle the case where t_mid is a vertical (rational) tangency
- // point, by checking whether the first vector of right is vertical.
- CGAL_assertion(CGAL::compare(right[0].x(), right[1].x()) != EQUAL);
+ if (! can_refine (cp, t_min, t_max))
+ {
+ // Failed in the bounded approach - stop the subdivision and indicate
+ // that the subdivision has failed.
+ ref_pt = in_pt;
+ ref_pt.bound.can_refine = false;
- vertical_tangency_points(left, t0, t_mid, tangency_points);
- vertical_tangency_points(right, t_mid, t1, tangency_points);
- return;
+ return (false);
}
- // Debug print
- //std::cout << "number of tangency points is " << tangency_points.size() << std::endl;
- }
+ // Apply de Casteljau's algorithm and bisect the bounding polyline.
+ Control_points cp_a, cp_b;
+ const NT t_mid = (t_min + t_max) / 2;
- // Refinement of vertical tangency point.
- // Assumes the tangency point is isolated (only a single point in the interval).
- void refine_tangency_point(const Control_point_vec& cp, const NT& t0, const NT& t1,
- std::pair<Bez_point_bound, Bez_point_bbox>& tangency_point)
- {
- // TODO - a better implmentation that just checks which of the two subdivided
- // curves is x-monotone (since they should both be y-monotone and convex).
+ bisect_control_polygon_2 (cp.begin(), cp.end(),
+ std::back_inserter(cp_a),
+ std::front_inserter(cp_b));
- Control_point_vec left, right;
- /*
- // iddo: the can_refine check is probably not needed here (the check
- // in the recursion will handle it).
- bool can_refine_pt = can_refine(cp, left, right);
- if (!can_refine)
+ // Handle the case where t_mid is a vertical (rational) tangency
+ // point, by checking whether the first vector of right is vertical.
+ if (f_compare_x (cp_b[0], cp_b[1]) == EQUAL)
{
- CGAL_assertion(false);
- }
- */
+ // The subdivision point at t_mid is a vertical tangency point with
+ // rational coordinates.
+ const typename Kernel::Point_2 &vpt = cp_b[0];
- bisect_control_polygon_2 (cp.begin(), cp.end(),
- std::back_inserter(left),
- std::front_inserter(right));
- //DeCasteljau(cp, 0.5, left, right);
- NT t_mid = NT(0.5) * (t0 + t1);
+ Bez_point_bound bound (Bez_point_bound::RATIONAL_PT,
+ cp_b, t_mid, t_mid, true);
+ Bez_point_bbox bbox (vpt.x(), vpt.x(), vpt.y(), vpt.y());
- // TODO - handle the case where t_mid is a vertical (rational) tangency point,
- // by checking whether the first vector of right is vertical.
- CGAL_assertion(CGAL::compare(right[0].x(), right[1].x()) != EQUAL);
+ ref_pt = Vertical_tangency_point (bound, bbox);
+ return (true);
+ }
- bool can_refine_pt_left = can_refine(left, t0, t_mid);
- bool can_refine_pt_right = can_refine(right, t_mid, t1);
+ const bool can_refine_pt_cp_a = can_refine (cp_a, t_min, t_mid);
+ const bool can_refine_pt_cp_b = can_refine (cp_b, t_mid, t_max);
- if (!can_refine_pt_left || !can_refine_pt_right)
+ if (! can_refine_pt_cp_a || ! can_refine_pt_cp_b)
{
- // Construct a Bez_point_bound that is inconclusive:
- // all parameters found so far + can_refine = false,
- // and add to output.
- Bez_point_bbox bbox;
- cp_bbox(cp, bbox);
- tangency_point =
- std::make_pair(Bez_point_bound(cp, t0, t1,
- Bez_point_bound::VERTICAL_TANGENCY_PT,
- false),
- bbox);
- return;
+ // Construct an inconclusive point bound, which includes all parameters
+ // found so far, with can_refine = false, and with the bounding
+ // box of the curve.
+ Bez_point_bound bound (Bez_point_bound::VERTICAL_TANGENCY_PT,
+ cp, t_min, t_max,
+ false); // Cannot refine further.
+ Bez_point_bbox ref_bbox;
+
+ construct_bbox (cp, ref_bbox);
+ ref_pt = Vertical_tangency_point (bound, ref_bbox);
+ return (true);
}
- std::list<std::pair<Bez_point_bound, Bez_point_bbox> > aux1, aux2;
- vertical_tangency_points(left, t0, t_mid, aux1);
- vertical_tangency_points(right, t_mid, t1, aux2);
+ // Compute the vertical tangency points of the two subcurves in order to
+ // refine the vertical tangency point.
+ std::list<Vertical_tangency_point> vpts;
- CGAL_assertion(aux1.size() + aux2.size() == 1);
- if (aux1.size() == 1)
- {
- tangency_point = *(aux1.begin());
- }
- else
- {
- tangency_point = *(aux2.begin());
- }
+ _compute_vertical_tangency_points(cp_a, t_min, t_mid, vpts);
+ _compute_vertical_tangency_points(cp_b, t_mid, t_max, vpts);
+
+ CGAL_assertion(vpts.size() == 1);
+
+ ref_pt = vpts.front();
+ return (true);
}
+ //@}
- ////////////////////////////////////////////////////////////////
- // Slope comparison.
- ///////////////////////////////////////////////////////////////
- Comparison_result compare_slopes_of_bounded_intersection_point(
- const Bez_point_bound& bound1, const Bez_point_bound& bound2)
+ /// \name Predicates.
+ //@{
+
+ /*!
+ * Check whether the current polyline can be refined.
+ * In this case we simply check the values that define the parametric
+ * range are well-separated.
+ */
+ bool can_refine (const Control_points& ,
+ const NT& t_min, const NT& t_max)
{
+ return (CGAL::compare (t_max - t_min, m_accuracy_bound) != SMALLER);
+ }
- // Since we assume the spans do not overlap, we can just compare any vector
- // of the hodograph spans.
- const Control_point_vec& cv1 = bound1.bounding_polyline;
- const Control_point_vec& cv2 = bound2.bounding_polyline;
+ /*!
+ * Compare the slopes of two Bezier curve at their intersection point,
+ * given the point-on-curve bounds at this point.
+ * \param bound1 The point-on-curve bound of the first curve.
+ * \param bound2 The point-on-curve bound of the second curve.
+ * \return The comparison result.
+ */
+ Comparison_result
+ compare_slopes_at_intersection_point (const Bez_point_bound& bound1,
+ const Bez_point_bound& bound2)
+ {
+ const Control_points& cp1 = bound1.ctrl;
+ const Control_points& cp2 = bound2.ctrl;
- // An (expensive) check for the precondition that the spans do not overlap.
- // Therefor it is commented out.
- /*
- Vector_2 dir1, dir2, leftMost1, leftMost2, rightMost1, rightMost2;
- bool spansOverlap;
- bool cv1SpanOK = _CrvTanAngularSpan(cv1, dir1, leftMost1, rightMost1);
- bool cv2SpanOK = _CrvTanAngularSpan(cv2, dir2, leftMost2, rightMost2);
- if (cv1SpanOK && cv2SpanOK)
- {
- spansOverlap = _AngularSpansOverlap(leftMost1, rightMost1, leftMost2, rightMost2);
- }
- else
- {
- spansOverlap = true;
- }
- CGAL_precondition(!spansOverlap);
- */
+ // An (expensive) check that the angular spans do not overlap.
+ CGAL_expensive_precondition_code (
+ Vector_2 v_min1;
+ Vector_2 v_max1;
+ const bool span_ok1 = _compute_angular_span (cp1,
+ v_min1, v_max1);
+ Vector_2 v_min2;
+ Vector_2 v_max2;
+ const bool span_ok2 = _compute_angular_span (cp2,
+ v_min2, v_max2);
+ bool spans_overlap = true;
- Vector_2 dir1 = cv1.back() - cv1.front();
- Vector_2 dir2 = cv2.back() - cv2.front();
+ if (span_ok1 && span_ok2)
+ {
+ spans_overlap = _angular_spans_overlap (v_min1, v_max1,
+ v_min2, v_max2);
+ }
+ );
- // RI: use compare_angle_with_x_axis (after doing: if (dir.x()<0) dir=-dir)
- // instead of division.
- NT slope1 = dir1.y()/dir1.x();
- NT slope2 = dir2.y()/dir2.x();
- return CGAL::compare(slope1, slope2);
+ CGAL_expensive_precondition (! spans_overlap);
+
+ // As the angular spans of the control polygons do not overlap, we can
+ // just compare any vector of the hodograph spans.
+ const Vector_2 dir1 = f_construct_vector (cp1.front(), cp1.back());
+ const Vector_2 dir2 = f_construct_vector (cp2.front(), cp2.back());
+
+ // The slopes are given by dir1.y() / dir1.x() and dir2.y()/dir2.x().
+ // However, to avoid potential division by zero, we compare:
+ // dir1.y()*dir2.x() and dir2.y()*dir1.x(), after considering the signs
+ // of the xs.
+ Comparison_result res;
+ bool swap_res = false;
+
+ if (CGAL::sign (dir1.x()) == CGAL::NEGATIVE)
+ swap_res = ! swap_res;
+
+ if (CGAL::sign (dir2.x()) == CGAL::NEGATIVE)
+ swap_res = ! swap_res;
+
+ res = CGAL::compare (dir1.y() * dir2.x(), dir2.y() * dir1.x());
+ if (swap_res)
+ res = CGAL::opposite (res);
+
+ return (res);
}
-
- // Can be made a template function <Control_point_vec, Kernel>
- void cp_bbox(const Control_point_vec& cp, Bez_point_bbox& bez_bbox) const
+ /*!
+ * Construct a bounding box for the given control polygon.
+ * \param cp A sequence of control point (the control polgon).
+ * \param bbox Output: The bounding box.
+ * \pre cp is not empty.
+ */
+ void construct_bbox (const Control_points& cp,
+ Bez_point_bbox& bez_bbox)
{
- CGAL_assertion (! cp.empty());
+ CGAL_precondition (! cp.empty());
- typename Kernel::Compare_x_2 comp_x = kernel.compare_x_2_object();
- typename Kernel::Compare_y_2 comp_y = kernel.compare_y_2_object();
+ // Go over the points, and locate the ones with extremal x and y values.
+ typename Control_points::const_iterator it = cp.begin();
+ typename Control_points::const_iterator it_end = cp.end();
+ typename Control_points::const_iterator min_x = it;
+ typename Control_points::const_iterator max_x = it;
+ typename Control_points::const_iterator min_y = it;
+ typename Control_points::const_iterator max_y = it;
- typename Control_point_vec::const_iterator it = cp.begin();
- typename Control_point_vec::const_iterator min_x = it;
- typename Control_point_vec::const_iterator max_x = it;
- typename Control_point_vec::const_iterator min_y = it;
- typename Control_point_vec::const_iterator max_y = it;
-
- while (it != cp.end())
+ while (it != it_end)
{
- if (comp_x (*it, *min_x) == SMALLER)
+ if (f_compare_x (*it, *min_x) == SMALLER)
min_x = it;
- else if (comp_x (*it, *max_x) == LARGER)
+ else if (f_compare_x (*it, *max_x) == LARGER)
max_x = it;
- if (comp_y (*it, *min_y) == SMALLER)
+ if (f_compare_y (*it, *min_y) == SMALLER)
min_y = it;
- else if (comp_y (*it, *max_y) == LARGER)
+ else if (f_compare_y (*it, *max_y) == LARGER)
max_y = it;
++it;
}
+ // Set the bounding box.
bez_bbox.min_x = min_x->x();
bez_bbox.max_x = max_x->x();
bez_bbox.min_y = min_y->y();
bez_bbox.max_y = max_y->y();
+ return;
}
+ //@}
private:
- /*
- void _construct_bboxes (const Control_point_vec& cv1,
- const Control_point_vec& cv2,
- Bez_point_bbox& bbox1, Bez_point_bbox& bbox2) const
+ /*!
+ * Check whether the given control polygon defines an x-monotone (or
+ * y-monotone) Bezier curve.
+ * \param cp The control points.
+ * \param check_x (true) to check x-monotonicity;
+ * (false) to check y-monotonicity.
+ * \return Is the curve x-monotone (or y-monotone).
+ */
+ bool _is_monotone (const Control_points& cp, bool check_x) const
{
- const double diff_bound = 0.1875; // radians, a bit more than 10 degrees.
- const double _pi = 3.14159265;
+ Comparison_result curr_res;
+ Comparison_result res = EQUAL;
- // Compute the angles theta1 and theta2 the lines that connect the first
- // and last control points in each curve form with the x-axis.
- const double dx1 = CGAL::to_double (cv1.back().x()) -
- CGAL::to_double (cv1.front().x());
- const double dy1 = CGAL::to_double (cv1.back().y()) -
- CGAL::to_double (cv1.front().y());
- const double theta1 = atan2 (dy1, dx1);
+ // Look for the first pair of consecutive points whose x-coordinate
+ // (or y-coordinate) are not equal. Their comparsion result will be
+ // set as the "reference" comparison result.
+ typename Control_points::const_iterator pt_curr = cp.begin();
+ typename Control_points::const_iterator pt_end = cp.end();
+ typename Control_points::const_iterator pt_next = pt_curr;
- const double dx2 = CGAL::to_double (cv2.back().x()) -
- CGAL::to_double (cv2.front().x());
- const double dy2 = CGAL::to_double (cv2.back().y()) -
- CGAL::to_double (cv2.front().y());
- const double theta2 = atan2 (dy2, dx2);
-
- // Check if the two angles are close, and compute half their average.
- const double diff_theta = fabs (theta2 - theta1);
-
- double half_theta;
-
- if (diff_theta < diff_bound)
+ ++pt_next;
+ while (pt_next != pt_end)
{
- // The two angles are close:
- half_theta = (theta1 + theta2) / 4;
- }
- else if (fabs (diff_theta - _pi) < diff_bound)
- {
- // We have theta2 ~= theta1 +/- PI, thus we take:
- if (theta2 > 0)
- half_theta = (theta1 + theta2 - _pi) / 4;
- else
- half_theta = (theta1 + theta2 + _pi) / 4;
- }
- else
- {
- // In this case the angles are not close, so we compute their
- // axis-alligned bounding boxes.
- cp_bbox(cv1, bbox1);
- cp_bbox(cv2, bbox2);
- return;
- }
+ curr_res = check_x ? f_compare_x (*pt_curr, *pt_next) :
+ f_compare_y (*pt_curr, *pt_next);
- // Check if theta/2 is close to 0 or to +/- PI/2.
- if (fabs(half_theta) < 0.0625 || _pi/2 - fabs(half_theta) < 0.0625)
- {
- // In this case, using the axis-alligned bounding box is fine.
- cp_bbox(cv1, bbox1);
- cp_bbox(cv2, bbox2);
- return;
- }
+ pt_curr = pt_next;
+ ++pt_next;
- // Compute a rational tau ~= tan(theta/2) and use it to find a rotation
- // angle close to theta with rational sine and cosine.
- const double tan_half_theta = tan (half_theta);
- const int denom = 1000;
- const int numer = static_cast<int> (denom * tan_half_theta + 0.5);
- const NT tau (numer, denom);
- const NT sqr_tau = tau*tau;
- const NT sin_theta = 2*tau / (1 + sqr_tau);
- const NT cos_theta = (1 - sqr_tau) / (1 + sqr_tau);
-
- // Rotate both control polygons by the minus the approximated angle.
- Control_point_vec rot_cv1;
- Control_point_vec rot_cv2;
- typename Control_point_vec::const_iterator it;
-
- for (it = cv1.begin(); it != cv1.end(); ++it)
- {
- rot_cv1.push_back (Point_2 (cos_theta * it->x() + sin_theta * it->y(),
- -sin_theta * it->x() + cos_theta * it->y()));
- }
-
- for (it = cv2.begin(); it != cv2.end(); ++it)
- {
- rot_cv2.push_back (Point_2 (cos_theta * it->x() + sin_theta * it->y(),
- -sin_theta * it->x() + cos_theta * it->y()));
- }
-
- // Compute the axis-alligned bounding boxes of the rotated curves, which
- // are equivalent of the rotated skewed bounding boxes of the curves.
- cp_bbox (rot_cv1, bbox1);
- cp_bbox (rot_cv2, bbox2);
-
- return;
- }
- */
-
- // Can be made a template function <Control_point_vec, Kernel>
- bool _is_monotone(const Control_point_vec& cp, bool check_x) const
- {
- CGAL_precondition(cp.size() > 1);
- Comparison_result curr_res;
- Comparison_result res = EQUAL;
-
- typename Control_point_vec::const_iterator cp_iter = cp.begin();
- typename Control_point_vec::const_iterator cp_iter_end = cp.end();
- typename Control_point_vec::const_iterator cp_iter_next = cp_iter;
- ++cp_iter_next;
-
- while (cp_iter_next != cp_iter_end)
- {
- curr_res = check_x ? compare_x(*cp_iter, *cp_iter_next) :
- compare_y(*cp_iter, *cp_iter_next);
if (curr_res != EQUAL)
{
res = curr_res;
- ++cp_iter;
- ++cp_iter_next;
break;
}
- ++cp_iter;
- ++cp_iter_next;
}
- if (cp_iter_next == cp_iter_end)
- // iso-parallel segment.
+
+ if (pt_next == pt_end || res == EQUAL)
+ {
+ // If we reached here, we have an iso-parallel segment.
// In our context, we consider a horizontal segment as non-y-monotone
- // but a vertical segment is considered x-monotone (is this right?).
- {
+ // but a vertical segment is considered x-monotone.
if (check_x)
- return true; // Vertical segment when checking x-monoticity
- return false; // Horizontal segment when checking y-monoticity
+ return (true); // Vertical segment when checking x-monoticity
+
+ return (false); // Horizontal segment when checking y-monoticity
}
- while (cp_iter_next != cp_iter_end)
+ // Go over the rest of the control polygons, and make sure that the
+ // comparison result between each pair of consecutive points is the
+ // same as the reference result.
+ while (pt_next != pt_end)
{
- curr_res = check_x ? compare_x(*cp_iter, *cp_iter_next) :
- compare_y(*cp_iter, *cp_iter_next);
- if (curr_res == EQUAL)
- {
- ++cp_iter;
- ++cp_iter_next;
- continue;
- }
- if (res != curr_res)
- break;
+ curr_res = check_x ? f_compare_x (*pt_curr, *pt_next) :
+ f_compare_y (*pt_curr, *pt_next);
- ++cp_iter;
- ++cp_iter_next;
+ if (curr_res != EQUAL && res != curr_res)
+ return (false);
+
+ pt_curr = pt_next;
+ ++pt_next;
}
- if (cp_iter_next != cp_iter_end)
- return false;
- return true;
+
+ // If we reached here, the control polygon is x-monotone (or y-monotone).
+ return (true);
}
- inline bool _is_x_monotone(const Control_point_vec& cp) const
+ /*!
+ * Check whether the control polygon defines an x-monotone Bezier curve.
+ * \param cp The control points.
+ * \return Is the curve x-monotone.
+ */
+ inline bool _is_x_monotone (const Control_points& cp) const
{
- return _is_monotone(cp, true);
+ return (_is_monotone (cp, true));
}
- inline bool _is_y_monotone(const Control_point_vec& cp) const
+ /*!
+ * Check whether the control polygon defines a y-monotone Bezier curve.
+ * \param cp The control points.
+ * \return Is the curve y-monotone.
+ */
+ inline bool _is_y_monotone (const Control_points& cp) const
{
- return _is_monotone(cp, false);
+ return (_is_monotone (cp, false));
}
-
- // Can be made a template function <Control_point_vec, Kernel>
- // Returns true if the span is less than 90 degrees (we don't want to work with larger spans).
- bool _CrvTanAngularSpan(const Control_point_vec& cp,
- Vector_2& dir,
- Vector_2& leftMost,
- Vector_2& rightMost)
+ /*!
+ * Compute the angular span of a Bezier curve, given by its control polygon.
+ * \param cp The control points.
+ * \param v_min Output: The minimal direction of the angular span.
+ * \param v_max Output: The maximal direction of the angular span.
+ * \return Whether the span is less than 90 degrees
+ * (as we don't want to work with larger spans).
+ */
+ bool _compute_angular_span (const Control_points& cp,
+ Vector_2& v_min,
+ Vector_2& v_max)
{
- Vector_2 v;
- const Point_2& O(ORIGIN);
+ // Initialize the output directions.
+ Vector_2 dir = f_construct_vector (cp.front(), cp.back());
+ Vector_2 v;
- dir = cp.back() - cp.front();
+ v_min = v_max = dir;
- // Initialize LeftMost and RightMost
- leftMost = dir;
- rightMost = dir;
+ // Go over the control points and examine the vectors defined by pairs
+ // of consecutive control points.
+ typename Control_points::const_iterator pt_curr = cp.begin();
+ typename Control_points::const_iterator pt_end = cp.end();
+ typename Control_points::const_iterator pt_next = pt_curr;
- typename Control_point_vec::const_iterator iterCP0 = cp.begin();
- typename Control_point_vec::const_iterator iterCP1 = iterCP0;
- ++iterCP1;
- for (; iterCP1!=cp.end(); ++iterCP1, ++iterCP0)
+ ++pt_next;
+ while (pt_next != pt_end)
{
- v = (*iterCP1) - (*iterCP0);
+ v = f_construct_vector (*pt_curr, *pt_next);
- // Test for 90 degree angle, we can loosen it to test for 180 degree between
- // left and right but we don't really want to work with larger spans.
- // TODO - try get rid of these ORIGINs
- if (angle(ORIGIN+v, O, ORIGIN+dir) != ACUTE)
- return false;
+ // If the current vector forms a right-turn with dir, it is a candidate
+ // for the minimal vector, otherwise it is a candidate for the maximal
+ // vector in the span.
+ bool updated_span = false;
- if (left_turn(ORIGIN+v, O, ORIGIN+dir))
+ if (f_orientation (f_construct_point (ORIGIN, v),
+ ORIGIN,
+ f_construct_point (ORIGIN, dir)) == RIGHT_TURN)
{
- if (left_turn(ORIGIN+v, O, ORIGIN+leftMost))
+ // Check if we need to update the minimal vector in the span.
+ if (f_orientation (f_construct_point (ORIGIN, v),
+ ORIGIN,
+ f_construct_point (ORIGIN, v_min)) == RIGHT_TURN)
{
- leftMost = v;
- }
- }
- else
+ v_min = v;
+ updated_span = true;
+ }
+ }
+ else
{
- if (right_turn(ORIGIN+v, O, ORIGIN+rightMost))
+ // Check if we need to update the maximal vector in the span.
+ if (f_orientation (f_construct_point (ORIGIN, v),
+ ORIGIN,
+ f_construct_point (ORIGIN, v_max)) == LEFT_TURN)
{
- rightMost = v;
- }
- }
- } //for
+ v_max = v;
+ updated_span = true;
+ }
+ }
- return true;
- }
+ // Before proceeding, check if the current span is larger than
+ // 180 degrees. We do that by checking that (v_min, 0, v_max)
+ // is still a right-turn, and did not become a left-turn.
+ if (updated_span &&
+ f_orientation (f_construct_point (ORIGIN, v_min),
+ ORIGIN,
+ f_construct_point (ORIGIN, v_max)) != RIGHT_TURN)
+ {
+ return (false);
+ }
- // TODO - there is already a kernel function that does this ccw?
- // Can be made a template function <Kernel>
- inline bool _VecIsBetween(const Vector_2& vec,
- const Vector_2& leftVec,
- const Vector_2& rightVec)
- {
- const Point_2& O(ORIGIN);
- return (right_turn(ORIGIN+rightVec, O, ORIGIN+vec) && left_turn(ORIGIN+leftVec, O, ORIGIN+vec));
+ // Move to the next pair of points.
+ pt_curr = pt_next;
+ ++pt_next;
+ }
+
+ // If we reached here, the angular span is less than 180 degrees.
+ return (true);
}
- // Can be made a template function <Kernel>
- bool _AngularSpansOverlap(
- const Vector_2& leftMost1,
- const Vector_2& rightMost1,
- const Vector_2& leftMost2,
- const Vector_2& rightMost2
- )
+ /*!
+ * Check whether to angular spans overlap.
+ * \param v_min1 The minimal vector in the first angular span.
+ * \param v_max1 The maximal vector in the first angular span.
+ * \param v_min2 The minimal vector in the second angular span.
+ * \param v_max2 The maximal vector in the second angular span.
+ * \return Do the two spans overlap.
+ */
+ bool _angular_spans_overlap (const Vector_2& v_min1,
+ const Vector_2& v_max1,
+ const Vector_2& v_min2,
+ const Vector_2& v_max2)
{
- Vector_2 negVec;
+ const Direction_2 dir_l1 = f_construct_direction (v_min1);
+ const Direction_2 dir_r1 = f_construct_direction (v_max1);
+ const Direction_2 dir_l2 = f_construct_direction (v_min2);
+ const Direction_2 dir_r2 = f_construct_direction (v_max2);
- // Check if any of span2 vectors is between span1.
- if (_VecIsBetween(leftMost2, leftMost1, rightMost1))
- return true;
+ // Check whether any of the vectors of the second span (or their
+ // opposite vectors) is between the vectors of the first span.
+ if (! f_equal (v_min1, v_max1) &&
+ (f_ccw_in_between (dir_l2, dir_l1, dir_r1) ||
+ f_ccw_in_between (f_opposite_direction (dir_l2), dir_l1, dir_r1) ||
+ f_ccw_in_between (dir_r2, dir_l1, dir_r1) ||
+ f_ccw_in_between (f_opposite_direction (dir_r2), dir_l1, dir_r1)))
+ {
+ return (true);
+ }
- negVec = -leftMost2;
- if (_VecIsBetween(negVec, leftMost1, rightMost1))
- return true;
+ // Check whether the left vector of the first span (or its opposite
+ // vector) is between the vectors of the second span. Note that at this
+ // point, the only possibility is that the first span is totally contained
+ // within the second, so we do not need to check both vectors.
+ if (! f_equal (v_min2, v_max2) &&
+ (f_ccw_in_between (dir_l1, dir_l2, dir_r2) ||
+ f_ccw_in_between (f_opposite_direction (dir_l1), dir_l2, dir_r2)))
+ {
+ return (true);
+ }
- if (_VecIsBetween(rightMost2, leftMost1, rightMost1))
- return true;
+ // If we reached here, the two angular spans do not overlap.
+ return (false);
+ }
- negVec = -rightMost2;
- if (_VecIsBetween(negVec, leftMost1, rightMost1))
- return true;
+ /*!
+ * Construct a skewed bounding box for the control polygon.
+ * The skewed bounding box is represented as a pair of lines, both parallel
+ * to the line connecting the first and last points of the control polygon,
+ * that bound the curve from above and below.
+ * \param cp The control points.
+ * \param l_min Output: The line with minimal signed distance.
+ * \param l_max Output: The line with maximal signed distance.
+ */
+ void _skewed_bbox (const Control_points& cp,
+ Line_2& l_min, Line_2& l_max)
+ {
+ // Construct the line that connects the first and the last control points.
+ const Line_2 l = f_construct_line (cp.front(), cp.back());
+ const Vector_2 v = f_construct_vector (cp.front(), cp.back());
+ // Select the points with minimum and maximum distance from the line l,
+ // and construct two lines parallel to l that pass through these points.
+ typename Control_points::const_iterator pt_curr = cp.begin();
+ typename Control_points::const_iterator pt_end = cp.end();
- // Note: If we got here, the only possibility for intersection
- // is that span1 is totally between span2
+ typename Control_points::const_iterator pt_min = pt_curr;
+ typename Control_points::const_iterator pt_max = pt_curr;
- // Check if span1 is between span2
- if (_VecIsBetween(leftMost1, leftMost2, rightMost2))
- return true;
+ ++pt_curr;
+ while (pt_curr != pt_end)
+ {
+ if (f_less_signed_distance (l, *pt_curr, *pt_min))
+ {
+ pt_min = pt_curr;
+ }
+ else if (f_less_signed_distance (l, *pt_max, *pt_curr))
+ {
+ pt_max = pt_curr;
+ }
- negVec = -leftMost1;
- if (_VecIsBetween(negVec, leftMost2, rightMost2))
- return true;
+ ++pt_curr;
+ }
-
- // The rest of the symmetric conditions below are redundant (see note above)
- return false;
-
- /*
- if (_VecIsBetween(rightMost1, leftMost2, rightMost2))
- return true;
-
- negVec = -rightMost1;
- if (_VecIsBetween(negVec, leftMost2, rightMost2))
- return true;
-
- return false;
- */
+ // Construct the output lines.
+ l_min = f_construct_line (*pt_min, v);
+ l_max = f_construct_line (*pt_max, v);
+ return;
}
- // Constructing the curve skewed bbox - actually returning
- // two parallel lines (parallel to the first-last line) that bound
- // the curve from above and below.
- // Using less_signed_distance_to_line_2 enables not to construct
- // new points (only take maximum and minimum signed dist points
- // and add the first-last vector).
- // Can be made a template function <Control_point_vec, Kernel>
- void _cvSkewBbox(const Control_point_vec& cp,
- Line_2& la,
- Line_2& lb)
+ /*!
+ * Check whether the endpoints of the two curves coincide.
+ * If so, it adds the Intersection_point to the list if the endpoint hasn't
+ * already been inserted there.
+ * Note that this function is not able to detect an endpoint that lies in
+ * the interior of the other curve. Such cases will not be approximated
+ * and we will have to resort to exact computation.
+ * \param cp1 The control points of the first curve.
+ * \param t_min1 The lower bound of the parameter range of the first curve.
+ * \param t_max1 The upper bound of the parameter range of the first curve.
+ * \param cp2 The control points of the second curve.
+ * \param t_min2 The lower bound of the parameter range of the second curve.
+ * \param t_max2 The upper bound of the parameter range of the second curve.
+ * \param iept_params Input/Output: The set of parameter values (for cp1)
+ * of known intersections at the endpoints.
+ * \param ipts Input/Output: The computed intersection points.
+ * \return Whether a pair of curve endpoints coincide.
+ */
+ bool _endpoints_coincide (const Control_points& cp1,
+ const NT& t_min1, const NT& t_max1,
+ const Control_points& cp2,
+ const NT& t_min2, const NT& t_max2,
+ std::set<NT>& iept_params,
+ std::list<Intersection_point>& ipts) const
{
- // Functor for comparing two points signed distance from a given line (used for min/max_elem).
- // Binder taken from <CGAL/functional.h>, if this doesn't compile, construct your own functor.
- typedef CGALi::Binder<typename Kernel::Less_signed_distance_to_line_2, Arity_tag< 3 >, Line_2, 1 > Less_signed_dist_to_line;
- // Explanation: Take the 3-ary kernel functor Less_signed_distance_to_line_2 and bind its
- // first parameter, the line, to be our given line l.
+ // Get the first and last control points of each curve.
+ const Point_2& s1 = cp1.front();
+ const Point_2& t1 = cp1.back();
+ const Point_2& s2 = cp2.front();
+ const Point_2& t2 = cp2.back();
- Line_2 l(cp.front(), cp.back());
- Vector_2 v(cp.front(), cp.back());
- typename Kernel::Less_signed_distance_to_line_2 lsdtl2 = kernel.less_signed_distance_to_line_2_object();
+ // Check whether any pair of these endpoints conincide.
+ NT x, y; // Coordinate of a common endpoint.
+ NT t_val1, t_val2; // Its respective parameters.
- typename Control_point_vec::const_iterator aux;
- aux = std::min_element(cp.begin(), cp.end(), Less_signed_dist_to_line(lsdtl2, l));
- la = Line_2(*aux, v);
-
- aux = std::max_element(cp.begin(), cp.end(), Less_signed_dist_to_line(lsdtl2, l));
- lb = Line_2(*aux, v);
- }
-
-
- // Returns true if endpoints coincide, and adds the Bound_pair
- // to the list if the endpoint hasn't already been inserted there.
- // Precondition !spansOverlap
- bool _endpoints_coincide(const Control_point_vec& cv1, const NT& l1, const NT& r1,
- const Control_point_vec& cv2, const NT& l2, const NT& r2,
- std::set<NT>& endpoint_intersections,
- Bound_pair_lst& intersection_pairs) const
- {
- const Point_2& s1 = cv1.front();
- const Point_2& t1 = cv1.back();
- const Point_2& s2 = cv2.front();
- const Point_2& t2 = cv2.back();
-
- // What happens when an endpoint of one is an inner point (e.g., 1/3) of another?
- // We will not be able to discover this endcase and therefore get to can_refine()==false.
- int endpoint_intersection_num = 0; // For sanity check and return value.
-
- NT x, y, t_param1, t_param2;
- if (s1==s2) // Using operator== of kernel Point_2
+ if (f_equal (s1, s2))
{
- ++endpoint_intersection_num;
x = s1.x();
y = s1.y();
- t_param1 = l1;
- t_param2 = l2;
+ t_val1 = t_min1;
+ t_val2 = t_min2;
}
- if (s1==t2)
+ else if (f_equal (s1, t2))
{
- ++endpoint_intersection_num;
x = s1.x();
y = s1.y();
- t_param1 = l1;
- t_param2 = r2;
+ t_val1 = t_min1;
+ t_val2 = t_max2;
}
- if (t1==s2)
+ else if (f_equal (t1, s2))
{
- ++endpoint_intersection_num;
x = t1.x();
y = t1.y();
- t_param1 = r1;
- t_param2 = l2;
+ t_val1 = t_max1;
+ t_val2 = t_min2;
}
- if (t1==t2)
+ else if (f_equal (t1, t2))
{
- ++endpoint_intersection_num;
x = t1.x();
y = t1.y();
- t_param1 = r1;
- t_param2 = r2;
+ t_val1 = t_max1;
+ t_val2 = t_max2;
}
+ else
+ {
+ // No common endpoint found:
+ return (false);
+ }
- CGAL_assertion(endpoint_intersection_num <= 1);
- if (endpoint_intersection_num == 0)
- return false;
+ // Try inserting the parameter value t1 into the set of parameters of
+ // already discovered intersecting endpoints.
+ std::pair<typename std::set<NT>::iterator,
+ bool> res = iept_params.insert (t_val1);
- // Construct a degenerate rational Bound_pair
- if (endpoint_intersections.find(t_param1) == endpoint_intersections.end())
+ if (res.second)
{
- // Found a new intersection point that wasn't inserted yet.
- endpoint_intersections.insert(t_param1);
- intersection_pairs.push_back(Bound_pair(
- Bez_point_bound(cv1, t_param1, t_param1, Bez_point_bound::RATIONAL_PT, true),
- Bez_point_bound(cv2, t_param2, t_param2, Bez_point_bound::RATIONAL_PT, true),
- Bez_point_bbox(x,x,y,y)
- )
- );
+ // In case the insertion has succeeded, report on a new rational
+ // intersection point.
+ Bez_point_bound bound1 (Bez_point_bound::RATIONAL_PT,
+ cp1, t_val1, t_val1, true);
+ Bez_point_bound bound2 (Bez_point_bound::RATIONAL_PT,
+ cp2, t_val2, t_val2, true);
+ Bez_point_bbox bbox (x, x, y, y);
+
+ ipts.push_back (Intersection_point (bound1, bound2, bbox));
}
- return true;
+
+ return (true);
}
-
- // Auxilary recursive function for intersecting two curves.
- void _CrvCrvInterAux
- (const Control_point_vec& cv1, const NT& l1, const NT& r1,
- const Control_point_vec& cv2, const NT& l2, const NT& r2,
- bool should_check_span,
- std::set<NT>& endpoint_intersections,
- Bound_pair_lst& intersection_pairs)
+ /*!
+ * An auxilary recursive function for computing the approximated
+ * intersection points between two Bezier curves.
+ * \param cp1 The control points of the first curve.
+ * \param t_min1 The lower bound of the parameter range of the first curve.
+ * \param t_max1 The upper bound of the parameter range of the first curve.
+ * \param cp2 The control points of the second curve.
+ * \param t_min2 The lower bound of the parameter range of the second curve.
+ * \param t_max2 The upper bound of the parameter range of the second curve.
+ * \param check_span Should we check the angular span of the curves.
+ * \param iept_params Input/Output: The set of parameter values (for cp1)
+ * of known intersections at the endpoints.
+ * \param ipts Input/Output: The computed intersection points.
+ */
+ void _compute_intersection_points (const Control_points& cp1,
+ const NT& t_min1, const NT& t_max1,
+ const Control_points& cp2,
+ const NT& t_min2, const NT& t_max2,
+ bool check_span,
+ std::set<NT>& iept_params,
+ std::list<Intersection_point>& ipts)
{
- //iddo debug
- //static unsigned int dbg = 0;
- //std::cout << ++dbg << ", " << recursion_level << "; ";
+ // Check if we got to subdivision termination criteria.
+ const bool can_refine1 = can_refine (cp1, t_min1, t_max1);
+ const bool can_refine2 = can_refine (cp2, t_min2, t_max2);;
- // Check if we got to subdivision termination criteria.
- bool can_refine1 = can_refine(cv1, l1, r1);
- bool can_refine2 = can_refine(cv2, l2, r2);;
- if (!can_refine1 || !can_refine2)
+ if (! can_refine1 || ! can_refine2)
{
- //Failed in the bounded approach - stop the subdivision.
- //Construct output with can_refine=false... We don't really need a
- // meaningful bbox since we can't refine, so we take cv1 bbox (which
- // certainly contains the point).
- Bez_point_bbox bbox1;
+ // It is not possible to further refine the approximation, so we stop
+ // the recursive subdivision, and construct an output intersection point
+ // with can_refine = false. Note that we do not need a tight bounding
+ // box here, since we cannot refine it anyway, so we take cp1's bbox
+ // (which certainly contains the intersection point).
+ Bez_point_bound bound1 (Bez_point_bound::INTERSECTION_PT,
+ cp1, t_min1, t_max1,
+ false); // Cannot refine further.
+ Bez_point_bound bound2 (Bez_point_bound::INTERSECTION_PT,
+ cp2, t_min2, t_max2,
+ false); // Cannot refine further.
+ Bez_point_bbox ipt_bbox;
- cp_bbox(cv1, bbox1);
- intersection_pairs.push_back
- (Bound_pair (Bez_point_bound (cv1, l1, r1,
- Bez_point_bound::INTERSECTION_PT, false),
- Bez_point_bound (cv2, l2, r2,
- Bez_point_bound::INTERSECTION_PT, false),
- bbox1));
+ construct_bbox (cp1, ipt_bbox);
+ ipts.push_back (Intersection_point (bound1, bound2, ipt_bbox));
+ m_active_nodes--;
return;
}
- Bez_point_bbox bbox1, bbox2;
+ // Consturct bounding boxes for the two curves and check whether they
+ // overlap.
+ Bez_point_bbox bbox1;
+ Bez_point_bbox bbox2;
- //_construct_bboxes (cv1, cv2,
- // bbox1, bbox2);
- cp_bbox(cv1, bbox1);
- cp_bbox(cv2, bbox2);
+ construct_bbox (cp1, bbox1);
+ construct_bbox (cp2, bbox2);
- // The bounding boxes do not overlap - return.
- if (!bbox1.Overlaps(bbox2))
+ if (! bbox1.overlaps (bbox2))
{
- //iddo debug
- //std::cout << "boxes do not overlap - curve do not intersect" << std::endl;
+ // The bounding boxes do not overlap, so the two input curves do not
+ // intersect.
+ m_active_nodes--;
return;
}
- else // Debug print
- {
- //std::cout << "boxes overlap - continuing..." << std::endl;
- }
- Vector_2 dir1, dir2, leftMost1, rightMost1, leftMost2, rightMost2;
+ // Check the angular spans, if necessary.
+ bool spans_overlap = false;
- bool cv1SpanOK = true;
- bool cv2SpanOK = true;
- bool spansOverlap = false;
+ if (check_span)
+ {
+ // Compute the angular spans of the two curves.
+ Vector_2 v_min1, v_max1;
+ const bool span_ok1 = _compute_angular_span (cp1,
+ v_min1, v_max1);
+ Vector_2 v_min2, v_max2;
+ const bool span_ok2 = _compute_angular_span (cp2,
+ v_min2, v_max2);
- if (should_check_span)
- {
- cv1SpanOK = _CrvTanAngularSpan(cv1, dir1, leftMost1, rightMost1);
- cv2SpanOK = _CrvTanAngularSpan(cv2, dir2, leftMost2, rightMost2);
- if (cv1SpanOK && cv2SpanOK)
+ if (span_ok1 && span_ok2)
{
- spansOverlap = _AngularSpansOverlap(leftMost1, rightMost1, leftMost2, rightMost2);
+ spans_overlap = _angular_spans_overlap (v_min1, v_max1,
+ v_min2, v_max2);
}
else
{
- spansOverlap = true;
+ // One of the spans is greater than 180 degrees, so we do not have to
+ // check for overlaps (we know that it must overlap the other span).
+ spans_overlap = true;
}
}
- if (!spansOverlap)
+ if (! spans_overlap)
{
- // Debug print.
- //std::cout << "Spans do NOT Overlap." << std::endl;
+ // In case the spans do not overlap, we potentially have a single
+ // intersection point.
- // Debug print
- /*
- std::cout << "CV1:\n";
- Control_point_vec::const_iterator debug_iter = cv1.begin();
- while (debug_iter != cv1.end())
- {
- std::cout << "[" << *debug_iter << "]";
- ++debug_iter;
- }
- std::cout << std::endl;
- std::cout << "CV2:\n";
- debug_iter = cv2.begin();
- while (debug_iter != cv2.end())
- {
- std::cout << "[" << *debug_iter << "]";
- ++debug_iter;
- }
- std::cout << std::endl;
- */
-
// Checking for endpoint intersections.
- bool endpoint_coincidence = _endpoints_coincide(cv1, l1, r1,
- cv2, l2, r2,
- endpoint_intersections,
- intersection_pairs);
-
- if (endpoint_coincidence)
+ if (_endpoints_coincide (cp1, t_min1, t_max1,
+ cp2, t_min2, t_max2,
+ iept_params,
+ ipts))
{
- // Debug print
- //std::cout << "endpoint coincidence intersection." << std::endl;
- //std::cout << "curves intersect! intersectionIntervals.size() == " << intersection_pairs.size() << std::endl;
-
- //Since spans do not overlap and there is an endpoint coincidence
- //there is no need for further checking (the endpoint is the intersection).
+ // As we have located the single intersection point (a common endpoint
+ // in this case), we can stop here.
+ m_active_nodes--;
return;
}
- const Point_2& s1 = cv1.front();
- const Point_2& t1 = cv1.back();
- const Point_2& s2 = cv2.front();
- const Point_2& t2 = cv2.back();
+ // Construct the skewed bounding boxes for the two curves and check
+ // whether the endpoints of the first curve lie on opposite side of
+ // the skewed bounding box of the second curve, and vice versa.
+ Line_2 skew1a, skew1b;
+ Line_2 skew2a, skew2b;
+ const Point_2& s1 = cp1.front();
+ const Point_2& t1 = cp1.back();
+ const Point_2& s2 = cp2.front();
+ const Point_2& t2 = cp2.back();
- Line_2 skew1a, skew1b;
- _cvSkewBbox(cv1, skew1a, skew1b);
+ _skewed_bbox(cp1, skew1a, skew1b);
+ _skewed_bbox(cp2, skew2a, skew2b);
- Line_2 skew2a, skew2b;
- _cvSkewBbox(cv2, skew2a, skew2b);
+ const Oriented_side or_2a_s1 = f_oriented_side (skew2a, s1);
+ const Oriented_side or_2a_t1 = f_oriented_side (skew2a, t1);
- // Check for opposite orientations between:
- // (skew2a, s1)(skew2a, t1)
- // (Skew2b, s1)(Skew2b, t1)
- // and vice versa
- //RI: TODO not use < 0 (count on numeric value) on Oriented_side
- Oriented_side or_2a_s1 = skew2a.oriented_side(s1);
- Oriented_side or_2a_t1 = skew2a.oriented_side(t1);
+ const Oriented_side or_2b_s1 = f_oriented_side (skew2b, s1);
+ const Oriented_side or_2b_t1 = f_oriented_side (skew2b, t1);
- Oriented_side or_2b_s1 = skew2b.oriented_side(s1);
- Oriented_side or_2b_t1 = skew2b.oriented_side(t1);
+ const bool s1_t1_are_opposite =
+ ((or_2a_s1 == CGAL::opposite (or_2a_t1)) &&
+ (or_2b_s1 == CGAL:: opposite (or_2b_t1)));
- bool s1_t1_are_opposite = ((or_2a_s1*or_2a_t1 < 0) &&
- (or_2b_s1*or_2b_t1 < 0));
+ const Oriented_side or_1a_s2 = f_oriented_side (skew1a, s2);
+ const Oriented_side or_1a_t2 = f_oriented_side (skew1a, t2);
+ const Oriented_side or_1b_s2 = f_oriented_side (skew1b, s2);
+ const Oriented_side or_1b_t2 = f_oriented_side (skew1b, t2);
- Oriented_side or_1a_s2 = skew1a.oriented_side(s2);
- Oriented_side or_1a_t2 = skew1a.oriented_side(t2);
+ const bool s2_t2_are_opposite =
+ ((or_1a_s2 == CGAL::opposite (or_1a_t2)) &&
+ (or_1b_s2 == CGAL::opposite (or_1b_t2)));
- Oriented_side or_1b_s2 = skew1b.oriented_side(s2);
- Oriented_side or_1b_t2 = skew1b.oriented_side(t2);
-
- bool s2_t2_are_opposite = ((or_1a_s2*or_1a_t2 < 0) &&
- (or_1b_s2*or_1b_t2 < 0));
-
if (s1_t1_are_opposite && s2_t2_are_opposite)
{
- // Construct the bbox from intersection of skew lines.
- Bez_point_bbox inter_bbox;
- Control_point_vec aux_vec;
- Object res;
- Point_2 p;
+ // Construct a finer bounding box for the intersection point from
+ // the intersection of the two skewed bounding boxes.
+ Bez_point_bbox ipt_bbox;
+ Control_points aux_vec;
+ CGAL::Object res;
+ Point_2 p;
- res = intersection (skew1a, skew2a);
- if (!assign(p, res))
+ res = f_intersect (skew1a, skew2a);
+ if (! assign (p, res))
{
CGAL_assertion(false);
}
aux_vec.push_back(p);
- res = intersection (skew1a, skew2b);
- if (!assign(p, res))
+ res = f_intersect (skew1a, skew2b);
+ if (! assign(p, res))
{
CGAL_assertion(false);
}
aux_vec.push_back(p);
- res = intersection (skew1b, skew2a);
- if (!assign(p, res))
+ res = f_intersect (skew1b, skew2a);
+ if (! assign(p, res))
{
CGAL_assertion(false);
}
aux_vec.push_back(p);
- res = intersection (skew1b, skew2b);
+ res = f_intersect (skew1b, skew2b);
if (!assign(p, res))
{
CGAL_assertion(false);
}
aux_vec.push_back(p);
- cp_bbox(aux_vec, inter_bbox);
+ construct_bbox (aux_vec, ipt_bbox);
- intersection_pairs.push_back(Bound_pair(
- Bez_point_bound(cv1, l1, r1, Bez_point_bound::INTERSECTION_PT, true),
- Bez_point_bound(cv2, l2, r2, Bez_point_bound::INTERSECTION_PT, true),
- inter_bbox
- )
- );
+ // Report on the intersection point we have managed to isolate.
+ Bez_point_bound bound1 (Bez_point_bound::INTERSECTION_PT,
+ cp1, t_min1, t_max1,
+ true); // We can further refine it.
+ Bez_point_bound bound2 (Bez_point_bound::INTERSECTION_PT,
+ cp2, t_min2, t_max2,
+ true); // We can further refine it.
- // Debug print
- //std::cout << "curves intersect! intersectionIntervals.size() == " << intersection_pairs.size() << std::endl;
- // gets(aux);
+ ipts.push_back (Intersection_point (bound1, bound2, ipt_bbox));
- return;
+ m_active_nodes--;
+ return;
}
}
- // Subdivide the two curves and recurse.
- Control_point_vec cv1a, cv1b;
+ // Apply de Casteljau's algorithm and bisect both bounding polylines.
+ Control_points cp1a, cp1b;
+ const NT t_mid1 = (t_min1 + t_max1) / 2;
- bisect_control_polygon_2 (cv1.begin(), cv1.end(),
- std::back_inserter(cv1a),
- std::front_inserter(cv1b));
- //DeCasteljau(cv1, 0.5, cv1a, cv1b);
- NT t_mid1 = NT(0.5) * (l1 + r1);
+ bisect_control_polygon_2 (cp1.begin(), cp1.end(),
+ std::back_inserter(cp1a),
+ std::front_inserter(cp1b));
- Control_point_vec cv2a, cv2b;
-
- bisect_control_polygon_2 (cv2.begin(), cv2.end(),
- std::back_inserter(cv2a),
- std::front_inserter(cv2b));
- //DeCasteljau(cv2, 0.5, cv2a, cv2b);
- NT t_mid2 = NT(0.5) * (l2 + r2);
+ Control_points cp2a, cp2b;
+ const NT t_mid2 = (t_min2 + t_max2) / 2;
- // Recursion:
- _CrvCrvInterAux (cv1a, l1, t_mid1, cv2a, l2, t_mid2,
- spansOverlap, endpoint_intersections,
- intersection_pairs);
- _CrvCrvInterAux (cv1a, l1, t_mid1, cv2b, t_mid2, r2,
- spansOverlap, endpoint_intersections,
- intersection_pairs);
- _CrvCrvInterAux (cv1b, t_mid1, r1, cv2a, l2, t_mid2,
- spansOverlap, endpoint_intersections,
- intersection_pairs);
- _CrvCrvInterAux (cv1b, t_mid1, r1, cv2b, t_mid2, r2,
- spansOverlap, endpoint_intersections,
- intersection_pairs);
+ bisect_control_polygon_2 (cp2.begin(), cp2.end(),
+ std::back_inserter(cp2a),
+ std::front_inserter(cp2b));
+
+ // Recursively compute the intersection points on all pairs of bisected
+ // curves.
+ m_active_nodes += 4;
+ if (m_active_nodes > 4 * cp1.size() * cp2.size())
+ {
+ // It is not possible to further refine the approximation, as the
+ // number of active nodes in the recursive search is too high. We stop
+ // the recursive subdivision, and construct an output intersection point
+ // with can_refine = false. Note that we do not need a tight bounding
+ // box here, since we cannot refine it anyway, so we take cp1's bbox
+ // (which certainly contains the intersection point).
+ Bez_point_bound bound1 (Bez_point_bound::INTERSECTION_PT,
+ cp1, t_min1, t_max1,
+ false); // Cannot refine further.
+ Bez_point_bound bound2 (Bez_point_bound::INTERSECTION_PT,
+ cp2, t_min2, t_max2,
+ false); // Cannot refine further.
+ Bez_point_bbox ipt_bbox;
+
+ construct_bbox (cp1, ipt_bbox);
+ ipts.push_back (Intersection_point (bound1, bound2, ipt_bbox));
+
+ return;
+ }
+
+ _compute_intersection_points (cp1a, t_min1, t_mid1,
+ cp2a, t_min2, t_mid2,
+ spans_overlap,
+ iept_params,
+ ipts);
+
+ _compute_intersection_points (cp1a, t_min1, t_mid1,
+ cp2b, t_mid2, t_max2,
+ spans_overlap,
+ iept_params,
+ ipts);
+
+ _compute_intersection_points (cp1b, t_mid1, t_max1,
+ cp2a, t_min2, t_mid2,
+ spans_overlap,
+ iept_params,
+ ipts);
+
+ _compute_intersection_points (cp1b, t_mid1, t_max1,
+ cp2b, t_mid2, t_max2,
+ spans_overlap,
+ iept_params,
+ ipts);
+
+ m_active_nodes--;
+ return;
}
-};
+ /*!
+ * An auxilary recursive function for computing the approximated vertical
+ * tangency points of a Bezier curves.
+ * \param cp The control points of the curve.
+ * \param t_min The lower bound of the parameter range of the curve.
+ * \param t_max The upper bound of the parameter range of the curve.
+ * \param vpts Input/Output: The computed vertical tangency points.
+ */
+ void _compute_vertical_tangency_points
+ (const Control_points& cp,
+ const NT& t_min, const NT& t_max,
+ std::list<Vertical_tangency_point>& vpts)
+ {
+ // \todo Handle the special case of degree two curves.
-//TODO - These classes can both go outside to another file.
-//if we decide to take some of the functionality to a template algo file.
-template <class _NT>
-struct _Bez_point_bound {
- enum Point_type {RATIONAL_PT, VERTICAL_TANGENCY_PT, INTERSECTION_PT, UNDEFINED};
+ // Check if we got to subdivision termination criteria.
+ if (! can_refine (cp, t_min, t_max))
+ {
+ // It is not possible to further refine the approximation, so we stop
+ // the recursive subdivision, and construct an output vertical tangency
+ // point with can_refine = false. Note that we do not need a tight
+ // bounding box here, since we cannot refine it anyway, so we take cp's
+ // bbox (which certainly contains the intersection point).
+ Bez_point_bound bound (Bez_point_bound::VERTICAL_TANGENCY_PT,
+ cp, t_min, t_max,
+ false);
+ Bez_point_bbox vpt_bbox;
- typedef _NT NT;
- typedef typename Cartesian<NT>::Point_2 Point_2;
- typedef std::deque<Point_2> Control_point_vec;
+ construct_bbox (cp, vpt_bbox);
+ vpts.push_back (Vertical_tangency_point (bound,
+ vpt_bbox));
- Control_point_vec bounding_polyline;
+ return;
+ }
- NT t_min; // iddo: in the future, maybe a better rep for these numbers
- NT t_max; // possibly a mantissa+exp representation.
- Point_type point_type;
- bool can_refine;
+ // If the control polygon is x-monotone, the curve does not contain any
+ // vertical tangency points.
+ if (_is_x_monotone (cp))
+ {
+ return;
+ }
- _Bez_point_bound() : bounding_polyline(),
- t_min(), t_max(), point_type(UNDEFINED), can_refine(false)
- {}
+ // Check whether the control polygon is y-monotone.
+ if (! _is_y_monotone (cp))
+ {
+ // Use de Casteljau's algorithm and subdivide the control polygon into
+ // two, in order to obtain y-monotone subcurves.
+ Control_points cp_a, cp_b;
+ const NT t_mid = (t_min + t_max) / 2;
- _Bez_point_bound(const _Bez_point_bound& o) : bounding_polyline(o.bounding_polyline),
- t_min(o.t_min), t_max(o.t_max), point_type(o.point_type), can_refine(o.can_refine)
- {}
+ bisect_control_polygon_2 (cp.begin(), cp.end(),
+ std::back_inserter (cp_a),
+ std::front_inserter (cp_b));
- _Bez_point_bound(const Control_point_vec& bp,
- const NT& min_t, const NT& max_t,
- Point_type pt, bool cr) : bounding_polyline(bp),
- t_min(min_t), t_max(max_t), point_type(pt), can_refine(cr)
- {}
+ // Check the case where t_mid is a vertical (rational) tangency
+ // point, by checking whether the first vector of right is vertical.
+ if (f_compare_x (cp_b[0], cp_b[1]) == EQUAL)
+ {
+ // The subdivision point at t_mid is a vertical tangency point with
+ // rational coordinates.
+ const typename Kernel::Point_2 &vpt = cp_b[0];
-};
+ Bez_point_bound bound (Bez_point_bound::RATIONAL_PT,
+ cp_b, t_mid, t_mid, true);
+ Bez_point_bbox bbox (vpt.x(), vpt.x(), vpt.y(), vpt.y());
-template <class NT>
-struct _Bez_point_bbox {
- NT min_x, max_x, min_y, max_y;
+ vpts.push_back (Vertical_tangency_point (bound,
+ bbox));
+ return;
+ }
- _Bez_point_bbox() : min_x(), max_x(), min_y(), max_y() {}
- _Bez_point_bbox(const _Bez_point_bbox& other) :
- min_x(other.min_x), max_x(other.max_x),
- min_y(other.min_y), max_y(other.max_y)
- {}
+ // Recursively compute the vertical tangency points of the two
+ // subcurves.
+ _compute_vertical_tangency_points (cp_a, t_min, t_mid,
+ vpts);
- _Bez_point_bbox(const NT& _min_x, const NT& _max_x,
- const NT& _min_y, const NT& _max_y) :
- min_x(_min_x), max_x(_max_x), min_y(_min_y), max_y(_max_y) {}
+ _compute_vertical_tangency_points (cp_b, t_mid, t_max,
+ vpts);
- void AddS(const _Bez_point_bbox& other) {
- min_x = std::min(min_x, other.min_x);
- min_y = std::min(min_y, other.min_y);
- max_x = std::max(max_x, other.max_x);
- max_y = std::max(max_y, other.max_y);
- }
+ return;
+ }
- _Bez_point_bbox Add(const _Bez_point_bbox& other) const {
- _Bez_point_bbox res(other);
- res.AddS(*this);
- return res;
- }
+ // If we reached here, the control polygon of the curve is y-monotone,
+ // but not x-monotone. Check whether it is convex. If it is, we know the
+ // curve contains exactly one vertical tangency point.
+ if (is_convex_2 (cp.begin(), cp.end(), m_kernel))
+ {
+ // We managed to isolate a single vertical tangency point.
+ Bez_point_bound bound (Bez_point_bound::VERTICAL_TANGENCY_PT,
+ cp, t_min, t_max,
+ true); // It is possible to refine further.
+ Bez_point_bbox vpt_bbox;
- bool Overlaps(const _Bez_point_bbox& other) const
- {
- if (max_x < other.min_x ||
- max_y < other.min_y ||
- other.max_x < min_x ||
- other.max_y < min_y)
- return false;
- return true;
- }
+ construct_bbox (cp, vpt_bbox);
+ vpts.push_back (Vertical_tangency_point (bound,
+ vpt_bbox));
- bool Overlaps_x(const _Bez_point_bbox& other) const
- {
- if (max_x < other.min_x ||
- other.max_x < min_x)
- return false;
- return true;
+ return;
+ }
+
+ // If we reached here, the curve contains more than a single vertical
+ // tangency point. We therefore bisect it and continue recursively.
+ Control_points cp_a, cp_b;
+ const NT t_mid = (t_min + t_max) / 2;
+
+ bisect_control_polygon_2 (cp.begin(), cp.end(),
+ std::back_inserter(cp_a),
+ std::front_inserter(cp_b));
+
+ // Check the case where t_mid is a vertical (rational) tangency
+ // point, by checking whether the first vector of right is vertical.
+ if (f_compare_x (cp_b[0], cp_b[1]) == EQUAL)
+ {
+ // The subdivision point at t_mid is a vertical tangency point with
+ // rational coordinates.
+ const typename Kernel::Point_2 &vpt = cp_b[0];
+
+ Bez_point_bound bound (Bez_point_bound::RATIONAL_PT,
+ cp_b, t_mid, t_mid, true);
+ Bez_point_bbox bbox (vpt.x(), vpt.x(), vpt.y(), vpt.y());
+
+ vpts.push_back (Vertical_tangency_point (bound,
+ bbox));
+ return;
+ }
+
+ // Recursively compute the vertical tangency points on the two subcurves.
+ _compute_vertical_tangency_points (cp_a, t_min, t_mid,
+ vpts);
+
+ _compute_vertical_tangency_points (cp_b, t_mid, t_max,
+ vpts);
+
+ return;
}
-
};
CGAL_END_NAMESPACE
Index: Arr_Bezier_curve_traits_2.h
===================================================================
--- Arr_Bezier_curve_traits_2.h (revision 40349)
+++ Arr_Bezier_curve_traits_2.h (working copy)
@@ -102,63 +102,60 @@
typedef typename X_monotone_curve_2::Intersection_map Intersection_map;
// Data members:
- Bezier_cache *_cache; // Caches vertical tangency points and
- // intersection points that have been
- // exactly computed.
+ Bezier_cache *p_cache; /*!< Caches vertical tangency points and
+ intersection points that have been
+ computed in an exact manner. */
+ Intersection_map *p_inter_map; /*!< Maps curve pairs to their intersection
+ points. */
+ bool m_owner; /*!< Does this instance own its cache amd
+ map structures. */
- Intersection_map *_inter_map; // Mapping curve pairs to their intersection
- // points.
-
- bool _owner; // Does the traits class own its structures.
-
public:
- /*!
- * Default constructor.
- */
+ /// \name Construction.
+ //@{
+
+ /*! Default constructor. */
Arr_Bezier_curve_traits_2 ()
{
- _cache = new Bezier_cache;
- _inter_map = new Intersection_map;
- _owner = true;
+ p_cache = new Bezier_cache;
+ p_inter_map = new Intersection_map;
+ m_owner = true;
}
-
- /*!
- * Copy constructor.
- */
+
+ /*! Copy constructor. */
Arr_Bezier_curve_traits_2 (const Self& tr) :
- _cache (tr._cache),
- _inter_map (tr._inter_map),
- _owner (false)
+ p_cache (tr.p_cache),
+ p_inter_map (tr.p_inter_map),
+ m_owner (false)
{}
- /*!
- * Assignment operator.
- */
+ /*! Assignmnet operator. */
Self& operator= (const Self& tr)
{
if (this == &tr)
return (*this);
- _cache = tr._cache;
- _inter_map = tr._inter_map;
- _owner = false;
+ p_cache = tr.p_cache;
+ p_inter_map = tr.p_inter_map;
+ m_owner = false;
return (*this);
}
- /*!
- * Destructor.
- */
+ /*! Destructor. */
~Arr_Bezier_curve_traits_2 ()
{
- if (_owner)
+ if (m_owner)
{
- delete _cache;
- delete _inter_map;
+ delete p_cache;
+ delete p_inter_map;
}
+ p_cache = NULL;
+ p_inter_map = NULL;
}
+ //@}
- /// \name Functor definitions.
+ /// \name Functor definitions for the arrangement traits.
//@{
/*! \class Compare_x_2
@@ -194,7 +191,7 @@
/*! Get a Compare_x_2 functor object. */
Compare_x_2 compare_x_2_object () const
{
- return (Compare_x_2 (_cache));
+ return (Compare_x_2 (p_cache));
}
/*! \class Compare_xy_2
@@ -230,7 +227,7 @@
/*! Get a Compare_xy_2 functor object. */
Compare_xy_2 compare_xy_2_object () const
{
- return (Compare_xy_2 (_cache));
+ return (Compare_xy_2 (p_cache));
}
/*! \class Construct_min_vertex_2
@@ -329,8 +326,6 @@
Comparison_result operator() (const Point_2& p,
const X_monotone_curve_2& cv) const
{
- // Iddo: Need to enhance point_position (maybe in the future we will
- // store subdivided subcurves).
return (cv.point_position (p,
const_cast<Bezier_cache&> (*p_cache)));
}
@@ -339,7 +334,7 @@
/*! Get a Compare_y_at_x_2 functor object. */
Compare_y_at_x_2 compare_y_at_x_2_object () const
{
- return (Compare_y_at_x_2 (_cache));
+ return (Compare_y_at_x_2 (p_cache));
}
/*! \class Compare_y_at_x_left_2
@@ -380,7 +375,7 @@
/*! Get a Compare_y_at_x_left_2 functor object. */
Compare_y_at_x_left_2 compare_y_at_x_left_2_object () const
{
- return (Compare_y_at_x_left_2 (_cache));
+ return (Compare_y_at_x_left_2 (p_cache));
}
/*! \class Compare_y_at_x_right_2
@@ -421,7 +416,7 @@
/*! Get a Compare_y_at_x_right_2 functor object. */
Compare_y_at_x_right_2 compare_y_at_x_right_2_object () const
{
- return (Compare_y_at_x_right_2 (_cache));
+ return (Compare_y_at_x_right_2 (p_cache));
}
/*! \class Equal_2
@@ -430,7 +425,7 @@
class Equal_2
{
private:
- const Bezier_cache *p_cache;
+ const Bezier_cache *p_cache;
public:
@@ -468,7 +463,7 @@
/*! Get an Equal_2 functor object. */
Equal_2 equal_2_object () const
{
- return (Equal_2 (_cache));
+ return (Equal_2 (p_cache));
}
/*! \class Make_x_monotone_2
@@ -497,53 +492,35 @@
template<class OutputIterator>
OutputIterator operator() (const Curve_2& B, OutputIterator oi)
{
+ typedef typename Bounding_traits::Vertical_tangency_point
+ Vertical_tangency_point;
+
// Try to compute the bounds of the vertical tangency points.
- Bounding_traits bound_tr;
- typename Bounding_traits::Control_point_vec cpts;
+ Bounding_traits bound_tr;
+ typename Bounding_traits::Control_points cpts;
+ std::list<Vertical_tangency_point> vpt_bounds;
std::copy (B.control_points_begin(), B.control_points_end(),
std::back_inserter(cpts));
- typedef std::pair<typename Bounding_traits::Bez_point_bound,
- typename Bounding_traits::Bez_point_bbox> Tang_pair;
+ bound_tr.compute_vertical_tangency_points
+ (cpts, std::back_inserter (vpt_bounds));
- std::list<Tang_pair> tang_bounds;
-
- bound_tr.vertical_tangency_points (cpts, 0, 1, tang_bounds);
-
- // Go over the computed bounds approximated for the vertical tangency
- // points in increasing order of their t parameters, and construct
- // Point_2 objects from the bounded tangency points.
- std::list<Point_2> tang_points;
+ // Construct Point_2 from bounded tangency points.
+ std::list<Point_2> vpts;
bool app_ok = true;
+ typename std::list<Vertical_tangency_point>::const_iterator iter;
- while (! tang_bounds.empty())
+ for (iter = vpt_bounds.begin(); iter != vpt_bounds.end(); ++iter)
{
- // Locate the bound with minimal t-value.
- typename std::list<Tang_pair>::iterator iter;
- typename std::list<Tang_pair>::iterator it_min;
+ const typename Bounding_traits::Bez_point_bound& bound = iter->bound;
+ const typename Bounding_traits::Bez_point_bbox& bbox = iter->bbox;
- it_min = iter = tang_bounds.begin();
- ++iter;
- while (iter != tang_bounds.end())
- {
- if (CGAL::compare (iter->first.t_min,
- it_min->first.t_min) == SMALLER)
- {
- it_min = iter;
- }
- ++iter;
- }
-
- // Continue with the bound for the point with minimal t-value.
- const typename Bounding_traits::Bez_point_bound& bound = it_min->first;
- const typename Bounding_traits::Bez_point_bbox& bbox = it_min->second;
-
if (! bound.can_refine)
{
// If we cannot refine the vertical-tangency bound anymore, then
// we failed to bound the vertical tangency point.
- // Iddo: In the future, we might want to use the info.
+ // \todo In the future, we might want to use the info.
app_ok = false;
break;
}
@@ -551,7 +528,7 @@
// Construct an approximate vertical tangency point.
Point_2 pt;
- if (bound.point_type == Bounding_traits::Bez_point_bound::RATIONAL_PT)
+ if (bound.type == Bounding_traits::Bez_point_bound::RATIONAL_PT)
{
CGAL_assertion (CGAL::compare (bound.t_min, bound.t_max) == EQUAL);
Rational t0 = bound.t_min;
@@ -564,10 +541,7 @@
}
pt.set_bbox (bbox);
- tang_points.push_back(pt);
-
- // Remove the tangency point with minimal t-value from the list.
- tang_bounds.erase (it_min);
+ vpts.push_back(pt);
}
// If bounding the approximated vertical-tangency points went fine,
@@ -576,20 +550,24 @@
{
// Create the x-monotone subcurves with approximate endpoints.
typename std::list<Point_2>::const_iterator pit;
- Point_2 p0(B, Rational(0)); // A rational start point.
+ Point_2 p0(B, Rational(0)); // A rational start point.
+ unsigned int xid = 1; // Serial number of the subcurve.
- for (pit = tang_points.begin(); pit != tang_points.end(); ++pit)
+ for (pit = vpts.begin(); pit != vpts.end(); ++pit)
{
- *oi = CGAL::make_object (X_monotone_curve_2 (B, p0, *pit,
+ *oi = CGAL::make_object (X_monotone_curve_2 (B, xid,
+ p0, *pit,
*p_cache));
++oi;
+ xid++;
p0 = *pit;
}
Point_2 p1(B, Rational(1)); // A rational end point.
- *oi = CGAL::make_object (X_monotone_curve_2 (B, p0, p1,
+ *oi = CGAL::make_object (X_monotone_curve_2 (B, xid,
+ p0, p1,
*p_cache));
return (oi);
@@ -607,19 +585,24 @@
Point_2 p0 (B, Rational(0));
Point_2 p1;
typename Bezier_cache::Vertical_tangency_iter it;
+ unsigned int xid = 1; // Serial number of the subcurve.
for (it = vt_list.begin(); it != vt_list.end(); ++it)
{
p1 = Point_2 (B, *it);
- *oi = CGAL::make_object (X_monotone_curve_2 (B, p0, p1,
+ *oi = CGAL::make_object (X_monotone_curve_2 (B, xid,
+ p0, p1,
*p_cache));
++oi;
+
+ xid++;
p0 = p1;
}
// Create the final subcurve.
p1 = Point_2 (B, Rational(1));
- *oi = CGAL::make_object (X_monotone_curve_2 (B, p0, p1,
+ *oi = CGAL::make_object (X_monotone_curve_2 (B, xid,
+ p0, p1,
*p_cache));
return (oi);
@@ -629,7 +612,7 @@
/*! Get a Make_x_monotone_2 functor object. */
Make_x_monotone_2 make_x_monotone_2_object ()
{
- return (Make_x_monotone_2 (_cache));
+ return (Make_x_monotone_2 (p_cache));
}
/*! \class Split_2
@@ -697,7 +680,7 @@
/*! Get an Intersect_2 functor object. */
Intersect_2 intersect_2_object ()
{
- return (Intersect_2 (_cache, _inter_map));
+ return (Intersect_2 (p_cache, p_inter_map));
}
/*! \class Are_mergeable_2
@@ -779,7 +762,7 @@
if (cv.is_directed_right())
return (SMALLER);
else
- return (LARGER);
+ return (LARGER);
}
};
- Arrangement_2: assertion violation in Arr_traits_2/Bezier_bounding_rational_traits.h, David Keller, 09/13/2007
- Re: [cgal-discuss] Arrangement_2: assertion violation in Arr_traits_2/Bezier_bounding_rational_traits.h, Efraim Fogel, 09/16/2007
- Re: [cgal-discuss] Arrangement_2: assertion violation in Arr_traits_2/Bezier_bounding_rational_traits.h, David Keller, 09/17/2007
- Re: [cgal-discuss] Arrangement_2: assertion violation in Arr_traits_2/Bezier_bounding_rational_traits.h, David Keller, 09/17/2007
- Re: [cgal-discuss] Arrangement_2: assertion violation in Arr_traits_2/Bezier_bounding_rational_traits.h, Efraim Fogel, 09/17/2007
- Re: [cgal-discuss] Arrangement_2: assertion violation in Arr_traits_2/Bezier_bounding_rational_traits.h, David Keller, 09/17/2007
- Re: [cgal-discuss] Arrangement_2: assertion violation in Arr_traits_2/Bezier_bounding_rational_traits.h, Efraim Fogel, 09/17/2007
- Re: [cgal-discuss] Arrangement_2: assertion violation in Arr_traits_2/Bezier_bounding_rational_traits.h, Efraim Fogel, 09/16/2007
Archive powered by MHonArc 2.6.16.