From: johanengelen Date: Mon, 1 Sep 2008 19:29:30 +0000 (+0000) Subject: update 2geom (rev. 1569) X-Git-Url: https://git.tokkee.org/?a=commitdiff_plain;h=6bc0b25077dcb0cce5dea357de5bab735babe891;p=inkscape.git update 2geom (rev. 1569) --- diff --git a/src/2geom/basic-intersection.cpp b/src/2geom/basic-intersection.cpp index 48b990445..334c23fea 100644 --- a/src/2geom/basic-intersection.cpp +++ b/src/2geom/basic-intersection.cpp @@ -1,8 +1,14 @@ + + + #include <2geom/basic-intersection.h> +#include <2geom/sbasis-to-bezier.h> #include <2geom/exception.h> + + #include #include - + unsigned intersect_steps = 0; @@ -16,10 +22,10 @@ public: } void split(double t, OldBezier &a, OldBezier &b) const; Point operator()(double t) const; - + ~OldBezier() {} - void bounds(double &minax, double &maxax, + 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 @@ -45,17 +51,17 @@ public: } } - + }; static std::vector > find_intersections( OldBezier a, OldBezier b); -static std::vector > +static std::vector > find_self_intersections(OldBezier const &Sb, D2 const & A); std::vector > -find_intersections( vector const & A, +find_intersections( vector const & A, vector const & B) { OldBezier a, b; a.p = A; @@ -63,23 +69,23 @@ find_intersections( vector const & A, return find_intersections(a,b); } -std::vector > +std::vector > find_self_intersections(OldBezier const &/*Sb*/) { THROW_NOTIMPLEMENTED(); } -std::vector > +std::vector > find_self_intersections(D2 const & A) { OldBezier Sb; - Sb.p = sbasis_to_bezier(A); + sbasis_to_bezier(Sb.p, A); return find_self_intersections(Sb, A); } -static std::vector > +static std::vector > find_self_intersections(OldBezier const &Sb, D2 const & A) { - + vector dr = roots(derivative(A[X])); { vector dyr = roots(derivative(A[Y])); @@ -90,9 +96,9 @@ find_self_intersections(OldBezier const &Sb, D2 const & A) { // We want to be sure that we have no empty segments sort(dr.begin(), dr.end()); unique(dr.begin(), dr.end()); - + std::vector > all_si; - + vector pieces; { OldBezier in = Sb, l, r; @@ -104,7 +110,7 @@ find_self_intersections(OldBezier const &Sb, D2 const & A) { } for(unsigned i = 0; i < dr.size()-1; i++) { for(unsigned j = i+1; j < dr.size()-1; j++) { - std::vector > section = + std::vector > section = find_intersections( pieces[i], pieces[j]); for(unsigned k = 0; k < section.size(); k++) { double l = section[k].first; @@ -118,17 +124,17 @@ find_self_intersections(OldBezier const &Sb, D2 const & A) { } } } - + // 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 @@ -142,12 +148,12 @@ void OldBezier::split(double t, OldBezier &left, OldBezier &right) const { std::copy(p.begin(), p.end(), Vtemp[0]); /* Triangle computation */ - for (unsigned i = 1; i < sz; i++) { + 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++) @@ -156,6 +162,7 @@ void OldBezier::split(double t, OldBezier &left, OldBezier &right) const { 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 @@ -169,13 +176,35 @@ Point OldBezier::operator()(double t) const { std::copy(p.begin(), p.end(), Vtemp[0]); /* Triangle computation */ - for (unsigned i = 1; i < sz; i++) { + 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 maxbx ) || ( minay > maxby ) || ( minbx > maxax ) || ( minby > maxay ) ); } - -/* - * Recursively intersect two curves keeping track of their real parameters + +/* + * 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 + * 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 + * 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 + * 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 + * 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 + * 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 @@ -335,7 +364,7 @@ void recursively_intersect( OldBezier a, double t0, double t1, int deptha, } 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. @@ -354,7 +383,7 @@ unsigned wangs_theorem(OldBezier a) { 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 ) + 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 ) ); @@ -367,63 +396,109 @@ struct rparams OldBezier &A; OldBezier &B; }; - + static int intersect_polish_f (const gsl_vector * x, void *params, gsl_vector * f) { const double x0 = gsl_vector_get (x, 0); const double x1 = gsl_vector_get (x, 1); - - Geom::Point dx = ((struct rparams *) params)->A(x0) - + + Geom::Point dx = ((struct rparams *) params)->A(x0) - ((struct rparams *) params)->B(x1); - + gsl_vector_set (f, 0, dx[0]); gsl_vector_set (f, 1, dx[1]); - + return GSL_SUCCESS; } +typedef union dbl_64{ + long long i64; + double d64; +}; + +static double EpsilonBy(double value, int eps) +{ + dbl_64 s; + s.d64 = value; + s.i64 += eps; + return s.d64; +} + + static void intersect_polish_root (OldBezier &A, double &s, OldBezier &B, double &t) { const gsl_multiroot_fsolver_type *T; gsl_multiroot_fsolver *sol; - + int status; size_t iter = 0; - + const size_t n = 2; struct rparams p = {A, B}; gsl_multiroot_function f = {&intersect_polish_f, n, &p}; - + double x_init[2] = {s, t}; gsl_vector *x = gsl_vector_alloc (n); - + gsl_vector_set (x, 0, x_init[0]); gsl_vector_set (x, 1, x_init[1]); - + T = gsl_multiroot_fsolver_hybrids; sol = gsl_multiroot_fsolver_alloc (T, 2); gsl_multiroot_fsolver_set (sol, &f, x); - + do { iter++; status = gsl_multiroot_fsolver_iterate (sol); - + if (status) /* check if solver is stuck */ break; - + status = gsl_multiroot_test_residual (sol->f, 1e-12); } while (status == GSL_CONTINUE && iter < 1000); - + s = gsl_vector_get (sol->x, 0); t = gsl_vector_get (sol->x, 1); - + gsl_multiroot_fsolver_free (sol); gsl_vector_free (x); + + // This code does a neighbourhood search for minor improvements. + double best_v = L1(A(s) - B(t)); + //std::cout << "------\n" << best_v << std::endl; + Point best(s,t); + while (true) { + Point trial = best; + double trial_v = best_v; + for(int nsi = -1; nsi < 2; nsi++) { + for(int nti = -1; nti < 2; nti++) { + Point n(EpsilonBy(best[0], nsi), + EpsilonBy(best[1], nti)); + double c = L1(A(n[0]) - B(n[1])); + //std::cout << c << "; "; + if (c < trial_v) { + trial = n; + trial_v = c; + } + } + } + if(trial == best) { + //std::cout << "\n" << s << " -> " << s - best[0] << std::endl; + //std::cout << t << " -> " << t - best[1] << std::endl; + //std::cout << best_v << std::endl; + s = best[0]; + t = best[1]; + return; + } else { + best = trial; + best_v = trial_v; + } + } } @@ -432,8 +507,8 @@ std::vector > find_intersections( OldBezier a, OldBezi std::vector > parameters; if( intersect_BB( a, b ) ) { - recursively_intersect( a, 0., 1., wangs_theorem(a), - b, 0., 1., wangs_theorem(b), + recursively_intersect( a, 0., 1., wangs_theorem(a), + b, 0., 1., wangs_theorem(b), parameters); } for(unsigned i = 0; i < parameters.size(); i++) diff --git a/src/2geom/basic-intersection.h b/src/2geom/basic-intersection.h index 0090b0305..d464891f9 100644 --- a/src/2geom/basic-intersection.h +++ b/src/2geom/basic-intersection.h @@ -1,25 +1,47 @@ + + +#include <2geom/point.h> #include <2geom/sbasis.h> -#include <2geom/bezier-to-sbasis.h> -#include <2geom/sbasis-to-bezier.h> #include <2geom/d2.h> +#include +#include + + namespace Geom { -std::vector > -find_intersections( D2 const & A, +std::vector > +find_intersections( D2 const & A, D2 const & B); -std::vector > +std::vector > find_self_intersections(D2 const & A); // Bezier form -std::vector > -find_intersections( std::vector const & A, +std::vector > +find_intersections( std::vector const & A, std::vector const & B); -std::vector > +std::vector > find_self_intersections(std::vector const & A); + +/* + * find_intersection + * + * input: A, B - set of control points of two Bezier curve + * input: precision - required precision of computation + * output: xs - set of pairs of parameter values + * at which crossing happens + * + * This routine is based on the Bezier Clipping Algorithm, + * see: Sederberg - Computer Aided Geometric Design + */ +void find_intersections (std::vector< std::pair > & xs, + std::vector const& A, + std::vector const& B, + double precision = 1e-5); + }; /* diff --git a/src/2geom/bezier-to-sbasis.h b/src/2geom/bezier-to-sbasis.h index 8b421f2e7..d5279a570 100644 --- a/src/2geom/bezier-to-sbasis.h +++ b/src/2geom/bezier-to-sbasis.h @@ -32,12 +32,14 @@ #define _BEZIER_TO_SBASIS #include <2geom/coord.h> - -#include <2geom/d2.h> #include <2geom/point.h> +#include <2geom/d2.h> +#include <2geom/sbasis-to-bezier.h> -namespace Geom{ +namespace Geom +{ +#if 0 inline SBasis bezier_to_sbasis(Coord const *handles, unsigned order) { if(order == 0) return Linear(handles[0]); @@ -50,7 +52,8 @@ inline SBasis bezier_to_sbasis(Coord const *handles, unsigned order) { template -inline D2 handles_to_sbasis(T const &handles, unsigned order) { +inline D2 handles_to_sbasis(T const &handles, unsigned order) +{ double v[2][order+1]; for(unsigned i = 0; i <= order; i++) for(unsigned j = 0; j < 2; j++) @@ -58,8 +61,29 @@ inline D2 handles_to_sbasis(T const &handles, unsigned order) { return D2(bezier_to_sbasis(v[0], order), bezier_to_sbasis(v[1], order)); } +#endif + + +template +inline +D2 handles_to_sbasis(T const& handles, unsigned order) +{ + D2 sbc; + size_t sz = order + 1; + std::vector v; + v.reserve(sz); + for (size_t i = 0; i < sz; ++i) + v.push_back(handles[i]); + bezier_to_sbasis(sbc, v); + return sbc; +} + + +} // end namespace Geom + + + -}; #endif /* Local Variables: diff --git a/src/2geom/bezier.h b/src/2geom/bezier.h index 79b15cd03..6e11ad58b 100644 --- a/src/2geom/bezier.h +++ b/src/2geom/bezier.h @@ -36,7 +36,6 @@ #include <2geom/coord.h> #include #include <2geom/isnan.h> -#include <2geom/bezier-to-sbasis.h> #include <2geom/d2.h> #include <2geom/solver.h> #include @@ -50,7 +49,7 @@ inline Coord subdivideArr(Coord t, Coord const *v, Coord *left, Coord *right, un std::copy(v, v+order+1, vtemp[0]); /* Triangle computation */ - for (unsigned i = 1; i <= order; i++) { + for (unsigned i = 1; i <= order; i++) { for (unsigned j = 0; j <= order - i; j++) { vtemp[i][j] = lerp(t, vtemp[i-1][j], vtemp[i-1][j+1]); } @@ -84,8 +83,8 @@ protected: public: unsigned int order() const { return c_.size()-1;} unsigned int size() const { return c_.size();} - - Bezier() :c_(0., 32) {} + + Bezier() {} Bezier(const Bezier& b) :c_(b.c_) {} Bezier &operator=(Bezier const &other) { if ( c_.size() != other.c_.size() ) { @@ -126,11 +125,21 @@ public: c_[0] = c0; c_[1] = c1; c_[2] = c2; c_[3] = c3; } + void resize (unsigned int n, Coord v = 0) + { + c_.resize (n, v); + } + + void clear() + { + c_.resize(0); + } + inline unsigned degree() const { return order(); } //IMPL: FragmentConcept typedef Coord output_type; - inline bool isZero() const { + inline bool isZero() const { for(unsigned i = 0; i <= order(); i++) { if(c_[i] != 0) return false; } @@ -151,23 +160,27 @@ public: inline Coord at0() const { return c_[0]; } inline Coord at1() const { return c_[order()]; } - inline Coord valueAt(double t) const { - return subdivideArr(t, &c_[0], NULL, NULL, order()); + inline Coord valueAt(double t) const { + return subdivideArr(t, &c_[0], NULL, NULL, order()); } inline Coord operator()(double t) const { return valueAt(t); } - inline SBasis toSBasis() const { - return bezier_to_sbasis(&c_[0], order()); - } + SBasis toSBasis() const; +// inline SBasis toSBasis() const { +// SBasis sb; +// bezier_to_sbasis(sb, (*this)); +// return sb; +// //return bezier_to_sbasis(&c_[0], order()); +// } //Only mutator inline Coord &operator[](unsigned ix) { return c_[ix]; } inline Coord const &operator[](unsigned ix) const { return c_[ix]; } inline void setPoint(unsigned ix, double val) { c_[ix] = val; } - + /* This is inelegant, as it uses several extra stores. I think there might be a way to * evaluate roughly in situ. */ - + std::vector valueAndDerivatives(Coord t, unsigned n_derivs) const { std::vector val_n_der; Coord d_[order()+1]; @@ -186,7 +199,7 @@ public: val_n_der.resize(n_derivs); return val_n_der; } - + std::pair subdivide(Coord t) const { Bezier a(Bezier::Order(*this)), b(Bezier::Order(*this)); subdivideArr(t, &c_[0], &a.c_[0], &b.c_[0], order()); @@ -200,6 +213,18 @@ public: } }; + +void bezier_to_sbasis (SBasis & sb, Bezier const& bz); + + +inline +SBasis Bezier::toSBasis() const { + SBasis sb; + bezier_to_sbasis(sb, (*this)); + return sb; + //return bezier_to_sbasis(&c_[0], order()); +} + //TODO: implement others inline Bezier operator+(const Bezier & a, double v) { Bezier result = Bezier(Bezier::Order(a)); @@ -266,7 +291,7 @@ inline Bezier derivative(const Bezier & a) { //if(a.order() == 1) return Bezier(0.0); if(a.order() == 1) return Bezier(a.c_[1]-a.c_[0]); Bezier der(Bezier::Order(a.order()-1)); - + for(unsigned i = 0; i < a.order(); i++) { der.c_[i] = a.order()*(a.c_[i+1] - a.c_[i]); } @@ -275,7 +300,7 @@ inline Bezier derivative(const Bezier & a) { inline Bezier integral(const Bezier & a) { Bezier inte(Bezier::Order(a.order()+1)); - + inte[0] = 0; for(unsigned i = 0; i < inte.order(); i++) { inte[i+1] = inte[i] + a[i]/(inte.order()); diff --git a/src/2geom/crossing.cpp b/src/2geom/crossing.cpp index a4d4f6ab1..d717a4ed5 100644 --- a/src/2geom/crossing.cpp +++ b/src/2geom/crossing.cpp @@ -27,7 +27,7 @@ double wrap_dist(double from, double to, double size, bool rev) { } } } - +/* CrossingGraph create_crossing_graph(std::vector const &p, Crossings const &crs) { std::vector locs; CrossingGraph ret; @@ -75,7 +75,7 @@ CrossingGraph create_crossing_graph(std::vector const &p, Crossings const } return ret; - + */ /* Various incoherent code bits // list of sets of edges, each set corresponding to those emanating from the path CrossingGraph ret; @@ -111,7 +111,7 @@ CrossingGraph create_crossing_graph(std::vector const &p, Crossings const } } */ -} +//} void merge_crossings(Crossings &a, Crossings &b, unsigned i) { Crossings n; diff --git a/src/2geom/crossing.h b/src/2geom/crossing.h index 546e33ebd..b16c7e46d 100644 --- a/src/2geom/crossing.h +++ b/src/2geom/crossing.h @@ -25,6 +25,7 @@ struct Crossing { }; +/* struct Edge { unsigned node, path; double time; @@ -49,7 +50,6 @@ struct CrossingNode { } }; -typedef std::vector Crossings; typedef std::vector CrossingGraph; @@ -61,6 +61,7 @@ struct TimeOrder { class Path; CrossingGraph create_crossing_graph(std::vector const &p, Crossings const &crs); +*/ /*inline bool are_near(Crossing a, Crossing b) { return are_near(a.ta, b.ta) && are_near(a.tb, b.tb); @@ -78,6 +79,7 @@ struct CrossingOrder { } }; +typedef std::vector Crossings; typedef std::vector CrossingSet; diff --git a/src/2geom/forward.h b/src/2geom/forward.h index 9bd5dedbe..d4ac1af1e 100644 --- a/src/2geom/forward.h +++ b/src/2geom/forward.h @@ -3,7 +3,7 @@ * * Authors: * Johan Engelen - * + * * Copyright 2008 authors * * This library is free software; you can redistribute it and/or @@ -37,6 +37,7 @@ namespace Geom { +class Bezier; template class BezierCurve; template<> class BezierCurve<0>; typedef BezierCurve<2> QuadraticBezier; diff --git a/src/2geom/interval.h b/src/2geom/interval.h index 74e03aec9..2d9f4abd8 100644 --- a/src/2geom/interval.h +++ b/src/2geom/interval.h @@ -84,8 +84,8 @@ public: return contains(val._b[0]) || contains(val._b[1]) || val.contains(*this); } - inline bool operator==(Interval other) { return _b[0] == other._b[0] && _b[1] == other._b[1]; } - inline bool operator!=(Interval other) { return _b[0] != other._b[0] || _b[1] != other._b[1]; } + inline bool operator==(Interval other) const { return _b[0] == other._b[0] && _b[1] == other._b[1]; } + inline bool operator!=(Interval other) const { return _b[0] != other._b[0] || _b[1] != other._b[1]; } //IMPL: OffsetableConcept //TODO: rename output_type to something else in the concept diff --git a/src/2geom/numeric/fitting-model.h b/src/2geom/numeric/fitting-model.h index de8f24760..564663cf7 100644 --- a/src/2geom/numeric/fitting-model.h +++ b/src/2geom/numeric/fitting-model.h @@ -47,6 +47,30 @@ namespace Geom { namespace NL { +/* + * A model is an abstraction for an expression dependent from a parameter where + * the coefficients of this expression are the unknowns of the fitting problem. + * For a ceratain number of parameter values we know the related values + * the expression evaluates to: from each parameter value we get a row of + * the matrix of the fitting problem, from each expression value we get + * the related constant term. + * Example: given the model a*x^2 + b*x + c = 0; from x = 1 we get + * the equation a + b + c = 0, in this example the constant term is always + * the same for each parameter value. + * + * A model is required to implement 3 methods: + * + * - size : returns the number of unknown coefficients that appear in + * the expression of the fitting problem; + * - feed : its input is a parameter value and the related expression value, + * it generates a matrix row and a new entry of the constant vector + * of the fitting problem; + * - instance : it has an input parameter represented by the raw vector + * solution of the fitting problem and an output parameter + * of type InstanceType that return a specific object that is + * generated using the fitting problem solution, in the example + * above the object could be a Poly type. + */ /* * completely unknown models must inherit from this template class; diff --git a/src/2geom/numeric/fitting-tool.h b/src/2geom/numeric/fitting-tool.h index d589d86e3..47c588cb2 100644 --- a/src/2geom/numeric/fitting-tool.h +++ b/src/2geom/numeric/fitting-tool.h @@ -43,11 +43,20 @@ #include +/* + * The least_square_fitter class represents a tool for solving a fitting + * problem with respect to a given model that represents an expression + * dependent from a parameter where the coefficients of this expression + * are the unknowns of the fitting problem. + * The minimizing solution is found by computing the pseudo-inverse + * of the problem matrix + */ + + namespace Geom { namespace NL { namespace detail { - template< typename ModelT> class lsf_base { diff --git a/src/2geom/path.h b/src/2geom/path.h index 9dfc92975..04b77862c 100644 --- a/src/2geom/path.h +++ b/src/2geom/path.h @@ -254,7 +254,7 @@ public: size_type max_size() const { return get_curves().max_size()-1; } - bool empty() const { return get_curves().size() == 1; } + bool empty() const { return (get_curves().size() == 1); } bool closed() const { return closed_; } void close(bool closed=true) { closed_ = closed; } diff --git a/src/2geom/rect.h b/src/2geom/rect.h index fa79201e8..21cdcf849 100644 --- a/src/2geom/rect.h +++ b/src/2geom/rect.h @@ -103,7 +103,7 @@ class D2 { inline double maxExtent() const { return std::max(f[X].extent(), f[Y].extent()); } inline bool isEmpty() const { - return f[X].isEmpty() && f[Y].isEmpty(); + return f[X].isEmpty() || f[Y].isEmpty(); } inline bool intersects(Rect const &r) const { return f[X].intersects(r[X]) && f[Y].intersects(r[Y]); diff --git a/src/2geom/sbasis-roots.cpp b/src/2geom/sbasis-roots.cpp index 861baf76d..09a84050b 100644 --- a/src/2geom/sbasis-roots.cpp +++ b/src/2geom/sbasis-roots.cpp @@ -8,7 +8,7 @@ * multi-roots using bernstein method, one approach would be: sort c take median and find roots of that - whenever a segment lies entirely on one side of the median, + whenever a segment lies entirely on one side of the median, find the median of the half and recurse. in essence we are implementing quicksort on a continuous function @@ -126,7 +126,7 @@ Interval bounds_local(const SBasis &sb, const Interval &i, int order) { going from f(a) up to C with slope M takes at least time (C-f(a))/M From this we conclude there are no roots before a'=a+min((f(a)-c)/m,(C-f(a))/M). Do the same for b: compute some b' such that there are no roots in (b',b]. - -if [a',b'] is not empty, repeat the process with [a',(a'+b')/2] and [(a'+b')/2,b']. + -if [a',b'] is not empty, repeat the process with [a',(a'+b')/2] and [(a'+b')/2,b']. unfortunately, extra care is needed about rounding errors, and also to avoid the repetition of roots, making things tricky and unpleasant... */ @@ -147,7 +147,7 @@ static void multi_roots_internal(SBasis const &f, double fa, double b, double fb){ - + if (f.size()==0){ int idx; idx=upper_level(levels,0,vtol); @@ -157,7 +157,7 @@ static void multi_roots_internal(SBasis const &f, } return; } -////usefull? +////usefull? // if (f.size()==1){ // int idxa=upper_level(levels,fa); // int idxb=upper_level(levels,fb); @@ -188,7 +188,7 @@ static void multi_roots_internal(SBasis const &f, } return; } - + int idxa=upper_level(levels,fa,vtol); int idxb=upper_level(levels,fb,vtol); @@ -219,9 +219,9 @@ static void multi_roots_internal(SBasis const &f, if (bs.max()>0 && idxb>0) tb_lo=b+(levels.at(idxb-1)-fb)/bs.max(); } - + double t0,t1; - t0=std::min(ta_hi,ta_lo); + t0=std::min(ta_hi,ta_lo); t1=std::max(tb_hi,tb_lo); //hum, rounding errors frighten me! so I add this +tol... if (t0>t1+htol) return;//no root here. @@ -256,7 +256,7 @@ std::vector > multi_roots(SBasis const &f, std::vector > roots(levels.size(), std::vector()); SBasis df=derivative(f); - multi_roots_internal(f,df,levels,roots,htol,vtol,a,f(a),b,f(b)); + multi_roots_internal(f,df,levels,roots,htol,vtol,a,f(a),b,f(b)); return(roots); } @@ -283,9 +283,9 @@ double Laguerre_internal(SBasis const & p, b = xk*b + p[j]; err = fabs(b) + abx*err; } - + err *= 1e-7; // magic epsilon for convergence, should be computed from tol - + double px = b; if(fabs(b) < err) return xk; @@ -293,7 +293,7 @@ double Laguerre_internal(SBasis const & p, // return xk; double G = d / px; double H = G*G - f / px; - + //std::cout << "G = " << G << "H = " << H; double radicand = (n - 1)*(n*H-G*G); //assert(radicand.real() > 0); @@ -315,7 +315,7 @@ double Laguerre_internal(SBasis const & p, #endif void subdiv_sbasis(SBasis const & s, - std::vector & roots, + std::vector & roots, double left, double right) { Interval bs = bounds_fast(s); if(bs.min() > 0 || bs.max() < 0) @@ -345,12 +345,16 @@ std::vector roots1(SBasis const & s) { std::vector roots(SBasis const & s) { switch(s.size()) { - case 0: - return std::vector(); - case 1: - return roots1(s); - default: - return sbasis_to_bezier(s).roots(); + case 0: + return std::vector(); + case 1: + return roots1(s); + default: + { + Bezier bz; + sbasis_to_bezier(bz, s); + return bz.roots(); + } } } diff --git a/src/2geom/sbasis-to-bezier.cpp b/src/2geom/sbasis-to-bezier.cpp index cbddccda8..27e3047fd 100644 --- a/src/2geom/sbasis-to-bezier.cpp +++ b/src/2geom/sbasis-to-bezier.cpp @@ -1,3 +1,260 @@ +/* + * Symmetric Power Basis - Bernstein Basis conversion routines + * + * Authors: + * Marco Cecchetti + * Nathan Hurst + * + * Copyright 2007-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 + * License version 2.1 as published by the Free Software Foundation + * (the "LGPL") or, at your option, under the terms of the Mozilla + * Public License Version 1.1 (the "MPL"). If you do not alter this + * notice, a recipient may use your version of this file under either + * the MPL or the LGPL. + * + * You should have received a copy of the LGPL along with this library + * in the file COPYING-LGPL-2.1; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * You should have received a copy of the MPL along with this library + * in the file COPYING-MPL-1.1 + * + * The contents of this file are subject to the Mozilla Public License + * Version 1.1 (the "License"); you may not use this file except in + * compliance with the License. You may obtain a copy of the License at + * http://www.mozilla.org/MPL/ + * + * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY + * OF ANY KIND, either express or implied. See the LGPL or the MPL for + * the specific language governing rights and limitations. + */ + + +#include <2geom/sbasis-to-bezier.h> +#include <2geom/choose.h> +#include <2geom/svg-path.h> +#include <2geom/exception.h> + +#include + + + + +namespace Geom +{ + +/* + * Symmetric Power Basis - Bernstein Basis conversion routines + * + * some remark about precision: + * interval [0,1], subdivisions: 10^3 + * - bezier_to_sbasis : up to degree ~72 precision is at least 10^-5 + * up to degree ~87 precision is at least 10^-3 + * - sbasis_to_bezier : up to order ~63 precision is at least 10^-15 + * precision is at least 10^-14 even beyond order 200 + * + * interval [-1,1], subdivisions: 10^3 + * - bezier_to_sbasis : up to degree ~21 precision is at least 10^-5 + * up to degree ~24 precision is at least 10^-3 + * - sbasis_to_bezier : up to order ~11 precision is at least 10^-5 + * up to order ~13 precision is at least 10^-3 + * + * interval [-10,10], subdivisions: 10^3 + * - bezier_to_sbasis : up to degree ~7 precision is at least 10^-5 + * up to degree ~8 precision is at least 10^-3 + * - sbasis_to_bezier : up to order ~3 precision is at least 10^-5 + * up to order ~4 precision is at least 10^-3 + * + * references: + * this implementation is based on the following article: + * J.Sanchez-Reyes - The Symmetric Analogue of the Polynomial Power Basis + */ + +inline +double binomial(unsigned int n, unsigned int k) +{ + return choose(n, k); +} + +inline +int sgn(unsigned int j, unsigned int k) +{ + assert (j >= k); + // we are sure that j >= k + return ((j-k) & 1u) ? -1 : 1; +} + + +void sbasis_to_bezier (Bezier & bz, SBasis const& sb, size_t sz) +{ + // if the degree is even q is the order in the symmetrical power basis, + // if the degree is odd q is the order + 1 + // n is always the polynomial degree, i. e. the Bezier order + size_t q, n; + bool even; + if (sz == 0) + { + q = sb.size(); + if (sb[q-1][0] == sb[q-1][1]) + { + even = true; + --q; + n = 2*q; + } + else + { + even = false; + n = 2*q-1; + } + } + else + { + q = (sz > sb.size()) ? sb.size() : sz; + n = 2*sz-1; + even = false; + } + bz.clear(); + bz.resize(n+1); + double Tjk; + for (size_t k = 0; k < q; ++k) + { + for (size_t j = k; j < n-k; ++j) // j <= n-k-1 + { + Tjk = binomial(n-2*k-1, j-k); + bz[j] += (Tjk * sb[k][0]); + bz[n-j] += (Tjk * sb[k][1]); // n-k <-> [k][1] + } + } + if (even) + { + bz[q] += sb[q][0]; + } + // the resulting coefficients are with respect to the scaled Bernstein + // basis so we need to divide them by (n, j) binomial coefficient + for (size_t j = 1; j < n; ++j) + { + bz[j] /= binomial(n, j); + } + bz[0] = sb[0][0]; + bz[n] = sb[0][1]; +} + +void sbasis_to_bezier (std::vector & bz, D2 const& sb, size_t sz) +{ + Bezier bzx, bzy; + sbasis_to_bezier(bzx, sb[X], sz); + sbasis_to_bezier(bzy, sb[Y], sz); + size_t n = (bzx.size() >= bzy.size()) ? bzx.size() : bzy.size(); + + bz.resize(n, Point(0,0)); + for (size_t i = 0; i < bzx.size(); ++i) + { + bz[i][X] = bzx[i]; + } + for (size_t i = 0; i < bzy.size(); ++i) + { + bz[i][Y] = bzy[i]; + } +} + + +void bezier_to_sbasis (SBasis & sb, Bezier const& bz) +{ + // if the degree is even q is the order in the symmetrical power basis, + // if the degree is odd q is the order + 1 + // n is always the polynomial degree, i. e. the Bezier order + size_t n = bz.order(); + size_t q = (n+1) / 2; + size_t even = (n & 1u) ? 0 : 1; + sb.clear(); + sb.resize(q + even, Linear(0, 0)); + double Tjk; + for (size_t k = 0; k < q; ++k) + { + for (size_t j = k; j < q; ++j) + { + Tjk = sgn(j, k) * binomial(n-j-k, j-k) * binomial(n, k); + sb[j][0] += (Tjk * bz[k]); + sb[j][1] += (Tjk * bz[n-k]); // n-j <-> [j][1] + } + for (size_t j = k+1; j < q; ++j) + { + Tjk = sgn(j, k) * binomial(n-j-k-1, j-k-1) * binomial(n, k); + sb[j][0] += (Tjk * bz[n-k]); + sb[j][1] += (Tjk * bz[k]); // n-j <-> [j][1] + } + } + if (even) + { + for (size_t k = 0; k < q; ++k) + { + Tjk = sgn(q,k) * binomial(n, k); + sb[q][0] += (Tjk * (bz[k] + bz[n-k])); + } + sb[q][0] += (binomial(n, q) * bz[q]); + sb[q][1] = sb[q][0]; + } + sb[0][0] = bz[0]; + sb[0][1] = bz[n]; +} + + +void bezier_to_sbasis (D2 & sb, std::vector const& bz) +{ + size_t n = bz.size() - 1; + size_t q = (n+1) / 2; + size_t even = (n & 1u) ? 0 : 1; + sb[X].clear(); + sb[Y].clear(); + sb[X].resize(q + even, Linear(0, 0)); + sb[Y].resize(q + even, Linear(0, 0)); + double Tjk; + for (size_t k = 0; k < q; ++k) + { + for (size_t j = k; j < q; ++j) + { + Tjk = sgn(j, k) * binomial(n-j-k, j-k) * binomial(n, k); + sb[X][j][0] += (Tjk * bz[k][X]); + sb[X][j][1] += (Tjk * bz[n-k][X]); + sb[Y][j][0] += (Tjk * bz[k][Y]); + sb[Y][j][1] += (Tjk * bz[n-k][Y]); + } + for (size_t j = k+1; j < q; ++j) + { + Tjk = sgn(j, k) * binomial(n-j-k-1, j-k-1) * binomial(n, k); + sb[X][j][0] += (Tjk * bz[n-k][X]); + sb[X][j][1] += (Tjk * bz[k][X]); + sb[Y][j][0] += (Tjk * bz[n-k][Y]); + sb[Y][j][1] += (Tjk * bz[k][Y]); + } + } + if (even) + { + for (size_t k = 0; k < q; ++k) + { + Tjk = sgn(q,k) * binomial(n, k); + sb[X][q][0] += (Tjk * (bz[k][X] + bz[n-k][X])); + sb[Y][q][0] += (Tjk * (bz[k][Y] + bz[n-k][Y])); + } + sb[X][q][0] += (binomial(n, q) * bz[q][X]); + sb[X][q][1] = sb[X][q][0]; + sb[Y][q][0] += (binomial(n, q) * bz[q][Y]); + sb[Y][q][1] = sb[Y][q][0]; + } + sb[X][0][0] = bz[0][X]; + sb[X][0][1] = bz[n][X]; + sb[Y][0][0] = bz[0][Y]; + sb[Y][0][1] = bz[n][Y]; +} + + +} // end namespace Geom + +namespace Geom{ +#if 0 + /* From Sanchez-Reyes 1997 W_{j,k} = W_{n0j, n-k} = choose(n-2k-1, j-k)choose(2k+1,k)/choose(n,j) for k=0,...,q-1; j = k, ...,n-k-1 @@ -9,14 +266,6 @@ This is wrong, it should read W_{q,q} = 1 (n even) */ -#include <2geom/sbasis-to-bezier.h> -#include <2geom/choose.h> -#include <2geom/svg-path.h> -#include -#include <2geom/exception.h> - -namespace Geom{ - double W(unsigned n, unsigned j, unsigned k) { unsigned q = (n+1)/2; if((n & 1) == 0 && j == q && k == q) @@ -32,6 +281,7 @@ double W(unsigned n, unsigned j, unsigned k) { choose(n,j); } + // this produces a degree 2q bezier from a degree k sbasis Bezier sbasis_to_bezier(SBasis const &B, unsigned q) { @@ -59,6 +309,7 @@ double mopi(int i) { return (i&1)?-1:1; } +// WARNING: this is wrong! // this produces a degree k sbasis from a degree 2q bezier SBasis bezier_to_sbasis(Bezier const &B) { @@ -106,6 +357,7 @@ D2 > sbasis_to_bezier(D2 const &B) { return D2 >(sbasis_to_bezier(B[0]), sbasis_to_bezier(B[1])); } */ +#endif #if 0 // using old path //std::vector @@ -135,7 +387,7 @@ subpath_from_sbasis(Geom::OldPathSetBuilder &pb, D2 const &B, double tol /* * This version works by inverting a reasonable upper bound on the error term after subdividing the * curve at $a$. We keep biting off pieces until there is no more curve left. -* +* * Derivation: The tail of the power series is $a_ks^k + a_{k+1}s^{k+1} + \ldots = e$. A * subdivision at $a$ results in a tail error of $e*A^k, A = (1-a)a$. Let this be the desired * tolerance tol $= e*A^k$ and invert getting $A = e^{1/k}$ and $a = 1/2 - \sqrt{1/4 - A}$ @@ -146,7 +398,7 @@ subpath_from_sbasis_incremental(Geom::OldPathSetBuilder &pb, D2 B, doubl double te = B.tail_error(k); assert(B[0].IS_FINITE()); assert(B[1].IS_FINITE()); - + //std::cout << "tol = " << tol << std::endl; while(1) { double A = std::sqrt(tol/te); // pow(te, 1./k) @@ -169,10 +421,10 @@ subpath_from_sbasis_incremental(Geom::OldPathSetBuilder &pb, D2 B, doubl initial = false; } pb.push_cubic(bez[1], bez[2], bez[3]); - + // move to next piece of curve if(a >= 1) break; - B = compose(B, Linear(a, 1)); + B = compose(B, Linear(a, 1)); te = B.tail_error(k); } } @@ -190,7 +442,8 @@ void build_from_sbasis(Geom::PathBuilder &pb, D2 const &B, double tol, b if( !only_cubicbeziers && (sbasis_size(B) <= 1) ) { pb.lineTo(B.at1()); } else { - std::vector bez = sbasis_to_bezier(B, 2); + std::vector bez; + sbasis_to_bezier(bez, B, 2); pb.curveTo(bez[1], bez[2], bez[3]); } } else { @@ -246,7 +499,7 @@ path_from_piecewise(Geom::Piecewise > const &B, double to return pb.peek(); } -}; +} /* Local Variables: diff --git a/src/2geom/sbasis-to-bezier.h b/src/2geom/sbasis-to-bezier.h index 2fedb83c7..c9d5cbbbc 100644 --- a/src/2geom/sbasis-to-bezier.h +++ b/src/2geom/sbasis-to-bezier.h @@ -4,7 +4,17 @@ #include <2geom/d2.h> #include <2geom/path.h> -namespace Geom{ +#include + +namespace Geom { + +void sbasis_to_bezier (Bezier & bz, SBasis const& sb, size_t sz = 0); +void sbasis_to_bezier (std::vector & bz, D2 const& sb, size_t sz = 0); +void bezier_to_sbasis (SBasis & sb, Bezier const& bz); +void bezier_to_sbasis (D2 & sb, std::vector const& bz); + + +#if 0 // this produces a degree k bezier from a degree k sbasis Bezier sbasis_to_bezier(SBasis const &B, unsigned q = 0); @@ -15,14 +25,19 @@ SBasis bezier_to_sbasis(Bezier const &B); std::vector sbasis_to_bezier(D2 const &B, unsigned q = 0); +#endif + std::vector path_from_piecewise(Piecewise > const &B, double tol, bool only_cubicbeziers = false); Path path_from_sbasis(D2 const &B, double tol, bool only_cubicbeziers = false); inline Path cubicbezierpath_from_sbasis(D2 const &B, double tol) - { return path_from_sbasis(B, tol, true); } ; + { return path_from_sbasis(B, tol, true); } + +} // end namespace Geom + + -}; #endif /* diff --git a/src/2geom/solve-bezier-one-d.cpp b/src/2geom/solve-bezier-one-d.cpp index b922e970e..9a7b5633f 100644 --- a/src/2geom/solve-bezier-one-d.cpp +++ b/src/2geom/solve-bezier-one-d.cpp @@ -14,8 +14,8 @@ static int SGN(t x) { return (x > 0 ? 1 : (x < 0 ? -1 : 0)); } const unsigned MAXDEPTH = 23; // Maximum depth for recursion. Using floats means 23 bits precision max -const double BEPSILON = ldexp(1.0,-MAXDEPTH-1); /*Flatness control value */ - +const double BEPSILON = ldexp(1.0,(-MAXDEPTH-1)); /*Flatness control value */ +const double SECANT_EPSILON = 1e-13; // secant method converges much faster, get a bit more precision /** * This function is called _a lot_. We have included various manual memory management stuff to reduce the amount of mallocing that goes on. In the future it is possible that this will hurt performance. **/ @@ -24,7 +24,8 @@ public: double *Vtemp; unsigned N,degree; std::vector &solutions; - Bernsteins(int degr, std::vector &so) : N(degr+1), degree(degr),solutions(so) { + bool use_secant; + Bernsteins(int degr, std::vector &so) : N(degr+1), degree(degr),solutions(so), use_secant(false) { Vtemp = new double[N*2]; } ~Bernsteins() { @@ -35,6 +36,8 @@ public: double *Left, double *Right); + double horner(const double *b, double t); + unsigned control_poly_flat_enough(double const *V); @@ -52,9 +55,10 @@ find_bernstein_roots(double const *w, /* The control points */ unsigned degree, /* The degree of the polynomial */ std::vector &solutions, /* RETURN candidate t-values */ unsigned depth, /* The depth of the recursion */ - double left_t, double right_t) + double left_t, double right_t, bool use_secant) { Bernsteins B(degree, solutions); + B.use_secant = use_secant; B.find_bernstein_roots(w, depth, left_t, right_t); } @@ -96,17 +100,37 @@ Bernsteins::find_bernstein_roots(double const *w, /* The control points */ } // I thought secant method would be faster here, but it'aint. -- njh - - // The original code went to some effort to try and terminate early when the curve is sufficiently flat. However, it seems that in practice it almost always bottoms out, so leaving this code out actually speeds things up - if(0) if (control_poly_flat_enough(w)) { - //printf("flatten out %d\n", depth); - const double Ax = right_t - left_t; - const double Ay = w[degree] - w[0]; - - solutions.push_back(left_t - Ax*w[0] / Ay); - return; + // Actually, it was, I just was using the wrong method for bezier evaluation. Horner's rule results in a very efficient algorithm - 10* faster (20080816) + // Future work: try using brent's method + if(use_secant) { // false position + double s = 0;double t = 1; + double e = 1e-10; + int n,side=0; + double r,fr,fs = w[0],ft = w[degree]; + + for (n = 1; n <= 100; n++) + { + r = (fs*t - ft*s) / (fs - ft); + if (fabs(t-s) < e*fabs(t+s)) break; + fr = horner(w, r); + + if (fr * ft > 0) + { + t = r; ft = fr; + if (side==-1) fs /= 2; + side = -1; + } + else if (fs * fr > 0) + { + s = r; fs = fr; + if (side==+1) ft /= 2; + side = +1; + } + else break; } - + solutions.push_back(r*right_t + (1-r)*left_t); + return; + } } /* Otherwise, solve recursively after subdividing control polygon */ @@ -195,6 +219,22 @@ Bernsteins::control_poly_flat_enough(double const *V) return error < BEPSILON * a; } +// suggested by Sederberg. +double Bernsteins::horner(const double *b, double t) { + int n = degree; + double u, bc, tn, tmp; + int i; + u = 1.0 - t; + bc = 1; + tn = 1; + tmp = b[0]*u; + for(i=1; i & solutions, /* RETURN candidate t-values */ unsigned depth, /* The depth of the recursion */ - double left_t=0, double right_t=1); + double left_t=0, double right_t=1, bool use_secant=true); }; #endif diff --git a/src/2geom/svg-elliptical-arc.cpp b/src/2geom/svg-elliptical-arc.cpp index 5ae92097e..23645fe11 100644 --- a/src/2geom/svg-elliptical-arc.cpp +++ b/src/2geom/svg-elliptical-arc.cpp @@ -335,7 +335,10 @@ SVGEllipticalArc::roots(double v, Dim2 d) const } -// D(E(t,C),t) = E(t+PI/2,O) +// D(E(t,C),t) = E(t+PI/2,O), where C is the ellipse center +// the derivative doesn't rotate the ellipse but there is a translation +// of the parameter t by an angle of PI/2 so the ellipse points are shifted +// of such an angle in the cw direction Curve* SVGEllipticalArc::derivative() const { if (isDegenerate() && is_svg_compliant()) @@ -411,10 +414,12 @@ bool SVGEllipticalArc::containsAngle(Coord angle) const Curve* SVGEllipticalArc::portion(double f, double t) const { + // fix input arguments if (f < 0) f = 0; if (f > 1) f = 1; if (t < 0) t = 0; if (t > 1) t = 1; + if ( are_near(f, t) ) { SVGEllipticalArc* arc = new SVGEllipticalArc(); @@ -424,6 +429,7 @@ Curve* SVGEllipticalArc::portion(double f, double t) const arc->m_sweep = m_sweep; arc->m_large_arc = m_large_arc; } + SVGEllipticalArc* arc = new SVGEllipticalArc( *this ); arc->m_initial_point = pointAt(f); arc->m_final_point = pointAt(t); @@ -909,8 +915,6 @@ D2 SVGEllipticalArc::toSBasis() const Coord cos_rot_angle = std::cos(rotation_angle()); Coord sin_rot_angle = std::sin(rotation_angle()); // order = 4 seems to be enough to get a perfect looking elliptical arc - // should it be choosen in function of the arc length anyway ? - // or maybe a user settable parameter: toSBasis(unsigned int order) ? SBasis arc_x = ray(X) * cos(param,4); SBasis arc_y = ray(Y) * sin(param,4); arc[0] = arc_x * cos_rot_angle - arc_y * sin_rot_angle + Linear(center(X),center(X)); @@ -928,8 +932,6 @@ D2 SVGEllipticalArc::toSBasis() const Curve* SVGEllipticalArc::transformed(Matrix const& m) const { - // return SBasisCurve(toSBasis()).transformed(m); - Ellipse e(center(X), center(Y), ray(X), ray(Y), rotation_angle()); Ellipse et = e.transformed(m); Point inner_point = pointAt(0.5); @@ -940,7 +942,11 @@ Curve* SVGEllipticalArc::transformed(Matrix const& m) const return ea.duplicate(); } - +/* + * helper routine to convert the parameter t value + * btw [0,1] and [0,2PI] domain and back + * + */ Coord SVGEllipticalArc::map_to_02PI(Coord t) const { Coord angle = start_angle(); @@ -966,7 +972,14 @@ Coord SVGEllipticalArc::map_to_01(Coord angle) const namespace detail { - +/* + * ellipse_equation + * + * this is an helper struct, it provides two routines: + * the first one evaluates the implicit form of an ellipse on a given point + * the second one computes the normal versor at a given point of an ellipse + * in implicit form + */ struct ellipse_equation { ellipse_equation(double a, double b, double c, double d, double e, double f) @@ -1017,6 +1030,10 @@ make_elliptical_arc( SVGEllipticalArc& _ea, { } +/* + * check that the coefficients computed by the fit method satisfy + * the tolerance parameters at the k-th sample point + */ bool make_elliptical_arc:: bound_exceeded( unsigned int k, detail::ellipse_equation const & ee, @@ -1024,39 +1041,50 @@ bound_exceeded( unsigned int k, detail::ellipse_equation const & ee, { dist_err = std::fabs( ee(p[k]) ); dist_bound = std::fabs( e1x * p[k][X] + e1y * p[k][Y] + e2 ); + // check that the angle btw the tangent versor to the input curve + // and the normal versor of the elliptical arc, both evaluate + // at the k-th sample point, are really othogonal angle_err = std::fabs( dot( dcurve(k/partitions), ee.normal(p[k]) ) ); //angle_err *= angle_err; return ( dist_err > dist_bound || angle_err > angle_tol ); } +/* + * check that the coefficients computed by the fit method satisfy + * the tolerance parameters at each sample point + */ bool make_elliptical_arc:: check_bound(double A, double B, double C, double D, double E, double F) { - // check error magnitude detail::ellipse_equation ee(A, B, C, D, E, F); + // check error magnitude at the end-points double e1x = (2*A + B) * tol_at_extr; double e1y = (B + 2*C) * tol_at_extr; double e2 = ((D + E) + (A + B + C) * tol_at_extr) * tol_at_extr; - if ( bound_exceeded(0, ee, e1x, e1y, e2) ) + if (bound_exceeded(0, ee, e1x, e1y, e2)) { print_bound_error(0); return false; } - if ( bound_exceeded(0, ee, e1x, e1y, e2) ) + if (bound_exceeded(0, ee, e1x, e1y, e2)) { print_bound_error(last); return false; } + // e1x = derivative((ee(x,y), x) | x->tolerance, y->tolerance e1x = (2*A + B) * tolerance; + // e1y = derivative((ee(x,y), y) | x->tolerance, y->tolerance e1y = (B + 2*C) * tolerance; + // e2 = ee(tolerance, tolerance) - F; e2 = ((D + E) + (A + B + C) * tolerance) * tolerance; // std::cerr << "e1x = " << e1x << std::endl; // std::cerr << "e1y = " << e1y << std::endl; // std::cerr << "e2 = " << e2 << std::endl; + // check error magnitude at sample points for ( unsigned int k = 1; k < last; ++k ) { if ( bound_exceeded(k, ee, e1x, e1y, e2) ) @@ -1069,6 +1097,12 @@ check_bound(double A, double B, double C, double D, double E, double F) return true; } +/* + * fit + * + * supply the samples to the fitter and compute + * the ellipse implicit equation coefficients + */ void make_elliptical_arc::fit() { for (unsigned int k = 0; k < N; ++k) diff --git a/src/2geom/svg-elliptical-arc.h b/src/2geom/svg-elliptical-arc.h index 7d5087c48..5c42e6e1d 100644 --- a/src/2geom/svg-elliptical-arc.h +++ b/src/2geom/svg-elliptical-arc.h @@ -60,12 +60,34 @@ class SVGEllipticalArc : public Curve : m_initial_point(Point(0,0)), m_final_point(Point(0,0)), m_rx(0), m_ry(0), m_rot_angle(0), m_large_arc(true), m_sweep(true), - m_svg_compliant(_svg_compliant) - { - m_start_angle = m_end_angle = 0; - m_center = Point(0,0); - } - + m_svg_compliant(_svg_compliant), + m_start_angle(0), m_end_angle(0), + m_center(Point(0,0)) + { + } + + /* + * constructor + * + * input parameters: + * _initial_point: initial arc end point; + * _rx: ellipse x-axis ray length + * _ry: ellipse y-axis ray length + * _rot_angle: ellipse x-axis rotation angle; + * _large_arc: if true the largest arc is chosen, + * if false the smallest arc is chosen; + * _sweep : if true the clockwise arc is chosen, + * if false the counter-clockwise arc is chosen; + * _final_point: final arc end point; + * _svg_compliant: if true the class behaviour follows the Standard + * SVG 1.1 implementation guidelines (see Appendix F.6) + * if false the class behavoiur is more strict + * on input parameter + * + * in case the initial and the final arc end-points overlaps + * a degenerate arc of zero length is generated + * + */ SVGEllipticalArc( Point _initial_point, double _rx, double _ry, double _rot_angle, bool _large_arc, bool _sweep, Point _final_point, @@ -196,9 +218,19 @@ class SVGEllipticalArc : public Curve std::vector roots(double v, Dim2 d) const; + /* + * find all the points on the curve portion between "from" and "to" + * at the same smallest distance from the point "p" the points are returned + * as their parameter t value; + */ std::vector allNearestPoints( Point const& p, double from = 0, double to = 1 ) const; + /* + * find a point on the curve portion between "from" and "to" + * at the same smallest distance from the point "p"; + * the point is returned as its parameter t value; + */ double nearestPoint( Point const& p, double from = 0, double to = 1 ) const { if ( are_near(ray(X), ray(Y)) && are_near(center(), p) ) @@ -225,10 +257,22 @@ class SVGEllipticalArc : public Curve D2 toSBasis() const; + /* + * return true if the angle argument (in radiants) is contained + * in the range [start_angle(), end_angle() ] + */ bool containsAngle(Coord angle) const; + /* + * return the value of the d-dimensional coordinate related to "t" + * here t belongs to the [0,2PI] domain + */ double valueAtAngle(Coord t, Dim2 d) const; + /* + * return the point related to the parameter value "t" + * here t belongs to the [0,2PI] domain + */ Point pointAtAngle(Coord t) const { double sin_rot_angle = std::sin(rotation_angle()); @@ -240,6 +284,10 @@ class SVGEllipticalArc : public Curve return p * m; } + /* + * return the value of the d-dimensional coordinate related to "t" + * here t belongs to the [0,1] domain + */ double valueAt(Coord t, Dim2 d) const { if (isDegenerate() && is_svg_compliant()) @@ -249,6 +297,10 @@ class SVGEllipticalArc : public Curve return valueAtAngle(tt, d); } + /* + * return the point related to the parameter value "t" + * here t belongs to the [0,1] domain + */ Point pointAt(Coord t) const { if (isDegenerate() && is_svg_compliant()) @@ -284,6 +336,7 @@ class SVGEllipticalArc : public Curve return rarc; } + double sweep_angle() const { Coord d = end_angle() - start_angle(); @@ -307,12 +360,16 @@ class SVGEllipticalArc : public Curve Point m_initial_point, m_final_point; double m_rx, m_ry, m_rot_angle; bool m_large_arc, m_sweep; + bool m_svg_compliant; double m_start_angle, m_end_angle; Point m_center; - bool m_svg_compliant; }; // end class SVGEllipticalArc + +/* + * useful for testing and debugging + */ template< class charT > inline std::basic_ostream & @@ -331,17 +388,44 @@ operator<< (std::basic_ostream & os, const SVGEllipticalArc & ea) +// forward declation namespace detail { struct ellipse_equation; } - +/* + * make_elliptical_arc + * + * convert a parametric polynomial curve given in symmetric power basis form + * into an SVGEllipticalArc type; in order to be successfull the input curve + * has to look like an actual elliptical arc even if a certain tolerance + * is allowed through an ad-hoc parameter. + * The conversion is performed through an interpolation on a certain amount of + * sample points computed on the input curve; + * the interpolation computes the coefficients of the general implicit equation + * of an ellipse (A*X^2 + B*XY + C*Y^2 + D*X + E*Y + F = 0), then from the + * implicit equation we compute the parametric form. + * + */ class make_elliptical_arc { public: typedef D2 curve_type; + /* + * constructor + * + * it doesn't execute the conversion but set the input and output parameters + * + * _ea: the output SVGEllipticalArc that will be generated; + * _curve: the input curve to be converted; + * _total_samples: the amount of sample points to be taken + * on the input curve for performing the conversion + * _tolerance: how much likelihood is required between the input curve + * and the generated elliptical arc; the smaller it is the + * the tolerance the higher it is the likelihood. + */ make_elliptical_arc( SVGEllipticalArc& _ea, curve_type const& _curve, unsigned int _total_samples, @@ -369,8 +453,13 @@ class make_elliptical_arc } public: + /* + * perform the actual conversion + * return true if the conversion is successfull, false on the contrary + */ bool operator()() { + // initialize the reference const NL::Vector & coeff = fitter.result(); fit(); if ( !check_bound(1, coeff[0], coeff[1], coeff[2], coeff[3], coeff[4]) ) @@ -379,6 +468,11 @@ class make_elliptical_arc return true; } + /* + * you can set a boolean parameter to tell the conversion routine + * if the output elliptical arc has to be svg compliant or not; + * the default value is true + */ bool svg_compliant_flag() const { return svg_compliant; @@ -390,17 +484,24 @@ class make_elliptical_arc } private: - SVGEllipticalArc& ea; - const curve_type & curve; - Piecewise > dcurve; - NL::LFMEllipse model; + SVGEllipticalArc& ea; // output elliptical arc + const curve_type & curve; // input curve + Piecewise > dcurve; // derivative of the input curve + NL::LFMEllipse model; // model used for fitting + // perform the actual fitting task NL::least_squeares_fitter fitter; + // tolerance: the user-defined tolerance parameter; + // tol_at_extr: the tolerance at end-points automatically computed + // on the value of "tolerance", and usually more strict; + // tol_at_center: tolerance at the center of the ellipse + // angle_tol: tolerance for the angle btw the input curve tangent + // versor and the ellipse normal versor at the sample points double tolerance, tol_at_extr, tol_at_center, angle_tol; - Point initial_point, final_point; - unsigned int N; + Point initial_point, final_point; // initial and final end-points + unsigned int N; // total samples unsigned int last; // N-1 double partitions; // N-1 - std::vector p; // sample points + std::vector p; // sample points double dist_err, dist_bound, angle_err; bool svg_compliant; }; diff --git a/src/2geom/svg-path.cpp b/src/2geom/svg-path.cpp index bed04d3ff..898c72bf5 100644 --- a/src/2geom/svg-path.cpp +++ b/src/2geom/svg-path.cpp @@ -35,7 +35,8 @@ namespace Geom { void output(Curve const &curve, SVGPathSink &sink) { - std::vector pts = sbasis_to_bezier(curve.toSBasis(), 2); //TODO: use something better! + std::vector pts; + sbasis_to_bezier(pts, curve.toSBasis(), 2); //TODO: use something better! sink.curveTo(pts[0], pts[1], pts[2]); } diff --git a/src/2geom/sweep.cpp b/src/2geom/sweep.cpp index 1b8813990..227c822bd 100644 --- a/src/2geom/sweep.cpp +++ b/src/2geom/sweep.cpp @@ -9,8 +9,10 @@ std::vector > sweep_bounds(std::vector rs) { std::vector > pairs(rs.size()); for(unsigned i = 0; i < rs.size(); i++) { - events.push_back(Event(rs[i].left(), i, false)); - events.push_back(Event(rs[i].right(), i, true)); + if(!rs[i].isEmpty()) { + events.push_back(Event(rs[i].left(), i, false)); + events.push_back(Event(rs[i].right(), i, true)); + } } std::sort(events.begin(), events.end()); @@ -45,8 +47,11 @@ std::vector > sweep_bounds(std::vector a, std::vecto unsigned sz = n ? b.size() : a.size(); events[n].reserve(sz*2); for(unsigned i = 0; i < sz; i++) { - events[n].push_back(Event(n ? b[i].left() : a[i].left(), i, false)); - events[n].push_back(Event(n ? b[i].right() : a[i].right(), i, true)); + Rect r = n ? b[i] : a[i]; + if(!r.isEmpty()) { + events[n].push_back(Event(r.left(), i, false)); + events[n].push_back(Event(r.right(), i, true)); + } } std::sort(events[n].begin(), events[n].end()); }