From: johanengelen Date: Thu, 30 Aug 2007 18:32:36 +0000 (+0000) Subject: Update to 2Geom rev. 1113 X-Git-Url: https://git.tokkee.org/?a=commitdiff_plain;h=29684a16b6c92bee28a94fdc2607bcc143950fa8;p=inkscape.git Update to 2Geom rev. 1113 --- diff --git a/src/2geom/Makefile_insert b/src/2geom/Makefile_insert index a100c418d..4a87c1a93 100644 --- a/src/2geom/Makefile_insert +++ b/src/2geom/Makefile_insert @@ -9,6 +9,8 @@ 2geom/bezier-to-sbasis.h \ 2geom/bezier-utils.cpp \ 2geom/bezier-utils.h \ + 2geom/chebyshev.cpp \ + 2geom/chebyshev.h \ 2geom/choose.h \ 2geom/circle-circle.cpp \ 2geom/circulator.h \ @@ -18,6 +20,10 @@ 2geom/convex-cover.cpp \ 2geom/convex-cover.h \ 2geom/coord.h \ + 2geom/crossing.cpp \ + 2geom/crossing.cpp \ + 2geom/d2-sbasis.cpp \ + 2geom/d2-sbasis.h \ 2geom/d2.cpp \ 2geom/d2.h \ 2geom/geom.cpp \ @@ -27,6 +33,9 @@ 2geom/linear.h \ 2geom/matrix.cpp \ 2geom/matrix.h \ + 2geom/ord.h \ + 2geom/path-intersection.cpp \ + 2geom/path-intersection.h \ 2geom/path.cpp \ 2geom/path.h \ 2geom/piecewise.cpp \ @@ -44,6 +53,8 @@ 2geom/quadtree.cpp \ 2geom/quadtree.h \ 2geom/rect.h \ + 2geom/region.cpp \ + 2geom/region.h \ 2geom/sbasis-2d.cpp \ 2geom/sbasis-2d.h \ 2geom/sbasis.cpp \ @@ -57,6 +68,8 @@ 2geom/sbasis-roots.cpp \ 2geom/sbasis-to-bezier.cpp \ 2geom/sbasis-to-bezier.h \ + 2geom/shape.cpp \ + 2geom/shape.h \ 2geom/solve-bezier-one-d.cpp \ 2geom/solve-bezier-parametric.cpp \ 2geom/solver.h \ @@ -65,6 +78,8 @@ 2geom/svg-path.h \ 2geom/svg-path-parser.cpp \ 2geom/svg-path-parser.h \ + 2geom/sweep.cpp \ + 2geom/sweep.h \ 2geom/transforms.cpp \ 2geom/transforms.h \ 2geom/utils.h \ diff --git a/src/2geom/basic-intersection.cpp b/src/2geom/basic-intersection.cpp index a5b827023..694760d5a 100644 --- a/src/2geom/basic-intersection.cpp +++ b/src/2geom/basic-intersection.cpp @@ -159,7 +159,7 @@ void OldBezier::split(double t, OldBezier &left, OldBezier &right) const { * the various edges of the bounding box of the first curve to test * for interference. * Second, after a few subdivisions it is highly probable that two corners - * of the bounding box of a given OldBezier curve are the first and last + * of the bounding box of a given Bezier curve are the first and last * control point. Once this happens once, it happens for all subsequent * subcurves. It might be worth putting in a test and then short-circuit * code for further subdivision levels. diff --git a/src/2geom/bezier.h b/src/2geom/bezier.h index 7949e006c..2c1b49213 100644 --- a/src/2geom/bezier.h +++ b/src/2geom/bezier.h @@ -3,6 +3,7 @@ * * Copyright 2007 MenTaLguY * Copyright 2007 Michael Sloan + * Copyright 2007 Nathan Hurst * * This library is free software; you can redistribute it and/or * modify it either under the terms of the GNU Lesser General Public @@ -35,107 +36,246 @@ #include "coord.h" #include "isnan.h" #include "bezier-to-sbasis.h" +#include "d2.h" +#include "solver.h" +#include namespace Geom { +template +Coord subdivideArr(Coord t, Coord const *v, Coord *left, Coord *right) { + Coord vtemp[order+1][order+1]; + + /* Copy control points */ + std::copy(v, v+order+1, vtemp[0]); + + /* Triangle computation */ + 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]); + } + } + if(left != NULL) + for (unsigned j = 0; j <= order; j++) + left[j] = vtemp[j][0]; + if(right != NULL) + for (unsigned j = 0; j <= order; j++) + right[j] = vtemp[order-j][j]; + + return (vtemp[order][0]); +} + + template class Bezier { private: - Coord c_[order+1]; + Coord c_[order+1]; + + template + friend Bezier portion(const Bezier & a, Coord from, Coord to); + + template + friend Interval bounds_fast(Bezier const & b); + + template + friend Bezier derivative(const Bezier & a); protected: - Bezier(Coord c[]) { - std::copy(c, c+order+1, c_); - } + Bezier(Coord const c[]) { + std::copy(c, c+order+1, c_); + } public: - template - static void assert_order(Bezier const *) {} - - Bezier() {} - - //Construct an order-0 bezier (constant Bézier) - explicit Bezier(Coord c0) { - assert_order<0>(this); - c_[0] = c0; - } - - //Construct an order-1 bezier (linear Bézier) - Bezier(Coord c0, Coord c1) { - assert_order<1>(this); - c_[0] = c0; c_[1] = c1; - } - - //Construct an order-2 bezier (quadratic Bézier) - Bezier(Coord c0, Coord c1, Coord c2) { - assert_order<2>(this); - c_[0] = c0; c_[1] = c1; c_[2] = c2; - } - - //Construct an order-3 bezier (cubic Bézier) - Bezier(Coord c0, Coord c1, Coord c2, Coord c3) { - assert_order<3>(this); - c_[0] = c0; c_[1] = c1; c_[2] = c2; c_[3] = c3; - } - - inline unsigned degree() const { return order; } - - //IMPL: FragmentConcept - typedef Coord output_type; - inline bool isZero() const { - for(int i = 0; i <= order; i++) { - if(c_[i] != 0) return false; - } - return true; - } - inline bool isFinite() const { - for(int i = 0; i <= order; i++) { - if(!is_finite(c_[i])) return false; + + //TODO: fix this assert - get it compiling! + //template + //static void assert_order(Bezier const *) {} + + Bezier() {} + + //Construct an order-0 bezier (constant Bézier) + explicit Bezier(Coord c0) { + //assert_order<0>(this); + c_[0] = c0; + } + + //Construct an order-1 bezier (linear Bézier) + Bezier(Coord c0, Coord c1) { + //assert_order<1>(this); + c_[0] = c0; c_[1] = c1; + } + + //Construct an order-2 bezier (quadratic Bézier) + Bezier(Coord c0, Coord c1, Coord c2) { + //assert_order<2>(this); + c_[0] = c0; c_[1] = c1; c_[2] = c2; + } + + //Construct an order-3 bezier (cubic Bézier) + Bezier(Coord c0, Coord c1, Coord c2, Coord c3) { + //assert_order<3>(this); + c_[0] = c0; c_[1] = c1; c_[2] = c2; c_[3] = c3; } - return true; - } - inline Coord at0() const { return c_[0]; } - inline Coord at1() const { return c_[order]; } - inline SBasis toSBasis() const { return bezier_to_sbasis(c_); } + inline unsigned degree() const { return order; } - inline Interval bounds_fast() const { return Interval::fromArray(c_, order+1); } - //TODO: better bounds exact - inline Interval bounds_exact() const { return toSBasis().bounds_exact(); } - inline Interval bounds_local(double u, double v) const { return toSBasis().bounds_local(u, v); } + //IMPL: FragmentConcept + typedef Coord output_type; + inline bool isZero() const { + for(unsigned i = 0; i <= order; i++) { + if(c_[i] != 0) return false; + } + return true; + } + inline bool isFinite() const { + for(unsigned i = 0; i <= order; i++) { + if(!is_finite(c_[i])) return false; + } + return true; + } + inline Coord at0() const { return c_[0]; } + inline Coord at1() const { return c_[order]; } - //Only mutator - inline Coord &operator[](int index) { return c_[index]; } - inline Coord const &operator[](int index) const { return c_[index]; } + inline Coord valueAt(double t) const { return subdivideArr(t, c_, NULL, NULL); } + inline Coord operator()(double t) const { return valueAt(t); } - Maybe winding(Point p) const { - return sbasis_winding(toSBasis(), p); - } + inline SBasis toSBasis() const { return bezier_to_sbasis(c_); } + + //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 { + throw NotImplemented(); + // Can't be implemented without a dynamic version of subdivide. + /*std::vector val_n_der; + Coord d_[order+1]; + for(int di = n_derivs; di > 0; di--) { + val_n_der.push_back(subdivideArr(t, d_, NULL, NULL)); + for(unsigned i = 0; i < di; i++) { + d[i] = order*(a._c[i+1] - a._c[i]); + } + }*/ + } - Point pointAndDerivativesAt(Coord t, unsigned n_derivs, Point *derivs) const { - //TODO - return Point(0,0); - } + std::pair, Bezier > subdivide(Coord t) const { + Bezier a, b; + subdivideArr(t, order, c_, a.c_, b.c_); + return std::pair, Bezier >(a, b); + } + + std::vector roots() const { + std::vector solutions; + find_bernstein_roots(c_, order, solutions, 0, 0.0, 1.0); + return solutions; + } }; +//TODO: implement others +template +Bezier operator+(const Bezier & a, double v) { + Bezier result; + for(unsigned i = 0; i <= order; i++) + result[i] = a[i] + v; + return result; +} +template +Bezier operator-(const Bezier & a, double v) { + Bezier result; + for(unsigned i = 0; i <= order; i++) + result[i] = a[i] - v; + return result; +} +template +Bezier operator*(const Bezier & a, double v) { + Bezier result; + for(unsigned i = 0; i <= order; i++) + result[i] = a[i] * v; + return result; +} +template +Bezier operator/(const Bezier & a, double v) { + Bezier result; + for(unsigned i = 0; i <= order; i++) + result[i] = a[i] / v; + return result; +} + template Bezier reverse(const Bezier & a) { - Bezier result; - for(int i = 0; i <= order; i++) - result[i] = a[order - i]; - return result; + Bezier result; + for(unsigned i = 0; i <= order; i++) + result[i] = a[order - i]; + return result; +} + +template +Bezier portion(const Bezier & a, double from, double to) { + //TODO: implement better? + Coord res[order+1]; + if(from == 0) { + if(to == 1) { return Bezier(a); } + subdivideArr(to, a.c_, res, NULL); + return Bezier(res); + } + subdivideArr(from, a.c_, NULL, res); + if(to == 1) return Bezier(res); + Coord res2[order+1]; + subdivideArr((to - from)/(1 - from), res, res2, NULL); + return Bezier(res2); +} + +template +std::vector bezier_points(const D2 > & a) { + std::vector result; + for(unsigned i = 0; i <= order; i++) { + Point p; + for(unsigned d = 0; d < 2; d++) p[d] = a[d][i]; + result.push_back(p); + } + return result; +} + +template +Bezier derivative(const Bezier & a) { + Bezier der; + + for(unsigned i = 0; i < order; i++) { + der.c_[i] = order*(a.c_[i+1] - a.c_[i]); + } + return der; } template -vector bezier_points(const D2 > & a) { - vector result; - for(int i = 0; i <= order; i++) { - Point p; - for(unsigned d = 0; d < 2; d++) p[d] = a[d][i]; - result[i] = p; - } - return result; +inline Interval bounds_fast(Bezier const & b) { + return Interval::fromArray(b.c_, order+1); +} + +//TODO: better bounds exact +template +inline Interval bounds_exact(Bezier const & b) { + return bounds_exact(b.toSBasis()); +} + +template +inline Interval bounds_local(Bezier const & b, Interval i) { + return bounds_fast(portion(b, i.min(), i.max())); + //return bounds_local(b.toSBasis(), i); } } #endif //SEEN_BEZIER_H +/* + Local Variables: + mode:c++ + c-file-style:"stroustrup" + c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) + indent-tabs-mode:nil + fill-column:99 + End: +*/ +// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 : diff --git a/src/2geom/chebyshev.cpp b/src/2geom/chebyshev.cpp new file mode 100644 index 000000000..b56c57eac --- /dev/null +++ b/src/2geom/chebyshev.cpp @@ -0,0 +1,126 @@ +#include "chebyshev.h" + +#include "sbasis.h" +#include "sbasis-poly.h" + +#include +using std::vector; + +#include +#include + +namespace Geom{ + +SBasis cheb(unsigned n) { + static std::vector basis; + if(basis.empty()) { + basis.push_back(Linear(1,1)); + basis.push_back(Linear(0,1)); + } + for(unsigned i = basis.size(); i <= n; i++) { + basis.push_back(Linear(0,2)*basis[i-1] - basis[i-2]); + } + + return basis[n]; +} + +SBasis cheb_series(unsigned n, double* cheb_coeff) { + SBasis r; + for(unsigned i = 0; i < n; i++) { + double cof = cheb_coeff[i]; + //if(i == 0) + //cof /= 2; + r += cheb(i)*cof; + } + + return r; +} + +SBasis clenshaw_series(unsigned m, double* cheb_coeff) { + /** b_n = a_n + b_n-1 = 2*x*b_n + a_n-1 + b_n-k = 2*x*b_{n-k+1} + a_{n-k} - b_{n - k + 2} + b_0 = x*b_1 + a_0 - b_2 + */ + + double a = -1, b = 1; + SBasis d, dd; + SBasis y = (Linear(0, 2) - (a+b)) / (b-a); + SBasis y2 = 2*y; + for(int j = m - 1; j >= 1; j--) { + SBasis sv = d; + d = y2*d - dd + cheb_coeff[j]; + dd = sv; + } + + return y*d - dd + 0.5*cheb_coeff[0]; +} + +SBasis chebyshev_approximant (double (*f)(double,void*), int order, Interval in, void* p) { + gsl_cheb_series *cs = gsl_cheb_alloc (order+2); + + gsl_function F; + + F.function = f; + F.params = p; + + gsl_cheb_init (cs, &F, in[0], in[1]); + + SBasis r = compose(clenshaw_series(order, cs->c), Linear(-1,1)); + + gsl_cheb_free (cs); + return r; +} + +struct wrap { + double (*f)(double,void*); + void* pp; + double fa, fb; + Interval in; +}; + +double f_interp(double x, void* p) { + struct wrap *wr = (struct wrap *)p; + double z = (x - wr->in[0]) / (wr->in[1] - wr->in[0]); + return (wr->f)(x, wr->pp) - ((1 - z)*wr->fa + z*wr->fb); +} + +SBasis chebyshev_approximant_interpolating (double (*f)(double,void*), + int order, Interval in, void* p) { + double fa = f(in[0], p); + double fb = f(in[1], p); + struct wrap wr; + wr.fa = fa; + wr.fb = fb; + wr.in = in; + printf("%f %f\n", fa, fb); + wr.f = f; + wr.pp = p; + return compose(Linear(in[0], in[1]), Linear(fa, fb)) + chebyshev_approximant(f_interp, order, in, &wr) + Linear(fa, fb); +} + +SBasis chebyshev(unsigned n) { + static std::vector basis; + if(basis.empty()) { + basis.push_back(Linear(1,1)); + basis.push_back(Linear(0,1)); + } + for(unsigned i = basis.size(); i <= n; i++) { + basis.push_back(Linear(0,2)*basis[i-1] - basis[i-2]); + } + + return basis[n]; +} + +}; + +/* + Local Variables: + mode:c++ + c-file-style:"stroustrup" + c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) + indent-tabs-mode:nil + fill-column:99 + End: +*/ +// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 : diff --git a/src/2geom/chebyshev.h b/src/2geom/chebyshev.h new file mode 100644 index 000000000..f7e82a5df --- /dev/null +++ b/src/2geom/chebyshev.h @@ -0,0 +1,30 @@ +#ifndef _CHEBYSHEV +#define _CHEBYSHEV + +#include "sbasis.h" +#include "interval.h" + +/*** Conversion between Chebyshev approximation and SBasis. + * + */ + +namespace Geom{ + +SBasis chebyshev_approximant (double (*f)(double,void*), int order, Interval in, void* p=0); +SBasis chebyshev_approximant_interpolating (double (*f)(double,void*), int order, Interval in, void* p=0); +SBasis chebyshev(unsigned n); + +}; + +/* + Local Variables: + mode:c++ + c-file-style:"stroustrup" + c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) + indent-tabs-mode:nil + fill-column:99 + End: +*/ +// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 : + +#endif diff --git a/src/2geom/circulator.h b/src/2geom/circulator.h index 4c7f83bad..512aac39f 100644 --- a/src/2geom/circulator.h +++ b/src/2geom/circulator.h @@ -38,18 +38,18 @@ namespace Geom { template class Circulator { public: - typedef std::random_access_iterator_tag std::iterator_category; - typedef std::iterator_traits::value_type value_type; - typedef std::iterator_traits::difference_type difference_type; - typedef std::iterator_traits::pointer pointer; - typedef std::iterator_traits::reference reference; + typedef std::random_access_iterator_tag iterator_category; + typedef typename std::iterator_traits::value_type value_type; + typedef typename std::iterator_traits::difference_type difference_type; + typedef typename std::iterator_traits::pointer pointer; + typedef typename std::iterator_traits::reference reference; Circulator(Iterator const &first, Iterator const &last, Iterator const &pos) : _first(first), _last(last), _pos(pos) { - match_random_access(std::iterator_category(first)); + match_random_access(iterator_category(first)); } reference operator*() const { @@ -101,12 +101,12 @@ public: return _pos - other._pos; } - reference operator[n] const { + reference operator[](int n) const { return *_offset(n); } private: - void match_random_access(random_access_iterator_tag) {} + void match_random_access(iterator_category) {} Iterator _offset(int n) { difference_type range=( _last - _first ); diff --git a/src/2geom/concepts.h b/src/2geom/concepts.h index 42192c445..bd020ab66 100644 --- a/src/2geom/concepts.h +++ b/src/2geom/concepts.h @@ -34,7 +34,7 @@ #include "sbasis.h" #include "interval.h" #include "point.h" - +#include #include namespace Geom { @@ -66,6 +66,9 @@ struct FragmentConcept { bool b; BoundsType i; Interval dom; + std::vector v; + unsigned u; + SbType sb; void constraints() { t = T(o); b = t.isZero(); @@ -74,7 +77,10 @@ struct FragmentConcept { o = t.at1(); o = t.valueAt(d); o = t(d); - SbType sb = t.toSBasis(); + v = t.valueAndDerivatives(d, u); + //Is a pure derivative (ignoring others) accessor ever much faster? + //u = number of values returned. first val is value. + sb = t.toSBasis(); t = reverse(t); i = bounds_fast(t); i = bounds_exact(t); diff --git a/src/2geom/crossing.cpp b/src/2geom/crossing.cpp new file mode 100644 index 000000000..8f29635e8 --- /dev/null +++ b/src/2geom/crossing.cpp @@ -0,0 +1,81 @@ +#include "crossing.h" + +namespace Geom { + +void merge_crossings(Crossings &a, Crossings &b, unsigned i) { + Crossings n; + sort_crossings(b, i); + n.resize(a.size() + b.size()); + std::merge(a.begin(), a.end(), b.begin(), b.end(), n.begin(), CrossingOrder(i)); + a = n; +} + +void offset_crossings(Crossings &cr, double a, double b) { + for(unsigned i = 0; i < cr.size(); i++) { + cr[i].ta += a; + cr[i].tb += b; + } +} + +Crossings reverse_ta(Crossings const &cr, std::vector max) { + Crossings ret; + for(Crossings::const_iterator i = cr.begin(); i != cr.end(); ++i) { + double mx = max[i->a]; + ret.push_back(Crossing(i->ta > mx+0.01 ? (1 - (i->ta - mx) + mx) : mx - i->ta, + i->tb, !i->dir)); + } + return ret; +} + +Crossings reverse_tb(Crossings const &cr, unsigned split, std::vector max) { + Crossings ret; + for(Crossings::const_iterator i = cr.begin(); i != cr.end(); ++i) { + double mx = max[i->b - split]; + std::cout << i->b << "\n"; + ret.push_back(Crossing(i->ta, i->tb > mx+0.01 ? (1 - (i->tb - mx) + mx) : mx - i->tb, + !i->dir)); + } + return ret; +} + +CrossingSet reverse_ta(CrossingSet const &cr, unsigned split, std::vector max) { + CrossingSet ret; + for(unsigned i = 0; i < cr.size(); i++) { + Crossings res = reverse_ta(cr[i], max); + if(i < split) std::reverse(res.begin(), res.end()); + ret.push_back(res); + } + return ret; +} + +CrossingSet reverse_tb(CrossingSet const &cr, unsigned split, std::vector max) { + CrossingSet ret; + for(unsigned i = 0; i < cr.size(); i++) { + Crossings res = reverse_tb(cr[i], split, max); + if(i >= split) std::reverse(res.begin(), res.end()); + ret.push_back(res); + } + return ret; +} + +void clean(Crossings &cr_a, Crossings &cr_b) { +/* if(cr_a.empty()) return; + + //Remove anything with dupes + + for(Eraser i(&cr_a); !i.ended(); i++) { + const Crossing cur = *i; + Eraser next(i); + next++; + if(near(cur, *next)) { + cr_b.erase(std::find(cr_b.begin(), cr_b.end(), cur)); + for(i = next; near(*i, cur); i++) { + cr_b.erase(std::find(cr_b.begin(), cr_b.end(), *i)); + } + continue; + } + } +*/ +} + +} diff --git a/src/2geom/crossing.h b/src/2geom/crossing.h new file mode 100644 index 000000000..72b2eea4b --- /dev/null +++ b/src/2geom/crossing.h @@ -0,0 +1,103 @@ +#ifndef __GEOM_CROSSING_H +#define __GEOM_CROSSING_H + +#include +#include +#include "rect.h" +#include "sweep.h" +namespace Geom { + +struct Crossing { + bool dir; //True: along a, a becomes outside. + double ta, tb; //time on a and b of crossing + unsigned a, b; //storage of indices + Crossing() : dir(false), ta(0), tb(1), a(0), b(1) {} + Crossing(double t_a, double t_b, bool direction) : dir(direction), ta(t_a), tb(t_b), a(0), b(1) {} + Crossing(double t_a, double t_b, unsigned ai, unsigned bi, bool direction) : dir(direction), ta(t_a), tb(t_b), a(ai), b(bi) {} + bool operator==(const Crossing & other) const { return a == other.a && b == other.b && dir == other.dir && ta == other.ta && tb == other.tb; } + bool operator!=(const Crossing & other) const { return !(*this == other); } + unsigned getOther(unsigned cur) { return a == cur ? b : a; } +}; + + +/*inline bool near(Crossing a, Crossing b) { + return near(a.ta, b.ta) && near(a.tb, b.tb); +} + +struct NearF { bool operator()(Crossing a, Crossing b) { return near(a, b); } }; +*/ + +struct CrossingOrder { + unsigned ix; + CrossingOrder(unsigned i) : ix(i) {} + bool operator()(Crossing a, Crossing b) { + return (ix == a.a ? a.ta : a.tb) < + (ix == b.a ? b.ta : b.tb); + } +}; + +typedef std::vector Crossings; +typedef std::vector CrossingSet; + +template +std::vector bounds(C const &a) { + std::vector rs; + for(unsigned i = 0; i < a.size(); i++) rs.push_back(a[i].boundsFast()); + return rs; +} + +inline void sort_crossings(Crossings &cr, unsigned ix) { std::sort(cr.begin(), cr.end(), CrossingOrder(ix)); } + +template +struct Crosser { + virtual Crossings crossings(T const &a, T const &b) { return crossings(std::vector(1,a), std::vector(1,b))[0]; } + virtual CrossingSet crossings(std::vector const &a, std::vector const &b) { + CrossingSet results(a.size() + b.size(), Crossings()); + + std::vector > cull = sweep_bounds(bounds(a), bounds(b)); + for(unsigned i = 0; i < cull.size(); i++) { + for(unsigned jx = 0; jx < cull[i].size(); jx++) { + unsigned j = cull[i][jx]; + unsigned jc = j + a.size(); + Crossings cr = crossings(a[i], b[j]); + for(unsigned k = 0; k < cr.size(); k++) { cr[k].a = i; cr[k].b = jc; } + + //Sort & add A-sorted crossings + sort_crossings(cr, i); + Crossings n(results[i].size() + cr.size()); + std::merge(results[i].begin(), results[i].end(), cr.begin(), cr.end(), n.begin(), CrossingOrder(i)); + results[i] = n; + + //Sort & add B-sorted crossings + sort_crossings(cr, jc); + n.resize(results[jc].size() + cr.size()); + std::merge(results[jc].begin(), results[jc].end(), cr.begin(), cr.end(), n.begin(), CrossingOrder(jc)); + results[jc] = n; + } + } + return results; + } +}; +void merge_crossings(Crossings &a, Crossings &b, unsigned i); +void offset_crossings(Crossings &cr, double a, double b); + +Crossings reverse_ta(Crossings const &cr, std::vector max); +Crossings reverse_tb(Crossings const &cr, unsigned split, std::vector max); +CrossingSet reverse_ta(CrossingSet const &cr, unsigned split, std::vector max); +CrossingSet reverse_tb(CrossingSet const &cr, unsigned split, std::vector max); + +void clean(Crossings &cr_a, Crossings &cr_b); + +} + +#endif +/* + Local Variables: + mode:c++ + c-file-style:"stroustrup" + c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) + indent-tabs-mode:nil + fill-column:99 + End: +*/ +// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 : diff --git a/src/2geom/d2-sbasis.cpp b/src/2geom/d2-sbasis.cpp new file mode 100644 index 000000000..b68e3e6c5 --- /dev/null +++ b/src/2geom/d2-sbasis.cpp @@ -0,0 +1,148 @@ +#include "d2.h" +/* One would think that we would include d2-sbasis.h, however, + * you cannot actually include it in anything - only d2 may import it. + * This is due to the trickinesses of template submatching. */ + +namespace Geom { + +SBasis L2(D2 const & a, unsigned k) { return sqrt(dot(a, a), k); } + +D2 multiply(Linear const & a, D2 const & b) { + return D2(multiply(a, b[X]), multiply(a, b[Y])); +} + +D2 multiply(SBasis const & a, D2 const & b) { + return D2(multiply(a, b[X]), multiply(a, b[Y])); +} + +D2 truncate(D2 const & a, unsigned terms) { + return D2(truncate(a[X], terms), truncate(a[Y], terms)); +} + +unsigned sbasis_size(D2 const & a) { + return std::max((unsigned) a[0].size(), (unsigned) a[1].size()); +} + +//TODO: Is this sensical? shouldn't it be like pythagorean or something? +double tail_error(D2 const & a, unsigned tail) { + return std::max(a[0].tailError(tail), a[1].tailError(tail)); +} + +Piecewise > sectionize(D2 > const &a) { + Piecewise x = partition(a[0], a[1].cuts), y = partition(a[1], a[0].cuts); + assert(x.size() == y.size()); + Piecewise > ret; + for(unsigned i = 0; i < x.size(); i++) + ret.push_seg(D2(x[i], y[i])); + ret.cuts.insert(ret.cuts.end(), x.cuts.begin(), x.cuts.end()); + return ret; +} + +D2 > make_cuts_independant(Piecewise > const &a) { + D2 > ret; + for(unsigned d = 0; d < 2; d++) { + for(unsigned i = 0; i < a.size(); i++) + ret[d].push_seg(a[i][d]); + ret[d].cuts.insert(ret[d].cuts.end(), a.cuts.begin(), a.cuts.end()); + } + return ret; +} + +Piecewise > rot90(Piecewise > const &M){ + Piecewise > result; + if (M.empty()) return M; + result.push_cut(M.cuts[0]); + for (unsigned i=0; i dot(Piecewise > const &a, + Piecewise > const &b){ + Piecewise result; + if (a.empty() || b.empty()) return result; + Piecewise > aa = partition(a,b.cuts); + Piecewise > bb = partition(b,a.cuts); + + result.push_cut(aa.cuts.front()); + for (unsigned i=0; i cross(Piecewise > const &a, + Piecewise > const &b){ + Piecewise result; + if (a.empty() || b.empty()) return result; + Piecewise > aa = partition(a,b.cuts); + Piecewise > bb = partition(b,a.cuts); + + result.push_cut(aa.cuts.front()); + for (unsigned i=0; i > remove_short_cuts(Piecewise > const &f, double tol){ + double min = tol; + unsigned idx = f.size(); + for(unsigned i=0; i f.cuts[i+1]-f.cuts[i]){ + min = f.cuts[i+1]-f.cuts[i]; + idx = int(i); + } + } + if (idx==f.size()){ + return f; + } + if (f.size()==1) { + //removing this seg would result in an empty pw>... + return f; + } + Piecewise > new_f=f; + for (int dim=0; dim<2; dim++){ + double v = Hat(f.segs.at(idx)[dim][0]); + //TODO: what about closed curves? + if (idx>0 && f.segs.at(idx-1).at1()==f.segs.at(idx).at0()) + new_f.segs.at(idx-1)[dim][0][1] = v; + if (idx0, only force continuity where the jump is smaller than tol. +Piecewise > force_continuity(Piecewise > const &f, + double tol, + bool closed){ + if (f.size()==0) return f; + Piecewise > result=f; + unsigned cur = (closed)? 0:1; + unsigned prev = (closed)? f.size()-1:0; + while(cur compose(D2 const & a, SBasis const & b) { + return D2(compose(a[X], b), compose(a[Y], b)); +} + +SBasis L2(D2 const & a, unsigned k); +double L2(D2 const & a); + +D2 multiply(Linear const & a, D2 const & b); +inline D2 operator*(Linear const & a, D2 const & b) { return multiply(a, b); } +D2 multiply(SBasis const & a, D2 const & b); +inline D2 operator*(SBasis const & a, D2 const & b) { return multiply(a, b); } +D2 truncate(D2 const & a, unsigned terms); + +unsigned sbasis_size(D2 const & a); +double tail_error(D2 const & a, unsigned tail); + +//Piecewise > specific decls: + +Piecewise > sectionize(D2 > const &a); +D2 > make_cuts_independant(Piecewise > const &a); +Piecewise > rot90(Piecewise > const &a); +Piecewise dot(Piecewise > const &a, Piecewise > const &b); +Piecewise cross(Piecewise > const &a, Piecewise > const &b); + +Piecewise > force_continuity(Piecewise > const &f, + double tol=0, + bool closed=false); + +class CoordIterator +: public std::iterator +{ +public: + CoordIterator(std::vector >::const_iterator const &iter, unsigned d) : impl_(iter), ix_(d) {} + + inline bool operator==(CoordIterator const &other) { return other.impl_ == impl_; } + inline bool operator!=(CoordIterator const &other) { return other.impl_ != impl_; } + + inline SBasis operator*() const { + return (*impl_)[ix_]; + } + + inline CoordIterator &operator++() { + ++impl_; + return *this; + } + inline CoordIterator operator++(int) { + CoordIterator old=*this; + ++(*this); + return old; + } + +private: + std::vector >::const_iterator impl_; + unsigned ix_; +}; + +inline CoordIterator iterateCoord(Piecewise > const &a, unsigned d) { + return CoordIterator(a.segs.begin(), d); +} + +//bounds specializations with order +inline Rect bounds_fast(D2 const & s, unsigned order=0) { + return Rect(bounds_fast(s[X], order), + bounds_fast(s[Y], order)); +} +inline Rect bounds_local(D2 const & s, Interval i, unsigned order=0) { + return Rect(bounds_local(s[X], i, order), + bounds_local(s[Y], i, order)); +} + +} + +#endif +#endif diff --git a/src/2geom/d2.h b/src/2geom/d2.h index 82070519c..696ad0191 100644 --- a/src/2geom/d2.h +++ b/src/2geom/d2.h @@ -83,6 +83,15 @@ class D2{ Point valueAt(double t) const { boost::function_requires >(); return (*this)(t); + } + std::vector valueAndDerivatives(double t, unsigned count) const { + std::vector x = f[X].valueAndDerivatives(t, count), + y = f[Y].valueAndDerivatives(t, count); + std::vector res; + for(unsigned i = 0; i < count; i++) { + res.push_back(Point(x[i], y[i])); + } + return res; } D2 toSBasis() const { boost::function_requires >(); @@ -91,13 +100,18 @@ class D2{ Point operator()(double t) const; Point operator()(double x, double y) const; -}; - +}; template -D2 reverse(const D2 &a) { +inline D2 reverse(const D2 &a) { boost::function_requires >(); return D2(reverse(a[X]), reverse(a[Y])); } + +template +inline D2 portion(const D2 &a, Coord f, Coord t) { + boost::function_requires >(); + return D2(portion(a[X], f, t), portion(a[Y], f, t)); +} //IMPL: boost::EqualityComparableConcept template @@ -347,24 +361,16 @@ D2::operator()(double x, double y) const { template D2 derivative(D2 const & a) { return D2(derivative(a[X]), derivative(a[Y])); -} - +} template D2 integral(D2 const & a) { return D2(integral(a[X]), integral(a[Y])); -} +} } //end namespace Geom - - -//TODO: implement intersect - - #include "rect.h" -#include "sbasis.h" -#include "sbasis-2d.h" -#include "piecewise.h" +#include "d2-sbasis.h" namespace Geom{ @@ -383,83 +389,7 @@ template Rect bounds_local(const D2 &a, const Interval &t) { boost::function_requires >(); return Rect(bounds_local(a[X], t), bounds_local(a[Y], t)); -} - -//D2 specific decls: - -inline D2 compose(D2 const & a, SBasis const & b) { - return D2(compose(a[X], b), compose(a[Y], b)); -} - -SBasis L2(D2 const & a, unsigned k); -double L2(D2 const & a); - -inline D2 portion(D2 const &M, double t0, double t1){ - return D2(portion(M[0],t0,t1),portion(M[1],t0,t1)); -} - -D2 multiply(Linear const & a, D2 const & b); -inline D2 operator*(Linear const & a, D2 const & b) { return multiply(a, b); } -D2 multiply(SBasis const & a, D2 const & b); -inline D2 operator*(SBasis const & a, D2 const & b) { return multiply(a, b); } -D2 truncate(D2 const & a, unsigned terms); - -unsigned sbasis_size(D2 const & a); -double tail_error(D2 const & a, unsigned tail); - -//Piecewise > specific decls: - -Piecewise > sectionize(D2 > const &a); -D2 > make_cuts_independant(Piecewise > const &a); -Piecewise > rot90(Piecewise > const &a); -Piecewise dot(Piecewise > const &a, Piecewise > const &b); -Piecewise cross(Piecewise > const &a, Piecewise > const &b); - -Piecewise > force_continuity(Piecewise > const &f, - double tol=0, - bool closed=false); - -class CoordIterator -: public std::iterator -{ -public: - CoordIterator(std::vector >::const_iterator const &iter, unsigned d) : impl_(iter), ix_(d) {} - - inline bool operator==(CoordIterator const &other) { return other.impl_ == impl_; } - inline bool operator!=(CoordIterator const &other) { return other.impl_ != impl_; } - - inline SBasis operator*() const { - return (*impl_)[ix_]; - } - - inline CoordIterator &operator++() { - ++impl_; - return *this; - } - inline CoordIterator operator++(int) { - CoordIterator old=*this; - ++(*this); - return old; - } - -private: - std::vector >::const_iterator impl_; - unsigned ix_; -}; - -inline CoordIterator iterateCoord(Piecewise > const &a, unsigned d) { - return CoordIterator(a.segs.begin(), d); -} - -//bounds specializations with order -inline Rect bounds_fast(D2 const & s, unsigned order=0) { - return Rect(bounds_fast(s[X], order), - bounds_fast(s[Y], order)); -} -inline Rect bounds_local(D2 const & s, Interval i, unsigned order=0) { - return Rect(bounds_local(s[X], i, order), - bounds_local(s[Y], i, order)); -} +} }; /* diff --git a/src/2geom/interval.h b/src/2geom/interval.h index 459f2cd49..19e08978c 100644 --- a/src/2geom/interval.h +++ b/src/2geom/interval.h @@ -59,57 +59,50 @@ public: _b[0] = v; _b[1] = u; } } - + double operator[](unsigned i) const { assert(i < 2); return _b[i]; } - double& operator[](unsigned i) { return _b[i]; } //Trust the user... - - Coord min() const { return _b[0]; } - Coord max() const { return _b[1]; } - Coord extent() const { return _b[1] - _b[0]; } - Coord middle() const { return (_b[1] + _b[0]) * 0.5; } - - bool isEmpty() const { return _b[0] == _b[1]; } - bool contains(Coord val) const { return _b[0] <= val && val <= _b[1]; } + inline double& operator[](unsigned i) { return _b[i]; } //Trust the user... + + inline Coord min() const { return _b[0]; } + inline Coord max() const { return _b[1]; } + inline Coord extent() const { return _b[1] - _b[0]; } + inline Coord middle() const { return (_b[1] + _b[0]) * 0.5; } + + inline bool isEmpty() const { return _b[0] == _b[1]; } + inline bool contains(Coord val) const { return _b[0] <= val && val <= _b[1]; } bool contains(const Interval & val) const { return _b[0] <= val._b[0] && val._b[1] <= _b[1]; } bool intersects(const Interval & val) const { return contains(val._b[0]) || contains(val._b[1]) || val.contains(*this); } - - static Interval fromArray(const Coord* c, int n) { - assert(n > 0); - Interval result(c[0]); - for(int i = 1; i < n; i++) result.extendTo(c[i]); - return result; - } - - bool operator==(Interval other) { return _b[0] == other._b[0] && _b[1] == other._b[1]; } - 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) { return _b[0] != other._b[0] || _b[1] != other._b[1]; } + //IMPL: OffsetableConcept //TODO: rename output_type to something else in the concept typedef Coord output_type; - Interval operator+(Coord amnt) { + inline Interval operator+(Coord amnt) { return Interval(_b[0] + amnt, _b[1] + amnt); } - Interval operator-(Coord amnt) { + inline Interval operator-(Coord amnt) { return Interval(_b[0] - amnt, _b[1] - amnt); } - Interval operator+=(Coord amnt) { + inline Interval operator+=(Coord amnt) { _b[0] += amnt; _b[1] += amnt; return *this; } - Interval operator-=(Coord amnt) { + inline Interval operator-=(Coord amnt) { _b[0] -= amnt; _b[1] -= amnt; return *this; } - + //IMPL: ScalableConcept - Interval operator-() const { return Interval(*this); } - Interval operator*(Coord s) const { return Interval(_b[0]*s, _b[1]*s); } - Interval operator/(Coord s) const { return Interval(_b[0]/s, _b[1]/s); } + inline Interval operator-() const { return Interval(*this); } + inline Interval operator*(Coord s) const { return Interval(_b[0]*s, _b[1]*s); } + inline Interval operator/(Coord s) const { return Interval(_b[0]/s, _b[1]/s); } Interval operator*=(Coord s) { if(s < 0) { Coord temp = _b[0]; @@ -133,7 +126,7 @@ public: } return *this; } - + //TODO: NaN handleage for the next two? //TODO: Evaluate if wrap behaviour is proper. //If val > max, then rather than becoming a min==max range, it 'wraps' over @@ -154,18 +147,25 @@ public: _b[1] = val; } } - - void extendTo(Coord val) { + + inline void extendTo(Coord val) { if(val < _b[0]) _b[0] = val; if(val > _b[1]) _b[1] = val; //no else, as we want to handle NaN } - - void expandBy(double amnt) { + + static Interval fromArray(const Coord* c, int n) { + assert(n > 0); + Interval result(c[0]); + for(int i = 1; i < n; i++) result.extendTo(c[i]); + return result; + } + + inline void expandBy(double amnt) { _b[0] -= amnt; _b[1] += amnt; } - - void unionWith(const Interval & a) { + + inline void unionWith(const Interval & a) { if(a._b[0] < _b[0]) _b[0] = a._b[0]; if(a._b[1] > _b[1]) _b[1] = a._b[1]; } diff --git a/src/2geom/matrix.cpp b/src/2geom/matrix.cpp index f53795943..c0e64ad33 100644 --- a/src/2geom/matrix.cpp +++ b/src/2geom/matrix.cpp @@ -94,7 +94,7 @@ void Matrix::setExpansionY(double val) { double exp_y = expansionY(); if(!near(exp_y, 0.0)) { //TODO: best way to deal with it is to skip op? double coef = val / expansionY(); - for(unsigned i=2;i<4;i++) _c[i] *= coef; + for(unsigned i=2; i<4; i++) _c[i] *= coef; } } @@ -105,6 +105,8 @@ void Matrix::setIdentity() { _c[4] = 0.0; _c[5] = 0.0; } +//TODO: use eps + bool Matrix::isIdentity(Coord const eps) const { return near(_c[0], 1.0) && near(_c[1], 0.0) && near(_c[2], 0.0) && near(_c[3], 1.0) && @@ -151,6 +153,14 @@ bool Matrix::isRotation(Coord const eps) const { near(_c[0]*_c[0] + _c[1]*_c[1], 1.0); } +bool Matrix::onlyScaleAndTranslation(Coord const eps) const { + return near(_c[0], _c[3]) && near(_c[1], 0) && near(_c[2], 0); +} + +bool Matrix::flips() const { + return cross(xAxis(), yAxis()) > 0; +} + /** Returns the Scale/Rotate/skew part of the matrix without the translation part. */ Matrix Matrix::without_translation() const { return Matrix(_c[0], _c[1], _c[2], _c[3], 0, 0); diff --git a/src/2geom/matrix.h b/src/2geom/matrix.h index 6a45b50c8..c39c99716 100644 --- a/src/2geom/matrix.h +++ b/src/2geom/matrix.h @@ -90,6 +90,9 @@ class Matrix { bool isRotation(double eps = EPSILON) const; bool isScale(double eps = EPSILON) const; bool isUniformScale(double eps = EPSILON) const; + bool onlyScaleAndTranslation(double eps = EPSILON) const; + + bool flips() const; Matrix without_translation() const; diff --git a/src/2geom/ord.h b/src/2geom/ord.h new file mode 100644 index 000000000..d0e348aec --- /dev/null +++ b/src/2geom/ord.h @@ -0,0 +1,37 @@ + +#ifndef __2GEOM_ORD__ +#define __2GEOM_ORD__ + +namespace { + +enum Cmp { + LESS_THAN=-1, + GREATER_THAN=1, + EQUAL_TO=0 +}; + +inline Cmp operator-(Cmp x) { + switch(x) { + case LESS_THAN: + return GREATER_THAN; + case GREATER_THAN: + return LESS_THAN; + case EQUAL_TO: + return EQUAL_TO; + } +} + +template +inline Cmp cmp(T1 const &a, T2 const &b) { + if ( a < b ) { + return LESS_THAN; + } else if ( b < a ) { + return GREATER_THAN; + } else { + return EQUAL_TO; + } +} + +} + +#endif diff --git a/src/2geom/path-intersection.cpp b/src/2geom/path-intersection.cpp new file mode 100644 index 000000000..d641fcc08 --- /dev/null +++ b/src/2geom/path-intersection.cpp @@ -0,0 +1,588 @@ +#include "path-intersection.h" + +#include "ord.h" + +//for path_direction: +#include "sbasis-geometric.h" + +namespace Geom { + +/* This function computes the winding of the path, given a reference point. + * Positive values correspond to counter-clockwise in the mathematical coordinate system, + * and clockwise in screen coordinates. This particular implementation casts a ray in + * the positive x direction. It iterates the path, checking for intersection with the + * bounding boxes. If an intersection is found, the initial/final Y value of the curve is + * used to derive a delta on the winding value. If the point is within the bounding box, + * the curve specific winding function is called. + */ +int winding(Path const &path, Point p) { + //start on a segment which is not a horizontal line with y = p[y] + Path::const_iterator start; + for(Path::const_iterator iter = path.begin(); ; ++iter) { + if(iter == path.end_closed()) { return 0; } + if(iter->initialPoint()[Y]!=p[Y]) { start = iter; break; } + if(iter->finalPoint()[Y]!=p[Y]) { start = iter; break; } + if(iter->boundsFast().height()!=0.){ start = iter; break; } + } + int wind = 0; + int cnt = 0; + bool starting = true; + for (Path::const_iterator iter = start; iter != start || starting + ; ++iter, iter = (iter == path.end_closed()) ? path.begin() : iter ) + { + cnt++; + if(cnt > path.size()) return wind; //some bug makes this required + starting = false; + Rect bounds = iter->boundsFast(); + Coord x = p[X], y = p[Y]; + + if(x > bounds.right() || !bounds[Y].contains(y)) continue; //ray doesn't intersect box + + Point final = iter->finalPoint(); + Point initial = iter->initialPoint(); + Cmp final_to_ray = cmp(final[Y], y); + Cmp initial_to_ray = cmp(initial[Y], y); + + // if y is included, these will have opposite values, giving order. + Cmp c = cmp(final_to_ray, initial_to_ray); + if(x < bounds.left()) { + // ray goes through bbox + // winding delta determined by position of endpoints + if(final_to_ray != EQUAL_TO) { + wind += int(c); // GT = counter-clockwise = 1; LT = clockwise = -1; EQ = not-included = 0 + //std::cout << int(c) << " "; + goto cont; + } + } else { + //inside bbox, use custom per-curve winding thingie + int delt = iter->winding(p); + wind += delt; + //std::cout << "n" << delt << " "; + } + //Handling the special case of an endpoint on the ray: + if(final[Y] == y) { + //Traverse segments until it breaks away from y + //99.9% of the time this will happen the first go + Path::const_iterator next = iter; + next++; + for(; ; next++) { + if(next == path.end_closed()) next = path.begin(); + Rect bnds = next->boundsFast(); + //TODO: X considerations + if(bnds.height() > 0) { + //It has diverged + if(bnds.contains(p)) { + const double fudge = 0.01; + if(cmp(y, next->valueAt(fudge, Y)) == initial_to_ray) { + wind += int(c); + std::cout << "!!!!!" << int(c) << " "; + } + iter = next; // No increment, as the rest of the thing hasn't been counted. + } else { + Coord ny = next->initialPoint()[Y]; + if(cmp(y, ny) == initial_to_ray) { + //Is a continuation through the ray, so counts windingwise + wind += int(c); + std::cout << "!!!!!" << int(c) << " "; + } + iter = ++next; + } + goto cont; + } + if(next==start) return wind; + } + //Looks like it looped, which means everything's flat + return 0; + } + + cont:(void)0; + } + return wind; +} + +/* This function should only be applied to simple paths (regions), as otherwise + * a boolean winding direction is undefined. It returns true for fill, false for + * hole. Defaults to using the sign of area when it reaches funny cases. + */ +bool path_direction(Path const &p) { + //could probably be more efficient, but this is a quick job + double y = p.initialPoint()[Y]; + double x = p.initialPoint()[X]; + Cmp res = cmp(p[0].finalPoint()[Y], y); + goto doh; + for(unsigned i = 1; i <= p.size(); i++) { + Cmp final_to_ray = cmp(p[i].finalPoint()[Y], y); + Cmp initial_to_ray = cmp(p[i].initialPoint()[Y], y); + // if y is included, these will have opposite values, giving order. + Cmp c = cmp(final_to_ray, initial_to_ray); + if(c != EQUAL_TO) { + std::vector rs = p[i].roots(y, Y); + for(unsigned j = 0; j < rs.size(); j++) { + double nx = p[i].valueAt(rs[j], X); + if(nx > x) { + x = nx; + res = c; + } + } + } else if(final_to_ray == EQUAL_TO) goto doh; + } + return res < 0; + + doh: + //Otherwise fallback on area + + Piecewise > pw = p.toPwSb(); + double area; + Point centre; + Geom::centroid(pw, centre, area); + return area > 0; +} + +//pair intersect code based on njh's pair-intersect + +// A little sugar for appending a list to another +template +void append(T &a, T const &b) { + a.insert(a.end(), b.begin(), b.end()); +} + +/* Finds the intersection between the lines defined by A0 & A1, and B0 & B1. + * Returns through the last 3 parameters, returning the t-values on the lines + * and the cross-product of the deltas (a useful byproduct). The return value + * indicates if the time values are within their proper range on the line segments. + */ +bool +linear_intersect(Point A0, Point A1, Point B0, Point B1, + double &tA, double &tB, double &det) { + // kramers rule as cross products + Point Ad = A1 - A0, + Bd = B1 - B0, + d = B0 - A0; + det = cross(Ad, Bd); + if( 1.0 + det == 1.0 ) + return false; + else + { + double detinv = 1.0 / det; + tA = cross(d, Bd) * detinv; + tB = cross(d, Ad) * detinv; + return tA >= 0. && tA <= 1. && tB >= 0. && tB <= 1.; + } +} + +/* This uses the local bounds functions of curves to generically intersect two. + * It passes in the curves, time intervals, and keeps track of depth, while + * returning the results through the Crossings parameter. + */ +void pair_intersect(Curve const & A, double Al, double Ah, + Curve const & B, double Bl, double Bh, + Crossings &ret, unsigned depth=0) { + // std::cout << depth << "(" << Al << ", " << Ah << ")\n"; + Rect Ar = A.boundsLocal(Interval(Al, Ah)); + if(Ar.isEmpty()) return; + + Rect Br = B.boundsLocal(Interval(Bl, Bh)); + if(Br.isEmpty()) return; + + if(!Ar.intersects(Br)) return; + + //Checks the general linearity of the function + if((depth > 12) || (A.boundsLocal(Interval(Al, Ah), 1).maxExtent() < 0.1 + && B.boundsLocal(Interval(Bl, Bh), 1).maxExtent() < 0.1)) { + double tA, tB, c; + if(linear_intersect(A.pointAt(Al), A.pointAt(Ah), + B.pointAt(Bl), B.pointAt(Bh), + tA, tB, c)) { + tA = tA * (Ah - Al) + Al; + tB = tB * (Bh - Bl) + Bl; + if(depth % 2) + ret.push_back(Crossing(tB, tA, c < 0)); + else + ret.push_back(Crossing(tA, tB, c > 0)); + return; + } + } + if(depth > 12) return; + double mid = (Bl + Bh)/2; + pair_intersect(B, Bl, mid, + A, Al, Ah, + ret, depth+1); + pair_intersect(B, mid, Bh, + A, Al, Ah, + ret, depth+1); +} + +// A simple wrapper around pair_intersect +Crossings SimpleCrosser::crossings(Curve const &a, Curve const &b) { + Crossings ret; + pair_intersect(a, 0, 1, b, 0, 1, ret); + return ret; +} + +/* Takes two paths and time ranges on them, with the invariant that the + * paths are monotonic on the range. Splits A when the linear intersection + * doesn't exist or is inaccurate. Uses the fact that it is monotonic to + * do very fast local bounds. + */ +void mono_pair(Path const &A, double Al, double Ah, + Path const &B, double Bl, double Bh, + Crossings &ret, double tol, unsigned depth = 0) { + if( Al >= Ah || Bl >= Bh) return; + std::cout << " " << depth << "[" << Al << ", " << Ah << "]" << "[" << Bl << ", " << Bh << "]"; + + Point A0 = A.pointAt(Al), A1 = A.pointAt(Ah), + B0 = B.pointAt(Bl), B1 = B.pointAt(Bh); + //inline code that this implies? (without rect/interval construction) + if(!Rect(A0, A1).intersects(Rect(B0, B1)) || A0 == A1 || B0 == B1) return; + + //Checks the general linearity of the function + //if((depth > 12) || (A.boundsLocal(Interval(Al, Ah), 1).maxExtent() < 0.1 + // && B.boundsLocal(Interval(Bl, Bh), 1).maxExtent() < 0.1)) { + double tA, tB, c; + if(linear_intersect(A0, A1, B0, B1, + tA, tB, c)) { + tA = tA * (Ah - Al) + Al; + tB = tB * (Bh - Bl) + Bl; + if(depth % 2) + ret.push_back(Crossing(tB, tA, c < 0)); + else + ret.push_back(Crossing(tA, tB, c > 0)); + return; + } + //} + if(depth > 12) return; + double mid = (Bl + Bh)/2; + mono_pair(B, Bl, mid, + A, Al, Ah, + ret, depth+1); + mono_pair(B, mid, Bh, + A, Al, Ah, + ret, depth+1); +} + +// This returns the times when the x or y derivative is 0 in the curve. +std::vector curve_mono_splits(Curve const &d) { + std::vector rs = d.roots(0, X); + append(rs, d.roots(0, Y)); + std::sort(rs.begin(), rs.end()); + return rs; +} + +// Convenience function to add a value to each entry in a vector of doubles. +std::vector offset_doubles(std::vector const &x, double offs) { + std::vector ret; + for(unsigned i = 0; i < x.size(); i++) { + ret.push_back(x[i] + offs); + } + return ret; +} + +/* Finds all the monotonic splits for a path. Only includes the split between + * curves if they switch derivative directions at that point. + */ +std::vector path_mono_splits(Path const &p) { + std::vector ret; + if(p.empty()) return ret; + ret.push_back(0); + + Curve* deriv = p[0].derivative(); + append(ret, curve_mono_splits(*deriv)); + delete deriv; + + bool pdx=2, pdy=2; //Previous derivative direction + for(unsigned i = 0; i <= p.size(); i++) { + deriv = p[i].derivative(); + std::vector spl = offset_doubles(curve_mono_splits(*deriv), i); + delete deriv; + bool dx = p[i].initialPoint()[X] > (spl.empty()? p[i].finalPoint()[X] : + p.valueAt(spl.front(), X)); + bool dy = p[i].initialPoint()[Y] > (spl.empty()? p[i].finalPoint()[Y] : + p.valueAt(spl.front(), Y)); + //The direction changed, include the split time + if(dx != pdx || dy != pdy) { + ret.push_back(i); + pdx = dx; pdy = dy; + } + append(ret, spl); + } + return ret; +} + +/* Applies path_mono_splits to multiple paths, and returns the results such that + * time-set i corresponds to Path i. + */ +std::vector > paths_mono_splits(std::vector const &ps) { + std::vector > ret; + for(unsigned i = 0; i < ps.size(); i++) + ret.push_back(path_mono_splits(ps[i])); + return ret; +} + +/* Processes the bounds for a list of paths and a list of splits on them, yielding a list of rects for each. + * Each entry i corresponds to path i of the input. The number of rects in each entry is guaranteed to be the + * number of splits for that path, subtracted by one. + */ +std::vector > split_bounds(std::vector const &p, std::vector > splits) { + std::vector > ret; + for(unsigned i = 0; i < p.size(); i++) { + std::vector res; + for(unsigned j = 1; j < splits[i].size(); j++) + res.push_back(Rect(p[i].pointAt(splits[i][j-1]), p[i].pointAt(splits[i][j]))); + ret.push_back(res); + } + return ret; +} + +/* This is the main routine of "MonoCrosser", and implements a monotonic strategy on multiple curves. + * Finds crossings between two sets of paths, yielding a CrossingSet. [0, a.size()) of the return correspond + * to the sorted crossings of a with paths of b. The rest of the return, [a.size(), a.size() + b.size()], + * corresponds to the sorted crossings of b with paths of a. + * + * This function does two sweeps, one on the bounds of each path, and after that cull, one on the curves within. + * This leads to a certain amount of code complexity, however, most of that is factored into the above functions + */ +CrossingSet MonoCrosser::crossings(std::vector const &a, std::vector const &b) { + if(b.empty()) return CrossingSet(a.size(), Crossings()); + CrossingSet results(a.size() + b.size(), Crossings()); + if(a.empty()) return results; + + std::vector > splits_a = paths_mono_splits(a), splits_b = paths_mono_splits(b); + std::vector > bounds_a = split_bounds(a, splits_a), bounds_b = split_bounds(b, splits_b); + + std::vector bounds_a_union, bounds_b_union; + for(unsigned i = 0; i < bounds_a.size(); i++) bounds_a_union.push_back(union_list(bounds_a[i])); + for(unsigned i = 0; i < bounds_b.size(); i++) bounds_b_union.push_back(union_list(bounds_b[i])); + + std::vector > cull = sweep_bounds(bounds_a_union, bounds_b_union); + Crossings n; + for(unsigned i = 0; i < cull.size(); i++) { + for(unsigned jx = 0; jx < cull[i].size(); jx++) { + unsigned j = cull[i][jx]; + unsigned jc = j + a.size(); + Crossings res; + + //Sweep of the monotonic portions + std::vector > cull2 = sweep_bounds(bounds_a[i], bounds_b[j]); + for(unsigned k = 0; k < cull2.size(); k++) { + for(unsigned lx = 0; lx < cull2[k].size(); lx++) { + unsigned l = cull2[k][lx]; + mono_pair(a[i], splits_a[i][k-1], splits_a[i][k], + b[j], splits_b[j][l-1], splits_b[j][l], + res, .1); + } + } + + for(unsigned k = 0; k < res.size(); k++) { res[k].a = i; res[k].b = jc; } + + merge_crossings(results[i], res, i); + merge_crossings(results[i], res, jc); + } + } + + return results; +} + +/* This function is similar codewise to the MonoCrosser, the main difference is that it deals with + * only one set of paths and includes self intersection +CrossingSet crossings_among(std::vector const &p) { + CrossingSet results(p.size(), Crossings()); + if(p.empty()) return results; + + std::vector > splits = paths_mono_splits(p); + std::vector > prs = split_bounds(p, splits); + std::vector rs; + for(unsigned i = 0; i < prs.size(); i++) rs.push_back(union_list(prs[i])); + + std::vector > cull = sweep_bounds(rs); + + //we actually want to do the self-intersections, so add em in: + for(unsigned i = 0; i < cull.size(); i++) cull[i].push_back(i); + + for(unsigned i = 0; i < cull.size(); i++) { + for(unsigned jx = 0; jx < cull[i].size(); jx++) { + unsigned j = cull[i][jx]; + Crossings res; + + //Sweep of the monotonic portions + std::vector > cull2 = sweep_bounds(prs[i], prs[j]); + for(unsigned k = 0; k < cull2.size(); k++) { + for(unsigned lx = 0; lx < cull2[k].size(); lx++) { + unsigned l = cull2[k][lx]; + mono_pair(p[i], splits[i][k-1], splits[i][k], + p[j], splits[j][l-1], splits[j][l], + res, .1); + } + } + + for(unsigned k = 0; k < res.size(); k++) { res[k].a = i; res[k].b = j; } + + merge_crossings(results[i], res, i); + merge_crossings(results[j], res, j); + } + } + + return results; +} +*/ + +Crossings curve_self_crossings(Curve const &a) { + Crossings res; + std::vector spl; + spl.push_back(0); + append(spl, curve_mono_splits(a)); + spl.push_back(1); + for(unsigned i = 1; i < spl.size(); i++) + for(unsigned j = i+1; j < spl.size(); j++) + pair_intersect(a, spl[i-1], spl[i], a, spl[j-1], spl[j], res); + return res; +} + +/* +void mono_curve_intersect(Curve const & A, double Al, double Ah, + Curve const & B, double Bl, double Bh, + Crossings &ret, unsigned depth=0) { + // std::cout << depth << "(" << Al << ", " << Ah << ")\n"; + Point A0 = A.pointAt(Al), A1 = A.pointAt(Ah), + B0 = B.pointAt(Bl), B1 = B.pointAt(Bh); + //inline code that this implies? (without rect/interval construction) + if(!Rect(A0, A1).intersects(Rect(B0, B1)) || A0 == A1 || B0 == B1) return; + + //Checks the general linearity of the function + if((depth > 12) || (A.boundsLocal(Interval(Al, Ah), 1).maxExtent() < 0.1 + && B.boundsLocal(Interval(Bl, Bh), 1).maxExtent() < 0.1)) { + double tA, tB, c; + if(linear_intersect(A0, A1, B0, B1, tA, tB, c)) { + tA = tA * (Ah - Al) + Al; + tB = tB * (Bh - Bl) + Bl; + if(depth % 2) + ret.push_back(Crossing(tB, tA, c < 0)); + else + ret.push_back(Crossing(tA, tB, c > 0)); + return; + } + } + if(depth > 12) return; + double mid = (Bl + Bh)/2; + mono_curve_intersect(B, Bl, mid, + A, Al, Ah, + ret, depth+1); + mono_curve_intersect(B, mid, Bh, + A, Al, Ah, + ret, depth+1); +} + +std::vector > curves_mono_splits(Path const &p) { + std::vector > ret; + for(unsigned i = 0; i <= p.size(); i++) { + std::vector spl; + spl.push_back(0); + append(spl, curve_mono_splits(p[i])); + spl.push_back(1); + ret.push_back(spl); + } +} + +std::vector > curves_split_bounds(Path const &p, std::vector > splits) { + std::vector > ret; + for(unsigned i = 0; i < splits.size(); i++) { + std::vector res; + for(unsigned j = 1; j < splits[i].size(); j++) + res.push_back(Rect(p.pointAt(splits[i][j-1]+i), p.pointAt(splits[i][j]+i))); + ret.push_back(res); + } + return ret; +} + + +Crossings path_self_crossings(Path const &p) { + Crossings ret; + std::vector > cull = sweep_bounds(bounds(p)); + std::vector > spl = curves_mono_splits(p); + std::vector > bnds = curves_split_bounds(p, spl); + for(unsigned i = 0; i < cull.size(); i++) { + Crossings res; + for(unsigned k = 1; k < spl[i].size(); k++) + for(unsigned l = k+1; l < spl[i].size(); l++) + mono_curve_intersect(p[i], spl[i][k-1], spl[i][k], p[i], spl[i][l-1], spl[i][l], res); + offset_crossings(res, i, i); + append(ret, res); + for(unsigned jx = 0; jx < cull[i].size(); jx++) { + unsigned j = cull[i][jx]; + res.clear(); + + std::vector > cull2 = sweep_bounds(bnds[i], bnds[j]); + for(unsigned k = 0; k < cull2.size(); k++) { + for(unsigned lx = 0; lx < cull2[k].size(); lx++) { + unsigned l = cull2[k][lx]; + mono_curve_intersect(p[i], spl[i][k-1], spl[i][k], p[j], spl[j][l-1], spl[j][l], res); + } + } + + //if(fabs(int(i)-j) == 1 || fabs(int(i)-j) == p.size()-1) { + Crossings res2; + for(unsigned k = 0; k < res.size(); k++) { + if(res[k].ta != 0 && res[k].ta != 1 && res[k].tb != 0 && res[k].tb != 1) { + res.push_back(res[k]); + } + } + res = res2; + //} + offset_crossings(res, i, j); + append(ret, res); + } + } + return ret; +} +*/ + +Crossings path_self_crossings(Path const &p) { + Crossings ret; + std::vector > cull = sweep_bounds(bounds(p)); + for(unsigned i = 0; i < cull.size(); i++) { + Crossings res = curve_self_crossings(p[i]); + offset_crossings(res, i, i); + append(ret, res); + for(unsigned jx = 0; jx < cull[i].size(); jx++) { + unsigned j = cull[i][jx]; + res.clear(); + pair_intersect(p[i], 0, 1, p[j], 0, 1, res); + + //if(fabs(int(i)-j) == 1 || fabs(int(i)-j) == p.size()-1) { + Crossings res2; + for(unsigned k = 0; k < res.size(); k++) { + if(res[k].ta != 0 && res[k].ta != 1 && res[k].tb != 0 && res[k].tb != 1) { + res.push_back(res[k]); + } + } + res = res2; + //} + offset_crossings(res, i, j); + append(ret, res); + } + } + return ret; +} + +CrossingSet crossings_among(std::vector const &p) { + CrossingSet results(p.size(), Crossings()); + if(p.empty()) return results; + + SimpleCrosser cc; + + std::vector > cull = sweep_bounds(bounds(p)); + for(unsigned i = 0; i < cull.size(); i++) { + Crossings res = path_self_crossings(p[i]); + for(unsigned k = 0; k < res.size(); k++) { res[k].a = res[k].b = i; } + merge_crossings(results[i], res, i); + for(unsigned jx = 0; jx < cull[i].size(); jx++) { + unsigned j = cull[i][jx]; + + Crossings res = cc.crossings(p[i], p[j]); + for(unsigned k = 0; k < res.size(); k++) { res[k].a = i; res[k].b = j; } + merge_crossings(results[i], res, i); + merge_crossings(results[j], res, j); + } + } +} + +} diff --git a/src/2geom/path-intersection.h b/src/2geom/path-intersection.h new file mode 100644 index 000000000..3401386e0 --- /dev/null +++ b/src/2geom/path-intersection.h @@ -0,0 +1,65 @@ +#ifndef __GEOM_PATH_INTERSECTION_H +#define __GEOM_PATH_INTERSECTION_H + +#include "path.h" + +#include "crossing.h" + +#include "sweep.h" + +namespace Geom { + +int winding(Path const &path, Point p); +bool path_direction(Path const &p); + +inline bool contains(Path const & p, Point i, bool evenodd = true) { + return (evenodd ? winding(p, i) % 2 : winding(p, i)) != 0; +} + +template +Crossings curve_sweep(Path const &a, Path const &b) { + T t; + Crossings ret; + std::vector bounds_a = bounds(a), bounds_b = bounds(b); + std::vector > ixs = sweep_bounds(bounds_a, bounds_b); + for(unsigned i = 0; i < a.size(); i++) { + for(std::vector::iterator jp = ixs[i].begin(); jp != ixs[i].end(); jp++) { + Crossings cc = t.crossings(a[i], b[*jp]); + offset_crossings(cc, i, *jp); + ret.insert(ret.end(), cc.begin(), cc.end()); + } + } + return ret; +} + +struct SimpleCrosser : public Crosser { + Crossings crossings(Curve const &a, Curve const &b); + Crossings crossings(Path const &a, Path const &b) { return curve_sweep(a, b); } + CrossingSet crossings(std::vector const &a, std::vector const &b) { return Crosser::crossings(a, b); } +}; + +struct MonoCrosser : public Crosser { + Crossings crossings(Path const &a, Path const &b) { return crossings(std::vector(1,a), std::vector(1,b))[0]; } + CrossingSet crossings(std::vector const &a, std::vector const &b); +}; + +typedef SimpleCrosser DefaultCrosser; + +std::vector path_mono_splits(Path const &p); + +CrossingSet crossings_among(std::vector const & p); +inline Crossings self_crossings(Path const & a) { return crossings_among(std::vector(1, a))[0]; } + +inline Crossings crossings(Path const & a, Path const & b) { + DefaultCrosser c = DefaultCrosser(); + return c.crossings(a, b); +} + +inline CrossingSet crossings(std::vector const & a, std::vector const & b) { + DefaultCrosser c = DefaultCrosser(); + return c.crossings(a, b); +} + +} + +#endif diff --git a/src/2geom/path.cpp b/src/2geom/path.cpp index 91868eb7e..f05d3b0cf 100644 --- a/src/2geom/path.cpp +++ b/src/2geom/path.cpp @@ -1,367 +1,258 @@ -/* - * Path - Series of continuous curves - * - * Copyright 2007 MenTaLguY - * - * 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 "path.h" - -namespace Geom { - -namespace { - -enum Cmp { - LESS_THAN=-1, - GREATER_THAN=1, - EQUAL_TO=0 -}; - -template -inline Cmp cmp(T1 const &a, T2 const &b) { - if ( a < b ) { - return LESS_THAN; - } else if ( b < a ) { - return GREATER_THAN; - } else { - return EQUAL_TO; - } -} - -} - -boost::optional CurveHelpers::sbasis_winding(D2 const &sb, Point p) { - Interval ix = bounds_fast(sb[X]); - - if ( p[X] > ix.max() ) { /* ray does not intersect bbox */ - return 0; - } - - SBasis fy = sb[Y]; - fy -= p[Y]; - - if (fy.empty()) { /* coincident horizontal segment */ - return boost::optional(); - } - - if ( p[X] < ix.min() ) { /* ray does not originate in bbox */ - double y = p[Y]; - /* winding determined by position of endpoints */ - Cmp initial_to_ray = cmp(fy[0][0], y); - Cmp final_to_ray = cmp(fy[0][1], y); - switch (cmp(final_to_ray, initial_to_ray)) { - case GREATER_THAN: - /* exclude final endpoint */ - return ( final_to_ray != EQUAL_TO ); - case LESS_THAN: - /* exclude final endpoint */ - return -( final_to_ray != EQUAL_TO ); - default: - /* any intersections cancel out */ - return 0; - } - } else { /* ray originates in bbox */ - std::vector ts = roots(fy); - - static const unsigned MAX_DERIVATIVES=8; - boost::optional ds[MAX_DERIVATIVES]; - ds[0] = derivative(fy); - - /* winding determined by summing signs of derivatives at intersections */ - int winding=0; - for ( std::vector::iterator ti = ts.begin() - ; ti != ts.end() - ; ++ti ) - { - double t = *ti; - if ( sb[X](t) >= p[X] ) { /* root is ray intersection */ - for ( boost::optional *di = ds - ; di != ( ds + MAX_DERIVATIVES ) - ; ++di ) - { - if (!*di) { - *di = derivative(**(di-1)); - } - switch (cmp((**di)(t), 0)) { - case GREATER_THAN: - if ( t < 1 ) { /* exclude final endpoint */ - winding += 1; - } - goto next_root; - case LESS_THAN: - if ( t < 1 ) { /* exclude final endpoint */ - winding -= 1; - } - goto next_root; - default: (void)0; - /* give up */ - }; - } - } -next_root: (void)0; - } - - return winding; - } -} - -Rect BezierHelpers::bounds(unsigned degree, Point const *points) { - Point min=points[0]; - Point max=points[0]; - for ( unsigned i = 1 ; i <= degree ; ++i ) { - for ( unsigned axis = 0 ; axis < 2 ; ++axis ) { - min[axis] = std::min(min[axis], points[i][axis]); - max[axis] = std::max(max[axis], points[i][axis]); - } - } - return Rect(min, max); -} - -Point BezierHelpers::point_and_derivatives_at(Coord t, - unsigned degree, - Point const *points, - unsigned n_derivs, - Point *derivs) -{ - return Point(0,0); // TODO -} - -Geom::Point -BezierHelpers::subdivideArr(Coord t, // Parameter value - unsigned degree, // Degree of bezier curve - Geom::Point const *V, // Control pts - Geom::Point *Left, // RETURN left half ctl pts - Geom::Point *Right) // RETURN right half ctl pts -{ - Geom::Point Vtemp[degree+1][degree+1]; - - /* Copy control points */ - std::copy(V, V+degree+1, Vtemp[0]); - - /* Triangle computation */ - for (unsigned i = 1; i <= degree; i++) { - for (unsigned j = 0; j <= degree - i; j++) { - Vtemp[i][j] = lerp(t, Vtemp[i-1][j], Vtemp[i-1][j+1]); - } - } - - for (unsigned j = 0; j <= degree; j++) - Left[j] = Vtemp[j][0]; - for (unsigned j = 0; j <= degree; j++) - Right[j] = Vtemp[degree-j][j]; - - return (Vtemp[degree][0]); -} - -void Path::swap(Path &other) { - std::swap(curves_, other.curves_); - std::swap(closed_, other.closed_); - std::swap(*final_, *other.final_); - curves_[curves_.size()-1] = final_; - other.curves_[other.curves_.size()-1] = other.final_; -} - -Rect Path::bounds_fast() const { - Rect bounds=front().bounds_fast(); - for ( const_iterator iter=++begin(); iter != end() ; ++iter ) { - bounds.unionWith(iter->bounds_fast()); - } - return bounds; -} - -Rect Path::bounds_exact() const { - Rect bounds=front().bounds_exact(); - for ( const_iterator iter=++begin(); iter != end() ; ++iter ) { - bounds.unionWith(iter->bounds_exact()); - } - return bounds; -} - -int Path::winding(Point p) const { - int winding = 0; - boost::optional ignore = boost::optional(); - for ( const_iterator iter = begin() - ; iter != end_closed() - ; ++iter ) - { - boost::optional w = iter->winding(p); - if (w) { - winding += *w; - ignore = boost::optional(); - } else { - Point initial = iter->initialPoint(); - Point final = iter->finalPoint(); - switch (cmp(initial[X], final[X])) { - case GREATER_THAN: - if ( !ignore || *ignore != GREATER_THAN ) { /* ignore repeated */ - winding += 1; - ignore = GREATER_THAN; - } - break; - case LESS_THAN: - if ( !ignore || *ignore != LESS_THAN ) { /* ignore repeated */ - if ( p[X] < final[X] ) { /* ignore final point */ - winding -= 1; - ignore = LESS_THAN; - } - } - break; - case EQUAL_TO: - /* always ignore null segments */ - break; - } - } - } - return winding; -} - -void Path::append(Curve const &curve) { - if ( curves_.front() != final_ && curve.initialPoint() != (*final_)[0] ) { - throw ContinuityError(); - } - do_append(curve.duplicate()); -} - -void Path::append(D2 const &curve) { - if ( curves_.front() != final_ ) { - for ( int i = 0 ; i < 2 ; ++i ) { - if ( curve[i][0][0] != (*final_)[0][i] ) { - throw ContinuityError(); - } - } - } - do_append(new SBasisCurve(curve)); -} - -void Path::do_update(Sequence::iterator first_replaced, - Sequence::iterator last_replaced, - Sequence::iterator first, - Sequence::iterator last) -{ - // note: modifies the contents of [first,last) - - check_continuity(first_replaced, last_replaced, first, last); - delete_range(first_replaced, last_replaced); - if ( ( last - first ) == ( last_replaced - first_replaced ) ) { - std::copy(first, last, first_replaced); - } else { - // this approach depends on std::vector's behavior WRT iterator stability - curves_.erase(first_replaced, last_replaced); - curves_.insert(first_replaced, first, last); - } - - if ( curves_.front() != final_ ) { - (*final_)[0] = back().finalPoint(); - (*final_)[1] = front().initialPoint(); - } -} - -void Path::do_append(Curve *curve) { - if ( curves_.front() == final_ ) { - (*final_)[1] = curve->initialPoint(); - } - curves_.insert(curves_.end()-1, curve); - (*final_)[0] = curve->finalPoint(); -} - -void Path::delete_range(Sequence::iterator first, Sequence::iterator last) { - for ( Sequence::iterator iter=first ; iter != last ; ++iter ) { - delete *iter; - } -} - -void Path::check_continuity(Sequence::iterator first_replaced, - Sequence::iterator last_replaced, - Sequence::iterator first, - Sequence::iterator last) -{ - if ( first != last ) { - if ( first_replaced != curves_.begin() ) { - if ( (*first_replaced)->initialPoint() != (*first)->initialPoint() ) { - throw ContinuityError(); - } - } - if ( last_replaced != (curves_.end()-1) ) { - if ( (*(last_replaced-1))->finalPoint() != (*(last-1))->finalPoint() ) { - throw ContinuityError(); - } - } - } else if ( first_replaced != last_replaced && first_replaced != curves_.begin() && last_replaced != curves_.end()-1) { - if ( (*first_replaced)->initialPoint() != - (*(last_replaced-1))->finalPoint() ) - { - throw ContinuityError(); - } - } -} - -Rect SBasisCurve::bounds_fast() const { - throw NotImplemented(); - return Rect(Point(0,0), Point(0,0)); -} - -Rect SBasisCurve::bounds_exact() const { - throw NotImplemented(); - return Rect(Point(0,0), Point(0,0)); -} - -Point SBasisCurve::pointAndDerivativesAt(Coord t, unsigned n_derivs, Point *derivs) const { - throw NotImplemented(); - return Point(0,0); -} - -Path const &SBasisCurve::subdivide(Coord t, Path &out) const { - throw NotImplemented(); -} - -Rect SVGEllipticalArc::bounds_fast() const { - throw NotImplemented(); -} -Rect SVGEllipticalArc::bounds_exact() const { - throw NotImplemented(); -} - -Point SVGEllipticalArc::pointAndDerivativesAt(Coord t, unsigned n_derivs, Point *derivs) const { - throw NotImplemented(); -} - -D2 SVGEllipticalArc::sbasis() const { - throw NotImplemented(); -} - -} - -/* - Local Variables: - mode:c++ - c-file-style:"stroustrup" - c-file-offsets:((innamespace . 0)(substatement-open . 0)) - indent-tabs-mode:nil - c-brace-offset:0 - fill-column:99 - End: - vim: filetype=cpp:expandtab:shiftwidth=2:tabstop=8:softtabstop=2 : -*/ - +/* + * Path - Series of continuous curves + * + * Copyright 2007 MenTaLguY + * + * 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 "path.h" + +#include "ord.h" + +namespace Geom { + +int CurveHelpers::root_winding(Curve const &c, Point p) { + std::vector ts = c.roots(p[Y], Y); + + if(ts.empty()) return 0; + + double const fudge = 0.01; //fudge factor used on first and last + + std::sort(ts.begin(), ts.end()); + + // winding determined by crossings at roots + int wind=0; + // previous time + double pt = ts.front() - fudge; + for ( std::vector::iterator ti = ts.begin() + ; ti != ts.end() + ; ++ti ) + { + double t = *ti; + if ( t <= 0. || t >= 1. ) continue; //skip endpoint roots + if ( c.valueAt(t, X) > p[X] ) { // root is ray intersection + // Get t of next: + std::vector::iterator next = ti; + next++; + double nt; + if(next == ts.end()) nt = t + fudge; else nt = *next; + + // Check before in time and after in time for positions + // Currently we're using the average times between next and previous segs + Cmp after_to_ray = cmp(c.valueAt((t + nt) / 2, Y), p[Y]); + Cmp before_to_ray = cmp(c.valueAt((t + pt) / 2, Y), p[Y]); + // if y is included, these will have opposite values, giving order. + Cmp dt = cmp(after_to_ray, before_to_ray); + if(dt != EQUAL_TO) //Should always be true, but yah never know.. + wind += dt; + pt = t; + } + } + + return wind; +} + +void Path::swap(Path &other) { + std::swap(curves_, other.curves_); + std::swap(closed_, other.closed_); + std::swap(*final_, *other.final_); + curves_[curves_.size()-1] = final_; + other.curves_[other.curves_.size()-1] = other.final_; +} + +Rect Path::boundsFast() const { + Rect bounds=front().boundsFast(); + for ( const_iterator iter=++begin(); iter != end() ; ++iter ) { + bounds.unionWith(iter->boundsFast()); + } + return bounds; +} + +Rect Path::boundsExact() const { + Rect bounds=front().boundsExact(); + for ( const_iterator iter=++begin(); iter != end() ; ++iter ) { + bounds.unionWith(iter->boundsExact()); + } + return bounds; +} + +template +iter inc(iter const &x, unsigned n) { + iter ret = x; + for(unsigned i = 0; i < n; i++) + ret++; + return ret; +} + +//This assumes that you can't be perfect in your t-vals, and as such, tweaks the start +void Path::appendPortionTo(Path &ret, double from, double to) const { + assert(from >= 0 && to >= 0); + if(to == 0) to = size()+0.999999; + if(from == to) { return; } + double fi, ti; + double ff = modf(from, &fi), tf = modf(to, &ti); + if(tf == 0) { ti--; tf = 1; } + const_iterator fromi = inc(begin(), (unsigned)fi); + if(fi == ti && from < to) { + Curve *v = fromi->portion(ff, tf); + ret.append(*v); + delete v; + return; + } + const_iterator toi = inc(begin(), (unsigned)ti); + if(ff != 1.) { + Curve *fromv = fromi->portion(ff, 1.); + //fromv->setInitial(ret.finalPoint()); + ret.append(*fromv); + delete fromv; + } + if(from > to) { + const_iterator ender = end(); + if(ender->initialPoint() == ender->finalPoint()) ender++; + ret.insert(ret.end(), ++fromi, ender); + ret.insert(ret.end(), begin(), toi); + } else { + ret.insert(ret.end(), ++fromi, toi); + } + Curve *tov = toi->portion(0., tf); + ret.append(*tov); + delete tov; +} + +const double eps = .1; + +void Path::append(Curve const &curve) { + if ( curves_.front() != final_ && !near(curve.initialPoint(), (*final_)[0], eps) ) { + throw ContinuityError(); + } + do_append(curve.duplicate()); +} + +void Path::append(D2 const &curve) { + if ( curves_.front() != final_ ) { + for ( int i = 0 ; i < 2 ; ++i ) { + if ( !near(curve[i][0][0], (*final_)[0][i], eps) ) { + throw ContinuityError(); + } + } + } + do_append(new SBasisCurve(curve)); +} + +void Path::do_update(Sequence::iterator first_replaced, + Sequence::iterator last_replaced, + Sequence::iterator first, + Sequence::iterator last) +{ + // note: modifies the contents of [first,last) + + check_continuity(first_replaced, last_replaced, first, last); + delete_range(first_replaced, last_replaced); + if ( ( last - first ) == ( last_replaced - first_replaced ) ) { + std::copy(first, last, first_replaced); + } else { + // this approach depends on std::vector's behavior WRT iterator stability + curves_.erase(first_replaced, last_replaced); + curves_.insert(first_replaced, first, last); + } + + if ( curves_.front() != final_ ) { + final_->setPoint(0, back().finalPoint()); + final_->setPoint(1, front().initialPoint()); + } +} + +void Path::do_append(Curve *curve) { + if ( curves_.front() == final_ ) { + final_->setPoint(1, curve->initialPoint()); + } + curves_.insert(curves_.end()-1, curve); + final_->setPoint(0, curve->finalPoint()); +} + +void Path::delete_range(Sequence::iterator first, Sequence::iterator last) { + for ( Sequence::iterator iter=first ; iter != last ; ++iter ) { + delete *iter; + } +} + +void Path::check_continuity(Sequence::iterator first_replaced, + Sequence::iterator last_replaced, + Sequence::iterator first, + Sequence::iterator last) +{ + if ( first != last ) { + if ( first_replaced != curves_.begin() ) { + if ( !near( (*first_replaced)->initialPoint(), (*first)->initialPoint(), eps ) ) { + throw ContinuityError(); + } + } + if ( last_replaced != (curves_.end()-1) ) { + if ( !near( (*(last_replaced-1))->finalPoint(), (*(last-1))->finalPoint(), eps ) ) { + throw ContinuityError(); + } + } + } else if ( first_replaced != last_replaced && first_replaced != curves_.begin() && last_replaced != curves_.end()-1) { + if ( !near((*first_replaced)->initialPoint(), (*(last_replaced-1))->finalPoint(), eps ) ) { + throw ContinuityError(); + } + } +} + +Rect SVGEllipticalArc::boundsFast() const { + throw NotImplemented(); +} +Rect SVGEllipticalArc::boundsExact() const { + throw NotImplemented(); +} +Rect SVGEllipticalArc::boundsLocal(Interval i, unsigned deg) const { + //throw NotImplemented(); +} + +std::vector SVGEllipticalArc::pointAndDerivatives(Coord t, unsigned n) const { + throw NotImplemented(); +} + +std::vector SVGEllipticalArc::roots(double v, Dim2 d) const { + //throw NotImplemented(); +} + +D2 SVGEllipticalArc::toSBasis() const { + return D2(Linear(initial_[X], final_[X]), Linear(initial_[Y], final_[Y])); +} + +} + +/* + Local Variables: + mode:c++ + c-file-style:"stroustrup" + c-file-offsets:((innamespace . 0)(substatement-open . 0)) + indent-tabs-mode:nil + c-brace-offset:0 + fill-column:99 + End: + vim: filetype=cpp:expandtab:shiftwidth=2:tabstop=8:softtabstop=2 : +*/ diff --git a/src/2geom/path.h b/src/2geom/path.h index 31a7173b7..6f39eb7ef 100644 --- a/src/2geom/path.h +++ b/src/2geom/path.h @@ -35,15 +35,21 @@ #include #include #include -#include #include "d2.h" -#include "bezier-to-sbasis.h" +#include "matrix.h" +#include "bezier.h" +#include "crossing.h" namespace Geom { -class Path; +class Curve; -class Curve { +struct CurveHelpers { +protected: + static int root_winding(Curve const &c, Point p); +}; + +class Curve : private CurveHelpers { public: virtual ~Curve() {} @@ -52,112 +58,211 @@ public: virtual Curve *duplicate() const = 0; - virtual Rect bounds_fast() const = 0; - virtual Rect bounds_exact() const = 0; - - virtual boost::optional winding(Point p) const = 0; + virtual Rect boundsFast() const = 0; + virtual Rect boundsExact() const = 0; + virtual Rect boundsLocal(Interval i, unsigned deg) const = 0; + Rect boundsLocal(Interval i) const { return boundsLocal(i, 0); } + + virtual std::vector roots(double v, Dim2 d) const = 0; - virtual Path const &subdivide(Coord t, Path &out) const = 0; + virtual int winding(Point p) const { return root_winding(*this, p); } - Point pointAt(Coord t) const { return pointAndDerivativesAt(t, 0, NULL); } - virtual Point pointAndDerivativesAt(Coord t, unsigned n, Point *ds) const = 0; - virtual D2 sbasis() const = 0; + //mental: review these + virtual Curve *portion(double f, double t) const = 0; + virtual Curve *reverse() const { return portion(1, 0); } + virtual Curve *derivative() const = 0; + + virtual void setInitial(Point v) = 0; + virtual void setFinal(Point v) = 0; + + virtual Curve *transformed(Matrix const &m) const = 0; + + virtual Point pointAt(Coord t) const { return pointAndDerivatives(t, 1).front(); } + virtual Coord valueAt(Coord t, Dim2 d) const { return pointAt(t)[d]; } + virtual std::vector pointAndDerivatives(Coord t, unsigned n) const = 0; + virtual D2 toSBasis() const = 0; }; -struct CurveHelpers { -protected: - static boost::optional sbasis_winding(D2 const &sbasis, Point p); -}; +class SBasisCurve : public Curve { +private: + SBasisCurve(); + D2 inner; +public: + explicit SBasisCurve(D2 const &sb) : inner(sb) {} + explicit SBasisCurve(Curve const &other) : inner(other.toSBasis()) {} + Curve *duplicate() const { return new SBasisCurve(*this); } -struct BezierHelpers { -protected: - static Rect bounds(unsigned degree, Point const *points); - static Point point_and_derivatives_at(Coord t, - unsigned degree, Point const *points, - unsigned n_derivs, Point *derivs); - static Point subdivideArr(Coord t, unsigned degree, Point const *V, - Point *Left, Point *Right); + Point initialPoint() const { return inner.at0(); } + Point finalPoint() const { return inner.at1(); } + Point pointAt(Coord t) const { return inner.valueAt(t); } + std::vector pointAndDerivatives(Coord t, unsigned n) const { + return inner.valueAndDerivatives(t, n); + } + double valueAt(Coord t, Dim2 d) const { return inner[d].valueAt(t); } + + void setInitial(Point v) { for(unsigned d = 0; d < 2; d++) { inner[d][0][0] = v[d]; } } + void setFinal(Point v) { for(unsigned d = 0; d < 2; d++) { inner[d][0][1] = v[d]; } } + + Rect boundsFast() const { return bounds_fast(inner); } + Rect boundsExact() const { return bounds_exact(inner); } + Rect boundsLocal(Interval i, unsigned deg) const { return bounds_local(inner, i, deg); } + + std::vector roots(double v, Dim2 d) const { return Geom::roots(inner[d] - v); } + + Curve *portion(double f, double t) const { + return new SBasisCurve(Geom::portion(inner, f, t)); + } + + Curve *transformed(Matrix const &m) const { + return new SBasisCurve(inner * m); + } + + Curve *derivative() const { + return new SBasisCurve(Geom::derivative(inner)); + } + + D2 toSBasis() const { return inner; } }; -template -class Bezier : public Curve, private CurveHelpers, private BezierHelpers { +template +class BezierCurve : public Curve { +private: + D2 > inner; public: template - static void assert_degree(Bezier const *) {} + static void assert_degree(BezierCurve const *) {} + + BezierCurve() {} - Bezier() {} + explicit BezierCurve(D2 > const &x) : inner(x) {} + + BezierCurve(Bezier x, Bezier y) : inner(x, y) {} // default copy // default assign - Bezier(Point c0, Point c1) { + BezierCurve(Point c0, Point c1) { assert_degree<1>(this); - c_[0] = c0; - c_[1] = c1; + for(unsigned d = 0; d < 2; d++) + inner[d] = Bezier(c0[d], c1[d]); } - Bezier(Point c0, Point c1, Point c2) { + BezierCurve(Point c0, Point c1, Point c2) { assert_degree<2>(this); - c_[0] = c0; - c_[1] = c1; - c_[2] = c2; + for(unsigned d = 0; d < 2; d++) + inner[d] = Bezier(c0[d], c1[d], c2[d]); } - Bezier(Point c0, Point c1, Point c2, Point c3) { + BezierCurve(Point c0, Point c1, Point c2, Point c3) { assert_degree<3>(this); - c_[0] = c0; - c_[1] = c1; - c_[2] = c2; - c_[3] = c3; + for(unsigned d = 0; d < 2; d++) + inner[d] = Bezier(c0[d], c1[d], c2[d], c3[d]); } - unsigned degree() const { return bezier_degree; } + unsigned degree() const { return order; } - Curve *duplicate() const { return new Bezier(*this); } + Curve *duplicate() const { return new BezierCurve(*this); } - Point initialPoint() const { return c_[0]; } - Point finalPoint() const { return c_[bezier_degree]; } + Point initialPoint() const { return inner.at0(); } + Point finalPoint() const { return inner.at1(); } - Point &operator[](int index) { return c_[index]; } - Point const &operator[](int index) const { return c_[index]; } + void setInitial(Point v) { setPoint(0, v); } + void setFinal(Point v) { setPoint(1, v); } - Rect bounds_fast() const { return bounds(bezier_degree, c_); } - Rect bounds_exact() const { return bounds(bezier_degree, c_); } + void setPoint(unsigned ix, Point v) { inner[X].setPoint(ix, v[X]); inner[Y].setPoint(ix, v[Y]); } + Point const operator[](unsigned ix) const { return Point(inner[X][ix], inner[Y][ix]); } + + Rect boundsFast() const { return bounds_fast(inner); } + Rect boundsExact() const { return bounds_exact(inner); } + Rect boundsLocal(Interval i, unsigned deg) const { + if(i.min() == 0 && i.max() == 1) return boundsFast(); + if(deg == 0) return bounds_local(inner, i); + // TODO: UUUUUUGGGLLY + if(deg == 1 && order > 1) return Rect(bounds_local(Geom::derivative(inner[X]), i), + bounds_local(Geom::derivative(inner[Y]), i)); + return Rect(Interval(0,0), Interval(0,0)); + } +//TODO: local - boost::optional winding(Point p) const { - return sbasis_winding(sbasis(), p); +//TODO: implement next 3 natively + int winding(Point p) const { + return SBasisCurve(toSBasis()).winding(p); } - Path const &subdivide(Coord t, Path &out) const; + std::vector + roots(double v, Dim2 d) const { + return (inner[d] - v).roots(); + } - Point pointAndDerivativesAt(Coord t, unsigned n_derivs, Point *derivs) - const - { - return point_and_derivatives_at(t, bezier_degree, c_, n_derivs, derivs); + void setPoints(std::vector ps) { + for(unsigned i = 0; i <= order; i++) { + setPoint(i, ps[i]); + } + } + std::vector points() const { return bezier_points(inner); } + + std::pair, BezierCurve > subdivide(Coord t) const { + std::pair, Bezier > sx = inner[X].subdivide(t), sy = inner[Y].subdivide(t); + return std::pair, BezierCurve >( + BezierCurve(sx.first, sy.first), + BezierCurve(sx.second, sy.second)); + } + + Curve *portion(double f, double t) const { + return new BezierCurve(Geom::portion(inner, f, t)); } - D2 sbasis() const { - return handles_to_sbasis(c_); + Curve *reverse() const { + return new BezierCurve(Geom::reverse(inner)); } -protected: - Bezier(Point c[]) { - std::copy(c, c+bezier_degree+1, c_); + Curve *transformed(Matrix const &m) const { + BezierCurve *ret = new BezierCurve(); + std::vector ps = points(); + for(unsigned i = 0; i <= order; i++) ps[i] = ps[i] * m; + ret->setPoints(ps); + return ret; } -private: - Point c_[bezier_degree+1]; + Curve *derivative() const { + if(order > 1) + return new BezierCurve(Geom::derivative(inner[X]), Geom::derivative(inner[Y])); + else if (order == 1) { + double dx = inner[X][1] - inner[X][0], dy = inner[Y][1] - inner[Y][0]; + if(dx == 0) return new BezierCurve<1>(Point(0,0), Point(0,0)); + double slope = dy / dx; + Geom::Point pnt; + if(slope == 0) pnt = Geom::Point(0, 0); else pnt = Geom::Point(slope, 1./slope); + return new BezierCurve<1>(pnt, pnt); + } + } + + Point pointAt(double t) const { return inner.valueAt(t); } + std::vector pointAndDerivatives(Coord t, unsigned n) const { return inner.valueAndDerivatives(t, n); } + + double valueAt(double t, Dim2 d) const { return inner[d].valueAt(t); } + + D2 toSBasis() const {return inner.toSBasis(); } + +protected: + BezierCurve(Point c[]) { + Coord x[order+1], y[order+1]; + for(unsigned i = 0; i <= order; i++) { + x[i] = c[i][X]; y[i] = c[i][Y]; + } + inner = Bezier(x, y); + } }; -// Bezier<0> is meaningless; specialize it out -template<> class Bezier<0> { Bezier(); }; +// BezierCurve<0> is meaningless; specialize it out +template<> class BezierCurve<0> : public BezierCurve<1> { public: BezierCurve(); BezierCurve(Bezier<0> x, Bezier<0> y); }; -typedef Bezier<1> LineSegment; -typedef Bezier<2> QuadraticBezier; -typedef Bezier<3> CubicBezier; +typedef BezierCurve<1> LineSegment; +typedef BezierCurve<2> QuadraticBezier; +typedef BezierCurve<3> CubicBezier; -class SVGEllipticalArc : public Curve, private CurveHelpers { +class SVGEllipticalArc : public Curve { public: SVGEllipticalArc() {} @@ -173,18 +278,55 @@ public: Point initialPoint() const { return initial_; } Point finalPoint() const { return final_; } - Rect bounds_fast() const; - Rect bounds_exact() const; + void setInitial(Point v) { initial_ = v; } + void setFinal(Point v) { final_ = v; } + + //TODO: implement funcs - boost::optional winding(Point p) const { - return sbasis_winding(sbasis(), p); + Rect boundsFast() const; + Rect boundsExact() const; + Rect boundsLocal(Interval i, unsigned deg) const; + + int winding(Point p) const { + return SBasisCurve(toSBasis()).winding(p); } + + std::vector roots(double v, Dim2 d) const; - Path const &subdivide(Coord t, Path &out) const; + inline std::pair + subdivide(Coord t) { + SVGEllipticalArc a(*this), b(*this); + a.final_ = b.initial_ = pointAt(t); + return std::pair(a, b); + } + + Curve *portion(double f, double t) const { + SVGEllipticalArc *ret = new SVGEllipticalArc (*this); + ret->initial_ = pointAt(f); + ret->final_ = pointAt(t); + return ret; + } + + Curve *reverse(double f, double t) const { + SVGEllipticalArc *ret = new SVGEllipticalArc (*this); + ret->initial_ = final_; + ret->final_ = initial_; + return ret; + } - Point pointAndDerivativesAt(Coord t, unsigned n_derivs, Point *derivs) const; + //TODO: this next def isn't right + Curve *transformed(Matrix const & m) const { + SVGEllipticalArc *ret = new SVGEllipticalArc (*this); + ret->initial_ = initial_ * m; + ret->final_ = final_ * m; + return ret; + } + + Curve *derivative() const { throw NotImplemented(); } + + std::vector pointAndDerivatives(Coord t, unsigned n) const; - D2 sbasis() const; + D2 toSBasis() const; private: Point initial_; @@ -196,39 +338,6 @@ private: Point final_; }; -class SBasisCurve : public Curve, private CurveHelpers { -private: - SBasisCurve(); -public: - explicit SBasisCurve(D2 const &coeffs) - : coeffs_(coeffs) {} - - Point initialPoint() const { - return Point(coeffs_[X][0][0], coeffs_[Y][0][0]); - } - Point finalPoint() const { - return Point(coeffs_[X][0][1], coeffs_[Y][0][1]); - } - - Curve *duplicate() const { return new SBasisCurve(*this); } - - Rect bounds_fast() const; - Rect bounds_exact() const; - - boost::optional winding(Point p) const { - return sbasis_winding(coeffs_, p); - } - - Path const &subdivide(Coord t, Path &out) const; - - Point pointAndDerivativesAt(Coord t, unsigned n_derivs, Point *derivs) const; - - D2 sbasis() const { return coeffs_; } - -private: - D2 coeffs_; -}; - template class BaseIterator : public std::iterator @@ -253,6 +362,7 @@ public: ++impl_; return *this; } + BaseIterator operator++(int) { BaseIterator old=*this; ++(*this); @@ -282,7 +392,7 @@ public: } Curve *operator*() const { return (*impl_)->duplicate(); } - + DuplicatingIterator &operator++() { ++impl_; return *this; @@ -378,20 +488,77 @@ public: bool closed() const { return closed_; } void close(bool closed=true) { closed_ = closed; } - int winding(Point p) const; - - Rect bounds_fast() const; - Rect bounds_exact() const; + Rect boundsFast() const; + Rect boundsExact() const; Piecewise > toPwSb() const { Piecewise > ret; ret.push_cut(0); - for(unsigned i = 0; i < size() + (closed_ ? 1 : 0); i++) { - ret.push(curves_[i]->sbasis(), i+1); + unsigned i = 1; + for(const_iterator it = begin(); it != end_default(); ++it, i++) { + ret.push(it->toSBasis(), i); + } + return ret; + } + + Path operator*(Matrix const &m) const { + Path ret; + for(const_iterator it = begin(); it != end(); ++it) { + Curve *temp = it->transformed(m); + //Possible point of discontinuity? + ret.append(*temp); + delete temp; } return ret; } + Point pointAt(double t) const { + if(empty()) return Point(0,0); + double i, f = modf(t, &i); + if(i == size() && f == 0) { i--; } + assert(i >= 0 && i <= size()); + return (*this)[unsigned(i)].pointAt(f); + } + + double valueAt(double t, Dim2 d) const { + if(empty()) return 0; + double i, f = modf(t, &i); + if(i == size() && f == 0) { i--; } + assert(i >= 0 && i <= size()); + return (*this)[unsigned(i)].valueAt(f, d); + } + + std::vector roots(double v, Dim2 d) const { + std::vector res; + for(const_iterator it = begin(); it != end_closed(); ++it) { + std::vector temp = it->roots(v, d); + res.insert(res.end(), temp.begin(), temp.end()); + } + return res; + } + + void appendPortionTo(Path &p, double f, double t) const; + + Path portion(double f, double t) const { + Path ret; + ret.close(false); + appendPortionTo(ret, f, t); + return ret; + } + Path portion(Interval i) const { return portion(i.min(), i.max()); } + + Path reverse() const { + Path ret; + ret.close(closed_); + for(int i = size() - (closed_ ? 0 : 1); i >= 0; i--) { + //TODO: do we really delete? + Curve *temp = (*this)[i].reverse(); + ret.append(*temp); + delete temp; + } + return ret; + } + void insert(iterator pos, Curve const &curve) { Sequence source(1, curve.duplicate()); try { @@ -481,14 +648,14 @@ public: void start(Point p) { clear(); - (*final_)[0] = (*final_)[1] = p; + final_->setPoint(0, p); + final_->setPoint(1, p); } Point initialPoint() const { return (*final_)[1]; } Point finalPoint() const { return (*final_)[0]; } void append(Curve const &curve); - void append(D2 const &curve); template @@ -573,37 +740,52 @@ inline static Piecewise > paths_to_pw(std::vector paths) { return ret; } -template -inline Path const &Bezier::subdivide(Coord t, Path &out) const { - Bezier a, b; - subdivideArr(t, bezier_degree, c_, a.c_, b.c_); - out.clear(); - out.close(false); - out.append(a); - out.append(b); - return out; -} +/* +class PathPortion : public Curve { + Path *source; + double f, t; + boost::optional result; + + public: + double from() const { return f; } + double to() const { return t; } + + explicit PathPortion(Path *s, double fp, double tp) : source(s), f(fp), t(tp) {} + Curve *duplicate() const { return new PathPortion(*this); } + + Point initialPoint() const { return source->pointAt(f); } + Point finalPoint() const { return source->pointAt(t); } + + Path actualPath() { + if(!result) *result = source->portion(f, t); + return *result; + } + + Rect boundsFast() const { return actualPath().boundsFast; } + Rect boundsExact() const { return actualPath().boundsFast; } + Rect boundsLocal(Interval i) const { throw NotImplemented(); } -inline Path const &SVGEllipticalArc::subdivide(Coord t, Path &out) const { - SVGEllipticalArc a; - SVGEllipticalArc b; - a.rx_ = b.rx_ = rx_; - a.ry_ = b.ry_ = ry_; - a.x_axis_rotation_ = b.x_axis_rotation_ = x_axis_rotation_; - a.initial_ = initial_; - a.final_ = b.initial_ = pointAt(t); - b.final_ = final_; - out.clear(); - out.close(false); - out.append(a); - out.append(b); - return out; -} + std::vector roots(double v, Dim2 d) const = 0; -class Set { -public: -private: + virtual int winding(Point p) const { return root_winding(*this, p); } + + virtual Curve *portion(double f, double t) const = 0; + virtual Curve *reverse() const { return portion(1, 0); } + + virtual Crossings crossingsWith(Curve const & other) const; + + virtual void setInitial(Point v) = 0; + virtual void setFinal(Point v) = 0; + + virtual Curve *transformed(Matrix const &m) const = 0; + + virtual Point pointAt(Coord t) const { return pointAndDerivatives(t, 1).front(); } + virtual Coord valueAt(Coord t, Dim2 d) const { return pointAt(t)[d]; } + virtual std::vector pointAndDerivatives(Coord t, unsigned n) const = 0; + virtual D2 toSBasis() const = 0; + }; +*/ } diff --git a/src/2geom/poly-dk-solve.cpp b/src/2geom/poly-dk-solve.cpp index dde230dc4..136687ae8 100644 --- a/src/2geom/poly-dk-solve.cpp +++ b/src/2geom/poly-dk-solve.cpp @@ -1,64 +1,64 @@ -#include "poly-dk-solve.h" -#include - -/*** implementation of the Durand-Kerner method. seems buggy*/ - -std::complex evalu(Poly const & p, std::complex x) { - std::complex result = 0; - std::complex xx = 1; - - for(unsigned i = 0; i < p.size(); i++) { - result += p[i]*xx; - xx *= x; - } - return result; -} - -std::vector > DK(Poly const & ply, const double tol) { - std::vector > roots; - const int N = ply.degree(); - - std::complex b(0.4, 0.9); - std::complex p = 1; - for(int i = 0; i < N; i++) { - roots.push_back(p); - p *= b; - } - assert(roots.size() == ply.degree()); - - double error = 0; - int i; - for( i = 0; i < 30; i++) { - error = 0; - for(int r_i = 0; r_i < N; r_i++) { - std::complex denom = 1; - std::complex R = roots[r_i]; - for(int d_i = 0; d_i < N; d_i++) { - if(r_i != d_i) - denom *= R-roots[d_i]; - } - assert(norm(denom) != 0); - std::complex dr = evalu(ply, R)/denom; - error += norm(dr); - roots[r_i] = R - dr; - } - /*std::copy(roots.begin(), roots.end(), std::ostream_iterator >(std::cout, ",\t")); - std::cout << std::endl;*/ - if(error < tol) - break; - } - //std::cout << error << ", " << i<< std::endl; - return roots; -} - - -/* - Local Variables: - mode:c++ - c-file-style:"stroustrup" - c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) - indent-tabs-mode:nil - fill-column:99 - End: -*/ -// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 : +#include "poly-dk-solve.h" +#include + +/*** implementation of the Durand-Kerner method. seems buggy*/ + +std::complex evalu(Poly const & p, std::complex x) { + std::complex result = 0; + std::complex xx = 1; + + for(unsigned i = 0; i < p.size(); i++) { + result += p[i]*xx; + xx *= x; + } + return result; +} + +std::vector > DK(Poly const & ply, const double tol) { + std::vector > roots; + const int N = ply.degree(); + + std::complex b(0.4, 0.9); + std::complex p = 1; + for(int i = 0; i < N; i++) { + roots.push_back(p); + p *= b; + } + assert(roots.size() == ply.degree()); + + double error = 0; + int i; + for( i = 0; i < 30; i++) { + error = 0; + for(int r_i = 0; r_i < N; r_i++) { + std::complex denom = 1; + std::complex R = roots[r_i]; + for(int d_i = 0; d_i < N; d_i++) { + if(r_i != d_i) + denom *= R-roots[d_i]; + } + assert(norm(denom) != 0); + std::complex dr = evalu(ply, R)/denom; + error += norm(dr); + roots[r_i] = R - dr; + } + /*std::copy(roots.begin(), roots.end(), std::ostream_iterator >(std::cout, ",\t")); + std::cout << std::endl;*/ + if(error < tol) + break; + } + //std::cout << error << ", " << i<< std::endl; + return roots; +} + + +/* + Local Variables: + mode:c++ + c-file-style:"stroustrup" + c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) + indent-tabs-mode:nil + fill-column:99 + End: +*/ +// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 : diff --git a/src/2geom/poly-laguerre-solve.cpp b/src/2geom/poly-laguerre-solve.cpp index 83f049286..fb3c2075b 100644 --- a/src/2geom/poly-laguerre-solve.cpp +++ b/src/2geom/poly-laguerre-solve.cpp @@ -1,147 +1,147 @@ -#include "poly-laguerre-solve.h" -#include - -typedef std::complex cdouble; - -cdouble laguerre_internal_complex(Poly const & p, - double x0, - double tol, - bool & quad_root) { - cdouble a = 2*tol; - cdouble xk = x0; - double n = p.degree(); - quad_root = false; - const unsigned shuffle_rate = 10; - static double shuffle[] = {0, 0.5, 0.25, 0.75, 0.125, 0.375, 0.625, 0.875, 1.0}; - unsigned shuffle_counter = 0; - while(std::norm(a) > (tol*tol)) { - //std::cout << "xk = " << xk << std::endl; - cdouble b = p.back(); - cdouble d = 0, f = 0; - double err = abs(b); - double abx = abs(xk); - for(int j = p.size()-2; j >= 0; j--) { - f = xk*f + d; - d = xk*d + b; - b = xk*b + p[j]; - err = abs(b) + abx*err; - } - - err *= 1e-7; // magic epsilon for convergence, should be computed from tol - - cdouble px = b; - if(abs(b) < err) - return xk; - //if(std::norm(px) < tol*tol) - // return xk; - cdouble G = d / px; - cdouble H = G*G - f / px; - - //std::cout << "G = " << G << "H = " << H; - cdouble radicand = (n - 1)*(n*H-G*G); - //assert(radicand.real() > 0); - if(radicand.real() < 0) - quad_root = true; - //std::cout << "radicand = " << radicand << std::endl; - if(G.real() < 0) // here we try to maximise the denominator avoiding cancellation - a = - sqrt(radicand); - else - a = sqrt(radicand); - //std::cout << "a = " << a << std::endl; - a = n / (a + G); - //std::cout << "a = " << a << std::endl; - if(shuffle_counter % shuffle_rate == 0) - ;//a *= shuffle[shuffle_counter / shuffle_rate]; - xk -= a; - shuffle_counter++; - if(shuffle_counter >= 90) - break; - } - //std::cout << "xk = " << xk << std::endl; - return xk; -} - -double laguerre_internal(Poly const & p, - Poly const & pp, - Poly const & ppp, - double x0, - double tol, - bool & quad_root) { - double a = 2*tol; - double xk = x0; - double n = p.degree(); - quad_root = false; - while(a*a > (tol*tol)) { - //std::cout << "xk = " << xk << std::endl; - double px = p(xk); - if(px*px < tol*tol) - return xk; - double G = pp(xk) / px; - double H = G*G - ppp(xk) / px; - - //std::cout << "G = " << G << "H = " << H; - double radicand = (n - 1)*(n*H-G*G); - assert(radicand > 0); - //std::cout << "radicand = " << radicand << std::endl; - if(G < 0) // here we try to maximise the denominator avoiding cancellation - a = - sqrt(radicand); - else - a = sqrt(radicand); - //std::cout << "a = " << a << std::endl; - a = n / (a + G); - //std::cout << "a = " << a << std::endl; - xk -= a; - } - //std::cout << "xk = " << xk << std::endl; - return xk; -} - - -std::vector -laguerre(Poly p, const double tol) { - std::vector solutions; - //std::cout << "p = " << p << " = "; - while(p.size() > 1) - { - double x0 = 0; - bool quad_root = false; - cdouble sol = laguerre_internal_complex(p, x0, tol, quad_root); - //if(abs(sol) > 1) break; - Poly dvs; - if(quad_root) { - dvs.push_back((sol*conj(sol)).real()); - dvs.push_back(-(sol + conj(sol)).real()); - dvs.push_back(1.0); - //std::cout << "(" << dvs << ")"; - //solutions.push_back(sol); - //solutions.push_back(conj(sol)); - } else { - //std::cout << sol << std::endl; - dvs.push_back(-sol.real()); - dvs.push_back(1.0); - solutions.push_back(sol); - //std::cout << "(" << dvs << ")"; - } - Poly r; - p = divide(p, dvs, r); - //std::cout << r << std::endl; - } - return solutions; -} - -std::vector -laguerre_real_interval(Poly const & ply, - const double lo, const double hi, - const double tol) { -} - -/* - Local Variables: - mode:c++ - c-file-style:"stroustrup" - c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) - indent-tabs-mode:nil - fill-column:99 - End: -*/ -// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 : +#include "poly-laguerre-solve.h" +#include + +typedef std::complex cdouble; + +cdouble laguerre_internal_complex(Poly const & p, + double x0, + double tol, + bool & quad_root) { + cdouble a = 2*tol; + cdouble xk = x0; + double n = p.degree(); + quad_root = false; + const unsigned shuffle_rate = 10; + static double shuffle[] = {0, 0.5, 0.25, 0.75, 0.125, 0.375, 0.625, 0.875, 1.0}; + unsigned shuffle_counter = 0; + while(std::norm(a) > (tol*tol)) { + //std::cout << "xk = " << xk << std::endl; + cdouble b = p.back(); + cdouble d = 0, f = 0; + double err = abs(b); + double abx = abs(xk); + for(int j = p.size()-2; j >= 0; j--) { + f = xk*f + d; + d = xk*d + b; + b = xk*b + p[j]; + err = abs(b) + abx*err; + } + + err *= 1e-7; // magic epsilon for convergence, should be computed from tol + + cdouble px = b; + if(abs(b) < err) + return xk; + //if(std::norm(px) < tol*tol) + // return xk; + cdouble G = d / px; + cdouble H = G*G - f / px; + + //std::cout << "G = " << G << "H = " << H; + cdouble radicand = (n - 1)*(n*H-G*G); + //assert(radicand.real() > 0); + if(radicand.real() < 0) + quad_root = true; + //std::cout << "radicand = " << radicand << std::endl; + if(G.real() < 0) // here we try to maximise the denominator avoiding cancellation + a = - sqrt(radicand); + else + a = sqrt(radicand); + //std::cout << "a = " << a << std::endl; + a = n / (a + G); + //std::cout << "a = " << a << std::endl; + if(shuffle_counter % shuffle_rate == 0) + ;//a *= shuffle[shuffle_counter / shuffle_rate]; + xk -= a; + shuffle_counter++; + if(shuffle_counter >= 90) + break; + } + //std::cout << "xk = " << xk << std::endl; + return xk; +} + +double laguerre_internal(Poly const & p, + Poly const & pp, + Poly const & ppp, + double x0, + double tol, + bool & quad_root) { + double a = 2*tol; + double xk = x0; + double n = p.degree(); + quad_root = false; + while(a*a > (tol*tol)) { + //std::cout << "xk = " << xk << std::endl; + double px = p(xk); + if(px*px < tol*tol) + return xk; + double G = pp(xk) / px; + double H = G*G - ppp(xk) / px; + + //std::cout << "G = " << G << "H = " << H; + double radicand = (n - 1)*(n*H-G*G); + assert(radicand > 0); + //std::cout << "radicand = " << radicand << std::endl; + if(G < 0) // here we try to maximise the denominator avoiding cancellation + a = - sqrt(radicand); + else + a = sqrt(radicand); + //std::cout << "a = " << a << std::endl; + a = n / (a + G); + //std::cout << "a = " << a << std::endl; + xk -= a; + } + //std::cout << "xk = " << xk << std::endl; + return xk; +} + + +std::vector +laguerre(Poly p, const double tol) { + std::vector solutions; + //std::cout << "p = " << p << " = "; + while(p.size() > 1) + { + double x0 = 0; + bool quad_root = false; + cdouble sol = laguerre_internal_complex(p, x0, tol, quad_root); + //if(abs(sol) > 1) break; + Poly dvs; + if(quad_root) { + dvs.push_back((sol*conj(sol)).real()); + dvs.push_back(-(sol + conj(sol)).real()); + dvs.push_back(1.0); + //std::cout << "(" << dvs << ")"; + //solutions.push_back(sol); + //solutions.push_back(conj(sol)); + } else { + //std::cout << sol << std::endl; + dvs.push_back(-sol.real()); + dvs.push_back(1.0); + solutions.push_back(sol); + //std::cout << "(" << dvs << ")"; + } + Poly r; + p = divide(p, dvs, r); + //std::cout << r << std::endl; + } + return solutions; +} + +std::vector +laguerre_real_interval(Poly const & ply, + const double lo, const double hi, + const double tol) { +} + +/* + Local Variables: + mode:c++ + c-file-style:"stroustrup" + c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) + indent-tabs-mode:nil + fill-column:99 + End: +*/ +// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 : diff --git a/src/2geom/poly.cpp b/src/2geom/poly.cpp index 5c3a30002..17f286494 100644 --- a/src/2geom/poly.cpp +++ b/src/2geom/poly.cpp @@ -112,8 +112,8 @@ Poly derivative(Poly const & p) { Poly compose(Poly const & a, Poly const & b) { Poly result; - for(unsigned i = a.size()-1; i >=0; i--) { - result = Poly(a[i]) + result * b; + for(unsigned i = a.size(); i > 0; i--) { + result = Poly(a[i-1]) + result * b; } return result; diff --git a/src/2geom/rect.h b/src/2geom/rect.h index d96802014..f0811baf6 100644 --- a/src/2geom/rect.h +++ b/src/2geom/rect.h @@ -1,91 +1,153 @@ -//D2 specialization to Rect: - - /* Authors of original rect class: - * Lauris Kaplinski - * Nathan Hurst - * bulia byak - * MenTaLguY - */ +/* + * rect.h - D2 specialization to Rect + * + * Copyright 2007 Michael Sloan + * + * 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, output 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. + * + */ + +/* Authors of original rect class: + * Lauris Kaplinski + * Nathan Hurst + * bulia byak + * MenTaLguY + */ #ifdef _2GEOM_D2 /*This is intentional: we don't actually want anyone to - include this, other than D2.h */ + include this, other than D2.h. If somone else tries, D2 + won't be defined. If it is, this will already be included. */ #ifndef _2GEOM_RECT #define _2GEOM_RECT - -#include "matrix.h" + +#include "matrix.h" #include - -namespace Geom { - -typedef D2 Rect; + +namespace Geom { + +typedef D2 Rect; Rect unify(const Rect &, const Rect &); - -template<> -class D2 { - private: - Interval f[2]; - D2();// { f[X] = f[Y] = Interval(0, 0); } - - public: - D2(Interval const &a, Interval const &b) { - f[X] = a; - f[Y] = b; - } - - D2(Point const & a, Point const & b) { - f[X] = Interval(a[X], b[X]); - f[Y] = Interval(a[Y], b[Y]); - } - - inline Interval& operator[](unsigned i) { return f[i]; } - inline Interval const & operator[](unsigned i) const { return f[i]; } - - inline Point min() const { return Point(f[X].min(), f[Y].min()); } - inline Point max() const { return Point(f[X].max(), f[Y].max()); } - - /** returns the four corners of the rectangle in order - * (clockwise if +Y is up, anticlockwise if +Y is down) */ - Point corner(unsigned i) const { - switch(i % 4) { - case 0: return Point(f[X].min(), f[Y].min()); - case 1: return Point(f[X].max(), f[Y].min()); - case 2: return Point(f[X].max(), f[Y].max()); case 3: return Point(f[X].min(), f[Y].max()); - } - } +template<> +class D2 { + private: + Interval f[2]; + D2();// { f[X] = f[Y] = Interval(0, 0); } + + public: + D2(Interval const &a, Interval const &b) { + f[X] = a; + f[Y] = b; + } + + D2(Point const & a, Point const & b) { + f[X] = Interval(a[X], b[X]); + f[Y] = Interval(a[Y], b[Y]); + } + + inline Interval& operator[](unsigned i) { return f[i]; } + inline Interval const & operator[](unsigned i) const { return f[i]; } + + inline Point min() const { return Point(f[X].min(), f[Y].min()); } + inline Point max() const { return Point(f[X].max(), f[Y].max()); } + + /** returns the four corners of the rectangle in positive order + * (clockwise if +Y is up, anticlockwise if +Y is down) */ + Point corner(unsigned i) const { + switch(i % 4) { + case 0: return Point(f[X].min(), f[Y].min()); + case 1: return Point(f[X].max(), f[Y].min()); + case 2: return Point(f[X].max(), f[Y].max()); + case 3: return Point(f[X].min(), f[Y].max()); + } + } + + //We should probably remove these - they're coord sys gnostic + inline double top() const { return f[Y].min(); } + inline double bottom() const { return f[Y].max(); } + inline double left() const { return f[X].min(); } + inline double right() const { return f[X].max(); } + inline double width() const { return f[X].extent(); } inline double height() const { return f[Y].extent(); } - - /** returns a vector from min to max. */ - inline Point dimensions() const { return Point(f[X].extent(), f[Y].extent()); } - inline Point midpoint() const { return Point(f[X].middle(), f[Y].middle()); } - - inline double area() const { return f[X].extent() * f[Y].extent(); } - inline double maxExtent() const { return std::max(f[X].extent(), f[Y].extent()); } - - inline bool isEmpty() const { 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]); } - inline bool contains(Rect const &r) const { return f[X].contains(r[X]) && f[Y].contains(r[Y]); } - inline bool contains(Point const &p) const { return f[X].contains(p[X]) && f[Y].contains(p[Y]); } - - inline void expandTo(Point p) { f[X].extendTo(p[X]); f[Y].extendTo(p[Y]); } - inline void unionWith(Rect const &b) { f[X].unionWith(b[X]); f[Y].unionWith(b[Y]); } - - inline void expandBy(double amnt) { f[X].expandBy(amnt); f[Y].expandBy(amnt); } - inline void expandBy(Point const p) { f[X].expandBy(p[X]); f[Y].expandBy(p[Y]); } - + + /** returns a vector from min to max. */ + inline Point dimensions() const { return Point(f[X].extent(), f[Y].extent()); } + inline Point midpoint() const { return Point(f[X].middle(), f[Y].middle()); } + + inline double area() const { return f[X].extent() * f[Y].extent(); } + inline double maxExtent() const { return std::max(f[X].extent(), f[Y].extent()); } + + inline bool isEmpty() const { + 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]); + } + inline bool contains(Rect const &r) const { + return f[X].contains(r[X]) && f[Y].contains(r[Y]); + } + inline bool contains(Point const &p) const { + return f[X].contains(p[X]) && f[Y].contains(p[Y]); + } + + inline void expandTo(Point p) { + f[X].extendTo(p[X]); f[Y].extendTo(p[Y]); + } + inline void unionWith(Rect const &b) { + f[X].unionWith(b[X]); f[Y].unionWith(b[Y]); + } + + inline void expandBy(double amnt) { + f[X].expandBy(amnt); f[Y].expandBy(amnt); + } + inline void expandBy(Point const p) { + f[X].expandBy(p[X]); f[Y].expandBy(p[Y]); + } + /** Transforms the rect by m. Note that it gives correct results only for scales and translates, - in the case of rotations, the area of the rect will grow as it cannot rotate. */ - inline Rect operator*(Matrix const m) const { return unify(Rect(corner(0) * m, corner(2) * m), - Rect(corner(1) * m, corner(3) * m)); } + in the case of rotations, the area of the rect will grow as it cannot rotate. */ + inline Rect operator*(Matrix const m) const { + return unify(Rect(corner(0) * m, corner(2) * m), + Rect(corner(1) * m, corner(3) * m)); + } }; -inline Rect unify(const Rect & a, const Rect & b) { +inline Rect unify(Rect const & a, Rect const & b) { return Rect(unify(a[X], b[X]), unify(a[Y], b[Y])); } -inline boost::optional intersect(const Rect & a, const Rect & b) { +inline Rect union_list(std::vector const &r) { + if(r.empty()) return Rect(Interval(0,0), Interval(0,0)); + Rect ret = r[0]; + for(unsigned i = 1; i < r.size(); i++) + ret.unionWith(r[i]); + return ret; +} + +inline boost::optional intersect(Rect const & a, Rect const & b) { boost::optional x = intersect(a[X], b[X]); boost::optional y = intersect(a[Y], b[Y]); return x && y ? boost::optional(Rect(*x, *y)) : boost::optional(); @@ -95,3 +157,13 @@ inline boost::optional intersect(const Rect & a, const Rect & b) { #endif //_2GEOM_RECT #endif //_2GEOM_D2 +/* + Local Variables: + mode:c++ + c-file-style:"stroustrup" + c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) + indent-tabs-mode:nil + fill-column:99 + End: +*/ +// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 : diff --git a/src/2geom/region.cpp b/src/2geom/region.cpp new file mode 100644 index 000000000..cfab3a35a --- /dev/null +++ b/src/2geom/region.cpp @@ -0,0 +1,42 @@ +#include "region.h" +#include "utils.h" + +#include "shape.h" + +namespace Geom { + +Regions sanitize_path(Path const &p) { + Regions results; + Crossings crs = self_crossings(p); + for(unsigned i = 0; i < crs.size(); i++) { + + } +} + +Region Region::operator*(Matrix const &m) const { + Region r((m.flips() ? boundary.reverse() : boundary) * m, fill); + if(box && m.onlyScaleAndTranslation()) r.box = (*box) * m; + return r; +} + +bool Region::invariants() const { + return self_crossings(boundary).empty(); +} + +unsigned outer_index(Regions const &ps) { + if(ps.size() <= 1 || ps[0].contains(ps[1])) { + return 0; + } else { + /* Since we've already shown that chunks[0] is not outside + it can be used as an exemplar inner. */ + Point exemplar = Path(ps[0]).initialPoint(); + for(unsigned i = 1; i < ps.size(); i++) { + if(ps[i].contains(exemplar)) { + return i; + } + } + } + return ps.size(); +} + +} diff --git a/src/2geom/region.h b/src/2geom/region.h new file mode 100644 index 000000000..4b434f1e5 --- /dev/null +++ b/src/2geom/region.h @@ -0,0 +1,80 @@ +#ifndef __2GEOM_REGION_H +#define __2GEOM_REGION_H + +#include "path.h" +#include "path-intersection.h" + +namespace Geom { + +class Shape; + +class Region { + friend Crossings crossings(Region const &a, Region const &b); + friend class Shape; + friend Shape shape_boolean(bool rev, Shape const & a, Shape const & b, CrossingSet const & crs); + + Path boundary; + mutable boost::optional box; + bool fill; + public: + Region() : fill(true) {} + explicit Region(Path const &p) : boundary(p) { fill = path_direction(p); } + Region(Path const &p, bool dir) : boundary(p), fill(dir) {} + Region(Path const &p, boost::optional const &b) : boundary(p), box(b) { fill = path_direction(p); } + Region(Path const &p, boost::optional const &b, bool dir) : boundary(p), box(b), fill(dir) {} + + unsigned size() const { return boundary.size(); } + + bool isFill() const { return fill; } + Region asFill() const { if(fill) return Region(*this); else return inverse(); } + Region asHole() const { if(fill) return inverse(); else return Region(*this); } + + operator Path() const { return boundary; } + Rect boundsFast() const { + if(!box) box = boost::optional(boundary.boundsFast()); + return *box; + } + bool contains(Point const &p) const { + if(box && !box->contains(p)) return false; + return Geom::contains(boundary, p); + } + bool contains(Region const &other) const { return contains(other.boundary.initialPoint()); } + + Region inverse() const { return Region(boundary.reverse(), box, !fill); } + + Region operator*(Matrix const &m) const; + + bool invariants() const; +}; + +typedef std::vector Regions; + +unsigned outer_index(Regions const &ps); + +//assumes they're already sanitized somewhat +inline Regions regions_from_paths(std::vector const &ps) { + Regions res; + for(unsigned i = 0; i < ps.size(); i++) + res.push_back(Region(ps[i])); + return res; +} + +inline std::vector paths_from_regions(Regions const &rs) { + std::vector res; + for(unsigned i = 0; i < rs.size(); i++) + res.push_back(rs[i]); + return res; +} + +Regions sanitize_path(Path const &p); + +Regions region_boolean(bool rev, Region const & a, Region const & b, Crossings const &cr); +Regions region_boolean(bool rev, Region const & a, Region const & b, Crossings const & cr_a, Crossings const & cr_b); + +inline Regions region_boolean(bool rev, Region const & a, Region const & b) { + return region_boolean(rev, a, b, crossings(a, b)); +} + +} + +#endif diff --git a/src/2geom/sbasis-geometric.cpp b/src/2geom/sbasis-geometric.cpp index 170323c6c..fde907857 100644 --- a/src/2geom/sbasis-geometric.cpp +++ b/src/2geom/sbasis-geometric.cpp @@ -350,7 +350,7 @@ unsigned Geom::centroid(Piecewise > const &p, Point& centroid, double } // join ends centroid_tmp *= 2; - Point final = p[p.size()].at1(), initial = p[0].at0(); + Point final = p[p.size()-1].at1(), initial = p[0].at0(); const double ai = cross(final, initial); atmp += ai; centroid_tmp += (final + initial)*ai; // first moment. diff --git a/src/2geom/sbasis-math.cpp b/src/2geom/sbasis-math.cpp index d88c26832..db4cf5e08 100644 --- a/src/2geom/sbasis-math.cpp +++ b/src/2geom/sbasis-math.cpp @@ -231,11 +231,11 @@ Piecewise reciprocalOnDomain(Interval range, double tol){ if (a<=tol){ reciprocal_fn.push_cut(0); - int i0=(int) floor(log(tol)/log(R)); + int i0=(int) floor(std::log(tol)/std::log(R)); a=pow(R,i0); reciprocal_fn.push(Linear(1/a),a); }else{ - int i0=(int) floor(log(a)/log(R)); + int i0=(int) floor(std::log(a)/std::log(R)); a=pow(R,i0); reciprocal_fn.cuts.push_back(a); } diff --git a/src/2geom/sbasis-math.h b/src/2geom/sbasis-math.h index 529641068..8f4e1612d 100644 --- a/src/2geom/sbasis-math.h +++ b/src/2geom/sbasis-math.h @@ -70,6 +70,9 @@ Piecewise cos( SBasis const &f, double tol=1e-3, int order=3); Piecewise cos(Piecewise const &f, double tol=1e-3, int order=3); Piecewise sin( SBasis const &f, double tol=1e-3, int order=3); Piecewise sin(Piecewise const &f, double tol=1e-3, int order=3); +//-Log--------------------------------------------------------------- +Piecewise log( SBasis const &f, double tol=1e-3, int order=3); +Piecewise log(Piecewiseconst &f, double tol=1e-3, int order=3); //--1/x------------------------------------------------------------ //TODO: change this... diff --git a/src/2geom/sbasis-poly.cpp b/src/2geom/sbasis-poly.cpp index e0fa828f9..ff797920f 100644 --- a/src/2geom/sbasis-poly.cpp +++ b/src/2geom/sbasis-poly.cpp @@ -15,6 +15,8 @@ SBasis poly_to_sbasis(Poly const & p) { } Poly sbasis_to_poly(SBasis const & sb) { + if(sb.isZero()) + return Poly(); Poly S; // (1-x)x = -1*x^2 + 1*x + 0 Poly A, B; B.push_back(0); diff --git a/src/2geom/sbasis-roots.cpp b/src/2geom/sbasis-roots.cpp index ad2dfbda4..c4ef3c16d 100644 --- a/src/2geom/sbasis-roots.cpp +++ b/src/2geom/sbasis-roots.cpp @@ -68,7 +68,7 @@ Interval bounds_fast(const SBasis &sb, int order) { double a=sb[j][0]; double b=sb[j][1]; - double t, v; + double v, t = 0; v = res[0]; if (v<0) t = ((b-a)/v+1)*0.5; if (v>=0 || t<0 || t>1) { @@ -95,7 +95,7 @@ Interval bounds_local(const SBasis &sb, const Interval &i, int order) { double a=sb[j][0]; double b=sb[j][1]; - double t; + double t = 0; if (lo<0) t = ((b-a)/lo+1)*0.5; if (lo>=0 || tt1) { lo = std::min(a*(1-t0)+b*t0+lo*t0*(1-t0),a*(1-t1)+b*t1+lo*t1*(1-t1)); diff --git a/src/2geom/sbasis-to-bezier.cpp b/src/2geom/sbasis-to-bezier.cpp index 613c39288..206f18931 100644 --- a/src/2geom/sbasis-to-bezier.cpp +++ b/src/2geom/sbasis-to-bezier.cpp @@ -184,34 +184,6 @@ path_from_sbasis(D2 const &B, double tol) { //TODO: some of this logic should be lifted into svg-path std::vector path_from_piecewise(Geom::Piecewise > const &B, double tol) { -/* - Geom::PathBuilder pb; - if(B.size() == 0) return pb.peek(); - Geom::Point start = B[0].at0(); - pb.moveTo(start); - for(unsigned i = 0; ; i++) { - if(i+1 == B.size() || !near(B[i+1].at0(), B[i].at1(), tol)) { - //start of a new path - if(near(start, B[i].at1())) { - //it's closed - pb.closePath(); - if(sbasis_size(B[i]) <= 1) { - //last line seg already there - goto no_add; - } - } - build_from_sbasis(pb, B[i], tol); - no_add: - if(i+1 >= B.size()) break; - start = B[i+1].at0(); - pb.moveTo(start); - i++; - } - build_from_sbasis(pb, B[i], tol); - } - pb.finish(); - return pb.peek(); -*/ Geom::PathBuilder pb; if(B.size() == 0) return pb.peek(); Geom::Point start = B[0].at0(); diff --git a/src/2geom/sbasis.h b/src/2geom/sbasis.h index a27d4bfce..b0df4435c 100644 --- a/src/2geom/sbasis.h +++ b/src/2geom/sbasis.h @@ -39,6 +39,7 @@ #include "linear.h" #include "interval.h" +#include "utils.h" namespace Geom { @@ -88,6 +89,12 @@ public: double operator()(double t) const { return valueAt(t); } + + std::vector valueAndDerivatives(double t, unsigned n) const { + //TODO + throw NotImplemented(); + } + SBasis toSBasis() const { return SBasis(*this); } double tailError(unsigned tail) const; diff --git a/src/2geom/shape.cpp b/src/2geom/shape.cpp new file mode 100644 index 000000000..670792521 --- /dev/null +++ b/src/2geom/shape.cpp @@ -0,0 +1,442 @@ +#include "shape.h" +#include "utils.h" +#include "sweep.h" + +#include +#include + +namespace Geom { + +// Utility funcs + +// Yes, xor is !=, but I'm pretty sure this is safer in the event of strange bools +bool logical_xor (bool a, bool b) { return (a || b) && !(a && b); } + +// A little sugar for appending a list to another +template +void append(T &a, T const &b) { + a.insert(a.end(), b.begin(), b.end()); +} + +/* Used within shape_boolean and related functions, as the name describes, finds the + * first false within the list of lists of booleans. + */ +void first_false(std::vector > visited, unsigned &i, unsigned &j) { + for(i = 0, j = 0; i < visited.size(); i++) { + std::vector::iterator unvisited = std::find(visited[i].begin(), visited[i].end(), false); + if(unvisited != visited[i].end()) { + j = unvisited - visited[i].begin(); + break; + } + } +} + +// Finds a crossing in a list of them, given the sorting index. +unsigned find_crossing(Crossings const &cr, Crossing x, unsigned i) { + return std::lower_bound(cr.begin(), cr.end(), x, CrossingOrder(i)) - cr.begin(); +} + +/* This function handles boolean ops on shapes. The first parameter is a bool + * which determines its behavior in each combination of cases. For proper + * fill information and noncrossing behavior, the fill data of the regions + * must be correct. The boolean parameter determines whether the operation + * is a union or a subtraction. Reversed paths represent inverse regions, + * where everything is included in the fill except for the insides. + * + * Here is a chart of the behavior under various circumstances: + * + * rev = false (union) + * A + * F H + * F A+B -> F A-B -> H + *B + * H B-A -> H AxB -> H + * + * rev = true (intersect) + * A + * F H + * F AxB -> F B-A -> F + *B + * H A-B -> F A+B -> H + * + * F/H = Fill outer / Hole outer + * A/B specify operands + * + = union, - = subtraction, x = intersection + * -> read as "produces" + * + * This is the main function of boolops, yet its operation isn't very complicated. + * It traverses the crossings, and uses the crossing direction to decide whether + * the next segment should be taken from A or from B. The second half of the + * function deals with figuring out what to do with bits that have no intersection. + */ +Shape shape_boolean(bool rev, Shape const & a, Shape const & b, CrossingSet const & crs) { + const Regions ac = a.content, bc = b.content; + + //Keep track of which crossings we've hit. + std::vector > visited; + for(unsigned i = 0; i < crs.size(); i++) + visited.push_back(std::vector(crs[i].size(), false)); + + //bool const exception = + + //Traverse the crossings, creating chunks + Regions chunks; + while(true) { + unsigned i, j; + first_false(visited, i, j); + if(i == visited.size()) break; + + Path res; + do { + Crossing cur = crs[i][j]; + visited[i][j] = true; + + //get indices of the dual: + unsigned io = cur.getOther(i), jo = find_crossing(crs[io], cur, io); + if(jo < visited[io].size()) visited[io][jo] = true; + + //main driving logic + if(logical_xor(cur.dir, rev)) { + if(i >= ac.size()) { i = io; j = jo; } + j++; + if(j >= crs[i].size()) j = 0; + Crossing next = crs[i][j]; + ac[next.a].boundary.appendPortionTo(res, cur.ta, next.ta); + } else { + if(i < ac.size()) { i = io; j = jo; } + j++; + if(j >= crs[i].size()) j = 0; + Crossing next = crs[i][j]; + bc[next.b - ac.size()].boundary.appendPortionTo(res, cur.tb, next.tb); + } + } while (!visited[i][j]); + if(res.size() > 0) chunks.push_back(Region(res)); + } + + //If true, then we are on the 'subtraction diagonal' + bool const on_sub = logical_xor(a.fill, b.fill); + //If true, then the hole must be inside the other to be included + bool const a_mode = logical_xor(logical_xor(!rev, a.fill), on_sub), + b_mode = logical_xor(logical_xor(!rev, b.fill), on_sub); + + //Handle unintersecting portions + for(unsigned i = 0; i < crs.size(); i++) { + if(crs[i].size() == 0) { + Region r(i < ac.size() ? ac[i] : bc[i - ac.size()]); + bool mode(i < ac.size() ? a_mode : b_mode); + + if(logical_xor(r.fill, i < ac.size() ? a.fill : b.fill)) { + //is an inner (fill is opposite the outside fill) + Point exemplar = r.boundary.initialPoint(); + Regions const & others = i < ac.size() ? bc : ac; + for(unsigned j = 0; j < others.size(); j++) { + if(others[j].contains(exemplar)) { + //contained in another + if(mode) chunks.push_back(r); + goto skip; + } + } + } + //disjoint + if(!mode) chunks.push_back(r); + skip: (void)0; + } + } + + return Shape(chunks); +} + +// Just a convenience wrapper for shape_boolean, which handles the crossings +Shape shape_boolean(bool rev, Shape const & a, Shape const & b) { + CrossingSet crs = crossings_between(a, b); + + return shape_boolean(rev, a, b, crs); +} + + +// Some utility functions for boolop: + +std::vector region_sizes(Shape const &a) { + std::vector ret; + for(unsigned i = 0; i < a.size(); i++) { + ret.push_back(double(a[i].size())); + } + return ret; +} + +Shape shape_boolean_ra(bool rev, Shape const &a, Shape const &b, CrossingSet const &crs) { + return shape_boolean(rev, a.inverse(), b, reverse_ta(crs, a.size(), region_sizes(a))); +} + +Shape shape_boolean_rb(bool rev, Shape const &a, Shape const &b, CrossingSet const &crs) { + return shape_boolean(rev, a, b.inverse(), reverse_tb(crs, a.size(), region_sizes(b))); +} + +/* This is a function based on shape_boolean which allows boolean operations + * to be specified as a logic table. This logic table is 4 bit-flags, which + * correspond to the elements of the 'truth table' for a particular operation. + * These flags are defined with the enums starting with BOOLOP_ . + */ +Shape boolop(Shape const &a, Shape const &b, unsigned flags, CrossingSet const &crs) { + flags &= 15; + if(flags <= BOOLOP_UNION) { + switch(flags) { + case BOOLOP_INTERSECT: return shape_boolean(true, a, b, crs); + case BOOLOP_SUBTRACT_A_B: return shape_boolean_rb(true, a, b, crs); + case BOOLOP_IDENTITY_A: return a; + case BOOLOP_SUBTRACT_B_A: return shape_boolean_ra(true, a, b, crs); + case BOOLOP_IDENTITY_B: return b; + case BOOLOP_EXCLUSION: { + Shape res = shape_boolean_rb(true, a, b, crs); + append(res.content, shape_boolean_ra(true, a, b, crs).content); + return res; + } + case BOOLOP_UNION: return shape_boolean(false, a, b); + } + } else { + switch(flags - BOOLOP_NEITHER) { + case BOOLOP_SUBTRACT_A_B: return shape_boolean_ra(false, a, b, crs); + case BOOLOP_SUBTRACT_B_A: return shape_boolean_rb(false, a, b, crs); + case BOOLOP_EXCLUSION: { + Shape res = shape_boolean_ra(false, a, b, CrossingSet(crs)); + append(res.content, shape_boolean_rb(false, a, b, CrossingSet(crs)).content); + return res; + } + } + return boolop(a, b, ~flags, crs).inverse(); + } + return Shape(); +} + +/* This version of the boolop function doesn't require a set of crossings, as + * it computes them for you. This is more efficient in some cases, as the + * shape can be inverted before finding crossings. In the special case of + * exclusion it uses the other version of boolop. + */ +Shape boolop(Shape const &a, Shape const &b, unsigned flags) { + flags &= 15; + if(flags <= BOOLOP_UNION) { + switch(flags) { + case BOOLOP_INTERSECT: return shape_boolean(true, a, b); + case BOOLOP_SUBTRACT_A_B: return shape_boolean(true, a, b.inverse()); + case BOOLOP_IDENTITY_A: return a; + case BOOLOP_SUBTRACT_B_A: return shape_boolean(true, b, a.inverse()); + case BOOLOP_IDENTITY_B: return b; + case BOOLOP_EXCLUSION: { + Shape res = shape_boolean(true, a, b.inverse()); + append(res.content, shape_boolean(true, b, a.inverse()).content); + return res; + } //return boolop(a, b, flags, crossings_between(a, b)); + case BOOLOP_UNION: return shape_boolean(false, a, b); + } + } else { + switch(flags - BOOLOP_NEITHER) { + case BOOLOP_SUBTRACT_A_B: return shape_boolean(false, b, a.inverse()); + case BOOLOP_SUBTRACT_B_A: return shape_boolean(false, a, b.inverse()); + case BOOLOP_EXCLUSION: { + Shape res = shape_boolean(false, a, b.inverse()); + append(res.content, shape_boolean(false, b, a.inverse()).content); + return res; + } //return boolop(a, b, flags, crossings_between(a, b)); + } + return boolop(a, b, ~flags).inverse(); + } + return Shape(); +} + + +int paths_winding(std::vector const &ps, Point p) { + int ret; + for(unsigned i = 0; i < ps.size(); i++) + ret += winding(ps[i], p); + return ret; +} + +std::vector y_of_roots(std::vector const & ps, double x) { + std::vector res; + for(unsigned i = 0; i < ps.size(); i++) { + std::vector temp = ps[i].roots(x, X); + for(unsigned i = 0; i < temp.size(); i++) + res.push_back(ps[i].valueAt(temp[i], Y)); + } + std::sort(res.begin(), res.end()); + return res; +} + +struct Edge { + unsigned ix; + double from, to; + bool rev; + int wind; + Edge(unsigned i, double ft, double tt, bool r, unsigned w) : ix(i), from(ft), to(tt), rev(r), wind(w) {} + Edge(unsigned i, double ft, double tt, bool r, std::vector const &ps) : ix(i), from(ft), to(tt), rev(r) { + //TODO: get the edge wind data some other way + Point p = ps[i].pointAt(ft); + std::vector rs = y_of_roots(ps, p[X]); + unsigned interv = std::lower_bound(rs.begin(), rs.end(), p[Y]) - rs.begin(); + wind = interv % 2; + } + double initial() { return rev ? to : from; } + double final() { return rev ? from : to; } + void addTo(Path &res, std::vector const &ps) { + if(rev) { + Path p = ps[ix].portion(to, from).reverse(); + for(unsigned i = 0; i < p.size(); i++) + res.append(p[i]); + } else { + ps[ix].appendPortionTo(res, from, to); + } + } +}; + +typedef std::vector Edges; + +double point_cosine(Point a, Point b, Point c) { + Point db = b - a, dc = c - a; + return cross(db, dc) / (db.length() * dc.length()); +} + +//sanitize +Regions regionize_paths(std::vector const &ps, bool evenodd) { + CrossingSet crs = crossings_among(ps); + + Edges es; + + for(unsigned i = 0; i < crs.size(); i++) { + for(unsigned j = 0; j < crs[i].size(); j++) { + Crossing cur = crs[i][j]; + int io = i, jo = j; + + jo++; + if(jo >= crs[io].size()) jo = 0; + //std::cout << io << ", " << jo << "\n"; + if(cur.a == i) + es.push_back(Edge(i, cur.ta, crs[io][jo].ta, false, ps)); + else + es.push_back(Edge(i, cur.tb, crs[io][jo].tb, false, ps)); + + jo-=2; + if(jo < 0) jo += crs[io].size(); + // std::cout << io << ", " << jo << "\n"; + if(cur.a == i) + es.push_back(Edge(i, crs[io][jo].ta, cur.ta, true, ps)); + else + es.push_back(Edge(i, crs[io][jo].tb, cur.tb, true, ps)); + } + } + for(unsigned i = 0; i used(es2.size(), false); + while(true) { + unsigned i = std::find(used.begin(), used.end(), false) - used.begin(); + if(i == used.size()) break; + Path res; + do { + es2[i].addTo(res, ps); + Point pnt = res.finalPoint(); + std::vector poss; + for(unsigned j = 0; j < es2.size(); j++) + if(near(pnt, ps[es2[j].ix].pointAt(es2[j].initial()))) poss.push_back(j); + if(poss.empty()) break; + unsigned best = 0; + if(poss.size() > 1) { + double crossval = 10; + Point along = ps[i].pointAt(es2[i].final()+0.1); + for(unsigned j = 0; j < poss.size(); j++) { + unsigned ix = poss[j]; + double val = point_cosine(pnt, along, ps[ix].pointAt(es2[ix].initial()+.01)); + if(val < crossval) { + crossval = val; + best = j; + } + } + } + i = poss[best]; + } while(!used[i]); + chunks.push_back(Region(res)); + } + return chunks; +} + +/* This transforms a shape by a matrix. In the case that the matrix flips + * the shape, it reverses the paths in order to preserve the fill. + */ +Shape Shape::operator*(Matrix const &m) const { + Shape ret; + for(unsigned i = 0; i < size(); i++) + ret.content.push_back(content[i] * m); + ret.fill = fill; + return ret; +} + +// Inverse is a boolean not, and simply reverses all the paths & fill flags +Shape Shape::inverse() const { + Shape ret; + for(unsigned i = 0; i < size(); i++) + ret.content.push_back(content[i].inverse()); + ret.fill = !fill; + return ret; +} + +struct ContainmentOrder { + std::vector const *rs; + explicit ContainmentOrder(std::vector const *r) : rs(r) {} + bool operator()(unsigned a, unsigned b) const { return (*rs)[b].contains((*rs)[a]); } +}; + +bool Shape::contains(Point const &p) const { + std::vector pnt; + pnt.push_back(Rect(p, p)); + std::vector > cull = sweep_bounds(pnt, bounds(*this)); + if(cull[0].size() == 0) return !fill; + return content[*min_element(cull[0].begin(), cull[0].end(), ContainmentOrder(&content))].isFill(); +} + +bool Shape::inside_invariants() const { //semi-slow & easy to violate + for(unsigned i = 0; i < size(); i++) + if( logical_xor(content[i].isFill(), contains(content[i].boundary.initialPoint())) ) return false; + return true; +} +bool Shape::region_invariants() const { //semi-slow + for(unsigned i = 0; i < size(); i++) + if(!content[i].invariants()) return false; + return true; +} +bool Shape::cross_invariants() const { //slow + CrossingSet crs; // = crossings_among(paths_from_regions(content)); + for(unsigned i = 0; i < crs.size(); i++) + if(!crs[i].empty()) return false; + return true; +} + +bool Shape::invariants() const { + return inside_invariants() && region_invariants() && cross_invariants(); +} + +} diff --git a/src/2geom/shape.h b/src/2geom/shape.h new file mode 100644 index 000000000..3700e9e5a --- /dev/null +++ b/src/2geom/shape.h @@ -0,0 +1,90 @@ +#ifndef __2GEOM_SHAPE_H +#define __2GEOM_SHAPE_H + +#include +#include + +#include "region.h" + +//TODO: BBOX optimizations + +namespace Geom { + +enum { + BOOLOP_JUST_A = 1, + BOOLOP_JUST_B = 2, + BOOLOP_BOTH = 4, + BOOLOP_NEITHER = 8 +}; + +enum { + BOOLOP_NULL = 0, + BOOLOP_INTERSECT = BOOLOP_BOTH, + BOOLOP_SUBTRACT_A_B = BOOLOP_JUST_B, + BOOLOP_IDENTITY_A = BOOLOP_JUST_A | BOOLOP_BOTH, + BOOLOP_SUBTRACT_B_A = BOOLOP_JUST_A, + BOOLOP_IDENTITY_B = BOOLOP_JUST_B | BOOLOP_BOTH, + BOOLOP_EXCLUSION = BOOLOP_JUST_A | BOOLOP_JUST_B, + BOOLOP_UNION = BOOLOP_JUST_A | BOOLOP_JUST_B | BOOLOP_BOTH +}; + +class Shape { + Regions content; + mutable bool fill; + //friend Shape shape_region_boolean(bool rev, Shape const & a, Region const & b); + friend CrossingSet crossings_between(Shape const &a, Shape const &b); + friend Shape shape_boolean(bool rev, Shape const &, Shape const &, CrossingSet const &); + friend Shape boolop(Shape const &a, Shape const &b, unsigned); + friend Shape boolop(Shape const &a, Shape const &b, unsigned, CrossingSet const &); + + public: + Shape() : fill(true) {} + explicit Shape(Region const & r) { + content = Regions(1, r); + fill = r.fill; + } + explicit Shape(Regions const & r) : content(r) { update_fill(); } + explicit Shape(bool f) : fill(f) {} + Shape(Regions const & r, bool f) : content(r), fill(f) {} + + Regions getContent() const { return content; } + bool isFill() const { return fill; } + + unsigned size() const { return content.size(); } + const Region &operator[](unsigned ix) const { return content[ix]; } + + Shape inverse() const; + Shape operator*(Matrix const &m) const; + + bool contains(Point const &p) const; + + bool inside_invariants() const; //semi-slow & easy to violate : checks that the insides are inside, the outsides are outside + bool region_invariants() const; //semi-slow : checks for self crossing + bool cross_invariants() const; //slow : checks that everything is disjoint + bool invariants() const; //vera slow (combo rombo, checks the above) + + private: + void update_fill() const { + unsigned ix = outer_index(content); + if(ix < size()) + fill = content[ix].fill; + else if(size() > 0) + fill = content.front().fill; + else + fill = true; + } +}; + +inline CrossingSet crossings_between(Shape const &a, Shape const &b) { return crossings(paths_from_regions(a.content), paths_from_regions(b.content)); } + +Shape shape_boolean(bool rev, Shape const &, Shape const &, CrossingSet const &); +Shape shape_boolean(bool rev, Shape const &, Shape const &); + +Shape boolop(Shape const &, Shape const &, unsigned flags); +Shape boolop(Shape const &, Shape const &, unsigned flags, CrossingSet &); + +Regions regionize_paths(std::vector const &ps, bool evenodd=true); + +} + +#endif diff --git a/src/2geom/solve-bezier-one-d.cpp b/src/2geom/solve-bezier-one-d.cpp index 558c10c0f..4f37cd049 100644 --- a/src/2geom/solve-bezier-one-d.cpp +++ b/src/2geom/solve-bezier-one-d.cpp @@ -35,7 +35,7 @@ const double BEPSILON = ldexp(1.0,-MAXDEPTH-1); /*Flatness control value */ * of the roots in the open interval (0, 1). Return the number of roots found. */ void -find_bernstein_roots(double *w, /* The control points */ +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 */ @@ -43,13 +43,11 @@ find_bernstein_roots(double *w, /* The control points */ { unsigned n_crossings = 0; /* Number of zero-crossings */ - double split = 0.5; int old_sign = SGN(w[0]); for (unsigned i = 1; i <= degree; i++) { int sign = SGN(w[i]); if (sign) { if (sign != old_sign && old_sign) { - split = (i-0.5)/degree; n_crossings++; } old_sign = sign; @@ -75,18 +73,16 @@ find_bernstein_roots(double *w, /* The control points */ const double Ax = right_t - left_t; const double Ay = w[degree] - w[0]; - solutions.push_back(left_t + Ax*w[0] / Ay); + solutions.push_back(left_t - Ax*w[0] / Ay); return; } break; } /* Otherwise, solve recursively after subdividing control polygon */ - - - /* This bit is very clever (if I say so myself) - rather than arbitrarily subdividing at the t = 0.5 point, we subdivide at the most likely point of change of direction. This buys us a factor of 10 speed up. We also avoid lots of stack frames by avoiding tail recursion. */ double Left[degree+1], /* New left and right */ Right[degree+1]; /* control polygons */ + const double split = 0.5; Bernstein(w, degree, split, Left, Right); double mid_t = left_t*(1-split) + right_t*split; @@ -100,73 +96,6 @@ find_bernstein_roots(double *w, /* The control points */ find_bernstein_roots(Right, degree, solutions, depth+1, mid_t, right_t); } -/* - * find_bernstein_roots : Given an equation in Bernstein-Bernstein form, find all - * of the roots in the interval [0, 1]. Return the number of roots found. - */ -void -find_bernstein_roots_buggy(double *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) -{ - unsigned n_crossings = 0; /* Number of zero-crossings */ - - int old_sign = SGN(w[0]); - for (unsigned i = 1; i <= degree; i++) { - int sign = SGN(w[i]); - if (sign) { - if (sign != old_sign && old_sign) - n_crossings++; - old_sign = sign; - } - } - - switch (n_crossings) { - case 0: /* No solutions here */ - return; - - case 1: - /* Unique solution */ - /* Stop recursion when the tree is deep enough */ - /* if deep enough, return 1 solution at midpoint */ - if (depth >= MAXDEPTH) { - solutions.push_back((left_t + right_t) / 2.0); - return; - } - - // I thought secant method would be faster here, but it'aint. -- njh - - if (control_poly_flat_enough(w, degree, left_t, right_t)) { - const double Ax = right_t - left_t; - const double Ay = w[degree] - w[0]; - - solutions.push_back(left_t + Ax*w[0] / Ay); - return; - } - break; - } - - /* Otherwise, solve recursively after subdividing control polygon */ - - double split = 0.5; - - double Left[degree+1], /* New left and right */ - Right[degree+1]; /* control polygons */ - Bernstein(w, degree, split, Left, Right); - - double mid_t = left_t*(1-split) + right_t*split; - - find_bernstein_roots_buggy(Left, degree, solutions, depth+1, left_t, mid_t); - - /* Solution is exactly on the subdivision point. */ - if (Right[0] == 0) - solutions.push_back(mid_t); - - find_bernstein_roots_buggy(Right, degree, solutions, depth+1, mid_t, right_t); -} - /* * control_poly_flat_enough : * Check if the control polygon of a Bernstein curve is flat enough diff --git a/src/2geom/solver.h b/src/2geom/solver.h index c2e290251..e0d66fc2a 100644 --- a/src/2geom/solver.h +++ b/src/2geom/solver.h @@ -5,6 +5,8 @@ namespace Geom{ + class Point; + unsigned crossing_count(Geom::Point const *V, /* Control pts of Bezier curve */ unsigned degree); /* Degree of Bezier curve */ @@ -21,14 +23,7 @@ crossing_count(double const *V, /* Control pts of Bezier curve */ double left_t, double right_t); void find_bernstein_roots( - double *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=0, double right_t=1); -void -find_bernstein_roots_buggy( - double *w, /* The control points */ + 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 */ diff --git a/src/2geom/svg-path-parser.cpp b/src/2geom/svg-path-parser.cpp index 0063fcfe9..0bd15e8b9 100644 --- a/src/2geom/svg-path-parser.cpp +++ b/src/2geom/svg-path-parser.cpp @@ -1,4 +1,4 @@ -#line 1 "/home/michael/2geom/src/svg-path-parser.rl" +#line 1 "/home/njh/svn/lib2geom/src/svg-path-parser.rl" /* * parse SVG path specifications * @@ -129,7 +129,7 @@ private: }; -#line 133 "/home/michael/2geom/src/svg-path-parser.cpp" +#line 133 "/home/njh/svn/lib2geom/src/svg-path-parser.cpp" static const char _svg_path_actions[] = { 0, 1, 0, 1, 1, 1, 2, 1, 3, 1, 4, 1, 5, 1, 15, 1, @@ -905,8 +905,8 @@ static const short _svg_path_indicies[] = { 185, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 194, 179, 180, 205, 0, 575, - 575, 576, 576, 577, 575, 578, 0, 761, - 761, 762, 762, 763, 761, 764, 0, 781, + 575, 576, 576, 577, 575, 578, 0, 769, + 769, 770, 770, 771, 769, 772, 0, 781, 781, 782, 782, 783, 781, 784, 0, 750, 741, 0, 714, 0, 699, 699, 701, 702, 715, 715, 699, 700, 714, 0, 738, 738, @@ -939,8 +939,8 @@ static const short _svg_path_indicies[] = { 294, 295, 295, 297, 320, 300, 301, 303, 304, 305, 306, 307, 308, 309, 310, 311, 312, 313, 314, 315, 316, 317, 318, 319, - 310, 295, 296, 321, 0, 757, 757, 758, - 758, 759, 757, 760, 0, 777, 777, 778, + 310, 295, 296, 321, 0, 765, 765, 766, + 766, 767, 765, 768, 0, 777, 777, 778, 778, 779, 777, 780, 0, 749, 740, 0, 712, 0, 694, 694, 696, 697, 713, 713, 694, 695, 712, 0, 737, 737, 728, 730, @@ -1054,9 +1054,9 @@ static const short _svg_path_indicies[] = { 123, 125, 148, 128, 129, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 138, 123, - 124, 149, 0, 769, 769, 770, 770, 771, - 769, 772, 0, 765, 765, 766, 766, 767, - 765, 768, 0, 599, 599, 600, 600, 601, + 124, 149, 0, 761, 761, 762, 762, 763, + 761, 764, 0, 757, 757, 758, 758, 759, + 757, 760, 0, 599, 599, 600, 600, 601, 599, 602, 0, 607, 607, 608, 608, 609, 607, 610, 0, 208, 209, 209, 211, 629, 214, 215, 628, 217, 218, 219, 220, 221, @@ -1345,9 +1345,9 @@ static const unsigned char _svg_path_trans_actions_wi[] = { 0, 1, 1, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, - 17, 3, 17, 0, 0, 9, 59, 59, - 59, 9, 59, 59, 59, 11, 62, 62, - 62, 11, 62, 62, 62, 0, 1, 1, + 17, 3, 17, 0, 0, 11, 62, 62, + 62, 11, 62, 62, 62, 9, 59, 59, + 59, 9, 59, 59, 59, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0 }; @@ -1356,7 +1356,7 @@ static const int svg_path_start = 0; static const int svg_path_first_final = 326; -#line 133 "/home/michael/2geom/src/svg-path-parser.rl" +#line 133 "/home/njh/svn/lib2geom/src/svg-path-parser.rl" void Parser::parse(char const *str) @@ -1369,7 +1369,7 @@ throw(SVGPathParseError) _reset(); -#line 1373 "/home/michael/2geom/src/svg-path-parser.cpp" +#line 1373 "/home/njh/svn/lib2geom/src/svg-path-parser.cpp" { cs = svg_path_start; } @@ -1445,13 +1445,13 @@ _match: switch ( *_acts++ ) { case 0: -#line 145 "/home/michael/2geom/src/svg-path-parser.rl" +#line 145 "/home/njh/svn/lib2geom/src/svg-path-parser.rl" { start = p; } break; case 1: -#line 149 "/home/michael/2geom/src/svg-path-parser.rl" +#line 149 "/home/njh/svn/lib2geom/src/svg-path-parser.rl" { char const *end=p; std::string buf(start, end); @@ -1460,55 +1460,55 @@ _match: } break; case 2: -#line 156 "/home/michael/2geom/src/svg-path-parser.rl" +#line 156 "/home/njh/svn/lib2geom/src/svg-path-parser.rl" { _push(1.0); } break; case 3: -#line 160 "/home/michael/2geom/src/svg-path-parser.rl" +#line 160 "/home/njh/svn/lib2geom/src/svg-path-parser.rl" { _push(0.0); } break; case 4: -#line 164 "/home/michael/2geom/src/svg-path-parser.rl" +#line 164 "/home/njh/svn/lib2geom/src/svg-path-parser.rl" { _absolute = true; } break; case 5: -#line 168 "/home/michael/2geom/src/svg-path-parser.rl" +#line 168 "/home/njh/svn/lib2geom/src/svg-path-parser.rl" { _absolute = false; } break; case 6: -#line 172 "/home/michael/2geom/src/svg-path-parser.rl" +#line 172 "/home/njh/svn/lib2geom/src/svg-path-parser.rl" { _moveTo(_pop_point()); } break; case 7: -#line 176 "/home/michael/2geom/src/svg-path-parser.rl" +#line 176 "/home/njh/svn/lib2geom/src/svg-path-parser.rl" { _lineTo(_pop_point()); } break; case 8: -#line 180 "/home/michael/2geom/src/svg-path-parser.rl" +#line 180 "/home/njh/svn/lib2geom/src/svg-path-parser.rl" { _lineTo(Point(_pop_coord(X), _current[Y])); } break; case 9: -#line 184 "/home/michael/2geom/src/svg-path-parser.rl" +#line 184 "/home/njh/svn/lib2geom/src/svg-path-parser.rl" { _lineTo(Point(_current[X], _pop_coord(Y))); } break; case 10: -#line 188 "/home/michael/2geom/src/svg-path-parser.rl" +#line 188 "/home/njh/svn/lib2geom/src/svg-path-parser.rl" { Point p = _pop_point(); Point c1 = _pop_point(); @@ -1517,7 +1517,7 @@ _match: } break; case 11: -#line 195 "/home/michael/2geom/src/svg-path-parser.rl" +#line 195 "/home/njh/svn/lib2geom/src/svg-path-parser.rl" { Point p = _pop_point(); Point c1 = _pop_point(); @@ -1525,7 +1525,7 @@ _match: } break; case 12: -#line 201 "/home/michael/2geom/src/svg-path-parser.rl" +#line 201 "/home/njh/svn/lib2geom/src/svg-path-parser.rl" { Point p = _pop_point(); Point c = _pop_point(); @@ -1533,14 +1533,14 @@ _match: } break; case 13: -#line 207 "/home/michael/2geom/src/svg-path-parser.rl" +#line 207 "/home/njh/svn/lib2geom/src/svg-path-parser.rl" { Point p = _pop_point(); _quadTo(_quad_tangent, p); } break; case 14: -#line 212 "/home/michael/2geom/src/svg-path-parser.rl" +#line 212 "/home/njh/svn/lib2geom/src/svg-path-parser.rl" { Point point = _pop_point(); bool sweep = _pop_flag(); @@ -1553,16 +1553,16 @@ _match: } break; case 15: -#line 223 "/home/michael/2geom/src/svg-path-parser.rl" +#line 223 "/home/njh/svn/lib2geom/src/svg-path-parser.rl" { _closePath(); } break; case 16: -#line 360 "/home/michael/2geom/src/svg-path-parser.rl" +#line 360 "/home/njh/svn/lib2geom/src/svg-path-parser.rl" {goto _out;} break; -#line 1566 "/home/michael/2geom/src/svg-path-parser.cpp" +#line 1566 "/home/njh/svn/lib2geom/src/svg-path-parser.cpp" } } @@ -1571,7 +1571,7 @@ _again: goto _resume; _out: {} } -#line 370 "/home/michael/2geom/src/svg-path-parser.rl" +#line 370 "/home/njh/svn/lib2geom/src/svg-path-parser.rl" if ( cs < svg_path_first_final ) { diff --git a/src/2geom/svg-path.cpp b/src/2geom/svg-path.cpp index 141ddbcf3..312db9d23 100644 --- a/src/2geom/svg-path.cpp +++ b/src/2geom/svg-path.cpp @@ -34,7 +34,7 @@ namespace Geom { void output(Curve const &curve, SVGPathSink &sink) { - std::vector pts = sbasis_to_bezier(curve.sbasis(), 2); //TODO: use something better! + std::vector pts = sbasis_to_bezier(curve.toSBasis(), 2); //TODO: use something better! sink.curveTo(pts[0], pts[1], pts[2]); } diff --git a/src/2geom/svg-path.h b/src/2geom/svg-path.h index d09002220..ebc23a9b5 100644 --- a/src/2geom/svg-path.h +++ b/src/2geom/svg-path.h @@ -46,6 +46,7 @@ public: bool large_arc, bool sweep, Point p) = 0; virtual void closePath() = 0; virtual void finish() = 0; + virtual ~SVGPathSink() {} }; void output_svg_path(Path &path, SVGPathSink &sink); @@ -66,14 +67,14 @@ public: _path.appendNew(p); } - void curveTo(Point c0, Point c1, Point p) { - _path.appendNew(c0, c1, p); - } - void quadTo(Point c, Point p) { _path.appendNew(c, p); } + void curveTo(Point c0, Point c1, Point p) { + _path.appendNew(c0, c1, p); + } + void arcTo(double rx, double ry, double angle, bool large_arc, bool sweep, Point p) { diff --git a/src/2geom/sweep.cpp b/src/2geom/sweep.cpp new file mode 100644 index 000000000..4329bc559 --- /dev/null +++ b/src/2geom/sweep.cpp @@ -0,0 +1,104 @@ +#include "sweep.h" + +#include + +namespace Geom { + +std::vector > sweep_bounds(std::vector rs) { + std::vector events; events.reserve(rs.size()*2); + 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)); + } + std::sort(events.begin(), events.end()); + + std::vector open; + for(unsigned i = 0; i < events.size(); i++) { + unsigned ix = events[i].ix; + if(events[i].closing) { + std::vector::iterator iter = std::find(open.begin(), open.end(), ix); + //if(iter != open.end()) + open.erase(iter); + } else { + for(unsigned j = 0; j < open.size(); j++) { + unsigned jx = open[j]; + if(rs[jx][Y].intersects(rs[ix][Y])) { + pairs[jx].push_back(ix); + } + } + open.push_back(ix); + } + } + return pairs; +} + +std::vector > sweep_bounds(std::vector a, std::vector b) { + std::vector > pairs(a.size()); + if(a.empty() || b.empty()) return pairs; + std::vector events[2]; + events[0].reserve(a.size()*2); + events[1].reserve(b.size()*2); + + for(unsigned n = 0; n < 2; n++) { + 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)); + } + std::sort(events[n].begin(), events[n].end()); + } + + std::vector open[2]; + bool n = events[1].front() < events[0].front(); + for(unsigned i[] = {0,0}; i[n] < events[n].size();) { + unsigned ix = events[n][i[n]].ix; + bool closing = events[n][i[n]].closing; + //std::cout << n << "[" << ix << "] - " << (closing ? "closer" : "opener") << "\n"; + if(closing) { + open[n].erase(std::find(open[n].begin(), open[n].end(), ix)); + } else { + if(n) { + //n = 1 + //opening a B, add to all open a + for(unsigned j = 0; j < open[0].size(); j++) { + unsigned jx = open[0][j]; + if(a[jx][Y].intersects(b[ix][Y])) { + pairs[jx].push_back(ix); + } + } + } else { + //n = 0 + //opening an A, add all open b + for(unsigned j = 0; j < open[1].size(); j++) { + unsigned jx = open[1][j]; + if(b[jx][Y].intersects(a[ix][Y])) { + pairs[ix].push_back(jx); + } + } + } + open[n].push_back(ix); + } + i[n]++; + n = (events[!n][i[!n]] < events[n][i[n]]) ? !n : n; + } + return pairs; +} + +//Fake cull, until the switch to the real sweep is made. +std::vector > fake_cull(unsigned a, unsigned b) { + std::vector > ret; + + std::vector all; + for(unsigned j = 0; j < b; j++) + all.push_back(j); + + for(unsigned i = 0; i < a; i++) + ret.push_back(all); + + return ret; +} + +} diff --git a/src/2geom/sweep.h b/src/2geom/sweep.h new file mode 100644 index 000000000..9587cec36 --- /dev/null +++ b/src/2geom/sweep.h @@ -0,0 +1,29 @@ +#ifndef __2GEOM_SWEEP_H__ +#define __2GEOM_SWEEP_H__ + +#include +#include "d2.h" + +namespace Geom { + +struct Event { + double x; + unsigned ix; + bool closing; + Event(double pos, unsigned i, bool c) : x(pos), ix(i), closing(c) {} +// Lexicographic ordering by x then closing + bool operator<(Event const &other) const { + if(x < other.x) return true; + if(x > other.x) return false; + return closing < other.closing; + } + +}; +std::vector > sweep_bounds(std::vector); +std::vector > sweep_bounds(std::vector, std::vector); + +std::vector > fake_cull(unsigned a, unsigned b); + +} + +#endif diff --git a/src/2geom/transforms.h b/src/2geom/transforms.h index 1cd1d3e5f..94058bff0 100644 --- a/src/2geom/transforms.h +++ b/src/2geom/transforms.h @@ -28,7 +28,7 @@ class Translate { public: explicit Translate(Point const &p) : vec(p) {} explicit Translate(Coord const x, Coord const y) : vec(x, y) {} - inline operator Matrix() const { return Matrix(0, 0, 0, 0, vec[X], vec[Y]); } + inline operator Matrix() const { return Matrix(1, 0, 0, 1, vec[X], vec[Y]); } inline Coord operator[](Dim2 const dim) const { return vec[dim]; } inline Coord operator[](unsigned const dim) const { return vec[dim]; } diff --git a/src/2geom/utils.h b/src/2geom/utils.h index a2f906ff4..ca880640d 100644 --- a/src/2geom/utils.h +++ b/src/2geom/utils.h @@ -38,6 +38,9 @@ public: NotImplemented() : std::logic_error("method not implemented") {} }; +// proper logical xor +inline bool logical_xor (bool a, bool b) { return (a || b) && !(a && b); } + /** Sign function - indicates the sign of a numeric type. -1 indicates negative, 1 indicates * positive, and 0 indicates, well, 0. Mathsy people will know this is basically the derivative * of abs, except for the fact that it is defined on 0. diff --git a/src/live_effects/n-art-bpath-2geom.cpp b/src/live_effects/n-art-bpath-2geom.cpp index f04c80f2b..9e5b966ea 100644 --- a/src/live_effects/n-art-bpath-2geom.cpp +++ b/src/live_effects/n-art-bpath-2geom.cpp @@ -45,7 +45,7 @@ static void curve_to_svgd(std::ostream & f, Geom::Curve const* c) { // } else { //this case handles sbasis as well as all other curve types - Geom::Path sbasis_path = path_from_sbasis(c->sbasis(), 0.1); + Geom::Path sbasis_path = Geom::path_from_sbasis(c->toSBasis(), 0.1); //recurse to convert the new path resulting from the sbasis to svgd for(Geom::Path::iterator iter = sbasis_path.begin(); iter != sbasis_path.end(); ++iter) {