From 9d9fc63aac9464b0b642f1818570cf6f6ca601e9 Mon Sep 17 00:00:00 2001 From: "Johan B. C. Engelen" Date: Sat, 13 Dec 2008 19:56:16 +0000 Subject: update to 2geom rev.1723 (bzr r6996) --- src/2geom/basic-intersection.cpp | 453 +++++++++------------------------------ 1 file changed, 96 insertions(+), 357 deletions(-) (limited to 'src/2geom/basic-intersection.cpp') diff --git a/src/2geom/basic-intersection.cpp b/src/2geom/basic-intersection.cpp index 2a40e7f45..16b5c0240 100644 --- a/src/2geom/basic-intersection.cpp +++ b/src/2geom/basic-intersection.cpp @@ -1,6 +1,3 @@ - - - #include <2geom/basic-intersection.h> #include <2geom/sbasis-to-bezier.h> #include <2geom/exception.h> @@ -10,82 +7,57 @@ #include -unsigned intersect_steps = 0; - using std::vector; namespace Geom { -class OldBezier { -public: - std::vector p; - OldBezier() { - } - void split(double t, OldBezier &a, OldBezier &b) const; - Point operator()(double t) const; - - ~OldBezier() {} - - void bounds(double &minax, double &maxax, - double &minay, double &maxay) { - // Compute bounding box for a - minax = p[0][X]; // These are the most likely to be extremal - maxax = p.back()[X]; - if( minax > maxax ) - std::swap(minax, maxax); - for(unsigned i = 1; i < p.size()-1; i++) { - if( p[i][X] < minax ) - minax = p[i][X]; - else if( p[i][X] > maxax ) - maxax = p[i][X]; - } - - minay = p[0][Y]; // These are the most likely to be extremal - maxay = p.back()[Y]; - if( minay > maxay ) - std::swap(minay, maxay); - for(unsigned i = 1; i < p.size()-1; i++) { - if( p[i][Y] < minay ) - minay = p[i][Y]; - else if( p[i][Y] > maxay ) - maxay = p[i][Y]; - } - - } - -}; - -static std::vector > -find_intersections( OldBezier a, OldBezier b); +namespace detail{ namespace bezier_clipping { +void portion (std::vector & B, Interval const& I); +}; }; -static std::vector > -find_self_intersections(OldBezier const &Sb, D2 const & A); - -std::vector > -find_intersections( vector const & A, - vector const & B) { - OldBezier a, b; - a.p = A; - b.p = B; - return find_intersections(a,b); +void find_intersections(std::vector > &xs, + D2 const & A, + D2 const & B) { + vector BezA, BezB; + sbasis_to_bezier(BezA, A); + sbasis_to_bezier(BezB, B); + + xs.clear(); + + find_intersections(xs, BezA, BezB); } -std::vector > -find_self_intersections(OldBezier const &/*Sb*/) { - THROW_NOTIMPLEMENTED(); -} +/* + * split the curve at the midpoint, returning an array with the two parts + * Temporary storage is minimized by using part of the storage for the result + * to hold an intermediate value until it is no longer needed. + */ +void split(vector const &p, double t, + vector &left, vector &right) { + const unsigned sz = p.size(); + Geom::Point Vtemp[sz][sz]; -std::vector > -find_self_intersections(D2 const & A) { - OldBezier Sb; - sbasis_to_bezier(Sb.p, A); - return find_self_intersections(Sb, A); -} + /* Copy control points */ + std::copy(p.begin(), p.end(), Vtemp[0]); + /* Triangle computation */ + for (unsigned i = 1; i < sz; i++) { + for (unsigned j = 0; j < sz - i; j++) { + Vtemp[i][j] = lerp(t, Vtemp[i-1][j], Vtemp[i-1][j+1]); + } + } -static std::vector > -find_self_intersections(OldBezier const &Sb, D2 const & A) { + left.resize(sz); + right.resize(sz); + for (unsigned j = 0; j < sz; j++) + left[j] = Vtemp[j][0]; + for (unsigned j = 0; j < sz; j++) + right[j] = Vtemp[sz-1-j][j]; +} +void +find_self_intersections(std::vector > &xs, + D2 const & A) { vector dr = roots(derivative(A[X])); { vector dyr = roots(derivative(A[Y])); @@ -97,21 +69,21 @@ find_self_intersections(OldBezier const &Sb, D2 const & A) { sort(dr.begin(), dr.end()); unique(dr.begin(), dr.end()); - std::vector > all_si; - - vector pieces; + vector > pieces; { - OldBezier in = Sb, l, r; + vector in, l, r; + sbasis_to_bezier(in, A); for(unsigned i = 0; i < dr.size()-1; i++) { - in.split((dr[i+1]-dr[i]) / (1 - dr[i]), l, r); + split(in, (dr[i+1]-dr[i]) / (1 - dr[i]), l, r); pieces.push_back(l); in = r; } } for(unsigned i = 0; i < dr.size()-1; i++) { for(unsigned j = i+1; j < dr.size()-1; j++) { - std::vector > section = - find_intersections( pieces[i], pieces[j]); + std::vector > section; + + find_intersections( section, pieces[i], pieces[j]); for(unsigned k = 0; k < section.size(); k++) { double l = section[k].first; double r = section[k].second; @@ -119,287 +91,31 @@ find_self_intersections(OldBezier const &Sb, D2 const & A) { if(j == i+1) if((l == 1) && (r == 0)) continue; - all_si.push_back(std::make_pair((1-l)*dr[i] + l*dr[i+1], + xs.push_back(std::make_pair((1-l)*dr[i] + l*dr[i+1], (1-r)*dr[j] + r*dr[j+1])); } } } - // Because i is in order, all_si should be roughly already in order? - //sort(all_si.begin(), all_si.end()); - //unique(all_si.begin(), all_si.end()); - - return all_si; -} - -/* The value of 1.0 / (1L<<14) is enough for most applications */ -const double INV_EPS = (1L<<14); - -/* - * split the curve at the midpoint, returning an array with the two parts - * Temporary storage is minimized by using part of the storage for the result - * to hold an intermediate value until it is no longer needed. - */ -void OldBezier::split(double t, OldBezier &left, OldBezier &right) const { - const unsigned sz = p.size(); - Geom::Point Vtemp[sz][sz]; - - /* Copy control points */ - std::copy(p.begin(), p.end(), Vtemp[0]); - - /* Triangle computation */ - for (unsigned i = 1; i < sz; i++) { - for (unsigned j = 0; j < sz - i; j++) { - Vtemp[i][j] = lerp(t, Vtemp[i-1][j], Vtemp[i-1][j+1]); - } - } - - left.p.resize(sz); - right.p.resize(sz); - for (unsigned j = 0; j < sz; j++) - left.p[j] = Vtemp[j][0]; - for (unsigned j = 0; j < sz; j++) - right.p[j] = Vtemp[sz-1-j][j]; -} - -#if 0 -/* - * split the curve at the midpoint, returning an array with the two parts - * Temporary storage is minimized by using part of the storage for the result - * to hold an intermediate value until it is no longer needed. - */ -Point OldBezier::operator()(double t) const { - const unsigned sz = p.size(); - Geom::Point Vtemp[sz][sz]; - - /* Copy control points */ - std::copy(p.begin(), p.end(), Vtemp[0]); - - /* Triangle computation */ - for (unsigned i = 1; i < sz; i++) { - for (unsigned j = 0; j < sz - i; j++) { - Vtemp[i][j] = lerp(t, Vtemp[i-1][j], Vtemp[i-1][j+1]); - } - } - return Vtemp[sz-1][0]; -} -#endif - -// suggested by Sederberg. -Point OldBezier::operator()(double t) const { - int n = p.size()-1; - double u, bc, tn, tmp; - int i; - Point r; - for(int dim = 0; dim < 2; dim++) { - u = 1.0 - t; - bc = 1; - tn = 1; - tmp = p[0][dim]*u; - for(i=1; i= : need boundary case - return not( ( minax > maxbx ) || ( minay > maxby ) - || ( minbx > maxax ) || ( minby > maxay ) ); -} - -/* - * Recursively intersect two curves keeping track of their real parameters - * and depths of intersection. - * The results are returned in a 2-D array of doubles indicating the parameters - * for which intersections are found. The parameters are in the order the - * intersections were found, which is probably not in sorted order. - * When an intersection is found, the parameter value for each of the two - * is stored in the index elements array, and the index is incremented. - * - * If either of the curves has subdivisions left before it is straight - * (depth > 0) - * that curve (possibly both) is (are) subdivided at its (their) midpoint(s). - * the depth(s) is (are) decremented, and the parameter value(s) corresponding - * to the midpoints(s) is (are) computed. - * Then each of the subcurves of one curve is intersected with each of the - * subcurves of the other curve, first by testing the bounding boxes for - * interference. If there is any bounding box interference, the corresponding - * subcurves are recursively intersected. - * - * If neither curve has subdivisions left, the line segments from the first - * to last control point of each segment are intersected. (Actually the - * only the parameter value corresponding to the intersection point is found). - * - * The apriori flatness test is probably more efficient than testing at each - * level of recursion, although a test after three or four levels would - * probably be worthwhile, since many curves become flat faster than their - * asymptotic rate for the first few levels of recursion. - * - * The bounding box test fails much more frequently than it succeeds, providing - * substantial pruning of the search space. - * - * Each (sub)curve is subdivided only once, hence it is not possible that for - * one final line intersection test the subdivision was at one level, while - * for another final line intersection test the subdivision (of the same curve) - * was at another. Since the line segments share endpoints, the intersection - * is robust: a near-tangential intersection will yield zero or two - * intersections. - */ -void recursively_intersect( OldBezier a, double t0, double t1, int deptha, - OldBezier b, double u0, double u1, int depthb, - std::vector > ¶meters) -{ - intersect_steps ++; - if( deptha > 0 ) - { - OldBezier A[2]; - a.split(0.5, A[0], A[1]); - double tmid = (t0+t1)*0.5; - deptha--; - if( depthb > 0 ) - { - OldBezier B[2]; - b.split(0.5, B[0], B[1]); - double umid = (u0+u1)*0.5; - depthb--; - if( intersect_BB( A[0], B[0] ) ) - recursively_intersect( A[0], t0, tmid, deptha, - B[0], u0, umid, depthb, - parameters ); - if( intersect_BB( A[1], B[0] ) ) - recursively_intersect( A[1], tmid, t1, deptha, - B[0], u0, umid, depthb, - parameters ); - if( intersect_BB( A[0], B[1] ) ) - recursively_intersect( A[0], t0, tmid, deptha, - B[1], umid, u1, depthb, - parameters ); - if( intersect_BB( A[1], B[1] ) ) - recursively_intersect( A[1], tmid, t1, deptha, - B[1], umid, u1, depthb, - parameters ); - } - else - { - if( intersect_BB( A[0], b ) ) - recursively_intersect( A[0], t0, tmid, deptha, - b, u0, u1, depthb, - parameters ); - if( intersect_BB( A[1], b ) ) - recursively_intersect( A[1], tmid, t1, deptha, - b, u0, u1, depthb, - parameters ); - } - } - else - if( depthb > 0 ) - { - OldBezier B[2]; - b.split(0.5, B[0], B[1]); - double umid = (u0 + u1)*0.5; - depthb--; - if( intersect_BB( a, B[0] ) ) - recursively_intersect( a, t0, t1, deptha, - B[0], u0, umid, depthb, - parameters ); - if( intersect_BB( a, B[1] ) ) - recursively_intersect( a, t0, t1, deptha, - B[0], umid, u1, depthb, - parameters ); - } - else // Both segments are fully subdivided; now do line segments - { - double xlk = a.p.back()[X] - a.p[0][X]; - double ylk = a.p.back()[Y] - a.p[0][Y]; - double xnm = b.p.back()[X] - b.p[0][X]; - double ynm = b.p.back()[Y] - b.p[0][Y]; - double xmk = b.p[0][X] - a.p[0][X]; - double ymk = b.p[0][Y] - a.p[0][Y]; - double det = xnm * ylk - ynm * xlk; - if( 1.0 + det == 1.0 ) - return; - else - { - double detinv = 1.0 / det; - double s = ( xnm * ymk - ynm *xmk ) * detinv; - double t = ( xlk * ymk - ylk * xmk ) * detinv; - if( ( s < 0.0 ) || ( s > 1.0 ) || ( t < 0.0 ) || ( t > 1.0 ) ) - return; - parameters.push_back(std::pair(t0 + s * ( t1 - t0 ), - u0 + t * ( u1 - u0 ))); - } - } -} - -inline double log4( double x ) { return log(x)/log(4.); } - -/* - * Wang's theorem is used to estimate the level of subdivision required, - * but only if the bounding boxes interfere at the top level. - * Assuming there is a possible intersection, recursively_intersect is - * used to find all the parameters corresponding to intersection points. - * these are then sorted and returned in an array. - */ - -double Lmax(Point p) { - return std::max(fabs(p[X]), fabs(p[Y])); -} - -unsigned wangs_theorem(OldBezier a) { - return 6; // seems a good approximation! - double la1 = Lmax( ( a.p[2] - a.p[1] ) - (a.p[1] - a.p[0]) ); - double la2 = Lmax( ( a.p[3] - a.p[2] ) - (a.p[2] - a.p[1]) ); - double l0 = std::max(la1, la2); - unsigned ra; - if( l0 * 0.75 * M_SQRT2 + 1.0 == 1.0 ) - ra = 0; - else - ra = (unsigned)ceil( log4( M_SQRT2 * 6.0 / 8.0 * INV_EPS * l0 ) ); - std::cout << ra << std::endl; - return ra; -} +#include + + struct rparams { - OldBezier &A; - OldBezier &B; + D2 const &A; + D2 const &B; }; static int intersect_polish_f (const gsl_vector * x, void *params, - gsl_vector * f) + gsl_vector * f) { const double x0 = gsl_vector_get (x, 0); const double x1 = gsl_vector_get (x, 1); @@ -413,7 +129,7 @@ intersect_polish_f (const gsl_vector * x, void *params, return GSL_SUCCESS; } -typedef union dbl_64{ +union dbl_64{ long long i64; double d64; }; @@ -427,8 +143,8 @@ static double EpsilonBy(double value, int eps) } -static void intersect_polish_root (OldBezier &A, double &s, - OldBezier &B, double &t) { +static void intersect_polish_root (D2 const &A, double &s, + D2 const &B, double &t) { const gsl_multiroot_fsolver_type *T; gsl_multiroot_fsolver *sol; @@ -502,23 +218,46 @@ static void intersect_polish_root (OldBezier &A, double &s, } -std::vector > find_intersections( OldBezier a, OldBezier b) +void polish_intersections(std::vector > &xs, + D2 const &A, D2 const &B) { - std::vector > parameters; - if( intersect_BB( a, b ) ) - { - recursively_intersect( a, 0., 1., wangs_theorem(a), - b, 0., 1., wangs_theorem(b), - parameters); - } - for(unsigned i = 0; i < parameters.size(); i++) - intersect_polish_root(a, parameters[i].first, - b, parameters[i].second); - std::sort(parameters.begin(), parameters.end()); - return parameters; + for(unsigned i = 0; i < xs.size(); i++) + intersect_polish_root(A, xs[i].first, + B, xs[i].second); } + /** + * Compute the Hausdorf distance from A to B only. + */ + + +#if 0 +/** Compute the value of a bezier + Todo: find a good palce for this. + */ +// suggested by Sederberg. +Point OldBezier::operator()(double t) const { + int n = p.size()-1; + double u, bc, tn, tmp; + int i; + Point r; + for(int dim = 0; dim < 2; dim++) { + u = 1.0 - t; + bc = 1; + tn = 1; + tmp = p[0][dim]*u; + for(i=1; i