From 8001ba81cb851b38d86650a2fef5817facffb763 Mon Sep 17 00:00:00 2001 From: johanengelen Date: Thu, 3 Jul 2008 20:06:40 +0000 Subject: [PATCH] update 2geom (svn rev1433) --- src/2geom/basic-intersection.cpp | 6 +- src/2geom/basic-intersection.h | 8 +- src/2geom/bezier-curve.h | 6 +- src/2geom/bezier-to-sbasis.h | 6 +- src/2geom/bezier-utils.cpp | 4 +- src/2geom/bezier-utils.h | 2 +- src/2geom/bezier.h | 10 +- src/2geom/choose.h | 2 +- src/2geom/circle-circle.cpp | 2 +- src/2geom/concepts.h | 6 +- src/2geom/conjugate_gradient.cpp | 7 +- src/2geom/convex-cover.cpp | 24 +- src/2geom/convex-cover.h | 2 +- src/2geom/crossing.cpp | 6 +- src/2geom/crossing.h | 4 +- src/2geom/curve-helpers.cpp | 4 +- src/2geom/curve.h | 25 +- src/2geom/curves.h | 13 +- src/2geom/d2-sbasis.cpp | 4 +- src/2geom/d2-sbasis.h | 10 +- src/2geom/d2.h | 12 +- src/2geom/ellipse.cpp | 206 ++++ src/2geom/ellipse.h | 133 +++ src/2geom/elliptical-arc.cpp | 6 +- src/2geom/elliptical-arc.h | 8 +- src/2geom/geom.cpp | 4 +- src/2geom/geom.h | 2 +- src/2geom/hvlinesegment.h | 2 +- src/2geom/interval.h | 2 +- src/2geom/isnan.h | 4 +- src/2geom/linear.h | 4 +- src/2geom/matrix.cpp | 6 +- src/2geom/matrix.h | 2 +- src/2geom/nearest-point.cpp | 8 +- src/2geom/nearest-point.h | 6 +- src/2geom/numeric/fitting-model.h | 423 ++++++++ src/2geom/numeric/fitting-tool.h | 532 ++++++++++ src/2geom/numeric/linear_system.h | 227 +++-- src/2geom/numeric/matrix.cpp | 115 +++ src/2geom/numeric/matrix.h | 762 ++++++++++---- src/2geom/numeric/vector.h | 750 ++++++++++---- src/2geom/path-intersection.cpp | 8 +- src/2geom/path-intersection.h | 6 +- src/2geom/path.cpp | 141 ++- src/2geom/path.h | 263 +++-- src/2geom/pathvector.cpp | 6 +- src/2geom/pathvector.h | 6 +- src/2geom/piecewise.cpp | 17 +- src/2geom/piecewise.h | 6 +- src/2geom/point-l.h | 2 +- src/2geom/point.cpp | 8 +- src/2geom/point.h | 12 +- src/2geom/poly-dk-solve.cpp | 2 +- src/2geom/poly-dk-solve.h | 2 +- src/2geom/poly-laguerre-solve.cpp | 12 +- src/2geom/poly-laguerre-solve.h | 2 +- src/2geom/poly.cpp | 4 +- src/2geom/poly.h | 2 +- src/2geom/quadtree.cpp | 2 +- src/2geom/quadtree.h | 2 +- src/2geom/rect.h | 2 +- src/2geom/region.cpp | 6 +- src/2geom/region.h | 4 +- src/2geom/sbasis-2d.cpp | 2 +- src/2geom/sbasis-2d.h | 4 +- src/2geom/sbasis-curve.h | 2 +- src/2geom/sbasis-geometric.cpp | 10 +- src/2geom/sbasis-geometric.h | 4 +- src/2geom/sbasis-math.cpp | 2 +- src/2geom/sbasis-math.h | 4 +- src/2geom/sbasis-poly.cpp | 2 +- src/2geom/sbasis-poly.h | 4 +- src/2geom/sbasis-roots.cpp | 6 +- src/2geom/sbasis-to-bezier.cpp | 8 +- src/2geom/sbasis-to-bezier.h | 4 +- src/2geom/sbasis.cpp | 7 +- src/2geom/sbasis.h | 8 +- src/2geom/shape.cpp | 8 +- src/2geom/shape.h | 2 +- src/2geom/solve-bezier-one-d.cpp | 4 +- src/2geom/solve-bezier-parametric.cpp | 4 +- src/2geom/solver.h | 4 +- src/2geom/sturm.h | 4 +- src/2geom/svg-elliptical-arc.cpp | 1124 +++++++++++++++++++++ src/2geom/svg-elliptical-arc.h | 435 ++++++++ src/2geom/svg-path-parser.cpp | 53 +- src/2geom/svg-path-parser.h | 6 +- src/2geom/svg-path.cpp | 6 +- src/2geom/svg-path.h | 2 +- src/2geom/sweep.cpp | 2 +- src/2geom/sweep.h | 2 +- src/2geom/transforms.cpp | 2 +- src/2geom/transforms.h | 2 +- src/2geom/utils.h | 4 + src/live_effects/lpe-bendpath.cpp | 2 +- src/live_effects/lpe-curvestitch.cpp | 4 +- src/live_effects/lpe-envelope.cpp | 2 +- src/live_effects/lpe-patternalongpath.cpp | 2 +- src/live_effects/lpe-perspective_path.cpp | 2 +- src/live_effects/lpe-vonkoch.cpp | 2 +- 100 files changed, 4734 insertions(+), 891 deletions(-) create mode 100644 src/2geom/ellipse.cpp create mode 100644 src/2geom/ellipse.h create mode 100644 src/2geom/numeric/fitting-model.h create mode 100644 src/2geom/numeric/fitting-tool.h create mode 100644 src/2geom/numeric/matrix.cpp create mode 100644 src/2geom/svg-elliptical-arc.cpp create mode 100644 src/2geom/svg-elliptical-arc.h 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 const & A, } std::vector > -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 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 #endif -#include "bezier-utils.h" +#include <2geom/bezier-utils.h> -#include "isnan.h" +#include <2geom/isnan.h> #include 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 -#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 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 #include -#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 #include 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 #include #include -#include "conjugate_gradient.h" +#include <2geom/conjugate_gradient.h> /* lifted wholely from wikipedia. */ @@ -55,9 +55,12 @@ matrix_times_vector(valarray const &matrix, /* m * n */ } } +/** +// only used in commented code below static double Linfty(valarray const &vec) { return std::max(vec.max(), -vec.min()); } +**/ double inner(valarray const &x, @@ -96,7 +99,7 @@ conjugate_gradient(valarray const &A, valarray &x, valarray const &b, unsigned n, double tol, - unsigned max_iterations, bool ortho1) { + unsigned max_iterations, bool /*ortho1*/) { valarray 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 #include /** 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 abridges; map 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 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 -#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 @@ -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 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 * Marco Cecchetti - * + * * 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 > sectionize(D2 > const &a) { return ret; } -D2 > make_cuts_independant(Piecewise > const &a) { +D2 > make_cuts_independent(Piecewise > const &a) { D2 > 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 const & a, unsigned tail); //Piecewise > specific decls: Piecewise > sectionize(D2 > const &a); -D2 > make_cuts_independant(Piecewise > const &a); +D2 > make_cuts_independent(Piecewise > const &a); Piecewise > rot90(Piecewise > const &a); Piecewise dot(Piecewise > const &a, Piecewise > const &b); Piecewise cross(Piecewise > const &a, Piecewise > 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 -#include "concepts.h" +#include <2geom/concepts.h> namespace Geom{ @@ -373,8 +373,8 @@ D2 integral(D2 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 + * + * 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 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 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 + * + * 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 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 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 #include 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 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 #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 -#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 -#include "coord.h" +#include <2geom/coord.h> #include 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 #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 -#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 all_nearest_points( Point const& p, - D2 const& c, - D2 const& dc, - double from, double to ) + D2 const& c, + D2 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 -#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 + * + * 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 +{ + 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 +{ + public: + LFMPoly(size_t degree) + : LFMPowerBasis(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 +{ + 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 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 +{ + 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 +{ + 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 +{ + 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 objects +class LFMD2SBasis + : public LinearFittingModel< double, Point, D2 > +{ + 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 & 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 +{ + 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 m_bc; +}; + + +// this model generates Bezier curves +template< unsigned int N > +class LFMBezierCurve + : public LinearFittingModel< double, Point, BezierCurve > +{ + 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 & 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(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 + * + * 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 + + +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 + : public lsf_base +{ +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 base_type; + + using base_type::m_model; + using base_type::m_psdinv_matrix; + using base_type::total_samples; + +public: + lsf_solution( 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 + + +// 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 + : public lsf_base +{ +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 base_type; + + using base_type::m_model; + using base_type::m_psdinv_matrix; + using base_type::total_samples; + +public: + lsf_solution( model_type const& _model, + size_t forecasted_samples ) + : base_type(_model, forecasted_samples), + m_solution(_model.size(), 2) + { + } + + solution_type& result(std::vector 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 const& old_sample_values, + std::vector 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 + + + + +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 + : public lsf_solution +{ + public: + typedef ModelT model_type; + typedef typename model_type::parameter_type parameter_type; + typedef typename model_type::value_type value_type; + typedef lsf_solution 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( 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 + + +// fitting tool for partially known models +template< typename ModelT> +class lsf_with_fixed_terms + : public lsf_solution +{ + public: + typedef ModelT model_type; + typedef typename model_type::parameter_type parameter_type; + typedef typename model_type::value_type value_type; + typedef lsf_solution 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( 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() + { + if (m_vector_view != NULL) + { + delete m_vector_view; + } + } + + protected: + Vector m_vector; + VectorView* m_vector_view; + +}; // end class lsf_with_fixed_terms + + +} // 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 + : public detail::lsf_with_fixed_terms +{ + public: + typedef ModelT model_type; + typedef detail::lsf_with_fixed_terms 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( model_type const& _model, + size_t forecasted_samples ) + : base_type(_model, forecasted_samples) + { + } +}; // end class least_squeares_fitter + + +template< typename ModelT> +class least_squeares_fitter + : public detail::lsf_with_fixed_terms +{ + public: + typedef ModelT model_type; + typedef detail::lsf_with_fixed_terms 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( 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 + + +template< typename ModelT> +class least_squeares_fitter + : public detail::lsf_with_fixed_terms +{ + public: + typedef ModelT model_type; + typedef detail::lsf_with_fixed_terms 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( model_type const& _model, + size_t forecasted_samples ) + : base_type(_model, forecasted_samples) + { + } + + solution_type& result(std::vector 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 + + +} // 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 - -#include - -#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 + * + * 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 + +#include + +#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 + * + * 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 +#include + + + +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 -#include -#include - -#include - - - -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 - max_index() const - { - std::pair indices; - gsl_matrix_max_index(m_matrix, &(indices.first), &(indices.second)); - return indices; - } - - std::pair - min_index() const - { - std::pair 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 + * + * 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 +#include // for std::pair +#include // for std::swap +#include +#include + +#include +#include + + +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 + max_index() const + { + std::pair indices; + gsl_matrix_max_index(m_matrix, &(indices.first), &(indices.second)); + return indices; + } + + std::pair + min_index() const + { + std::pair 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 & +operator<< (std::basic_ostream & 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( base_type::transpose() ); + } + + Matrix & scale(double x) + { + return static_cast( base_type::scale(x) ); + } + + Matrix & translate(double x) + { + return static_cast( base_type::translate(x) ); + } + + Matrix & operator+=(base_type::base_type const& _matrix) + { + return static_cast( base_type::operator+=(_matrix) ); + } + + Matrix & operator-=(base_type::base_type const& _matrix) + { + return static_cast( 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( &(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( &(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( &(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( base_type::transpose() ); + } + + MatrixView & scale(double x) + { + return static_cast( base_type::scale(x) ); + } + + MatrixView & translate(double x) + { + return static_cast( base_type::translate(x) ); + } + + MatrixView & operator+=(base_type::base_type const& _matrix) + { + return static_cast( base_type::operator+=(_matrix) ); + } + + MatrixView & operator-=(base_type::base_type const& _matrix) + { + return static_cast( 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 -#include - -#include - - -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 + * + * 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 +#include // for std::swap +#include +#include +#include + + +#include +#include + + +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 & +operator<< (std::basic_ostream & 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( base_type::scale(x) ); + } + + Vector & translate(double x) + { + return static_cast( base_type::translate(x) ); + } + + Vector & operator+=(base_type::base_type const& _vector) + { + return static_cast( base_type::operator+=(_vector) ); + } + + Vector & operator-=(base_type::base_type const& _vector) + { + return static_cast( 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( &(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( &(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( &(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( &(m_vector_view.vector) ); + } + + explicit + ConstVectorView(gsl_vector_const_view _gsl_vector_view) + : m_vector_view(_gsl_vector_view) + { + m_vector = const_cast( &(m_vector_view.vector) ); + m_size = m_vector->size; + } + + explicit + ConstVectorView(const std::vector& _vector) + : m_vector_view( gsl_vector_const_view_array(&(_vector[0]), _vector.size()) ) + { + m_vector = const_cast( &(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( &(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( &(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 & _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( base_type::scale(x) ); + } + + VectorView & translate(double x) + { + return static_cast( base_type::translate(x) ); + } + + VectorView & operator+=(base_type::base_type const& _vector) + { + return static_cast( base_type::operator+=(_vector) ); + } + + VectorView & operator-=(base_type::base_type const& _vector) + { + return static_cast( 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 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 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 @@ -46,6 +46,16 @@ namespace Geom { +// Conditional expression for types. If true, first, if false, second. +template + struct __conditional_type + { typedef _Iftrue __type; }; + +template + struct __conditional_type + { typedef _Iffalse __type; }; + + template class BaseIterator : public std::iterator @@ -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 Sequence; + BaseIterator ( typename __conditional_type< + (std::__are_same::__value), // check if this instantiation is of const_iterator type + const BaseIterator< Sequence::iterator >, // if true: accept iterator in const_iterator instantiation + const BaseIterator > ::__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 Path(BaseIterator first, BaseIterator 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 - void insert(iterator pos, BaseIterator first, BaseIterator last) + void insert(iterator pos, BaseIterator first, BaseIterator last, Stitching stitching=NO_STITCHING) { Sequence source(DuplicatingIterator(first.impl_), DuplicatingIterator(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 void replace(iterator replaced, - BaseIterator first, BaseIterator last) + BaseIterator first, BaseIterator last, + Stitching stitching=NO_STITCHING) { Sequence source(DuplicatingIterator(first.impl_), DuplicatingIterator(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 void replace(iterator first_replaced, iterator last_replaced, - BaseIterator first, BaseIterator last) + BaseIterator first, BaseIterator 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 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 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 void appendNew(A a) { - do_append(new CurveType((*final_)[0], a)); + do_append(new CurveType(finalPoint(), a)); } template void appendNew(A a, B b) { - do_append(new CurveType((*final_)[0], a, b)); + do_append(new CurveType(finalPoint(), a, b)); } template 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 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 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 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 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 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 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 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 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 pointAndDerivatives(Coord t, unsigned n) const = 0; - virtual D2 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 #include @@ -167,21 +167,6 @@ std::vector roots(Piecewise const &f){ return result; } -Piecewise reverse(Piecewise const &f) { - Piecewise ret = Piecewise(); - 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 #include -#include "concepts.h" -#include "isnan.h" +#include <2geom/concepts.h> +#include <2geom/isnan.h> #include 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 -#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 -#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 -#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 /*** 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 std::vector > 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 typedef std::complex 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 -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 std::vector > 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 #include #include -#include "utils.h" +#include <2geom/utils.h> class Poly : public std::vector{ 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 #include -#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 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 #include #include -#include "d2.h" -#include "sbasis.h" +#include <2geom/d2.h> +#include <2geom/sbasis.h> #include 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 (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 /** 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. //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 #include -#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 -#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 -#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 #include -#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 #include 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 #include -#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 /*** 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 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 + * + * 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 +#include + +#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 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::vectorarc_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 +SVGEllipticalArc::roots(double v, Dim2 d) const +{ + if ( d > Y ) + { + THROW_RANGEERROR("dimention out of range"); + } + + std::vector 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 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 +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 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 SVGEllipticalArc:: +allNearestPoints( Point const& p, double from, double to ) const +{ + std::vector 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 solX = roots(np[X],X); +// std::vector 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 == 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 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 -> 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::max(); + double mindistsq2 = std::numeric_limits::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 SVGEllipticalArc::toSBasis() const +{ + + if (isDegenerate() && is_svg_compliant()) + return chord().toSBasis(); + + D2 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 + * Marco Cecchetti + * + * 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 + + + +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 roots(double v, Dim2 d) const; + + std::vector + 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 pointAndDerivatives(Coord t, unsigned int n) const; + + D2 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 + subdivide(Coord t) const + { + SVGEllipticalArc* arc1 = static_cast(portion(0, t)); + SVGEllipticalArc* arc2 = static_cast(portion(t, 1)); + assert( arc1 != NULL && arc2 != NULL); + std::pair 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 & +operator<< (std::basic_ostream & 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 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 > dcurve; + NL::LFMEllipse model; + NL::least_squeares_fitter 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 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 #include -#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 #include #include -#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 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 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 -#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 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 +#include 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& 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 > const & pwd Piecewise > n = rot90(derivative(uskeleton)); n = force_continuity(remove_short_cuts(n,.1)); - D2 > patternd2 = make_cuts_independant(pwd2_in); + D2 > patternd2 = make_cuts_independent(pwd2_in); Piecewise x = vertical_pattern.get_value() ? Piecewise(patternd2[1]) : Piecewise(patternd2[0]); Piecewise y = vertical_pattern.get_value() ? Piecewise(patternd2[0]) : Piecewise(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 const & path_in) startpoint_spacing_variation.resetRandomizer(); endpoint_spacing_variation.resetRandomizer(); - D2 > stroke = make_cuts_independant(strokepath.get_pwd2()); + D2 > 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 > d2pw = make_cuts_independant(pwd2); + D2 > 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 > const & pwd n4 = force_continuity(remove_short_cuts(n4,.1)); - D2 > patternd2 = make_cuts_independant(pwd2_in); + D2 > patternd2 = make_cuts_independent(pwd2_in); Piecewise x = Piecewise(patternd2[0]); Piecewise y = Piecewise(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 > con PAPCopyType type = copytype.get_value(); - D2 > patternd2 = make_cuts_independant(pattern.get_pwd2()); + D2 > patternd2 = make_cuts_independent(pattern.get_pwd2()); Piecewise x0 = vertical_pattern.get_value() ? Piecewise(patternd2[1]) : Piecewise(patternd2[0]); Piecewise y0 = vertical_pattern.get_value() ? Piecewise(patternd2[0]) : Piecewise(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 > cons // remove this once we have unified coordinate systems path_a_pw = path_a_pw + Geom::Point(offsetx, -offsety); - D2 > B = make_cuts_independant(path_a_pw); + D2 > B = make_cuts_independent(path_a_pw); Piecewise 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 > d2pw = make_cuts_independant(pwd2); + D2 > 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); -- 2.30.2