![[tokkee]](http://tokkee.org/images/avatar.png){kind=avatar}
summary | shortlog | log | commit | commitdiff | tree
raw | patch | inline | side by side (parent: cfc4305)
raw | patch | inline | side by side (parent: cfc4305)
| author | johanengelen <johanengelen@users.sourceforge.net> | |
| Thu, 30 Aug 2007 18:32:36 +0000 (18:32 +0000) | ||
| committer | johanengelen <johanengelen@users.sourceforge.net> | |
| Thu, 30 Aug 2007 18:32:36 +0000 (18:32 +0000) |
45 files changed:
index a100c418dfecf71116388d6dac8dee58e25a19a0..4a87c1a93d3a21a69bce3f0115d71516ec149c3d 100644 (file)
2geom/bezier-to-sbasis.h \\r
2geom/bezier-utils.cpp \\r
2geom/bezier-utils.h \\r
+ 2geom/chebyshev.cpp \\r
+ 2geom/chebyshev.h \\r
2geom/choose.h \\r
2geom/circle-circle.cpp \\r
2geom/circulator.h \\r
2geom/convex-cover.cpp \\r
2geom/convex-cover.h \\r
2geom/coord.h \\r
+ 2geom/crossing.cpp \\r
+ 2geom/crossing.cpp \\r
+ 2geom/d2-sbasis.cpp \\r
+ 2geom/d2-sbasis.h \\r
2geom/d2.cpp \\r
2geom/d2.h \\r
2geom/geom.cpp \\r
2geom/linear.h \\r
2geom/matrix.cpp \\r
2geom/matrix.h \\r
+ 2geom/ord.h \\r
+ 2geom/path-intersection.cpp \\r
+ 2geom/path-intersection.h \\r
2geom/path.cpp \\r
2geom/path.h \\r
2geom/piecewise.cpp \\r
2geom/quadtree.cpp \\r
2geom/quadtree.h \\r
2geom/rect.h \\r
+ 2geom/region.cpp \\r
+ 2geom/region.h \\r
2geom/sbasis-2d.cpp \\r
2geom/sbasis-2d.h \\r
2geom/sbasis.cpp \\r
2geom/sbasis-roots.cpp \\r
2geom/sbasis-to-bezier.cpp \\r
2geom/sbasis-to-bezier.h \\r
+ 2geom/shape.cpp \\r
+ 2geom/shape.h \\r
2geom/solve-bezier-one-d.cpp \\r
2geom/solve-bezier-parametric.cpp \\r
2geom/solver.h \\r
2geom/svg-path.h \\r
2geom/svg-path-parser.cpp \\r
2geom/svg-path-parser.h \\r
+ 2geom/sweep.cpp \\r
+ 2geom/sweep.h \\r
2geom/transforms.cpp \\r
2geom/transforms.h \\r
2geom/utils.h \\r
index a5b82702356f006a91a895b3c3eca2c1bc4ea297..694760d5abcdc9b000f96c8a1fa4954d7434bc1f 100644 (file)
* 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 7949e006c957752ded62cc98d0fc24ffc6a6697e..2c1b49213acd4fd40b73b24eca18054ecac81b6c 100644 (file)
--- a/src/2geom/bezier.h
+++ b/src/2geom/bezier.h
*
* Copyright 2007 MenTaLguY <mental@rydia.net>
* Copyright 2007 Michael Sloan <mgsloan@gmail.com>
+ * Copyright 2007 Nathan Hurst <njh@njhurst.com>
*
* This library is free software; you can redistribute it and/or
* modify it either under the terms of the GNU Lesser General Public
#include "coord.h"
#include "isnan.h"
#include "bezier-to-sbasis.h"
+#include "d2.h"
+#include "solver.h"
+#include <boost/optional/optional.hpp>
namespace Geom {
+template<unsigned order>
+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 <unsigned order>
class Bezier {
private:
- Coord c_[order+1];
+ Coord c_[order+1];
+
+ template<unsigned ord>
+ friend Bezier<ord> portion(const Bezier<ord> & a, Coord from, Coord to);
+
+ template<unsigned ord>
+ friend Interval bounds_fast(Bezier<ord> const & b);
+
+ template<unsigned ord>
+ friend Bezier<ord-1> derivative(const Bezier<ord> & 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 <unsigned required_order>
- static void assert_order(Bezier<required_order> 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 <unsigned required_order>
+ //static void assert_order(Bezier<required_order> 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<order>(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<order>(t, c_, NULL, NULL); }
+ inline Coord operator()(double t) const { return valueAt(t); }
- Maybe<int> winding(Point p) const {
- return sbasis_winding(toSBasis(), p);
- }
+ inline SBasis toSBasis() const { return bezier_to_sbasis<order>(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<Coord> valueAndDerivatives(Coord t, unsigned n_derivs) const {
+ throw NotImplemented();
+ // Can't be implemented without a dynamic version of subdivide.
+ /*std::vector<Coord> val_n_der;
+ Coord d_[order+1];
+ for(int di = n_derivs; di > 0; di--) {
+ val_n_der.push_back(subdivideArr<di>(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<order>, Bezier<order> > subdivide(Coord t) const {
+ Bezier<order> a, b;
+ subdivideArr(t, order, c_, a.c_, b.c_);
+ return std::pair<Bezier<order>, Bezier<order> >(a, b);
+ }
+
+ std::vector<double> roots() const {
+ std::vector<double> solutions;
+ find_bernstein_roots(c_, order, solutions, 0, 0.0, 1.0);
+ return solutions;
+ }
};
+//TODO: implement others
+template<unsigned order>
+Bezier<order> operator+(const Bezier<order> & a, double v) {
+ Bezier<order> result;
+ for(unsigned i = 0; i <= order; i++)
+ result[i] = a[i] + v;
+ return result;
+}
+template<unsigned order>
+Bezier<order> operator-(const Bezier<order> & a, double v) {
+ Bezier<order> result;
+ for(unsigned i = 0; i <= order; i++)
+ result[i] = a[i] - v;
+ return result;
+}
+template<unsigned order>
+Bezier<order> operator*(const Bezier<order> & a, double v) {
+ Bezier<order> result;
+ for(unsigned i = 0; i <= order; i++)
+ result[i] = a[i] * v;
+ return result;
+}
+template<unsigned order>
+Bezier<order> operator/(const Bezier<order> & a, double v) {
+ Bezier<order> result;
+ for(unsigned i = 0; i <= order; i++)
+ result[i] = a[i] / v;
+ return result;
+}
+
template<unsigned order>
Bezier<order> reverse(const Bezier<order> & a) {
- Bezier<order> result;
- for(int i = 0; i <= order; i++)
- result[i] = a[order - i];
- return result;
+ Bezier<order> result;
+ for(unsigned i = 0; i <= order; i++)
+ result[i] = a[order - i];
+ return result;
+}
+
+template<unsigned order>
+Bezier<order> portion(const Bezier<order> & a, double from, double to) {
+ //TODO: implement better?
+ Coord res[order+1];
+ if(from == 0) {
+ if(to == 1) { return Bezier<order>(a); }
+ subdivideArr<order>(to, a.c_, res, NULL);
+ return Bezier<order>(res);
+ }
+ subdivideArr<order>(from, a.c_, NULL, res);
+ if(to == 1) return Bezier<order>(res);
+ Coord res2[order+1];
+ subdivideArr<order>((to - from)/(1 - from), res, res2, NULL);
+ return Bezier<order>(res2);
+}
+
+template<unsigned order>
+std::vector<Point> bezier_points(const D2<Bezier<order> > & a) {
+ std::vector<Point> 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<unsigned order>
+Bezier<order-1> derivative(const Bezier<order> & a) {
+ Bezier<order-1> der;
+
+ for(unsigned i = 0; i < order; i++) {
+ der.c_[i] = order*(a.c_[i+1] - a.c_[i]);
+ }
+ return der;
}
template<unsigned order>
-vector<Point> bezier_points(const D2<Bezier<order> > & a) {
- vector<Point> 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<order> const & b) {
+ return Interval::fromArray(b.c_, order+1);
+}
+
+//TODO: better bounds exact
+template<unsigned order>
+inline Interval bounds_exact(Bezier<order> const & b) {
+ return bounds_exact(b.toSBasis());
+}
+
+template<unsigned order>
+inline Interval bounds_local(Bezier<order> 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
--- /dev/null
+++ b/src/2geom/chebyshev.cpp
@@ -0,0 +1,126 @@
+#include "chebyshev.h"
+
+#include "sbasis.h"
+#include "sbasis-poly.h"
+
+#include <vector>
+using std::vector;
+
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_chebyshev.h>
+
+namespace Geom{
+
+SBasis cheb(unsigned n) {
+ static std::vector<SBasis> 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<SBasis> 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
--- /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 4c7f83baddeb1bf432d5dd604e5db9edb6623df2..512aac39f3cdd781bbd1097a4a8cb0b92d0a07c2 100644 (file)
--- a/src/2geom/circulator.h
+++ b/src/2geom/circulator.h
template <typename Iterator>
class Circulator {
public:
- typedef std::random_access_iterator_tag std::iterator_category;
- typedef std::iterator_traits<Iterator>::value_type value_type;
- typedef std::iterator_traits<Iterator>::difference_type difference_type;
- typedef std::iterator_traits<Iterator>::pointer pointer;
- typedef std::iterator_traits<Iterator>::reference reference;
+ typedef std::random_access_iterator_tag iterator_category;
+ typedef typename std::iterator_traits<Iterator>::value_type value_type;
+ typedef typename std::iterator_traits<Iterator>::difference_type difference_type;
+ typedef typename std::iterator_traits<Iterator>::pointer pointer;
+ typedef typename std::iterator_traits<Iterator>::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 {
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 42192c4458178444b6970aab08d3c6a40704652e..bd020ab6625d9fab3ae7bcc7fe83cc4a120dc47e 100644 (file)
--- a/src/2geom/concepts.h
+++ b/src/2geom/concepts.h
#include "sbasis.h"
#include "interval.h"
#include "point.h"
-
+#include <vector>
#include <boost/concept_check.hpp>
namespace Geom {
bool b;
BoundsType i;
Interval dom;
+ std::vector<OutputType> v;
+ unsigned u;
+ SbType sb;
void constraints() {
t = T(o);
b = t.isZero();
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
--- /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<double> 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<double> 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<double> 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<double> 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<Crossings> i(&cr_a); !i.ended(); i++) {
+ const Crossing cur = *i;
+ Eraser<Crossings> 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
--- /dev/null
+++ b/src/2geom/crossing.h
@@ -0,0 +1,103 @@
+#ifndef __GEOM_CROSSING_H
+#define __GEOM_CROSSING_H
+
+#include <vector>
+#include <set>
+#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<Crossing> Crossings;
+typedef std::vector<Crossings> CrossingSet;
+
+template<typename C>
+std::vector<Rect> bounds(C const &a) {
+ std::vector<Rect> 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<typename T>
+struct Crosser {
+ virtual Crossings crossings(T const &a, T const &b) { return crossings(std::vector<T>(1,a), std::vector<T>(1,b))[0]; }
+ virtual CrossingSet crossings(std::vector<T> const &a, std::vector<T> const &b) {
+ CrossingSet results(a.size() + b.size(), Crossings());
+
+ std::vector<std::vector<unsigned> > 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<double> max);
+Crossings reverse_tb(Crossings const &cr, unsigned split, std::vector<double> max);
+CrossingSet reverse_ta(CrossingSet const &cr, unsigned split, std::vector<double> max);
+CrossingSet reverse_tb(CrossingSet const &cr, unsigned split, std::vector<double> 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
--- /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 {
+\r
+SBasis L2(D2<SBasis> const & a, unsigned k) { return sqrt(dot(a, a), k); }\r
+\r
+D2<SBasis> multiply(Linear const & a, D2<SBasis> const & b) {\r
+ return D2<SBasis>(multiply(a, b[X]), multiply(a, b[Y]));\r
+}\r
+\r
+D2<SBasis> multiply(SBasis const & a, D2<SBasis> const & b) {\r
+ return D2<SBasis>(multiply(a, b[X]), multiply(a, b[Y]));\r
+}\r
+\r
+D2<SBasis> truncate(D2<SBasis> const & a, unsigned terms) {\r
+ return D2<SBasis>(truncate(a[X], terms), truncate(a[Y], terms));\r
+}\r
+\r
+unsigned sbasis_size(D2<SBasis> const & a) {\r
+ return std::max((unsigned) a[0].size(), (unsigned) a[1].size());\r
+}\r
+\r
+//TODO: Is this sensical? shouldn't it be like pythagorean or something?\r
+double tail_error(D2<SBasis> const & a, unsigned tail) {\r
+ return std::max(a[0].tailError(tail), a[1].tailError(tail));\r
+}\r
+\r
+Piecewise<D2<SBasis> > sectionize(D2<Piecewise<SBasis> > const &a) {\r
+ Piecewise<SBasis> x = partition(a[0], a[1].cuts), y = partition(a[1], a[0].cuts);\r
+ assert(x.size() == y.size());\r
+ Piecewise<D2<SBasis> > ret;\r
+ for(unsigned i = 0; i < x.size(); i++)\r
+ ret.push_seg(D2<SBasis>(x[i], y[i]));\r
+ ret.cuts.insert(ret.cuts.end(), x.cuts.begin(), x.cuts.end());\r
+ return ret;\r
+}\r
+\r
+D2<Piecewise<SBasis> > make_cuts_independant(Piecewise<D2<SBasis> > const &a) {\r
+ D2<Piecewise<SBasis> > ret;\r
+ for(unsigned d = 0; d < 2; d++) {\r
+ for(unsigned i = 0; i < a.size(); i++)\r
+ ret[d].push_seg(a[i][d]);\r
+ ret[d].cuts.insert(ret[d].cuts.end(), a.cuts.begin(), a.cuts.end());\r
+ }\r
+ return ret;\r
+}\r
+\r
+Piecewise<D2<SBasis> > rot90(Piecewise<D2<SBasis> > const &M){\r
+ Piecewise<D2<SBasis> > result;\r
+ if (M.empty()) return M;\r
+ result.push_cut(M.cuts[0]);\r
+ for (unsigned i=0; i<M.size(); i++){\r
+ result.push(rot90(M[i]),M.cuts[i+1]);\r
+ }\r
+ return result;\r
+}\r
+\r
+Piecewise<SBasis> dot(Piecewise<D2<SBasis> > const &a, \r
+ Piecewise<D2<SBasis> > const &b){\r
+ Piecewise<SBasis > result;\r
+ if (a.empty() || b.empty()) return result;\r
+ Piecewise<D2<SBasis> > aa = partition(a,b.cuts);\r
+ Piecewise<D2<SBasis> > bb = partition(b,a.cuts);\r
+\r
+ result.push_cut(aa.cuts.front());\r
+ for (unsigned i=0; i<aa.size(); i++){\r
+ result.push(dot(aa.segs[i],bb.segs[i]),aa.cuts[i+1]);\r
+ }\r
+ return result;\r
+}\r
+\r
+Piecewise<SBasis> cross(Piecewise<D2<SBasis> > const &a, \r
+ Piecewise<D2<SBasis> > const &b){\r
+ Piecewise<SBasis > result;\r
+ if (a.empty() || b.empty()) return result;\r
+ Piecewise<D2<SBasis> > aa = partition(a,b.cuts);\r
+ Piecewise<D2<SBasis> > bb = partition(b,a.cuts);\r
+\r
+ result.push_cut(aa.cuts.front());\r
+ for (unsigned i=0; i<a.size(); i++){\r
+ result.push(cross(aa.segs[i],bb.segs[i]),aa.cuts[i+1]);\r
+ }\r
+ return result;\r
+}\r
+\r
+/* Replaced by remove_short_cuts in piecewise.h
+//this recursively removes the shortest cut interval until none is shorter than tol.\r
+//TODO: code this in a more efficient way!\r
+Piecewise<D2<SBasis> > remove_short_cuts(Piecewise<D2<SBasis> > const &f, double tol){\r
+ double min = tol;\r
+ unsigned idx = f.size();\r
+ for(unsigned i=0; i<f.size(); i++){\r
+ if (min > f.cuts[i+1]-f.cuts[i]){\r
+ min = f.cuts[i+1]-f.cuts[i];\r
+ idx = int(i);\r
+ }\r
+ }\r
+ if (idx==f.size()){\r
+ return f;\r
+ }\r
+ if (f.size()==1) {\r
+ //removing this seg would result in an empty pw<d2<sb>>...\r
+ return f;\r
+ }\r
+ Piecewise<D2<SBasis> > new_f=f;\r
+ for (int dim=0; dim<2; dim++){\r
+ double v = Hat(f.segs.at(idx)[dim][0]);\r
+ //TODO: what about closed curves?\r
+ if (idx>0 && f.segs.at(idx-1).at1()==f.segs.at(idx).at0()) \r
+ new_f.segs.at(idx-1)[dim][0][1] = v;\r
+ if (idx<f.size() && f.segs.at(idx+1).at0()==f.segs.at(idx).at1()) \r
+ new_f.segs.at(idx+1)[dim][0][0] = v;\r
+ }\r
+ double t = (f.cuts.at(idx)+f.cuts.at(idx+1))/2;\r
+ new_f.cuts.at(idx+1) = t; \r
+ \r
+ new_f.segs.erase(new_f.segs.begin()+idx);\r
+ new_f.cuts.erase(new_f.cuts.begin()+idx); \r
+ return remove_short_cuts(new_f, tol);\r
+}\r
+*/
+\r
+//if tol>0, only force continuity where the jump is smaller than tol.\r
+Piecewise<D2<SBasis> > force_continuity(Piecewise<D2<SBasis> > const &f, \r
+ double tol,\r
+ bool closed){\r
+ if (f.size()==0) return f;\r
+ Piecewise<D2<SBasis> > result=f;\r
+ unsigned cur = (closed)? 0:1;\r
+ unsigned prev = (closed)? f.size()-1:0;\r
+ while(cur<f.size()){\r
+ Point pt0 = f.segs[prev].at1();\r
+ Point pt1 = f.segs[cur ].at0();\r
+ if (tol<=0 || L2sq(pt0-pt1)<tol*tol){\r
+ pt0 = (pt0+pt1)/2;\r
+ for (unsigned dim=0; dim<2; dim++){\r
+ result.segs[prev][dim][0][1]=pt0[dim];\r
+ result.segs[cur ][dim][0][0]=pt0[dim];\r
+ }\r
+ }\r
+ prev = cur++;\r
+ }\r
+ return result;\r
+}\r\r
+}
diff --git a/src/2geom/d2-sbasis.h b/src/2geom/d2-sbasis.h
--- /dev/null
+++ b/src/2geom/d2-sbasis.h
@@ -0,0 +1,88 @@
+#ifdef _2GEOM_D2 /*This is intentional: we don't actually want anyone to
+ 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_SBASIS_CURVE_H
+#define __2GEOM_SBASIS_CURVE_H
+
+#include "sbasis.h"\r
+#include "sbasis-2d.h"\r
+#include "piecewise.h"
+
+//TODO: implement intersect
+
+namespace Geom {
+\r
+inline D2<SBasis> compose(D2<SBasis> const & a, SBasis const & b) {\r
+ return D2<SBasis>(compose(a[X], b), compose(a[Y], b));\r
+}\r
+\r
+SBasis L2(D2<SBasis> const & a, unsigned k);\r
+double L2(D2<double> const & a);\r
+\r
+D2<SBasis> multiply(Linear const & a, D2<SBasis> const & b);\r
+inline D2<SBasis> operator*(Linear const & a, D2<SBasis> const & b) { return multiply(a, b); }\r
+D2<SBasis> multiply(SBasis const & a, D2<SBasis> const & b);\r
+inline D2<SBasis> operator*(SBasis const & a, D2<SBasis> const & b) { return multiply(a, b); }\r
+D2<SBasis> truncate(D2<SBasis> const & a, unsigned terms);\r
+\r
+unsigned sbasis_size(D2<SBasis> const & a);\r
+double tail_error(D2<SBasis> const & a, unsigned tail);\r
+\r
+//Piecewise<D2<SBasis> > specific decls:\r
+\r
+Piecewise<D2<SBasis> > sectionize(D2<Piecewise<SBasis> > const &a);\r
+D2<Piecewise<SBasis> > make_cuts_independant(Piecewise<D2<SBasis> > const &a);\r
+Piecewise<D2<SBasis> > rot90(Piecewise<D2<SBasis> > const &a);\r
+Piecewise<SBasis> dot(Piecewise<D2<SBasis> > const &a, Piecewise<D2<SBasis> > const &b);\r
+Piecewise<SBasis> cross(Piecewise<D2<SBasis> > const &a, Piecewise<D2<SBasis> > const &b);\r
+\r
+Piecewise<D2<SBasis> > force_continuity(Piecewise<D2<SBasis> > const &f, \r
+ double tol=0,\r
+ bool closed=false);\r
+\r
+class CoordIterator\r
+: public std::iterator<std::input_iterator_tag, SBasis const>\r
+{\r
+public:\r
+ CoordIterator(std::vector<D2<SBasis> >::const_iterator const &iter, unsigned d) : impl_(iter), ix_(d) {}\r
+\r
+ inline bool operator==(CoordIterator const &other) { return other.impl_ == impl_; }\r
+ inline bool operator!=(CoordIterator const &other) { return other.impl_ != impl_; }\r
+\r
+ inline SBasis operator*() const {\r
+ return (*impl_)[ix_];\r
+ }\r
+\r
+ inline CoordIterator &operator++() {\r
+ ++impl_;\r
+ return *this;\r
+ }\r
+ inline CoordIterator operator++(int) {\r
+ CoordIterator old=*this;\r
+ ++(*this);\r
+ return old;\r
+ }\r
+\r
+private:\r
+ std::vector<D2<SBasis> >::const_iterator impl_;\r
+ unsigned ix_;\r
+};\r
+\r
+inline CoordIterator iterateCoord(Piecewise<D2<SBasis> > const &a, unsigned d) {\r
+ return CoordIterator(a.segs.begin(), d);\r
+}\r
+\r
+//bounds specializations with order\r
+inline Rect bounds_fast(D2<SBasis> const & s, unsigned order=0) {\r
+ return Rect(bounds_fast(s[X], order),\r
+ bounds_fast(s[Y], order));\r
+}\r
+inline Rect bounds_local(D2<SBasis> const & s, Interval i, unsigned order=0) {\r
+ return Rect(bounds_local(s[X], i, order),\r
+ bounds_local(s[Y], i, order));\r
+}
+
+}
+
+#endif
+#endif
diff --git a/src/2geom/d2.h b/src/2geom/d2.h
index 82070519c349f2955b4436e9865a0053a2101244..696ad0191411b04f5ddc5f2b7e1ef79cb61b7111 100644 (file)
--- a/src/2geom/d2.h
+++ b/src/2geom/d2.h
Point valueAt(double t) const {\r
boost::function_requires<FragmentConcept<T> >();\r
return (*this)(t);\r
+ }
+ std::vector<Point > valueAndDerivatives(double t, unsigned count) const {
+ std::vector<Coord> x = f[X].valueAndDerivatives(t, count),
+ y = f[Y].valueAndDerivatives(t, count);
+ std::vector<Point> res;
+ for(unsigned i = 0; i < count; i++) {
+ res.push_back(Point(x[i], y[i]));
+ }
+ return res;
}\r
D2<SBasis> toSBasis() const {\r
boost::function_requires<FragmentConcept<T> >();\r
\r
Point operator()(double t) const;\r
Point operator()(double x, double y) const;\r
-};\r
-\r
+};\r\r
template <typename T>\r
-D2<T> reverse(const D2<T> &a) {\r
+inline D2<T> reverse(const D2<T> &a) {\r
boost::function_requires<FragmentConcept<T> >();\r
return D2<T>(reverse(a[X]), reverse(a[Y]));\r
}\r
+
+template <typename T>
+inline D2<T> portion(const D2<T> &a, Coord f, Coord t) {
+ boost::function_requires<FragmentConcept<T> >();
+ return D2<T>(portion(a[X], f, t), portion(a[Y], f, t));
+}
\r
//IMPL: boost::EqualityComparableConcept\r
template <typename T>\r
template<typename T>\r
D2<T> derivative(D2<T> const & a) {\r
return D2<T>(derivative(a[X]), derivative(a[Y]));\r
-}\r
-\r
+}\r\r
template<typename T>\r
D2<T> integral(D2<T> const & a) {\r
return D2<T>(integral(a[X]), integral(a[Y]));\r
-}\r
+}
\r
} //end namespace Geom\r
\r
-\r
-\r
-//TODO: implement intersect\r
-\r
-\r
#include "rect.h"\r
-#include "sbasis.h"\r
-#include "sbasis-2d.h"\r
-#include "piecewise.h"\r
+#include "d2-sbasis.h"\r
\r
namespace Geom{\r
\r
Rect bounds_local(const D2<T> &a, const Interval &t) {\r
boost::function_requires<FragmentConcept<T> >(); \r
return Rect(bounds_local(a[X], t), bounds_local(a[Y], t));\r
-}\r
-\r
-//D2<SBasis> specific decls:\r
-\r
-inline D2<SBasis> compose(D2<SBasis> const & a, SBasis const & b) {\r
- return D2<SBasis>(compose(a[X], b), compose(a[Y], b));\r
-}\r
-\r
-SBasis L2(D2<SBasis> const & a, unsigned k);\r
-double L2(D2<double> const & a);\r
-\r
-inline D2<SBasis> portion(D2<SBasis> const &M, double t0, double t1){\r
- return D2<SBasis>(portion(M[0],t0,t1),portion(M[1],t0,t1));\r
-}\r
-\r
-D2<SBasis> multiply(Linear const & a, D2<SBasis> const & b);\r
-inline D2<SBasis> operator*(Linear const & a, D2<SBasis> const & b) { return multiply(a, b); }\r
-D2<SBasis> multiply(SBasis const & a, D2<SBasis> const & b);\r
-inline D2<SBasis> operator*(SBasis const & a, D2<SBasis> const & b) { return multiply(a, b); }\r
-D2<SBasis> truncate(D2<SBasis> const & a, unsigned terms);\r
-\r
-unsigned sbasis_size(D2<SBasis> const & a);\r
-double tail_error(D2<SBasis> const & a, unsigned tail);\r
-\r
-//Piecewise<D2<SBasis> > specific decls:\r
-\r
-Piecewise<D2<SBasis> > sectionize(D2<Piecewise<SBasis> > const &a);\r
-D2<Piecewise<SBasis> > make_cuts_independant(Piecewise<D2<SBasis> > const &a);\r
-Piecewise<D2<SBasis> > rot90(Piecewise<D2<SBasis> > const &a);\r
-Piecewise<SBasis> dot(Piecewise<D2<SBasis> > const &a, Piecewise<D2<SBasis> > const &b);\r
-Piecewise<SBasis> cross(Piecewise<D2<SBasis> > const &a, Piecewise<D2<SBasis> > const &b);\r
-\r
-Piecewise<D2<SBasis> > force_continuity(Piecewise<D2<SBasis> > const &f, \r
- double tol=0,\r
- bool closed=false);\r
-\r
-class CoordIterator\r
-: public std::iterator<std::input_iterator_tag, SBasis const>\r
-{\r
-public:\r
- CoordIterator(std::vector<D2<SBasis> >::const_iterator const &iter, unsigned d) : impl_(iter), ix_(d) {}\r
-\r
- inline bool operator==(CoordIterator const &other) { return other.impl_ == impl_; }\r
- inline bool operator!=(CoordIterator const &other) { return other.impl_ != impl_; }\r
-\r
- inline SBasis operator*() const {\r
- return (*impl_)[ix_];\r
- }\r
-\r
- inline CoordIterator &operator++() {\r
- ++impl_;\r
- return *this;\r
- }\r
- inline CoordIterator operator++(int) {\r
- CoordIterator old=*this;\r
- ++(*this);\r
- return old;\r
- }\r
-\r
-private:\r
- std::vector<D2<SBasis> >::const_iterator impl_;\r
- unsigned ix_;\r
-};\r
-\r
-inline CoordIterator iterateCoord(Piecewise<D2<SBasis> > const &a, unsigned d) {\r
- return CoordIterator(a.segs.begin(), d);\r
-}\r
-\r
-//bounds specializations with order\r
-inline Rect bounds_fast(D2<SBasis> const & s, unsigned order=0) {\r
- return Rect(bounds_fast(s[X], order),\r
- bounds_fast(s[Y], order));\r
-}\r
-inline Rect bounds_local(D2<SBasis> const & s, Interval i, unsigned order=0) {\r
- return Rect(bounds_local(s[X], i, order),\r
- bounds_local(s[Y], i, order));\r
-}\r
+}\r\r
};\r
\r
/*\r
diff --git a/src/2geom/interval.h b/src/2geom/interval.h
index 459f2cd495e8f337a191c508d24568daf4180318..19e08978cf343c3947b43b098b4264f94bd2575b 100644 (file)
--- a/src/2geom/interval.h
+++ b/src/2geom/interval.h
_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];
}
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
_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 f537959434c26c7bfa905f34fe10127a011bcbac..c0e64ad33a0525f31b923faa6ef216f4252b8886 100644 (file)
--- a/src/2geom/matrix.cpp
+++ b/src/2geom/matrix.cpp
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;
}
}
_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) &&
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 6a45b50c84283d1e1f9d407976abccbb9c666b4e..c39c997166a51ef4a1a6a2c71f4b6abddecfe061 100644 (file)
--- a/src/2geom/matrix.h
+++ b/src/2geom/matrix.h
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
--- /dev/null
+++ b/src/2geom/ord.h
@@ -0,0 +1,37 @@
+
+#ifndef __2GEOM_ORD__
+#define __2GEOM_ORD__
+
+namespace {
+
+enum Cmp {\r
+ LESS_THAN=-1,\r
+ GREATER_THAN=1,\r
+ EQUAL_TO=0\r
+};\r
+
+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;
+ }
+}
+\r
+template <typename T1, typename T2>\r
+inline Cmp cmp(T1 const &a, T2 const &b) {\r
+ if ( a < b ) {\r
+ return LESS_THAN;\r
+ } else if ( b < a ) {\r
+ return GREATER_THAN;\r
+ } else {\r
+ return EQUAL_TO;\r
+ }\r
+}
+
+}
+
+#endif
diff --git a/src/2geom/path-intersection.cpp b/src/2geom/path-intersection.cpp
--- /dev/null
@@ -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<double> 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<D2<SBasis> > 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<typename T>
+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<double> curve_mono_splits(Curve const &d) {
+ std::vector<double> 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<double> offset_doubles(std::vector<double> const &x, double offs) {
+ std::vector<double> 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<double> path_mono_splits(Path const &p) {
+ std::vector<double> 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<double> 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<std::vector<double> > paths_mono_splits(std::vector<Path> const &ps) {
+ std::vector<std::vector<double> > 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<std::vector<Rect> > split_bounds(std::vector<Path> const &p, std::vector<std::vector<double> > splits) {
+ std::vector<std::vector<Rect> > ret;
+ for(unsigned i = 0; i < p.size(); i++) {
+ std::vector<Rect> 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<Path> const &a, std::vector<Path> const &b) {
+ if(b.empty()) return CrossingSet(a.size(), Crossings());
+ CrossingSet results(a.size() + b.size(), Crossings());
+ if(a.empty()) return results;
+
+ std::vector<std::vector<double> > splits_a = paths_mono_splits(a), splits_b = paths_mono_splits(b);
+ std::vector<std::vector<Rect> > bounds_a = split_bounds(a, splits_a), bounds_b = split_bounds(b, splits_b);
+
+ std::vector<Rect> 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<std::vector<unsigned> > 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<std::vector<unsigned> > 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<Path> const &p) {
+ CrossingSet results(p.size(), Crossings());
+ if(p.empty()) return results;
+
+ std::vector<std::vector<double> > splits = paths_mono_splits(p);
+ std::vector<std::vector<Rect> > prs = split_bounds(p, splits);
+ std::vector<Rect> rs;
+ for(unsigned i = 0; i < prs.size(); i++) rs.push_back(union_list(prs[i]));
+
+ std::vector<std::vector<unsigned> > 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<std::vector<unsigned> > 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<double> 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<std::vector<double> > curves_mono_splits(Path const &p) {
+ std::vector<std::vector<double> > ret;
+ for(unsigned i = 0; i <= p.size(); i++) {
+ std::vector<double> spl;
+ spl.push_back(0);
+ append(spl, curve_mono_splits(p[i]));
+ spl.push_back(1);
+ ret.push_back(spl);
+ }
+}
+
+std::vector<std::vector<Rect> > curves_split_bounds(Path const &p, std::vector<std::vector<double> > splits) {
+ std::vector<std::vector<Rect> > ret;
+ for(unsigned i = 0; i < splits.size(); i++) {
+ std::vector<Rect> 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<std::vector<unsigned> > cull = sweep_bounds(bounds(p));
+ std::vector<std::vector<double> > spl = curves_mono_splits(p);
+ std::vector<std::vector<Rect> > 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<std::vector<unsigned> > 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<std::vector<unsigned> > 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<Path> const &p) {
+ CrossingSet results(p.size(), Crossings());
+ if(p.empty()) return results;
+
+ SimpleCrosser cc;
+
+ std::vector<std::vector<unsigned> > 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
--- /dev/null
@@ -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<typename T>
+Crossings curve_sweep(Path const &a, Path const &b) {
+ T t;
+ Crossings ret;
+ std::vector<Rect> bounds_a = bounds(a), bounds_b = bounds(b);
+ std::vector<std::vector<unsigned> > ixs = sweep_bounds(bounds_a, bounds_b);
+ for(unsigned i = 0; i < a.size(); i++) {
+ for(std::vector<unsigned>::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<Path> {
+ Crossings crossings(Curve const &a, Curve const &b);
+ Crossings crossings(Path const &a, Path const &b) { return curve_sweep<SimpleCrosser>(a, b); }
+ CrossingSet crossings(std::vector<Path> const &a, std::vector<Path> const &b) { return Crosser<Path>::crossings(a, b); }
+};
+
+struct MonoCrosser : public Crosser<Path> {
+ Crossings crossings(Path const &a, Path const &b) { return crossings(std::vector<Path>(1,a), std::vector<Path>(1,b))[0]; }
+ CrossingSet crossings(std::vector<Path> const &a, std::vector<Path> const &b);
+};
+
+typedef SimpleCrosser DefaultCrosser;
+
+std::vector<double> path_mono_splits(Path const &p);
+
+CrossingSet crossings_among(std::vector<Path> const & p);
+inline Crossings self_crossings(Path const & a) { return crossings_among(std::vector<Path>(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<Path> const & a, std::vector<Path> const & b) {
+ DefaultCrosser c = DefaultCrosser();
+ return c.crossings(a, b);
+}
+
+}
+
+#endif
diff --git a/src/2geom/path.cpp b/src/2geom/path.cpp
index 91868eb7e349d09ad1bae2b9233384ff0846f6c5..f05d3b0cf7153c5d8c0cb5b81f0ca0199a4321d9 100644 (file)
--- a/src/2geom/path.cpp
+++ b/src/2geom/path.cpp
-/*\r
- * Path - Series of continuous curves\r
- * \r
- * Copyright 2007 MenTaLguY <mental@rydia.net>\r
- * \r
- * This library is free software; you can redistribute it and/or \r
- * modify it either under the terms of the GNU Lesser General Public\r
- * License version 2.1 as published by the Free Software Foundation\r
- * (the "LGPL") or, at your option, under the terms of the Mozilla\r
- * Public License Version 1.1 (the "MPL"). If you do not alter this \r
- * notice, a recipient may use your version of this file under either\r
- * the MPL or the LGPL.\r
- * \r
- * You should have received a copy of the LGPL along with this library\r
- * in the file COPYING-LGPL-2.1; if not, write to the Free Software \r
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
- * You should have received a copy of the MPL along with this library \r
- * in the file COPYING-MPL-1.1\r
- * \r
- * The contents of this file are subject to the Mozilla Public License\r
- * Version 1.1 (the "License"); you may not use this file except in\r
- * compliance with the License. You may obtain a copy of the License at\r
- * http://www.mozilla.org/MPL/\r
- * \r
- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY\r
- * OF ANY KIND, either express or implied. See the LGPL or the MPL for\r
- * the specific language governing rights and limitations.\r
- */\r
-\r
-#include "path.h"\r
-\r
-namespace Geom {\r
-\r
-namespace {\r
-\r
-enum Cmp {\r
- LESS_THAN=-1,\r
- GREATER_THAN=1,\r
- EQUAL_TO=0\r
-};\r
-\r
-template <typename T1, typename T2>\r
-inline Cmp cmp(T1 const &a, T2 const &b) {\r
- if ( a < b ) {\r
- return LESS_THAN;\r
- } else if ( b < a ) {\r
- return GREATER_THAN;\r
- } else {\r
- return EQUAL_TO;\r
- }\r
-}\r
-\r
-}\r
-\r
-boost::optional<int> CurveHelpers::sbasis_winding(D2<SBasis> const &sb, Point p) {\r
- Interval ix = bounds_fast(sb[X]);\r
-\r
- if ( p[X] > ix.max() ) { /* ray does not intersect bbox */\r
- return 0;\r
- }\r
-\r
- SBasis fy = sb[Y];\r
- fy -= p[Y];\r
-\r
- if (fy.empty()) { /* coincident horizontal segment */\r
- return boost::optional<int>();\r
- }\r
-\r
- if ( p[X] < ix.min() ) { /* ray does not originate in bbox */\r
- double y = p[Y];\r
- /* winding determined by position of endpoints */\r
- Cmp initial_to_ray = cmp(fy[0][0], y);\r
- Cmp final_to_ray = cmp(fy[0][1], y);\r
- switch (cmp(final_to_ray, initial_to_ray)) {\r
- case GREATER_THAN:\r
- /* exclude final endpoint */\r
- return ( final_to_ray != EQUAL_TO );\r
- case LESS_THAN:\r
- /* exclude final endpoint */\r
- return -( final_to_ray != EQUAL_TO );\r
- default:\r
- /* any intersections cancel out */\r
- return 0;\r
- }\r
- } else { /* ray originates in bbox */\r
- std::vector<double> ts = roots(fy);\r
-\r
- static const unsigned MAX_DERIVATIVES=8;\r
- boost::optional<SBasis> ds[MAX_DERIVATIVES];\r
- ds[0] = derivative(fy);\r
-\r
- /* winding determined by summing signs of derivatives at intersections */\r
- int winding=0;\r
- for ( std::vector<double>::iterator ti = ts.begin()\r
- ; ti != ts.end()\r
- ; ++ti )\r
- { \r
- double t = *ti;\r
- if ( sb[X](t) >= p[X] ) { /* root is ray intersection */\r
- for ( boost::optional<SBasis> *di = ds\r
- ; di != ( ds + MAX_DERIVATIVES )\r
- ; ++di )\r
- {\r
- if (!*di) {\r
- *di = derivative(**(di-1));\r
- }\r
- switch (cmp((**di)(t), 0)) {\r
- case GREATER_THAN:\r
- if ( t < 1 ) { /* exclude final endpoint */\r
- winding += 1;\r
- }\r
- goto next_root;\r
- case LESS_THAN:\r
- if ( t < 1 ) { /* exclude final endpoint */\r
- winding -= 1;\r
- }\r
- goto next_root;\r
- default: (void)0;\r
- /* give up */\r
- };\r
- }\r
- } \r
-next_root: (void)0;\r
- }\r
- \r
- return winding;\r
- }\r
-}\r
-\r
-Rect BezierHelpers::bounds(unsigned degree, Point const *points) {\r
- Point min=points[0];\r
- Point max=points[0];\r
- for ( unsigned i = 1 ; i <= degree ; ++i ) {\r
- for ( unsigned axis = 0 ; axis < 2 ; ++axis ) {\r
- min[axis] = std::min(min[axis], points[i][axis]);\r
- max[axis] = std::max(max[axis], points[i][axis]);\r
- }\r
- }\r
- return Rect(min, max);\r
-}\r
-\r
-Point BezierHelpers::point_and_derivatives_at(Coord t,\r
- unsigned degree,\r
- Point const *points,\r
- unsigned n_derivs,\r
- Point *derivs)\r
-{\r
- return Point(0,0); // TODO\r
-}\r
-\r
-Geom::Point\r
-BezierHelpers::subdivideArr(Coord t, // Parameter value\r
- unsigned degree, // Degree of bezier curve\r
- Geom::Point const *V, // Control pts\r
- Geom::Point *Left, // RETURN left half ctl pts\r
- Geom::Point *Right) // RETURN right half ctl pts\r
-{\r
- Geom::Point Vtemp[degree+1][degree+1];\r
-\r
- /* Copy control points */\r
- std::copy(V, V+degree+1, Vtemp[0]);\r
-\r
- /* Triangle computation */\r
- for (unsigned i = 1; i <= degree; i++) { \r
- for (unsigned j = 0; j <= degree - i; j++) {\r
- Vtemp[i][j] = lerp(t, Vtemp[i-1][j], Vtemp[i-1][j+1]);\r
- }\r
- }\r
- \r
- for (unsigned j = 0; j <= degree; j++)\r
- Left[j] = Vtemp[j][0];\r
- for (unsigned j = 0; j <= degree; j++)\r
- Right[j] = Vtemp[degree-j][j];\r
-\r
- return (Vtemp[degree][0]);\r
-}\r
-\r
-void Path::swap(Path &other) {\r
- std::swap(curves_, other.curves_);\r
- std::swap(closed_, other.closed_);\r
- std::swap(*final_, *other.final_);\r
- curves_[curves_.size()-1] = final_;\r
- other.curves_[other.curves_.size()-1] = other.final_;\r
-}\r
-\r
-Rect Path::bounds_fast() const {\r
- Rect bounds=front().bounds_fast();\r
- for ( const_iterator iter=++begin(); iter != end() ; ++iter ) {\r
- bounds.unionWith(iter->bounds_fast());\r
- }\r
- return bounds;\r
-}\r
-\r
-Rect Path::bounds_exact() const {\r
- Rect bounds=front().bounds_exact();\r
- for ( const_iterator iter=++begin(); iter != end() ; ++iter ) {\r
- bounds.unionWith(iter->bounds_exact());\r
- }\r
- return bounds;\r
-}\r
-\r
-int Path::winding(Point p) const {\r
- int winding = 0;\r
- boost::optional<Cmp> ignore = boost::optional<Cmp>();\r
- for ( const_iterator iter = begin()\r
- ; iter != end_closed()\r
- ; ++iter )\r
- {\r
- boost::optional<int> w = iter->winding(p);\r
- if (w) {\r
- winding += *w;\r
- ignore = boost::optional<Cmp>();\r
- } else {\r
- Point initial = iter->initialPoint();\r
- Point final = iter->finalPoint();\r
- switch (cmp(initial[X], final[X])) {\r
- case GREATER_THAN:\r
- if ( !ignore || *ignore != GREATER_THAN ) { /* ignore repeated */\r
- winding += 1;\r
- ignore = GREATER_THAN;\r
- }\r
- break;\r
- case LESS_THAN:\r
- if ( !ignore || *ignore != LESS_THAN ) { /* ignore repeated */\r
- if ( p[X] < final[X] ) { /* ignore final point */\r
- winding -= 1;\r
- ignore = LESS_THAN;\r
- }\r
- }\r
- break;\r
- case EQUAL_TO:\r
- /* always ignore null segments */\r
- break;\r
- }\r
- }\r
- }\r
- return winding;\r
-}\r
-\r
-void Path::append(Curve const &curve) {\r
- if ( curves_.front() != final_ && curve.initialPoint() != (*final_)[0] ) {\r
- throw ContinuityError();\r
- }\r
- do_append(curve.duplicate());\r
-}\r
-\r
-void Path::append(D2<SBasis> const &curve) {\r
- if ( curves_.front() != final_ ) {\r
- for ( int i = 0 ; i < 2 ; ++i ) {\r
- if ( curve[i][0][0] != (*final_)[0][i] ) {\r
- throw ContinuityError();\r
- }\r
- }\r
- }\r
- do_append(new SBasisCurve(curve));\r
-}\r
-\r
-void Path::do_update(Sequence::iterator first_replaced,\r
- Sequence::iterator last_replaced,\r
- Sequence::iterator first,\r
- Sequence::iterator last)\r
-{\r
- // note: modifies the contents of [first,last)\r
-\r
- check_continuity(first_replaced, last_replaced, first, last);\r
- delete_range(first_replaced, last_replaced);\r
- if ( ( last - first ) == ( last_replaced - first_replaced ) ) {\r
- std::copy(first, last, first_replaced);\r
- } else {\r
- // this approach depends on std::vector's behavior WRT iterator stability\r
- curves_.erase(first_replaced, last_replaced);\r
- curves_.insert(first_replaced, first, last);\r
- }\r
-\r
- if ( curves_.front() != final_ ) {\r
- (*final_)[0] = back().finalPoint();\r
- (*final_)[1] = front().initialPoint();\r
- }\r
-}\r
-\r
-void Path::do_append(Curve *curve) {\r
- if ( curves_.front() == final_ ) {\r
- (*final_)[1] = curve->initialPoint();\r
- }\r
- curves_.insert(curves_.end()-1, curve);\r
- (*final_)[0] = curve->finalPoint();\r
-}\r
-\r
-void Path::delete_range(Sequence::iterator first, Sequence::iterator last) {\r
- for ( Sequence::iterator iter=first ; iter != last ; ++iter ) {\r
- delete *iter;\r
- }\r
-}\r
-\r
-void Path::check_continuity(Sequence::iterator first_replaced,\r
- Sequence::iterator last_replaced,\r
- Sequence::iterator first,\r
- Sequence::iterator last)\r
-{\r
- if ( first != last ) {\r
- if ( first_replaced != curves_.begin() ) {\r
- if ( (*first_replaced)->initialPoint() != (*first)->initialPoint() ) {\r
- throw ContinuityError();\r
- }\r
- }\r
- if ( last_replaced != (curves_.end()-1) ) {\r
- if ( (*(last_replaced-1))->finalPoint() != (*(last-1))->finalPoint() ) {\r
- throw ContinuityError();\r
- }\r
- }\r
- } else if ( first_replaced != last_replaced && first_replaced != curves_.begin() && last_replaced != curves_.end()-1) {\r
- if ( (*first_replaced)->initialPoint() !=\r
- (*(last_replaced-1))->finalPoint() )\r
- {\r
- throw ContinuityError();\r
- }\r
- }\r
-}\r
-\r
-Rect SBasisCurve::bounds_fast() const {\r
- throw NotImplemented();\r
- return Rect(Point(0,0), Point(0,0));\r
-}\r
-\r
-Rect SBasisCurve::bounds_exact() const {\r
- throw NotImplemented();\r
- return Rect(Point(0,0), Point(0,0));\r
-}\r
-\r
-Point SBasisCurve::pointAndDerivativesAt(Coord t, unsigned n_derivs, Point *derivs) const {\r
- throw NotImplemented();\r
- return Point(0,0);\r
-}\r
-\r
-Path const &SBasisCurve::subdivide(Coord t, Path &out) const {\r
- throw NotImplemented();\r
-}\r
-\r
-Rect SVGEllipticalArc::bounds_fast() const {\r
- throw NotImplemented();\r
-}\r
-Rect SVGEllipticalArc::bounds_exact() const {\r
- throw NotImplemented();\r
-}\r
-\r
-Point SVGEllipticalArc::pointAndDerivativesAt(Coord t, unsigned n_derivs, Point *derivs) const {\r
- throw NotImplemented();\r
-}\r
-\r
-D2<SBasis> SVGEllipticalArc::sbasis() const {\r
- throw NotImplemented();\r
-}\r
-\r
-}\r
-\r
-/*\r
- Local Variables:\r
- mode:c++\r
- c-file-style:"stroustrup"\r
- c-file-offsets:((innamespace . 0)(substatement-open . 0))\r
- indent-tabs-mode:nil\r
- c-brace-offset:0\r
- fill-column:99\r
- End:\r
- vim: filetype=cpp:expandtab:shiftwidth=2:tabstop=8:softtabstop=2 :\r
-*/\r
-\r
+/*
+ * Path - Series of continuous curves
+ *
+ * Copyright 2007 MenTaLguY <mental@rydia.net>
+ *
+ * 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<double> 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<double>::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<double>::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<typename iter>
+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<SBasis> 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<Point> SVGEllipticalArc::pointAndDerivatives(Coord t, unsigned n) const {
+ throw NotImplemented();
+}
+
+std::vector<double> SVGEllipticalArc::roots(double v, Dim2 d) const {
+ //throw NotImplemented();
+}
+
+D2<SBasis> SVGEllipticalArc::toSBasis() const {
+ return D2<SBasis>(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 31a7173b79c7119b44e20c9388a75d7850e236b4..6f39eb7ef7c1b59437fbb3836260c522da74ac39 100644 (file)
--- a/src/2geom/path.h
+++ b/src/2geom/path.h
#include <algorithm>
#include <exception>
#include <stdexcept>
-#include <boost/optional/optional.hpp>
#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() {}
virtual Curve *duplicate() const = 0;
- virtual Rect bounds_fast() const = 0;
- virtual Rect bounds_exact() const = 0;
-
- virtual boost::optional<int> 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<double> 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> 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<Point> pointAndDerivatives(Coord t, unsigned n) const = 0;
+ virtual D2<SBasis> toSBasis() const = 0;
};
-struct CurveHelpers {
-protected:
- static boost::optional<int> sbasis_winding(D2<SBasis> const &sbasis, Point p);
-};
+class SBasisCurve : public Curve {
+private:
+ SBasisCurve();
+ D2<SBasis> inner;
+public:
+ explicit SBasisCurve(D2<SBasis> 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<Point> 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<double> 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<SBasis> toSBasis() const { return inner; }
};
-template <unsigned bezier_degree>
-class Bezier : public Curve, private CurveHelpers, private BezierHelpers {
+template <unsigned order>
+class BezierCurve : public Curve {
+private:
+ D2<Bezier<order> > inner;
public:
template <unsigned required_degree>
- static void assert_degree(Bezier<required_degree> const *) {}
+ static void assert_degree(BezierCurve<required_degree> const *) {}
+
+ BezierCurve() {}
- Bezier() {}
+ explicit BezierCurve(D2<Bezier<order> > const &x) : inner(x) {}
+
+ BezierCurve(Bezier<order> x, Bezier<order> 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<order>(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<order>(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<order>(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<int> 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<double>
+ 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<Point> ps) {
+ for(unsigned i = 0; i <= order; i++) {
+ setPoint(i, ps[i]);
+ }
+ }
+ std::vector<Point> points() const { return bezier_points(inner); }
+
+ std::pair<BezierCurve<order>, BezierCurve<order> > subdivide(Coord t) const {
+ std::pair<Bezier<order>, Bezier<order> > sx = inner[X].subdivide(t), sy = inner[Y].subdivide(t);
+ return std::pair<BezierCurve<order>, BezierCurve<order> >(
+ BezierCurve<order>(sx.first, sy.first),
+ BezierCurve<order>(sx.second, sy.second));
+ }
+
+ Curve *portion(double f, double t) const {
+ return new BezierCurve(Geom::portion(inner, f, t));
}
- D2<SBasis> sbasis() const {
- return handles_to_sbasis<bezier_degree>(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<Point> 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<order-1>(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<Point> pointAndDerivatives(Coord t, unsigned n) const { return inner.valueAndDerivatives(t, n); }
+
+ double valueAt(double t, Dim2 d) const { return inner[d].valueAt(t); }
+
+ D2<SBasis> 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<order>(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() {}
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<int> 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<double> roots(double v, Dim2 d) const;
- Path const &subdivide(Coord t, Path &out) const;
+ inline std::pair<SVGEllipticalArc, SVGEllipticalArc>
+ subdivide(Coord t) {
+ SVGEllipticalArc a(*this), b(*this);
+ a.final_ = b.initial_ = pointAt(t);
+ return std::pair<SVGEllipticalArc, SVGEllipticalArc>(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<Point> pointAndDerivatives(Coord t, unsigned n) const;
- D2<SBasis> sbasis() const;
+ D2<SBasis> toSBasis() const;
private:
Point initial_;
Point final_;
};
-class SBasisCurve : public Curve, private CurveHelpers {
-private:
- SBasisCurve();
-public:
- explicit SBasisCurve(D2<SBasis> 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<int> 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> sbasis() const { return coeffs_; }
-
-private:
- D2<SBasis> coeffs_;
-};
-
template <typename IteratorImpl>
class BaseIterator
: public std::iterator<std::forward_iterator_tag, Curve const>
++impl_;
return *this;
}
+
BaseIterator operator++(int) {
BaseIterator old=*this;
++(*this);
}
Curve *operator*() const { return (*impl_)->duplicate(); }
-
+
DuplicatingIterator &operator++() {
++impl_;
return *this;
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<D2<SBasis> > toPwSb() const {
Piecewise<D2<SBasis> > 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<double> roots(double v, Dim2 d) const {
+ std::vector<double> res;
+ for(const_iterator it = begin(); it != end_closed(); ++it) {
+ std::vector<double> 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 {
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<SBasis> const &curve);
template <typename CurveType, typename A>
return ret;
}
-template <unsigned bezier_degree>
-inline Path const &Bezier<bezier_degree>::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<Path> 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<double> 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<Point> pointAndDerivatives(Coord t, unsigned n) const = 0;
+ virtual D2<SBasis> toSBasis() const = 0;
+
};
+*/
}
index dde230dc49c2d0e1fcec27ae5694065ec9b464be..136687ae87c8fa518080f5945240e65652284521 100644 (file)
-#include "poly-dk-solve.h"
-#include <iterator>
-
-/*** implementation of the Durand-Kerner method. seems buggy*/
-
-std::complex<double> evalu(Poly const & p, std::complex<double> x) {
- std::complex<double> result = 0;
- std::complex<double> xx = 1;
-
- for(unsigned i = 0; i < p.size(); i++) {
- result += p[i]*xx;
- xx *= x;
- }
- return result;
-}
-
-std::vector<std::complex<double> > DK(Poly const & ply, const double tol) {
- std::vector<std::complex<double> > roots;
- const int N = ply.degree();
-
- std::complex<double> b(0.4, 0.9);
- std::complex<double> 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<double> denom = 1;
- std::complex<double> 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<double> dr = evalu(ply, R)/denom;
- error += norm(dr);
- roots[r_i] = R - dr;
- }
- /*std::copy(roots.begin(), roots.end(), std::ostream_iterator<std::complex<double> >(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"\r
+#include <iterator>\r
+\r
+/*** implementation of the Durand-Kerner method. seems buggy*/\r
+\r
+std::complex<double> evalu(Poly const & p, std::complex<double> x) {\r
+ std::complex<double> result = 0;\r
+ std::complex<double> xx = 1;\r
+ \r
+ for(unsigned i = 0; i < p.size(); i++) {\r
+ result += p[i]*xx;\r
+ xx *= x;\r
+ }\r
+ return result;\r
+}\r
+\r
+std::vector<std::complex<double> > DK(Poly const & ply, const double tol) {\r
+ std::vector<std::complex<double> > roots;\r
+ const int N = ply.degree();\r
+ \r
+ std::complex<double> b(0.4, 0.9);\r
+ std::complex<double> p = 1;\r
+ for(int i = 0; i < N; i++) {\r
+ roots.push_back(p);\r
+ p *= b;\r
+ }\r
+ assert(roots.size() == ply.degree());\r
+ \r
+ double error = 0;\r
+ int i;\r
+ for( i = 0; i < 30; i++) {\r
+ error = 0;\r
+ for(int r_i = 0; r_i < N; r_i++) {\r
+ std::complex<double> denom = 1;\r
+ std::complex<double> R = roots[r_i];\r
+ for(int d_i = 0; d_i < N; d_i++) {\r
+ if(r_i != d_i)\r
+ denom *= R-roots[d_i];\r
+ }\r
+ assert(norm(denom) != 0);\r
+ std::complex<double> dr = evalu(ply, R)/denom;\r
+ error += norm(dr);\r
+ roots[r_i] = R - dr;\r
+ }\r
+ /*std::copy(roots.begin(), roots.end(), std::ostream_iterator<std::complex<double> >(std::cout, ",\t"));\r
+ std::cout << std::endl;*/\r
+ if(error < tol)\r
+ break;\r
+ }\r
+ //std::cout << error << ", " << i<< std::endl;\r
+ return roots;\r
+}\r
+\r
+\r
+/*\r
+ Local Variables:\r
+ mode:c++\r
+ c-file-style:"stroustrup"\r
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))\r
+ indent-tabs-mode:nil\r
+ fill-column:99\r
+ End:\r
+*/\r
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :\r
index 83f049286df25b78167ff6bae9cd6e8d717d9d7b..fb3c2075b213da13d21d01a9793d675d3e9934f5 100644 (file)
-#include "poly-laguerre-solve.h"
-#include <iterator>
-
-typedef std::complex<double> 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<cdouble >
-laguerre(Poly p, const double tol) {
- std::vector<cdouble > 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<double>
-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"\r
+#include <iterator>\r
+\r
+typedef std::complex<double> cdouble;\r
+\r
+cdouble laguerre_internal_complex(Poly const & p,\r
+ double x0,\r
+ double tol,\r
+ bool & quad_root) {\r
+ cdouble a = 2*tol;\r
+ cdouble xk = x0;\r
+ double n = p.degree();\r
+ quad_root = false;\r
+ const unsigned shuffle_rate = 10;\r
+ static double shuffle[] = {0, 0.5, 0.25, 0.75, 0.125, 0.375, 0.625, 0.875, 1.0};\r
+ unsigned shuffle_counter = 0;\r
+ while(std::norm(a) > (tol*tol)) {\r
+ //std::cout << "xk = " << xk << std::endl;\r
+ cdouble b = p.back();\r
+ cdouble d = 0, f = 0;\r
+ double err = abs(b);\r
+ double abx = abs(xk);\r
+ for(int j = p.size()-2; j >= 0; j--) {\r
+ f = xk*f + d;\r
+ d = xk*d + b;\r
+ b = xk*b + p[j];\r
+ err = abs(b) + abx*err;\r
+ }\r
+ \r
+ err *= 1e-7; // magic epsilon for convergence, should be computed from tol\r
+ \r
+ cdouble px = b;\r
+ if(abs(b) < err)\r
+ return xk;\r
+ //if(std::norm(px) < tol*tol)\r
+ // return xk;\r
+ cdouble G = d / px;\r
+ cdouble H = G*G - f / px;\r
+ \r
+ //std::cout << "G = " << G << "H = " << H;\r
+ cdouble radicand = (n - 1)*(n*H-G*G);\r
+ //assert(radicand.real() > 0);\r
+ if(radicand.real() < 0)\r
+ quad_root = true;\r
+ //std::cout << "radicand = " << radicand << std::endl;\r
+ if(G.real() < 0) // here we try to maximise the denominator avoiding cancellation\r
+ a = - sqrt(radicand);\r
+ else\r
+ a = sqrt(radicand);\r
+ //std::cout << "a = " << a << std::endl;\r
+ a = n / (a + G);\r
+ //std::cout << "a = " << a << std::endl;\r
+ if(shuffle_counter % shuffle_rate == 0)\r
+ ;//a *= shuffle[shuffle_counter / shuffle_rate];\r
+ xk -= a;\r
+ shuffle_counter++;\r
+ if(shuffle_counter >= 90)\r
+ break;\r
+ }\r
+ //std::cout << "xk = " << xk << std::endl;\r
+ return xk;\r
+}\r
+\r
+double laguerre_internal(Poly const & p,\r
+ Poly const & pp,\r
+ Poly const & ppp,\r
+ double x0,\r
+ double tol,\r
+ bool & quad_root) {\r
+ double a = 2*tol;\r
+ double xk = x0;\r
+ double n = p.degree();\r
+ quad_root = false;\r
+ while(a*a > (tol*tol)) {\r
+ //std::cout << "xk = " << xk << std::endl;\r
+ double px = p(xk);\r
+ if(px*px < tol*tol)\r
+ return xk;\r
+ double G = pp(xk) / px;\r
+ double H = G*G - ppp(xk) / px;\r
+ \r
+ //std::cout << "G = " << G << "H = " << H;\r
+ double radicand = (n - 1)*(n*H-G*G);\r
+ assert(radicand > 0);\r
+ //std::cout << "radicand = " << radicand << std::endl;\r
+ if(G < 0) // here we try to maximise the denominator avoiding cancellation\r
+ a = - sqrt(radicand);\r
+ else\r
+ a = sqrt(radicand);\r
+ //std::cout << "a = " << a << std::endl;\r
+ a = n / (a + G);\r
+ //std::cout << "a = " << a << std::endl;\r
+ xk -= a;\r
+ }\r
+ //std::cout << "xk = " << xk << std::endl;\r
+ return xk;\r
+}\r
+\r
+\r
+std::vector<cdouble > \r
+laguerre(Poly p, const double tol) {\r
+ std::vector<cdouble > solutions;\r
+ //std::cout << "p = " << p << " = ";\r
+ while(p.size() > 1)\r
+ {\r
+ double x0 = 0;\r
+ bool quad_root = false;\r
+ cdouble sol = laguerre_internal_complex(p, x0, tol, quad_root);\r
+ //if(abs(sol) > 1) break;\r
+ Poly dvs;\r
+ if(quad_root) {\r
+ dvs.push_back((sol*conj(sol)).real());\r
+ dvs.push_back(-(sol + conj(sol)).real());\r
+ dvs.push_back(1.0);\r
+ //std::cout << "(" << dvs << ")";\r
+ //solutions.push_back(sol);\r
+ //solutions.push_back(conj(sol));\r
+ } else {\r
+ //std::cout << sol << std::endl;\r
+ dvs.push_back(-sol.real());\r
+ dvs.push_back(1.0);\r
+ solutions.push_back(sol);\r
+ //std::cout << "(" << dvs << ")";\r
+ }\r
+ Poly r;\r
+ p = divide(p, dvs, r);\r
+ //std::cout << r << std::endl;\r
+ }\r
+ return solutions;\r
+}\r
+\r
+std::vector<double> \r
+laguerre_real_interval(Poly const & ply, \r
+ const double lo, const double hi,\r
+ const double tol) {\r
+}\r
+\r
+/*\r
+ Local Variables:\r
+ mode:c++\r
+ c-file-style:"stroustrup"\r
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))\r
+ indent-tabs-mode:nil\r
+ fill-column:99\r
+ End:\r
+*/\r
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :\r
diff --git a/src/2geom/poly.cpp b/src/2geom/poly.cpp
index 5c3a300029b1b4ab3e02b8004ececc6ddc24c33c..17f28649401d43bce4939de5b45321371666b333 100644 (file)
--- a/src/2geom/poly.cpp
+++ b/src/2geom/poly.cpp
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 d96802014fd1af43e53aee91f6fdb0417749767e..f0811baf6f7359a238cf81c337c9b9147fdefe8b 100644 (file)
--- a/src/2geom/rect.h
+++ b/src/2geom/rect.h
-//D2<Interval> specialization to Rect:\r
-\r
- /* Authors of original rect class:\r
- * Lauris Kaplinski <lauris@kaplinski.com>\r
- * Nathan Hurst <njh@mail.csse.monash.edu.au>\r
- * bulia byak <buliabyak@users.sf.net>\r
- * MenTaLguY <mental@rydia.net>\r
- */\r
+/*
+ * rect.h - D2<Interval> specialization to Rect
+ *
+ * Copyright 2007 Michael Sloan <mgsloan@gmail.com>
+ *
+ * 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 <lauris@kaplinski.com>
+ * Nathan Hurst <njh@mail.csse.monash.edu.au>
+ * bulia byak <buliabyak@users.sf.net>
+ * MenTaLguY <mental@rydia.net>
+ */
#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
-\r
-#include "matrix.h"\r
+
+#include "matrix.h"
#include <boost/optional/optional.hpp>
-\r
-namespace Geom {\r
-\r
-typedef D2<Interval> Rect;\r
+
+namespace Geom {
+
+typedef D2<Interval> Rect;
Rect unify(const Rect &, const Rect &);
-\r
-template<>\r
-class D2<Interval> {\r
- private:\r
- Interval f[2];\r
- D2<Interval>();// { f[X] = f[Y] = Interval(0, 0); }\r
-\r
- public:\r
- D2<Interval>(Interval const &a, Interval const &b) {\r
- f[X] = a;\r
- f[Y] = b;\r
- }\r
-\r
- D2<Interval>(Point const & a, Point const & b) {\r
- f[X] = Interval(a[X], b[X]);\r
- f[Y] = Interval(a[Y], b[Y]);\r
- }\r
-\r
- inline Interval& operator[](unsigned i) { return f[i]; }\r
- inline Interval const & operator[](unsigned i) const { return f[i]; }\r
-\r
- inline Point min() const { return Point(f[X].min(), f[Y].min()); }\r
- inline Point max() const { return Point(f[X].max(), f[Y].max()); }\r
-\r
- /** returns the four corners of the rectangle in order\r
- * (clockwise if +Y is up, anticlockwise if +Y is down) */\r
- Point corner(unsigned i) const {\r
- switch(i % 4) {\r
- case 0: return Point(f[X].min(), f[Y].min());\r
- case 1: return Point(f[X].max(), f[Y].min());\r
- case 2: return Point(f[X].max(), f[Y].max());\r case 3: return Point(f[X].min(), f[Y].max());\r
- }\r
- }\r
+template<>
+class D2<Interval> {
+ private:
+ Interval f[2];
+ D2<Interval>();// { f[X] = f[Y] = Interval(0, 0); }
+
+ public:
+ D2<Interval>(Interval const &a, Interval const &b) {
+ f[X] = a;
+ f[Y] = b;
+ }
+
+ D2<Interval>(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(); }
-\r
- /** returns a vector from min to max. */\r
- inline Point dimensions() const { return Point(f[X].extent(), f[Y].extent()); }\r
- inline Point midpoint() const { return Point(f[X].middle(), f[Y].middle()); }\r
-\r
- inline double area() const { return f[X].extent() * f[Y].extent(); }\r
- inline double maxExtent() const { return std::max(f[X].extent(), f[Y].extent()); }\r
-\r
- inline bool isEmpty() const { return f[X].isEmpty() && f[Y].isEmpty(); }\r
- inline bool intersects(Rect const &r) const { return f[X].intersects(r[X]) && f[Y].intersects(r[Y]); }\r
- inline bool contains(Rect const &r) const { return f[X].contains(r[X]) && f[Y].contains(r[Y]); }\r
- inline bool contains(Point const &p) const { return f[X].contains(p[X]) && f[Y].contains(p[Y]); }\r
-\r
- inline void expandTo(Point p) { f[X].extendTo(p[X]); f[Y].extendTo(p[Y]); }\r
- inline void unionWith(Rect const &b) { f[X].unionWith(b[X]); f[Y].unionWith(b[Y]); }\r
-\r
- inline void expandBy(double amnt) { f[X].expandBy(amnt); f[Y].expandBy(amnt); }\r
- inline void expandBy(Point const p) { f[X].expandBy(p[X]); f[Y].expandBy(p[Y]); }\r
-\r
+
+ /** 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. */\r
- inline Rect operator*(Matrix const m) const { return unify(Rect(corner(0) * m, corner(2) * m),
- Rect(corner(1) * m, corner(3) * m)); }\r
+ 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<Rect> intersect(const Rect & a, const Rect & b) {
+inline Rect union_list(std::vector<Rect> 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<Rect> intersect(Rect const & a, Rect const & b) {
boost::optional<Interval> x = intersect(a[X], b[X]);
boost::optional<Interval> y = intersect(a[Y], b[Y]);
return x && y ? boost::optional<Rect>(Rect(*x, *y)) : boost::optional<Rect>();
#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
--- /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
--- /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<Rect> 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<Rect> const &b) : boundary(p), box(b) { fill = path_direction(p); }
+ Region(Path const &p, boost::optional<Rect> 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<Rect>(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<Region> Regions;
+
+unsigned outer_index(Regions const &ps);
+
+//assumes they're already sanitized somewhat
+inline Regions regions_from_paths(std::vector<Path> const &ps) {
+ Regions res;
+ for(unsigned i = 0; i < ps.size(); i++)
+ res.push_back(Region(ps[i]));
+ return res;
+}
+
+inline std::vector<Path> paths_from_regions(Regions const &rs) {
+ std::vector<Path> 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
index 170323c6ca19214d15669308757275ea10f1cd81..fde90785705c7bc13dc3c9f679614643fd32d744 100644 (file)
@@ -350,7 +350,7 @@ unsigned Geom::centroid(Piecewise<D2<SBasis> > 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.
index d88c2683261cd39535f2cb1584f5da1e48233107..db4cf5e08dee315c62eee2d354dbfe0a6f09d5f9 100644 (file)
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);
}
index 5296410686d0a018d2f57cb93b9ba90c989d2e63..8f4e1612dacd02a6211e4d8f2cdae3443070375a 100644 (file)
--- a/src/2geom/sbasis-math.h
+++ b/src/2geom/sbasis-math.h
Piecewise<SBasis> cos(Piecewise<SBasis> const &f, double tol=1e-3, int order=3);
Piecewise<SBasis> sin( SBasis const &f, double tol=1e-3, int order=3);
Piecewise<SBasis> sin(Piecewise<SBasis> const &f, double tol=1e-3, int order=3);
+//-Log---------------------------------------------------------------
+Piecewise<SBasis> log( SBasis const &f, double tol=1e-3, int order=3);
+Piecewise<SBasis> log(Piecewise<SBasis>const &f, double tol=1e-3, int order=3);
//--1/x------------------------------------------------------------
//TODO: change this...
index e0fa828f965efc333ed4dce5d22b3ab6929a2299..ff797920f7d2b376600a41169e7a7c7b2ee68ce6 100644 (file)
}
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);
index ad2dfbda4f868caa2c7bb18afcf52178af5dcd04..c4ef3c16d921050b33ec83f7ec05289fce06e6af 100644 (file)
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) {
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 || t<t0 || t>t1) {
lo = std::min(a*(1-t0)+b*t0+lo*t0*(1-t0),a*(1-t1)+b*t1+lo*t1*(1-t1));
index 613c39288b817bf15f936389d4aecf3b5b8ad0a3..206f18931430abdda1860f918419585919e08be6 100644 (file)
//TODO: some of this logic should be lifted into svg-path
std::vector<Geom::Path>
path_from_piecewise(Geom::Piecewise<Geom::D2<Geom::SBasis> > 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 a27d4bfcee4512c1a58b3a4a4cffdd43d021dde7..b0df4435c3653bc96be35c6913f11595330ffd46 100644 (file)
--- a/src/2geom/sbasis.h
+++ b/src/2geom/sbasis.h
#include "linear.h"
#include "interval.h"
+#include "utils.h"
namespace Geom {
double operator()(double t) const {
return valueAt(t);
}
+
+ std::vector<double> 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
--- /dev/null
+++ b/src/2geom/shape.cpp
@@ -0,0 +1,442 @@
+#include "shape.h"
+#include "utils.h"
+#include "sweep.h"
+
+#include <iostream>
+#include <algorithm>
+
+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<typename T>
+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<std::vector<bool> > visited, unsigned &i, unsigned &j) {
+ for(i = 0, j = 0; i < visited.size(); i++) {
+ std::vector<bool>::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<std::vector<bool> > visited;
+ for(unsigned i = 0; i < crs.size(); i++)
+ visited.push_back(std::vector<bool>(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<double> region_sizes(Shape const &a) {
+ std::vector<double> 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<Path> const &ps, Point p) {
+ int ret;
+ for(unsigned i = 0; i < ps.size(); i++)
+ ret += winding(ps[i], p);
+ return ret;
+}
+
+std::vector<double> y_of_roots(std::vector<Path> const & ps, double x) {
+ std::vector<double> res;
+ for(unsigned i = 0; i < ps.size(); i++) {
+ std::vector<double> 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<Path> 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<double> 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<Path> 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<Edge> 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<Path> 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<crs.size(); i++) {
+ if(crs[i].empty()) {
+ es.push_back(Edge(i, 0, ps[i].size(), false, ps));
+ es.push_back(Edge(i, ps[i].size(), 0, true, ps));
+ }
+ }
+
+ Edges es2;
+ //filter
+ for(unsigned i = 0; i < es.size(); i++) {
+ if(true) //(evenodd && es[i].wind % 2 == 0) || (!evenodd && es[i].wind == 0))
+ es2.push_back(es[i]);
+ }
+ es = es2;
+
+ std::cout << es.size() << " edges\n";
+
+ Regions chunks;
+ for(unsigned i = 0; i < es.size(); i++) {
+ Edge cur = es[i];
+ if(cur.rev)
+ chunks.push_back(Region(ps[cur.ix].portion(cur.from, cur.to).reverse(), cur.wind % 2 != 0));
+ else
+ chunks.push_back(Region(ps[cur.ix].portion(cur.from, cur.to), cur.wind % 2 != 0));
+ }
+ return chunks;
+
+ //Regions chunks;
+ std::vector<bool> 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<unsigned> 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<Region> const *rs;
+ explicit ContainmentOrder(std::vector<Region> 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<Rect> pnt;
+ pnt.push_back(Rect(p, p));
+ std::vector<std::vector<unsigned> > 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
--- /dev/null
+++ b/src/2geom/shape.h
@@ -0,0 +1,90 @@
+#ifndef __2GEOM_SHAPE_H
+#define __2GEOM_SHAPE_H
+
+#include <vector>
+#include <set>
+
+#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<Path> const &ps, bool evenodd=true);
+
+}
+
+#endif
index 558c10c0fbb4b6fe3028400822429564949f45d0..4f37cd0492866c4c5063a0b3e705567631961401 100644 (file)
* 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<double> &solutions, /* RETURN candidate t-values */
unsigned depth, /* The depth of the recursion */
{
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;
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;
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<double> &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 c2e2902516a7c8642bbeb009e6d0aec9a56ea28b..e0d66fc2a04329b3f56774cd51d9df36a552aaed 100644 (file)
--- a/src/2geom/solver.h
+++ b/src/2geom/solver.h
namespace Geom{
+ class Point;
+
unsigned
crossing_count(Geom::Point const *V, /* Control pts of Bezier curve */
unsigned degree); /* Degree 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<double> & 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<double> & solutions, /* RETURN candidate t-values */
unsigned depth, /* The depth of the recursion */
index 0063fcfe964cf6cad871599884cc3088013726c6..0bd15e8b9e55f94b0f0ae37279f152a3c0c88fea 100644 (file)
-#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
*
};
-#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,
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,
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,
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,
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
};
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)
_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;
}
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);
}
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();
}
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();
}
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();
}
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();
}
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"
}
}
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 141ddbcf38b931988c51fc07a5077e03250c7e79..312db9d2320fbd635f063f4ffdea77899f838c8d 100644 (file)
--- a/src/2geom/svg-path.cpp
+++ b/src/2geom/svg-path.cpp
namespace Geom {
void output(Curve const &curve, SVGPathSink &sink) {
- std::vector<Point> pts = sbasis_to_bezier(curve.sbasis(), 2); //TODO: use something better!
+ std::vector<Point> 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 d09002220df0630151ba24138737a4d881953fa5..ebc23a9b5383a0067dc77f17c144a7ca9da9f104 100644 (file)
--- a/src/2geom/svg-path.h
+++ b/src/2geom/svg-path.h
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);
_path.appendNew<LineSegment>(p);
}
- void curveTo(Point c0, Point c1, Point p) {
- _path.appendNew<CubicBezier>(c0, c1, p);
- }
-
void quadTo(Point c, Point p) {
_path.appendNew<QuadraticBezier>(c, p);
}
+ void curveTo(Point c0, Point c1, Point p) {
+ _path.appendNew<CubicBezier>(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
--- /dev/null
+++ b/src/2geom/sweep.cpp
@@ -0,0 +1,104 @@
+#include "sweep.h"\r
+\r
+#include <algorithm>\r
+\r
+namespace Geom {\r
+\r
+std::vector<std::vector<unsigned> > sweep_bounds(std::vector<Rect> rs) {\r
+ std::vector<Event> events; events.reserve(rs.size()*2);\r
+ std::vector<std::vector<unsigned> > pairs(rs.size());\r
+ \r
+ for(unsigned i = 0; i < rs.size(); i++) {\r
+ events.push_back(Event(rs[i].left(), i, false));\r
+ events.push_back(Event(rs[i].right(), i, true));\r
+ }\r
+ std::sort(events.begin(), events.end());\r
+\r
+ std::vector<unsigned> open;\r
+ for(unsigned i = 0; i < events.size(); i++) {\r
+ unsigned ix = events[i].ix;\r
+ if(events[i].closing) {\r
+ std::vector<unsigned>::iterator iter = std::find(open.begin(), open.end(), ix);\r
+ //if(iter != open.end())\r
+ open.erase(iter);\r
+ } else {\r
+ for(unsigned j = 0; j < open.size(); j++) {\r
+ unsigned jx = open[j];\r
+ if(rs[jx][Y].intersects(rs[ix][Y])) {\r
+ pairs[jx].push_back(ix);\r
+ }\r
+ }\r
+ open.push_back(ix);\r
+ }\r
+ }\r
+ return pairs;\r
+}\r
+\r
+std::vector<std::vector<unsigned> > sweep_bounds(std::vector<Rect> a, std::vector<Rect> b) {\r
+ std::vector<std::vector<unsigned> > pairs(a.size());\r
+ if(a.empty() || b.empty()) return pairs;\r
+ std::vector<Event> events[2];\r
+ events[0].reserve(a.size()*2);\r
+ events[1].reserve(b.size()*2);\r
+ \r
+ for(unsigned n = 0; n < 2; n++) {\r
+ unsigned sz = n ? b.size() : a.size();\r
+ events[n].reserve(sz*2);\r
+ for(unsigned i = 0; i < sz; i++) {\r
+ events[n].push_back(Event(n ? b[i].left() : a[i].left(), i, false));\r
+ events[n].push_back(Event(n ? b[i].right() : a[i].right(), i, true));\r
+ }\r
+ std::sort(events[n].begin(), events[n].end());\r
+ }\r
+\r
+ std::vector<unsigned> open[2];\r
+ bool n = events[1].front() < events[0].front();\r
+ for(unsigned i[] = {0,0}; i[n] < events[n].size();) {\r
+ unsigned ix = events[n][i[n]].ix;\r
+ bool closing = events[n][i[n]].closing;\r
+ //std::cout << n << "[" << ix << "] - " << (closing ? "closer" : "opener") << "\n";\r
+ if(closing) {\r
+ open[n].erase(std::find(open[n].begin(), open[n].end(), ix));\r
+ } else {\r
+ if(n) {\r
+ //n = 1\r
+ //opening a B, add to all open a\r
+ for(unsigned j = 0; j < open[0].size(); j++) {\r
+ unsigned jx = open[0][j];\r
+ if(a[jx][Y].intersects(b[ix][Y])) {\r
+ pairs[jx].push_back(ix);\r
+ }\r
+ }\r
+ } else {\r
+ //n = 0\r
+ //opening an A, add all open b\r
+ for(unsigned j = 0; j < open[1].size(); j++) {\r
+ unsigned jx = open[1][j];\r
+ if(b[jx][Y].intersects(a[ix][Y])) {\r
+ pairs[ix].push_back(jx);\r
+ }\r
+ }\r
+ }\r
+ open[n].push_back(ix);\r
+ }\r
+ i[n]++;\r
+ n = (events[!n][i[!n]] < events[n][i[n]]) ? !n : n;\r
+ }\r
+ return pairs;\r
+}\r
+\r
+//Fake cull, until the switch to the real sweep is made.\r
+std::vector<std::vector<unsigned> > fake_cull(unsigned a, unsigned b) {\r
+ std::vector<std::vector<unsigned> > ret;\r
+ \r
+ std::vector<unsigned> all;\r
+ for(unsigned j = 0; j < b; j++)\r
+ all.push_back(j);\r
+ \r
+ for(unsigned i = 0; i < a; i++)\r
+ ret.push_back(all);\r
+ \r
+ return ret;\r
+}\r
+\r
+}\r
diff --git a/src/2geom/sweep.h b/src/2geom/sweep.h
--- /dev/null
+++ b/src/2geom/sweep.h
@@ -0,0 +1,29 @@
+#ifndef __2GEOM_SWEEP_H__
+#define __2GEOM_SWEEP_H__
+
+#include <vector>
+#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<std::vector<unsigned> > sweep_bounds(std::vector<Rect>);
+std::vector<std::vector<unsigned> > sweep_bounds(std::vector<Rect>, std::vector<Rect>);
+
+std::vector<std::vector<unsigned> > fake_cull(unsigned a, unsigned b);
+
+}
+
+#endif
diff --git a/src/2geom/transforms.h b/src/2geom/transforms.h
index 1cd1d3e5f605e842373355349ec250411123c4c4..94058bff00f0d067f4035a290bcb23e78185a941 100644 (file)
--- a/src/2geom/transforms.h
+++ b/src/2geom/transforms.h
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 a2f906ff494ac533ab756dcd4a1efe5147fec4de..ca880640d49d41e588891fe57a61066e8f579d81 100644 (file)
--- a/src/2geom/utils.h
+++ b/src/2geom/utils.h
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.
index f04c80f2b14c3802ea90b84c2fd2184fb5744aad..9e5b966eaffaf22a4074ca83f0dd6342d99aba94 100644 (file)
// }\r
else { \r
//this case handles sbasis as well as all other curve types\r
- Geom::Path sbasis_path = path_from_sbasis(c->sbasis(), 0.1);\r
+ Geom::Path sbasis_path = Geom::path_from_sbasis(c->toSBasis(), 0.1);\r
\r
//recurse to convert the new path resulting from the sbasis to svgd\r
for(Geom::Path::iterator iter = sbasis_path.begin(); iter != sbasis_path.end(); ++iter) {\r