diff options
100 files changed, 4734 insertions, 891 deletions
diff --git a/src/2geom/basic-intersection.cpp b/src/2geom/basic-intersection.cpp index 97c4c6e5c..5375e5b58 100644 --- a/src/2geom/basic-intersection.cpp +++ b/src/2geom/basic-intersection.cpp @@ -1,5 +1,5 @@ -#include "basic-intersection.h" -#include "exception.h" +#include <2geom/basic-intersection.h> +#include <2geom/exception.h> unsigned intersect_steps = 0; @@ -60,7 +60,7 @@ find_intersections( vector<Geom::Point> const & A, } std::vector<std::pair<double, double> > -find_self_intersections(OldBezier const &Sb) { +find_self_intersections(OldBezier const &/*Sb*/) { THROW_NOTIMPLEMENTED(); } diff --git a/src/2geom/basic-intersection.h b/src/2geom/basic-intersection.h index 76abcce2a..0090b0305 100644 --- a/src/2geom/basic-intersection.h +++ b/src/2geom/basic-intersection.h @@ -1,7 +1,7 @@ -#include "sbasis.h" -#include "bezier-to-sbasis.h" -#include "sbasis-to-bezier.h" -#include "d2.h" +#include <2geom/sbasis.h> +#include <2geom/bezier-to-sbasis.h> +#include <2geom/sbasis-to-bezier.h> +#include <2geom/d2.h> namespace Geom { diff --git a/src/2geom/bezier-curve.h b/src/2geom/bezier-curve.h index ceab3cc11..244328f8a 100644 --- a/src/2geom/bezier-curve.h +++ b/src/2geom/bezier-curve.h @@ -38,9 +38,9 @@ #define _2GEOM_BEZIER_CURVE_H_ -#include "curve.h" -#include "sbasis-curve.h" // for non-native winding method -#include "bezier.h" +#include <2geom/curve.h> +#include <2geom/sbasis-curve.h> // for non-native winding method +#include <2geom/bezier.h> #include <algorithm> diff --git a/src/2geom/bezier-to-sbasis.h b/src/2geom/bezier-to-sbasis.h index f88913584..8b421f2e7 100644 --- a/src/2geom/bezier-to-sbasis.h +++ b/src/2geom/bezier-to-sbasis.h @@ -31,10 +31,10 @@ #ifndef _BEZIER_TO_SBASIS #define _BEZIER_TO_SBASIS -#include "coord.h" +#include <2geom/coord.h> -#include "d2.h" -#include "point.h" +#include <2geom/d2.h> +#include <2geom/point.h> namespace Geom{ diff --git a/src/2geom/bezier-utils.cpp b/src/2geom/bezier-utils.cpp index 59aac8951..83767af98 100644 --- a/src/2geom/bezier-utils.cpp +++ b/src/2geom/bezier-utils.cpp @@ -51,9 +51,9 @@ # include <ieefp.h> #endif -#include "bezier-utils.h" +#include <2geom/bezier-utils.h> -#include "isnan.h" +#include <2geom/isnan.h> #include <assert.h> namespace Geom{ diff --git a/src/2geom/bezier-utils.h b/src/2geom/bezier-utils.h index fba4333ff..68ae8c0e7 100644 --- a/src/2geom/bezier-utils.h +++ b/src/2geom/bezier-utils.h @@ -38,7 +38,7 @@ * */ -#include "point.h" +#include <2geom/point.h> namespace Geom{ diff --git a/src/2geom/bezier.h b/src/2geom/bezier.h index bc9d6032e..79b15cd03 100644 --- a/src/2geom/bezier.h +++ b/src/2geom/bezier.h @@ -33,12 +33,12 @@ #ifndef SEEN_BEZIER_H #define SEEN_BEZIER_H -#include "coord.h" +#include <2geom/coord.h> #include <valarray> -#include "isnan.h" -#include "bezier-to-sbasis.h" -#include "d2.h" -#include "solver.h" +#include <2geom/isnan.h> +#include <2geom/bezier-to-sbasis.h> +#include <2geom/d2.h> +#include <2geom/solver.h> #include <boost/optional/optional.hpp> namespace Geom { diff --git a/src/2geom/choose.h b/src/2geom/choose.h index f3a8b0f4f..bfd66e8a9 100644 --- a/src/2geom/choose.h +++ b/src/2geom/choose.h @@ -42,7 +42,7 @@ T choose(unsigned n, unsigned k) { static unsigned rows_done = 0; // indexing is (0,0,), (1,0), (1,1), (2, 0)... // to get (i, j) i*(i+1)/2 + j - if(k < 0 || k > n) return 0; + if(/*k < 0 ||*/ k > n) return 0; if(rows_done <= n) {// we haven't got there yet if(rows_done == 0) { pascals_triangle.push_back(1); diff --git a/src/2geom/circle-circle.cpp b/src/2geom/circle-circle.cpp index 024864091..25385180b 100644 --- a/src/2geom/circle-circle.cpp +++ b/src/2geom/circle-circle.cpp @@ -42,7 +42,7 @@ */ #include <stdio.h> #include <math.h> -#include "point.h" +#include <2geom/point.h> namespace Geom{ diff --git a/src/2geom/concepts.h b/src/2geom/concepts.h index 3b6fd3577..6f10c8bb0 100644 --- a/src/2geom/concepts.h +++ b/src/2geom/concepts.h @@ -31,9 +31,9 @@ #ifndef SEEN_CONCEPTS_H #define SEEN_CONCEPTS_H -#include "sbasis.h" -#include "interval.h" -#include "point.h" +#include <2geom/sbasis.h> +#include <2geom/interval.h> +#include <2geom/point.h> #include <vector> #include <boost/concept_check.hpp> diff --git a/src/2geom/conjugate_gradient.cpp b/src/2geom/conjugate_gradient.cpp index b98bb314c..f5a0f9cd8 100644 --- a/src/2geom/conjugate_gradient.cpp +++ b/src/2geom/conjugate_gradient.cpp @@ -32,7 +32,7 @@ #include <stdlib.h> #include <valarray> #include <cassert> -#include "conjugate_gradient.h" +#include <2geom/conjugate_gradient.h> /* lifted wholely from wikipedia. */ @@ -55,9 +55,12 @@ matrix_times_vector(valarray<double> const &matrix, /* m * n */ } } +/** +// only used in commented code below static double Linfty(valarray<double> const &vec) { return std::max(vec.max(), -vec.min()); } +**/ double inner(valarray<double> const &x, @@ -96,7 +99,7 @@ conjugate_gradient(valarray<double> const &A, valarray<double> &x, valarray<double> const &b, unsigned n, double tol, - unsigned max_iterations, bool ortho1) { + unsigned max_iterations, bool /*ortho1*/) { valarray<double> Ap(n), p(n), r(n); matrix_times_vector(A,x,Ap); r=b-Ap; diff --git a/src/2geom/convex-cover.cpp b/src/2geom/convex-cover.cpp index 50c43e6d1..7127a7c09 100644 --- a/src/2geom/convex-cover.cpp +++ b/src/2geom/convex-cover.cpp @@ -29,7 +29,7 @@ * */ -#include "convex-cover.h" +#include <2geom/convex-cover.h> #include <algorithm> #include <map> /** Todo: @@ -264,7 +264,7 @@ proposed algorithm: We must be very careful about rounding here. */ bool ConvexHull::no_colinear_points() const { - + // XXX: implement me! } bool @@ -292,8 +292,8 @@ bridges(ConvexHull a, ConvexHull b) { map<int, int> abridges; map<int, int> bbridges; - for(int ia = 0; ia < a.boundary.size(); ia++) { - for(int ib = 0; ib < b.boundary.size(); ib++) { + for(unsigned ia = 0; ia < a.boundary.size(); ia++) { + for(unsigned ib = 0; ib < b.boundary.size(); ib++) { Point d = b[ib] - a[ia]; Geom::Coord e = cross(d, a[ia - 1] - a[ia]), f = cross(d, a[ia + 1] - a[ia]); Geom::Coord g = cross(d, b[ib - 1] - a[ia]), h = cross(d, b[ib + 1] - a[ia]); @@ -336,8 +336,8 @@ unsigned find_bottom_right(ConvexHull const &a) { ConvexHull sweepline_intersection(ConvexHull const &a, ConvexHull const &b) { ConvexHull ret; - int al = 0; - int bl = 0; + unsigned al = 0; + unsigned bl = 0; while(al+1 < a.boundary.size() && (a.boundary[al+1][Y] > b.boundary[bl][Y])) { @@ -348,8 +348,8 @@ ConvexHull sweepline_intersection(ConvexHull const &a, ConvexHull const &b) { bl++; } // al and bl now point to the top of the first pair of edges that overlap in y value - double sweep_y = std::min(a.boundary[al][Y], - b.boundary[bl][Y]); + //double sweep_y = std::min(a.boundary[al][Y], + // b.boundary[bl][Y]); return ret; } @@ -358,7 +358,7 @@ ConvexHull sweepline_intersection(ConvexHull const &a, ConvexHull const &b) { * (Proof: take any two points both in a and in b. Any point between them is in a by convexity, * and in b by convexity, thus in both. Need to prove still finite bounds.) */ -ConvexHull intersection(ConvexHull a, ConvexHull b) { +ConvexHull intersection(ConvexHull /*a*/, ConvexHull /*b*/) { ConvexHull ret; /* int ai = 0, bi = 0; @@ -384,13 +384,13 @@ ConvexHull merge(ConvexHull a, ConvexHull b) { ab[-1] = 0; bb[-1] = 0; - int i = -1; + int i = -1; // XXX: i is int but refers to vector indices if(a.boundary[0][1] > b.boundary[0][1]) goto start_b; while(true) { for(; ab.count(i) == 0; i++) { ret.boundary.push_back(a[i]); - if(i >= a.boundary.size()) return ret; + if(i >= (int)a.boundary.size()) return ret; } if(ab[i] == 0 && i != -1) break; i = ab[i]; @@ -398,7 +398,7 @@ ConvexHull merge(ConvexHull a, ConvexHull b) { for(; bb.count(i) == 0; i++) { ret.boundary.push_back(b[i]); - if(i >= b.boundary.size()) return ret; + if(i >= (int)b.boundary.size()) return ret; } if(bb[i] == 0 && i != -1) break; i = bb[i]; diff --git a/src/2geom/convex-cover.h b/src/2geom/convex-cover.h index d99e07b95..1fe86e40d 100644 --- a/src/2geom/convex-cover.h +++ b/src/2geom/convex-cover.h @@ -36,7 +36,7 @@ * convex hull class is included here (the convex-hull header is wrong) */ -#include "point.h" +#include <2geom/point.h> #include <vector> namespace Geom{ diff --git a/src/2geom/crossing.cpp b/src/2geom/crossing.cpp index 880b99e1a..a4d4f6ab1 100644 --- a/src/2geom/crossing.cpp +++ b/src/2geom/crossing.cpp @@ -1,5 +1,5 @@ -#include "crossing.h" -#include "path.h" +#include <2geom/crossing.h> +#include <2geom/path.h> namespace Geom { @@ -169,7 +169,7 @@ CrossingSet reverse_tb(CrossingSet const &cr, unsigned split, std::vector<double return ret; } -void clean(Crossings &cr_a, Crossings &cr_b) { +void clean(Crossings &/*cr_a*/, Crossings &/*cr_b*/) { /* if(cr_a.empty()) return; //Remove anything with dupes diff --git a/src/2geom/crossing.h b/src/2geom/crossing.h index 8e3a45142..546e33ebd 100644 --- a/src/2geom/crossing.h +++ b/src/2geom/crossing.h @@ -2,8 +2,8 @@ #define __GEOM_CROSSING_H #include <vector> -#include "rect.h" -#include "sweep.h" +#include <2geom/rect.h> +#include <2geom/sweep.h> namespace Geom { diff --git a/src/2geom/curve-helpers.cpp b/src/2geom/curve-helpers.cpp index c0e46bc06..c767af54f 100644 --- a/src/2geom/curve-helpers.cpp +++ b/src/2geom/curve-helpers.cpp @@ -31,8 +31,8 @@ */ -#include "curve.h" -#include "ord.h" +#include <2geom/curve.h> +#include <2geom/ord.h> namespace Geom diff --git a/src/2geom/curve.h b/src/2geom/curve.h index 22d2ec556..8eaf73fcf 100644 --- a/src/2geom/curve.h +++ b/src/2geom/curve.h @@ -38,14 +38,14 @@ #define _2GEOM_CURVE_H_ -#include "coord.h" -#include "point.h" -#include "interval.h" -#include "nearest-point.h" -#include "sbasis.h" -#include "d2.h" -#include "matrix.h" -#include "exception.h" +#include <2geom/coord.h> +#include <2geom/point.h> +#include <2geom/interval.h> +#include <2geom/nearest-point.h> +#include <2geom/sbasis.h> +#include <2geom/d2.h> +#include <2geom/matrix.h> +#include <2geom/exception.h> #include <vector> @@ -101,6 +101,14 @@ public: return all_nearest_points(p, toSBasis(), from, to); } + /* + Path operator*=(Matrix) + This is not possible, because: + A Curve can be many things, for example a HLineSegment. + Such a segment cannot be transformed and stay a HLineSegment in general (take for example rotations). + This means that these curves become a different type of curve, hence one should use "transformed(Matrix). + */ + virtual Curve *transformed(Matrix const &m) const = 0; virtual Point pointAt(Coord t) const { return pointAndDerivatives(t, 0).front(); } @@ -129,6 +137,7 @@ public: }; virtual D2<SBasis> toSBasis() const = 0; + virtual bool operator==(Curve const &c) const { return this == &c;} }; inline diff --git a/src/2geom/curves.h b/src/2geom/curves.h index a4065930a..3bf7d9b59 100644 --- a/src/2geom/curves.h +++ b/src/2geom/curves.h @@ -4,7 +4,7 @@ * Authors: * MenTaLguY <mental@rydia.net> * Marco Cecchetti <mrcekets at gmail.com> - * + * * Copyright 2007-2008 authors * * This library is free software; you can redistribute it and/or @@ -38,11 +38,12 @@ #define _2GEOM_CURVES_H_ -#include "curve.h" -#include "sbasis-curve.h" -#include "bezier-curve.h" -#include "hvlinesegment.h" -#include "elliptical-arc.h" +#include <2geom/curve.h> +#include <2geom/sbasis-curve.h> +#include <2geom/bezier-curve.h> +#include <2geom/hvlinesegment.h> +#include <2geom/elliptical-arc.h> +#include <2geom/svg-elliptical-arc.h> #endif // _2GEOM_CURVES_H_ diff --git a/src/2geom/d2-sbasis.cpp b/src/2geom/d2-sbasis.cpp index 9b6ca269c..01f83bf5e 100644 --- a/src/2geom/d2-sbasis.cpp +++ b/src/2geom/d2-sbasis.cpp @@ -1,4 +1,4 @@ -#include "d2.h" +#include <2geom/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. */ @@ -38,7 +38,7 @@ Piecewise<D2<SBasis> > sectionize(D2<Piecewise<SBasis> > const &a) { return ret; } -D2<Piecewise<SBasis> > make_cuts_independant(Piecewise<D2<SBasis> > const &a) { +D2<Piecewise<SBasis> > make_cuts_independent(Piecewise<D2<SBasis> > const &a) { D2<Piecewise<SBasis> > ret; for(unsigned d = 0; d < 2; d++) { for(unsigned i = 0; i < a.size(); i++) diff --git a/src/2geom/d2-sbasis.h b/src/2geom/d2-sbasis.h index e921896f5..91ceb31c9 100644 --- a/src/2geom/d2-sbasis.h +++ b/src/2geom/d2-sbasis.h @@ -4,10 +4,10 @@ #ifndef __2GEOM_SBASIS_CURVE_H #define __2GEOM_SBASIS_CURVE_H -#include "sbasis.h" -#include "sbasis-2d.h" -#include "piecewise.h" -#include "matrix.h" +#include <2geom/sbasis.h> +#include <2geom/sbasis-2d.h> +#include <2geom/piecewise.h> +#include <2geom/matrix.h> //TODO: implement intersect @@ -32,7 +32,7 @@ double tail_error(D2<SBasis> const & a, unsigned tail); //Piecewise<D2<SBasis> > specific decls: Piecewise<D2<SBasis> > sectionize(D2<Piecewise<SBasis> > const &a); -D2<Piecewise<SBasis> > make_cuts_independant(Piecewise<D2<SBasis> > const &a); +D2<Piecewise<SBasis> > make_cuts_independent(Piecewise<D2<SBasis> > const &a); Piecewise<D2<SBasis> > rot90(Piecewise<D2<SBasis> > const &a); Piecewise<SBasis> dot(Piecewise<D2<SBasis> > const &a, Piecewise<D2<SBasis> > const &b); Piecewise<SBasis> cross(Piecewise<D2<SBasis> > const &a, Piecewise<D2<SBasis> > const &b); diff --git a/src/2geom/d2.h b/src/2geom/d2.h index db8cf68c4..731a084a1 100644 --- a/src/2geom/d2.h +++ b/src/2geom/d2.h @@ -31,12 +31,12 @@ #ifndef _2GEOM_D2 //If this is change, change the guard in rect.h as well. #define _2GEOM_D2 -#include "point.h" -#include "interval.h" -#include "matrix.h" +#include <2geom/point.h> +#include <2geom/interval.h> +#include <2geom/matrix.h> #include <boost/concept_check.hpp> -#include "concepts.h" +#include <2geom/concepts.h> namespace Geom{ @@ -373,8 +373,8 @@ D2<T> integral(D2<T> const & a) { } //end namespace Geom -#include "rect.h" -#include "d2-sbasis.h" +#include <2geom/rect.h> +#include <2geom/d2-sbasis.h> namespace Geom{ diff --git a/src/2geom/ellipse.cpp b/src/2geom/ellipse.cpp new file mode 100644 index 000000000..c9c5b9ec4 --- /dev/null +++ b/src/2geom/ellipse.cpp @@ -0,0 +1,206 @@ +/*
+ * Ellipse Curve
+ *
+ * Authors:
+ * Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2008 authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+
+#include <2geom/ellipse.h>
+#include <2geom/svg-elliptical-arc.h>
+#include <2geom/numeric/fitting-tool.h>
+#include <2geom/numeric/fitting-model.h>
+
+
+namespace Geom
+{
+
+void Ellipse::set(double A, double B, double C, double D, double E, double F)
+{
+ double den = 4*A*C - B*B;
+ if ( den == 0 )
+ {
+ THROW_LOGICALERROR("den == 0, while computing ellipse centre");
+ }
+ m_centre[X] = (B*E - 2*C*D) / den;
+ m_centre[Y] = (B*D - 2*A*E) / den;
+
+ // evaluate the a coefficient of the ellipse equation in normal form
+ // E(x,y) = a*(x-cx)^2 + b*(x-cx)*(y-cy) + c*(y-cy)^2 = 1
+ // where b = a*B , c = a*C, (cx,cy) == centre
+ double num = A * sqr(m_centre[X])
+ + B * m_centre[X] * m_centre[Y]
+ + C * sqr(m_centre[Y])
+ - A * F;
+
+
+ //evaluate ellipse rotation angle
+ double rot = std::atan2( -B, -(A - C) )/2;
+// std::cerr << "rot = " << rot << std::endl;
+ bool swap_axes = false;
+ if ( are_near(rot, 0) ) rot = 0;
+ if ( are_near(rot, M_PI/2) || rot < 0 )
+ {
+ swap_axes = true;
+ }
+
+ // evaluate the length of the ellipse rays
+ double cosrot = std::cos(rot);
+ double sinrot = std::sin(rot);
+ double cos2 = cosrot * cosrot;
+ double sin2 = sinrot * sinrot;
+ double cossin = cosrot * sinrot;
+
+ den = A * cos2 + B * cossin + C * sin2;
+ if ( den == 0 )
+ {
+ THROW_LOGICALERROR("den == 0, while computing 'rx' coefficient");
+ }
+ double rx2 = num/den;
+ if ( rx2 < 0 )
+ {
+ THROW_LOGICALERROR("rx2 < 0, while computing 'rx' coefficient");
+ }
+ double rx = std::sqrt(rx2);
+
+ den = C * cos2 - B * cossin + A * sin2;
+ if ( den == 0 )
+ {
+ THROW_LOGICALERROR("den == 0, while computing 'ry' coefficient");
+ }
+ double ry2 = num/den;
+ if ( ry2 < 0 )
+ {
+ THROW_LOGICALERROR("ry2 < 0, while computing 'rx' coefficient");
+ }
+ double ry = std::sqrt(ry2);
+
+ // the solution is not unique so we choose always the ellipse
+ // with a rotation angle between 0 and PI/2
+ if ( swap_axes ) std::swap(rx, ry);
+ if ( are_near(rot, M_PI/2)
+ || are_near(rot, -M_PI/2)
+ || are_near(rx, ry) )
+ {
+ rot = 0;
+ }
+ else if ( rot < 0 )
+ {
+ rot += M_PI/2;
+ }
+
+ m_ray[X] = rx;
+ m_ray[Y] = ry;
+ m_angle = rot;
+}
+
+
+void Ellipse::set(std::vector<Point> const& points)
+{
+ size_t sz = points.size();
+ if (sz < 5)
+ {
+ THROW_RANGEERROR("fitting error: too few points passed");
+ }
+ NL::LFMEllipse model;
+ NL::least_squeares_fitter<NL::LFMEllipse> fitter(model, sz);
+
+ for (size_t i = 0; i < sz; ++i)
+ {
+ fitter.append(points[i]);
+ }
+ fitter.update();
+
+ NL::Vector z(sz, 0.0);
+ model.instance(*this, fitter.result(z));
+}
+
+
+SVGEllipticalArc
+Ellipse::arc(Point const& initial, Point const& inner, Point const& final,
+ bool _svg_compliant)
+{
+ Point sp_cp = initial - center();
+ Point ep_cp = final - center();
+ Point ip_cp = inner - center();
+
+ double angle1 = angle_between(sp_cp, ep_cp);
+ double angle2 = angle_between(sp_cp, ip_cp);
+ double angle3 = angle_between(ip_cp, ep_cp);
+
+ bool large_arc_flag = true;
+ bool sweep_flag = true;
+
+ if ( angle1 > 0 )
+ {
+ if ( angle2 > 0 && angle3 > 0 )
+ {
+ large_arc_flag = false;
+ sweep_flag = true;
+ }
+ else
+ {
+ large_arc_flag = true;
+ sweep_flag = false;
+ }
+ }
+ else
+ {
+ if ( angle2 < 0 && angle3 < 0 )
+ {
+ large_arc_flag = false;
+ sweep_flag = false;
+ }
+ else
+ {
+ large_arc_flag = true;
+ sweep_flag = true;
+ }
+ }
+
+ SVGEllipticalArc ea( initial, ray(X), ray(Y), rot_angle(),
+ large_arc_flag, sweep_flag, final, _svg_compliant);
+ return ea;
+}
+
+
+} // end namespace Geom
+
+
+/*
+ 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/ellipse.h b/src/2geom/ellipse.h new file mode 100644 index 000000000..af8b01e78 --- /dev/null +++ b/src/2geom/ellipse.h @@ -0,0 +1,133 @@ +/*
+ * Ellipse Curve
+ *
+ * Authors:
+ * Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2008 authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+
+#ifndef _2GEOM_ELLIPSE_H_
+#define _2GEOM_ELLIPSE_H_
+
+
+#include <2geom/point.h>
+#include <2geom/exception.h>
+
+
+namespace Geom
+{
+
+class SVGEllipticalArc;
+
+class Ellipse
+{
+ public:
+ Ellipse()
+ {}
+
+ Ellipse(double cx, double cy, double rx, double ry, double a)
+ : m_centre(cx, cy), m_ray(rx, ry), m_angle(a)
+ {
+ }
+
+ // build an ellipse by its implicit equation:
+ // Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0
+ Ellipse(double A, double B, double C, double D, double E, double F)
+ {
+ set(A, B, C, D, E, F);
+ }
+
+ Ellipse(std::vector<Point> const& points)
+ {
+ set(points);
+ }
+
+ void set(double cx, double cy, double rx, double ry, double a)
+ {
+ m_centre[X] = cx;
+ m_centre[Y] = cy;
+ m_ray[X] = rx;
+ m_ray[Y] = ry;
+ m_angle = a;
+ }
+
+ // build an ellipse by its implicit equation:
+ // Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0
+ void set(double A, double B, double C, double D, double E, double F);
+
+ // biuld up the best fitting ellipse wrt the passed points
+ // prerequisite: at least 5 points must be passed
+ void set(std::vector<Point> const& points);
+
+ SVGEllipticalArc
+ arc(Point const& initial, Point const& inner, Point const& final,
+ bool _svg_compliant = true);
+
+ Point center() const
+ {
+ return m_centre;
+ }
+
+ Coord center(Dim2 d) const
+ {
+ return m_centre[d];
+ }
+
+ Coord ray(Dim2 d) const
+ {
+ return m_ray[d];
+ }
+
+ Coord rot_angle() const
+ {
+ return m_angle;
+ }
+
+ private:
+ Point m_centre, m_ray;
+ double m_angle;
+};
+
+
+} // end namespace Geom
+
+
+
+#endif // _2GEOM_ELLIPSE_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/elliptical-arc.cpp b/src/2geom/elliptical-arc.cpp index b30dfde5d..552cdf5b3 100644 --- a/src/2geom/elliptical-arc.cpp +++ b/src/2geom/elliptical-arc.cpp @@ -28,9 +28,9 @@ */ -#include "elliptical-arc.h" -#include "bezier-curve.h" -#include "poly.h" +#include <2geom/elliptical-arc.h> +#include <2geom/bezier-curve.h> +#include <2geom/poly.h> #include <cfloat> #include <limits> diff --git a/src/2geom/elliptical-arc.h b/src/2geom/elliptical-arc.h index 3a3b90670..46b94c0df 100644 --- a/src/2geom/elliptical-arc.h +++ b/src/2geom/elliptical-arc.h @@ -38,10 +38,10 @@ #define _2GEOM_ELLIPTICAL_ARC_H_ -#include "curve.h" -#include "angle.h" -#include "utils.h" -#include "sbasis-curve.h" // for non-native methods +#include <2geom/curve.h> +#include <2geom/angle.h> +#include <2geom/utils.h> +#include <2geom/sbasis-curve.h> // for non-native methods #include <algorithm> diff --git a/src/2geom/geom.cpp b/src/2geom/geom.cpp index 8845cefa8..64bf21ee6 100644 --- a/src/2geom/geom.cpp +++ b/src/2geom/geom.cpp @@ -6,8 +6,8 @@ #ifdef HAVE_CONFIG_H # include <config.h> #endif -#include "geom.h" -#include "point.h" +#include <2geom/geom.h> +#include <2geom/point.h> namespace Geom { diff --git a/src/2geom/geom.h b/src/2geom/geom.h index 49156482c..207d609b0 100644 --- a/src/2geom/geom.h +++ b/src/2geom/geom.h @@ -37,7 +37,7 @@ //TODO: move somewhere else #include <vector> -#include "point.h" +#include <2geom/point.h> namespace Geom { diff --git a/src/2geom/hvlinesegment.h b/src/2geom/hvlinesegment.h index 732c29938..a34c5a962 100644 --- a/src/2geom/hvlinesegment.h +++ b/src/2geom/hvlinesegment.h @@ -32,7 +32,7 @@ #define _2GEOM_HVLINESEGMENT_H_ -#include "bezier-curve.h" +#include <2geom/bezier-curve.h> namespace Geom diff --git a/src/2geom/interval.h b/src/2geom/interval.h index eb506dc1f..f042294ff 100644 --- a/src/2geom/interval.h +++ b/src/2geom/interval.h @@ -37,7 +37,7 @@ #define SEEN_INTERVAL_H #include <assert.h> -#include "coord.h" +#include <2geom/coord.h> #include <boost/optional/optional.hpp> diff --git a/src/2geom/isnan.h b/src/2geom/isnan.h index d95a45f10..6b94daa6e 100644 --- a/src/2geom/isnan.h +++ b/src/2geom/isnan.h @@ -34,7 +34,7 @@ # define IS_NAN(_a) (_isnan(_a)) /* Win32 definition */ #elif defined(isnan) || defined(__FreeBSD__) || defined(__osf__) # define IS_NAN(_a) (isnan(_a)) /* GNU definition */ -#elif defined (SOLARIS) +#elif defined (SOLARIS_2_8) && __GNUC__ == 3 && __GNUC_MINOR__ == 2 # define IS_NAN(_a) (isnan(_a)) /* GNU definition */ #else # define IS_NAN(_a) (std::isnan(_a)) @@ -55,7 +55,7 @@ # define IS_FINITE(_a) (isfinite(_a)) #elif defined(__osf__) # define IS_FINITE(_a) (finite(_a) && !IS_NAN(_a)) -#elif defined (SOLARIS) +#elif defined (SOLARIS_2_8) && __GNUC__ == 3 && __GNUC_MINOR__ == 2 #include <ieeefp.h> #define IS_FINITE(_a) (finite(_a) && !IS_NAN(_a)) #else diff --git a/src/2geom/linear.h b/src/2geom/linear.h index 271e87be4..bc7af564c 100644 --- a/src/2geom/linear.h +++ b/src/2geom/linear.h @@ -33,8 +33,8 @@ #ifndef SEEN_LINEAR_H #define SEEN_LINEAR_H -#include "interval.h" -#include "isnan.h" +#include <2geom/interval.h> +#include <2geom/isnan.h> namespace Geom{ diff --git a/src/2geom/matrix.cpp b/src/2geom/matrix.cpp index f90bb6d42..d25cc8f7e 100644 --- a/src/2geom/matrix.cpp +++ b/src/2geom/matrix.cpp @@ -12,9 +12,9 @@ * This code is in public domain */ -#include "utils.h" -#include "matrix.h" -#include "point.h" +#include <2geom/utils.h> +#include <2geom/matrix.h> +#include <2geom/point.h> namespace Geom { diff --git a/src/2geom/matrix.h b/src/2geom/matrix.h index ba4451265..dcf765c50 100644 --- a/src/2geom/matrix.h +++ b/src/2geom/matrix.h @@ -19,7 +19,7 @@ //#include <glib/gmessages.h> -#include "point.h" +#include <2geom/point.h> namespace Geom { diff --git a/src/2geom/nearest-point.cpp b/src/2geom/nearest-point.cpp index 074de1cb3..59d55ba32 100644 --- a/src/2geom/nearest-point.cpp +++ b/src/2geom/nearest-point.cpp @@ -32,7 +32,7 @@ */
-#include "nearest-point.h"
+#include <2geom/nearest-point.h>
namespace Geom
{
@@ -87,9 +87,9 @@ double nearest_point( Point const& p, std::vector<double>
all_nearest_points( Point const& p,
- D2<SBasis> const& c,
- D2<SBasis> const& dc,
- double from, double to )
+ D2<SBasis> const& c,
+ D2<SBasis> const& /*dc*/,
+ double from, double to )
{
std::swap(from, to);
if ( from > to ) std::swap(from, to);
diff --git a/src/2geom/nearest-point.h b/src/2geom/nearest-point.h index 73ac0c3ce..0b43ce67b 100644 --- a/src/2geom/nearest-point.h +++ b/src/2geom/nearest-point.h @@ -39,9 +39,9 @@ #include <vector>
-#include "d2.h"
-#include "piecewise.h"
-#include "exception.h"
+#include <2geom/d2.h>
+#include <2geom/piecewise.h>
+#include <2geom/exception.h>
diff --git a/src/2geom/numeric/fitting-model.h b/src/2geom/numeric/fitting-model.h new file mode 100644 index 000000000..145be40e4 --- /dev/null +++ b/src/2geom/numeric/fitting-model.h @@ -0,0 +1,423 @@ +/*
+ * Fitting Models for Geom Types
+ *
+ * Authors:
+ * Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2008 authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+
+#ifndef _NL_FITTING_MODEL_H_
+#define _NL_FITTING_MODEL_H_
+
+
+#include <2geom/d2.h>
+#include <2geom/sbasis.h>
+#include <2geom/bezier.h>
+#include <2geom/bezier-curve.h>
+#include <2geom/poly.h>
+#include <2geom/ellipse.h>
+#include <2geom/utils.h>
+
+
+namespace Geom { namespace NL {
+
+
+/*
+ * completely unknown models must inherit from this template class;
+ * example: the model a*x^2 + b*x + c = 0 to be solved wrt a, b, c;
+ * example: the model A(t) = known_sample_value_at(t) to be solved wrt
+ * the coefficients of the curve A(t) expressed in S-Basis form;
+ * parameter type: the type of x and t variable in the examples above;
+ * value type: the type of the known sample values (in the first example
+ * is constant )
+ * instance type: the type of the objects produced by using
+ * the fitting raw data solution
+ */
+template< typename ParameterType, typename ValueType, typename InstanceType >
+class LinearFittingModel
+{
+ public:
+ typedef ParameterType parameter_type;
+ typedef ValueType value_type;
+ typedef InstanceType instance_type;
+
+ static const bool WITH_FIXED_TERMS = false;
+
+ /*
+ * a LinearFittingModel must implement the following methods:
+ *
+ * void feed( VectorView & vector,
+ * parameter_type const& sample_parameter ) const;
+ *
+ * size_t size() const;
+ *
+ * void instance(instance_type &, raw_type const& raw_data) const;
+ *
+ */
+};
+
+
+/*
+ * partially known models must inherit from this template class
+ * example: the model a*x^2 + 2*x + c = 0 to be solved wrt a and c
+ */
+template< typename ParameterType, typename ValueType, typename InstanceType >
+class LinearFittingModelWithFixedTerms
+{
+ public:
+ typedef ParameterType parameter_type;
+ typedef ValueType value_type;
+ typedef InstanceType instance_type;
+
+ static const bool WITH_FIXED_TERMS = true;
+
+ /*
+ * a LinearFittingModelWithFixedTerms must implement the following methods:
+ *
+ * void feed( VectorView & vector,
+ * value_type & fixed_term,
+ * parameter_type const& sample_parameter ) const;
+ *
+ * size_t size() const;
+ *
+ * void instance(instance_type &, raw_type const& raw_data) const;
+ *
+ */
+
+
+};
+
+
+// incomplete model, it can be inherited to make up different kinds of
+// instance type; the raw data is a vector of coefficients of a polynomial
+// rapresented in standard power basis
+template< typename InstanceType >
+class LFMPowerBasis
+ : public LinearFittingModel<double, double, InstanceType>
+{
+ public:
+ LFMPowerBasis(size_t degree)
+ : m_size(degree + 1)
+ {
+ }
+
+ void feed( VectorView & coeff, double sample_parameter ) const
+ {
+ coeff[0] = 1;
+ double x_i = 1;
+ for (size_t i = 1; i < coeff.size(); ++i)
+ {
+ x_i *= sample_parameter;
+ coeff[i] = x_i;
+ }
+ }
+
+ size_t size() const
+ {
+ return m_size;
+ }
+
+ private:
+ size_t m_size;
+};
+
+
+// this model generates Geom::Poly objects
+class LFMPoly
+ : public LFMPowerBasis<Poly>
+{
+ public:
+ LFMPoly(size_t degree)
+ : LFMPowerBasis<Poly>(degree)
+ {
+ }
+
+ void instance(Poly & poly, ConstVectorView const& raw_data) const
+ {
+ poly.clear();
+ poly.resize(size());
+ for (size_t i = 0; i < raw_data.size(); ++i)
+ {
+ poly[i] = raw_data[i];
+ }
+ }
+};
+
+
+// incomplete model, it can be inherited to make up different kinds of
+// instance type; the raw data is a vector of coefficients of a polynomial
+// rapresented in standard power basis with leading term coefficient equal to 1
+template< typename InstanceType >
+class LFMNormalizedPowerBasis
+ : public LinearFittingModelWithFixedTerms<double, double, InstanceType>
+{
+ public:
+ LFMNormalizedPowerBasis(size_t _degree)
+ : m_model( _degree - 1)
+ {
+ assert(_degree > 0);
+ }
+
+
+ void feed( VectorView & coeff,
+ double & known_term,
+ double sample_parameter ) const
+ {
+ m_model.feed(coeff, sample_parameter);
+ known_term = coeff[m_model.size()-1] * sample_parameter;
+ }
+
+ size_t size() const
+ {
+ return m_model.size();
+ }
+
+ private:
+ LFMPowerBasis<InstanceType> m_model;
+};
+
+
+// incomplete model, it can be inherited to make up different kinds of
+// instance type; the raw data is a vector of coefficients of the equation
+// of an ellipse curve
+template< typename InstanceType >
+class LFMEllipseEquation
+ : public LinearFittingModelWithFixedTerms<Point, double, InstanceType>
+{
+ public:
+ void feed( VectorView & coeff, double & fixed_term, Point const& p ) const
+ {
+ coeff[0] = p[X] * p[Y];
+ coeff[1] = p[Y] * p[Y];
+ coeff[2] = p[X];
+ coeff[3] = p[Y];
+ coeff[4] = 1;
+ fixed_term = p[X] * p[X];
+ }
+
+ size_t size() const
+ {
+ return 5;
+ }
+};
+
+
+// this model generates Ellipse curves
+class LFMEllipse
+ : public LFMEllipseEquation<Ellipse>
+{
+ public:
+ void instance(Ellipse & e, ConstVectorView const& coeff) const
+ {
+ e.set(1, coeff[0], coeff[1], coeff[2], coeff[3], coeff[4]);
+ }
+};
+
+
+// this model generates SBasis objects
+class LFMSBasis
+ : public LinearFittingModel<double, double, SBasis>
+{
+ public:
+ LFMSBasis( size_t _order )
+ : m_size( 2*(_order+1) ),
+ m_order(_order)
+ {
+ }
+
+ void feed( VectorView & coeff, double t ) const
+ {
+ double u0 = 1-t;
+ double u1 = t;
+ double s = u0 * u1;
+ coeff[0] = u0;
+ coeff[1] = u1;
+ for (size_t i = 2; i < size(); i+=2)
+ {
+ u0 *= s;
+ u1 *= s;
+ coeff[i] = u0;
+ coeff[i+1] = u1;
+ }
+ }
+
+ size_t size() const
+ {
+ return m_size;
+ }
+
+ void instance(SBasis & sb, ConstVectorView const& raw_data) const
+ {
+ sb.clear();
+ sb.resize(m_order+1);
+ for (unsigned int i = 0, k = 0; i < raw_data.size(); i+=2, ++k)
+ {
+ sb[k][0] = raw_data[i];
+ sb[k][1] = raw_data[i+1];
+ }
+ }
+
+ private:
+ size_t m_size;
+ size_t m_order;
+};
+
+
+// this model generates D2<SBasis> objects
+class LFMD2SBasis
+ : public LinearFittingModel< double, Point, D2<SBasis> >
+{
+ public:
+ LFMD2SBasis( size_t _order )
+ : mosb(_order)
+ {
+ }
+
+ void feed( VectorView & coeff, double t ) const
+ {
+ mosb.feed(coeff, t);
+ }
+
+ size_t size() const
+ {
+ return mosb.size();
+ }
+
+ void instance(D2<SBasis> & d2sb, ConstMatrixView const& raw_data) const
+ {
+ mosb.instance(d2sb[X], raw_data.column_const_view(X));
+ mosb.instance(d2sb[Y], raw_data.column_const_view(Y));
+ }
+
+ private:
+ LFMSBasis mosb;
+};
+
+
+// this model generates Bezier objects
+class LFMBezier
+ : public LinearFittingModel<double, double, Bezier>
+{
+ public:
+ LFMBezier( size_t _order )
+ : m_size(_order + 1),
+ m_order(_order)
+ {
+ binomial_coefficients(m_bc, m_order);
+ }
+
+ void feed( VectorView & coeff, double t ) const
+ {
+ double s = 1;
+ for (size_t i = 0; i < size(); ++i)
+ {
+ coeff[i] = s * m_bc[i];
+ s *= t;
+ }
+ double u = 1-t;
+ s = 1;
+ for (size_t i = size()-1; i > 0; --i)
+ {
+ coeff[i] *= s;
+ s *= u;
+ }
+ coeff[0] *= s;
+ }
+
+ size_t size() const
+ {
+ return m_size;
+ }
+
+ void instance(Bezier & b, ConstVectorView const& raw_data) const
+ {
+ assert(b.size() == raw_data.size());
+ for (unsigned int i = 0; i < raw_data.size(); ++i)
+ {
+ b[i] = raw_data[i];
+ }
+ }
+
+ private:
+ size_t m_size;
+ size_t m_order;
+ std::vector<size_t> m_bc;
+};
+
+
+// this model generates Bezier curves
+template< unsigned int N >
+class LFMBezierCurve
+ : public LinearFittingModel< double, Point, BezierCurve<N> >
+{
+ public:
+ LFMBezierCurve( size_t _order )
+ : mob(_order)
+ {
+ }
+
+ void feed( VectorView & coeff, double t ) const
+ {
+ mob.feed(coeff, t);
+ }
+
+ size_t size() const
+ {
+ return mob.size();
+ }
+
+ void instance(BezierCurve<N> & bc, ConstMatrixView const& raw_data) const
+ {
+ Bezier bx(size()-1);
+ Bezier by(size()-1);
+ mob.instance(bx, raw_data.column_const_view(X));
+ mob.instance(by, raw_data.column_const_view(Y));
+ bc = BezierCurve<N>(bx, by);
+ }
+
+ private:
+ LFMBezier mob;
+};
+
+} // end namespace NL
+} // end namespace Geom
+
+
+#endif // _NL_FITTING_MODEL_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/numeric/fitting-tool.h b/src/2geom/numeric/fitting-tool.h new file mode 100644 index 000000000..edacc663a --- /dev/null +++ b/src/2geom/numeric/fitting-tool.h @@ -0,0 +1,532 @@ +/*
+ * Fitting Tools
+ *
+ * Authors:
+ * Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2008 authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+
+#ifndef _NL_FITTING_TOOL_H_
+#define _NL_FITTING_TOOL_H_
+
+
+#include <2geom/numeric/vector.h>
+#include <2geom/numeric/matrix.h>
+
+#include <2geom/point.h>
+
+#include <vector>
+
+
+namespace Geom { namespace NL {
+
+namespace detail {
+
+
+template< typename ModelT>
+class lsf_base
+{
+ public:
+ typedef ModelT model_type;
+ typedef typename model_type::parameter_type parameter_type;
+ typedef typename model_type::value_type value_type;
+
+ lsf_base( model_type const& _model, size_t forecasted_samples )
+ : m_model(_model),
+ m_total_samples(0),
+ m_matrix(forecasted_samples, m_model.size()),
+ m_psdinv_matrix(NULL)
+ {
+ }
+
+ // compute pseudo inverse
+ void update()
+ {
+ if (total_samples() == 0) return;
+ if (m_psdinv_matrix != NULL)
+ {
+ delete m_psdinv_matrix;
+ }
+ MatrixView mv(m_matrix, 0, 0, total_samples(), m_matrix.columns());
+ m_psdinv_matrix = new Matrix( pseudo_inverse(mv) );
+ assert(m_psdinv_matrix != NULL);
+ }
+
+ size_t total_samples() const
+ {
+ return m_total_samples;
+ }
+
+ bool is_full() const
+ {
+ return (total_samples() == m_matrix.rows());
+ }
+
+ void clear()
+ {
+ m_total_samples = 0;
+ }
+
+ virtual
+ ~lsf_base()
+ {
+ if (m_psdinv_matrix != NULL)
+ {
+ delete m_psdinv_matrix;
+ }
+ }
+
+ protected:
+ const model_type & m_model;
+ size_t m_total_samples;
+ Matrix m_matrix;
+ Matrix* m_psdinv_matrix;
+
+}; // end class lsf_base
+
+
+
+
+template< typename ModelT, typename ValueType = typename ModelT::value_type>
+class lsf_solution
+{
+};
+
+// a fitting process on samples with value of type double
+// produces a solution of type Vector
+template< typename ModelT>
+class lsf_solution<ModelT, double>
+ : public lsf_base<ModelT>
+{
+public:
+ typedef ModelT model_type;
+ typedef typename model_type::parameter_type parameter_type;
+ typedef typename model_type::value_type value_type;
+ typedef Vector solution_type;
+ typedef lsf_base<model_type> base_type;
+
+ using base_type::m_model;
+ using base_type::m_psdinv_matrix;
+ using base_type::total_samples;
+
+public:
+ lsf_solution<ModelT, double>( model_type const& _model,
+ size_t forecasted_samples )
+ : base_type(_model, forecasted_samples),
+ m_solution(_model.size())
+ {
+ }
+
+ template< typename VectorT >
+ solution_type& result(VectorT const& sample_values)
+ {
+ assert(sample_values.size() == total_samples());
+ ConstVectorView sv(sample_values);
+ m_solution = (*m_psdinv_matrix) * sv;
+ return m_solution;
+ }
+
+ // a comparison between old sample values and the new ones is performed
+ // in order to minimize computation
+ // prerequisite:
+ // old_sample_values.size() == new_sample_values.size()
+ // no update() call can be performed between two result invocations
+ template< typename VectorT >
+ solution_type& result( VectorT const& old_sample_values,
+ VectorT const& new_sample_values )
+ {
+ assert(old_sample_values.size() == total_samples());
+ assert(new_sample_values.size() == total_samples());
+ Vector diff(total_samples());
+ for (size_t i = 0; i < diff.size(); ++i)
+ {
+ diff[i] = new_sample_values[i] - old_sample_values[i];
+ }
+ Vector column(m_model.size());
+ Vector delta(m_model.size(), 0.0);
+ for (size_t i = 0; i < diff.size(); ++i)
+ {
+ if (diff[i] != 0)
+ {
+ column = m_psdinv_matrix->column_view(i);
+ column.scale(diff[i]);
+ delta += column;
+ }
+ }
+ m_solution += delta;
+ return m_solution;
+ }
+
+ solution_type& result()
+ {
+ return m_solution;
+ }
+
+private:
+ solution_type m_solution;
+
+}; // end class lsf_solution<ModelT, double>
+
+
+// a fitting process on samples with value of type Point
+// produces a solution of type Matrix (with 2 columns)
+template< typename ModelT>
+class lsf_solution<ModelT, Point>
+ : public lsf_base<ModelT>
+{
+public:
+ typedef ModelT model_type;
+ typedef typename model_type::parameter_type parameter_type;
+ typedef typename model_type::value_type value_type;
+ typedef Matrix solution_type;
+ typedef lsf_base<model_type> base_type;
+
+ using base_type::m_model;
+ using base_type::m_psdinv_matrix;
+ using base_type::total_samples;
+
+public:
+ lsf_solution<ModelT, Point>( model_type const& _model,
+ size_t forecasted_samples )
+ : base_type(_model, forecasted_samples),
+ m_solution(_model.size(), 2)
+ {
+ }
+
+ solution_type& result(std::vector<Point> const& sample_values)
+ {
+ assert(sample_values.size() == total_samples());
+ Matrix svm(total_samples(), 2);
+ for (size_t i = 0; i < total_samples(); ++i)
+ {
+ svm(i, X) = sample_values[i][X];
+ svm(i, Y) = sample_values[i][Y];
+ }
+ m_solution = (*m_psdinv_matrix) * svm;
+ return m_solution;
+ }
+
+ // a comparison between old sample values and the new ones is performed
+ // in order to minimize computation
+ // prerequisite:
+ // old_sample_values.size() == new_sample_values.size()
+ // no update() call can to be performed between two result invocations
+ solution_type& result( std::vector<Point> const& old_sample_values,
+ std::vector<Point> const& new_sample_values )
+ {
+ assert(old_sample_values.size() == total_samples());
+ assert(new_sample_values.size() == total_samples());
+ Matrix diff(total_samples(), 2);
+ for (size_t i = 0; i < total_samples(); ++i)
+ {
+ diff(i, X) = new_sample_values[i][X] - old_sample_values[i][X];
+ diff(i, Y) = new_sample_values[i][Y] - old_sample_values[i][Y];
+ }
+ Vector column(m_model.size());
+ Matrix delta(m_model.size(), 2, 0.0);
+ VectorView deltax = delta.column_view(X);
+ VectorView deltay = delta.column_view(Y);
+ for (size_t i = 0; i < total_samples(); ++i)
+ {
+ if (diff(i, X) != 0)
+ {
+ column = m_psdinv_matrix->column_view(i);
+ column.scale(diff(i, X));
+ deltax += column;
+ }
+ if (diff(i, Y) != 0)
+ {
+ column = m_psdinv_matrix->column_view(i);
+ column.scale(diff(i, Y));
+ deltay += column;
+ }
+ }
+ m_solution += delta;
+ return m_solution;
+ }
+
+ solution_type& result()
+ {
+ return m_solution;
+ }
+
+private:
+ solution_type m_solution;
+
+}; // end class lsf_solution<ModelT, Point>
+
+
+
+
+template< typename ModelT,
+ bool WITH_FIXED_TERMS = ModelT::WITH_FIXED_TERMS >
+class lsf_with_fixed_terms
+{
+};
+
+
+// fitting tool for completely unknown models
+template< typename ModelT>
+class lsf_with_fixed_terms<ModelT, false>
+ : public lsf_solution<ModelT>
+{
+ public:
+ typedef ModelT model_type;
+ typedef typename model_type::parameter_type parameter_type;
+ typedef typename model_type::value_type value_type;
+ typedef lsf_solution<model_type> base_type;
+ typedef typename base_type::solution_type solution_type;
+
+ using base_type::total_samples;
+ using base_type::is_full;
+ using base_type::m_matrix;
+ using base_type::m_total_samples;
+ using base_type::m_model;
+
+ public:
+ lsf_with_fixed_terms<ModelT, false>( model_type const& _model,
+ size_t forecasted_samples )
+ : base_type(_model, forecasted_samples)
+ {
+ }
+
+ void append(parameter_type const& sample_parameter)
+ {
+ assert(!is_full());
+ VectorView row = m_matrix.row_view(total_samples());
+ m_model.feed(row, sample_parameter);
+ ++m_total_samples;
+ }
+
+ void append_copy(size_t sample_index)
+ {
+ assert(!is_full());
+ assert(sample_index < total_samples());
+ VectorView dest_row = m_matrix.row_view(total_samples());
+ VectorView source_row = m_matrix.row_view(sample_index);
+ dest_row = source_row;
+ ++m_total_samples;
+ }
+
+}; // end class lsf_with_fixed_terms<ModelT, false>
+
+
+// fitting tool for partially known models
+template< typename ModelT>
+class lsf_with_fixed_terms<ModelT, true>
+ : public lsf_solution<ModelT>
+{
+ public:
+ typedef ModelT model_type;
+ typedef typename model_type::parameter_type parameter_type;
+ typedef typename model_type::value_type value_type;
+ typedef lsf_solution<model_type> base_type;
+ typedef typename base_type::solution_type solution_type;
+
+ using base_type::total_samples;
+ using base_type::is_full;
+ using base_type::m_matrix;
+ using base_type::m_total_samples;
+ using base_type::m_model;
+
+ public:
+ lsf_with_fixed_terms<ModelT, true>( model_type const& _model,
+ size_t forecasted_samples )
+ : base_type(_model, forecasted_samples),
+ m_vector(forecasted_samples),
+ m_vector_view(NULL)
+ {
+ }
+ void append(parameter_type const& sample_parameter)
+ {
+ assert(!is_full());
+ VectorView row = m_matrix.row_view(total_samples());
+ m_model.feed(row, m_vector[total_samples()], sample_parameter);
+ ++m_total_samples;
+ }
+
+ void append_copy(size_t sample_index)
+ {
+ assert(!is_full());
+ assert(sample_index < total_samples());
+ VectorView dest_row = m_matrix.row_view(total_samples());
+ VectorView source_row = m_matrix.row_view(sample_index);
+ dest_row = source_row;
+ m_vector[total_samples()] = m_vector[sample_index];
+ ++m_total_samples;
+ }
+
+ void update()
+ {
+ base_type::update();
+ if (total_samples() == 0) return;
+ if (m_vector_view != NULL)
+ {
+ delete m_vector_view;
+ }
+ m_vector_view = new VectorView(m_vector, base_type::total_samples());
+ assert(m_vector_view != NULL);
+ }
+
+ virtual
+ ~lsf_with_fixed_terms<model_type, true>()
+ {
+ if (m_vector_view != NULL)
+ {
+ delete m_vector_view;
+ }
+ }
+
+ protected:
+ Vector m_vector;
+ VectorView* m_vector_view;
+
+}; // end class lsf_with_fixed_terms<ModelT, true>
+
+
+} // end namespace detail
+
+
+
+
+template< typename ModelT,
+ typename ValueType = typename ModelT::value_type,
+ bool WITH_FIXED_TERMS = ModelT::WITH_FIXED_TERMS >
+class least_squeares_fitter
+{
+};
+
+
+template< typename ModelT, typename ValueType >
+class least_squeares_fitter<ModelT, ValueType, false>
+ : public detail::lsf_with_fixed_terms<ModelT>
+{
+ public:
+ typedef ModelT model_type;
+ typedef detail::lsf_with_fixed_terms<model_type> base_type;
+ typedef typename base_type::parameter_type parameter_type;
+ typedef typename base_type::value_type value_type;
+ typedef typename base_type::solution_type solution_type;
+
+ public:
+ least_squeares_fitter<ModelT, ValueType, false>( model_type const& _model,
+ size_t forecasted_samples )
+ : base_type(_model, forecasted_samples)
+ {
+ }
+}; // end class least_squeares_fitter<ModelT, ValueType, true>
+
+
+template< typename ModelT>
+class least_squeares_fitter<ModelT, double, true>
+ : public detail::lsf_with_fixed_terms<ModelT>
+{
+ public:
+ typedef ModelT model_type;
+ typedef detail::lsf_with_fixed_terms<model_type> base_type;
+ typedef typename base_type::parameter_type parameter_type;
+ typedef typename base_type::value_type value_type;
+ typedef typename base_type::solution_type solution_type;
+
+ using base_type::m_vector_view;
+ using base_type::result;
+
+ public:
+ least_squeares_fitter<ModelT, double, true>( model_type const& _model,
+ size_t forecasted_samples )
+ : base_type(_model, forecasted_samples)
+ {
+ }
+
+ template< typename VectorT >
+ solution_type& result(VectorT const& sample_values)
+ {
+ assert(sample_values.size() == m_vector_view->size());
+ Vector sv(sample_values.size());
+ for (size_t i = 0; i < sv.size(); ++i)
+ sv[i] = sample_values[i] - (*m_vector_view)[i];
+ return base_type::result(sv);
+ }
+
+}; // end class least_squeares_fitter<ModelT, double, true>
+
+
+template< typename ModelT>
+class least_squeares_fitter<ModelT, Point, true>
+ : public detail::lsf_with_fixed_terms<ModelT>
+{
+ public:
+ typedef ModelT model_type;
+ typedef detail::lsf_with_fixed_terms<model_type> base_type;
+ typedef typename base_type::parameter_type parameter_type;
+ typedef typename base_type::value_type value_type;
+ typedef typename base_type::solution_type solution_type;
+
+ using base_type::m_vector_view;
+ using base_type::result;
+
+ public:
+ least_squeares_fitter<ModelT, Point, true>( model_type const& _model,
+ size_t forecasted_samples )
+ : base_type(_model, forecasted_samples)
+ {
+ }
+
+ solution_type& result(std::vector<Point> const& sample_values)
+ {
+ assert(sample_values.size() == m_vector_view->size());
+ NL::Matrix sv(sample_values.size(), 2);
+ for (size_t i = 0; i < sample_values.size(); ++i)
+ {
+ sv(i, X) = sample_values[i][X] - (*m_vector_view)[i];
+ sv(i, Y) = sample_values[i][Y] - (*m_vector_view)[i];
+ }
+ return base_type::result(sv);
+ }
+
+}; // end class least_squeares_fitter<ModelT, Point, true>
+
+
+} // end namespace NL
+} // end namespace Geom
+
+
+
+#endif // _NL_FITTING_TOOL_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/numeric/linear_system.h b/src/2geom/numeric/linear_system.h index 7699c5224..5b516c9e6 100644 --- a/src/2geom/numeric/linear_system.h +++ b/src/2geom/numeric/linear_system.h @@ -1,89 +1,138 @@ -#ifndef _NL_LINEAR_SYSTEM_H_ -#define _NL_LINEAR_SYSTEM_H_ - - -#include <cassert> - -#include <gsl/gsl_linalg.h> - -#include "2geom/numeric/matrix.h" -#include "2geom/numeric/vector.h" - - -namespace Geom { namespace NL { - - -class LinearSystem -{ -public: - LinearSystem(Matrix & _matrix, Vector & _vector) - : m_matrix(_matrix), m_vector(_vector), m_solution(_matrix.columns()) - { - } - - const Vector & LU_solve() - { - assert( matrix().rows() == matrix().columns() - && matrix().rows() == vector().size() ); - int s; - gsl_permutation * p = gsl_permutation_alloc(matrix().rows()); - gsl_linalg_LU_decomp (matrix().get_gsl_matrix(), p, &s); - gsl_linalg_LU_solve( matrix().get_gsl_matrix(), - p, - vector().get_gsl_vector(), - m_solution.get_gsl_vector() - ); - gsl_permutation_free(p); - return solution(); - } - - const Vector & SV_solve() - { - assert( matrix().rows() >= matrix().columns() - && matrix().rows() == vector().size() ); - - gsl_matrix* U = matrix().get_gsl_matrix(); - gsl_matrix* V = gsl_matrix_alloc(matrix().columns(), matrix().columns()); - gsl_vector* S = gsl_vector_alloc(matrix().columns()); - gsl_vector* work = gsl_vector_alloc(matrix().columns()); - - gsl_linalg_SV_decomp( U, V, S, work ); - - gsl_vector* b = vector().get_gsl_vector(); - gsl_vector* x = m_solution.get_gsl_vector(); - - gsl_linalg_SV_solve( U, V, S, b, x); - - gsl_matrix_free(V); - gsl_vector_free(S); - gsl_vector_free(work); - - return solution(); - } - - Matrix & matrix() - { - return m_matrix; - } - - Vector & vector() - { - return m_vector; - } - - const Vector & solution() const - { - return m_solution; - } - -private: - Matrix & m_matrix; - Vector & m_vector; - Vector m_solution; -}; - - -} } // end namespaces - - -#endif /*_NL_LINEAR_SYSTEM_H_*/ +/*
+ * LinearSystem class wraps some gsl routines for solving linear systems
+ *
+ * Authors:
+ * Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2008 authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+
+#ifndef _NL_LINEAR_SYSTEM_H_
+#define _NL_LINEAR_SYSTEM_H_
+
+
+#include <cassert>
+
+#include <gsl/gsl_linalg.h>
+
+#include <2geom/numeric/matrix.h>
+#include <2geom/numeric/vector.h>
+
+
+namespace Geom { namespace NL {
+
+
+class LinearSystem
+{
+public:
+ LinearSystem(MatrixView & _matrix, VectorView & _vector)
+ : m_matrix(_matrix), m_vector(_vector), m_solution(_matrix.columns())
+ {
+ }
+
+ LinearSystem(Matrix & _matrix, Vector & _vector)
+ : m_matrix(_matrix), m_vector(_vector), m_solution(_matrix.columns())
+ {
+ }
+
+ const Vector & LU_solve()
+ {
+ assert( matrix().rows() == matrix().columns()
+ && matrix().rows() == vector().size() );
+ int s;
+ gsl_permutation * p = gsl_permutation_alloc(matrix().rows());
+ gsl_linalg_LU_decomp (matrix().get_gsl_matrix(), p, &s);
+ gsl_linalg_LU_solve( matrix().get_gsl_matrix(),
+ p,
+ vector().get_gsl_vector(),
+ m_solution.get_gsl_vector()
+ );
+ gsl_permutation_free(p);
+ return solution();
+ }
+
+ const Vector & SV_solve()
+ {
+ assert( matrix().rows() >= matrix().columns()
+ && matrix().rows() == vector().size() );
+
+ gsl_matrix* U = matrix().get_gsl_matrix();
+ gsl_matrix* V = gsl_matrix_alloc(matrix().columns(), matrix().columns());
+ gsl_vector* S = gsl_vector_alloc(matrix().columns());
+ gsl_vector* work = gsl_vector_alloc(matrix().columns());
+
+ gsl_linalg_SV_decomp( U, V, S, work );
+
+ gsl_vector* b = vector().get_gsl_vector();
+ gsl_vector* x = m_solution.get_gsl_vector();
+
+ gsl_linalg_SV_solve( U, V, S, b, x);
+
+ gsl_matrix_free(V);
+ gsl_vector_free(S);
+ gsl_vector_free(work);
+
+ return solution();
+ }
+
+ MatrixView & matrix()
+ {
+ return m_matrix;
+ }
+
+ VectorView & vector()
+ {
+ return m_vector;
+ }
+
+ const Vector & solution() const
+ {
+ return m_solution;
+ }
+
+private:
+ MatrixView m_matrix;
+ VectorView m_vector;
+ Vector m_solution;
+};
+
+
+} } // end namespaces
+
+
+#endif /*_NL_LINEAR_SYSTEM_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/numeric/matrix.cpp b/src/2geom/numeric/matrix.cpp new file mode 100644 index 000000000..dde80dfb6 --- /dev/null +++ b/src/2geom/numeric/matrix.cpp @@ -0,0 +1,115 @@ +/* + * Matrix, MatrixView, ConstMatrixView classes wrap the gsl matrix routines; + * "views" mimic the semantic of C++ references: any operation performed + * on a "view" is actually performed on the "viewed object" + * + * Authors: + * Marco Cecchetti <mrcekets at gmail.com> + * + * Copyright 2008 authors + * + * This library is free software; you can redistribute it and/or + * modify it either under the terms of the GNU Lesser General Public + * License version 2.1 as published by the Free Software Foundation + * (the "LGPL") or, at your option, under the terms of the Mozilla + * Public License Version 1.1 (the "MPL"). If you do not alter this + * notice, a recipient may use your version of this file under either + * the MPL or the LGPL. + * + * You should have received a copy of the LGPL along with this library + * in the file COPYING-LGPL-2.1; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * You should have received a copy of the MPL along with this library + * in the file COPYING-MPL-1.1 + * + * The contents of this file are subject to the Mozilla Public License + * Version 1.1 (the "License"); you may not use this file except in + * compliance with the License. You may obtain a copy of the License at + * http://www.mozilla.org/MPL/ + * + * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY + * OF ANY KIND, either express or implied. See the LGPL or the MPL for + * the specific language governing rights and limitations. + */ + + +#include <numeric/matrix.h> +#include <numeric/vector.h> + + + +namespace Geom { namespace NL { + + +Vector operator*( detail::BaseMatrixImpl const& A, + detail::BaseVectorImpl const& v ) +{ + assert(A.columns() == v.size()); + + Vector result(A.rows(), 0.0); + for (size_t i = 0; i < A.rows(); ++i) + for (size_t j = 0; j < A.columns(); ++j) + result[i] += A(i,j) * v[j]; + + return result; +} + +Matrix operator*( detail::BaseMatrixImpl const& A, + detail::BaseMatrixImpl const& B ) +{ + assert(A.columns() == B.rows()); + + Matrix C(A.rows(), B.columns(), 0.0); + for (size_t i = 0; i < C.rows(); ++i) + for (size_t j = 0; j < C.columns(); ++j) + for (size_t k = 0; k < A.columns(); ++k) + C(i,j) += A(i,k) * B(k, j); + + return C; +} + +Matrix pseudo_inverse(detail::BaseMatrixImpl const& A) +{ + + Matrix U(A); + Matrix V(A.columns(), A.columns()); + Vector s(A.columns()); + gsl_vector* work = gsl_vector_alloc(A.columns()); + + gsl_linalg_SV_decomp( U.get_gsl_matrix(), + V.get_gsl_matrix(), + s.get_gsl_vector(), + work ); + + Matrix P(A.columns(), A.rows(), 0.0); + + int sz = s.size(); + while ( sz-- > 0 && s[sz] == 0 ) {} + ++sz; + if (sz == 0) return P; + VectorView sv(s, sz); + + for (size_t i = 0; i < sv.size(); ++i) + { + VectorView v = V.column_view(i); + v.scale(1/sv[i]); + for (size_t h = 0; h < P.rows(); ++h) + for (size_t k = 0; k < P.columns(); ++k) + P(h,k) += V(h,i) * U(k,i); + } + + return P; +} + +} } // end namespaces + +/* + 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/numeric/matrix.h b/src/2geom/numeric/matrix.h index bdd647e53..64557a6f1 100644 --- a/src/2geom/numeric/matrix.h +++ b/src/2geom/numeric/matrix.h @@ -1,207 +1,555 @@ -#ifndef _NL_MATRIX_H_ -#define _NL_MATRIX_H_ - - -#include <cassert> -#include <utility> -#include <algorithm> - -#include <gsl/gsl_matrix.h> - - - -namespace Geom { namespace NL { - -class Matrix; -void swap(Matrix & m1, Matrix & m2); - - -class Matrix -{ -public: - // the matrix is not inizialized - Matrix(size_t n1, size_t n2) - : m_rows(n1), m_columns(n2) - { - m_matrix = gsl_matrix_alloc(n1, n2); - } - - Matrix(size_t n1, size_t n2, double x) - :m_rows(n1), m_columns(n2) - { - m_matrix = gsl_matrix_alloc(n1, n2); - gsl_matrix_set_all(m_matrix, x); - } - - Matrix( Matrix const& _matrix ) - : m_rows(_matrix.rows()), m_columns(_matrix.columns()) - { - m_matrix = gsl_matrix_alloc(rows(), columns()); - gsl_matrix_memcpy(m_matrix, _matrix.m_matrix); - } - -// explicit -// Matrix( gsl_matrix* m, size_t n1, size_t n2) -// : m_rows(n1), m_columns(n2), m_matrix(m) -// { -// } - - Matrix & operator=(Matrix const& _matrix) - { - assert( rows() == _matrix.rows() && columns() == _matrix.columns() ); - gsl_matrix_memcpy(m_matrix, _matrix.m_matrix); - - return *this; - } - - virtual ~Matrix() - { - gsl_matrix_free(m_matrix); - } - - void set_all( double x = 0) - { - gsl_matrix_set_all(m_matrix, x); - } - - void set_identity() - { - gsl_matrix_set_identity(m_matrix); - } - - const double & operator() (size_t i, size_t j) const - { - return *gsl_matrix_const_ptr(m_matrix, i, j); - } - - double & operator() (size_t i, size_t j) - { - return *gsl_matrix_ptr(m_matrix, i, j); - } - - gsl_matrix* get_gsl_matrix() - { - return m_matrix; - } - - void swap_rows(size_t i, size_t j) - { - gsl_matrix_swap_rows(m_matrix, i, j); - } - - void swap_columns(size_t i, size_t j) - { - gsl_matrix_swap_columns(m_matrix, i, j); - } - - Matrix & transpose() - { - gsl_matrix_transpose(m_matrix); - return (*this); - } - - Matrix & scale(double x) - { - gsl_matrix_scale(m_matrix, x); - return (*this); - } - - Matrix & translate(double x) - { - gsl_matrix_add_constant(m_matrix, x); - return (*this); - } - - Matrix & operator+=(Matrix const& _matrix) - { - gsl_matrix_add(m_matrix, _matrix.m_matrix); - return (*this); - } - - Matrix & operator-=(Matrix const& _matrix) - { - gsl_matrix_sub(m_matrix, _matrix.m_matrix); - return (*this); - } - - bool is_zero() const - { - return gsl_matrix_isnull(m_matrix); - } - - bool is_positive() const - { - return gsl_matrix_ispos(m_matrix); - } - - bool is_negative() const - { - return gsl_matrix_isneg(m_matrix); - } - - bool is_non_negative() const - { - for ( unsigned int i = 0; i < rows(); ++i ) - { - for ( unsigned int j = 0; j < columns(); ++j ) - { - if ( (*this)(i,j) < 0 ) return false; - } - } - return true; - } - - double max() const - { - return gsl_matrix_max(m_matrix); - } - - double min() const - { - return gsl_matrix_min(m_matrix); - } - - std::pair<size_t, size_t> - max_index() const - { - std::pair<size_t, size_t> indices; - gsl_matrix_max_index(m_matrix, &(indices.first), &(indices.second)); - return indices; - } - - std::pair<size_t, size_t> - min_index() const - { - std::pair<size_t, size_t> indices; - gsl_matrix_min_index(m_matrix, &(indices.first), &(indices.second)); - return indices; - } - - friend - void swap(Matrix & m1, Matrix & m2); - - size_t rows() const - { - return m_rows; - } - - size_t columns() const - { - return m_columns; - } - -private: - size_t m_rows, m_columns; - gsl_matrix* m_matrix; -}; - -void swap(Matrix & m1, Matrix & m2) -{ - assert( m1.rows() == m2.rows() && m1.columns() == m2.columns() ); - std::swap(m1.m_matrix, m2.m_matrix); -} - - -} } // end namespaces - -#endif /*_NL_MATRIX_H_*/ +/*
+ * Matrix, MatrixView, ConstMatrixView classes wrap the gsl matrix routines;
+ * "views" mimic the semantic of C++ references: any operation performed
+ * on a "view" is actually performed on the "viewed object"
+ *
+ * Authors:
+ * Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2008 authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+
+
+
+#ifndef _NL_MATRIX_H_
+#define _NL_MATRIX_H_
+
+#include <2geom/numeric/vector.h>
+
+#include <cassert>
+#include <utility> // for std::pair
+#include <algorithm> // for std::swap
+#include <sstream>
+#include <string>
+
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_linalg.h>
+
+
+namespace Geom { namespace NL {
+
+namespace detail
+{
+
+class BaseMatrixImpl
+{
+ public:
+ virtual ~BaseMatrixImpl()
+ {
+ }
+
+ ConstVectorView row_const_view(size_t i) const
+ {
+ return ConstVectorView(gsl_matrix_const_row(m_matrix, i));
+ }
+
+ ConstVectorView column_const_view(size_t i) const
+ {
+ return ConstVectorView(gsl_matrix_const_column(m_matrix, i));
+ }
+
+ const double & operator() (size_t i, size_t j) const
+ {
+ return *gsl_matrix_const_ptr(m_matrix, i, j);
+ }
+
+ const gsl_matrix* get_gsl_matrix() const
+ {
+ return m_matrix;
+ }
+
+ bool is_zero() const
+ {
+ return gsl_matrix_isnull(m_matrix);
+ }
+
+ bool is_positive() const
+ {
+ return gsl_matrix_ispos(m_matrix);
+ }
+
+ bool is_negative() const
+ {
+ return gsl_matrix_isneg(m_matrix);
+ }
+
+ bool is_non_negative() const
+ {
+ for ( unsigned int i = 0; i < rows(); ++i )
+ {
+ for ( unsigned int j = 0; j < columns(); ++j )
+ {
+ if ( (*this)(i,j) < 0 ) return false;
+ }
+ }
+ return true;
+ }
+
+ double max() const
+ {
+ return gsl_matrix_max(m_matrix);
+ }
+
+ double min() const
+ {
+ return gsl_matrix_min(m_matrix);
+ }
+
+ std::pair<size_t, size_t>
+ max_index() const
+ {
+ std::pair<size_t, size_t> indices;
+ gsl_matrix_max_index(m_matrix, &(indices.first), &(indices.second));
+ return indices;
+ }
+
+ std::pair<size_t, size_t>
+ min_index() const
+ {
+ std::pair<size_t, size_t> indices;
+ gsl_matrix_min_index(m_matrix, &(indices.first), &(indices.second));
+ return indices;
+ }
+
+ size_t rows() const
+ {
+ return m_rows;
+ }
+
+ size_t columns() const
+ {
+ return m_columns;
+ }
+
+ std::string str() const;
+
+ protected:
+ size_t m_rows, m_columns;
+ gsl_matrix* m_matrix;
+
+}; // end class BaseMatrixImpl
+
+
+inline
+bool operator== (BaseMatrixImpl const& m1, BaseMatrixImpl const& m2)
+{
+ if (m1.rows() != m2.rows() || m1.columns() != m2.columns()) return false;
+
+ for (size_t i = 0; i < m1.rows(); ++i)
+ for (size_t j = 0; j < m1.columns(); ++j)
+ if (m1(i,j) != m2(i,j)) return false;
+
+ return true;
+}
+
+template< class charT >
+inline
+std::basic_ostream<charT> &
+operator<< (std::basic_ostream<charT> & os, const BaseMatrixImpl & _matrix)
+{
+ if (_matrix.rows() == 0 || _matrix.columns() == 0) return os;
+
+ os << "[[" << _matrix(0,0);
+ for (size_t j = 1; j < _matrix.columns(); ++j)
+ {
+ os << ", " << _matrix(0,j);
+ }
+ os << "]";
+
+ for (size_t i = 1; i < _matrix.rows(); ++i)
+ {
+ os << ", [" << _matrix(i,0);
+ for (size_t j = 1; j < _matrix.columns(); ++j)
+ {
+ os << ", " << _matrix(i,j);
+ }
+ os << "]";
+ }
+ os << "]";
+ return os;
+}
+
+inline
+std::string BaseMatrixImpl::str() const
+{
+ std::ostringstream oss;
+ oss << (*this);
+ return oss.str();
+}
+
+
+class MatrixImpl : public BaseMatrixImpl
+{
+ public:
+
+ typedef BaseMatrixImpl base_type;
+
+ void set_all( double x )
+ {
+ gsl_matrix_set_all(m_matrix, x);
+ }
+
+ void set_identity()
+ {
+ gsl_matrix_set_identity(m_matrix);
+ }
+
+ using base_type::operator();
+
+ double & operator() (size_t i, size_t j)
+ {
+ return *gsl_matrix_ptr(m_matrix, i, j);
+ }
+
+ using base_type::get_gsl_matrix;
+
+ gsl_matrix* get_gsl_matrix()
+ {
+ return m_matrix;
+ }
+
+ VectorView row_view(size_t i)
+ {
+ return VectorView(gsl_matrix_row(m_matrix, i));
+ }
+
+ VectorView column_view(size_t i)
+ {
+ return VectorView(gsl_matrix_column(m_matrix, i));
+ }
+
+ void swap_rows(size_t i, size_t j)
+ {
+ gsl_matrix_swap_rows(m_matrix, i, j);
+ }
+
+ void swap_columns(size_t i, size_t j)
+ {
+ gsl_matrix_swap_columns(m_matrix, i, j);
+ }
+
+ MatrixImpl & transpose()
+ {
+ assert(columns() == rows());
+ gsl_matrix_transpose(m_matrix);
+ return (*this);
+ }
+
+ MatrixImpl & scale(double x)
+ {
+ gsl_matrix_scale(m_matrix, x);
+ return (*this);
+ }
+
+ MatrixImpl & translate(double x)
+ {
+ gsl_matrix_add_constant(m_matrix, x);
+ return (*this);
+ }
+
+ MatrixImpl & operator+=(base_type const& _matrix)
+ {
+ gsl_matrix_add(m_matrix, _matrix.get_gsl_matrix());
+ return (*this);
+ }
+
+ MatrixImpl & operator-=(base_type const& _matrix)
+ {
+ gsl_matrix_sub(m_matrix, _matrix.get_gsl_matrix());
+ return (*this);
+ }
+
+}; // end class MatrixImpl
+
+} // end namespace detail
+
+
+using detail::operator==;
+using detail::operator<<;
+
+
+
+
+class Matrix: public detail::MatrixImpl
+{
+ public:
+ typedef detail::MatrixImpl base_type;
+
+ public:
+ // the matrix is not inizialized
+ Matrix(size_t n1, size_t n2)
+ {
+ m_rows = n1;
+ m_columns = n2;
+ m_matrix = gsl_matrix_alloc(n1, n2);
+ }
+
+ Matrix(size_t n1, size_t n2, double x)
+ {
+ m_rows = n1;
+ m_columns = n2;
+ m_matrix = gsl_matrix_alloc(n1, n2);
+ gsl_matrix_set_all(m_matrix, x);
+ }
+
+ Matrix(Matrix const& _matrix)
+ : base_type()
+ {
+ m_rows = _matrix.rows();
+ m_columns = _matrix.columns();
+ m_matrix = gsl_matrix_alloc(rows(), columns());
+ gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());
+ }
+
+ explicit
+ Matrix(base_type::base_type const& _matrix)
+ {
+ m_rows = _matrix.rows();
+ m_columns = _matrix.columns();
+ m_matrix = gsl_matrix_alloc(rows(), columns());
+ gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());
+ }
+
+ Matrix & operator=(Matrix const& _matrix)
+ {
+ assert( rows() == _matrix.rows() && columns() == _matrix.columns() );
+ gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());
+ return *this;
+ }
+
+ Matrix & operator=(base_type::base_type const& _matrix)
+ {
+ assert( rows() == _matrix.rows() && columns() == _matrix.columns() );
+ gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());
+ return *this;
+ }
+
+ virtual ~Matrix()
+ {
+ gsl_matrix_free(m_matrix);
+ }
+
+ Matrix & transpose()
+ {
+ return static_cast<Matrix &>( base_type::transpose() );
+ }
+
+ Matrix & scale(double x)
+ {
+ return static_cast<Matrix &>( base_type::scale(x) );
+ }
+
+ Matrix & translate(double x)
+ {
+ return static_cast<Matrix &>( base_type::translate(x) );
+ }
+
+ Matrix & operator+=(base_type::base_type const& _matrix)
+ {
+ return static_cast<Matrix &>( base_type::operator+=(_matrix) );
+ }
+
+ Matrix & operator-=(base_type::base_type const& _matrix)
+ {
+ return static_cast<Matrix &>( base_type::operator-=(_matrix) );
+ }
+
+ friend
+ void swap(Matrix & m1, Matrix & m2);
+ friend
+ void swap_any(Matrix & m1, Matrix & m2);
+
+}; // end class Matrix
+
+
+// warning! this operation invalidates any view of the passed matrix objects
+inline
+void swap(Matrix & m1, Matrix & m2)
+{
+ assert( m1.rows() == m2.rows() && m1.columns() == m2.columns() );
+ std::swap(m1.m_matrix, m2.m_matrix);
+}
+
+inline
+void swap_any(Matrix & m1, Matrix & m2)
+{
+ std::swap(m1.m_matrix, m2.m_matrix);
+ std::swap(m1.m_rows, m2.m_rows);
+ std::swap(m1.m_columns, m2.m_columns);
+}
+
+
+
+class ConstMatrixView : public detail::BaseMatrixImpl
+{
+ public:
+ typedef detail::BaseMatrixImpl base_type;
+
+ public:
+ ConstMatrixView(const base_type & _matrix, size_t k1, size_t k2, size_t n1, size_t n2)
+ : m_matrix_view( gsl_matrix_const_submatrix(_matrix.get_gsl_matrix(), k1, k2, n1, n2) )
+ {
+ m_rows = n1;
+ m_columns = n2;
+ m_matrix = const_cast<gsl_matrix*>( &(m_matrix_view.matrix) );
+ }
+
+ ConstMatrixView(const ConstMatrixView & _matrix)
+ : base_type(),
+ m_matrix_view(_matrix.m_matrix_view)
+ {
+ m_rows = _matrix.rows();
+ m_columns = _matrix.columns();
+ m_matrix = const_cast<gsl_matrix*>( &(m_matrix_view.matrix) );
+ }
+
+ ConstMatrixView(const base_type & _matrix)
+ : m_matrix_view(gsl_matrix_const_submatrix(_matrix.get_gsl_matrix(), 0, 0, _matrix.rows(), _matrix.columns()))
+ {
+ m_rows = _matrix.rows();
+ m_columns = _matrix.columns();
+ m_matrix = const_cast<gsl_matrix*>( &(m_matrix_view.matrix) );
+ }
+
+ private:
+ gsl_matrix_const_view m_matrix_view;
+
+}; // end class ConstMatrixView
+
+
+
+
+class MatrixView : public detail::MatrixImpl
+{
+ public:
+ typedef detail::MatrixImpl base_type;
+
+ public:
+ MatrixView(base_type & _matrix, size_t k1, size_t k2, size_t n1, size_t n2)
+ {
+ m_rows = n1;
+ m_columns = n2;
+ m_matrix_view
+ = gsl_matrix_submatrix(_matrix.get_gsl_matrix(), k1, k2, n1, n2);
+ m_matrix = &(m_matrix_view.matrix);
+ }
+
+ MatrixView(const MatrixView & _matrix)
+ : base_type()
+ {
+ m_rows = _matrix.rows();
+ m_columns = _matrix.columns();
+ m_matrix_view = _matrix.m_matrix_view;
+ m_matrix = &(m_matrix_view.matrix);
+ }
+
+ MatrixView(Matrix & _matrix)
+ {
+ m_rows = _matrix.rows();
+ m_columns = _matrix.columns();
+ m_matrix_view
+ = gsl_matrix_submatrix(_matrix.get_gsl_matrix(), 0, 0, rows(), columns());
+ m_matrix = &(m_matrix_view.matrix);
+ }
+
+ MatrixView & operator=(MatrixView const& _matrix)
+ {
+ assert( rows() == _matrix.rows() && columns() == _matrix.columns() );
+ gsl_matrix_memcpy(m_matrix, _matrix.m_matrix);
+ return *this;
+ }
+
+ MatrixView & operator=(base_type::base_type const& _matrix)
+ {
+ assert( rows() == _matrix.rows() && columns() == _matrix.columns() );
+ gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());
+ return *this;
+ }
+
+ MatrixView & transpose()
+ {
+ return static_cast<MatrixView &>( base_type::transpose() );
+ }
+
+ MatrixView & scale(double x)
+ {
+ return static_cast<MatrixView &>( base_type::scale(x) );
+ }
+
+ MatrixView & translate(double x)
+ {
+ return static_cast<MatrixView &>( base_type::translate(x) );
+ }
+
+ MatrixView & operator+=(base_type::base_type const& _matrix)
+ {
+ return static_cast<MatrixView &>( base_type::operator+=(_matrix) );
+ }
+
+ MatrixView & operator-=(base_type::base_type const& _matrix)
+ {
+ return static_cast<MatrixView &>( base_type::operator-=(_matrix) );
+ }
+
+ friend
+ void swap_view(MatrixView & m1, MatrixView & m2);
+
+ private:
+ gsl_matrix_view m_matrix_view;
+
+}; // end class MatrixView
+
+
+inline
+void swap_view(MatrixView & m1, MatrixView & m2)
+{
+ assert( m1.rows() == m2.rows() && m1.columns() == m2.columns() );
+ std::swap(m1.m_matrix_view, m2.m_matrix_view);
+}
+
+Vector operator*( detail::BaseMatrixImpl const& A,
+ detail::BaseVectorImpl const& v );
+
+Matrix operator*( detail::BaseMatrixImpl const& A,
+ detail::BaseMatrixImpl const& B );
+
+Matrix pseudo_inverse(detail::BaseMatrixImpl const& A);
+
+} } // end namespaces
+
+#endif /*_NL_MATRIX_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/numeric/vector.h b/src/2geom/numeric/vector.h index b25861e22..43a39a1ac 100644 --- a/src/2geom/numeric/vector.h +++ b/src/2geom/numeric/vector.h @@ -1,185 +1,565 @@ -#ifndef _NL_VECTOR_H_ -#define _NL_VECTOR_H_ - -#include <cassert> -#include <utility> - -#include <gsl/gsl_vector.h> - - -namespace Geom { namespace NL { - -class Vector; -void swap(Vector & v1, Vector & v2); - - -class Vector -{ -public: - Vector(size_t n) - : m_size(n) - { - m_vector = gsl_vector_alloc(n); - } - - Vector(size_t n, double x) - : m_size(n) - { - m_vector = gsl_vector_alloc(n); - gsl_vector_set_all(m_vector, x); - } - - // create a vector with n elements all set to zero - // but the i-th that is set to 1 - Vector(size_t n, size_t i) - : m_size(n) - { - m_vector = gsl_vector_alloc(n); - gsl_vector_set_basis(m_vector, i); - } - - Vector(Vector const& _vector) - : m_size(_vector.size()) - { - m_vector = gsl_vector_alloc(size()); - gsl_vector_memcpy(m_vector, _vector.m_vector); - } - - virtual ~Vector() - { - gsl_vector_free(m_vector); - } - - Vector & operator=(Vector const& _vector) - { - assert( size() == _vector.size() ); - gsl_vector_memcpy(m_vector, _vector.m_vector); - return (*this); - } - - void set_all(double x = 0) - { - gsl_vector_set_all(m_vector, x); - } - - void set_basis(size_t i) - { - gsl_vector_set_basis(m_vector, i); - } - - double const& operator[](size_t i) const - { - return *gsl_vector_const_ptr(m_vector, i); - } - - double & operator[](size_t i) - { - return *gsl_vector_ptr(m_vector, i); - } - - gsl_vector* get_gsl_vector() - { - return m_vector; - } - - void swap_elements(size_t i, size_t j) - { - gsl_vector_swap_elements(m_vector, i, j); - } - - void reverse() - { - gsl_vector_reverse(m_vector); - } - - Vector & scale(double x) - { - gsl_vector_scale(m_vector, x); - return (*this); - } - - Vector & translate(double x) - { - gsl_vector_add_constant(m_vector, x); - return (*this); - } - - Vector & operator+=(Vector const& _vector) - { - gsl_vector_add(m_vector, _vector.m_vector); - return (*this); - } - - Vector & operator-=(Vector const& _vector) - { - gsl_vector_sub(m_vector, _vector.m_vector); - return (*this); - } - - bool is_zero() const - { - return gsl_vector_isnull(m_vector); - } - - bool is_positive() const - { - return gsl_vector_ispos(m_vector); - } - - bool is_negative() const - { - return gsl_vector_isneg(m_vector); - } - - bool is_non_negative() const - { - for ( size_t i = 0; i < size(); ++i ) - { - if ( (*this)[i] < 0 ) return false; - } - return true; - } - - double max() const - { - return gsl_vector_max(m_vector); - } - - double min() const - { - return gsl_vector_min(m_vector); - } - - size_t max_index() const - { - return gsl_vector_max_index(m_vector); - } - - size_t min_index() const - { - return gsl_vector_min_index(m_vector); - } - - friend - void swap(Vector & v1, Vector & v2); - - size_t size() const - { - return m_size; - } - -private: - size_t m_size; - gsl_vector* m_vector; -}; - -void swap(Vector & v1, Vector & v2) -{ - assert( v1.size() == v2.size() ); - std::swap(v1.m_vector, v2.m_vector); -} - -} } // end namespaces - - -#endif /*_NL_VECTOR_H_*/ +/*
+ * Vector, VectorView, ConstVectorView classes wrap the gsl vector routines;
+ * "views" mimic the semantic of C++ references: any operation performed
+ * on a "view" is actually performed on the "viewed object"
+ *
+ * Authors:
+ * Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2008 authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+
+
+
+#ifndef _NL_VECTOR_H_
+#define _NL_VECTOR_H_
+
+#include <cassert>
+#include <algorithm> // for std::swap
+#include <vector>
+#include <sstream>
+#include <string>
+
+
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_blas.h>
+
+
+namespace Geom { namespace NL {
+
+namespace detail
+{
+
+class BaseVectorImpl
+{
+ public:
+ double const& operator[](size_t i) const
+ {
+ return *gsl_vector_const_ptr(m_vector, i);
+ }
+
+ const gsl_vector* get_gsl_vector() const
+ {
+ return m_vector;
+ }
+ bool is_zero() const
+ {
+ return gsl_vector_isnull(m_vector);
+ }
+
+ bool is_positive() const
+ {
+ return gsl_vector_ispos(m_vector);
+ }
+
+ bool is_negative() const
+ {
+ return gsl_vector_isneg(m_vector);
+ }
+
+ bool is_non_negative() const
+ {
+ for ( size_t i = 0; i < size(); ++i )
+ {
+ if ( (*this)[i] < 0 ) return false;
+ }
+ return true;
+ }
+
+ double max() const
+ {
+ return gsl_vector_max(m_vector);
+ }
+
+ double min() const
+ {
+ return gsl_vector_min(m_vector);
+ }
+
+ size_t max_index() const
+ {
+ return gsl_vector_max_index(m_vector);
+ }
+
+ size_t min_index() const
+ {
+ return gsl_vector_min_index(m_vector);
+ }
+
+ size_t size() const
+ {
+ return m_size;
+ }
+
+ std::string str() const;
+
+ virtual ~BaseVectorImpl()
+ {
+ }
+
+ protected:
+ size_t m_size;
+ gsl_vector* m_vector;
+
+}; // end class BaseVectorImpl
+
+
+inline
+bool operator== (BaseVectorImpl const& v1, BaseVectorImpl const& v2)
+{
+ if (v1.size() != v2.size()) return false;
+
+ for (size_t i = 0; i < v1.size(); ++i)
+ {
+ if (v1[i] != v2[i]) return false;
+ }
+ return true;
+}
+
+template< class charT >
+inline
+std::basic_ostream<charT> &
+operator<< (std::basic_ostream<charT> & os, const BaseVectorImpl & _vector)
+{
+ if (_vector.size() == 0 ) return os;
+ os << "[" << _vector[0];
+ for (unsigned int i = 1; i < _vector.size(); ++i)
+ {
+ os << ", " << _vector[i];
+ }
+ os << "]";
+ return os;
+}
+
+inline
+std::string BaseVectorImpl::str() const
+{
+ std::ostringstream oss;
+ oss << (*this);
+ return oss.str();
+}
+
+inline
+double dot(BaseVectorImpl const& v1, BaseVectorImpl const& v2)
+{
+ double result;
+ gsl_blas_ddot(v1.get_gsl_vector(), v2.get_gsl_vector(), &result);
+ return result;
+}
+
+
+class VectorImpl : public BaseVectorImpl
+{
+ public:
+ typedef BaseVectorImpl base_type;
+
+ public:
+ void set_all(double x)
+ {
+ gsl_vector_set_all(m_vector, x);
+ }
+
+ void set_basis(size_t i)
+ {
+ gsl_vector_set_basis(m_vector, i);
+ }
+
+ using base_type::operator[];
+
+ double & operator[](size_t i)
+ {
+ return *gsl_vector_ptr(m_vector, i);
+ }
+
+ using base_type::get_gsl_vector;
+
+ gsl_vector* get_gsl_vector()
+ {
+ return m_vector;
+ }
+
+ void swap_elements(size_t i, size_t j)
+ {
+ gsl_vector_swap_elements(m_vector, i, j);
+ }
+
+ void reverse()
+ {
+ gsl_vector_reverse(m_vector);
+ }
+
+ VectorImpl & scale(double x)
+ {
+ gsl_vector_scale(m_vector, x);
+ return (*this);
+ }
+
+ VectorImpl & translate(double x)
+ {
+ gsl_vector_add_constant(m_vector, x);
+ return (*this);
+ }
+
+ VectorImpl & operator+=(base_type const& _vector)
+ {
+ gsl_vector_add(m_vector, _vector.get_gsl_vector());
+ return (*this);
+ }
+
+ VectorImpl & operator-=(base_type const& _vector)
+ {
+ gsl_vector_sub(m_vector, _vector.get_gsl_vector());
+ return (*this);
+ }
+
+}; // end class VectorImpl
+
+} // end namespace detail
+
+
+using detail::operator==;
+using detail::operator<<;
+
+class Vector : public detail::VectorImpl
+{
+ public:
+ typedef detail::VectorImpl base_type;
+
+ public:
+ Vector(size_t n)
+ {
+ m_size = n;
+ m_vector = gsl_vector_alloc(n);
+ }
+
+ Vector(size_t n, double x)
+ {
+ m_size = n;
+ m_vector = gsl_vector_alloc(n);
+ gsl_vector_set_all(m_vector, x);
+ }
+
+ // create a vector with n elements all set to zero
+ // but the i-th that is set to 1
+ Vector(size_t n, size_t i)
+ {
+ m_size = n;
+ m_vector = gsl_vector_alloc(n);
+ gsl_vector_set_basis(m_vector, i);
+ }
+
+ Vector(Vector const& _vector)
+ : base_type()
+ {
+ m_size = _vector.size();
+ m_vector = gsl_vector_alloc(size());
+ gsl_vector_memcpy(m_vector, _vector.m_vector);
+ }
+
+ explicit
+ Vector(base_type::base_type const& _vector)
+ {
+ m_size = _vector.size();
+ m_vector = gsl_vector_alloc(size());
+ gsl_vector_memcpy(m_vector, _vector.get_gsl_vector());
+ }
+
+ virtual ~Vector()
+ {
+ gsl_vector_free(m_vector);
+ }
+
+
+ Vector & operator=(Vector const& _vector)
+ {
+ assert( size() == _vector.size() );
+ gsl_vector_memcpy(m_vector, _vector.m_vector);
+ return (*this);
+ }
+
+ Vector & operator=(base_type::base_type const& _vector)
+ {
+ assert( size() == _vector.size() );
+ gsl_vector_memcpy(m_vector, _vector.get_gsl_vector());
+ return (*this);
+ }
+
+ Vector & scale(double x)
+ {
+ return static_cast<Vector&>( base_type::scale(x) );
+ }
+
+ Vector & translate(double x)
+ {
+ return static_cast<Vector&>( base_type::translate(x) );
+ }
+
+ Vector & operator+=(base_type::base_type const& _vector)
+ {
+ return static_cast<Vector&>( base_type::operator+=(_vector) );
+ }
+
+ Vector & operator-=(base_type::base_type const& _vector)
+ {
+ return static_cast<Vector&>( base_type::operator-=(_vector) );
+ }
+
+ friend
+ void swap(Vector & v1, Vector & v2);
+ friend
+ void swap_any(Vector & v1, Vector & v2);
+
+}; // end class Vector
+
+
+// warning! these operations invalidate any view of the passed vector objects
+inline
+void swap(Vector & v1, Vector & v2)
+{
+ assert( v1.size() == v2.size() );
+ std::swap(v1.m_vector, v2.m_vector);
+}
+
+inline
+void swap_any(Vector & v1, Vector & v2)
+{
+ std::swap(v1.m_vector, v2.m_vector);
+ std::swap(v1.m_size, v2.m_size);
+}
+
+
+class ConstVectorView : public detail::BaseVectorImpl
+{
+ public:
+ typedef detail::BaseVectorImpl base_type;
+
+ public:
+ ConstVectorView(const base_type & _vector, size_t n, size_t offset = 0)
+ : m_vector_view( gsl_vector_const_subvector(_vector.get_gsl_vector(), offset, n) )
+ {
+ m_size = n;
+ m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );
+ }
+
+ ConstVectorView(const base_type & _vector, size_t n, size_t offset , size_t stride)
+ : m_vector_view( gsl_vector_const_subvector_with_stride(_vector.get_gsl_vector(), offset, stride, n) )
+ {
+ m_size = n;
+ m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );
+ }
+
+ ConstVectorView(const double* _vector, size_t n, size_t offset = 0)
+ : m_vector_view( gsl_vector_const_view_array(_vector + offset, n) )
+ {
+ m_size = n;
+ m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );
+ }
+
+ ConstVectorView(const double* _vector, size_t n, size_t offset, size_t stride)
+ : m_vector_view( gsl_vector_const_view_array_with_stride(_vector + offset, stride, n) )
+ {
+ m_size = n;
+ m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );
+ }
+
+ explicit
+ ConstVectorView(gsl_vector_const_view _gsl_vector_view)
+ : m_vector_view(_gsl_vector_view)
+ {
+ m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );
+ m_size = m_vector->size;
+ }
+
+ explicit
+ ConstVectorView(const std::vector<double>& _vector)
+ : m_vector_view( gsl_vector_const_view_array(&(_vector[0]), _vector.size()) )
+ {
+ m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );
+ m_size = _vector.size();
+ }
+
+ ConstVectorView(const ConstVectorView & _vector)
+ : base_type(),
+ m_vector_view(_vector.m_vector_view)
+ {
+ m_size = _vector.size();
+ m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );
+ }
+
+ ConstVectorView(const base_type & _vector)
+ : m_vector_view(gsl_vector_const_subvector(_vector.get_gsl_vector(), 0, _vector.size()))
+ {
+ m_size = _vector.size();
+ m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );
+ }
+
+ private:
+ gsl_vector_const_view m_vector_view;
+
+}; // end class ConstVectorView
+
+
+
+
+class VectorView : public detail::VectorImpl
+{
+ public:
+ typedef detail::VectorImpl base_type;
+
+ public:
+ VectorView(base_type & _vector, size_t n, size_t offset = 0, size_t stride = 1)
+ {
+ m_size = n;
+ if (stride == 1)
+ {
+ m_vector_view
+ = gsl_vector_subvector(_vector.get_gsl_vector(), offset, n);
+ m_vector = &(m_vector_view.vector);
+ }
+ else
+ {
+ m_vector_view
+ = gsl_vector_subvector_with_stride(_vector.get_gsl_vector(), offset, stride, n);
+ m_vector = &(m_vector_view.vector);
+ }
+ }
+
+ VectorView(double* _vector, size_t n, size_t offset = 0, size_t stride = 1)
+ {
+ m_size = n;
+ if (stride == 1)
+ {
+ m_vector_view
+ = gsl_vector_view_array(_vector + offset, n);
+ m_vector = &(m_vector_view.vector);
+ }
+ else
+ {
+ m_vector_view
+ = gsl_vector_view_array_with_stride(_vector + offset, stride, n);
+ m_vector = &(m_vector_view.vector);
+ }
+
+ }
+
+ VectorView(const VectorView & _vector)
+ : base_type()
+ {
+ m_size = _vector.size();
+ m_vector_view = _vector.m_vector_view;
+ m_vector = &(m_vector_view.vector);
+ }
+
+ VectorView(Vector & _vector)
+ {
+ m_size = _vector.size();
+ m_vector_view = gsl_vector_subvector(_vector.get_gsl_vector(), 0, size());
+ m_vector = &(m_vector_view.vector);
+ }
+
+ explicit
+ VectorView(gsl_vector_view _gsl_vector_view)
+ : m_vector_view(_gsl_vector_view)
+ {
+ m_vector = &(m_vector_view.vector);
+ m_size = m_vector->size;
+ }
+
+ explicit
+ VectorView(std::vector<double> & _vector)
+ {
+ m_size = _vector.size();
+ m_vector_view = gsl_vector_view_array(&(_vector[0]), _vector.size());
+ m_vector = &(m_vector_view.vector);
+ }
+
+ VectorView & operator=(VectorView const& _vector)
+ {
+ assert( size() == _vector.size() );
+ gsl_vector_memcpy(m_vector, _vector.get_gsl_vector());
+ return (*this);
+ }
+
+ VectorView & operator=(base_type::base_type const& _vector)
+ {
+ assert( size() == _vector.size() );
+ gsl_vector_memcpy(m_vector, _vector.get_gsl_vector());
+ return (*this);
+ }
+
+ VectorView & scale(double x)
+ {
+ return static_cast<VectorView&>( base_type::scale(x) );
+ }
+
+ VectorView & translate(double x)
+ {
+ return static_cast<VectorView&>( base_type::translate(x) );
+ }
+
+ VectorView & operator+=(base_type::base_type const& _vector)
+ {
+ return static_cast<VectorView&>( base_type::operator+=(_vector) );
+ }
+
+ VectorView & operator-=(base_type::base_type const& _vector)
+ {
+ return static_cast<VectorView&>( base_type::operator-=(_vector) );
+ }
+
+ friend
+ void swap_view(VectorView & v1, VectorView & v2);
+
+ private:
+ gsl_vector_view m_vector_view;
+
+}; // end class VectorView
+
+
+inline
+void swap_view(VectorView & v1, VectorView & v2)
+{
+ assert( v1.size() == v2.size() );
+ std::swap(v1.m_vector_view, v2.m_vector_view); // not swap m_vector too
+}
+
+
+} } // end namespaces
+
+
+#endif /*_NL_VECTOR_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/path-intersection.cpp b/src/2geom/path-intersection.cpp index 2c32c3764..635996e51 100644 --- a/src/2geom/path-intersection.cpp +++ b/src/2geom/path-intersection.cpp @@ -1,9 +1,9 @@ -#include "path-intersection.h" +#include <2geom/path-intersection.h> -#include "ord.h" +#include <2geom/ord.h> //for path_direction: -#include "sbasis-geometric.h" +#include <2geom/sbasis-geometric.h> namespace Geom { @@ -227,7 +227,7 @@ Crossings SimpleCrosser::crossings(Curve const &a, Curve const &b) { */ void mono_pair(Path const &A, double Al, double Ah, Path const &B, double Bl, double Bh, - Crossings &ret, double tol, unsigned depth = 0) { + Crossings &ret, double /*tol*/, unsigned depth = 0) { if( Al >= Ah || Bl >= Bh) return; std::cout << " " << depth << "[" << Al << ", " << Ah << "]" << "[" << Bl << ", " << Bh << "]"; diff --git a/src/2geom/path-intersection.h b/src/2geom/path-intersection.h index 5cdee5509..2094c4beb 100644 --- a/src/2geom/path-intersection.h +++ b/src/2geom/path-intersection.h @@ -1,11 +1,11 @@ #ifndef __GEOM_PATH_INTERSECTION_H #define __GEOM_PATH_INTERSECTION_H -#include "path.h" +#include <2geom/path.h> -#include "crossing.h" +#include <2geom/crossing.h> -#include "sweep.h" +#include <2geom/sweep.h> namespace Geom { diff --git a/src/2geom/path.cpp b/src/2geom/path.cpp index 8bfccf87a..2e5b789d5 100644 --- a/src/2geom/path.cpp +++ b/src/2geom/path.cpp @@ -34,7 +34,7 @@ -#include "path.h" +#include <2geom/path.h> namespace Geom @@ -49,7 +49,9 @@ void Path::swap(Path &other) { } Rect Path::boundsFast() const { - Rect bounds=front().boundsFast(); + Rect bounds; + if (empty()) return bounds; + bounds=front().boundsFast(); const_iterator iter = begin(); if ( iter != end() ) { for ( ++iter; iter != end() ; ++iter ) { @@ -60,7 +62,9 @@ Rect Path::boundsFast() const { } Rect Path::boundsExact() const { - Rect bounds=front().boundsExact(); + Rect bounds; + if (empty()) return bounds; + bounds=front().boundsExact(); const_iterator iter = begin(); if ( iter != end() ) { for ( ++iter; iter != end() ; ++iter ) { @@ -70,6 +74,33 @@ Rect Path::boundsExact() const { return bounds; } +/* +Rect Path::boundsFast() +{ + Rect bound; + if (empty()) return bound; + + bound = begin()->boundsFast(); + for (iterator it = ++begin(); it != end(); ++it) + { + bound.unionWith(it->boundsFast()); + } + return bound; +} + +Rect Path::boundsExact() +{ + Rect bound; + if (empty()) return bound; + + bound = begin()->boundsExact(); + for (iterator it = ++begin(); it != end(); ++it) + { + bound.unionWith(it->boundsExact()); + } + return bound; + }*/ + template<typename iter> iter inc(iter const &x, unsigned n) { iter ret = x; @@ -239,33 +270,6 @@ double Path::nearestPoint(Point const& _point, double from, double to) const return ni + nearest; } -Rect Path::boundsFast() -{ - Rect bound; - if (empty()) return bound; - - bound = begin()->boundsFast(); - for (iterator it = ++begin(); it != end(); ++it) - { - bound.unionWith(it->boundsFast()); - } - return bound; -} - -Rect Path::boundsExact() -{ - Rect bound; - if (empty()) return bound; - - bound = begin()->boundsExact(); - for (iterator it = ++begin(); it != end(); ++it) - { - bound.unionWith(it->boundsExact()); - } - return bound; -} - -//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; @@ -276,7 +280,7 @@ void Path::appendPortionTo(Path &ret, double from, double to) const { const_iterator fromi = inc(begin(), (unsigned)fi); if(fi == ti && from < to) { Curve *v = fromi->portion(ff, tf); - ret.append(*v); + ret.append(*v, STITCH_DISCONTINUOUS); delete v; return; } @@ -284,59 +288,28 @@ void Path::appendPortionTo(Path &ret, double from, double to) const { if(ff != 1.) { Curve *fromv = fromi->portion(ff, 1.); //fromv->setInitial(ret.finalPoint()); - ret.append(*fromv); + ret.append(*fromv, STITCH_DISCONTINUOUS); 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); + ret.insert(ret.end(), ++fromi, ender, STITCH_DISCONTINUOUS); + ret.insert(ret.end(), begin(), toi, STITCH_DISCONTINUOUS); } else { - ret.insert(ret.end(), ++fromi, toi); + ret.insert(ret.end(), ++fromi, toi, STITCH_DISCONTINUOUS); } Curve *tov = toi->portion(0., tf); - ret.append(*tov); + ret.append(*tov, STITCH_DISCONTINUOUS); delete tov; } -const double eps = .1; - -void Path::append(Curve const &curve) { - if ( curves_.front() != final_ && !are_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 ( !are_near(curve[i][0][0], (*final_)[0][i], eps) ) { - THROW_CONTINUITYERROR(); - } - } - } - do_append(new SBasisCurve(curve)); -} - -void Path::append(Path const &other) -{ - // Check that path stays continuous: - if ( !are_near( finalPoint(), other.initialPoint() ) ) { - THROW_CONTINUITYERROR(); - } - - insert(begin(), other.begin(), other.end()); -} - void Path::do_update(Sequence::iterator first_replaced, Sequence::iterator last_replaced, Sequence::iterator first, - Sequence::iterator last) + 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 ) ) { @@ -356,6 +329,10 @@ void Path::do_update(Sequence::iterator first_replaced, void Path::do_append(Curve *curve) { if ( curves_.front() == final_ ) { final_->setPoint(1, curve->initialPoint()); + } else { + if (curve->initialPoint() != finalPoint()) { + THROW_CONTINUITYERROR(); + } } curves_.insert(curves_.end()-1, curve); final_->setPoint(0, curve->finalPoint()); @@ -367,6 +344,28 @@ void Path::delete_range(Sequence::iterator first, Sequence::iterator last) { } } +void Path::stitch(Sequence::iterator first_replaced, + Sequence::iterator last_replaced, + Sequence &source) +{ + if (!source.empty()) { + if ( first_replaced != curves_.begin() ) { + if ( (*first_replaced)->initialPoint() != source.front()->initialPoint() ) { + source.insert(source.begin(), new StitchSegment((*first_replaced)->initialPoint(), source.front()->initialPoint())); + } + } + if ( last_replaced != (curves_.end()-1) ) { + if ( (*last_replaced)->finalPoint() != source.back()->finalPoint() ) { + source.insert(source.end(), new StitchSegment(source.back()->finalPoint(), (*last_replaced)->finalPoint())); + } + } + } else if ( first_replaced != last_replaced && first_replaced != curves_.begin() && last_replaced != curves_.end()-1) { + if ( (*first_replaced)->initialPoint() != (*(last_replaced-1))->finalPoint() ) { + source.insert(source.begin(), new StitchSegment((*(last_replaced-1))->finalPoint(), (*first_replaced)->initialPoint())); + } + } +} + void Path::check_continuity(Sequence::iterator first_replaced, Sequence::iterator last_replaced, Sequence::iterator first, @@ -374,17 +373,17 @@ void Path::check_continuity(Sequence::iterator first_replaced, { if ( first != last ) { if ( first_replaced != curves_.begin() ) { - if ( !are_near( (*first_replaced)->initialPoint(), (*first)->initialPoint(), eps ) ) { + if ( (*first_replaced)->initialPoint() != (*first)->initialPoint() ) { THROW_CONTINUITYERROR(); } } if ( last_replaced != (curves_.end()-1) ) { - if ( !are_near( (*(last_replaced-1))->finalPoint(), (*(last-1))->finalPoint(), eps ) ) { + if ( (*(last_replaced-1))->finalPoint() != (*(last-1))->finalPoint() ) { THROW_CONTINUITYERROR(); } } } else if ( first_replaced != last_replaced && first_replaced != curves_.begin() && last_replaced != curves_.end()-1) { - if ( !are_near((*first_replaced)->initialPoint(), (*(last_replaced-1))->finalPoint(), eps ) ) { + if ( (*first_replaced)->initialPoint() != (*(last_replaced-1))->finalPoint() ) { THROW_CONTINUITYERROR(); } } diff --git a/src/2geom/path.h b/src/2geom/path.h index 5f2c86b95..1c142f98d 100644 --- a/src/2geom/path.h +++ b/src/2geom/path.h @@ -38,7 +38,7 @@ #define SEEN_GEOM_PATH_H -#include "curves.h" +#include <2geom/curves.h> #include <iterator> @@ -46,6 +46,16 @@ namespace Geom { +// Conditional expression for types. If true, first, if false, second. +template<bool _Cond, typename _Iftrue, typename _Iffalse> + struct __conditional_type + { typedef _Iftrue __type; }; + +template<typename _Iftrue, typename _Iffalse> + struct __conditional_type<false, _Iftrue, _Iffalse> + { typedef _Iffalse __type; }; + + template <typename IteratorImpl> class BaseIterator : public std::iterator<std::forward_iterator_tag, Curve const> @@ -56,6 +66,20 @@ public: // default construct // default copy + // Allow Sequence::iterator to Sequence::const_iterator conversion + // unfortunately I do not know how to imitate the way __normal_iterator + // does it, because I don't see a way to get the typename of the container + // IteratorImpl is pointing at... + typedef std::vector<Curve *> Sequence; + BaseIterator ( typename __conditional_type< + (std::__are_same<IteratorImpl, Sequence::const_iterator >::__value), // check if this instantiation is of const_iterator type + const BaseIterator< Sequence::iterator >, // if true: accept iterator in const_iterator instantiation + const BaseIterator<IteratorImpl> > ::__type // if false: default to standard copy constructor + & __other) + : impl_(__other.impl_) { } + friend class BaseIterator< Sequence::const_iterator >; + + bool operator==(BaseIterator const &other) { return other.impl_ == impl_; } @@ -151,28 +175,47 @@ public: typedef Sequence::size_type size_type; typedef Sequence::difference_type difference_type; + class ClosingSegment : public LineSegment { + public: + ClosingSegment() : LineSegment() {} + ClosingSegment(Point const &p1, Point const &p2) : LineSegment(p1, p2) {} + virtual Curve *duplicate() const { return new ClosingSegment(*this); } + }; + + enum Stitching { + NO_STITCHING=0, + STITCH_DISCONTINUOUS + }; + + class StitchSegment : public LineSegment { + public: + StitchSegment() : LineSegment() {} + StitchSegment(Point const &p1, Point const &p2) : LineSegment(p1, p2) {} + virtual Curve *duplicate() const { return new StitchSegment(*this); } + }; + Path() - : final_(new LineSegment()), closed_(false) + : final_(new ClosingSegment()), closed_(false) { curves_.push_back(final_); } Path(Path const &other) - : final_(new LineSegment()), closed_(other.closed_) + : final_(new ClosingSegment()), closed_(other.closed_) { curves_.push_back(final_); insert(begin(), other.begin(), other.end()); } explicit Path(Point p) - : final_(new LineSegment(p, p)), closed_(false) + : final_(new ClosingSegment(p, p)), closed_(false) { curves_.push_back(final_); } template <typename Impl> Path(BaseIterator<Impl> first, BaseIterator<Impl> last, bool closed=false) - : closed_(closed), final_(new LineSegment()) + : closed_(closed), final_(new ClosingSegment()) { curves_.push_back(final_); insert(begin(), first, last); @@ -236,18 +279,67 @@ public: return ret; } + bool operator==(Path const &m) const { + if (size() != m.size() || closed() != m.closed()) + return false; + const_iterator it2 = m.curves_.begin(); + for(const_iterator it = curves_.begin(); it != curves_.end(); ++it) { + const Curve& a = (*it); + const Curve& b = (*it2); + if(!(a == b)) + return false; + ++it2; + } + return true; + } + + /* + Path operator*=(Matrix) + This is not possible without at least partly regenerating the curves of + the path, because a path can consist of many types of curves, + e.g. a HLineSegment. + Such a segment cannot be transformed and stay a HLineSegment in general + (take for example rotations). + This means that these curves of the path have to be replaced with + LineSegments: new Curves. + So an implementation of this method should check the curve's type to see + whether operator*= is doable for that curve type, ... + */ Path operator*(Matrix const &m) const { Path ret; + ret.curves_.reserve(curves_.size()); for(const_iterator it = begin(); it != end(); ++it) { - Curve *temp = it->transformed(m); - //Possible point of discontinuity? - ret.append(*temp); - delete temp; + Curve *curve = it->transformed(m); + ret.do_append(curve); } - ret.closed_ = closed_; + ret.close(closed_); return ret; } + /* + // this should be even quickier but it works at low level + Path operator*(Matrix const &m) const + { + Path result; + size_t sz = curves_.size() - 1; + if (sz == 0) return result; + result.curves_.resize(curves_.size()); + result.curves_.back() = result.final_; + result.curves_[0] = (curves_[0])->transformed(m); + for (size_t i = 1; i < sz; ++i) + { + result.curves_[i] = (curves_[i])->transformed(m); + if ( result.curves_[i]->initialPoint() != result.curves_[i-1]->finalPoint() ) { + THROW_CONTINUITYERROR(); + } + } + result.final_->setInitial( (result.curves_[sz])->finalPoint() ); + result.final_->setFinal( (result.curves_[0])->initialPoint() ); + result.closed_ = closed_; + return result; + } + */ + Point pointAt(double t) const { unsigned int sz = size(); @@ -323,9 +415,6 @@ public: return nearestPoint(_point, 0, sz); } - Rect boundsFast(); - Rect boundsExact(); - void appendPortionTo(Path &p, double f, double t) const; Path portion(double f, double t) const { @@ -340,17 +429,18 @@ public: 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 since append makes a copy delete temp; } return ret; } - - void insert(iterator pos, Curve const &curve) { + + void insert(iterator pos, Curve const &curve, Stitching stitching=NO_STITCHING) { Sequence source(1, curve.duplicate()); try { + if (stitching) stitch(pos.impl_, pos.impl_, source); do_update(pos.impl_, pos.impl_, source.begin(), source.end()); } catch (...) { delete_range(source.begin(), source.end()); @@ -359,11 +449,12 @@ public: } template <typename Impl> - void insert(iterator pos, BaseIterator<Impl> first, BaseIterator<Impl> last) + void insert(iterator pos, BaseIterator<Impl> first, BaseIterator<Impl> last, Stitching stitching=NO_STITCHING) { Sequence source(DuplicatingIterator<Impl>(first.impl_), DuplicatingIterator<Impl>(last.impl_)); try { + if (stitching) stitch(pos.impl_, pos.impl_, source); do_update(pos.impl_, pos.impl_, source.begin(), source.end()); } catch (...) { delete_range(source.begin(), source.end()); @@ -376,12 +467,34 @@ public: curves_.begin(), curves_.begin()); } - void erase(iterator pos) { - do_update(pos.impl_, pos.impl_+1, curves_.begin(), curves_.begin()); + void erase(iterator pos, Stitching stitching=NO_STITCHING) { + if (stitching) { + Sequence stitched; + stitch(pos.impl_, pos.impl_+1, stitched); + try { + do_update(pos.impl_, pos.impl_+1, stitched.begin(), stitched.end()); + } catch (...) { + delete_range(stitched.begin(), stitched.end()); + throw; + } + } else { + do_update(pos.impl_, pos.impl_+1, curves_.begin(), curves_.begin()); + } } - void erase(iterator first, iterator last) { - do_update(first.impl_, last.impl_, curves_.begin(), curves_.begin()); + void erase(iterator first, iterator last, Stitching stitching=NO_STITCHING) { + if (stitching) { + Sequence stitched; + stitch(first.impl_, last.impl_, stitched); + try { + do_update(first.impl_, last.impl_, stitched.begin(), stitched.end()); + } catch (...) { + delete_range(stitched.begin(), stitched.end()); + throw; + } + } else { + do_update(first.impl_, last.impl_, curves_.begin(), curves_.begin()); + } } // erase last segment of path @@ -389,9 +502,10 @@ public: erase(curves_.end()-2); } - void replace(iterator replaced, Curve const &curve) { + void replace(iterator replaced, Curve const &curve, Stitching stitching=NO_STITCHING) { Sequence source(1, curve.duplicate()); try { + if (stitching) stitch(replaced.impl_, replaced.impl_+1, source); do_update(replaced.impl_, replaced.impl_+1, source.begin(), source.end()); } catch (...) { delete_range(source.begin(), source.end()); @@ -400,10 +514,11 @@ public: } void replace(iterator first_replaced, iterator last_replaced, - Curve const &curve) + Curve const &curve, Stitching stitching=NO_STITCHING) { Sequence source(1, curve.duplicate()); try { + if (stitching) stitch(first_replaced.impl_, last_replaced.impl_, source); do_update(first_replaced.impl_, last_replaced.impl_, source.begin(), source.end()); } catch (...) { @@ -414,11 +529,13 @@ public: template <typename Impl> void replace(iterator replaced, - BaseIterator<Impl> first, BaseIterator<Impl> last) + BaseIterator<Impl> first, BaseIterator<Impl> last, + Stitching stitching=NO_STITCHING) { Sequence source(DuplicatingIterator<Impl>(first.impl_), DuplicatingIterator<Impl>(last.impl_)); try { + if (stitching) stitch(replaced.impl_, replaced.impl_+1, source); do_update(replaced.impl_, replaced.impl_+1, source.begin(), source.end()); } catch (...) { delete_range(source.begin(), source.end()); @@ -428,10 +545,12 @@ public: template <typename Impl> void replace(iterator first_replaced, iterator last_replaced, - BaseIterator<Impl> first, BaseIterator<Impl> last) + BaseIterator<Impl> first, BaseIterator<Impl> last, + Stitching stitching=NO_STITCHING) { Sequence source(first.impl_, last.impl_); try { + if (stitching) stitch(first_replaced.impl_, last_replaced.impl_, source); do_update(first_replaced.impl_, last_replaced.impl_, source.begin(), source.end()); } catch (...) { @@ -485,65 +604,83 @@ public: } } - void append(Curve const &curve); - void append(D2<SBasis> const &curve); - void append(Path const &other); + void append(Curve const &curve, Stitching stitching=NO_STITCHING) { + if (stitching) stitchTo(curve.initialPoint()); + do_append(curve.duplicate()); + } + void append(D2<SBasis> const &curve, Stitching stitching=NO_STITCHING) { + if (stitching) stitchTo(Point(curve[X][0][0], curve[Y][0][0])); + do_append(new SBasisCurve(curve)); + } + void append(Path const &other, Stitching stitching=NO_STITCHING) { + insert(end(), other.begin(), other.end(), stitching); + } + + void stitchTo(Point const &p) { + if (!empty() && finalPoint() != p) { + do_append(new StitchSegment(finalPoint(), p)); + } + } template <typename CurveType, typename A> void appendNew(A a) { - do_append(new CurveType((*final_)[0], a)); + do_append(new CurveType(finalPoint(), a)); } template <typename CurveType, typename A, typename B> void appendNew(A a, B b) { - do_append(new CurveType((*final_)[0], a, b)); + do_append(new CurveType(finalPoint(), a, b)); } template <typename CurveType, typename A, typename B, typename C> void appendNew(A a, B b, C c) { - do_append(new CurveType((*final_)[0], a, b, c)); + do_append(new CurveType(finalPoint(), a, b, c)); } template <typename CurveType, typename A, typename B, typename C, typename D> void appendNew(A a, B b, C c, D d) { - do_append(new CurveType((*final_)[0], a, b, c, d)); + do_append(new CurveType(finalPoint(), a, b, c, d)); } template <typename CurveType, typename A, typename B, typename C, typename D, typename E> void appendNew(A a, B b, C c, D d, E e) { - do_append(new CurveType((*final_)[0], a, b, c, d, e)); + do_append(new CurveType(finalPoint(), a, b, c, d, e)); } template <typename CurveType, typename A, typename B, typename C, typename D, typename E, typename F> void appendNew(A a, B b, C c, D d, E e, F f) { - do_append(new CurveType((*final_)[0], a, b, c, d, e, f)); + do_append(new CurveType(finalPoint(), a, b, c, d, e, f)); } template <typename CurveType, typename A, typename B, typename C, typename D, typename E, typename F, typename G> void appendNew(A a, B b, C c, D d, E e, F f, G g) { - do_append(new CurveType((*final_)[0], a, b, c, d, e, f, g)); + do_append(new CurveType(finalPoint(), a, b, c, d, e, f, g)); } template <typename CurveType, typename A, typename B, typename C, typename D, typename E, typename F, typename G, typename H> void appendNew(A a, B b, C c, D d, E e, F f, G g, H h) { - do_append(new CurveType((*final_)[0], a, b, c, d, e, f, g, h)); + do_append(new CurveType(finalPoint(), a, b, c, d, e, f, g, h)); } template <typename CurveType, typename A, typename B, typename C, typename D, typename E, typename F, typename G, typename H, typename I> void appendNew(A a, B b, C c, D d, E e, F f, G g, H h, I i) { - do_append(new CurveType((*final_)[0], a, b, c, d, e, f, g, h, i)); + do_append(new CurveType(finalPoint(), a, b, c, d, e, f, g, h, i)); } private: + void stitch(Sequence::iterator first_replaced, + Sequence::iterator last_replaced, + Sequence &sequence); + void do_update(Sequence::iterator first_replaced, Sequence::iterator last_replaced, Sequence::iterator first, @@ -551,7 +688,7 @@ private: void do_append(Curve *curve); - void delete_range(Sequence::iterator first, Sequence::iterator last); + static void delete_range(Sequence::iterator first, Sequence::iterator last); void check_continuity(Sequence::iterator first_replaced, Sequence::iterator last_replaced, @@ -559,7 +696,7 @@ private: Sequence::iterator last); Sequence curves_; - LineSegment *final_; + ClosingSegment *final_; bool closed_; }; // end class Path @@ -578,54 +715,6 @@ Coord nearest_point(Point const& p, Path const& c) return c.nearestPoint(p); } - -/* -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(); } - - std::vector<double> roots(double v, Dim2 d) const = 0; - - 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, 0).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; - -}; -*/ - } // end namespace Geom namespace std { diff --git a/src/2geom/pathvector.cpp b/src/2geom/pathvector.cpp index fe2f1976c..9befc1f17 100644 --- a/src/2geom/pathvector.cpp +++ b/src/2geom/pathvector.cpp @@ -35,10 +35,10 @@ #ifndef SEEN_GEOM_PATHVECTOR_CPP #define SEEN_GEOM_PATHVECTOR_CPP -#include "pathvector.h" +#include <2geom/pathvector.h> -#include "path.h" -#include "matrix.h" +#include <2geom/path.h> +#include <2geom/matrix.h> namespace Geom { diff --git a/src/2geom/pathvector.h b/src/2geom/pathvector.h index 6a759e406..17f6b6212 100644 --- a/src/2geom/pathvector.h +++ b/src/2geom/pathvector.h @@ -35,9 +35,9 @@ #ifndef SEEN_GEOM_PATHVECTOR_H #define SEEN_GEOM_PATHVECTOR_H -#include "forward.h" -#include "path.h" -#include "rect.h" +#include <2geom/forward.h> +#include <2geom/path.h> +#include <2geom/rect.h> namespace Geom { diff --git a/src/2geom/piecewise.cpp b/src/2geom/piecewise.cpp index 222152aa3..2d8638818 100644 --- a/src/2geom/piecewise.cpp +++ b/src/2geom/piecewise.cpp @@ -29,7 +29,7 @@ * */ -#include "piecewise.h" +#include <2geom/piecewise.h> #include <iterator> #include <map> @@ -167,21 +167,6 @@ std::vector<double> roots(Piecewise<SBasis> const &f){ return result; } -Piecewise<SBasis> reverse(Piecewise<SBasis> const &f) { - Piecewise<SBasis> ret = Piecewise<SBasis>(); - ret.cuts.resize(f.cuts.size()); - ret.segs.resize(f.segs.size()); - double start = f.cuts[0]; - double end = f.cuts.back(); - for (unsigned i = 0; i < f.cuts.size(); i++) { - double x = f.cuts[f.cuts.size() - 1 - i]; - ret.cuts[i] = end - (x - start); - } - for (unsigned i = 0; i < f.segs.size(); i++) - ret.segs[i] = reverse(f[f.segs.size() - i - 1]); - return ret; -} - } /* diff --git a/src/2geom/piecewise.h b/src/2geom/piecewise.h index 96e64bcf9..73fc76e2f 100644 --- a/src/2geom/piecewise.h +++ b/src/2geom/piecewise.h @@ -31,12 +31,12 @@ #ifndef SEEN_GEOM_PW_SB_H #define SEEN_GEOM_PW_SB_H -#include "sbasis.h" +#include <2geom/sbasis.h> #include <vector> #include <map> -#include "concepts.h" -#include "isnan.h" +#include <2geom/concepts.h> +#include <2geom/isnan.h> #include <boost/concept_check.hpp> namespace Geom { diff --git a/src/2geom/point-l.h b/src/2geom/point-l.h index cd4768db0..88c4af948 100644 --- a/src/2geom/point-l.h +++ b/src/2geom/point-l.h @@ -2,7 +2,7 @@ #define SEEN_Geom_POINT_L_H #include <stdexcept> -#include "point.h" +#include <2geom/point.h> namespace Geom { diff --git a/src/2geom/point.cpp b/src/2geom/point.cpp index 66eaf8ca6..42a7ecef3 100644 --- a/src/2geom/point.cpp +++ b/src/2geom/point.cpp @@ -1,8 +1,8 @@ -#include "point.h" +#include <2geom/point.h> #include <assert.h> -#include "coord.h" -#include "isnan.h" //temporary fix for isnan() -#include "matrix.h" +#include <2geom/coord.h> +#include <2geom/isnan.h> //temporary fix for isnan() +#include <2geom/matrix.h> namespace Geom { diff --git a/src/2geom/point.h b/src/2geom/point.h index 91df159ca..2cab3d7fe 100644 --- a/src/2geom/point.h +++ b/src/2geom/point.h @@ -7,8 +7,8 @@ #include <iostream> -#include "coord.h" -#include "utils.h" +#include <2geom/coord.h> +#include <2geom/utils.h> namespace Geom { @@ -91,7 +91,7 @@ class Point { _pt[i] += o._pt[i]; } return *this; - } + } inline Point &operator-=(Point const &o) { for ( unsigned i = 0 ; i < 2 ; ++i ) { _pt[i] -= o._pt[i]; @@ -179,6 +179,12 @@ inline bool are_near(Point const &a, Point const &b, double const eps=EPSILON) { return ( are_near(a[X],b[X],eps) && are_near(a[Y],b[Y],eps) ); } +inline +Point middle_point(Point const& P1, Point const& P2) +{ + return (P1 + P2) / 2; +} + /** Returns p * Geom::rotate_degrees(90), but more efficient. * * Angle direction in Inkscape code: If you use the traditional mathematics convention that y diff --git a/src/2geom/poly-dk-solve.cpp b/src/2geom/poly-dk-solve.cpp index fdc1cefe5..2787c4478 100644 --- a/src/2geom/poly-dk-solve.cpp +++ b/src/2geom/poly-dk-solve.cpp @@ -1,4 +1,4 @@ -#include "poly-dk-solve.h" +#include <2geom/poly-dk-solve.h> #include <iterator> /*** implementation of the Durand-Kerner method. seems buggy*/ diff --git a/src/2geom/poly-dk-solve.h b/src/2geom/poly-dk-solve.h index f82caf394..dadd12110 100644 --- a/src/2geom/poly-dk-solve.h +++ b/src/2geom/poly-dk-solve.h @@ -1,4 +1,4 @@ -#include "poly.h" +#include <2geom/poly.h> #include <complex> std::vector<std::complex<double> > diff --git a/src/2geom/poly-laguerre-solve.cpp b/src/2geom/poly-laguerre-solve.cpp index 266e24c2b..bd2ddf305 100644 --- a/src/2geom/poly-laguerre-solve.cpp +++ b/src/2geom/poly-laguerre-solve.cpp @@ -1,4 +1,4 @@ -#include "poly-laguerre-solve.h" +#include <2geom/poly-laguerre-solve.h> #include <iterator> typedef std::complex<double> cdouble; @@ -50,8 +50,8 @@ cdouble laguerre_internal_complex(Poly const & p, //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]; + //if(shuffle_counter % shuffle_rate == 0) + //a *= shuffle[shuffle_counter / shuffle_rate]; xk -= a; shuffle_counter++; if(shuffle_counter >= 90) @@ -130,9 +130,9 @@ laguerre(Poly p, const double tol) { } std::vector<double> -laguerre_real_interval(Poly const & ply, - const double lo, const double hi, - const double tol) +laguerre_real_interval(Poly const & /*ply*/, + const double /*lo*/, const double /*hi*/, + const double /*tol*/) { /* not implemented*/ assert(false); diff --git a/src/2geom/poly-laguerre-solve.h b/src/2geom/poly-laguerre-solve.h index c7bfef245..6f7cb0aee 100644 --- a/src/2geom/poly-laguerre-solve.h +++ b/src/2geom/poly-laguerre-solve.h @@ -1,4 +1,4 @@ -#include "poly.h" +#include <2geom/poly.h> #include <complex> std::vector<std::complex<double> > diff --git a/src/2geom/poly.cpp b/src/2geom/poly.cpp index 27f98596b..5f7663957 100644 --- a/src/2geom/poly.cpp +++ b/src/2geom/poly.cpp @@ -1,4 +1,4 @@ -#include "poly.h" +#include <2geom/poly.h> Poly Poly::operator*(const Poly& p) const { Poly result; @@ -167,7 +167,7 @@ Poly divide(Poly const &a, Poly const &b, Poly &r) { return c; } -Poly gcd(Poly const &a, Poly const &b, const double tol) { +Poly gcd(Poly const &a, Poly const &b, const double /*tol*/) { if(a.size() < b.size()) return gcd(b, a); if(b.size() <= 0) diff --git a/src/2geom/poly.h b/src/2geom/poly.h index f653e80a6..2af24d25f 100644 --- a/src/2geom/poly.h +++ b/src/2geom/poly.h @@ -5,7 +5,7 @@ #include <iostream> #include <algorithm> #include <complex> -#include "utils.h" +#include <2geom/utils.h> class Poly : public std::vector<double>{ public: diff --git a/src/2geom/quadtree.cpp b/src/2geom/quadtree.cpp index bb041edbe..e322a091b 100644 --- a/src/2geom/quadtree.cpp +++ b/src/2geom/quadtree.cpp @@ -1,4 +1,4 @@ -#include "quadtree.h" +#include <2geom/quadtree.h> namespace Geom{ Quad* QuadTree::search(Rect const &r) { diff --git a/src/2geom/quadtree.h b/src/2geom/quadtree.h index 716c56673..0770cd5f8 100644 --- a/src/2geom/quadtree.h +++ b/src/2geom/quadtree.h @@ -1,7 +1,7 @@ #include <vector> #include <cassert> -#include "d2.h" +#include <2geom/d2.h> namespace Geom{ diff --git a/src/2geom/rect.h b/src/2geom/rect.h index c89946606..27a475cf3 100644 --- a/src/2geom/rect.h +++ b/src/2geom/rect.h @@ -41,7 +41,7 @@ #ifndef _2GEOM_RECT #define _2GEOM_RECT -#include "matrix.h" +#include <2geom/matrix.h> #include <boost/optional/optional.hpp> namespace Geom { diff --git a/src/2geom/region.cpp b/src/2geom/region.cpp index 1bc4a68e7..065f3f418 100644 --- a/src/2geom/region.cpp +++ b/src/2geom/region.cpp @@ -1,7 +1,7 @@ -#include "region.h" -#include "utils.h" +#include <2geom/region.h> +#include <2geom/utils.h> -#include "shape.h" +#include <2geom/shape.h> namespace Geom { diff --git a/src/2geom/region.h b/src/2geom/region.h index 561f3184d..960a570a4 100644 --- a/src/2geom/region.h +++ b/src/2geom/region.h @@ -1,8 +1,8 @@ #ifndef __2GEOM_REGION_H #define __2GEOM_REGION_H -#include "path.h" -#include "path-intersection.h" +#include <2geom/path.h> +#include <2geom/path-intersection.h> namespace Geom { diff --git a/src/2geom/sbasis-2d.cpp b/src/2geom/sbasis-2d.cpp index a4decd704..79e99519a 100644 --- a/src/2geom/sbasis-2d.cpp +++ b/src/2geom/sbasis-2d.cpp @@ -1,4 +1,4 @@ -#include "sbasis-2d.h" +#include <2geom/sbasis-2d.h> namespace Geom{ diff --git a/src/2geom/sbasis-2d.h b/src/2geom/sbasis-2d.h index 921a740b9..1845b5c99 100644 --- a/src/2geom/sbasis-2d.h +++ b/src/2geom/sbasis-2d.h @@ -3,8 +3,8 @@ #include <vector> #include <cassert> #include <algorithm> -#include "d2.h" -#include "sbasis.h" +#include <2geom/d2.h> +#include <2geom/sbasis.h> #include <iostream> namespace Geom{ diff --git a/src/2geom/sbasis-curve.h b/src/2geom/sbasis-curve.h index 1612f46f7..6b3888594 100644 --- a/src/2geom/sbasis-curve.h +++ b/src/2geom/sbasis-curve.h @@ -38,7 +38,7 @@ #define _2GEOM_SBASIS_CURVE_H_ -#include "curve.h" +#include <2geom/curve.h> namespace Geom diff --git a/src/2geom/sbasis-geometric.cpp b/src/2geom/sbasis-geometric.cpp index 655830389..887ca9995 100644 --- a/src/2geom/sbasis-geometric.cpp +++ b/src/2geom/sbasis-geometric.cpp @@ -1,8 +1,8 @@ -#include "sbasis-geometric.h" -#include "sbasis.h" -#include "sbasis-math.h" -//#include "solver.h" -#include "sbasis-geometric.h" +#include <2geom/sbasis-geometric.h> +#include <2geom/sbasis.h> +#include <2geom/sbasis-math.h> +//#include <2geom/solver.h> +#include <2geom/sbasis-geometric.h> /** Geometric operators on D2<SBasis> (1D->2D). * Copyright 2007 JF Barraud diff --git a/src/2geom/sbasis-geometric.h b/src/2geom/sbasis-geometric.h index 6e21a6395..7e067d801 100644 --- a/src/2geom/sbasis-geometric.h +++ b/src/2geom/sbasis-geometric.h @@ -1,7 +1,7 @@ #ifndef _SBASIS_GEOMETRIC #define _SBASIS_GEOMETRIC -#include "d2.h" -#include "piecewise.h" +#include <2geom/d2.h> +#include <2geom/piecewise.h> #include <vector> /** two-dimensional geometric operators. diff --git a/src/2geom/sbasis-math.cpp b/src/2geom/sbasis-math.cpp index b04e15bcc..f5a8ab7a1 100644 --- a/src/2geom/sbasis-math.cpp +++ b/src/2geom/sbasis-math.cpp @@ -34,7 +34,7 @@ //TODO: define a truncated compose(sb,sb, order) and extend it to pw<sb>. //TODO: in all these functions, compute 'order' according to 'tol'. -#include "sbasis-math.h" +#include <2geom/sbasis-math.h> //#define ZERO 1e-3 diff --git a/src/2geom/sbasis-math.h b/src/2geom/sbasis-math.h index 72428bc96..2a51bffb3 100644 --- a/src/2geom/sbasis-math.h +++ b/src/2geom/sbasis-math.h @@ -39,8 +39,8 @@ #define SEEN_GEOM_SB_CALCULS_H -#include "sbasis.h" -#include "piecewise.h" +#include <2geom/sbasis.h> +#include <2geom/piecewise.h> namespace Geom{ //-|x|--------------------------------------------------------------- diff --git a/src/2geom/sbasis-poly.cpp b/src/2geom/sbasis-poly.cpp index ff797920f..3fee0afe4 100644 --- a/src/2geom/sbasis-poly.cpp +++ b/src/2geom/sbasis-poly.cpp @@ -1,4 +1,4 @@ -#include "sbasis-poly.h" +#include <2geom/sbasis-poly.h> namespace Geom{ diff --git a/src/2geom/sbasis-poly.h b/src/2geom/sbasis-poly.h index 4a8aa1cd8..09c8e5487 100644 --- a/src/2geom/sbasis-poly.h +++ b/src/2geom/sbasis-poly.h @@ -1,8 +1,8 @@ #ifndef _SBASIS_TO_POLY #define _SBASIS_TO_POLY -#include "poly.h" -#include "sbasis.h" +#include <2geom/poly.h> +#include <2geom/sbasis.h> /*** Conversion between SBasis and Poly. Not recommended for general * use due to instability. diff --git a/src/2geom/sbasis-roots.cpp b/src/2geom/sbasis-roots.cpp index 51a4ec97d..30f5fb22f 100644 --- a/src/2geom/sbasis-roots.cpp +++ b/src/2geom/sbasis-roots.cpp @@ -44,9 +44,9 @@ also allow you to find intersections of multiple curves but require solving n*m #include <cmath> #include <map> -#include "sbasis.h" -#include "sbasis-to-bezier.h" -#include "solver.h" +#include <2geom/sbasis.h> +#include <2geom/sbasis-to-bezier.h> +#include <2geom/solver.h> using namespace std; diff --git a/src/2geom/sbasis-to-bezier.cpp b/src/2geom/sbasis-to-bezier.cpp index 9843b78d7..cca30291d 100644 --- a/src/2geom/sbasis-to-bezier.cpp +++ b/src/2geom/sbasis-to-bezier.cpp @@ -9,11 +9,11 @@ This is wrong, it should read W_{q,q} = 1 (n even) */ -#include "sbasis-to-bezier.h" -#include "choose.h" -#include "svg-path.h" +#include <2geom/sbasis-to-bezier.h> +#include <2geom/choose.h> +#include <2geom/svg-path.h> #include <iostream> -#include "exception.h" +#include <2geom/exception.h> namespace Geom{ diff --git a/src/2geom/sbasis-to-bezier.h b/src/2geom/sbasis-to-bezier.h index b515ebf8b..046e82996 100644 --- a/src/2geom/sbasis-to-bezier.h +++ b/src/2geom/sbasis-to-bezier.h @@ -1,8 +1,8 @@ #ifndef _SBASIS_TO_BEZIER #define _SBASIS_TO_BEZIER -#include "d2.h" -#include "path.h" +#include <2geom/d2.h> +#include <2geom/path.h> namespace Geom{ // this produces a degree k bezier from a degree k sbasis diff --git a/src/2geom/sbasis.cpp b/src/2geom/sbasis.cpp index d7a3972e1..377238d92 100644 --- a/src/2geom/sbasis.cpp +++ b/src/2geom/sbasis.cpp @@ -33,8 +33,8 @@ #include <cmath> -#include "sbasis.h" -#include "isnan.h" +#include <2geom/sbasis.h> +#include <2geom/isnan.h> namespace Geom{ @@ -258,6 +258,8 @@ SBasis integral(SBasis const &c) { SBasis derivative(SBasis const &a) { SBasis c; c.resize(a.size(), Linear(0,0)); + if(a.isZero()) + return c; for(unsigned k = 0; k < a.size()-1; k++) { double d = (2*k+1)*(a[k][1] - a[k][0]); @@ -278,6 +280,7 @@ SBasis derivative(SBasis const &a) { } void SBasis::derive() { // in place version + if(isZero()) return; for(unsigned k = 0; k < size()-1; k++) { double d = (2*k+1)*((*this)[k][1] - (*this)[k][0]); diff --git a/src/2geom/sbasis.h b/src/2geom/sbasis.h index d6f55c8d8..a9939f92a 100644 --- a/src/2geom/sbasis.h +++ b/src/2geom/sbasis.h @@ -37,10 +37,10 @@ #include <cassert> #include <iostream> -#include "linear.h" -#include "interval.h" -#include "utils.h" -#include "exception.h" +#include <2geom/linear.h> +#include <2geom/interval.h> +#include <2geom/utils.h> +#include <2geom/exception.h> namespace Geom { diff --git a/src/2geom/shape.cpp b/src/2geom/shape.cpp index a7a40066e..936dfaf70 100644 --- a/src/2geom/shape.cpp +++ b/src/2geom/shape.cpp @@ -1,7 +1,7 @@ -#include "shape.h" -#include "utils.h" -#include "sweep.h" -#include "ord.h" +#include <2geom/shape.h> +#include <2geom/utils.h> +#include <2geom/sweep.h> +#include <2geom/ord.h> #include <iostream> #include <algorithm> diff --git a/src/2geom/shape.h b/src/2geom/shape.h index a55fe9ad3..147aa9a50 100644 --- a/src/2geom/shape.h +++ b/src/2geom/shape.h @@ -4,7 +4,7 @@ #include <vector> #include <set> -#include "region.h" +#include <2geom/region.h> //TODO: BBOX optimizations diff --git a/src/2geom/solve-bezier-one-d.cpp b/src/2geom/solve-bezier-one-d.cpp index 3ec738df1..b922e970e 100644 --- a/src/2geom/solve-bezier-one-d.cpp +++ b/src/2geom/solve-bezier-one-d.cpp @@ -1,5 +1,5 @@ -#include "solver.h" -#include "point.h" +#include <2geom/solver.h> +#include <2geom/point.h> #include <algorithm> /*** Find the zeros of the bernstein function. The code subdivides until it is happy with the diff --git a/src/2geom/solve-bezier-parametric.cpp b/src/2geom/solve-bezier-parametric.cpp index fc49dcdb2..ad017f596 100644 --- a/src/2geom/solve-bezier-parametric.cpp +++ b/src/2geom/solve-bezier-parametric.cpp @@ -1,5 +1,5 @@ -#include "solver.h" -#include "point.h" +#include <2geom/solver.h> +#include <2geom/point.h> #include <algorithm> namespace Geom{ diff --git a/src/2geom/solver.h b/src/2geom/solver.h index 81ab431a2..87365ab33 100644 --- a/src/2geom/solver.h +++ b/src/2geom/solver.h @@ -1,7 +1,7 @@ #ifndef _SOLVE_SBASIS_H #define _SOLVE_SBASIS_H -#include "point.h" -#include "sbasis.h" +#include <2geom/point.h> +#include <2geom/sbasis.h> namespace Geom{ diff --git a/src/2geom/sturm.h b/src/2geom/sturm.h index 9937a6f15..097a5120a 100644 --- a/src/2geom/sturm.h +++ b/src/2geom/sturm.h @@ -1,8 +1,8 @@ #ifndef LIB2GEOM_STURM_HEADER #define LIB2GEOM_STURM_HEADER -#include "poly.h" -#include "utils.h" +#include <2geom/poly.h> +#include <2geom/utils.h> namespace Geom { diff --git a/src/2geom/svg-elliptical-arc.cpp b/src/2geom/svg-elliptical-arc.cpp new file mode 100644 index 000000000..86b919fd2 --- /dev/null +++ b/src/2geom/svg-elliptical-arc.cpp @@ -0,0 +1,1124 @@ +/*
+ * SVG Elliptical Arc Class
+ *
+ * Copyright 2008 Marco Cecchetti <mrcekets at 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, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+
+#include <2geom/svg-elliptical-arc.h>
+#include <2geom/ellipse.h>
+#include <2geom/sbasis-geometric.h>
+#include <2geom/bezier-curve.h>
+#include <2geom/poly.h>
+
+#include <cfloat>
+#include <limits>
+
+#include <2geom/numeric/vector.h>
+#include <2geom/numeric/fitting-tool.h>
+#include <2geom/numeric/fitting-model.h>
+
+
+
+namespace Geom
+{
+
+
+Rect SVGEllipticalArc::boundsExact() const
+{
+ if (isDegenerate() && is_svg_compliant())
+ return chord().boundsExact();
+
+ std::vector<double> extremes(4);
+ double cosrot = std::cos(rotation_angle());
+ double sinrot = std::sin(rotation_angle());
+ extremes[0] = std::atan2( -ray(Y) * sinrot, ray(X) * cosrot );
+ extremes[1] = extremes[0] + M_PI;
+ if ( extremes[0] < 0 ) extremes[0] += 2*M_PI;
+ extremes[2] = std::atan2( ray(Y) * cosrot, ray(X) * sinrot );
+ extremes[3] = extremes[2] + M_PI;
+ if ( extremes[2] < 0 ) extremes[2] += 2*M_PI;
+
+
+ std::vector<double>arc_extremes(4);
+ arc_extremes[0] = initialPoint()[X];
+ arc_extremes[1] = finalPoint()[X];
+ if ( arc_extremes[0] < arc_extremes[1] )
+ std::swap(arc_extremes[0], arc_extremes[1]);
+ arc_extremes[2] = initialPoint()[Y];
+ arc_extremes[3] = finalPoint()[Y];
+ if ( arc_extremes[2] < arc_extremes[3] )
+ std::swap(arc_extremes[2], arc_extremes[3]);
+
+
+ if ( start_angle() < end_angle() )
+ {
+ if ( sweep_flag() )
+ {
+ for ( unsigned int i = 0; i < extremes.size(); ++i )
+ {
+ if ( start_angle() < extremes[i] && extremes[i] < end_angle() )
+ {
+ arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];
+ }
+ }
+ }
+ else
+ {
+ for ( unsigned int i = 0; i < extremes.size(); ++i )
+ {
+ if ( start_angle() > extremes[i] || extremes[i] > end_angle() )
+ {
+ arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];
+ }
+ }
+ }
+ }
+ else
+ {
+ if ( sweep_flag() )
+ {
+ for ( unsigned int i = 0; i < extremes.size(); ++i )
+ {
+ if ( start_angle() < extremes[i] || extremes[i] < end_angle() )
+ {
+ arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];
+ }
+ }
+ }
+ else
+ {
+ for ( unsigned int i = 0; i < extremes.size(); ++i )
+ {
+ if ( start_angle() > extremes[i] && extremes[i] > end_angle() )
+ {
+ arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];
+ }
+ }
+ }
+ }
+
+ return Rect( Point(arc_extremes[1], arc_extremes[3]) ,
+ Point(arc_extremes[0], arc_extremes[2]) );
+}
+
+
+double SVGEllipticalArc::valueAtAngle(Coord t, Dim2 d) const
+{
+ double sin_rot_angle = std::sin(rotation_angle());
+ double cos_rot_angle = std::cos(rotation_angle());
+ if ( d == X )
+ {
+ return ray(X) * cos_rot_angle * std::cos(t)
+ - ray(Y) * sin_rot_angle * std::sin(t)
+ + center(X);
+ }
+ else if ( d == Y )
+ {
+ return ray(X) * sin_rot_angle * std::cos(t)
+ + ray(Y) * cos_rot_angle * std::sin(t)
+ + center(Y);
+ }
+ THROW_RANGEERROR("dimension parameter out of range");
+}
+
+
+std::vector<double>
+SVGEllipticalArc::roots(double v, Dim2 d) const
+{
+ if ( d > Y )
+ {
+ THROW_RANGEERROR("dimention out of range");
+ }
+
+ std::vector<double> sol;
+
+ if (isDegenerate() && is_svg_compliant())
+ {
+ return chord().roots(v, d);
+ }
+ else
+ {
+ if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )
+ {
+ if ( center(d) == v )
+ sol.push_back(0);
+ return sol;
+ }
+
+ const char* msg[2][2] =
+ {
+ { "d == X; ray(X) == 0; "
+ "s = (v - center(X)) / ( -ray(Y) * std::sin(rotation_angle()) ); "
+ "s should be contained in [-1,1]",
+ "d == X; ray(Y) == 0; "
+ "s = (v - center(X)) / ( ray(X) * std::cos(rotation_angle()) ); "
+ "s should be contained in [-1,1]"
+ },
+ { "d == Y; ray(X) == 0; "
+ "s = (v - center(X)) / ( ray(Y) * std::cos(rotation_angle()) ); "
+ "s should be contained in [-1,1]",
+ "d == Y; ray(Y) == 0; "
+ "s = (v - center(X)) / ( ray(X) * std::sin(rotation_angle()) ); "
+ "s should be contained in [-1,1]"
+ },
+ };
+
+ for ( unsigned int dim = 0; dim < 2; ++dim )
+ {
+ if ( are_near(ray(dim), 0) )
+ {
+ if ( initialPoint()[d] == v && finalPoint()[d] == v )
+ {
+ THROW_INFINITESOLUTIONS(0);
+ }
+ if ( (initialPoint()[d] < finalPoint()[d])
+ && (initialPoint()[d] > v || finalPoint()[d] < v) )
+ {
+ return sol;
+ }
+ if ( (initialPoint()[d] > finalPoint()[d])
+ && (finalPoint()[d] > v || initialPoint()[d] < v) )
+ {
+ return sol;
+ }
+ double ray_prj;
+ switch(d)
+ {
+ case X:
+ switch(dim)
+ {
+ case X: ray_prj = -ray(Y) * std::sin(rotation_angle());
+ break;
+ case Y: ray_prj = ray(X) * std::cos(rotation_angle());
+ break;
+ }
+ break;
+ case Y:
+ switch(dim)
+ {
+ case X: ray_prj = ray(Y) * std::cos(rotation_angle());
+ break;
+ case Y: ray_prj = ray(X) * std::sin(rotation_angle());
+ break;
+ }
+ break;
+ }
+
+ double s = (v - center(d)) / ray_prj;
+ if ( s < -1 || s > 1 )
+ {
+ THROW_LOGICALERROR(msg[d][dim]);
+ }
+ switch(dim)
+ {
+ case X:
+ s = std::asin(s); // return a value in [-PI/2,PI/2]
+ if ( logical_xor( sweep_flag(), are_near(start_angle(), M_PI/2) ) )
+ {
+ if ( s < 0 ) s += 2*M_PI;
+ }
+ else
+ {
+ s = M_PI - s;
+ if (!(s < 2*M_PI) ) s -= 2*M_PI;
+ }
+ break;
+ case Y:
+ s = std::acos(s); // return a value in [0,PI]
+ if ( logical_xor( sweep_flag(), are_near(start_angle(), 0) ) )
+ {
+ s = 2*M_PI - s;
+ if ( !(s < 2*M_PI) ) s -= 2*M_PI;
+ }
+ break;
+ }
+
+ //std::cerr << "s = " << rad_to_deg(s);
+ s = map_to_01(s);
+ //std::cerr << " -> t: " << s << std::endl;
+ if ( !(s < 0 || s > 1) )
+ sol.push_back(s);
+ return sol;
+ }
+ }
+
+ }
+
+ double rotx, roty;
+ switch(d)
+ {
+ case X:
+ rotx = std::cos(rotation_angle());
+ roty = -std::sin(rotation_angle());
+ break;
+ case Y:
+ rotx = std::sin(rotation_angle());
+ roty = std::cos(rotation_angle());
+ break;
+ }
+ double rxrotx = ray(X) * rotx;
+ double c_v = center(d) - v;
+
+ double a = -rxrotx + c_v;
+ double b = ray(Y) * roty;
+ double c = rxrotx + c_v;
+ //std::cerr << "a = " << a << std::endl;
+ //std::cerr << "b = " << b << std::endl;
+ //std::cerr << "c = " << c << std::endl;
+
+ if ( are_near(a,0) )
+ {
+ sol.push_back(M_PI);
+ if ( !are_near(b,0) )
+ {
+ double s = 2 * std::atan(-c/(2*b));
+ if ( s < 0 ) s += 2*M_PI;
+ sol.push_back(s);
+ }
+ }
+ else
+ {
+ double delta = b * b - a * c;
+ //std::cerr << "delta = " << delta << std::endl;
+ if ( are_near(delta, 0) )
+ {
+ double s = 2 * std::atan(-b/a);
+ if ( s < 0 ) s += 2*M_PI;
+ sol.push_back(s);
+ }
+ else if ( delta > 0 )
+ {
+ double sq = std::sqrt(delta);
+ double s = 2 * std::atan( (-b - sq) / a );
+ if ( s < 0 ) s += 2*M_PI;
+ sol.push_back(s);
+ s = 2 * std::atan( (-b + sq) / a );
+ if ( s < 0 ) s += 2*M_PI;
+ sol.push_back(s);
+ }
+ }
+
+ std::vector<double> arc_sol;
+ for (unsigned int i = 0; i < sol.size(); ++i )
+ {
+ //std::cerr << "s = " << rad_to_deg(sol[i]);
+ sol[i] = map_to_01(sol[i]);
+ //std::cerr << " -> t: " << sol[i] << std::endl;
+ if ( !(sol[i] < 0 || sol[i] > 1) )
+ arc_sol.push_back(sol[i]);
+ }
+ return arc_sol;
+}
+
+
+// D(E(t,C),t) = E(t+PI/2,O)
+Curve* SVGEllipticalArc::derivative() const
+{
+ if (isDegenerate() && is_svg_compliant())
+ return chord().derivative();
+
+ SVGEllipticalArc* result = new SVGEllipticalArc(*this);
+ result->m_center[X] = result->m_center[Y] = 0;
+ result->m_start_angle += M_PI/2;
+ if( !( result->m_start_angle < 2*M_PI ) )
+ {
+ result->m_start_angle -= 2*M_PI;
+ }
+ result->m_end_angle += M_PI/2;
+ if( !( result->m_end_angle < 2*M_PI ) )
+ {
+ result->m_end_angle -= 2*M_PI;
+ }
+ result->m_initial_point = result->pointAtAngle( result->start_angle() );
+ result->m_final_point = result->pointAtAngle( result->end_angle() );
+ return result;
+}
+
+
+std::vector<Point>
+SVGEllipticalArc::pointAndDerivatives(Coord t, unsigned int n) const
+{
+ if (isDegenerate() && is_svg_compliant())
+ return chord().pointAndDerivatives(t, n);
+
+ unsigned int nn = n+1; // nn represents the size of the result vector.
+ std::vector<Point> result;
+ result.reserve(nn);
+ double angle = map_unit_interval_on_circular_arc(t, start_angle(),
+ end_angle(), sweep_flag());
+ SVGEllipticalArc ea(*this);
+ ea.m_center = Point(0,0);
+ unsigned int m = std::min(nn, 4u);
+ for ( unsigned int i = 0; i < m; ++i )
+ {
+ result.push_back( ea.pointAtAngle(angle) );
+ angle += M_PI/2;
+ if ( !(angle < 2*M_PI) ) angle -= 2*M_PI;
+ }
+ m = nn / 4;
+ for ( unsigned int i = 1; i < m; ++i )
+ {
+ for ( unsigned int j = 0; j < 4; ++j )
+ result.push_back( result[j] );
+ }
+ m = nn - 4 * m;
+ for ( unsigned int i = 0; i < m; ++i )
+ {
+ result.push_back( result[i] );
+ }
+ if ( !result.empty() ) // nn != 0
+ result[0] = pointAtAngle(angle);
+ return result;
+}
+
+bool SVGEllipticalArc::containsAngle(Coord angle) const
+{
+ if ( sweep_flag() )
+ if ( start_angle() < end_angle() )
+ return ( !( angle < start_angle() || angle > end_angle() ) );
+ else
+ return ( !( angle < start_angle() && angle > end_angle() ) );
+ else
+ if ( start_angle() > end_angle() )
+ return ( !( angle > start_angle() || angle < end_angle() ) );
+ else
+ return ( !( angle > start_angle() && angle < end_angle() ) );
+}
+
+Curve* SVGEllipticalArc::portion(double f, double t) const
+{
+ if (f < 0) f = 0;
+ if (f > 1) f = 1;
+ if (t < 0) t = 0;
+ if (t > 1) t = 1;
+ if ( are_near(f, t) )
+ {
+ SVGEllipticalArc* arc = new SVGEllipticalArc();
+ arc->m_center = arc->m_initial_point = arc->m_final_point = pointAt(f);
+ arc->m_start_angle = arc->m_end_angle = m_start_angle;
+ arc->m_rot_angle = m_rot_angle;
+ arc->m_sweep = m_sweep;
+ arc->m_large_arc = m_large_arc;
+ }
+ SVGEllipticalArc* arc = new SVGEllipticalArc( *this );
+ arc->m_initial_point = pointAt(f);
+ arc->m_final_point = pointAt(t);
+ double sa = sweep_flag() ? sweep_angle() : -sweep_angle();
+ arc->m_start_angle = m_start_angle + sa * f;
+ if ( !(arc->m_start_angle < 2*M_PI) )
+ arc->m_start_angle -= 2*M_PI;
+ if ( arc->m_start_angle < 0 )
+ arc->m_start_angle += 2*M_PI;
+ arc->m_end_angle = m_start_angle + sa * t;
+ if ( !(arc->m_end_angle < 2*M_PI) )
+ arc->m_end_angle -= 2*M_PI;
+ if ( arc->m_end_angle < 0 )
+ arc->m_end_angle += 2*M_PI;
+ if ( f > t ) arc->m_sweep = !sweep_flag();
+ if ( large_arc_flag() && (arc->sweep_angle() < M_PI) )
+ arc->m_large_arc = false;
+ return arc;
+}
+
+
+std::vector<double> SVGEllipticalArc::
+allNearestPoints( Point const& p, double from, double to ) const
+{
+ std::vector<double> result;
+ if (isDegenerate() && is_svg_compliant())
+ {
+ result.push_back( chord().nearestPoint(p, from, to) );
+ return result;
+ }
+
+ if ( from > to ) std::swap(from, to);
+ if ( from < 0 || to > 1 )
+ {
+ THROW_RANGEERROR("[from,to] interval out of range");
+ }
+
+ if ( ( are_near(ray(X), 0) && are_near(ray(Y), 0) ) || are_near(from, to) )
+ {
+ result.push_back(from);
+ return result;
+ }
+ else if ( are_near(ray(X), 0) || are_near(ray(Y), 0) )
+ {
+ LineSegment seg(pointAt(from), pointAt(to));
+ Point np = seg.pointAt( seg.nearestPoint(p) );
+ if ( are_near(ray(Y), 0) )
+ {
+ if ( are_near(rotation_angle(), M_PI/2)
+ || are_near(rotation_angle(), 3*M_PI/2) )
+ {
+ result = roots(np[Y], Y);
+ }
+ else
+ {
+ result = roots(np[X], X);
+ }
+ }
+ else
+ {
+ if ( are_near(rotation_angle(), M_PI/2)
+ || are_near(rotation_angle(), 3*M_PI/2) )
+ {
+ result = roots(np[X], X);
+ }
+ else
+ {
+ result = roots(np[Y], Y);
+ }
+ }
+ return result;
+ }
+ else if ( are_near(ray(X), ray(Y)) )
+ {
+ Point r = p - center();
+ if ( are_near(r, Point(0,0)) )
+ {
+ THROW_INFINITESOLUTIONS(0);
+ }
+ // TODO: implement case r != 0
+// Point np = ray(X) * unit_vector(r);
+// std::vector<double> solX = roots(np[X],X);
+// std::vector<double> solY = roots(np[Y],Y);
+// double t;
+// if ( are_near(solX[0], solY[0]) || are_near(solX[0], solY[1]))
+// {
+// t = solX[0];
+// }
+// else
+// {
+// t = solX[1];
+// }
+// if ( !(t < from || t > to) )
+// {
+// result.push_back(t);
+// }
+// else
+// {
+//
+// }
+ }
+
+ // solve the equation <D(E(t),t)|E(t)-p> == 0
+ // that provides min and max distance points
+ // on the ellipse E wrt the point p
+ // after the substitutions:
+ // cos(t) = (1 - s^2) / (1 + s^2)
+ // sin(t) = 2t / (1 + s^2)
+ // where s = tan(t/2)
+ // we get a 4th degree equation in s
+ /*
+ * ry s^4 ((-cy + py) Cos[Phi] + (cx - px) Sin[Phi]) +
+ * ry ((cy - py) Cos[Phi] + (-cx + px) Sin[Phi]) +
+ * 2 s^3 (rx^2 - ry^2 + (-cx + px) rx Cos[Phi] + (-cy + py) rx Sin[Phi]) +
+ * 2 s (-rx^2 + ry^2 + (-cx + px) rx Cos[Phi] + (-cy + py) rx Sin[Phi])
+ */
+
+ Point p_c = p - center();
+ double rx2_ry2 = (ray(X) - ray(Y)) * (ray(X) + ray(Y));
+ double cosrot = std::cos( rotation_angle() );
+ double sinrot = std::sin( rotation_angle() );
+ double expr1 = ray(X) * (p_c[X] * cosrot + p_c[Y] * sinrot);
+ Poly coeff;
+ coeff.resize(5);
+ coeff[4] = ray(Y) * ( p_c[Y] * cosrot - p_c[X] * sinrot );
+ coeff[3] = 2 * ( rx2_ry2 + expr1 );
+ coeff[2] = 0;
+ coeff[1] = 2 * ( -rx2_ry2 + expr1 );
+ coeff[0] = -coeff[4];
+
+// for ( unsigned int i = 0; i < 5; ++i )
+// std::cerr << "c[" << i << "] = " << coeff[i] << std::endl;
+
+ std::vector<double> real_sol;
+ // gsl_poly_complex_solve raises an error
+ // if the leading coefficient is zero
+ if ( are_near(coeff[4], 0) )
+ {
+ real_sol.push_back(0);
+ if ( !are_near(coeff[3], 0) )
+ {
+ double sq = -coeff[1] / coeff[3];
+ if ( sq > 0 )
+ {
+ double s = std::sqrt(sq);
+ real_sol.push_back(s);
+ real_sol.push_back(-s);
+ }
+ }
+ }
+ else
+ {
+ real_sol = solve_reals(coeff);
+ }
+
+ for ( unsigned int i = 0; i < real_sol.size(); ++i )
+ {
+ real_sol[i] = 2 * std::atan(real_sol[i]);
+ if ( real_sol[i] < 0 ) real_sol[i] += 2*M_PI;
+ }
+ // when s -> Infinity then <D(E)|E-p> -> 0 iff coeff[4] == 0
+ // so we add M_PI to the solutions being lim arctan(s) = PI when s->Infinity
+ if ( (real_sol.size() % 2) != 0 )
+ {
+ real_sol.push_back(M_PI);
+ }
+
+ double mindistsq1 = std::numeric_limits<double>::max();
+ double mindistsq2 = std::numeric_limits<double>::max();
+ double dsq;
+ unsigned int mi1, mi2;
+ for ( unsigned int i = 0; i < real_sol.size(); ++i )
+ {
+ dsq = distanceSq(p, pointAtAngle(real_sol[i]));
+ if ( mindistsq1 > dsq )
+ {
+ mindistsq2 = mindistsq1;
+ mi2 = mi1;
+ mindistsq1 = dsq;
+ mi1 = i;
+ }
+ else if ( mindistsq2 > dsq )
+ {
+ mindistsq2 = dsq;
+ mi2 = i;
+ }
+ }
+
+ double t = map_to_01( real_sol[mi1] );
+ if ( !(t < from || t > to) )
+ {
+ result.push_back(t);
+ }
+
+ bool second_sol = false;
+ t = map_to_01( real_sol[mi2] );
+ if ( real_sol.size() == 4 && !(t < from || t > to) )
+ {
+ if ( result.empty() || are_near(mindistsq1, mindistsq2) )
+ {
+ result.push_back(t);
+ second_sol = true;
+ }
+ }
+
+ // we need to test extreme points too
+ double dsq1 = distanceSq(p, pointAt(from));
+ double dsq2 = distanceSq(p, pointAt(to));
+ if ( second_sol )
+ {
+ if ( mindistsq2 > dsq1 )
+ {
+ result.clear();
+ result.push_back(from);
+ mindistsq2 = dsq1;
+ }
+ else if ( are_near(mindistsq2, dsq) )
+ {
+ result.push_back(from);
+ }
+ if ( mindistsq2 > dsq2 )
+ {
+ result.clear();
+ result.push_back(to);
+ }
+ else if ( are_near(mindistsq2, dsq2) )
+ {
+ result.push_back(to);
+ }
+
+ }
+ else
+ {
+ if ( result.empty() )
+ {
+ if ( are_near(dsq1, dsq2) )
+ {
+ result.push_back(from);
+ result.push_back(to);
+ }
+ else if ( dsq2 > dsq1 )
+ {
+ result.push_back(from);
+ }
+ else
+ {
+ result.push_back(to);
+ }
+ }
+ }
+
+ return result;
+}
+
+
+/*
+ * NOTE: this implementation follows Standard SVG 1.1 implementation guidelines
+ * for elliptical arc curves. See Appendix F.6.
+ */
+void SVGEllipticalArc::calculate_center_and_extreme_angles()
+{
+ Point d = initialPoint() - finalPoint();
+
+ if (is_svg_compliant())
+ {
+ if ( initialPoint() == finalPoint() )
+ {
+ m_rx = m_ry = m_rot_angle = m_start_angle = m_end_angle = 0;
+ m_center = initialPoint();
+ m_large_arc = m_sweep = false;
+ return;
+ }
+
+ m_rx = std::fabs(m_rx);
+ m_ry = std::fabs(m_ry);
+
+ if ( are_near(ray(X), 0) || are_near(ray(Y), 0) )
+ {
+ m_rx = L2(d) / 2;
+ m_ry = 0;
+ m_rot_angle = std::atan2(d[Y], d[X]);
+ if (m_rot_angle < 0) m_rot_angle += 2*M_PI;
+ m_start_angle = 0;
+ m_end_angle = M_PI;
+ m_center = middle_point(initialPoint(), finalPoint());
+ m_large_arc = false;
+ m_sweep = false;
+ return;
+ }
+ }
+ else
+ {
+ if ( are_near(initialPoint(), finalPoint()) )
+ {
+ if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )
+ {
+ m_start_angle = m_end_angle = 0;
+ m_center = initialPoint();
+ return;
+ }
+ else
+ {
+ THROW_RANGEERROR("initial and final point are the same");
+ }
+ }
+ if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )
+ { // but initialPoint != finalPoint
+ THROW_RANGEERROR(
+ "there is no ellipse that satisfies the given constraints: "
+ "ray(X) == 0 && ray(Y) == 0 but initialPoint != finalPoint"
+ );
+ }
+ if ( are_near(ray(Y), 0) )
+ {
+ Point v = initialPoint() - finalPoint();
+ if ( are_near(L2sq(v), 4*ray(X)*ray(X)) )
+ {
+ double angle = std::atan2(v[Y], v[X]);
+ if (angle < 0) angle += 2*M_PI;
+ if ( are_near( angle, rotation_angle() ) )
+ {
+ m_start_angle = 0;
+ m_end_angle = M_PI;
+ m_center = v/2 + finalPoint();
+ return;
+ }
+ angle -= M_PI;
+ if ( angle < 0 ) angle += 2*M_PI;
+ if ( are_near( angle, rotation_angle() ) )
+ {
+ m_start_angle = M_PI;
+ m_end_angle = 0;
+ m_center = v/2 + finalPoint();
+ return;
+ }
+ THROW_RANGEERROR(
+ "there is no ellipse that satisfies the given constraints: "
+ "ray(Y) == 0 "
+ "and slope(initialPoint - finalPoint) != rotation_angle "
+ "and != rotation_angle + PI"
+ );
+ }
+ if ( L2sq(v) > 4*ray(X)*ray(X) )
+ {
+ THROW_RANGEERROR(
+ "there is no ellipse that satisfies the given constraints: "
+ "ray(Y) == 0 and distance(initialPoint, finalPoint) > 2*ray(X)"
+ );
+ }
+ else
+ {
+ THROW_RANGEERROR(
+ "there is infinite ellipses that satisfy the given constraints: "
+ "ray(Y) == 0 and distance(initialPoint, finalPoint) < 2*ray(X)"
+ );
+ }
+
+ }
+
+ if ( are_near(ray(X), 0) )
+ {
+ Point v = initialPoint() - finalPoint();
+ if ( are_near(L2sq(v), 4*ray(Y)*ray(Y)) )
+ {
+ double angle = std::atan2(v[Y], v[X]);
+ if (angle < 0) angle += 2*M_PI;
+ double rot_angle = rotation_angle() + M_PI/2;
+ if ( !(rot_angle < 2*M_PI) ) rot_angle -= 2*M_PI;
+ if ( are_near( angle, rot_angle ) )
+ {
+ m_start_angle = M_PI/2;
+ m_end_angle = 3*M_PI/2;
+ m_center = v/2 + finalPoint();
+ return;
+ }
+ angle -= M_PI;
+ if ( angle < 0 ) angle += 2*M_PI;
+ if ( are_near( angle, rot_angle ) )
+ {
+ m_start_angle = 3*M_PI/2;
+ m_end_angle = M_PI/2;
+ m_center = v/2 + finalPoint();
+ return;
+ }
+ THROW_RANGEERROR(
+ "there is no ellipse that satisfies the given constraints: "
+ "ray(X) == 0 "
+ "and slope(initialPoint - finalPoint) != rotation_angle + PI/2 "
+ "and != rotation_angle + (3/2)*PI"
+ );
+ }
+ if ( L2sq(v) > 4*ray(Y)*ray(Y) )
+ {
+ THROW_RANGEERROR(
+ "there is no ellipse that satisfies the given constraints: "
+ "ray(X) == 0 and distance(initialPoint, finalPoint) > 2*ray(Y)"
+ );
+ }
+ else
+ {
+ THROW_RANGEERROR(
+ "there is infinite ellipses that satisfy the given constraints: "
+ "ray(X) == 0 and distance(initialPoint, finalPoint) < 2*ray(Y)"
+ );
+ }
+
+ }
+
+ }
+
+ double sin_rot_angle = std::sin(rotation_angle());
+ double cos_rot_angle = std::cos(rotation_angle());
+
+
+ Matrix m( cos_rot_angle, -sin_rot_angle,
+ sin_rot_angle, cos_rot_angle,
+ 0, 0 );
+
+ Point p = (d / 2) * m;
+ double rx2 = m_rx * m_rx;
+ double ry2 = m_ry * m_ry;
+ double rxpy = m_rx * p[Y];
+ double rypx = m_ry * p[X];
+ double rx2py2 = rxpy * rxpy;
+ double ry2px2 = rypx * rypx;
+ double num = rx2 * ry2;
+ double den = rx2py2 + ry2px2;
+ assert(den != 0);
+ double rad = num / den;
+ Point c(0,0);
+ if (rad > 1)
+ {
+ rad -= 1;
+ rad = std::sqrt(rad);
+
+ if (m_large_arc == m_sweep) rad = -rad;
+ c = rad * Point(rxpy / m_ry, -rypx / m_rx);
+
+ m[1] = -m[1];
+ m[2] = -m[2];
+
+ m_center = c * m + middle_point(initialPoint(), finalPoint());
+ }
+ else if (rad == 1 || is_svg_compliant())
+ {
+ double lamda = std::sqrt(1 / rad);
+ m_rx *= lamda;
+ m_ry *= lamda;
+ m_center = middle_point(initialPoint(), finalPoint());
+ }
+ else
+ {
+ THROW_RANGEERROR(
+ "there is no ellipse that satisfies the given constraints"
+ );
+ }
+
+ Point sp((p[X] - c[X]) / m_rx, (p[Y] - c[Y]) / m_ry);
+ Point ep((-p[X] - c[X]) / m_rx, (-p[Y] - c[Y]) / m_ry);
+ Point v(1, 0);
+ m_start_angle = angle_between(v, sp);
+ double sweep_angle = angle_between(sp, ep);
+ if (!m_sweep && sweep_angle > 0) sweep_angle -= 2*M_PI;
+ if (m_sweep && sweep_angle < 0) sweep_angle += 2*M_PI;
+
+ if (m_start_angle < 0) m_start_angle += 2*M_PI;
+ m_end_angle = m_start_angle + sweep_angle;
+ if (m_end_angle < 0) m_end_angle += 2*M_PI;
+ if (m_end_angle >= 2*M_PI) m_end_angle -= 2*M_PI;
+}
+
+
+D2<SBasis> SVGEllipticalArc::toSBasis() const
+{
+
+ if (isDegenerate() && is_svg_compliant())
+ return chord().toSBasis();
+
+ D2<SBasis> arc;
+ // the interval of parametrization has to be [0,1]
+ Coord et = start_angle() + ( sweep_flag() ? sweep_angle() : -sweep_angle() );
+ Linear param(start_angle(), et);
+ Coord cos_rot_angle = std::cos(rotation_angle());
+ Coord sin_rot_angle = std::sin(rotation_angle());
+ // order = 4 seems to be enough to get a perfect looking elliptical arc
+ // should it be choosen in function of the arc length anyway ?
+ // or maybe a user settable parameter: toSBasis(unsigned int order) ?
+ SBasis arc_x = ray(X) * cos(param,4);
+ SBasis arc_y = ray(Y) * sin(param,4);
+ arc[0] = arc_x * cos_rot_angle - arc_y * sin_rot_angle + Linear(center(X),center(X));
+ arc[1] = arc_x * sin_rot_angle + arc_y * cos_rot_angle + Linear(center(Y),center(Y));
+ return arc;
+}
+
+
+Coord SVGEllipticalArc::map_to_02PI(Coord t) const
+{
+ if ( sweep_flag() )
+ {
+ Coord angle = start_angle() + sweep_angle() * t;
+ if ( !(angle < 2*M_PI) )
+ angle -= 2*M_PI;
+ return angle;
+ }
+ else
+ {
+ Coord angle = start_angle() - sweep_angle() * t;
+ if ( angle < 0 ) angle += 2*M_PI;
+ return angle;
+ }
+}
+
+Coord SVGEllipticalArc::map_to_01(Coord angle) const
+{
+ return map_circular_arc_on_unit_interval(angle, start_angle(),
+ end_angle(), sweep_flag());
+}
+
+
+namespace detail
+{
+
+struct ellipse_equation
+{
+ ellipse_equation(double a, double b, double c, double d, double e, double f)
+ : A(a), B(b), C(c), D(d), E(e), F(f)
+ {
+ }
+
+ double operator()(double x, double y) const
+ {
+ // A * x * x + B * x * y + C * y * y + D * x + E * y + F;
+ return (A * x + B * y + D) * x + (C * y + E) * y + F;
+ }
+
+ double operator()(Point const& p) const
+ {
+ return (*this)(p[X], p[Y]);
+ }
+
+ Point normal(double x, double y) const
+ {
+ Point n( 2 * A * x + B * y + D, 2 * C * y + B * x + E );
+ return unit_vector(n);
+ }
+
+ Point normal(Point const& p) const
+ {
+ return normal(p[X], p[Y]);
+ }
+
+ double A, B, C, D, E, F;
+};
+
+}
+
+make_elliptical_arc::
+make_elliptical_arc( SVGEllipticalArc& _ea,
+ curve_type const& _curve,
+ unsigned int _total_samples,
+ double _tolerance )
+ : ea(_ea), curve(_curve),
+ dcurve( unitVector(derivative(curve)) ),
+ model(), fitter(model, _total_samples),
+ tolerance(_tolerance), tol_at_extr(tolerance/2),
+ tol_at_center(0.1), angle_tol(0.1),
+ initial_point(curve.at0()), final_point(curve.at1()),
+ N(_total_samples), last(N-1), partitions(N-1), p(N),
+ svg_compliant(true)
+{
+}
+
+bool
+make_elliptical_arc::
+bound_exceeded( unsigned int k, detail::ellipse_equation const & ee,
+ double e1x, double e1y, double e2 )
+{
+ dist_err = std::fabs( ee(p[k]) );
+ dist_bound = std::fabs( e1x * p[k][X] + e1y * p[k][Y] + e2 );
+ angle_err = std::fabs( dot( dcurve(k/partitions), ee.normal(p[k]) ) );
+ //angle_err *= angle_err;
+ return ( dist_err > dist_bound || angle_err > angle_tol );
+}
+
+bool
+make_elliptical_arc::
+check_bound(double A, double B, double C, double D, double E, double F)
+{
+ // check error magnitude
+ detail::ellipse_equation ee(A, B, C, D, E, F);
+
+ double e1x = (2*A + B) * tol_at_extr;
+ double e1y = (B + 2*C) * tol_at_extr;
+ double e2 = ((D + E) + (A + B + C) * tol_at_extr) * tol_at_extr;
+ if ( bound_exceeded(0, ee, e1x, e1y, e2) )
+ {
+ print_bound_error(0);
+ return false;
+ }
+ if ( bound_exceeded(0, ee, e1x, e1y, e2) )
+ {
+ print_bound_error(last);
+ return false;
+ }
+
+ e1x = (2*A + B) * tolerance;
+ e1y = (B + 2*C) * tolerance;
+ e2 = ((D + E) + (A + B + C) * tolerance) * tolerance;
+// std::cerr << "e1x = " << e1x << std::endl;
+// std::cerr << "e1y = " << e1y << std::endl;
+// std::cerr << "e2 = " << e2 << std::endl;
+
+ for ( unsigned int k = 1; k < last; ++k )
+ {
+ if ( bound_exceeded(k, ee, e1x, e1y, e2) )
+ {
+ print_bound_error(k);
+ return false;
+ }
+ }
+
+ return true;
+}
+
+void make_elliptical_arc::fit()
+{
+ for (unsigned int k = 0; k < N; ++k)
+ {
+ p[k] = curve( k / partitions );
+ fitter.append(p[k]);
+ }
+ fitter.update();
+
+ NL::Vector z(N, 0.0);
+ fitter.result(z);
+}
+
+bool make_elliptical_arc::make_elliptiarc()
+{
+ const NL::Vector & coeff = fitter.result();
+ Ellipse e;
+ try
+ {
+ e.set(1, coeff[0], coeff[1], coeff[2], coeff[3], coeff[4]);
+ }
+ catch(LogicalError exc)
+ {
+ return false;
+ }
+
+ Point inner_point = curve(0.5);
+
+ if (svg_compliant_flag())
+ {
+ ea = e.arc(initial_point, inner_point, final_point);
+ }
+ else
+ {
+ try
+ {
+ ea = e.arc(initial_point, inner_point, final_point, false);
+ }
+ catch(RangeError exc)
+ {
+ return false;
+ }
+ }
+
+
+ if ( !are_near( e.center(),
+ ea.center(),
+ tol_at_center * std::min(e.ray(X),e.ray(Y))
+ )
+ )
+ {
+ return false;
+ }
+ return true;
+}
+
+
+
+} // end namespace Geom
+
+
+
+
+/*
+ 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/svg-elliptical-arc.h b/src/2geom/svg-elliptical-arc.h new file mode 100644 index 000000000..ae6d2e254 --- /dev/null +++ b/src/2geom/svg-elliptical-arc.h @@ -0,0 +1,435 @@ +/*
+ * Elliptical Arc - implementation of the svg elliptical arc path element
+ *
+ * Authors:
+ * MenTaLguY <mental@rydia.net>
+ * Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2007-2008 authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+
+#ifndef _2GEOM_SVG_ELLIPTICAL_ARC_H_
+#define _2GEOM_SVG_ELLIPTICAL_ARC_H_
+
+
+#include <2geom/curve.h>
+#include <2geom/angle.h>
+#include <2geom/utils.h>
+#include <2geom/bezier-curve.h>
+#include <2geom/sbasis-curve.h> // for non-native methods
+#include <2geom/numeric/vector.h>
+#include <2geom/numeric/fitting-tool.h>
+#include <2geom/numeric/fitting-model.h>
+
+
+#include <algorithm>
+
+
+
+namespace Geom
+{
+
+class SVGEllipticalArc : public Curve
+{
+ public:
+ SVGEllipticalArc(bool _svg_compliant = true)
+ : m_initial_point(Point(0,0)), m_final_point(Point(0,0)),
+ m_rx(0), m_ry(0), m_rot_angle(0),
+ m_large_arc(true), m_sweep(true),
+ m_svg_compliant(_svg_compliant)
+ {
+ m_start_angle = m_end_angle = 0;
+ m_center = Point(0,0);
+ }
+
+ SVGEllipticalArc( Point _initial_point, double _rx, double _ry,
+ double _rot_angle, bool _large_arc, bool _sweep,
+ Point _final_point,
+ bool _svg_compliant = true
+ )
+ : m_initial_point(_initial_point), m_final_point(_final_point),
+ m_rx(_rx), m_ry(_ry), m_rot_angle(_rot_angle),
+ m_large_arc(_large_arc), m_sweep(_sweep),
+ m_svg_compliant(_svg_compliant)
+ {
+ calculate_center_and_extreme_angles();
+ }
+
+ void set( Point _initial_point, double _rx, double _ry,
+ double _rot_angle, bool _large_arc, bool _sweep,
+ Point _final_point
+ )
+ {
+ m_initial_point = _initial_point;
+ m_final_point = _final_point;
+ m_rx = _rx;
+ m_ry = _ry;
+ m_rot_angle = _rot_angle;
+ m_large_arc = _large_arc;
+ m_sweep = _sweep;
+ calculate_center_and_extreme_angles();
+ }
+
+ Curve* duplicate() const
+ {
+ return new SVGEllipticalArc(*this);
+ }
+
+ double center(unsigned int i) const
+ {
+ return m_center[i];
+ }
+
+ Point center() const
+ {
+ return m_center;
+ }
+
+ Point initialPoint() const
+ {
+ return m_initial_point;
+ }
+
+ Point finalPoint() const
+ {
+ return m_final_point;
+ }
+
+ double start_angle() const
+ {
+ return m_start_angle;
+ }
+
+ double end_angle() const
+ {
+ return m_end_angle;
+ }
+
+ double ray(unsigned int i) const
+ {
+ return (i == 0) ? m_rx : m_ry;
+ }
+
+ bool large_arc_flag() const
+ {
+ return m_large_arc;
+ }
+
+ bool sweep_flag() const
+ {
+ return m_sweep;
+ }
+
+ double rotation_angle() const
+ {
+ return m_rot_angle;
+ }
+
+ void setInitial( const Point _point)
+ {
+ m_initial_point = _point;
+ calculate_center_and_extreme_angles();
+ }
+
+ void setFinal( const Point _point)
+ {
+ m_final_point = _point;
+ calculate_center_and_extreme_angles();
+ }
+
+ void setExtremes( const Point& _initial_point, const Point& _final_point )
+ {
+ m_initial_point = _initial_point;
+ m_final_point = _final_point;
+ calculate_center_and_extreme_angles();
+ }
+
+ bool isDegenerate() const
+ {
+ return ( are_near(ray(X), 0) || are_near(ray(Y), 0) );
+ }
+
+ bool is_svg_compliant() const
+ {
+ return m_svg_compliant;
+ }
+
+ Rect boundsFast() const
+ {
+ return boundsExact();
+ }
+
+ Rect boundsExact() const;
+
+ // TODO: native implementation of the following methods
+ Rect boundsLocal(Interval i, unsigned int deg) const
+ {
+ if (isDegenerate() && is_svg_compliant())
+ return chord().boundsLocal(i, deg);
+ else
+ return SBasisCurve(toSBasis()).boundsLocal(i, deg);
+ }
+
+ std::vector<double> roots(double v, Dim2 d) const;
+
+ std::vector<double>
+ allNearestPoints( Point const& p, double from = 0, double to = 1 ) const;
+
+ double nearestPoint( Point const& p, double from = 0, double to = 1 ) const
+ {
+ if ( are_near(ray(X), ray(Y)) && are_near(center(), p) )
+ {
+ return from;
+ }
+ return allNearestPoints(p, from, to).front();
+ }
+
+ // TODO: native implementation of the following methods
+ int winding(Point p) const
+ {
+ if (isDegenerate() && is_svg_compliant())
+ return chord().winding(p);
+ else
+ return SBasisCurve(toSBasis()).winding(p);
+ }
+
+ Curve *derivative() const;
+
+ // TODO: native implementation of the following methods
+ Curve *transformed(Matrix const &m) const
+ {
+ return SBasisCurve(toSBasis()).transformed(m);
+ }
+
+ std::vector<Point> pointAndDerivatives(Coord t, unsigned int n) const;
+
+ D2<SBasis> toSBasis() const;
+
+ bool containsAngle(Coord angle) const;
+
+ double valueAtAngle(Coord t, Dim2 d) const;
+
+ Point pointAtAngle(Coord t) const
+ {
+ double sin_rot_angle = std::sin(rotation_angle());
+ double cos_rot_angle = std::cos(rotation_angle());
+ Matrix m( ray(X) * cos_rot_angle, ray(X) * sin_rot_angle,
+ -ray(Y) * sin_rot_angle, ray(Y) * cos_rot_angle,
+ center(X), center(Y) );
+ Point p( std::cos(t), std::sin(t) );
+ return p * m;
+ }
+
+ double valueAt(Coord t, Dim2 d) const
+ {
+ if (isDegenerate() && is_svg_compliant())
+ return chord().valueAt(t, d);
+
+ Coord tt = map_to_02PI(t);
+ return valueAtAngle(tt, d);
+ }
+
+ Point pointAt(Coord t) const
+ {
+ if (isDegenerate() && is_svg_compliant())
+ return chord().pointAt(t);
+
+ Coord tt = map_to_02PI(t);
+ return pointAtAngle(tt);
+ }
+
+ std::pair<SVGEllipticalArc, SVGEllipticalArc>
+ subdivide(Coord t) const
+ {
+ SVGEllipticalArc* arc1 = static_cast<SVGEllipticalArc*>(portion(0, t));
+ SVGEllipticalArc* arc2 = static_cast<SVGEllipticalArc*>(portion(t, 1));
+ assert( arc1 != NULL && arc2 != NULL);
+ std::pair<SVGEllipticalArc, SVGEllipticalArc> arc_pair(*arc1, *arc2);
+ delete arc1;
+ delete arc2;
+ return arc_pair;
+ }
+
+ Curve* portion(double f, double t) const;
+
+ // the arc is the same but traversed in the opposite direction
+ Curve* reverse() const
+ {
+ SVGEllipticalArc* rarc = new SVGEllipticalArc( *this );
+ rarc->m_sweep = !m_sweep;
+ rarc->m_initial_point = m_final_point;
+ rarc->m_final_point = m_initial_point;
+ rarc->m_start_angle = m_end_angle;
+ rarc->m_end_angle = m_start_angle;
+ return rarc;
+ }
+
+ double sweep_angle() const
+ {
+ Coord d = end_angle() - start_angle();
+ if ( !sweep_flag() ) d = -d;
+ if ( d < 0 )
+ d += 2*M_PI;
+ return d;
+ }
+
+ LineSegment chord() const
+ {
+ return LineSegment(initialPoint(), finalPoint());
+ }
+
+ private:
+ Coord map_to_02PI(Coord t) const;
+ Coord map_to_01(Coord angle) const;
+ void calculate_center_and_extreme_angles();
+
+ private:
+ Point m_initial_point, m_final_point;
+ double m_rx, m_ry, m_rot_angle;
+ bool m_large_arc, m_sweep;
+ double m_start_angle, m_end_angle;
+ Point m_center;
+ bool m_svg_compliant;
+
+}; // end class SVGEllipticalArc
+
+template< class charT >
+inline
+std::basic_ostream<charT> &
+operator<< (std::basic_ostream<charT> & os, const SVGEllipticalArc & ea)
+{
+ os << "{ cx: " << ea.center(X) << ", cy: " << ea.center(Y)
+ << ", rx: " << ea.ray(X) << ", ry: " << ea.ray(Y)
+ << ", rot angle: " << decimal_round(rad_to_deg(ea.rotation_angle()),2)
+ << ", start angle: " << decimal_round(rad_to_deg(ea.start_angle()),2)
+ << ", end angle: " << decimal_round(rad_to_deg(ea.end_angle()),2)
+ << " }";
+
+ return os;
+}
+
+
+
+
+namespace detail
+{
+ struct ellipse_equation;
+}
+
+
+class make_elliptical_arc
+{
+ public:
+ typedef D2<SBasis> curve_type;
+
+ make_elliptical_arc( SVGEllipticalArc& _ea,
+ curve_type const& _curve,
+ unsigned int _total_samples,
+ double _tolerance );
+
+ private:
+ bool bound_exceeded( unsigned int k, detail::ellipse_equation const & ee,
+ double e1x, double e1y, double e2 );
+
+ bool check_bound(double A, double B, double C, double D, double E, double F);
+
+ void fit();
+
+ bool make_elliptiarc();
+
+ void print_bound_error(unsigned int k)
+ {
+ std::cerr
+ << "tolerance error" << std::endl
+ << "at point: " << k << std::endl
+ << "error value: "<< dist_err << std::endl
+ << "bound: " << dist_bound << std::endl
+ << "angle error: " << angle_err
+ << " (" << angle_tol << ")" << std::endl;
+ }
+
+ public:
+ bool operator()()
+ {
+ const NL::Vector & coeff = fitter.result();
+ fit();
+ if ( !check_bound(1, coeff[0], coeff[1], coeff[2], coeff[3], coeff[4]) )
+ return false;
+ if ( !(make_elliptiarc()) ) return false;
+ return true;
+ }
+
+ bool svg_compliant_flag() const
+ {
+ return svg_compliant;
+ }
+
+ void svg_compliant_on()
+ {
+ svg_compliant = true;
+ }
+
+ void svg_compliant_off()
+ {
+ svg_compliant = false;
+ }
+
+ private:
+ SVGEllipticalArc& ea;
+ const curve_type & curve;
+ Piecewise<D2<SBasis> > dcurve;
+ NL::LFMEllipse model;
+ NL::least_squeares_fitter<NL::LFMEllipse> fitter;
+ double tolerance, tol_at_extr, tol_at_center, angle_tol;
+ Point initial_point, final_point;
+ unsigned int N;
+ unsigned int last; // N-1
+ double partitions; // N-1
+ std::vector<Point> p; // sample points
+ double dist_err, dist_bound, angle_err;
+ bool svg_compliant;
+};
+
+
+} // end namespace Geom
+
+
+
+
+#endif /* _2GEOM_SVG_ELLIPTICAL_ARC_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/svg-path-parser.cpp b/src/2geom/svg-path-parser.cpp index fe996cef1..2a0ee9b1f 100644 --- a/src/2geom/svg-path-parser.cpp +++ b/src/2geom/svg-path-parser.cpp @@ -1,4 +1,4 @@ -#line 1 "/opt/shared/work/programming/eclipse/eclipse_3.3/lib2geom/src/svg-path-parser.rl" +#line 1 "/home/mental/trees/lib2geom/src/2geom/svg-path-parser.rl" /* * parse SVG path specifications * @@ -35,9 +35,8 @@ #include <vector> #include <glib.h> -#include "point.h" - -#include "svg-path-parser.h" +#include <2geom/point.h> +#include <2geom/svg-path-parser.h> namespace Geom { @@ -139,7 +138,7 @@ private: }; -#line 143 "/opt/shared/work/programming/eclipse/eclipse_3.3/lib2geom/src/svg-path-parser.cpp" +#line 142 "/home/mental/trees/lib2geom/src/2geom/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, @@ -1144,7 +1143,7 @@ static const int svg_path_first_final = 270; static const int svg_path_en_main = 1; -#line 143 "/opt/shared/work/programming/eclipse/eclipse_3.3/lib2geom/src/svg-path-parser.rl" +#line 142 "/home/mental/trees/lib2geom/src/2geom/svg-path-parser.rl" void Parser::parse(char const *str) @@ -1157,12 +1156,12 @@ throw(SVGPathParseError) _reset(); -#line 1161 "/opt/shared/work/programming/eclipse/eclipse_3.3/lib2geom/src/svg-path-parser.cpp" +#line 1160 "/home/mental/trees/lib2geom/src/2geom/svg-path-parser.cpp" { cs = svg_path_start; } -#line 1166 "/opt/shared/work/programming/eclipse/eclipse_3.3/lib2geom/src/svg-path-parser.cpp" +#line 1165 "/home/mental/trees/lib2geom/src/2geom/svg-path-parser.cpp" { int _klen; unsigned int _trans; @@ -1235,13 +1234,13 @@ _match: switch ( *_acts++ ) { case 0: -#line 155 "/opt/shared/work/programming/eclipse/eclipse_3.3/lib2geom/src/svg-path-parser.rl" +#line 154 "/home/mental/trees/lib2geom/src/2geom/svg-path-parser.rl" { start = p; } break; case 1: -#line 159 "/opt/shared/work/programming/eclipse/eclipse_3.3/lib2geom/src/svg-path-parser.rl" +#line 158 "/home/mental/trees/lib2geom/src/2geom/svg-path-parser.rl" { char const *end=p; std::string buf(start, end); @@ -1250,55 +1249,55 @@ _match: } break; case 2: -#line 166 "/opt/shared/work/programming/eclipse/eclipse_3.3/lib2geom/src/svg-path-parser.rl" +#line 165 "/home/mental/trees/lib2geom/src/2geom/svg-path-parser.rl" { _push(1.0); } break; case 3: -#line 170 "/opt/shared/work/programming/eclipse/eclipse_3.3/lib2geom/src/svg-path-parser.rl" +#line 169 "/home/mental/trees/lib2geom/src/2geom/svg-path-parser.rl" { _push(0.0); } break; case 4: -#line 174 "/opt/shared/work/programming/eclipse/eclipse_3.3/lib2geom/src/svg-path-parser.rl" +#line 173 "/home/mental/trees/lib2geom/src/2geom/svg-path-parser.rl" { _absolute = true; } break; case 5: -#line 178 "/opt/shared/work/programming/eclipse/eclipse_3.3/lib2geom/src/svg-path-parser.rl" +#line 177 "/home/mental/trees/lib2geom/src/2geom/svg-path-parser.rl" { _absolute = false; } break; case 6: -#line 182 "/opt/shared/work/programming/eclipse/eclipse_3.3/lib2geom/src/svg-path-parser.rl" +#line 181 "/home/mental/trees/lib2geom/src/2geom/svg-path-parser.rl" { _moveTo(_pop_point()); } break; case 7: -#line 186 "/opt/shared/work/programming/eclipse/eclipse_3.3/lib2geom/src/svg-path-parser.rl" +#line 185 "/home/mental/trees/lib2geom/src/2geom/svg-path-parser.rl" { _lineTo(_pop_point()); } break; case 8: -#line 190 "/opt/shared/work/programming/eclipse/eclipse_3.3/lib2geom/src/svg-path-parser.rl" +#line 189 "/home/mental/trees/lib2geom/src/2geom/svg-path-parser.rl" { _hlineTo(Point(_pop_coord(X), _current[Y])); } break; case 9: -#line 194 "/opt/shared/work/programming/eclipse/eclipse_3.3/lib2geom/src/svg-path-parser.rl" +#line 193 "/home/mental/trees/lib2geom/src/2geom/svg-path-parser.rl" { _vlineTo(Point(_current[X], _pop_coord(Y))); } break; case 10: -#line 198 "/opt/shared/work/programming/eclipse/eclipse_3.3/lib2geom/src/svg-path-parser.rl" +#line 197 "/home/mental/trees/lib2geom/src/2geom/svg-path-parser.rl" { Point p = _pop_point(); Point c1 = _pop_point(); @@ -1307,7 +1306,7 @@ _match: } break; case 11: -#line 205 "/opt/shared/work/programming/eclipse/eclipse_3.3/lib2geom/src/svg-path-parser.rl" +#line 204 "/home/mental/trees/lib2geom/src/2geom/svg-path-parser.rl" { Point p = _pop_point(); Point c1 = _pop_point(); @@ -1315,7 +1314,7 @@ _match: } break; case 12: -#line 211 "/opt/shared/work/programming/eclipse/eclipse_3.3/lib2geom/src/svg-path-parser.rl" +#line 210 "/home/mental/trees/lib2geom/src/2geom/svg-path-parser.rl" { Point p = _pop_point(); Point c = _pop_point(); @@ -1323,14 +1322,14 @@ _match: } break; case 13: -#line 217 "/opt/shared/work/programming/eclipse/eclipse_3.3/lib2geom/src/svg-path-parser.rl" +#line 216 "/home/mental/trees/lib2geom/src/2geom/svg-path-parser.rl" { Point p = _pop_point(); _quadTo(_quad_tangent, p); } break; case 14: -#line 222 "/opt/shared/work/programming/eclipse/eclipse_3.3/lib2geom/src/svg-path-parser.rl" +#line 221 "/home/mental/trees/lib2geom/src/2geom/svg-path-parser.rl" { Point point = _pop_point(); bool sweep = _pop_flag(); @@ -1343,16 +1342,16 @@ _match: } break; case 15: -#line 233 "/opt/shared/work/programming/eclipse/eclipse_3.3/lib2geom/src/svg-path-parser.rl" +#line 232 "/home/mental/trees/lib2geom/src/2geom/svg-path-parser.rl" { _closePath(); } break; case 16: -#line 369 "/opt/shared/work/programming/eclipse/eclipse_3.3/lib2geom/src/svg-path-parser.rl" +#line 368 "/home/mental/trees/lib2geom/src/2geom/svg-path-parser.rl" {goto _out;} break; -#line 1356 "/opt/shared/work/programming/eclipse/eclipse_3.3/lib2geom/src/svg-path-parser.cpp" +#line 1355 "/home/mental/trees/lib2geom/src/2geom/svg-path-parser.cpp" } } @@ -1363,7 +1362,7 @@ _again: goto _resume; _out: {} } -#line 379 "/opt/shared/work/programming/eclipse/eclipse_3.3/lib2geom/src/svg-path-parser.rl" +#line 378 "/home/mental/trees/lib2geom/src/2geom/svg-path-parser.rl" if ( cs < svg_path_first_final ) { diff --git a/src/2geom/svg-path-parser.h b/src/2geom/svg-path-parser.h index 98a885361..d24d25551 100644 --- a/src/2geom/svg-path-parser.h +++ b/src/2geom/svg-path-parser.h @@ -35,9 +35,9 @@ #include <vector> #include <iterator> #include <stdexcept> -#include "exception.h" -#include "point.h" -#include "svg-path.h" +#include <2geom/exception.h> +#include <2geom/point.h> +#include <2geom/svg-path.h> namespace Geom { diff --git a/src/2geom/svg-path.cpp b/src/2geom/svg-path.cpp index 819b53aac..c178b92da 100644 --- a/src/2geom/svg-path.cpp +++ b/src/2geom/svg-path.cpp @@ -28,9 +28,9 @@ * */ -#include "sbasis-to-bezier.h" -#include "svg-path.h" -#include "exception.h" +#include <2geom/sbasis-to-bezier.h> +#include <2geom/svg-path.h> +#include <2geom/exception.h> namespace Geom { diff --git a/src/2geom/svg-path.h b/src/2geom/svg-path.h index ad78680de..1ae4b054c 100644 --- a/src/2geom/svg-path.h +++ b/src/2geom/svg-path.h @@ -31,7 +31,7 @@ #ifndef SEEN_SVG_PATH_H #define SEEN_SVG_PATH_H -#include "path.h" +#include <2geom/path.h> #include <iterator> namespace Geom { diff --git a/src/2geom/sweep.cpp b/src/2geom/sweep.cpp index 92da524ce..1b8813990 100644 --- a/src/2geom/sweep.cpp +++ b/src/2geom/sweep.cpp @@ -1,4 +1,4 @@ -#include "sweep.h" +#include <2geom/sweep.h> #include <algorithm> diff --git a/src/2geom/sweep.h b/src/2geom/sweep.h index 94ea6855c..6500ccbc2 100644 --- a/src/2geom/sweep.h +++ b/src/2geom/sweep.h @@ -2,7 +2,7 @@ #define __2GEOM_SWEEP_H__ #include <vector> -#include "d2.h" +#include <2geom/d2.h> namespace Geom { diff --git a/src/2geom/transforms.cpp b/src/2geom/transforms.cpp index 62b340221..a6426fe81 100644 --- a/src/2geom/transforms.cpp +++ b/src/2geom/transforms.cpp @@ -1,4 +1,4 @@ -#include "transforms.h" +#include <2geom/transforms.h> namespace Geom { diff --git a/src/2geom/transforms.h b/src/2geom/transforms.h index fb077a910..ca0a81ac7 100644 --- a/src/2geom/transforms.h +++ b/src/2geom/transforms.h @@ -1,7 +1,7 @@ #ifndef SEEN_Geom_TRANSFORMS_H #define SEEN_Geom_TRANSFORMS_H -#include "matrix.h" +#include <2geom/matrix.h> #include <cmath> namespace Geom { diff --git a/src/2geom/utils.h b/src/2geom/utils.h index 96c63fddd..ac99b8b00 100644 --- a/src/2geom/utils.h +++ b/src/2geom/utils.h @@ -32,6 +32,7 @@ */ #include <cmath> +#include <vector> namespace Geom { @@ -76,6 +77,9 @@ inline double decimal_round(double const x, int const places) { return round( x * multiplier ) / multiplier; } + +void binomial_coefficients(std::vector<size_t>& bc, size_t n); + } #endif diff --git a/src/live_effects/lpe-bendpath.cpp b/src/live_effects/lpe-bendpath.cpp index 282ea8270..dc98561e1 100644 --- a/src/live_effects/lpe-bendpath.cpp +++ b/src/live_effects/lpe-bendpath.cpp @@ -94,7 +94,7 @@ LPEBendPath::doEffect_pwd2 (Geom::Piecewise<Geom::D2<Geom::SBasis> > const & pwd Piecewise<D2<SBasis> > n = rot90(derivative(uskeleton)); n = force_continuity(remove_short_cuts(n,.1)); - D2<Piecewise<SBasis> > patternd2 = make_cuts_independant(pwd2_in); + D2<Piecewise<SBasis> > patternd2 = make_cuts_independent(pwd2_in); Piecewise<SBasis> x = vertical_pattern.get_value() ? Piecewise<SBasis>(patternd2[1]) : Piecewise<SBasis>(patternd2[0]); Piecewise<SBasis> y = vertical_pattern.get_value() ? Piecewise<SBasis>(patternd2[0]) : Piecewise<SBasis>(patternd2[1]); diff --git a/src/live_effects/lpe-curvestitch.cpp b/src/live_effects/lpe-curvestitch.cpp index 814a381af..34fb60a31 100644 --- a/src/live_effects/lpe-curvestitch.cpp +++ b/src/live_effects/lpe-curvestitch.cpp @@ -78,7 +78,7 @@ LPECurveStitch::doEffect_path (std::vector<Geom::Path> const & path_in) startpoint_spacing_variation.resetRandomizer(); endpoint_spacing_variation.resetRandomizer(); - D2<Piecewise<SBasis> > stroke = make_cuts_independant(strokepath.get_pwd2()); + D2<Piecewise<SBasis> > stroke = make_cuts_independent(strokepath.get_pwd2()); Interval bndsStroke = bounds_exact(stroke[0]); gdouble scaling = bndsStroke.max() - bndsStroke.min(); Interval bndsStrokeY = bounds_exact(stroke[1]); @@ -162,7 +162,7 @@ LPECurveStitch::resetDefaults(SPItem * item) for (unsigned int i=0; i < temppath.size(); i++) { pwd2.concat( temppath[i].toPwSb() ); } - D2<Piecewise<SBasis> > d2pw = make_cuts_independant(pwd2); + D2<Piecewise<SBasis> > d2pw = make_cuts_independent(pwd2); Interval bndsX = bounds_exact(d2pw[0]); Interval bndsY = bounds_exact(d2pw[1]); Point start(bndsX.min(), (bndsY.max()+bndsY.min())/2); diff --git a/src/live_effects/lpe-envelope.cpp b/src/live_effects/lpe-envelope.cpp index 9ebe97b9f..7d2045d80 100755 --- a/src/live_effects/lpe-envelope.cpp +++ b/src/live_effects/lpe-envelope.cpp @@ -98,7 +98,7 @@ LPEEnvelope::doEffect_pwd2 (Geom::Piecewise<Geom::D2<Geom::SBasis> > const & pwd n4 = force_continuity(remove_short_cuts(n4,.1)); - D2<Piecewise<SBasis> > patternd2 = make_cuts_independant(pwd2_in); + D2<Piecewise<SBasis> > patternd2 = make_cuts_independent(pwd2_in); Piecewise<SBasis> x = Piecewise<SBasis>(patternd2[0]); Piecewise<SBasis> y = Piecewise<SBasis>(patternd2[1]); diff --git a/src/live_effects/lpe-patternalongpath.cpp b/src/live_effects/lpe-patternalongpath.cpp index 85de609fb..2e10efe7a 100644 --- a/src/live_effects/lpe-patternalongpath.cpp +++ b/src/live_effects/lpe-patternalongpath.cpp @@ -98,7 +98,7 @@ LPEPatternAlongPath::doEffect_pwd2 (Geom::Piecewise<Geom::D2<Geom::SBasis> > con PAPCopyType type = copytype.get_value(); - D2<Piecewise<SBasis> > patternd2 = make_cuts_independant(pattern.get_pwd2()); + D2<Piecewise<SBasis> > patternd2 = make_cuts_independent(pattern.get_pwd2()); Piecewise<SBasis> x0 = vertical_pattern.get_value() ? Piecewise<SBasis>(patternd2[1]) : Piecewise<SBasis>(patternd2[0]); Piecewise<SBasis> y0 = vertical_pattern.get_value() ? Piecewise<SBasis>(patternd2[0]) : Piecewise<SBasis>(patternd2[1]); Interval pattBndsX = bounds_exact(x0); diff --git a/src/live_effects/lpe-perspective_path.cpp b/src/live_effects/lpe-perspective_path.cpp index 785d7a245..9aa26e05d 100644 --- a/src/live_effects/lpe-perspective_path.cpp +++ b/src/live_effects/lpe-perspective_path.cpp @@ -97,7 +97,7 @@ LPEPerspectivePath::doEffect_pwd2 (Geom::Piecewise<Geom::D2<Geom::SBasis> > cons // remove this once we have unified coordinate systems path_a_pw = path_a_pw + Geom::Point(offsetx, -offsety); - D2<Piecewise<SBasis> > B = make_cuts_independant(path_a_pw); + D2<Piecewise<SBasis> > B = make_cuts_independent(path_a_pw); Piecewise<SBasis> preimage[4]; //Geom::Point orig = Geom::Point(bounds_X.min(), bounds_Y.middle()); diff --git a/src/live_effects/lpe-vonkoch.cpp b/src/live_effects/lpe-vonkoch.cpp index 114dd7263..e44868290 100644 --- a/src/live_effects/lpe-vonkoch.cpp +++ b/src/live_effects/lpe-vonkoch.cpp @@ -164,7 +164,7 @@ LPEVonKoch::resetDefaults(SPItem * item) pwd2.concat( temppath[i].toPwSb() ); } - D2<Piecewise<SBasis> > d2pw = make_cuts_independant(pwd2); + D2<Piecewise<SBasis> > d2pw = make_cuts_independent(pwd2); Interval bndsX = bounds_exact(d2pw[0]); Interval bndsY = bounds_exact(d2pw[1]); Point start(bndsX.min(), (bndsY.max()+bndsY.min())/2); |