From c3a8ad9235ff81909bd472707550aef5b91daf7b Mon Sep 17 00:00:00 2001 From: johanengelen Date: Sat, 13 Dec 2008 19:56:16 +0000 Subject: [PATCH] update to 2geom rev.1723 --- src/2geom/basic-intersection.cpp | 453 +++++++------------------------ src/2geom/basic-intersection.h | 54 +++- src/2geom/bezier-curve.h | 4 + src/2geom/circle.h | 1 - src/2geom/convex-cover.cpp | 33 ++- src/2geom/convex-cover.h | 3 + src/2geom/crossing.h | 8 +- src/2geom/curve.h | 6 +- src/2geom/d2-sbasis.cpp | 8 +- src/2geom/d2.h | 6 + src/2geom/ellipse.cpp | 6 + src/2geom/ellipse.h | 3 + src/2geom/elliptical-arc.h | 2 + src/2geom/geom.cpp | 32 ++- src/2geom/geom.h | 21 +- src/2geom/hvlinesegment.h | 4 + src/2geom/interval.h | 10 +- src/2geom/linear.h | 30 +- src/2geom/path-intersection.cpp | 6 +- src/2geom/pathvector.h | 1 - src/2geom/piecewise.h | 13 + src/2geom/sbasis-2d.cpp | 16 +- src/2geom/sbasis-curve.h | 2 + src/2geom/sbasis-geometric.cpp | 57 ++-- src/2geom/sbasis.cpp | 106 ++++---- src/2geom/sbasis.h | 74 +++-- src/2geom/svg-elliptical-arc.h | 2 + src/2geom/sweep.cpp | 18 ++ 28 files changed, 452 insertions(+), 527 deletions(-) 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 #include +#define USE_RECURSIVE_INTERSECTOR 0 namespace Geom { -std::vector > -find_intersections( D2 const & A, +#ifdef USE_RECURSIVE_INTERSECTOR +void +find_intersections( std::vector > &xs, + D2 const & A, D2 const & B); -std::vector > -find_self_intersections(D2 const & A); +void +find_self_intersections(std::vector > &xs, + D2 const & A); // Bezier form -std::vector > -find_intersections( std::vector const & A, - std::vector const & B); +void +find_intersections( std::vector > &xs, + std::vector const & A, + std::vector const & B, + double precision = 1e-5); -std::vector > -find_self_intersections(std::vector const & A); +// Bezier form +void +find_intersections_bezier_clipping( std::vector > &xs, + std::vector const & A, + std::vector const & B, + double precision = 1e-5); + +void +find_self_intersections(std::vector > &xs, + std::vector const & A); +#else /* * find_intersection @@ -72,11 +87,11 @@ find_self_intersections(std::vector const & A); * This routine is based on the Bezier Clipping Algorithm, * see: Sederberg, Nishita, 1990 - Curve intersection using Bezier clipping */ -void find_intersections (std::vector< std::pair > & xs, +void find_intersections_clipping (std::vector< std::pair > & xs, std::vector const& A, std::vector const& B, double precision = 1e-5); - +#endif /* * find_collinear_normal * @@ -93,6 +108,23 @@ void find_collinear_normal (std::vector< std::pair >& xs, std::vector const& B, double precision = 1e-5); +void polish_intersections(std::vector > &xs, + D2 const &A, + D2 const &B); + +void find_intersections(std::vector >& xs, + D2 const & A, + D2 const & B); + +void find_self_intersections(std::vector >& xs, + D2 const & A); + + +void find_self_intersections(std::vector >& xs, + std::vector const & A); + + + /** * Compute the Hausdorf distance from A to B only. */ diff --git a/src/2geom/bezier-curve.h b/src/2geom/bezier-curve.h index 589335d22..acf878ee8 100644 --- a/src/2geom/bezier-curve.h +++ b/src/2geom/bezier-curve.h @@ -121,6 +121,10 @@ public: int winding(Point p) const { return SBasisCurve(toSBasis()).winding(p); } + + virtual int degreesOfFreedom() const { + return 2*order; + } std::vector roots(double v, Dim2 d) const { diff --git a/src/2geom/circle.h b/src/2geom/circle.h index 1f9871276..27d4fcc3f 100644 --- a/src/2geom/circle.h +++ b/src/2geom/circle.h @@ -105,7 +105,6 @@ class Circle private: Point m_centre; Coord m_ray; - Coord m_angle; }; diff --git a/src/2geom/convex-cover.cpp b/src/2geom/convex-cover.cpp index c4d8df338..e8ea2280d 100644 --- a/src/2geom/convex-cover.cpp +++ b/src/2geom/convex-cover.cpp @@ -179,6 +179,13 @@ int mod(int i, int l) { * Tests if a point is left (outside) of a particular segment, n. */ bool ConvexHull::is_left(Point p, int n) { + return SignedTriangleArea((*this)[n], (*this)[n+1], p) >= 0; +} + +/*** ConvexHull::strict_left + * Tests if a point is left (outside) of a particular segment, n. */ +bool +ConvexHull::is_strict_left(Point p, int n) { return SignedTriangleArea((*this)[n], (*this)[n+1], p) > 0; } @@ -193,6 +200,20 @@ ConvexHull::find_left(Point p) { } return -1; } + + +/*** ConvexHull::find_positive + * May return any number n where the segment n -> n + 1 (possibly looped around) in the hull such + * that the point is on the wrong side to be within the hull. Returns -1 if it is within the hull.*/ +int +ConvexHull::find_strict_left(Point p) { + int l = boundary.size(); //Who knows if C++ is smart enough to optimize this? + for(int i = 0; i < l; i++) { + if(is_strict_left(p, i)) return i; + } + return -1; +} + //OPT: do a spread iteration - quasi-random with no repeats and full coverage. /*** ConvexHull::contains_point @@ -203,6 +224,14 @@ ConvexHull::contains_point(Point p) { return find_left(p) == -1; } +/*** ConvexHull::strict_contains_point + * In order to test whether a point is strictly inside (not on the boundary) a convex hull we can travel once around the outside making + * sure that each triangle made from an edge and the point has positive area. */ +bool +ConvexHull::strict_contains_point(Point p) { + return find_strict_left(p) == -1; +} + /*** ConvexHull::add_point * to add a point we need to find whether the new point extends the boundary, and if so, what it * obscures. Tarjan? Jarvis?*/ @@ -219,9 +248,9 @@ ConvexHull::merge(Point p) { bool pushed = false; - bool pre = is_left(p, -1); + bool pre = is_strict_left(p, -1); for(int i = 0; i < l; i++) { - bool cur = is_left(p, i); + bool cur = is_strict_left(p, i); if(pre) { if(cur) { if(!pushed) { diff --git a/src/2geom/convex-cover.h b/src/2geom/convex-cover.h index f22177746..524108965 100644 --- a/src/2geom/convex-cover.h +++ b/src/2geom/convex-cover.h @@ -65,6 +65,7 @@ public: void merge(Point p); bool contains_point(Point p); + bool strict_contains_point(Point p); inline Point operator[](int i) const { @@ -122,7 +123,9 @@ public: Point const * furthest(Point direction) const; bool is_left(Point p, int n); + bool is_strict_left(Point p, int n); int find_left(Point p); + int find_strict_left(Point p); double narrowest_diameter(Point &a, Point &b, Point &c); }; diff --git a/src/2geom/crossing.h b/src/2geom/crossing.h index 3ad034ac9..7737b24a9 100644 --- a/src/2geom/crossing.h +++ b/src/2geom/crossing.h @@ -3,9 +3,10 @@ * \brief \todo brief description * * Authors: - * ? + * Michael Sloane + * Marco * - * Copyright ?-? authors + * Copyright 2006-2008 authors * * This library is free software; you can redistribute it and/or * modify it either under the terms of the GNU Lesser General Public @@ -38,6 +39,7 @@ #include #include <2geom/rect.h> #include <2geom/sweep.h> +#include namespace Geom { @@ -58,6 +60,8 @@ struct Crossing { bool onIx(unsigned ix) const { return a == ix || b == ix; } }; +typedef boost::optional OptCrossing; + /* struct Edge { diff --git a/src/2geom/curve.h b/src/2geom/curve.h index 08cc2380c..ce1fec3a6 100644 --- a/src/2geom/curve.h +++ b/src/2geom/curve.h @@ -83,6 +83,8 @@ public: virtual int winding(Point p) const { return root_winding(*this, p); } + virtual int degreesOfFreedom() const { return 0;} + //mental: review these virtual Curve *portion(double f, double t) const = 0; virtual Curve *reverse() const { return portion(1, 0); } @@ -133,8 +135,8 @@ public: */ virtual Point unitTangentAt(Coord t, unsigned n = 3) const { - std::vector derivs = pointAndDerivatives(t, n); - for (unsigned deriv_n = 1; deriv_n < derivs.size(); deriv_n++) { + std::vector derivs = pointAndDerivatives(t, n); + for (unsigned deriv_n = 1; deriv_n < derivs.size(); deriv_n++) { Coord length = derivs[deriv_n].length(); if ( ! are_near(length, 0) ) { // length of derivative is non-zero, so return unit vector diff --git a/src/2geom/d2-sbasis.cpp b/src/2geom/d2-sbasis.cpp index 55c7ef55e..2c52e4782 100644 --- a/src/2geom/d2-sbasis.cpp +++ b/src/2geom/d2-sbasis.cpp @@ -113,12 +113,12 @@ Piecewise > force_continuity(Piecewise > const &f, double SBasis &cur_sb =result.segs[cur][dim]; Coord const c=pt0[dim]; if (prev_sb.empty()) { - prev_sb.push_back(Linear(0.0, c)); + prev_sb = SBasis(Linear(0.0, c)); } else { prev_sb[0][1] = c; } if (cur_sb.empty()) { - cur_sb.push_back(Linear(c, 0.0)); + cur_sb = SBasis(Linear(c, 0.0)); } else { cur_sb[0][0] = c; } @@ -157,7 +157,7 @@ static void set_first_point(Piecewise > &f, Point a){ } for (unsigned dim=0; dim<2; dim++){ if (f.segs.front()[dim].size() == 0){ - f.segs.front()[dim].push_back(Linear(a[dim],0)); + f.segs.front()[dim] = SBasis(Linear(a[dim],0)); }else{ f.segs.front()[dim][0][0] = a[dim]; } @@ -170,7 +170,7 @@ static void set_last_point(Piecewise > &f, Point a){ } for (unsigned dim=0; dim<2; dim++){ if (f.segs.back()[dim].size() == 0){ - f.segs.back()[dim].push_back(Linear(0,a[dim])); + f.segs.back()[dim] = SBasis(Linear(0,a[dim])); }else{ f.segs.back()[dim][0][1] = a[dim]; } diff --git a/src/2geom/d2.h b/src/2geom/d2.h index a14e3b0eb..f55f6a596 100644 --- a/src/2geom/d2.h +++ b/src/2geom/d2.h @@ -124,6 +124,12 @@ inline D2 portion(const D2 &a, Coord f, Coord t) { return D2(portion(a[X], f, t), portion(a[Y], f, t)); } +template +inline D2 portion(const D2 &a, Interval i) { + boost::function_requires >(); + return D2(portion(a[X], i), portion(a[Y], i)); +} + //IMPL: boost::EqualityComparableConcept template inline bool diff --git a/src/2geom/ellipse.cpp b/src/2geom/ellipse.cpp index 2690a091b..10071d09a 100644 --- a/src/2geom/ellipse.cpp +++ b/src/2geom/ellipse.cpp @@ -264,6 +264,12 @@ Ellipse Ellipse::transformed(Matrix const& m) const return e; } +Ellipse::Ellipse(Geom::Circle const &c) +{ + m_centre = c.center(); + m_ray[X] = m_ray[Y] = c.ray(); +} + } // end namespace Geom diff --git a/src/2geom/ellipse.h b/src/2geom/ellipse.h index d5b882bc4..7ed04e51b 100644 --- a/src/2geom/ellipse.h +++ b/src/2geom/ellipse.h @@ -44,6 +44,7 @@ namespace Geom { class SVGEllipticalArc; +class Circle; class Ellipse { @@ -67,6 +68,8 @@ class Ellipse { set(points); } + + Ellipse(Geom::Circle const &c); void set(double cx, double cy, double rx, double ry, double a) { diff --git a/src/2geom/elliptical-arc.h b/src/2geom/elliptical-arc.h index 25c79a363..b0c0bd9df 100644 --- a/src/2geom/elliptical-arc.h +++ b/src/2geom/elliptical-arc.h @@ -200,6 +200,8 @@ class EllipticalArc : public Curve { return SBasisCurve(toSBasis()).winding(p); } + + int degreesOfFreedom() const { return 7;} Curve *derivative() const; diff --git a/src/2geom/geom.cpp b/src/2geom/geom.cpp index f9b1a664b..5eade57f2 100644 --- a/src/2geom/geom.cpp +++ b/src/2geom/geom.cpp @@ -1,5 +1,4 @@ /** - * \file geom.cpp * \brief Various geometrical calculations. */ @@ -85,7 +84,6 @@ line_intersection(Geom::Point const &n0, double const d0, - /* ccw exists as a building block */ int intersector_ccw(const Geom::Point& p0, const Geom::Point& p1, @@ -314,13 +312,35 @@ rect_line_intersect(Geom::Point const &c0, Geom::Point const &c1, return results; } -std::vector +/** Determine whether \& where an (infinite) line intersects a rectangle. + * + * \a c0, \a c1 are diagonal corners of the rectangle and + * \a p1, \a p1 are distinct points on the line + * + * \return A list (possibly empty) of points of intersection. If two such points (say \a r0 and \a + * r1) then it is guaranteed that the order of \a r0, \a r1 along the line is the same as the that + * of \a c0, \a c1 (i.e., the vectors \a r1 - \a r0 and \a p1 - \a p0 point into the same + * direction). + */ +boost::optional rect_line_intersect(Geom::Rect &r, - Geom::Point const &p0, Geom::Point const &p1) + Geom::LineSegment ls) { - return rect_line_intersect(r.min(), r.max(), p0, p1); + std::vector results; + + results = rect_line_intersect(r.min(), r.max(), ls[0], ls[1]); + if(results.size() == 2) { + return LineSegment(results[0], results[1]); + } + return boost::optional(); } +boost::optional +rect_line_intersect(Geom::Rect &r, + Geom::Line l) +{ + return rect_line_intersect(r, l.segment(0, 1)); +} /** * polyCentroid: Calculates the centroid (xCentroid, yCentroid) and area of a polygon, given its @@ -336,7 +356,7 @@ rect_line_intersect(Geom::Rect &r, * for now we require the path to be a polyline and assume it is closed. **/ -int centroid(std::vector p, Geom::Point& centroid, double &area) { +int centroid(std::vector const &p, Geom::Point& centroid, double &area) { const unsigned n = p.size(); if (n < 3) return 1; diff --git a/src/2geom/geom.h b/src/2geom/geom.h index d0af7d7d2..9233696d7 100644 --- a/src/2geom/geom.h +++ b/src/2geom/geom.h @@ -37,8 +37,10 @@ //TODO: move somewhere else #include -#include <2geom/point.h> -#include <2geom/rect.h> +#include <2geom/forward.h> +#include +#include <2geom/bezier-curve.h> +#include <2geom/line.h> namespace Geom { @@ -55,6 +57,9 @@ intersector_ccw(const Geom::Point& p0, const Geom::Point& p1, /* intersectors */ +#if 0 +// Use the new routines provided in line.h + IntersectorKind line_intersection(Geom::Point const &n0, double const d0, Geom::Point const &n1, double const d1, @@ -69,17 +74,23 @@ IntersectorKind line_twopoint_intersect(Geom::Point const &p00, Geom::Point const &p01, Geom::Point const &p10, Geom::Point const &p11, Geom::Point &result); +#endif +#if 0 std::vector rect_line_intersect(Geom::Point const &E, Geom::Point const &F, Geom::Point const &p0, Geom::Point const &p1); +#endif +boost::optional +rect_line_intersect(Geom::Rect &r, + Geom::LineSegment ls); -std::vector +boost::optional rect_line_intersect(Geom::Rect &r, - Geom::Point const &p0, Geom::Point const &p1); + Geom::Line l); -int centroid(std::vector p, Geom::Point& centroid, double &area); +int centroid(std::vector const &p, Geom::Point& centroid, double &area); } diff --git a/src/2geom/hvlinesegment.h b/src/2geom/hvlinesegment.h index ccb95afdb..ac91ec80a 100644 --- a/src/2geom/hvlinesegment.h +++ b/src/2geom/hvlinesegment.h @@ -136,6 +136,8 @@ class HLineSegment : public Curve return m_line_seg.winding(p); } + int degreesOfFreedom() const { return 3;} + std::vector roots(double v, Dim2 d) const { @@ -375,6 +377,8 @@ class VLineSegment : public Curve return m_line_seg.winding(p); } + int degreesOfFreedom() const { return 3;} + std::vector roots(double v, Dim2 d) const { diff --git a/src/2geom/interval.h b/src/2geom/interval.h index 009528586..d4dae41b4 100644 --- a/src/2geom/interval.h +++ b/src/2geom/interval.h @@ -258,10 +258,18 @@ inline OptInterval intersect(const Interval & a, const Interval & b) { Coord u = std::max(a.min(), b.min()), v = std::min(a.max(), b.max()); //technically >= might be incorrect, but singulars suck - return u >= v ? OptInterval() + return u > v ? OptInterval() : OptInterval(Interval(u, v)); } +#ifdef _GLIBCXX_IOSTREAM +inline std::ostream &operator<< (std::ostream &os, + const Geom::Interval &I) { + os << "Interval("<= 0); @@ -113,10 +89,10 @@ public: inline OptInterval bounds_fast() const { return bounds_exact(); } inline OptInterval bounds_local(double u, double v) const { return Interval(valueAt(u), valueAt(v)); } - operator Tri() const { + double tri() const { return a[1] - a[0]; } - operator Hat() const { + double hat() const { return (a[1] + a[0])/2; } }; diff --git a/src/2geom/path-intersection.cpp b/src/2geom/path-intersection.cpp index 5a4e76020..bfb1fb96c 100644 --- a/src/2geom/path-intersection.cpp +++ b/src/2geom/path-intersection.cpp @@ -78,7 +78,7 @@ int winding(Path const &path, Point p) { const double fudge = 0.01; if(cmp(y, next->valueAt(fudge, Y)) == initial_to_ray) { wind += int(c); - std::cout << "!!!!!" << int(c) << " "; + //std::cout << "!!!!!" << int(c) << " "; } iter = next; // No increment, as the rest of the thing hasn't been counted. } else { @@ -86,7 +86,7 @@ int winding(Path const &path, Point p) { if(cmp(y, ny) == initial_to_ray) { //Is a continuation through the ray, so counts windingwise wind += int(c); - std::cout << "!!!!!" << int(c) << " "; + //std::cout << "!!!!!" << int(c) << " "; } iter = ++next; } @@ -177,6 +177,7 @@ linear_intersect(Point A0, Point A1, Point B0, Point B1, } +#if 0 typedef union dbl_64{ long long i64; double d64; @@ -196,6 +197,7 @@ static double EpsilonOf(double value) else return value - s.d64; } +#endif struct rparams { Curve const &A; diff --git a/src/2geom/pathvector.h b/src/2geom/pathvector.h index 9efae7c73..d1d785a07 100644 --- a/src/2geom/pathvector.h +++ b/src/2geom/pathvector.h @@ -38,7 +38,6 @@ #include <2geom/forward.h> #include <2geom/path.h> -#include <2geom/rect.h> #include <2geom/transforms.h> namespace Geom { diff --git a/src/2geom/piecewise.h b/src/2geom/piecewise.h index bbd1f054a..31bf6872a 100644 --- a/src/2geom/piecewise.h +++ b/src/2geom/piecewise.h @@ -754,8 +754,21 @@ Piecewise reverse(Piecewise const &f) { } +/** + * Interpolates between a and b. + * \return a if t = 0, b if t = 1, or an interpolation between a and b for t in [0,1] + */ +template +Piecewise lerp(Piecewise const &a, Piecewise b, double t) { + // Make sure both paths have the same number of segments and cuts at the same locations + b.setDomain(a.domain()); + Piecewise pA = partition(a, b.cuts); + Piecewise pB = partition(b, a.cuts); + + return (pA*(1-t) + pB*t); } +} #endif //SEEN_GEOM_PW_SB_H /* Local Variables: diff --git a/src/2geom/sbasis-2d.cpp b/src/2geom/sbasis-2d.cpp index e3a407094..399fb8595 100644 --- a/src/2geom/sbasis-2d.cpp +++ b/src/2geom/sbasis-2d.cpp @@ -4,7 +4,7 @@ namespace Geom{ SBasis extract_u(SBasis2d const &a, double u) { - SBasis sb; + SBasis sb(a.vs, Linear()); double s = u*(1-u); for(unsigned vi = 0; vi < a.vs; vi++) { @@ -14,14 +14,14 @@ SBasis extract_u(SBasis2d const &a, double u) { bo += (extract_u(a.index(ui, vi), u))*sk; sk *= s; } - sb.push_back(bo); + sb[vi] = bo; } return sb; } SBasis extract_v(SBasis2d const &a, double v) { - SBasis sb; + SBasis sb(a.us, Linear()); double s = v*(1-v); for(unsigned ui = 0; ui < a.us; ui++) { @@ -31,7 +31,7 @@ SBasis extract_v(SBasis2d const &a, double v) { bo += (extract_v(a.index(ui, vi), v))*sk; sk *= s; } - sb.push_back(bo); + sb[ui] = bo; } return sb; @@ -105,7 +105,6 @@ SBasis2d partial_derivative(SBasis2d const &f, int dim) { */ D2 sb2dsolve(SBasis2d const &f, Geom::Point const &A, Geom::Point const &B, unsigned degmax){ - D2result(Linear(A[X],B[X]),Linear(A[Y],B[Y])); //g_warning("check f(A)= %f = f(B) = %f =0!", f.apply(A[X],A[Y]), f.apply(B[X],B[Y])); SBasis2d dfdu = partial_derivative(f, 0); @@ -115,8 +114,11 @@ sb2dsolve(SBasis2d const &f, Geom::Point const &A, Geom::Point const &B, unsigne Geom::Point nA = dfA/(dfA[X]*dfA[X]+dfA[Y]*dfA[Y]); Geom::Point nB = dfB/(dfB[X]*dfB[X]+dfB[Y]*dfB[Y]); + D2result(SBasis(degmax, Linear()), SBasis(degmax, Linear())); double fact_k=1; double sign = 1.; + for(int dim = 0; dim < 2; dim++) + result[dim][0] = Linear(A[dim],B[dim]); for(unsigned k=1; k toSBasis() const { return inner; } + virtual int degreesOfFreedom() const { return inner[0].degreesOfFreedom() + inner[1].degreesOfFreedom(); + } }; diff --git a/src/2geom/sbasis-geometric.cpp b/src/2geom/sbasis-geometric.cpp index 96a5ed0ce..d27255749 100644 --- a/src/2geom/sbasis-geometric.cpp +++ b/src/2geom/sbasis-geometric.cpp @@ -27,6 +27,11 @@ using namespace std; //Some utils first. //TODO: remove this!! +/** + * Return a list of doubles that appear in both a and b to within error tol + * a, b, vector of double + * tol tolerance + */ static vector vect_intersect(vector const &a, vector const &b, double tol=0.){ vector inter; @@ -72,21 +77,7 @@ static SBasis divide_by_t0k(SBasis const &a, int k) { } static SBasis divide_by_t1k(SBasis const &a, int k) { - if(k < 0) { - SBasis c = Linear(1,0); - for (int i=2; i<=-k; i++){ - c*=c; - } - c*=a; - return(c); - }else{ - SBasis c = Linear(0,1); - for (int i=2; i<=k; i++){ - c*=c; - } - c*=a; - return(divide_by_sk(c,k)); - } + return divide_by_t0k(a, -k); } static D2 RescaleForNonVanishingEnds(D2 const &MM, double ZERO=1.e-4){ @@ -110,6 +101,14 @@ static D2 RescaleForNonVanishingEnds(D2 const &MM, double ZERO=1 return M; } +/*static D2 RescaleForNonVanishing(D2 const &MM, double ZERO=1.e-4){ + std::vector levels; + levels.push_back(-ZERO); + levels.push_back(ZERO); + //std::vector > mr = multi_roots(MM, levels); + }*/ + + //================================================================= //TODO: what's this for?!?! Piecewise > @@ -204,13 +203,15 @@ Geom::unitVector(D2 const &V_in, double tol, unsigned order){ D2 V = RescaleForNonVanishingEnds(V_in); if (V[0].empty() && V[1].empty()) return Piecewise >(D2(Linear(1),SBasis())); - SBasis x = V[0], y = V[1], a, b; + SBasis x = V[0], y = V[1]; SBasis r_eqn1, r_eqn2; Point v0 = unit_vector(V.at0()); Point v1 = unit_vector(V.at1()); - a.push_back(Linear(-v0[1],-v1[1])); - b.push_back(Linear( v0[0], v1[0])); + SBasis a(order+1, Linear()); + a[0] = Linear(-v0[1],-v1[1]); + SBasis b(order+1, Linear()); + b[0] = Linear( v0[0], v1[0]); r_eqn1 = -(a*x+b*y); r_eqn2 = Linear(1.)-(a*a+b*b); @@ -231,8 +232,8 @@ Geom::unitVector(D2 const &V_in, double tol, unsigned order){ a1 = r1/dot(v1,V(1))*v1[0]-rr1/2*v1[1]; b1 = r1/dot(v1,V(1))*v1[1]+rr1/2*v1[0]; - a.push_back(Linear(a0,a1)); - b.push_back(Linear(b0,b1)); + a[k] = Linear(a0,a1); + b[k] = Linear(b0,b1); //TODO: use "incremental" rather than explicit formulas. r_eqn1 = -(a*x+b*y); r_eqn2 = Linear(1)-(a*a+b*b); @@ -561,10 +562,10 @@ std::vector solve_lambda0(double a0,double a1,double c0,double c1, int insist_on_speeds_signs){ - SBasis p; - p.push_back(Linear( a1*c0*c0+c1, a1*a0*(a0+ 2*c0) +a1*c0*c0 +c1 -1 )); - p.push_back(Linear( -a1*a0*(a0+2*c0), -a1*a0*(3*a0+2*c0) )); - p.push_back(Linear( a1*a0*a0 )); + SBasis p(3, Linear()); + p[0] = Linear( a1*c0*c0+c1, a1*a0*(a0+ 2*c0) +a1*c0*c0 +c1 -1 ); + p[1] = Linear( -a1*a0*(a0+2*c0), -a1*a0*(3*a0+2*c0) ); + p[2] = Linear( a1*a0*a0 ); OptInterval domain = find_bounds_for_lambda0(a0,a1,c0,c1,insist_on_speeds_signs); if ( !domain ) @@ -671,9 +672,11 @@ Geom::cubics_fitting_curvature(Point const &M0, Point const &M1, Point V1 = lambda1[i]*dM1; D2 cubic; for(unsigned dim=0;dim<2;dim++){ - cubic[dim] = Linear(M0[dim],M1[dim]); - cubic[dim].push_back(Linear( M0[dim]-M1[dim]+V0[dim], - -M0[dim]+M1[dim]-V1[dim])); + SBasis c(2, Linear()); + c[0] = Linear(M0[dim],M1[dim]); + c[1] = Linear( M0[dim]-M1[dim]+V0[dim], + -M0[dim]+M1[dim]-V1[dim]); + cubic[dim] = c; } #if 0 Piecewise k = curvature(result); diff --git a/src/2geom/sbasis.cpp b/src/2geom/sbasis.cpp index 0bd672c15..bdc40c936 100644 --- a/src/2geom/sbasis.cpp +++ b/src/2geom/sbasis.cpp @@ -65,10 +65,10 @@ bool SBasis::isFinite() const { There is an elegant way to compute the value and n derivatives for a polynomial using a variant of horner's rule. Someone will someday work out how for sbasis. */ std::vector SBasis::valueAndDerivatives(double t, unsigned n) const { - std::vector ret(n); - ret.push_back(valueAt(t)); + std::vector ret(n+1); + ret[0] = valueAt(t); SBasis tmp = *this; - for(unsigned i = 0; i < n; i++) { + for(unsigned i = 1; i < n+1; i++) { tmp.derive(); ret[i+1] = tmp.valueAt(t); } @@ -82,18 +82,17 @@ std::vector SBasis::valueAndDerivatives(double t, unsigned n) const { */ SBasis operator+(const SBasis& a, const SBasis& b) { - SBasis result; const unsigned out_size = std::max(a.size(), b.size()); const unsigned min_size = std::min(a.size(), b.size()); - result.reserve(out_size); + SBasis result(out_size, Linear()); for(unsigned i = 0; i < min_size; i++) { - result.push_back(a[i] + b[i]); + result[i] = a[i] + b[i]; } for(unsigned i = min_size; i < a.size(); i++) - result.push_back(a[i]); + result[i] = a[i]; for(unsigned i = min_size; i < b.size(); i++) - result.push_back(b[i]); + result[i] = b[i]; assert(result.size() == out_size); return result; @@ -105,18 +104,17 @@ SBasis operator+(const SBasis& a, const SBasis& b) { */ SBasis operator-(const SBasis& a, const SBasis& b) { - SBasis result; const unsigned out_size = std::max(a.size(), b.size()); const unsigned min_size = std::min(a.size(), b.size()); - result.reserve(out_size); + SBasis result(out_size, Linear()); for(unsigned i = 0; i < min_size; i++) { - result.push_back(a[i] - b[i]); + result[i] = a[i] - b[i]; } for(unsigned i = min_size; i < a.size(); i++) - result.push_back(a[i]); + result[i] = a[i]; for(unsigned i = min_size; i < b.size(); i++) - result.push_back(-b[i]); + result[i] = -b[i]; assert(result.size() == out_size); return result; @@ -130,12 +128,12 @@ SBasis operator-(const SBasis& a, const SBasis& b) { SBasis& operator+=(SBasis& a, const SBasis& b) { const unsigned out_size = std::max(a.size(), b.size()); const unsigned min_size = std::min(a.size(), b.size()); - a.reserve(out_size); + a.resize(out_size); for(unsigned i = 0; i < min_size; i++) a[i] += b[i]; for(unsigned i = min_size; i < b.size(); i++) - a.push_back(b[i]); + a[i] = b[i]; assert(a.size() == out_size); return a; @@ -149,12 +147,12 @@ SBasis& operator+=(SBasis& a, const SBasis& b) { SBasis& operator-=(SBasis& a, const SBasis& b) { const unsigned out_size = std::max(a.size(), b.size()); const unsigned min_size = std::min(a.size(), b.size()); - a.reserve(out_size); + a.resize(out_size); for(unsigned i = 0; i < min_size; i++) a[i] -= b[i]; for(unsigned i = min_size; i < b.size(); i++) - a.push_back(-b[i]); + a[i] = -b[i]; assert(a.size() == out_size); return a; @@ -166,10 +164,9 @@ SBasis& operator-=(SBasis& a, const SBasis& b) { */ SBasis operator*(SBasis const &a, double k) { - SBasis c; - c.reserve(a.size()); + SBasis c(a.size(), Linear()); for(unsigned i = 0; i < a.size(); i++) - c.push_back(a[i] * k); + c[i] = a[i] * k; return c; } @@ -195,12 +192,14 @@ SBasis& operator*=(SBasis& a, double b) { */ SBasis shift(SBasis const &a, int sh) { - SBasis c = a; - if(sh > 0) { - c.insert(c.begin(), sh, Linear(0,0)); - } else { - //TODO: truncate - } + size_t n = a.size()+sh; + SBasis c(n, Linear()); + size_t m = std::max(0, sh); + + for(int i = 0; i < sh; i++) + c[i] = Linear(0,0); + for(size_t i = m, j = 0; i < n; i++, j++) + c[i] = a[j]; return c; } @@ -211,11 +210,13 @@ SBasis shift(SBasis const &a, int sh) { */ SBasis shift(Linear const &a, int sh) { - SBasis c; - if(sh >= 0) { - c.insert(c.begin(), sh, Linear(0,0)); - c.push_back(a); - } + size_t n = 1+sh; + SBasis c(n, Linear()); + + for(int i = 0; i < sh; i++) + c[i] = Linear(0,0); + if(sh >= 0) + c[sh] = a; return c; } @@ -230,8 +231,8 @@ SBasis multiply(SBasis const &a, SBasis const &b) { c.resize(a.size() + b.size(), Linear(0,0)); for(unsigned j = 0; j < b.size(); j++) { for(unsigned i = j; i < a.size()+j; i++) { - double tri = Tri(b[j])*Tri(a[i-j]); - c[i+1/*shift*/] += Linear(Hat(-tri)); + double tri = b[j].tri()*a[i-j].tri(); + c[i+1/*shift*/] += Linear(-tri); } } for(unsigned j = 0; j < b.size(); j++) { @@ -258,8 +259,8 @@ SBasis multiply_add(SBasis const &a, SBasis const &b, SBasis c) { c.resize(a.size() + b.size(), Linear(0,0)); for(unsigned j = 0; j < b.size(); j++) { for(unsigned i = j; i < a.size()+j; i++) { - double tri = Tri(b[j])*Tri(a[i-j]); - c[i+1/*shift*/] += Linear(Hat(-tri)); + double tri = b[j].tri()*a[i-j].tri(); + c[i+1/*shift*/] += Linear(-tri); } } for(unsigned j = 0; j < b.size(); j++) { @@ -296,12 +297,12 @@ SBasis integral(SBasis const &c) { a[0] = Linear(0,0); for(unsigned k = 1; k < c.size() + 1; k++) { - double ahat = -Tri(c[k-1])/(2*k); - a[k] = Hat(ahat); + double ahat = -c[k-1].tri()/(2*k); + a[k][0] = a[k][1] = ahat; } double aTri = 0; for(int k = c.size()-1; k >= 0; k--) { - aTri = (Hat(c[k]).d + (k+1)*aTri/2)/(2*k+1); + aTri = (c[k].hat() + (k+1)*aTri/2)/(2*k+1); a[k][0] -= aTri/2; a[k][1] += aTri/2; } @@ -397,7 +398,7 @@ SBasis reciprocal(Linear const &a, int k) { SBasis c; assert(!a.isZero()); c.resize(k, Linear(0,0)); - double r_s0 = (Tri(a)*Tri(a))/(-a[0]*a[1]); + double r_s0 = (a.tri()*a.tri())/(-a[0]*a[1]); double r_s0k = 1; for(unsigned i = 0; i < (unsigned)k; i++) { c[i] = Linear(r_s0k/a[0], r_s0k/a[1]); @@ -444,7 +445,7 @@ SBasis compose(SBasis const &a, SBasis const &b) { SBasis r; for(int i = a.size()-1; i >= 0; i--) { - r = multiply_add(r, s, SBasis(Linear(Hat(a[i][0]))) - b*a[i][0] + b*a[i][1]); + r = multiply_add(r, s, SBasis(Linear(a[i][0])) - b*a[i][0] + b*a[i][1]); } return r; } @@ -460,7 +461,7 @@ SBasis compose(SBasis const &a, SBasis const &b, unsigned k) { SBasis r; for(int i = a.size()-1; i >= 0; i--) { - r = multiply_add(r, s, SBasis(Linear(Hat(a[i][0]))) - b*a[i][0] + b*a[i][1]); + r = multiply_add(r, s, SBasis(Linear(a[i][0])) - b*a[i][0] + b*a[i][1]); } r.truncate(k); return r; @@ -499,11 +500,11 @@ SBasis inverse(SBasis a, int k) { if(a1 != 1) { a /= a1; } - SBasis c; // c(v) := 0 + SBasis c(k, Linear()); // c(v) := 0 if(a.size() >= 2 && k == 2) { - c.push_back(Linear(0,1)); + c[0] = Linear(0,1); Linear t1(1+a[1][0], 1-a[1][1]); // t_1 - c.push_back(Linear(-a[1][0]/t1[0], -a[1][1]/t1[1])); + c[1] = Linear(-a[1][0]/t1[0], -a[1][1]/t1[1]); } else if(a.size() >= 2) { // non linear SBasis r = Linear(0,1); // r(u) := r_0(u) := u Linear t1(1./(1+a[1][0]), 1./(1-a[1][1])); // 1./t_1 @@ -519,7 +520,7 @@ SBasis inverse(SBasis a, int k) { //assert(t1 == t[1]); #endif - c.resize(k+1, Linear(0,0)); + //c.resize(k+1, Linear(0,0)); for(unsigned i = 0; i < (unsigned)k; i++) { // for i:=0 to k do #ifdef DEBUG_INVERSION std::cout << "-------" << i << ": ---------" < #include <2geom/exception.h> +//#define USE_SBASISN 1 -#ifdef USE_SBASIS_OF + +#if defined(USE_SBASIS_OF) #include "sbasis-of.h" +#elif defined(USE_SBASISN) + +#include "sbasisN.h" +namespace Geom{ + +/*** An empty SBasis is identically 0. */ +class SBasis : public SBasisN<1>; + +}; #else namespace Geom{ /*** An empty SBasis is identically 0. */ -class SBasis : public std::vector{ +class SBasis{ + std::vector d; + void push_back(Linear const&l) { d.push_back(l); } + public: + size_t size() const {return d.size();} + Linear operator[](unsigned i) const { + assert(i < size()); + return d[i]; + } + Linear& operator[](unsigned i) { return d.at(i); } + Linear const* begin() const { return (Linear const*)&*d.begin();} + Linear const* end() const { return (Linear const*)&*d.end();} + Linear* begin() { return (Linear*)&*d.begin();} + Linear* end() { return (Linear*)&*d.end();} + bool empty() const {return d.empty();} + Linear &back() {return d.back();} + Linear const &back() const {return d.back();} + void pop_back() { d.pop_back();} + void resize(unsigned n) { d.resize(n);} + void resize(unsigned n, Linear const& l) { d.resize(n, l);} + void reserve(unsigned n) { d.reserve(n);} + void clear() {d.clear();} + void insert(Linear* before, const Linear* src_begin, const Linear* src_end) { d.insert(std::vector::iterator(before), src_begin, src_end);} + //void insert(Linear* aa, Linear* bb, Linear* cc} { d.insert(aa, bb, cc);} + Linear& at(unsigned i) { return d.at(i);} + //void insert(Linear* before, int& n, Linear const &l) { d.insert(std::vector::iterator(before), n, l);} + + SBasis() {} explicit SBasis(double a) { push_back(Linear(a,a)); } SBasis(SBasis const & a) : - std::vector(a) + d(a.d) {} SBasis(Linear const & bo) { push_back(bo); @@ -68,6 +106,7 @@ public: SBasis(Linear* bo) { push_back(*bo); } + explicit SBasis(size_t n, Linear const&l) : d(n, l) {} //IMPL: FragmentConcept typedef double output_type; @@ -93,6 +132,8 @@ public: inline double at1() const{ if(empty()) return 0; else return (*this)[0][1]; } + + int degreesOfFreedom() const { return size()*2;} double valueAt(double t) const { double s = t*(1-t); @@ -119,14 +160,7 @@ public: // compute f(g) SBasis operator()(SBasis const & g) const; - Linear operator[](unsigned i) const { - assert(i < size()); - return std::vector::operator[](i); - } - //MUTATOR PRISON - Linear& operator[](unsigned i) { return this->at(i); } - //remove extra zeros void normalize() { while(!empty() && 0 == back()[0] && 0 == back()[1]) @@ -152,21 +186,20 @@ OptInterval bounds_local(SBasis const &a, const OptInterval &t, int order = 0); useful for reversing a parameteric curve. */ inline SBasis reverse(SBasis const &a) { - SBasis result; - result.reserve(a.size()); + SBasis result(a.size(), Linear()); + for(unsigned k = 0; k < a.size(); k++) - result.push_back(reverse(a[k])); + result[k] = reverse(a[k]); return result; } //IMPL: ScalableConcept inline SBasis operator-(const SBasis& p) { if(p.isZero()) return SBasis(); - SBasis result; - result.reserve(p.size()); + SBasis result(p.size(), Linear()); for(unsigned i = 0; i < p.size(); i++) { - result.push_back(-p[i]); + result[i] = -p[i]; } return result; } @@ -183,7 +216,7 @@ SBasis& operator+=(SBasis& a, const SBasis& b); SBasis& operator-=(SBasis& a, const SBasis& b); //TODO: remove? -inline SBasis operator+(const SBasis & a, Linear const & b) { +/*inline SBasis operator+(const SBasis & a, Linear const & b) { if(b.isZero()) return a; if(a.isZero()) return b; SBasis result(a); @@ -209,7 +242,7 @@ inline SBasis& operator-=(SBasis& a, const Linear& b) { else a[0] -= b; return a; -} + }*/ //IMPL: OffsetableConcept inline SBasis operator+(const SBasis & a, double b) { @@ -226,14 +259,14 @@ inline SBasis operator-(const SBasis & a, double b) { } inline SBasis& operator+=(SBasis& a, double b) { if(a.isZero()) - a.push_back(Linear(b,b)); + a = SBasis(Linear(b,b)); else a[0] += b; return a; } inline SBasis& operator-=(SBasis& a, double b) { if(a.isZero()) - a.push_back(Linear(-b,-b)); + a = SBasis(Linear(-b,-b)); else a[0] -= b; return a; @@ -300,6 +333,7 @@ SBasis compose_inverse(SBasis const &f, SBasis const &g, unsigned order=2, doubl */ inline SBasis portion(const SBasis &t, double from, double to) { return compose(t, Linear(from, to)); } +inline SBasis portion(const SBasis &t, Interval ivl) { return compose(t, Linear(ivl[0], ivl[1])); } // compute f(g) inline SBasis diff --git a/src/2geom/svg-elliptical-arc.h b/src/2geom/svg-elliptical-arc.h index 92ec51b49..dad9000c1 100644 --- a/src/2geom/svg-elliptical-arc.h +++ b/src/2geom/svg-elliptical-arc.h @@ -248,6 +248,8 @@ class SVGEllipticalArc : public Curve else return SBasisCurve(toSBasis()).winding(p); } + + int degreesOfFreedom() const { return 5;} Curve *derivative() const; diff --git a/src/2geom/sweep.cpp b/src/2geom/sweep.cpp index 53b953255..227c47e09 100644 --- a/src/2geom/sweep.cpp +++ b/src/2geom/sweep.cpp @@ -4,6 +4,15 @@ namespace Geom { +/** + * \brief Make a list of pairs of self intersections in a list of Rects. + * + * \param rs: vector of Rect. + * + * [(A = rs[i], B = rs[j]) for i,J in enumerate(pairs) for j in J] + * then A.left <= B.left + */ + std::vector > sweep_bounds(std::vector rs) { std::vector events; events.reserve(rs.size()*2); std::vector > pairs(rs.size()); @@ -34,6 +43,15 @@ std::vector > sweep_bounds(std::vector rs) { return pairs; } +/** + * \brief Make a list of pairs of red-blue intersections between two lists of Rects. + * + * \param a: vector of Rect. + * \param b: vector of Rect. + * + * [(A = rs[i], B = rs[j]) for i,J in enumerate(pairs) for j in J] + * then A.left <= B.left, A in a, B in b + */ std::vector > sweep_bounds(std::vector a, std::vector b) { std::vector > pairs(a.size()); if(a.empty() || b.empty()) return pairs; -- 2.30.2