![[tokkee]](http://tokkee.org/images/avatar.png){kind=avatar}
summary | shortlog | log | commit | commitdiff | tree
raw | patch | inline | side by side (parent: aecee8c)
raw | patch | inline | side by side (parent: aecee8c)
| author | johanengelen <johanengelen@users.sourceforge.net> | |
| Thu, 3 Jul 2008 20:06:40 +0000 (20:06 +0000) | ||
| committer | johanengelen <johanengelen@users.sourceforge.net> | |
| Thu, 3 Jul 2008 20:06:40 +0000 (20:06 +0000) |
100 files changed:
index 97c4c6e5ce78d92a4bc2a5dfb123b299c3cc7eae..5375e5b58296c80009e83701c25bf009c9a5bb53 100644 (file)
-#include "basic-intersection.h"
-#include "exception.h"
+#include <2geom/basic-intersection.h>
+#include <2geom/exception.h>
unsigned intersect_steps = 0;
}
std::vector<std::pair<double, double> >
-find_self_intersections(OldBezier const &Sb) {
+find_self_intersections(OldBezier const &/*Sb*/) {
THROW_NOTIMPLEMENTED();
}
index 76abcce2a7055ab93e10ccf778b2e515877731ac..0090b0305e72e0675401000bfb1f2f8cbfa08891 100644 (file)
-#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 {
index ceab3cc111c655257f3bce18a17c7b27dd02c56c..244328f8af7a8ec80b5592021063d709e5fed590 100644 (file)
--- a/src/2geom/bezier-curve.h
+++ b/src/2geom/bezier-curve.h
#define _2GEOM_BEZIER_CURVE_H_
-#include "curve.h"
-#include "sbasis-curve.h" // for non-native winding method
-#include "bezier.h"
+#include <2geom/curve.h>
+#include <2geom/sbasis-curve.h> // for non-native winding method
+#include <2geom/bezier.h>
#include <algorithm>
index f88913584b3d54e7830c5bb1f0dbded05d983e4e..8b421f2e7e872299ce7c2b246f11a243f51f6892 100644 (file)
#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{
index 59aac895117fe4ab53665636dc655be62a46320d..83767af985008bead21fa00bd5f9c50a4ad6db09 100644 (file)
# include <ieefp.h>
#endif
-#include "bezier-utils.h"
+#include <2geom/bezier-utils.h>
-#include "isnan.h"
+#include <2geom/isnan.h>
#include <assert.h>
namespace Geom{
index fba4333ff989c9ab5e80c23bed71e949cb181a05..68ae8c0e7fe73fb77a608587b8588d9212831a09 100644 (file)
--- a/src/2geom/bezier-utils.h
+++ b/src/2geom/bezier-utils.h
*
*/
-#include "point.h"
+#include <2geom/point.h>
namespace Geom{
diff --git a/src/2geom/bezier.h b/src/2geom/bezier.h
index bc9d6032e63fe0c8befdacdb2ce28e0f7764be28..79b15cd0364259270c4fcc2c1da3c78e5e983992 100644 (file)
--- a/src/2geom/bezier.h
+++ b/src/2geom/bezier.h
#ifndef SEEN_BEZIER_H
#define SEEN_BEZIER_H
-#include "coord.h"
+#include <2geom/coord.h>
#include <valarray>
-#include "isnan.h"
-#include "bezier-to-sbasis.h"
-#include "d2.h"
-#include "solver.h"
+#include <2geom/isnan.h>
+#include <2geom/bezier-to-sbasis.h>
+#include <2geom/d2.h>
+#include <2geom/solver.h>
#include <boost/optional/optional.hpp>
namespace Geom {
diff --git a/src/2geom/choose.h b/src/2geom/choose.h
index f3a8b0f4fcda87fd796c06080bd179cf62be3031..bfd66e8a9091b233452a9917c7a89dc3dd8b2101 100644 (file)
--- a/src/2geom/choose.h
+++ b/src/2geom/choose.h
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);
index 024864091d0c42f8c0a27bbdec35bdaf29886065..25385180b531e8f2b2cd1f1d0b8fe791cc05efda 100644 (file)
*/
#include <stdio.h>
#include <math.h>
-#include "point.h"
+#include <2geom/point.h>
namespace Geom{
diff --git a/src/2geom/concepts.h b/src/2geom/concepts.h
index 3b6fd3577bd800d877c1ed6b38ae2dca33b5fa71..6f10c8bb0c31ac6dbdec82f85213dd85b94ff459 100644 (file)
--- a/src/2geom/concepts.h
+++ b/src/2geom/concepts.h
#ifndef SEEN_CONCEPTS_H
#define SEEN_CONCEPTS_H
-#include "sbasis.h"
-#include "interval.h"
-#include "point.h"
+#include <2geom/sbasis.h>
+#include <2geom/interval.h>
+#include <2geom/point.h>
#include <vector>
#include <boost/concept_check.hpp>
index b98bb314c6c56aa9ce397e90fbfd3c9a2646f750..f5a0f9cd82dffcd00b0c3bebb002a280ab8e1c7b 100644 (file)
#include <stdlib.h>
#include <valarray>
#include <cassert>
-#include "conjugate_gradient.h"
+#include <2geom/conjugate_gradient.h>
/* lifted wholely from wikipedia. */
}
}
+/**
+// only used in commented code below
static double Linfty(valarray<double> const &vec) {
return std::max(vec.max(), -vec.min());
}
+**/
double
inner(valarray<double> const &x,
valarray<double> &x,
valarray<double> const &b,
unsigned n, double tol,
- unsigned max_iterations, bool ortho1) {
+ unsigned max_iterations, bool /*ortho1*/) {
valarray<double> Ap(n), p(n), r(n);
matrix_times_vector(A,x,Ap);
r=b-Ap;
index 50c43e6d144de46e67904a21b7cf27156533dca5..7127a7c097ec5a6132d3a4f446f57716b3dcf6d0 100644 (file)
*
*/
-#include "convex-cover.h"
+#include <2geom/convex-cover.h>
#include <algorithm>
#include <map>
/** Todo:
*/
bool
ConvexHull::no_colinear_points() const {
-
+ // XXX: implement me!
}
bool
map<int, int> abridges;
map<int, int> bbridges;
- for(int ia = 0; ia < a.boundary.size(); ia++) {
- for(int ib = 0; ib < b.boundary.size(); ib++) {
+ for(unsigned ia = 0; ia < a.boundary.size(); ia++) {
+ for(unsigned ib = 0; ib < b.boundary.size(); ib++) {
Point d = b[ib] - a[ia];
Geom::Coord e = cross(d, a[ia - 1] - a[ia]), f = cross(d, a[ia + 1] - a[ia]);
Geom::Coord g = cross(d, b[ib - 1] - a[ia]), h = cross(d, b[ib + 1] - a[ia]);
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])) {
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;
}
* (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;
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];
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];
index d99e07b950d41c64802566dfedb034a463dd586b..1fe86e40d85b4215f7ba60ae70d536b4a91133c4 100644 (file)
--- a/src/2geom/convex-cover.h
+++ b/src/2geom/convex-cover.h
* convex hull class is included here (the convex-hull header is wrong)
*/
-#include "point.h"
+#include <2geom/point.h>
#include <vector>
namespace Geom{
diff --git a/src/2geom/crossing.cpp b/src/2geom/crossing.cpp
index 880b99e1a390010aafe40a1c1a49c3223fa5f975..a4d4f6ab1ac2e95afef6405a4ca022671ed74c90 100644 (file)
--- a/src/2geom/crossing.cpp
+++ b/src/2geom/crossing.cpp
-#include "crossing.h"
-#include "path.h"
+#include <2geom/crossing.h>
+#include <2geom/path.h>
namespace Geom {
@@ -169,7 +169,7 @@ CrossingSet reverse_tb(CrossingSet const &cr, unsigned split, std::vector<double
return ret;
}
-void clean(Crossings &cr_a, Crossings &cr_b) {
+void clean(Crossings &/*cr_a*/, Crossings &/*cr_b*/) {
/* if(cr_a.empty()) return;
//Remove anything with dupes
diff --git a/src/2geom/crossing.h b/src/2geom/crossing.h
index 8e3a45142b9b86054c4dd06b6606f5985bf02b47..546e33ebd41df747051a91379d3c11f83eee8d1e 100644 (file)
--- a/src/2geom/crossing.h
+++ b/src/2geom/crossing.h
#define __GEOM_CROSSING_H
#include <vector>
-#include "rect.h"
-#include "sweep.h"
+#include <2geom/rect.h>
+#include <2geom/sweep.h>
namespace Geom {
index c0e46bc0675988cd9769f4a71e19ed7541c1f2bc..c767af54f9e2f1ef28751e5edf4b38b6ecc3ff81 100644 (file)
*/
-#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 22d2ec55609cc4194b8374a1b0b10388419a4ba2..8eaf73fcf4d733e52bf31fef697783203b435cb1 100644 (file)
--- a/src/2geom/curve.h
+++ b/src/2geom/curve.h
#define _2GEOM_CURVE_H_
-#include "coord.h"
-#include "point.h"
-#include "interval.h"
-#include "nearest-point.h"
-#include "sbasis.h"
-#include "d2.h"
-#include "matrix.h"
-#include "exception.h"
+#include <2geom/coord.h>
+#include <2geom/point.h>
+#include <2geom/interval.h>
+#include <2geom/nearest-point.h>
+#include <2geom/sbasis.h>
+#include <2geom/d2.h>
+#include <2geom/matrix.h>
+#include <2geom/exception.h>
#include <vector>
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(); }
};
virtual D2<SBasis> toSBasis() const = 0;
+ virtual bool operator==(Curve const &c) const { return this == &c;}
};
inline
diff --git a/src/2geom/curves.h b/src/2geom/curves.h
index a4065930a558db880b34fd5860cd46564e0bc127..3bf7d9b598f37b4384793604a5815c70e8960805 100644 (file)
--- a/src/2geom/curves.h
+++ b/src/2geom/curves.h
* Authors:
* MenTaLguY <mental@rydia.net>
* Marco Cecchetti <mrcekets at gmail.com>
- *
+ *
* Copyright 2007-2008 authors
*
* This library is free software; you can redistribute it and/or
#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_
index 9b6ca269ca9fc0329d021d003490dd1fa7e0a20d..01f83bf5ecf86e26e7320d5a232c488b74c153c5 100644 (file)
--- a/src/2geom/d2-sbasis.cpp
+++ b/src/2geom/d2-sbasis.cpp
-#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. */
return ret;
}
-D2<Piecewise<SBasis> > make_cuts_independant(Piecewise<D2<SBasis> > const &a) {
+D2<Piecewise<SBasis> > make_cuts_independent(Piecewise<D2<SBasis> > const &a) {
D2<Piecewise<SBasis> > ret;
for(unsigned d = 0; d < 2; d++) {
for(unsigned i = 0; i < a.size(); i++)
diff --git a/src/2geom/d2-sbasis.h b/src/2geom/d2-sbasis.h
index e921896f5e0e433f65275a22ad3b461e224b52b0..91ceb31c98c916eb9f73809f068a67962947507f 100644 (file)
--- a/src/2geom/d2-sbasis.h
+++ b/src/2geom/d2-sbasis.h
#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
//Piecewise<D2<SBasis> > specific decls:
Piecewise<D2<SBasis> > sectionize(D2<Piecewise<SBasis> > const &a);
-D2<Piecewise<SBasis> > make_cuts_independant(Piecewise<D2<SBasis> > const &a);
+D2<Piecewise<SBasis> > make_cuts_independent(Piecewise<D2<SBasis> > const &a);
Piecewise<D2<SBasis> > rot90(Piecewise<D2<SBasis> > const &a);
Piecewise<SBasis> dot(Piecewise<D2<SBasis> > const &a, Piecewise<D2<SBasis> > const &b);
Piecewise<SBasis> cross(Piecewise<D2<SBasis> > const &a, Piecewise<D2<SBasis> > const &b);
diff --git a/src/2geom/d2.h b/src/2geom/d2.h
index db8cf68c4e503ca012a0f6c0c00edc30ca311e74..731a084a1d495dc9a99cb2041f70c8dcefa0eec8 100644 (file)
--- a/src/2geom/d2.h
+++ b/src/2geom/d2.h
#ifndef _2GEOM_D2 //If this is change, change the guard in rect.h as well.
#define _2GEOM_D2
-#include "point.h"
-#include "interval.h"
-#include "matrix.h"
+#include <2geom/point.h>
+#include <2geom/interval.h>
+#include <2geom/matrix.h>
#include <boost/concept_check.hpp>
-#include "concepts.h"
+#include <2geom/concepts.h>
namespace Geom{
} //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
--- /dev/null
+++ b/src/2geom/ellipse.cpp
@@ -0,0 +1,206 @@
+/*\r
+ * Ellipse Curve\r
+ *\r
+ * Authors:\r
+ * Marco Cecchetti <mrcekets at gmail.com>\r
+ *\r
+ * Copyright 2008 authors\r
+ *\r
+ * This library is free software; you can redistribute it and/or\r
+ * modify it either under the terms of the GNU Lesser General Public\r
+ * License version 2.1 as published by the Free Software Foundation\r
+ * (the "LGPL") or, at your option, under the terms of the Mozilla\r
+ * Public License Version 1.1 (the "MPL"). If you do not alter this\r
+ * notice, a recipient may use your version of this file under either\r
+ * the MPL or the LGPL.\r
+ *\r
+ * You should have received a copy of the LGPL along with this library\r
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software\r
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
+ * You should have received a copy of the MPL along with this library\r
+ * in the file COPYING-MPL-1.1\r
+ *\r
+ * The contents of this file are subject to the Mozilla Public License\r
+ * Version 1.1 (the "License"); you may not use this file except in\r
+ * compliance with the License. You may obtain a copy of the License at\r
+ * http://www.mozilla.org/MPL/\r
+ *\r
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY\r
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for\r
+ * the specific language governing rights and limitations.\r
+ */\r
+\r
+\r
+#include <2geom/ellipse.h>\r
+#include <2geom/svg-elliptical-arc.h>\r
+#include <2geom/numeric/fitting-tool.h>\r
+#include <2geom/numeric/fitting-model.h>\r
+\r
+\r
+namespace Geom\r
+{\r
+\r
+void Ellipse::set(double A, double B, double C, double D, double E, double F)\r
+{\r
+ double den = 4*A*C - B*B;\r
+ if ( den == 0 )\r
+ {\r
+ THROW_LOGICALERROR("den == 0, while computing ellipse centre");\r
+ }\r
+ m_centre[X] = (B*E - 2*C*D) / den;\r
+ m_centre[Y] = (B*D - 2*A*E) / den;\r
+\r
+ // evaluate the a coefficient of the ellipse equation in normal form\r
+ // E(x,y) = a*(x-cx)^2 + b*(x-cx)*(y-cy) + c*(y-cy)^2 = 1\r
+ // where b = a*B , c = a*C, (cx,cy) == centre\r
+ double num = A * sqr(m_centre[X])\r
+ + B * m_centre[X] * m_centre[Y]\r
+ + C * sqr(m_centre[Y])\r
+ - A * F;\r
+\r
+\r
+ //evaluate ellipse rotation angle\r
+ double rot = std::atan2( -B, -(A - C) )/2;\r
+// std::cerr << "rot = " << rot << std::endl;\r
+ bool swap_axes = false;\r
+ if ( are_near(rot, 0) ) rot = 0;\r
+ if ( are_near(rot, M_PI/2) || rot < 0 )\r
+ {\r
+ swap_axes = true;\r
+ }\r
+\r
+ // evaluate the length of the ellipse rays\r
+ double cosrot = std::cos(rot);\r
+ double sinrot = std::sin(rot);\r
+ double cos2 = cosrot * cosrot;\r
+ double sin2 = sinrot * sinrot;\r
+ double cossin = cosrot * sinrot;\r
+\r
+ den = A * cos2 + B * cossin + C * sin2;\r
+ if ( den == 0 )\r
+ {\r
+ THROW_LOGICALERROR("den == 0, while computing 'rx' coefficient");\r
+ }\r
+ double rx2 = num/den;\r
+ if ( rx2 < 0 )\r
+ {\r
+ THROW_LOGICALERROR("rx2 < 0, while computing 'rx' coefficient");\r
+ }\r
+ double rx = std::sqrt(rx2);\r
+\r
+ den = C * cos2 - B * cossin + A * sin2;\r
+ if ( den == 0 )\r
+ {\r
+ THROW_LOGICALERROR("den == 0, while computing 'ry' coefficient");\r
+ }\r
+ double ry2 = num/den;\r
+ if ( ry2 < 0 )\r
+ {\r
+ THROW_LOGICALERROR("ry2 < 0, while computing 'rx' coefficient");\r
+ }\r
+ double ry = std::sqrt(ry2);\r
+\r
+ // the solution is not unique so we choose always the ellipse\r
+ // with a rotation angle between 0 and PI/2\r
+ if ( swap_axes ) std::swap(rx, ry);\r
+ if ( are_near(rot, M_PI/2)\r
+ || are_near(rot, -M_PI/2)\r
+ || are_near(rx, ry) )\r
+ {\r
+ rot = 0;\r
+ }\r
+ else if ( rot < 0 )\r
+ {\r
+ rot += M_PI/2;\r
+ }\r
+\r
+ m_ray[X] = rx;\r
+ m_ray[Y] = ry;\r
+ m_angle = rot;\r
+}\r
+\r
+\r
+void Ellipse::set(std::vector<Point> const& points)\r
+{\r
+ size_t sz = points.size();\r
+ if (sz < 5)\r
+ {\r
+ THROW_RANGEERROR("fitting error: too few points passed");\r
+ }\r
+ NL::LFMEllipse model;\r
+ NL::least_squeares_fitter<NL::LFMEllipse> fitter(model, sz);\r
+\r
+ for (size_t i = 0; i < sz; ++i)\r
+ {\r
+ fitter.append(points[i]);\r
+ }\r
+ fitter.update();\r
+\r
+ NL::Vector z(sz, 0.0);\r
+ model.instance(*this, fitter.result(z));\r
+}\r
+\r
+\r
+SVGEllipticalArc\r
+Ellipse::arc(Point const& initial, Point const& inner, Point const& final,\r
+ bool _svg_compliant)\r
+{\r
+ Point sp_cp = initial - center();\r
+ Point ep_cp = final - center();\r
+ Point ip_cp = inner - center();\r
+\r
+ double angle1 = angle_between(sp_cp, ep_cp);\r
+ double angle2 = angle_between(sp_cp, ip_cp);\r
+ double angle3 = angle_between(ip_cp, ep_cp);\r
+\r
+ bool large_arc_flag = true;\r
+ bool sweep_flag = true;\r
+\r
+ if ( angle1 > 0 )\r
+ {\r
+ if ( angle2 > 0 && angle3 > 0 )\r
+ {\r
+ large_arc_flag = false;\r
+ sweep_flag = true;\r
+ }\r
+ else\r
+ {\r
+ large_arc_flag = true;\r
+ sweep_flag = false;\r
+ }\r
+ }\r
+ else\r
+ {\r
+ if ( angle2 < 0 && angle3 < 0 )\r
+ {\r
+ large_arc_flag = false;\r
+ sweep_flag = false;\r
+ }\r
+ else\r
+ {\r
+ large_arc_flag = true;\r
+ sweep_flag = true;\r
+ }\r
+ }\r
+\r
+ SVGEllipticalArc ea( initial, ray(X), ray(Y), rot_angle(),\r
+ large_arc_flag, sweep_flag, final, _svg_compliant);\r
+ return ea;\r
+}\r
+\r
+\r
+} // end namespace Geom\r
+\r
+\r
+/*\r
+ Local Variables:\r
+ mode:c++\r
+ c-file-style:"stroustrup"\r
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))\r
+ indent-tabs-mode:nil\r
+ fill-column:99\r
+ End:\r
+*/\r
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :\r
+\r
+\r
diff --git a/src/2geom/ellipse.h b/src/2geom/ellipse.h
--- /dev/null
+++ b/src/2geom/ellipse.h
@@ -0,0 +1,133 @@
+/*\r
+ * Ellipse Curve\r
+ *\r
+ * Authors:\r
+ * Marco Cecchetti <mrcekets at gmail.com>\r
+ *\r
+ * Copyright 2008 authors\r
+ *\r
+ * This library is free software; you can redistribute it and/or\r
+ * modify it either under the terms of the GNU Lesser General Public\r
+ * License version 2.1 as published by the Free Software Foundation\r
+ * (the "LGPL") or, at your option, under the terms of the Mozilla\r
+ * Public License Version 1.1 (the "MPL"). If you do not alter this\r
+ * notice, a recipient may use your version of this file under either\r
+ * the MPL or the LGPL.\r
+ *\r
+ * You should have received a copy of the LGPL along with this library\r
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software\r
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
+ * You should have received a copy of the MPL along with this library\r
+ * in the file COPYING-MPL-1.1\r
+ *\r
+ * The contents of this file are subject to the Mozilla Public License\r
+ * Version 1.1 (the "License"); you may not use this file except in\r
+ * compliance with the License. You may obtain a copy of the License at\r
+ * http://www.mozilla.org/MPL/\r
+ *\r
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY\r
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for\r
+ * the specific language governing rights and limitations.\r
+ */\r
+\r
+\r
+#ifndef _2GEOM_ELLIPSE_H_\r
+#define _2GEOM_ELLIPSE_H_\r
+\r
+\r
+#include <2geom/point.h>\r
+#include <2geom/exception.h>\r
+\r
+\r
+namespace Geom\r
+{\r
+\r
+class SVGEllipticalArc;\r
+\r
+class Ellipse\r
+{\r
+ public:\r
+ Ellipse()\r
+ {}\r
+\r
+ Ellipse(double cx, double cy, double rx, double ry, double a)\r
+ : m_centre(cx, cy), m_ray(rx, ry), m_angle(a)\r
+ {\r
+ }\r
+\r
+ // build an ellipse by its implicit equation:\r
+ // Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0\r
+ Ellipse(double A, double B, double C, double D, double E, double F)\r
+ {\r
+ set(A, B, C, D, E, F);\r
+ }\r
+\r
+ Ellipse(std::vector<Point> const& points)\r
+ {\r
+ set(points);\r
+ }\r
+\r
+ void set(double cx, double cy, double rx, double ry, double a)\r
+ {\r
+ m_centre[X] = cx;\r
+ m_centre[Y] = cy;\r
+ m_ray[X] = rx;\r
+ m_ray[Y] = ry;\r
+ m_angle = a;\r
+ }\r
+\r
+ // build an ellipse by its implicit equation:\r
+ // Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0\r
+ void set(double A, double B, double C, double D, double E, double F);\r
+\r
+ // biuld up the best fitting ellipse wrt the passed points\r
+ // prerequisite: at least 5 points must be passed\r
+ void set(std::vector<Point> const& points);\r
+\r
+ SVGEllipticalArc\r
+ arc(Point const& initial, Point const& inner, Point const& final,\r
+ bool _svg_compliant = true);\r
+\r
+ Point center() const\r
+ {\r
+ return m_centre;\r
+ }\r
+\r
+ Coord center(Dim2 d) const\r
+ {\r
+ return m_centre[d];\r
+ }\r
+\r
+ Coord ray(Dim2 d) const\r
+ {\r
+ return m_ray[d];\r
+ }\r
+\r
+ Coord rot_angle() const\r
+ {\r
+ return m_angle;\r
+ }\r
+\r
+ private:\r
+ Point m_centre, m_ray;\r
+ double m_angle;\r
+};\r
+\r
+\r
+} // end namespace Geom\r
+\r
+\r
+\r
+#endif // _2GEOM_ELLIPSE_H_\r
+\r
+\r
+/*\r
+ Local Variables:\r
+ mode:c++\r
+ c-file-style:"stroustrup"\r
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))\r
+ indent-tabs-mode:nil\r
+ fill-column:99\r
+ End:\r
+*/\r
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :\r
index b30dfde5d11d85da06163951a6984962d6f3a8a5..552cdf5b3188f23c0b8c3e53ebb32b2cf26b252b 100644 (file)
*/
-#include "elliptical-arc.h"
-#include "bezier-curve.h"
-#include "poly.h"
+#include <2geom/elliptical-arc.h>
+#include <2geom/bezier-curve.h>
+#include <2geom/poly.h>
#include <cfloat>
#include <limits>
index 3a3b90670199cdcda2db4bab342799d532cf280f..46b94c0df453dd13a07fe27835a4fface428a8c0 100644 (file)
#define _2GEOM_ELLIPTICAL_ARC_H_
-#include "curve.h"
-#include "angle.h"
-#include "utils.h"
-#include "sbasis-curve.h" // for non-native methods
+#include <2geom/curve.h>
+#include <2geom/angle.h>
+#include <2geom/utils.h>
+#include <2geom/sbasis-curve.h> // for non-native methods
#include <algorithm>
diff --git a/src/2geom/geom.cpp b/src/2geom/geom.cpp
index 8845cefa89fe443ae321fe022a7ee4be8c22c574..64bf21ee6b4b20bd48dce366d4db0ddbbe2c6e8e 100644 (file)
--- a/src/2geom/geom.cpp
+++ b/src/2geom/geom.cpp
#ifdef HAVE_CONFIG_H
# include <config.h>
#endif
-#include "geom.h"
-#include "point.h"
+#include <2geom/geom.h>
+#include <2geom/point.h>
namespace Geom {
diff --git a/src/2geom/geom.h b/src/2geom/geom.h
index 49156482c990ef54a25033b4d421857837b05514..207d609b09b8568c888d0966e6a16b9b722680a5 100644 (file)
--- a/src/2geom/geom.h
+++ b/src/2geom/geom.h
//TODO: move somewhere else
#include <vector>
-#include "point.h"
+#include <2geom/point.h>
namespace Geom {
index 732c299387efdac185bed092da9df3de9284d7b2..a34c5a962dc0060dac8492c0b0d0094ba3be0711 100644 (file)
#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 eb506dc1f060fff14f730e1d00de49e381a112e0..f042294ff605f1f276ac0ad40d1b74e043faf617 100644 (file)
--- a/src/2geom/interval.h
+++ b/src/2geom/interval.h
#define SEEN_INTERVAL_H
#include <assert.h>
-#include "coord.h"
+#include <2geom/coord.h>
#include <boost/optional/optional.hpp>
diff --git a/src/2geom/isnan.h b/src/2geom/isnan.h
index d95a45f10d721d655d4a4e958a67ee321858b9fe..6b94daa6e5c61826784c2799af54ee1226a37df5 100644 (file)
--- a/src/2geom/isnan.h
+++ b/src/2geom/isnan.h
# 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))
# define IS_FINITE(_a) (isfinite(_a))
#elif defined(__osf__)
# define IS_FINITE(_a) (finite(_a) && !IS_NAN(_a))
-#elif defined (SOLARIS)
+#elif defined (SOLARIS_2_8) && __GNUC__ == 3 && __GNUC_MINOR__ == 2
#include <ieeefp.h>
#define IS_FINITE(_a) (finite(_a) && !IS_NAN(_a))
#else
diff --git a/src/2geom/linear.h b/src/2geom/linear.h
index 271e87be4642a4a29d84483b95c4c07146e75220..bc7af564c9ad5dacf1037378bb3f88619f11150a 100644 (file)
--- a/src/2geom/linear.h
+++ b/src/2geom/linear.h
#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 f90bb6d420543a5c9f17b77ec663d60226b446a2..d25cc8f7e2151c313072ba1fda5793e31e5e5719 100644 (file)
--- a/src/2geom/matrix.cpp
+++ b/src/2geom/matrix.cpp
* 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 ba4451265a25275eb62453800fc2dfe46182b218..dcf765c50ae9dde84cc1ab1911049c0be81bd8bc 100644 (file)
--- a/src/2geom/matrix.h
+++ b/src/2geom/matrix.h
//#include <glib/gmessages.h>
-#include "point.h"
+#include <2geom/point.h>
namespace Geom {
index 074de1cb3364c54780731ffb0588c23e33eb9be3..59d55ba326383f5669f0a65b1aa323a83ebf871b 100644 (file)
*/\r
\r
\r
-#include "nearest-point.h"\r
+#include <2geom/nearest-point.h>\r
\r
namespace Geom\r
{\r
\r
std::vector<double> \r
all_nearest_points( Point const& p, \r
- D2<SBasis> const& c, \r
- D2<SBasis> const& dc, \r
- double from, double to )\r
+ D2<SBasis> const& c, \r
+ D2<SBasis> const& /*dc*/, \r
+ double from, double to )\r
{\r
std::swap(from, to);\r
if ( from > to ) std::swap(from, to);\r
index 73ac0c3ce73f9785a68d19d3afb9f83bdf310f79..0b43ce67b51ee195cce8b5edb4dc266dc2bcbab4 100644 (file)
\r
#include <vector>\r
\r
-#include "d2.h"\r
-#include "piecewise.h"\r
-#include "exception.h"\r
+#include <2geom/d2.h>\r
+#include <2geom/piecewise.h>\r
+#include <2geom/exception.h>\r
\r
\r
\r
diff --git a/src/2geom/numeric/fitting-model.h b/src/2geom/numeric/fitting-model.h
--- /dev/null
@@ -0,0 +1,423 @@
+/*\r
+ * Fitting Models for Geom Types\r
+ *\r
+ * Authors:\r
+ * Marco Cecchetti <mrcekets at gmail.com>\r
+ *\r
+ * Copyright 2008 authors\r
+ *\r
+ * This library is free software; you can redistribute it and/or\r
+ * modify it either under the terms of the GNU Lesser General Public\r
+ * License version 2.1 as published by the Free Software Foundation\r
+ * (the "LGPL") or, at your option, under the terms of the Mozilla\r
+ * Public License Version 1.1 (the "MPL"). If you do not alter this\r
+ * notice, a recipient may use your version of this file under either\r
+ * the MPL or the LGPL.\r
+ *\r
+ * You should have received a copy of the LGPL along with this library\r
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software\r
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
+ * You should have received a copy of the MPL along with this library\r
+ * in the file COPYING-MPL-1.1\r
+ *\r
+ * The contents of this file are subject to the Mozilla Public License\r
+ * Version 1.1 (the "License"); you may not use this file except in\r
+ * compliance with the License. You may obtain a copy of the License at\r
+ * http://www.mozilla.org/MPL/\r
+ *\r
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY\r
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for\r
+ * the specific language governing rights and limitations.\r
+ */\r
+\r
+\r
+#ifndef _NL_FITTING_MODEL_H_\r
+#define _NL_FITTING_MODEL_H_\r
+\r
+\r
+#include <2geom/d2.h>\r
+#include <2geom/sbasis.h>\r
+#include <2geom/bezier.h>\r
+#include <2geom/bezier-curve.h>\r
+#include <2geom/poly.h>\r
+#include <2geom/ellipse.h>\r
+#include <2geom/utils.h>\r
+\r
+\r
+namespace Geom { namespace NL {\r
+\r
+\r
+/*\r
+ * completely unknown models must inherit from this template class;\r
+ * example: the model a*x^2 + b*x + c = 0 to be solved wrt a, b, c;\r
+ * example: the model A(t) = known_sample_value_at(t) to be solved wrt\r
+ * the coefficients of the curve A(t) expressed in S-Basis form;\r
+ * parameter type: the type of x and t variable in the examples above;\r
+ * value type: the type of the known sample values (in the first example\r
+ * is constant )\r
+ * instance type: the type of the objects produced by using\r
+ * the fitting raw data solution\r
+ */\r
+template< typename ParameterType, typename ValueType, typename InstanceType >\r
+class LinearFittingModel\r
+{\r
+ public:\r
+ typedef ParameterType parameter_type;\r
+ typedef ValueType value_type;\r
+ typedef InstanceType instance_type;\r
+\r
+ static const bool WITH_FIXED_TERMS = false;\r
+\r
+ /*\r
+ * a LinearFittingModel must implement the following methods:\r
+ *\r
+ * void feed( VectorView & vector,\r
+ * parameter_type const& sample_parameter ) const;\r
+ *\r
+ * size_t size() const;\r
+ *\r
+ * void instance(instance_type &, raw_type const& raw_data) const;\r
+ *\r
+ */\r
+};\r
+\r
+\r
+/*\r
+ * partially known models must inherit from this template class\r
+ * example: the model a*x^2 + 2*x + c = 0 to be solved wrt a and c\r
+ */\r
+template< typename ParameterType, typename ValueType, typename InstanceType >\r
+class LinearFittingModelWithFixedTerms\r
+{\r
+ public:\r
+ typedef ParameterType parameter_type;\r
+ typedef ValueType value_type;\r
+ typedef InstanceType instance_type;\r
+\r
+ static const bool WITH_FIXED_TERMS = true;\r
+\r
+ /*\r
+ * a LinearFittingModelWithFixedTerms must implement the following methods:\r
+ *\r
+ * void feed( VectorView & vector,\r
+ * value_type & fixed_term,\r
+ * parameter_type const& sample_parameter ) const;\r
+ *\r
+ * size_t size() const;\r
+ *\r
+ * void instance(instance_type &, raw_type const& raw_data) const;\r
+ *\r
+ */\r
+\r
+\r
+};\r
+\r
+\r
+// incomplete model, it can be inherited to make up different kinds of\r
+// instance type; the raw data is a vector of coefficients of a polynomial\r
+// rapresented in standard power basis\r
+template< typename InstanceType >\r
+class LFMPowerBasis\r
+ : public LinearFittingModel<double, double, InstanceType>\r
+{\r
+ public:\r
+ LFMPowerBasis(size_t degree)\r
+ : m_size(degree + 1)\r
+ {\r
+ }\r
+\r
+ void feed( VectorView & coeff, double sample_parameter ) const\r
+ {\r
+ coeff[0] = 1;\r
+ double x_i = 1;\r
+ for (size_t i = 1; i < coeff.size(); ++i)\r
+ {\r
+ x_i *= sample_parameter;\r
+ coeff[i] = x_i;\r
+ }\r
+ }\r
+\r
+ size_t size() const\r
+ {\r
+ return m_size;\r
+ }\r
+\r
+ private:\r
+ size_t m_size;\r
+};\r
+\r
+\r
+// this model generates Geom::Poly objects\r
+class LFMPoly\r
+ : public LFMPowerBasis<Poly>\r
+{\r
+ public:\r
+ LFMPoly(size_t degree)\r
+ : LFMPowerBasis<Poly>(degree)\r
+ {\r
+ }\r
+\r
+ void instance(Poly & poly, ConstVectorView const& raw_data) const\r
+ {\r
+ poly.clear();\r
+ poly.resize(size());\r
+ for (size_t i = 0; i < raw_data.size(); ++i)\r
+ {\r
+ poly[i] = raw_data[i];\r
+ }\r
+ }\r
+};\r
+\r
+\r
+// incomplete model, it can be inherited to make up different kinds of\r
+// instance type; the raw data is a vector of coefficients of a polynomial\r
+// rapresented in standard power basis with leading term coefficient equal to 1\r
+template< typename InstanceType >\r
+class LFMNormalizedPowerBasis\r
+ : public LinearFittingModelWithFixedTerms<double, double, InstanceType>\r
+{\r
+ public:\r
+ LFMNormalizedPowerBasis(size_t _degree)\r
+ : m_model( _degree - 1)\r
+ {\r
+ assert(_degree > 0);\r
+ }\r
+\r
+\r
+ void feed( VectorView & coeff,\r
+ double & known_term,\r
+ double sample_parameter ) const\r
+ {\r
+ m_model.feed(coeff, sample_parameter);\r
+ known_term = coeff[m_model.size()-1] * sample_parameter;\r
+ }\r
+\r
+ size_t size() const\r
+ {\r
+ return m_model.size();\r
+ }\r
+\r
+ private:\r
+ LFMPowerBasis<InstanceType> m_model;\r
+};\r
+\r
+\r
+// incomplete model, it can be inherited to make up different kinds of\r
+// instance type; the raw data is a vector of coefficients of the equation\r
+// of an ellipse curve\r
+template< typename InstanceType >\r
+class LFMEllipseEquation\r
+ : public LinearFittingModelWithFixedTerms<Point, double, InstanceType>\r
+{\r
+ public:\r
+ void feed( VectorView & coeff, double & fixed_term, Point const& p ) const\r
+ {\r
+ coeff[0] = p[X] * p[Y];\r
+ coeff[1] = p[Y] * p[Y];\r
+ coeff[2] = p[X];\r
+ coeff[3] = p[Y];\r
+ coeff[4] = 1;\r
+ fixed_term = p[X] * p[X];\r
+ }\r
+\r
+ size_t size() const\r
+ {\r
+ return 5;\r
+ }\r
+};\r
+\r
+\r
+// this model generates Ellipse curves\r
+class LFMEllipse\r
+ : public LFMEllipseEquation<Ellipse>\r
+{\r
+ public:\r
+ void instance(Ellipse & e, ConstVectorView const& coeff) const\r
+ {\r
+ e.set(1, coeff[0], coeff[1], coeff[2], coeff[3], coeff[4]);\r
+ }\r
+};\r
+\r
+\r
+// this model generates SBasis objects\r
+class LFMSBasis\r
+ : public LinearFittingModel<double, double, SBasis>\r
+{\r
+ public:\r
+ LFMSBasis( size_t _order )\r
+ : m_size( 2*(_order+1) ),\r
+ m_order(_order)\r
+ {\r
+ }\r
+\r
+ void feed( VectorView & coeff, double t ) const\r
+ {\r
+ double u0 = 1-t;\r
+ double u1 = t;\r
+ double s = u0 * u1;\r
+ coeff[0] = u0;\r
+ coeff[1] = u1;\r
+ for (size_t i = 2; i < size(); i+=2)\r
+ {\r
+ u0 *= s;\r
+ u1 *= s;\r
+ coeff[i] = u0;\r
+ coeff[i+1] = u1;\r
+ }\r
+ }\r
+\r
+ size_t size() const\r
+ {\r
+ return m_size;\r
+ }\r
+\r
+ void instance(SBasis & sb, ConstVectorView const& raw_data) const\r
+ {\r
+ sb.clear();\r
+ sb.resize(m_order+1);\r
+ for (unsigned int i = 0, k = 0; i < raw_data.size(); i+=2, ++k)\r
+ {\r
+ sb[k][0] = raw_data[i];\r
+ sb[k][1] = raw_data[i+1];\r
+ }\r
+ }\r
+\r
+ private:\r
+ size_t m_size;\r
+ size_t m_order;\r
+};\r
+\r
+\r
+// this model generates D2<SBasis> objects\r
+class LFMD2SBasis\r
+ : public LinearFittingModel< double, Point, D2<SBasis> >\r
+{\r
+ public:\r
+ LFMD2SBasis( size_t _order )\r
+ : mosb(_order)\r
+ {\r
+ }\r
+\r
+ void feed( VectorView & coeff, double t ) const\r
+ {\r
+ mosb.feed(coeff, t);\r
+ }\r
+\r
+ size_t size() const\r
+ {\r
+ return mosb.size();\r
+ }\r
+\r
+ void instance(D2<SBasis> & d2sb, ConstMatrixView const& raw_data) const\r
+ {\r
+ mosb.instance(d2sb[X], raw_data.column_const_view(X));\r
+ mosb.instance(d2sb[Y], raw_data.column_const_view(Y));\r
+ }\r
+\r
+ private:\r
+ LFMSBasis mosb;\r
+};\r
+\r
+\r
+// this model generates Bezier objects\r
+class LFMBezier\r
+ : public LinearFittingModel<double, double, Bezier>\r
+{\r
+ public:\r
+ LFMBezier( size_t _order )\r
+ : m_size(_order + 1),\r
+ m_order(_order)\r
+ {\r
+ binomial_coefficients(m_bc, m_order);\r
+ }\r
+\r
+ void feed( VectorView & coeff, double t ) const\r
+ {\r
+ double s = 1;\r
+ for (size_t i = 0; i < size(); ++i)\r
+ {\r
+ coeff[i] = s * m_bc[i];\r
+ s *= t;\r
+ }\r
+ double u = 1-t;\r
+ s = 1;\r
+ for (size_t i = size()-1; i > 0; --i)\r
+ {\r
+ coeff[i] *= s;\r
+ s *= u;\r
+ }\r
+ coeff[0] *= s;\r
+ }\r
+\r
+ size_t size() const\r
+ {\r
+ return m_size;\r
+ }\r
+\r
+ void instance(Bezier & b, ConstVectorView const& raw_data) const\r
+ {\r
+ assert(b.size() == raw_data.size());\r
+ for (unsigned int i = 0; i < raw_data.size(); ++i)\r
+ {\r
+ b[i] = raw_data[i];\r
+ }\r
+ }\r
+\r
+ private:\r
+ size_t m_size;\r
+ size_t m_order;\r
+ std::vector<size_t> m_bc;\r
+};\r
+\r
+\r
+// this model generates Bezier curves\r
+template< unsigned int N >\r
+class LFMBezierCurve\r
+ : public LinearFittingModel< double, Point, BezierCurve<N> >\r
+{\r
+ public:\r
+ LFMBezierCurve( size_t _order )\r
+ : mob(_order)\r
+ {\r
+ }\r
+\r
+ void feed( VectorView & coeff, double t ) const\r
+ {\r
+ mob.feed(coeff, t);\r
+ }\r
+\r
+ size_t size() const\r
+ {\r
+ return mob.size();\r
+ }\r
+\r
+ void instance(BezierCurve<N> & bc, ConstMatrixView const& raw_data) const\r
+ {\r
+ Bezier bx(size()-1);\r
+ Bezier by(size()-1);\r
+ mob.instance(bx, raw_data.column_const_view(X));\r
+ mob.instance(by, raw_data.column_const_view(Y));\r
+ bc = BezierCurve<N>(bx, by);\r
+ }\r
+\r
+ private:\r
+ LFMBezier mob;\r
+};\r
+\r
+} // end namespace NL\r
+} // end namespace Geom\r
+\r
+\r
+#endif // _NL_FITTING_MODEL_H_\r
+\r
+\r
+/*\r
+ Local Variables:\r
+ mode:c++\r
+ c-file-style:"stroustrup"\r
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))\r
+ indent-tabs-mode:nil\r
+ fill-column:99\r
+ End:\r
+*/\r
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :\r
diff --git a/src/2geom/numeric/fitting-tool.h b/src/2geom/numeric/fitting-tool.h
--- /dev/null
@@ -0,0 +1,532 @@
+/*\r
+ * Fitting Tools\r
+ *\r
+ * Authors:\r
+ * Marco Cecchetti <mrcekets at gmail.com>\r
+ *\r
+ * Copyright 2008 authors\r
+ *\r
+ * This library is free software; you can redistribute it and/or\r
+ * modify it either under the terms of the GNU Lesser General Public\r
+ * License version 2.1 as published by the Free Software Foundation\r
+ * (the "LGPL") or, at your option, under the terms of the Mozilla\r
+ * Public License Version 1.1 (the "MPL"). If you do not alter this\r
+ * notice, a recipient may use your version of this file under either\r
+ * the MPL or the LGPL.\r
+ *\r
+ * You should have received a copy of the LGPL along with this library\r
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software\r
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
+ * You should have received a copy of the MPL along with this library\r
+ * in the file COPYING-MPL-1.1\r
+ *\r
+ * The contents of this file are subject to the Mozilla Public License\r
+ * Version 1.1 (the "License"); you may not use this file except in\r
+ * compliance with the License. You may obtain a copy of the License at\r
+ * http://www.mozilla.org/MPL/\r
+ *\r
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY\r
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for\r
+ * the specific language governing rights and limitations.\r
+ */\r
+\r
+\r
+#ifndef _NL_FITTING_TOOL_H_\r
+#define _NL_FITTING_TOOL_H_\r
+\r
+\r
+#include <2geom/numeric/vector.h>\r
+#include <2geom/numeric/matrix.h>\r
+\r
+#include <2geom/point.h>\r
+\r
+#include <vector>\r
+\r
+\r
+namespace Geom { namespace NL {\r
+\r
+namespace detail {\r
+\r
+\r
+template< typename ModelT>\r
+class lsf_base\r
+{\r
+ public:\r
+ typedef ModelT model_type;\r
+ typedef typename model_type::parameter_type parameter_type;\r
+ typedef typename model_type::value_type value_type;\r
+\r
+ lsf_base( model_type const& _model, size_t forecasted_samples )\r
+ : m_model(_model),\r
+ m_total_samples(0),\r
+ m_matrix(forecasted_samples, m_model.size()),\r
+ m_psdinv_matrix(NULL)\r
+ {\r
+ }\r
+\r
+ // compute pseudo inverse\r
+ void update()\r
+ {\r
+ if (total_samples() == 0) return;\r
+ if (m_psdinv_matrix != NULL)\r
+ {\r
+ delete m_psdinv_matrix;\r
+ }\r
+ MatrixView mv(m_matrix, 0, 0, total_samples(), m_matrix.columns());\r
+ m_psdinv_matrix = new Matrix( pseudo_inverse(mv) );\r
+ assert(m_psdinv_matrix != NULL);\r
+ }\r
+\r
+ size_t total_samples() const\r
+ {\r
+ return m_total_samples;\r
+ }\r
+\r
+ bool is_full() const\r
+ {\r
+ return (total_samples() == m_matrix.rows());\r
+ }\r
+\r
+ void clear()\r
+ {\r
+ m_total_samples = 0;\r
+ }\r
+\r
+ virtual\r
+ ~lsf_base()\r
+ {\r
+ if (m_psdinv_matrix != NULL)\r
+ {\r
+ delete m_psdinv_matrix;\r
+ }\r
+ }\r
+\r
+ protected:\r
+ const model_type & m_model;\r
+ size_t m_total_samples;\r
+ Matrix m_matrix;\r
+ Matrix* m_psdinv_matrix;\r
+\r
+}; // end class lsf_base\r
+\r
+\r
+\r
+\r
+template< typename ModelT, typename ValueType = typename ModelT::value_type>\r
+class lsf_solution\r
+{\r
+};\r
+\r
+// a fitting process on samples with value of type double\r
+// produces a solution of type Vector\r
+template< typename ModelT>\r
+class lsf_solution<ModelT, double>\r
+ : public lsf_base<ModelT>\r
+{\r
+public:\r
+ typedef ModelT model_type;\r
+ typedef typename model_type::parameter_type parameter_type;\r
+ typedef typename model_type::value_type value_type;\r
+ typedef Vector solution_type;\r
+ typedef lsf_base<model_type> base_type;\r
+\r
+ using base_type::m_model;\r
+ using base_type::m_psdinv_matrix;\r
+ using base_type::total_samples;\r
+\r
+public:\r
+ lsf_solution<ModelT, double>( model_type const& _model,\r
+ size_t forecasted_samples )\r
+ : base_type(_model, forecasted_samples),\r
+ m_solution(_model.size())\r
+ {\r
+ }\r
+\r
+ template< typename VectorT >\r
+ solution_type& result(VectorT const& sample_values)\r
+ {\r
+ assert(sample_values.size() == total_samples());\r
+ ConstVectorView sv(sample_values);\r
+ m_solution = (*m_psdinv_matrix) * sv;\r
+ return m_solution;\r
+ }\r
+\r
+ // a comparison between old sample values and the new ones is performed\r
+ // in order to minimize computation\r
+ // prerequisite:\r
+ // old_sample_values.size() == new_sample_values.size()\r
+ // no update() call can be performed between two result invocations\r
+ template< typename VectorT >\r
+ solution_type& result( VectorT const& old_sample_values,\r
+ VectorT const& new_sample_values )\r
+ {\r
+ assert(old_sample_values.size() == total_samples());\r
+ assert(new_sample_values.size() == total_samples());\r
+ Vector diff(total_samples());\r
+ for (size_t i = 0; i < diff.size(); ++i)\r
+ {\r
+ diff[i] = new_sample_values[i] - old_sample_values[i];\r
+ }\r
+ Vector column(m_model.size());\r
+ Vector delta(m_model.size(), 0.0);\r
+ for (size_t i = 0; i < diff.size(); ++i)\r
+ {\r
+ if (diff[i] != 0)\r
+ {\r
+ column = m_psdinv_matrix->column_view(i);\r
+ column.scale(diff[i]);\r
+ delta += column;\r
+ }\r
+ }\r
+ m_solution += delta;\r
+ return m_solution;\r
+ }\r
+\r
+ solution_type& result()\r
+ {\r
+ return m_solution;\r
+ }\r
+\r
+private:\r
+ solution_type m_solution;\r
+\r
+}; // end class lsf_solution<ModelT, double>\r
+\r
+\r
+// a fitting process on samples with value of type Point\r
+// produces a solution of type Matrix (with 2 columns)\r
+template< typename ModelT>\r
+class lsf_solution<ModelT, Point>\r
+ : public lsf_base<ModelT>\r
+{\r
+public:\r
+ typedef ModelT model_type;\r
+ typedef typename model_type::parameter_type parameter_type;\r
+ typedef typename model_type::value_type value_type;\r
+ typedef Matrix solution_type;\r
+ typedef lsf_base<model_type> base_type;\r
+\r
+ using base_type::m_model;\r
+ using base_type::m_psdinv_matrix;\r
+ using base_type::total_samples;\r
+\r
+public:\r
+ lsf_solution<ModelT, Point>( model_type const& _model,\r
+ size_t forecasted_samples )\r
+ : base_type(_model, forecasted_samples),\r
+ m_solution(_model.size(), 2)\r
+ {\r
+ }\r
+\r
+ solution_type& result(std::vector<Point> const& sample_values)\r
+ {\r
+ assert(sample_values.size() == total_samples());\r
+ Matrix svm(total_samples(), 2);\r
+ for (size_t i = 0; i < total_samples(); ++i)\r
+ {\r
+ svm(i, X) = sample_values[i][X];\r
+ svm(i, Y) = sample_values[i][Y];\r
+ }\r
+ m_solution = (*m_psdinv_matrix) * svm;\r
+ return m_solution;\r
+ }\r
+\r
+ // a comparison between old sample values and the new ones is performed\r
+ // in order to minimize computation\r
+ // prerequisite:\r
+ // old_sample_values.size() == new_sample_values.size()\r
+ // no update() call can to be performed between two result invocations\r
+ solution_type& result( std::vector<Point> const& old_sample_values,\r
+ std::vector<Point> const& new_sample_values )\r
+ {\r
+ assert(old_sample_values.size() == total_samples());\r
+ assert(new_sample_values.size() == total_samples());\r
+ Matrix diff(total_samples(), 2);\r
+ for (size_t i = 0; i < total_samples(); ++i)\r
+ {\r
+ diff(i, X) = new_sample_values[i][X] - old_sample_values[i][X];\r
+ diff(i, Y) = new_sample_values[i][Y] - old_sample_values[i][Y];\r
+ }\r
+ Vector column(m_model.size());\r
+ Matrix delta(m_model.size(), 2, 0.0);\r
+ VectorView deltax = delta.column_view(X);\r
+ VectorView deltay = delta.column_view(Y);\r
+ for (size_t i = 0; i < total_samples(); ++i)\r
+ {\r
+ if (diff(i, X) != 0)\r
+ {\r
+ column = m_psdinv_matrix->column_view(i);\r
+ column.scale(diff(i, X));\r
+ deltax += column;\r
+ }\r
+ if (diff(i, Y) != 0)\r
+ {\r
+ column = m_psdinv_matrix->column_view(i);\r
+ column.scale(diff(i, Y));\r
+ deltay += column;\r
+ }\r
+ }\r
+ m_solution += delta;\r
+ return m_solution;\r
+ }\r
+\r
+ solution_type& result()\r
+ {\r
+ return m_solution;\r
+ }\r
+\r
+private:\r
+ solution_type m_solution;\r
+\r
+}; // end class lsf_solution<ModelT, Point>\r
+\r
+\r
+\r
+\r
+template< typename ModelT,\r
+ bool WITH_FIXED_TERMS = ModelT::WITH_FIXED_TERMS >\r
+class lsf_with_fixed_terms\r
+{\r
+};\r
+\r
+\r
+// fitting tool for completely unknown models\r
+template< typename ModelT>\r
+class lsf_with_fixed_terms<ModelT, false>\r
+ : public lsf_solution<ModelT>\r
+{\r
+ public:\r
+ typedef ModelT model_type;\r
+ typedef typename model_type::parameter_type parameter_type;\r
+ typedef typename model_type::value_type value_type;\r
+ typedef lsf_solution<model_type> base_type;\r
+ typedef typename base_type::solution_type solution_type;\r
+\r
+ using base_type::total_samples;\r
+ using base_type::is_full;\r
+ using base_type::m_matrix;\r
+ using base_type::m_total_samples;\r
+ using base_type::m_model;\r
+\r
+ public:\r
+ lsf_with_fixed_terms<ModelT, false>( model_type const& _model,\r
+ size_t forecasted_samples )\r
+ : base_type(_model, forecasted_samples)\r
+ {\r
+ }\r
+\r
+ void append(parameter_type const& sample_parameter)\r
+ {\r
+ assert(!is_full());\r
+ VectorView row = m_matrix.row_view(total_samples());\r
+ m_model.feed(row, sample_parameter);\r
+ ++m_total_samples;\r
+ }\r
+\r
+ void append_copy(size_t sample_index)\r
+ {\r
+ assert(!is_full());\r
+ assert(sample_index < total_samples());\r
+ VectorView dest_row = m_matrix.row_view(total_samples());\r
+ VectorView source_row = m_matrix.row_view(sample_index);\r
+ dest_row = source_row;\r
+ ++m_total_samples;\r
+ }\r
+\r
+}; // end class lsf_with_fixed_terms<ModelT, false>\r
+\r
+\r
+// fitting tool for partially known models\r
+template< typename ModelT>\r
+class lsf_with_fixed_terms<ModelT, true>\r
+ : public lsf_solution<ModelT>\r
+{\r
+ public:\r
+ typedef ModelT model_type;\r
+ typedef typename model_type::parameter_type parameter_type;\r
+ typedef typename model_type::value_type value_type;\r
+ typedef lsf_solution<model_type> base_type;\r
+ typedef typename base_type::solution_type solution_type;\r
+\r
+ using base_type::total_samples;\r
+ using base_type::is_full;\r
+ using base_type::m_matrix;\r
+ using base_type::m_total_samples;\r
+ using base_type::m_model;\r
+\r
+ public:\r
+ lsf_with_fixed_terms<ModelT, true>( model_type const& _model,\r
+ size_t forecasted_samples )\r
+ : base_type(_model, forecasted_samples),\r
+ m_vector(forecasted_samples),\r
+ m_vector_view(NULL)\r
+ {\r
+ }\r
+ void append(parameter_type const& sample_parameter)\r
+ {\r
+ assert(!is_full());\r
+ VectorView row = m_matrix.row_view(total_samples());\r
+ m_model.feed(row, m_vector[total_samples()], sample_parameter);\r
+ ++m_total_samples;\r
+ }\r
+\r
+ void append_copy(size_t sample_index)\r
+ {\r
+ assert(!is_full());\r
+ assert(sample_index < total_samples());\r
+ VectorView dest_row = m_matrix.row_view(total_samples());\r
+ VectorView source_row = m_matrix.row_view(sample_index);\r
+ dest_row = source_row;\r
+ m_vector[total_samples()] = m_vector[sample_index];\r
+ ++m_total_samples;\r
+ }\r
+\r
+ void update()\r
+ {\r
+ base_type::update();\r
+ if (total_samples() == 0) return;\r
+ if (m_vector_view != NULL)\r
+ {\r
+ delete m_vector_view;\r
+ }\r
+ m_vector_view = new VectorView(m_vector, base_type::total_samples());\r
+ assert(m_vector_view != NULL);\r
+ }\r
+\r
+ virtual\r
+ ~lsf_with_fixed_terms<model_type, true>()\r
+ {\r
+ if (m_vector_view != NULL)\r
+ {\r
+ delete m_vector_view;\r
+ }\r
+ }\r
+\r
+ protected:\r
+ Vector m_vector;\r
+ VectorView* m_vector_view;\r
+\r
+}; // end class lsf_with_fixed_terms<ModelT, true>\r
+\r
+\r
+} // end namespace detail\r
+\r
+\r
+\r
+\r
+template< typename ModelT,\r
+ typename ValueType = typename ModelT::value_type,\r
+ bool WITH_FIXED_TERMS = ModelT::WITH_FIXED_TERMS >\r
+class least_squeares_fitter\r
+{\r
+};\r
+\r
+\r
+template< typename ModelT, typename ValueType >\r
+class least_squeares_fitter<ModelT, ValueType, false>\r
+ : public detail::lsf_with_fixed_terms<ModelT>\r
+{\r
+ public:\r
+ typedef ModelT model_type;\r
+ typedef detail::lsf_with_fixed_terms<model_type> base_type;\r
+ typedef typename base_type::parameter_type parameter_type;\r
+ typedef typename base_type::value_type value_type;\r
+ typedef typename base_type::solution_type solution_type;\r
+\r
+ public:\r
+ least_squeares_fitter<ModelT, ValueType, false>( model_type const& _model,\r
+ size_t forecasted_samples )\r
+ : base_type(_model, forecasted_samples)\r
+ {\r
+ }\r
+}; // end class least_squeares_fitter<ModelT, ValueType, true>\r
+\r
+\r
+template< typename ModelT>\r
+class least_squeares_fitter<ModelT, double, true>\r
+ : public detail::lsf_with_fixed_terms<ModelT>\r
+{\r
+ public:\r
+ typedef ModelT model_type;\r
+ typedef detail::lsf_with_fixed_terms<model_type> base_type;\r
+ typedef typename base_type::parameter_type parameter_type;\r
+ typedef typename base_type::value_type value_type;\r
+ typedef typename base_type::solution_type solution_type;\r
+\r
+ using base_type::m_vector_view;\r
+ using base_type::result;\r
+\r
+ public:\r
+ least_squeares_fitter<ModelT, double, true>( model_type const& _model,\r
+ size_t forecasted_samples )\r
+ : base_type(_model, forecasted_samples)\r
+ {\r
+ }\r
+\r
+ template< typename VectorT >\r
+ solution_type& result(VectorT const& sample_values)\r
+ {\r
+ assert(sample_values.size() == m_vector_view->size());\r
+ Vector sv(sample_values.size());\r
+ for (size_t i = 0; i < sv.size(); ++i)\r
+ sv[i] = sample_values[i] - (*m_vector_view)[i];\r
+ return base_type::result(sv);\r
+ }\r
+\r
+}; // end class least_squeares_fitter<ModelT, double, true>\r
+\r
+\r
+template< typename ModelT>\r
+class least_squeares_fitter<ModelT, Point, true>\r
+ : public detail::lsf_with_fixed_terms<ModelT>\r
+{\r
+ public:\r
+ typedef ModelT model_type;\r
+ typedef detail::lsf_with_fixed_terms<model_type> base_type;\r
+ typedef typename base_type::parameter_type parameter_type;\r
+ typedef typename base_type::value_type value_type;\r
+ typedef typename base_type::solution_type solution_type;\r
+\r
+ using base_type::m_vector_view;\r
+ using base_type::result;\r
+\r
+ public:\r
+ least_squeares_fitter<ModelT, Point, true>( model_type const& _model,\r
+ size_t forecasted_samples )\r
+ : base_type(_model, forecasted_samples)\r
+ {\r
+ }\r
+\r
+ solution_type& result(std::vector<Point> const& sample_values)\r
+ {\r
+ assert(sample_values.size() == m_vector_view->size());\r
+ NL::Matrix sv(sample_values.size(), 2);\r
+ for (size_t i = 0; i < sample_values.size(); ++i)\r
+ {\r
+ sv(i, X) = sample_values[i][X] - (*m_vector_view)[i];\r
+ sv(i, Y) = sample_values[i][Y] - (*m_vector_view)[i];\r
+ }\r
+ return base_type::result(sv);\r
+ }\r
+\r
+}; // end class least_squeares_fitter<ModelT, Point, true>\r
+\r
+\r
+} // end namespace NL\r
+} // end namespace Geom\r
+\r
+\r
+\r
+#endif // _NL_FITTING_TOOL_H_\r
+\r
+\r
+/*\r
+ Local Variables:\r
+ mode:c++\r
+ c-file-style:"stroustrup"\r
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))\r
+ indent-tabs-mode:nil\r
+ fill-column:99\r
+ End:\r
+*/\r
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :\r
index 7699c5224ead519e34c4bbe9ea7b7fb0b8a864ec..5b516c9e6ba6bf5ff0554ee3297f466e4cf80c34 100644 (file)
-#ifndef _NL_LINEAR_SYSTEM_H_
-#define _NL_LINEAR_SYSTEM_H_
-
-
-#include <cassert>
-
-#include <gsl/gsl_linalg.h>
-
-#include "2geom/numeric/matrix.h"
-#include "2geom/numeric/vector.h"
-
-
-namespace Geom { namespace NL {
-
-
-class LinearSystem
-{
-public:
- LinearSystem(Matrix & _matrix, Vector & _vector)
- : m_matrix(_matrix), m_vector(_vector), m_solution(_matrix.columns())
- {
- }
-
- const Vector & LU_solve()
- {
- assert( matrix().rows() == matrix().columns()
- && matrix().rows() == vector().size() );
- int s;
- gsl_permutation * p = gsl_permutation_alloc(matrix().rows());
- gsl_linalg_LU_decomp (matrix().get_gsl_matrix(), p, &s);
- gsl_linalg_LU_solve( matrix().get_gsl_matrix(),
- p,
- vector().get_gsl_vector(),
- m_solution.get_gsl_vector()
- );
- gsl_permutation_free(p);
- return solution();
- }
-
- const Vector & SV_solve()
- {
- assert( matrix().rows() >= matrix().columns()
- && matrix().rows() == vector().size() );
-
- gsl_matrix* U = matrix().get_gsl_matrix();
- gsl_matrix* V = gsl_matrix_alloc(matrix().columns(), matrix().columns());
- gsl_vector* S = gsl_vector_alloc(matrix().columns());
- gsl_vector* work = gsl_vector_alloc(matrix().columns());
-
- gsl_linalg_SV_decomp( U, V, S, work );
-
- gsl_vector* b = vector().get_gsl_vector();
- gsl_vector* x = m_solution.get_gsl_vector();
-
- gsl_linalg_SV_solve( U, V, S, b, x);
-
- gsl_matrix_free(V);
- gsl_vector_free(S);
- gsl_vector_free(work);
-
- return solution();
- }
-
- Matrix & matrix()
- {
- return m_matrix;
- }
-
- Vector & vector()
- {
- return m_vector;
- }
-
- const Vector & solution() const
- {
- return m_solution;
- }
-
-private:
- Matrix & m_matrix;
- Vector & m_vector;
- Vector m_solution;
-};
-
-
-} } // end namespaces
-
-
-#endif /*_NL_LINEAR_SYSTEM_H_*/
+/*\r
+ * LinearSystem class wraps some gsl routines for solving linear systems\r
+ *\r
+ * Authors:\r
+ * Marco Cecchetti <mrcekets at gmail.com>\r
+ * \r
+ * Copyright 2008 authors\r
+ *\r
+ * This library is free software; you can redistribute it and/or\r
+ * modify it either under the terms of the GNU Lesser General Public\r
+ * License version 2.1 as published by the Free Software Foundation\r
+ * (the "LGPL") or, at your option, under the terms of the Mozilla\r
+ * Public License Version 1.1 (the "MPL"). If you do not alter this\r
+ * notice, a recipient may use your version of this file under either\r
+ * the MPL or the LGPL.\r
+ *\r
+ * You should have received a copy of the LGPL along with this library\r
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software\r
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
+ * You should have received a copy of the MPL along with this library\r
+ * in the file COPYING-MPL-1.1\r
+ *\r
+ * The contents of this file are subject to the Mozilla Public License\r
+ * Version 1.1 (the "License"); you may not use this file except in\r
+ * compliance with the License. You may obtain a copy of the License at\r
+ * http://www.mozilla.org/MPL/\r
+ *\r
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY\r
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for\r
+ * the specific language governing rights and limitations.\r
+ */\r
+\r
+\r
+#ifndef _NL_LINEAR_SYSTEM_H_\r
+#define _NL_LINEAR_SYSTEM_H_\r
+\r
+\r
+#include <cassert>\r
+\r
+#include <gsl/gsl_linalg.h>\r
+\r
+#include <2geom/numeric/matrix.h>\r
+#include <2geom/numeric/vector.h>\r
+\r
+\r
+namespace Geom { namespace NL {\r
+\r
+\r
+class LinearSystem\r
+{\r
+public:\r
+ LinearSystem(MatrixView & _matrix, VectorView & _vector)\r
+ : m_matrix(_matrix), m_vector(_vector), m_solution(_matrix.columns())\r
+ {\r
+ }\r
+ \r
+ LinearSystem(Matrix & _matrix, Vector & _vector)\r
+ : m_matrix(_matrix), m_vector(_vector), m_solution(_matrix.columns())\r
+ {\r
+ }\r
+ \r
+ const Vector & LU_solve()\r
+ {\r
+ assert( matrix().rows() == matrix().columns() \r
+ && matrix().rows() == vector().size() );\r
+ int s;\r
+ gsl_permutation * p = gsl_permutation_alloc(matrix().rows());\r
+ gsl_linalg_LU_decomp (matrix().get_gsl_matrix(), p, &s);\r
+ gsl_linalg_LU_solve( matrix().get_gsl_matrix(), \r
+ p, \r
+ vector().get_gsl_vector(), \r
+ m_solution.get_gsl_vector()\r
+ );\r
+ gsl_permutation_free(p);\r
+ return solution();\r
+ }\r
+ \r
+ const Vector & SV_solve()\r
+ {\r
+ assert( matrix().rows() >= matrix().columns()\r
+ && matrix().rows() == vector().size() );\r
+ \r
+ gsl_matrix* U = matrix().get_gsl_matrix();\r
+ gsl_matrix* V = gsl_matrix_alloc(matrix().columns(), matrix().columns());\r
+ gsl_vector* S = gsl_vector_alloc(matrix().columns());\r
+ gsl_vector* work = gsl_vector_alloc(matrix().columns());\r
+ \r
+ gsl_linalg_SV_decomp( U, V, S, work );\r
+ \r
+ gsl_vector* b = vector().get_gsl_vector();\r
+ gsl_vector* x = m_solution.get_gsl_vector();\r
+ \r
+ gsl_linalg_SV_solve( U, V, S, b, x);\r
+ \r
+ gsl_matrix_free(V);\r
+ gsl_vector_free(S);\r
+ gsl_vector_free(work);\r
+ \r
+ return solution(); \r
+ }\r
+ \r
+ MatrixView & matrix()\r
+ {\r
+ return m_matrix;\r
+ }\r
+ \r
+ VectorView & vector()\r
+ {\r
+ return m_vector;\r
+ }\r
+ \r
+ const Vector & solution() const\r
+ {\r
+ return m_solution;\r
+ }\r
+ \r
+private:\r
+ MatrixView m_matrix;\r
+ VectorView m_vector;\r
+ Vector m_solution;\r
+};\r
+\r
+\r
+} } // end namespaces\r
+\r
+\r
+#endif /*_NL_LINEAR_SYSTEM_H_*/\r
+\r
+/*\r
+ Local Variables:\r
+ mode:c++\r
+ c-file-style:"stroustrup"\r
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))\r
+ indent-tabs-mode:nil\r
+ fill-column:99\r
+ End:\r
+*/\r
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :\r
diff --git a/src/2geom/numeric/matrix.cpp b/src/2geom/numeric/matrix.cpp
--- /dev/null
@@ -0,0 +1,115 @@
+/*
+ * Matrix, MatrixView, ConstMatrixView classes wrap the gsl matrix routines;
+ * "views" mimic the semantic of C++ references: any operation performed
+ * on a "view" is actually performed on the "viewed object"
+ *
+ * Authors:
+ * Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2008 authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+
+#include <numeric/matrix.h>
+#include <numeric/vector.h>
+
+
+
+namespace Geom { namespace NL {
+
+
+Vector operator*( detail::BaseMatrixImpl const& A,
+ detail::BaseVectorImpl const& v )
+{
+ assert(A.columns() == v.size());
+
+ Vector result(A.rows(), 0.0);
+ for (size_t i = 0; i < A.rows(); ++i)
+ for (size_t j = 0; j < A.columns(); ++j)
+ result[i] += A(i,j) * v[j];
+
+ return result;
+}
+
+Matrix operator*( detail::BaseMatrixImpl const& A,
+ detail::BaseMatrixImpl const& B )
+{
+ assert(A.columns() == B.rows());
+
+ Matrix C(A.rows(), B.columns(), 0.0);
+ for (size_t i = 0; i < C.rows(); ++i)
+ for (size_t j = 0; j < C.columns(); ++j)
+ for (size_t k = 0; k < A.columns(); ++k)
+ C(i,j) += A(i,k) * B(k, j);
+
+ return C;
+}
+
+Matrix pseudo_inverse(detail::BaseMatrixImpl const& A)
+{
+
+ Matrix U(A);
+ Matrix V(A.columns(), A.columns());
+ Vector s(A.columns());
+ gsl_vector* work = gsl_vector_alloc(A.columns());
+
+ gsl_linalg_SV_decomp( U.get_gsl_matrix(),
+ V.get_gsl_matrix(),
+ s.get_gsl_vector(),
+ work );
+
+ Matrix P(A.columns(), A.rows(), 0.0);
+
+ int sz = s.size();
+ while ( sz-- > 0 && s[sz] == 0 ) {}
+ ++sz;
+ if (sz == 0) return P;
+ VectorView sv(s, sz);
+
+ for (size_t i = 0; i < sv.size(); ++i)
+ {
+ VectorView v = V.column_view(i);
+ v.scale(1/sv[i]);
+ for (size_t h = 0; h < P.rows(); ++h)
+ for (size_t k = 0; k < P.columns(); ++k)
+ P(h,k) += V(h,i) * U(k,i);
+ }
+
+ return P;
+}
+
+} } // end namespaces
+
+/*
+ Local Variables:
+ mode:c++
+ c-file-style:"stroustrup"
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+ indent-tabs-mode:nil
+ fill-column:99
+ End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
index bdd647e53167fad3d34d143122d9497ebfdcd45c..64557a6f12ed64ffc28260871136adc3641501b9 100644 (file)
-#ifndef _NL_MATRIX_H_
-#define _NL_MATRIX_H_
-
-
-#include <cassert>
-#include <utility>
-#include <algorithm>
-
-#include <gsl/gsl_matrix.h>
-
-
-
-namespace Geom { namespace NL {
-
-class Matrix;
-void swap(Matrix & m1, Matrix & m2);
-
-
-class Matrix
-{
-public:
- // the matrix is not inizialized
- Matrix(size_t n1, size_t n2)
- : m_rows(n1), m_columns(n2)
- {
- m_matrix = gsl_matrix_alloc(n1, n2);
- }
-
- Matrix(size_t n1, size_t n2, double x)
- :m_rows(n1), m_columns(n2)
- {
- m_matrix = gsl_matrix_alloc(n1, n2);
- gsl_matrix_set_all(m_matrix, x);
- }
-
- Matrix( Matrix const& _matrix )
- : m_rows(_matrix.rows()), m_columns(_matrix.columns())
- {
- m_matrix = gsl_matrix_alloc(rows(), columns());
- gsl_matrix_memcpy(m_matrix, _matrix.m_matrix);
- }
-
-// explicit
-// Matrix( gsl_matrix* m, size_t n1, size_t n2)
-// : m_rows(n1), m_columns(n2), m_matrix(m)
-// {
-// }
-
- Matrix & operator=(Matrix const& _matrix)
- {
- assert( rows() == _matrix.rows() && columns() == _matrix.columns() );
- gsl_matrix_memcpy(m_matrix, _matrix.m_matrix);
-
- return *this;
- }
-
- virtual ~Matrix()
- {
- gsl_matrix_free(m_matrix);
- }
-
- void set_all( double x = 0)
- {
- gsl_matrix_set_all(m_matrix, x);
- }
-
- void set_identity()
- {
- gsl_matrix_set_identity(m_matrix);
- }
-
- const double & operator() (size_t i, size_t j) const
- {
- return *gsl_matrix_const_ptr(m_matrix, i, j);
- }
-
- double & operator() (size_t i, size_t j)
- {
- return *gsl_matrix_ptr(m_matrix, i, j);
- }
-
- gsl_matrix* get_gsl_matrix()
- {
- return m_matrix;
- }
-
- void swap_rows(size_t i, size_t j)
- {
- gsl_matrix_swap_rows(m_matrix, i, j);
- }
-
- void swap_columns(size_t i, size_t j)
- {
- gsl_matrix_swap_columns(m_matrix, i, j);
- }
-
- Matrix & transpose()
- {
- gsl_matrix_transpose(m_matrix);
- return (*this);
- }
-
- Matrix & scale(double x)
- {
- gsl_matrix_scale(m_matrix, x);
- return (*this);
- }
-
- Matrix & translate(double x)
- {
- gsl_matrix_add_constant(m_matrix, x);
- return (*this);
- }
-
- Matrix & operator+=(Matrix const& _matrix)
- {
- gsl_matrix_add(m_matrix, _matrix.m_matrix);
- return (*this);
- }
-
- Matrix & operator-=(Matrix const& _matrix)
- {
- gsl_matrix_sub(m_matrix, _matrix.m_matrix);
- return (*this);
- }
-
- bool is_zero() const
- {
- return gsl_matrix_isnull(m_matrix);
- }
-
- bool is_positive() const
- {
- return gsl_matrix_ispos(m_matrix);
- }
-
- bool is_negative() const
- {
- return gsl_matrix_isneg(m_matrix);
- }
-
- bool is_non_negative() const
- {
- for ( unsigned int i = 0; i < rows(); ++i )
- {
- for ( unsigned int j = 0; j < columns(); ++j )
- {
- if ( (*this)(i,j) < 0 ) return false;
- }
- }
- return true;
- }
-
- double max() const
- {
- return gsl_matrix_max(m_matrix);
- }
-
- double min() const
- {
- return gsl_matrix_min(m_matrix);
- }
-
- std::pair<size_t, size_t>
- max_index() const
- {
- std::pair<size_t, size_t> indices;
- gsl_matrix_max_index(m_matrix, &(indices.first), &(indices.second));
- return indices;
- }
-
- std::pair<size_t, size_t>
- min_index() const
- {
- std::pair<size_t, size_t> indices;
- gsl_matrix_min_index(m_matrix, &(indices.first), &(indices.second));
- return indices;
- }
-
- friend
- void swap(Matrix & m1, Matrix & m2);
-
- size_t rows() const
- {
- return m_rows;
- }
-
- size_t columns() const
- {
- return m_columns;
- }
-
-private:
- size_t m_rows, m_columns;
- gsl_matrix* m_matrix;
-};
-
-void swap(Matrix & m1, Matrix & m2)
-{
- assert( m1.rows() == m2.rows() && m1.columns() == m2.columns() );
- std::swap(m1.m_matrix, m2.m_matrix);
-}
-
-
-} } // end namespaces
-
-#endif /*_NL_MATRIX_H_*/
+/*\r
+ * Matrix, MatrixView, ConstMatrixView classes wrap the gsl matrix routines;\r
+ * "views" mimic the semantic of C++ references: any operation performed\r
+ * on a "view" is actually performed on the "viewed object"\r
+ *\r
+ * Authors:\r
+ * Marco Cecchetti <mrcekets at gmail.com>\r
+ *\r
+ * Copyright 2008 authors\r
+ *\r
+ * This library is free software; you can redistribute it and/or\r
+ * modify it either under the terms of the GNU Lesser General Public\r
+ * License version 2.1 as published by the Free Software Foundation\r
+ * (the "LGPL") or, at your option, under the terms of the Mozilla\r
+ * Public License Version 1.1 (the "MPL"). If you do not alter this\r
+ * notice, a recipient may use your version of this file under either\r
+ * the MPL or the LGPL.\r
+ *\r
+ * You should have received a copy of the LGPL along with this library\r
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software\r
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
+ * You should have received a copy of the MPL along with this library\r
+ * in the file COPYING-MPL-1.1\r
+ *\r
+ * The contents of this file are subject to the Mozilla Public License\r
+ * Version 1.1 (the "License"); you may not use this file except in\r
+ * compliance with the License. You may obtain a copy of the License at\r
+ * http://www.mozilla.org/MPL/\r
+ *\r
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY\r
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for\r
+ * the specific language governing rights and limitations.\r
+ */\r
+\r
+\r
+\r
+\r
+#ifndef _NL_MATRIX_H_\r
+#define _NL_MATRIX_H_\r
+\r
+#include <2geom/numeric/vector.h>\r
+\r
+#include <cassert>\r
+#include <utility> // for std::pair\r
+#include <algorithm> // for std::swap\r
+#include <sstream>\r
+#include <string>\r
+\r
+#include <gsl/gsl_matrix.h>\r
+#include <gsl/gsl_linalg.h>\r
+\r
+\r
+namespace Geom { namespace NL {\r
+\r
+namespace detail\r
+{\r
+\r
+class BaseMatrixImpl\r
+{\r
+ public:\r
+ virtual ~BaseMatrixImpl()\r
+ {\r
+ }\r
+\r
+ ConstVectorView row_const_view(size_t i) const\r
+ {\r
+ return ConstVectorView(gsl_matrix_const_row(m_matrix, i));\r
+ }\r
+\r
+ ConstVectorView column_const_view(size_t i) const\r
+ {\r
+ return ConstVectorView(gsl_matrix_const_column(m_matrix, i));\r
+ }\r
+\r
+ const double & operator() (size_t i, size_t j) const\r
+ {\r
+ return *gsl_matrix_const_ptr(m_matrix, i, j);\r
+ }\r
+\r
+ const gsl_matrix* get_gsl_matrix() const\r
+ {\r
+ return m_matrix;\r
+ }\r
+\r
+ bool is_zero() const\r
+ {\r
+ return gsl_matrix_isnull(m_matrix);\r
+ }\r
+\r
+ bool is_positive() const\r
+ {\r
+ return gsl_matrix_ispos(m_matrix);\r
+ }\r
+\r
+ bool is_negative() const\r
+ {\r
+ return gsl_matrix_isneg(m_matrix);\r
+ }\r
+\r
+ bool is_non_negative() const\r
+ {\r
+ for ( unsigned int i = 0; i < rows(); ++i )\r
+ {\r
+ for ( unsigned int j = 0; j < columns(); ++j )\r
+ {\r
+ if ( (*this)(i,j) < 0 ) return false;\r
+ }\r
+ }\r
+ return true;\r
+ }\r
+\r
+ double max() const\r
+ {\r
+ return gsl_matrix_max(m_matrix);\r
+ }\r
+\r
+ double min() const\r
+ {\r
+ return gsl_matrix_min(m_matrix);\r
+ }\r
+\r
+ std::pair<size_t, size_t>\r
+ max_index() const\r
+ {\r
+ std::pair<size_t, size_t> indices;\r
+ gsl_matrix_max_index(m_matrix, &(indices.first), &(indices.second));\r
+ return indices;\r
+ }\r
+\r
+ std::pair<size_t, size_t>\r
+ min_index() const\r
+ {\r
+ std::pair<size_t, size_t> indices;\r
+ gsl_matrix_min_index(m_matrix, &(indices.first), &(indices.second));\r
+ return indices;\r
+ }\r
+\r
+ size_t rows() const\r
+ {\r
+ return m_rows;\r
+ }\r
+\r
+ size_t columns() const\r
+ {\r
+ return m_columns;\r
+ }\r
+\r
+ std::string str() const;\r
+\r
+ protected:\r
+ size_t m_rows, m_columns;\r
+ gsl_matrix* m_matrix;\r
+\r
+}; // end class BaseMatrixImpl\r
+\r
+\r
+inline\r
+bool operator== (BaseMatrixImpl const& m1, BaseMatrixImpl const& m2)\r
+{\r
+ if (m1.rows() != m2.rows() || m1.columns() != m2.columns()) return false;\r
+\r
+ for (size_t i = 0; i < m1.rows(); ++i)\r
+ for (size_t j = 0; j < m1.columns(); ++j)\r
+ if (m1(i,j) != m2(i,j)) return false;\r
+\r
+ return true;\r
+}\r
+\r
+template< class charT >\r
+inline\r
+std::basic_ostream<charT> &\r
+operator<< (std::basic_ostream<charT> & os, const BaseMatrixImpl & _matrix)\r
+{\r
+ if (_matrix.rows() == 0 || _matrix.columns() == 0) return os;\r
+\r
+ os << "[[" << _matrix(0,0);\r
+ for (size_t j = 1; j < _matrix.columns(); ++j)\r
+ {\r
+ os << ", " << _matrix(0,j);\r
+ }\r
+ os << "]";\r
+\r
+ for (size_t i = 1; i < _matrix.rows(); ++i)\r
+ {\r
+ os << ", [" << _matrix(i,0);\r
+ for (size_t j = 1; j < _matrix.columns(); ++j)\r
+ {\r
+ os << ", " << _matrix(i,j);\r
+ }\r
+ os << "]";\r
+ }\r
+ os << "]";\r
+ return os;\r
+}\r
+\r
+inline\r
+std::string BaseMatrixImpl::str() const\r
+{\r
+ std::ostringstream oss;\r
+ oss << (*this);\r
+ return oss.str();\r
+}\r
+\r
+\r
+class MatrixImpl : public BaseMatrixImpl\r
+{\r
+ public:\r
+\r
+ typedef BaseMatrixImpl base_type;\r
+\r
+ void set_all( double x )\r
+ {\r
+ gsl_matrix_set_all(m_matrix, x);\r
+ }\r
+\r
+ void set_identity()\r
+ {\r
+ gsl_matrix_set_identity(m_matrix);\r
+ }\r
+\r
+ using base_type::operator();\r
+\r
+ double & operator() (size_t i, size_t j)\r
+ {\r
+ return *gsl_matrix_ptr(m_matrix, i, j);\r
+ }\r
+\r
+ using base_type::get_gsl_matrix;\r
+\r
+ gsl_matrix* get_gsl_matrix()\r
+ {\r
+ return m_matrix;\r
+ }\r
+\r
+ VectorView row_view(size_t i)\r
+ {\r
+ return VectorView(gsl_matrix_row(m_matrix, i));\r
+ }\r
+\r
+ VectorView column_view(size_t i)\r
+ {\r
+ return VectorView(gsl_matrix_column(m_matrix, i));\r
+ }\r
+\r
+ void swap_rows(size_t i, size_t j)\r
+ {\r
+ gsl_matrix_swap_rows(m_matrix, i, j);\r
+ }\r
+\r
+ void swap_columns(size_t i, size_t j)\r
+ {\r
+ gsl_matrix_swap_columns(m_matrix, i, j);\r
+ }\r
+\r
+ MatrixImpl & transpose()\r
+ {\r
+ assert(columns() == rows());\r
+ gsl_matrix_transpose(m_matrix);\r
+ return (*this);\r
+ }\r
+\r
+ MatrixImpl & scale(double x)\r
+ {\r
+ gsl_matrix_scale(m_matrix, x);\r
+ return (*this);\r
+ }\r
+\r
+ MatrixImpl & translate(double x)\r
+ {\r
+ gsl_matrix_add_constant(m_matrix, x);\r
+ return (*this);\r
+ }\r
+\r
+ MatrixImpl & operator+=(base_type const& _matrix)\r
+ {\r
+ gsl_matrix_add(m_matrix, _matrix.get_gsl_matrix());\r
+ return (*this);\r
+ }\r
+\r
+ MatrixImpl & operator-=(base_type const& _matrix)\r
+ {\r
+ gsl_matrix_sub(m_matrix, _matrix.get_gsl_matrix());\r
+ return (*this);\r
+ }\r
+\r
+}; // end class MatrixImpl\r
+\r
+} // end namespace detail\r
+\r
+\r
+using detail::operator==;\r
+using detail::operator<<;\r
+\r
+\r
+\r
+\r
+class Matrix: public detail::MatrixImpl\r
+{\r
+ public:\r
+ typedef detail::MatrixImpl base_type;\r
+\r
+ public:\r
+ // the matrix is not inizialized\r
+ Matrix(size_t n1, size_t n2)\r
+ {\r
+ m_rows = n1;\r
+ m_columns = n2;\r
+ m_matrix = gsl_matrix_alloc(n1, n2);\r
+ }\r
+\r
+ Matrix(size_t n1, size_t n2, double x)\r
+ {\r
+ m_rows = n1;\r
+ m_columns = n2;\r
+ m_matrix = gsl_matrix_alloc(n1, n2);\r
+ gsl_matrix_set_all(m_matrix, x);\r
+ }\r
+\r
+ Matrix(Matrix const& _matrix)\r
+ : base_type()\r
+ {\r
+ m_rows = _matrix.rows();\r
+ m_columns = _matrix.columns();\r
+ m_matrix = gsl_matrix_alloc(rows(), columns());\r
+ gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());\r
+ }\r
+\r
+ explicit\r
+ Matrix(base_type::base_type const& _matrix)\r
+ {\r
+ m_rows = _matrix.rows();\r
+ m_columns = _matrix.columns();\r
+ m_matrix = gsl_matrix_alloc(rows(), columns());\r
+ gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());\r
+ }\r
+\r
+ Matrix & operator=(Matrix const& _matrix)\r
+ {\r
+ assert( rows() == _matrix.rows() && columns() == _matrix.columns() );\r
+ gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());\r
+ return *this;\r
+ }\r
+\r
+ Matrix & operator=(base_type::base_type const& _matrix)\r
+ {\r
+ assert( rows() == _matrix.rows() && columns() == _matrix.columns() );\r
+ gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());\r
+ return *this;\r
+ }\r
+\r
+ virtual ~Matrix()\r
+ {\r
+ gsl_matrix_free(m_matrix);\r
+ }\r
+\r
+ Matrix & transpose()\r
+ {\r
+ return static_cast<Matrix &>( base_type::transpose() );\r
+ }\r
+\r
+ Matrix & scale(double x)\r
+ {\r
+ return static_cast<Matrix &>( base_type::scale(x) );\r
+ }\r
+\r
+ Matrix & translate(double x)\r
+ {\r
+ return static_cast<Matrix &>( base_type::translate(x) );\r
+ }\r
+\r
+ Matrix & operator+=(base_type::base_type const& _matrix)\r
+ {\r
+ return static_cast<Matrix &>( base_type::operator+=(_matrix) );\r
+ }\r
+\r
+ Matrix & operator-=(base_type::base_type const& _matrix)\r
+ {\r
+ return static_cast<Matrix &>( base_type::operator-=(_matrix) );\r
+ }\r
+\r
+ friend\r
+ void swap(Matrix & m1, Matrix & m2);\r
+ friend\r
+ void swap_any(Matrix & m1, Matrix & m2);\r
+\r
+}; // end class Matrix\r
+\r
+\r
+// warning! this operation invalidates any view of the passed matrix objects\r
+inline\r
+void swap(Matrix & m1, Matrix & m2)\r
+{\r
+ assert( m1.rows() == m2.rows() && m1.columns() == m2.columns() );\r
+ std::swap(m1.m_matrix, m2.m_matrix);\r
+}\r
+\r
+inline\r
+void swap_any(Matrix & m1, Matrix & m2)\r
+{\r
+ std::swap(m1.m_matrix, m2.m_matrix);\r
+ std::swap(m1.m_rows, m2.m_rows);\r
+ std::swap(m1.m_columns, m2.m_columns);\r
+}\r
+\r
+\r
+\r
+class ConstMatrixView : public detail::BaseMatrixImpl\r
+{\r
+ public:\r
+ typedef detail::BaseMatrixImpl base_type;\r
+\r
+ public:\r
+ ConstMatrixView(const base_type & _matrix, size_t k1, size_t k2, size_t n1, size_t n2)\r
+ : m_matrix_view( gsl_matrix_const_submatrix(_matrix.get_gsl_matrix(), k1, k2, n1, n2) )\r
+ {\r
+ m_rows = n1;\r
+ m_columns = n2;\r
+ m_matrix = const_cast<gsl_matrix*>( &(m_matrix_view.matrix) );\r
+ }\r
+\r
+ ConstMatrixView(const ConstMatrixView & _matrix)\r
+ : base_type(),\r
+ m_matrix_view(_matrix.m_matrix_view)\r
+ {\r
+ m_rows = _matrix.rows();\r
+ m_columns = _matrix.columns();\r
+ m_matrix = const_cast<gsl_matrix*>( &(m_matrix_view.matrix) );\r
+ }\r
+\r
+ ConstMatrixView(const base_type & _matrix)\r
+ : m_matrix_view(gsl_matrix_const_submatrix(_matrix.get_gsl_matrix(), 0, 0, _matrix.rows(), _matrix.columns()))\r
+ {\r
+ m_rows = _matrix.rows();\r
+ m_columns = _matrix.columns();\r
+ m_matrix = const_cast<gsl_matrix*>( &(m_matrix_view.matrix) );\r
+ }\r
+\r
+ private:\r
+ gsl_matrix_const_view m_matrix_view;\r
+\r
+}; // end class ConstMatrixView\r
+\r
+\r
+\r
+\r
+class MatrixView : public detail::MatrixImpl\r
+{\r
+ public:\r
+ typedef detail::MatrixImpl base_type;\r
+\r
+ public:\r
+ MatrixView(base_type & _matrix, size_t k1, size_t k2, size_t n1, size_t n2)\r
+ {\r
+ m_rows = n1;\r
+ m_columns = n2;\r
+ m_matrix_view\r
+ = gsl_matrix_submatrix(_matrix.get_gsl_matrix(), k1, k2, n1, n2);\r
+ m_matrix = &(m_matrix_view.matrix);\r
+ }\r
+\r
+ MatrixView(const MatrixView & _matrix)\r
+ : base_type()\r
+ {\r
+ m_rows = _matrix.rows();\r
+ m_columns = _matrix.columns();\r
+ m_matrix_view = _matrix.m_matrix_view;\r
+ m_matrix = &(m_matrix_view.matrix);\r
+ }\r
+\r
+ MatrixView(Matrix & _matrix)\r
+ {\r
+ m_rows = _matrix.rows();\r
+ m_columns = _matrix.columns();\r
+ m_matrix_view\r
+ = gsl_matrix_submatrix(_matrix.get_gsl_matrix(), 0, 0, rows(), columns());\r
+ m_matrix = &(m_matrix_view.matrix);\r
+ }\r
+\r
+ MatrixView & operator=(MatrixView const& _matrix)\r
+ {\r
+ assert( rows() == _matrix.rows() && columns() == _matrix.columns() );\r
+ gsl_matrix_memcpy(m_matrix, _matrix.m_matrix);\r
+ return *this;\r
+ }\r
+\r
+ MatrixView & operator=(base_type::base_type const& _matrix)\r
+ {\r
+ assert( rows() == _matrix.rows() && columns() == _matrix.columns() );\r
+ gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());\r
+ return *this;\r
+ }\r
+\r
+ MatrixView & transpose()\r
+ {\r
+ return static_cast<MatrixView &>( base_type::transpose() );\r
+ }\r
+\r
+ MatrixView & scale(double x)\r
+ {\r
+ return static_cast<MatrixView &>( base_type::scale(x) );\r
+ }\r
+\r
+ MatrixView & translate(double x)\r
+ {\r
+ return static_cast<MatrixView &>( base_type::translate(x) );\r
+ }\r
+\r
+ MatrixView & operator+=(base_type::base_type const& _matrix)\r
+ {\r
+ return static_cast<MatrixView &>( base_type::operator+=(_matrix) );\r
+ }\r
+\r
+ MatrixView & operator-=(base_type::base_type const& _matrix)\r
+ {\r
+ return static_cast<MatrixView &>( base_type::operator-=(_matrix) );\r
+ }\r
+\r
+ friend\r
+ void swap_view(MatrixView & m1, MatrixView & m2);\r
+\r
+ private:\r
+ gsl_matrix_view m_matrix_view;\r
+\r
+}; // end class MatrixView\r
+\r
+\r
+inline\r
+void swap_view(MatrixView & m1, MatrixView & m2)\r
+{\r
+ assert( m1.rows() == m2.rows() && m1.columns() == m2.columns() );\r
+ std::swap(m1.m_matrix_view, m2.m_matrix_view);\r
+}\r
+\r
+Vector operator*( detail::BaseMatrixImpl const& A,\r
+ detail::BaseVectorImpl const& v );\r
+\r
+Matrix operator*( detail::BaseMatrixImpl const& A,\r
+ detail::BaseMatrixImpl const& B );\r
+\r
+Matrix pseudo_inverse(detail::BaseMatrixImpl const& A);\r
+\r
+} } // end namespaces\r
+\r
+#endif /*_NL_MATRIX_H_*/\r
+\r
+/*\r
+ Local Variables:\r
+ mode:c++\r
+ c-file-style:"stroustrup"\r
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))\r
+ indent-tabs-mode:nil\r
+ fill-column:99\r
+ End:\r
+*/\r
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :\r
index b25861e22c7d353fac06412d462dd8e1869856ec..43a39a1ac5684ba5eba76c751e8d7336e456ef07 100644 (file)
-#ifndef _NL_VECTOR_H_
-#define _NL_VECTOR_H_
-
-#include <cassert>
-#include <utility>
-
-#include <gsl/gsl_vector.h>
-
-
-namespace Geom { namespace NL {
-
-class Vector;
-void swap(Vector & v1, Vector & v2);
-
-
-class Vector
-{
-public:
- Vector(size_t n)
- : m_size(n)
- {
- m_vector = gsl_vector_alloc(n);
- }
-
- Vector(size_t n, double x)
- : m_size(n)
- {
- m_vector = gsl_vector_alloc(n);
- gsl_vector_set_all(m_vector, x);
- }
-
- // create a vector with n elements all set to zero
- // but the i-th that is set to 1
- Vector(size_t n, size_t i)
- : m_size(n)
- {
- m_vector = gsl_vector_alloc(n);
- gsl_vector_set_basis(m_vector, i);
- }
-
- Vector(Vector const& _vector)
- : m_size(_vector.size())
- {
- m_vector = gsl_vector_alloc(size());
- gsl_vector_memcpy(m_vector, _vector.m_vector);
- }
-
- virtual ~Vector()
- {
- gsl_vector_free(m_vector);
- }
-
- Vector & operator=(Vector const& _vector)
- {
- assert( size() == _vector.size() );
- gsl_vector_memcpy(m_vector, _vector.m_vector);
- return (*this);
- }
-
- void set_all(double x = 0)
- {
- gsl_vector_set_all(m_vector, x);
- }
-
- void set_basis(size_t i)
- {
- gsl_vector_set_basis(m_vector, i);
- }
-
- double const& operator[](size_t i) const
- {
- return *gsl_vector_const_ptr(m_vector, i);
- }
-
- double & operator[](size_t i)
- {
- return *gsl_vector_ptr(m_vector, i);
- }
-
- gsl_vector* get_gsl_vector()
- {
- return m_vector;
- }
-
- void swap_elements(size_t i, size_t j)
- {
- gsl_vector_swap_elements(m_vector, i, j);
- }
-
- void reverse()
- {
- gsl_vector_reverse(m_vector);
- }
-
- Vector & scale(double x)
- {
- gsl_vector_scale(m_vector, x);
- return (*this);
- }
-
- Vector & translate(double x)
- {
- gsl_vector_add_constant(m_vector, x);
- return (*this);
- }
-
- Vector & operator+=(Vector const& _vector)
- {
- gsl_vector_add(m_vector, _vector.m_vector);
- return (*this);
- }
-
- Vector & operator-=(Vector const& _vector)
- {
- gsl_vector_sub(m_vector, _vector.m_vector);
- return (*this);
- }
-
- bool is_zero() const
- {
- return gsl_vector_isnull(m_vector);
- }
-
- bool is_positive() const
- {
- return gsl_vector_ispos(m_vector);
- }
-
- bool is_negative() const
- {
- return gsl_vector_isneg(m_vector);
- }
-
- bool is_non_negative() const
- {
- for ( size_t i = 0; i < size(); ++i )
- {
- if ( (*this)[i] < 0 ) return false;
- }
- return true;
- }
-
- double max() const
- {
- return gsl_vector_max(m_vector);
- }
-
- double min() const
- {
- return gsl_vector_min(m_vector);
- }
-
- size_t max_index() const
- {
- return gsl_vector_max_index(m_vector);
- }
-
- size_t min_index() const
- {
- return gsl_vector_min_index(m_vector);
- }
-
- friend
- void swap(Vector & v1, Vector & v2);
-
- size_t size() const
- {
- return m_size;
- }
-
-private:
- size_t m_size;
- gsl_vector* m_vector;
-};
-
-void swap(Vector & v1, Vector & v2)
-{
- assert( v1.size() == v2.size() );
- std::swap(v1.m_vector, v2.m_vector);
-}
-
-} } // end namespaces
-
-
-#endif /*_NL_VECTOR_H_*/
+/*\r
+ * Vector, VectorView, ConstVectorView classes wrap the gsl vector routines;\r
+ * "views" mimic the semantic of C++ references: any operation performed\r
+ * on a "view" is actually performed on the "viewed object"\r
+ *\r
+ * Authors:\r
+ * Marco Cecchetti <mrcekets at gmail.com>\r
+ *\r
+ * Copyright 2008 authors\r
+ *\r
+ * This library is free software; you can redistribute it and/or\r
+ * modify it either under the terms of the GNU Lesser General Public\r
+ * License version 2.1 as published by the Free Software Foundation\r
+ * (the "LGPL") or, at your option, under the terms of the Mozilla\r
+ * Public License Version 1.1 (the "MPL"). If you do not alter this\r
+ * notice, a recipient may use your version of this file under either\r
+ * the MPL or the LGPL.\r
+ *\r
+ * You should have received a copy of the LGPL along with this library\r
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software\r
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
+ * You should have received a copy of the MPL along with this library\r
+ * in the file COPYING-MPL-1.1\r
+ *\r
+ * The contents of this file are subject to the Mozilla Public License\r
+ * Version 1.1 (the "License"); you may not use this file except in\r
+ * compliance with the License. You may obtain a copy of the License at\r
+ * http://www.mozilla.org/MPL/\r
+ *\r
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY\r
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for\r
+ * the specific language governing rights and limitations.\r
+ */\r
+\r
+\r
+\r
+\r
+#ifndef _NL_VECTOR_H_\r
+#define _NL_VECTOR_H_\r
+\r
+#include <cassert>\r
+#include <algorithm> // for std::swap\r
+#include <vector>\r
+#include <sstream>\r
+#include <string>\r
+\r
+\r
+#include <gsl/gsl_vector.h>\r
+#include <gsl/gsl_blas.h>\r
+\r
+\r
+namespace Geom { namespace NL {\r
+\r
+namespace detail\r
+{\r
+\r
+class BaseVectorImpl\r
+{\r
+ public:\r
+ double const& operator[](size_t i) const\r
+ {\r
+ return *gsl_vector_const_ptr(m_vector, i);\r
+ }\r
+\r
+ const gsl_vector* get_gsl_vector() const\r
+ {\r
+ return m_vector;\r
+ }\r
+ bool is_zero() const\r
+ {\r
+ return gsl_vector_isnull(m_vector);\r
+ }\r
+\r
+ bool is_positive() const\r
+ {\r
+ return gsl_vector_ispos(m_vector);\r
+ }\r
+\r
+ bool is_negative() const\r
+ {\r
+ return gsl_vector_isneg(m_vector);\r
+ }\r
+\r
+ bool is_non_negative() const\r
+ {\r
+ for ( size_t i = 0; i < size(); ++i )\r
+ {\r
+ if ( (*this)[i] < 0 ) return false;\r
+ }\r
+ return true;\r
+ }\r
+\r
+ double max() const\r
+ {\r
+ return gsl_vector_max(m_vector);\r
+ }\r
+\r
+ double min() const\r
+ {\r
+ return gsl_vector_min(m_vector);\r
+ }\r
+\r
+ size_t max_index() const\r
+ {\r
+ return gsl_vector_max_index(m_vector);\r
+ }\r
+\r
+ size_t min_index() const\r
+ {\r
+ return gsl_vector_min_index(m_vector);\r
+ }\r
+\r
+ size_t size() const\r
+ {\r
+ return m_size;\r
+ }\r
+\r
+ std::string str() const;\r
+\r
+ virtual ~BaseVectorImpl()\r
+ {\r
+ }\r
+\r
+ protected:\r
+ size_t m_size;\r
+ gsl_vector* m_vector;\r
+\r
+}; // end class BaseVectorImpl\r
+\r
+\r
+inline\r
+bool operator== (BaseVectorImpl const& v1, BaseVectorImpl const& v2)\r
+{\r
+ if (v1.size() != v2.size()) return false;\r
+\r
+ for (size_t i = 0; i < v1.size(); ++i)\r
+ {\r
+ if (v1[i] != v2[i]) return false;\r
+ }\r
+ return true;\r
+}\r
+\r
+template< class charT >\r
+inline\r
+std::basic_ostream<charT> &\r
+operator<< (std::basic_ostream<charT> & os, const BaseVectorImpl & _vector)\r
+{\r
+ if (_vector.size() == 0 ) return os;\r
+ os << "[" << _vector[0];\r
+ for (unsigned int i = 1; i < _vector.size(); ++i)\r
+ {\r
+ os << ", " << _vector[i];\r
+ }\r
+ os << "]";\r
+ return os;\r
+}\r
+\r
+inline\r
+std::string BaseVectorImpl::str() const\r
+{\r
+ std::ostringstream oss;\r
+ oss << (*this);\r
+ return oss.str();\r
+}\r
+\r
+inline\r
+double dot(BaseVectorImpl const& v1, BaseVectorImpl const& v2)\r
+{\r
+ double result;\r
+ gsl_blas_ddot(v1.get_gsl_vector(), v2.get_gsl_vector(), &result);\r
+ return result;\r
+}\r
+\r
+\r
+class VectorImpl : public BaseVectorImpl\r
+{\r
+ public:\r
+ typedef BaseVectorImpl base_type;\r
+\r
+ public:\r
+ void set_all(double x)\r
+ {\r
+ gsl_vector_set_all(m_vector, x);\r
+ }\r
+\r
+ void set_basis(size_t i)\r
+ {\r
+ gsl_vector_set_basis(m_vector, i);\r
+ }\r
+\r
+ using base_type::operator[];\r
+\r
+ double & operator[](size_t i)\r
+ {\r
+ return *gsl_vector_ptr(m_vector, i);\r
+ }\r
+\r
+ using base_type::get_gsl_vector;\r
+\r
+ gsl_vector* get_gsl_vector()\r
+ {\r
+ return m_vector;\r
+ }\r
+\r
+ void swap_elements(size_t i, size_t j)\r
+ {\r
+ gsl_vector_swap_elements(m_vector, i, j);\r
+ }\r
+\r
+ void reverse()\r
+ {\r
+ gsl_vector_reverse(m_vector);\r
+ }\r
+\r
+ VectorImpl & scale(double x)\r
+ {\r
+ gsl_vector_scale(m_vector, x);\r
+ return (*this);\r
+ }\r
+\r
+ VectorImpl & translate(double x)\r
+ {\r
+ gsl_vector_add_constant(m_vector, x);\r
+ return (*this);\r
+ }\r
+\r
+ VectorImpl & operator+=(base_type const& _vector)\r
+ {\r
+ gsl_vector_add(m_vector, _vector.get_gsl_vector());\r
+ return (*this);\r
+ }\r
+\r
+ VectorImpl & operator-=(base_type const& _vector)\r
+ {\r
+ gsl_vector_sub(m_vector, _vector.get_gsl_vector());\r
+ return (*this);\r
+ }\r
+\r
+}; // end class VectorImpl\r
+\r
+} // end namespace detail\r
+\r
+\r
+using detail::operator==;\r
+using detail::operator<<;\r
+\r
+class Vector : public detail::VectorImpl\r
+{\r
+ public:\r
+ typedef detail::VectorImpl base_type;\r
+\r
+ public:\r
+ Vector(size_t n)\r
+ {\r
+ m_size = n;\r
+ m_vector = gsl_vector_alloc(n);\r
+ }\r
+\r
+ Vector(size_t n, double x)\r
+ {\r
+ m_size = n;\r
+ m_vector = gsl_vector_alloc(n);\r
+ gsl_vector_set_all(m_vector, x);\r
+ }\r
+\r
+ // create a vector with n elements all set to zero\r
+ // but the i-th that is set to 1\r
+ Vector(size_t n, size_t i)\r
+ {\r
+ m_size = n;\r
+ m_vector = gsl_vector_alloc(n);\r
+ gsl_vector_set_basis(m_vector, i);\r
+ }\r
+\r
+ Vector(Vector const& _vector)\r
+ : base_type()\r
+ {\r
+ m_size = _vector.size();\r
+ m_vector = gsl_vector_alloc(size());\r
+ gsl_vector_memcpy(m_vector, _vector.m_vector);\r
+ }\r
+\r
+ explicit\r
+ Vector(base_type::base_type const& _vector)\r
+ {\r
+ m_size = _vector.size();\r
+ m_vector = gsl_vector_alloc(size());\r
+ gsl_vector_memcpy(m_vector, _vector.get_gsl_vector());\r
+ }\r
+\r
+ virtual ~Vector()\r
+ {\r
+ gsl_vector_free(m_vector);\r
+ }\r
+\r
+\r
+ Vector & operator=(Vector const& _vector)\r
+ {\r
+ assert( size() == _vector.size() );\r
+ gsl_vector_memcpy(m_vector, _vector.m_vector);\r
+ return (*this);\r
+ }\r
+\r
+ Vector & operator=(base_type::base_type const& _vector)\r
+ {\r
+ assert( size() == _vector.size() );\r
+ gsl_vector_memcpy(m_vector, _vector.get_gsl_vector());\r
+ return (*this);\r
+ }\r
+\r
+ Vector & scale(double x)\r
+ {\r
+ return static_cast<Vector&>( base_type::scale(x) );\r
+ }\r
+\r
+ Vector & translate(double x)\r
+ {\r
+ return static_cast<Vector&>( base_type::translate(x) );\r
+ }\r
+\r
+ Vector & operator+=(base_type::base_type const& _vector)\r
+ {\r
+ return static_cast<Vector&>( base_type::operator+=(_vector) );\r
+ }\r
+\r
+ Vector & operator-=(base_type::base_type const& _vector)\r
+ {\r
+ return static_cast<Vector&>( base_type::operator-=(_vector) );\r
+ }\r
+\r
+ friend\r
+ void swap(Vector & v1, Vector & v2);\r
+ friend\r
+ void swap_any(Vector & v1, Vector & v2);\r
+\r
+}; // end class Vector\r
+\r
+\r
+// warning! these operations invalidate any view of the passed vector objects\r
+inline\r
+void swap(Vector & v1, Vector & v2)\r
+{\r
+ assert( v1.size() == v2.size() );\r
+ std::swap(v1.m_vector, v2.m_vector);\r
+}\r
+\r
+inline\r
+void swap_any(Vector & v1, Vector & v2)\r
+{\r
+ std::swap(v1.m_vector, v2.m_vector);\r
+ std::swap(v1.m_size, v2.m_size);\r
+}\r
+\r
+\r
+class ConstVectorView : public detail::BaseVectorImpl\r
+{\r
+ public:\r
+ typedef detail::BaseVectorImpl base_type;\r
+\r
+ public:\r
+ ConstVectorView(const base_type & _vector, size_t n, size_t offset = 0)\r
+ : m_vector_view( gsl_vector_const_subvector(_vector.get_gsl_vector(), offset, n) )\r
+ {\r
+ m_size = n;\r
+ m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );\r
+ }\r
+\r
+ ConstVectorView(const base_type & _vector, size_t n, size_t offset , size_t stride)\r
+ : m_vector_view( gsl_vector_const_subvector_with_stride(_vector.get_gsl_vector(), offset, stride, n) )\r
+ {\r
+ m_size = n;\r
+ m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );\r
+ }\r
+\r
+ ConstVectorView(const double* _vector, size_t n, size_t offset = 0)\r
+ : m_vector_view( gsl_vector_const_view_array(_vector + offset, n) )\r
+ {\r
+ m_size = n;\r
+ m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );\r
+ }\r
+\r
+ ConstVectorView(const double* _vector, size_t n, size_t offset, size_t stride)\r
+ : m_vector_view( gsl_vector_const_view_array_with_stride(_vector + offset, stride, n) )\r
+ {\r
+ m_size = n;\r
+ m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );\r
+ }\r
+\r
+ explicit\r
+ ConstVectorView(gsl_vector_const_view _gsl_vector_view)\r
+ : m_vector_view(_gsl_vector_view)\r
+ {\r
+ m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );\r
+ m_size = m_vector->size;\r
+ }\r
+\r
+ explicit\r
+ ConstVectorView(const std::vector<double>& _vector)\r
+ : m_vector_view( gsl_vector_const_view_array(&(_vector[0]), _vector.size()) )\r
+ {\r
+ m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );\r
+ m_size = _vector.size();\r
+ }\r
+\r
+ ConstVectorView(const ConstVectorView & _vector)\r
+ : base_type(),\r
+ m_vector_view(_vector.m_vector_view)\r
+ {\r
+ m_size = _vector.size();\r
+ m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );\r
+ }\r
+\r
+ ConstVectorView(const base_type & _vector)\r
+ : m_vector_view(gsl_vector_const_subvector(_vector.get_gsl_vector(), 0, _vector.size()))\r
+ {\r
+ m_size = _vector.size();\r
+ m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );\r
+ }\r
+\r
+ private:\r
+ gsl_vector_const_view m_vector_view;\r
+\r
+}; // end class ConstVectorView\r
+\r
+\r
+\r
+\r
+class VectorView : public detail::VectorImpl\r
+{\r
+ public:\r
+ typedef detail::VectorImpl base_type;\r
+\r
+ public:\r
+ VectorView(base_type & _vector, size_t n, size_t offset = 0, size_t stride = 1)\r
+ {\r
+ m_size = n;\r
+ if (stride == 1)\r
+ {\r
+ m_vector_view\r
+ = gsl_vector_subvector(_vector.get_gsl_vector(), offset, n);\r
+ m_vector = &(m_vector_view.vector);\r
+ }\r
+ else\r
+ {\r
+ m_vector_view\r
+ = gsl_vector_subvector_with_stride(_vector.get_gsl_vector(), offset, stride, n);\r
+ m_vector = &(m_vector_view.vector);\r
+ }\r
+ }\r
+\r
+ VectorView(double* _vector, size_t n, size_t offset = 0, size_t stride = 1)\r
+ {\r
+ m_size = n;\r
+ if (stride == 1)\r
+ {\r
+ m_vector_view\r
+ = gsl_vector_view_array(_vector + offset, n);\r
+ m_vector = &(m_vector_view.vector);\r
+ }\r
+ else\r
+ {\r
+ m_vector_view\r
+ = gsl_vector_view_array_with_stride(_vector + offset, stride, n);\r
+ m_vector = &(m_vector_view.vector);\r
+ }\r
+\r
+ }\r
+\r
+ VectorView(const VectorView & _vector)\r
+ : base_type()\r
+ {\r
+ m_size = _vector.size();\r
+ m_vector_view = _vector.m_vector_view;\r
+ m_vector = &(m_vector_view.vector);\r
+ }\r
+\r
+ VectorView(Vector & _vector)\r
+ {\r
+ m_size = _vector.size();\r
+ m_vector_view = gsl_vector_subvector(_vector.get_gsl_vector(), 0, size());\r
+ m_vector = &(m_vector_view.vector);\r
+ }\r
+\r
+ explicit\r
+ VectorView(gsl_vector_view _gsl_vector_view)\r
+ : m_vector_view(_gsl_vector_view)\r
+ {\r
+ m_vector = &(m_vector_view.vector);\r
+ m_size = m_vector->size;\r
+ }\r
+\r
+ explicit\r
+ VectorView(std::vector<double> & _vector)\r
+ {\r
+ m_size = _vector.size();\r
+ m_vector_view = gsl_vector_view_array(&(_vector[0]), _vector.size());\r
+ m_vector = &(m_vector_view.vector);\r
+ }\r
+\r
+ VectorView & operator=(VectorView const& _vector)\r
+ {\r
+ assert( size() == _vector.size() );\r
+ gsl_vector_memcpy(m_vector, _vector.get_gsl_vector());\r
+ return (*this);\r
+ }\r
+\r
+ VectorView & operator=(base_type::base_type const& _vector)\r
+ {\r
+ assert( size() == _vector.size() );\r
+ gsl_vector_memcpy(m_vector, _vector.get_gsl_vector());\r
+ return (*this);\r
+ }\r
+\r
+ VectorView & scale(double x)\r
+ {\r
+ return static_cast<VectorView&>( base_type::scale(x) );\r
+ }\r
+\r
+ VectorView & translate(double x)\r
+ {\r
+ return static_cast<VectorView&>( base_type::translate(x) );\r
+ }\r
+\r
+ VectorView & operator+=(base_type::base_type const& _vector)\r
+ {\r
+ return static_cast<VectorView&>( base_type::operator+=(_vector) );\r
+ }\r
+\r
+ VectorView & operator-=(base_type::base_type const& _vector)\r
+ {\r
+ return static_cast<VectorView&>( base_type::operator-=(_vector) );\r
+ }\r
+\r
+ friend\r
+ void swap_view(VectorView & v1, VectorView & v2);\r
+\r
+ private:\r
+ gsl_vector_view m_vector_view;\r
+\r
+}; // end class VectorView\r
+\r
+\r
+inline\r
+void swap_view(VectorView & v1, VectorView & v2)\r
+{\r
+ assert( v1.size() == v2.size() );\r
+ std::swap(v1.m_vector_view, v2.m_vector_view); // not swap m_vector too\r
+}\r
+\r
+\r
+} } // end namespaces\r
+\r
+\r
+#endif /*_NL_VECTOR_H_*/\r
+\r
+/*\r
+ Local Variables:\r
+ mode:c++\r
+ c-file-style:"stroustrup"\r
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))\r
+ indent-tabs-mode:nil\r
+ fill-column:99\r
+ End:\r
+*/\r
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :\r
index 2c32c376415ffc4fe7820cf9c069905229a87180..635996e515fde2ff8c12e732959db78d27c777dc 100644 (file)
-#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 {
*/
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 << "]";
index 5cdee5509a2701b9b2e204d09eb4c5f8eddbd120..2094c4beb501c13879a89562b8fb58b873491b88 100644 (file)
#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 8bfccf87a33dd2a82c428ff3b74e01925bbbc926..2e5b789d5936252895e72c07ac37c308232589b5 100644 (file)
--- a/src/2geom/path.cpp
+++ b/src/2geom/path.cpp
-#include "path.h"
+#include <2geom/path.h>
namespace Geom
}
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 ) {
}
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 ) {
return bounds;
}
+/*
+Rect Path::boundsFast()
+{
+ Rect bound;
+ if (empty()) return bound;
+
+ bound = begin()->boundsFast();
+ for (iterator it = ++begin(); it != end(); ++it)
+ {
+ bound.unionWith(it->boundsFast());
+ }
+ return bound;
+}
+
+Rect Path::boundsExact()
+{
+ Rect bound;
+ if (empty()) return bound;
+
+ bound = begin()->boundsExact();
+ for (iterator it = ++begin(); it != end(); ++it)
+ {
+ bound.unionWith(it->boundsExact());
+ }
+ return bound;
+ }*/
+
template<typename iter>
iter inc(iter const &x, unsigned n) {
iter ret = x;
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;
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;
}
if(ff != 1.) {
Curve *fromv = fromi->portion(ff, 1.);
//fromv->setInitial(ret.finalPoint());
- ret.append(*fromv);
+ ret.append(*fromv, STITCH_DISCONTINUOUS);
delete fromv;
}
if(from >= to) {
const_iterator ender = end();
if(ender->initialPoint() == ender->finalPoint()) ender++;
- ret.insert(ret.end(), ++fromi, ender);
- ret.insert(ret.end(), begin(), toi);
+ ret.insert(ret.end(), ++fromi, ender, STITCH_DISCONTINUOUS);
+ ret.insert(ret.end(), begin(), toi, STITCH_DISCONTINUOUS);
} else {
- ret.insert(ret.end(), ++fromi, toi);
+ ret.insert(ret.end(), ++fromi, toi, STITCH_DISCONTINUOUS);
}
Curve *tov = toi->portion(0., tf);
- ret.append(*tov);
+ ret.append(*tov, STITCH_DISCONTINUOUS);
delete tov;
}
-const double eps = .1;
-
-void Path::append(Curve const &curve) {
- if ( curves_.front() != final_ && !are_near(curve.initialPoint(), (*final_)[0], eps) ) {
- THROW_CONTINUITYERROR();
- }
- do_append(curve.duplicate());
-}
-
-void Path::append(D2<SBasis> const &curve) {
- if ( curves_.front() != final_ ) {
- for ( int i = 0 ; i < 2 ; ++i ) {
- if ( !are_near(curve[i][0][0], (*final_)[0][i], eps) ) {
- THROW_CONTINUITYERROR();
- }
- }
- }
- do_append(new SBasisCurve(curve));
-}
-
-void Path::append(Path const &other)
-{
- // Check that path stays continuous:
- if ( !are_near( finalPoint(), other.initialPoint() ) ) {
- THROW_CONTINUITYERROR();
- }
-
- insert(begin(), other.begin(), other.end());
-}
-
void Path::do_update(Sequence::iterator first_replaced,
Sequence::iterator last_replaced,
Sequence::iterator first,
- Sequence::iterator last)
+ Sequence::iterator last)
{
// note: modifies the contents of [first,last)
-
check_continuity(first_replaced, last_replaced, first, last);
delete_range(first_replaced, last_replaced);
if ( ( last - first ) == ( last_replaced - first_replaced ) ) {
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());
}
}
+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,
{
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 5f2c86b957dba17a155da2f9725d428dc1f7e766..1c142f98d542e42d74d2a9fc079df629136af09b 100644 (file)
--- a/src/2geom/path.h
+++ b/src/2geom/path.h
#define SEEN_GEOM_PATH_H
-#include "curves.h"
+#include <2geom/curves.h>
#include <iterator>
namespace Geom
{
+// Conditional expression for types. If true, first, if false, second.
+template<bool _Cond, typename _Iftrue, typename _Iffalse>
+ struct __conditional_type
+ { typedef _Iftrue __type; };
+
+template<typename _Iftrue, typename _Iffalse>
+ struct __conditional_type<false, _Iftrue, _Iffalse>
+ { typedef _Iffalse __type; };
+
+
template <typename IteratorImpl>
class BaseIterator
: public std::iterator<std::forward_iterator_tag, Curve const>
// default construct
// default copy
+ // Allow Sequence::iterator to Sequence::const_iterator conversion
+ // unfortunately I do not know how to imitate the way __normal_iterator
+ // does it, because I don't see a way to get the typename of the container
+ // IteratorImpl is pointing at...
+ typedef std::vector<Curve *> Sequence;
+ BaseIterator ( typename __conditional_type<
+ (std::__are_same<IteratorImpl, Sequence::const_iterator >::__value), // check if this instantiation is of const_iterator type
+ const BaseIterator< Sequence::iterator >, // if true: accept iterator in const_iterator instantiation
+ const BaseIterator<IteratorImpl> > ::__type // if false: default to standard copy constructor
+ & __other)
+ : impl_(__other.impl_) { }
+ friend class BaseIterator< Sequence::const_iterator >;
+
+
bool operator==(BaseIterator const &other) {
return other.impl_ == impl_;
}
typedef Sequence::size_type size_type;
typedef Sequence::difference_type difference_type;
+ class ClosingSegment : public LineSegment {
+ public:
+ ClosingSegment() : LineSegment() {}
+ ClosingSegment(Point const &p1, Point const &p2) : LineSegment(p1, p2) {}
+ virtual Curve *duplicate() const { return new ClosingSegment(*this); }
+ };
+
+ enum Stitching {
+ NO_STITCHING=0,
+ STITCH_DISCONTINUOUS
+ };
+
+ class StitchSegment : public LineSegment {
+ public:
+ StitchSegment() : LineSegment() {}
+ StitchSegment(Point const &p1, Point const &p2) : LineSegment(p1, p2) {}
+ virtual Curve *duplicate() const { return new StitchSegment(*this); }
+ };
+
Path()
- : final_(new LineSegment()), closed_(false)
+ : final_(new ClosingSegment()), closed_(false)
{
curves_.push_back(final_);
}
Path(Path const &other)
- : final_(new LineSegment()), closed_(other.closed_)
+ : final_(new ClosingSegment()), closed_(other.closed_)
{
curves_.push_back(final_);
insert(begin(), other.begin(), other.end());
}
explicit Path(Point p)
- : final_(new LineSegment(p, p)), closed_(false)
+ : final_(new ClosingSegment(p, p)), closed_(false)
{
curves_.push_back(final_);
}
template <typename Impl>
Path(BaseIterator<Impl> first, BaseIterator<Impl> last, bool closed=false)
- : closed_(closed), final_(new LineSegment())
+ : closed_(closed), final_(new ClosingSegment())
{
curves_.push_back(final_);
insert(begin(), first, last);
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();
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 {
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());
}
template <typename Impl>
- void insert(iterator pos, BaseIterator<Impl> first, BaseIterator<Impl> last)
+ void insert(iterator pos, BaseIterator<Impl> first, BaseIterator<Impl> last, Stitching stitching=NO_STITCHING)
{
Sequence source(DuplicatingIterator<Impl>(first.impl_),
DuplicatingIterator<Impl>(last.impl_));
try {
+ if (stitching) stitch(pos.impl_, pos.impl_, source);
do_update(pos.impl_, pos.impl_, source.begin(), source.end());
} catch (...) {
delete_range(source.begin(), source.end());
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
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());
}
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 (...) {
template <typename Impl>
void replace(iterator replaced,
- BaseIterator<Impl> first, BaseIterator<Impl> last)
+ BaseIterator<Impl> first, BaseIterator<Impl> last,
+ Stitching stitching=NO_STITCHING)
{
Sequence source(DuplicatingIterator<Impl>(first.impl_),
DuplicatingIterator<Impl>(last.impl_));
try {
+ if (stitching) stitch(replaced.impl_, replaced.impl_+1, source);
do_update(replaced.impl_, replaced.impl_+1, source.begin(), source.end());
} catch (...) {
delete_range(source.begin(), source.end());
template <typename Impl>
void replace(iterator first_replaced, iterator last_replaced,
- BaseIterator<Impl> first, BaseIterator<Impl> last)
+ BaseIterator<Impl> first, BaseIterator<Impl> last,
+ Stitching stitching=NO_STITCHING)
{
Sequence source(first.impl_, last.impl_);
try {
+ if (stitching) stitch(first_replaced.impl_, last_replaced.impl_, source);
do_update(first_replaced.impl_, last_replaced.impl_,
source.begin(), source.end());
} catch (...) {
}
}
- void append(Curve const &curve);
- void append(D2<SBasis> const &curve);
- void append(Path const &other);
+ void append(Curve const &curve, Stitching stitching=NO_STITCHING) {
+ if (stitching) stitchTo(curve.initialPoint());
+ do_append(curve.duplicate());
+ }
+ void append(D2<SBasis> const &curve, Stitching stitching=NO_STITCHING) {
+ if (stitching) stitchTo(Point(curve[X][0][0], curve[Y][0][0]));
+ do_append(new SBasisCurve(curve));
+ }
+ void append(Path const &other, Stitching stitching=NO_STITCHING) {
+ insert(end(), other.begin(), other.end(), stitching);
+ }
+
+ void stitchTo(Point const &p) {
+ if (!empty() && finalPoint() != p) {
+ do_append(new StitchSegment(finalPoint(), p));
+ }
+ }
template <typename CurveType, typename A>
void appendNew(A a) {
- do_append(new CurveType((*final_)[0], a));
+ do_append(new CurveType(finalPoint(), a));
}
template <typename CurveType, typename A, typename B>
void appendNew(A a, B b) {
- do_append(new CurveType((*final_)[0], a, b));
+ do_append(new CurveType(finalPoint(), a, b));
}
template <typename CurveType, typename A, typename B, typename C>
void appendNew(A a, B b, C c) {
- do_append(new CurveType((*final_)[0], a, b, c));
+ do_append(new CurveType(finalPoint(), a, b, c));
}
template <typename CurveType, typename A, typename B, typename C,
typename D>
void appendNew(A a, B b, C c, D d) {
- do_append(new CurveType((*final_)[0], a, b, c, d));
+ do_append(new CurveType(finalPoint(), a, b, c, d));
}
template <typename CurveType, typename A, typename B, typename C,
typename D, typename E>
void appendNew(A a, B b, C c, D d, E e) {
- do_append(new CurveType((*final_)[0], a, b, c, d, e));
+ do_append(new CurveType(finalPoint(), a, b, c, d, e));
}
template <typename CurveType, typename A, typename B, typename C,
typename D, typename E, typename F>
void appendNew(A a, B b, C c, D d, E e, F f) {
- do_append(new CurveType((*final_)[0], a, b, c, d, e, f));
+ do_append(new CurveType(finalPoint(), a, b, c, d, e, f));
}
template <typename CurveType, typename A, typename B, typename C,
typename D, typename E, typename F,
typename G>
void appendNew(A a, B b, C c, D d, E e, F f, G g) {
- do_append(new CurveType((*final_)[0], a, b, c, d, e, f, g));
+ do_append(new CurveType(finalPoint(), a, b, c, d, e, f, g));
}
template <typename CurveType, typename A, typename B, typename C,
typename D, typename E, typename F,
typename G, typename H>
void appendNew(A a, B b, C c, D d, E e, F f, G g, H h) {
- do_append(new CurveType((*final_)[0], a, b, c, d, e, f, g, h));
+ do_append(new CurveType(finalPoint(), a, b, c, d, e, f, g, h));
}
template <typename CurveType, typename A, typename B, typename C,
typename D, typename E, typename F,
typename G, typename H, typename I>
void appendNew(A a, B b, C c, D d, E e, F f, G g, H h, I i) {
- do_append(new CurveType((*final_)[0], a, b, c, d, e, f, g, h, i));
+ do_append(new CurveType(finalPoint(), a, b, c, d, e, f, g, h, i));
}
private:
+ void stitch(Sequence::iterator first_replaced,
+ Sequence::iterator last_replaced,
+ Sequence &sequence);
+
void do_update(Sequence::iterator first_replaced,
Sequence::iterator last_replaced,
Sequence::iterator first,
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,
Sequence::iterator last);
Sequence curves_;
- LineSegment *final_;
+ ClosingSegment *final_;
bool closed_;
}; // end class Path
return c.nearestPoint(p);
}
-
-/*
-class PathPortion : public Curve {
- Path *source;
- double f, t;
- boost::optional<Path> result;
-
- public:
- double from() const { return f; }
- double to() const { return t; }
-
- explicit PathPortion(Path *s, double fp, double tp) : source(s), f(fp), t(tp) {}
- Curve *duplicate() const { return new PathPortion(*this); }
-
- Point initialPoint() const { return source->pointAt(f); }
- Point finalPoint() const { return source->pointAt(t); }
-
- Path actualPath() {
- if(!result) *result = source->portion(f, t);
- return *result;
- }
-
- Rect boundsFast() const { return actualPath().boundsFast; }
- Rect boundsExact() const { return actualPath().boundsFast; }
- Rect boundsLocal(Interval i) const { THROW_NOTIMPLEMENTED(); }
-
- std::vector<double> roots(double v, Dim2 d) const = 0;
-
- virtual int winding(Point p) const { return root_winding(*this, p); }
-
- virtual Curve *portion(double f, double t) const = 0;
- virtual Curve *reverse() const { return portion(1, 0); }
-
- virtual Crossings crossingsWith(Curve const & other) const;
-
- virtual void setInitial(Point v) = 0;
- virtual void setFinal(Point v) = 0;
-
- virtual Curve *transformed(Matrix const &m) const = 0;
-
- virtual Point pointAt(Coord t) const { return pointAndDerivatives(t, 0).front(); }
- virtual Coord valueAt(Coord t, Dim2 d) const { return pointAt(t)[d]; }
- virtual std::vector<Point> pointAndDerivatives(Coord t, unsigned n) const = 0;
- virtual D2<SBasis> toSBasis() const = 0;
-
-};
-*/
-
} // end namespace Geom
namespace std {
index fe2f1976cdf73d755e7d9ff5b4339317c6b888e2..9befc1f170688edfee205f5ddd12bf70b459e3ee 100644 (file)
--- a/src/2geom/pathvector.cpp
+++ b/src/2geom/pathvector.cpp
#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 6a759e40638d9957b297e5f7b2240c3f38b7e757..17f6b62127e2afe976a225e89b843982766ece63 100644 (file)
--- a/src/2geom/pathvector.h
+++ b/src/2geom/pathvector.h
#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 {
index 222152aa36fd8209c39f7f9fb2bbc56ca22a9ca1..2d8638818c18868e50997ba962d925af5db99d29 100644 (file)
--- a/src/2geom/piecewise.cpp
+++ b/src/2geom/piecewise.cpp
*
*/
-#include "piecewise.h"
+#include <2geom/piecewise.h>
#include <iterator>
#include <map>
return result;
}
-Piecewise<SBasis> reverse(Piecewise<SBasis> const &f) {
- Piecewise<SBasis> ret = Piecewise<SBasis>();
- ret.cuts.resize(f.cuts.size());
- ret.segs.resize(f.segs.size());
- double start = f.cuts[0];
- double end = f.cuts.back();
- for (unsigned i = 0; i < f.cuts.size(); i++) {
- double x = f.cuts[f.cuts.size() - 1 - i];
- ret.cuts[i] = end - (x - start);
- }
- for (unsigned i = 0; i < f.segs.size(); i++)
- ret.segs[i] = reverse(f[f.segs.size() - i - 1]);
- return ret;
-}
-
}
/*
diff --git a/src/2geom/piecewise.h b/src/2geom/piecewise.h
index 96e64bcf9e40ba3ab53e9dedbf7e2d32c6ce85d4..73fc76e2ffcdd72f46a89ebe1592b2d25813cc01 100644 (file)
--- a/src/2geom/piecewise.h
+++ b/src/2geom/piecewise.h
#ifndef SEEN_GEOM_PW_SB_H
#define SEEN_GEOM_PW_SB_H
-#include "sbasis.h"
+#include <2geom/sbasis.h>
#include <vector>
#include <map>
-#include "concepts.h"
-#include "isnan.h"
+#include <2geom/concepts.h>
+#include <2geom/isnan.h>
#include <boost/concept_check.hpp>
namespace Geom {
diff --git a/src/2geom/point-l.h b/src/2geom/point-l.h
index cd4768db08d7792383d95688d3c4428f728c8911..88c4af9482a8fcdb6b93b75e143d754e8f027364 100644 (file)
--- a/src/2geom/point-l.h
+++ b/src/2geom/point-l.h
#define SEEN_Geom_POINT_L_H
#include <stdexcept>
-#include "point.h"
+#include <2geom/point.h>
namespace Geom {
diff --git a/src/2geom/point.cpp b/src/2geom/point.cpp
index 66eaf8ca671f1a06e49ce518be323167beecea46..42a7ecef305bf100fb168fb803a3726c82488882 100644 (file)
--- a/src/2geom/point.cpp
+++ b/src/2geom/point.cpp
-#include "point.h"
+#include <2geom/point.h>
#include <assert.h>
-#include "coord.h"
-#include "isnan.h" //temporary fix for isnan()
-#include "matrix.h"
+#include <2geom/coord.h>
+#include <2geom/isnan.h> //temporary fix for isnan()
+#include <2geom/matrix.h>
namespace Geom {
diff --git a/src/2geom/point.h b/src/2geom/point.h
index 91df159cad90f5364b3cdda88e8cbdaac5c5ee95..2cab3d7fe056356f095177cb6495bcda286b2d4f 100644 (file)
--- a/src/2geom/point.h
+++ b/src/2geom/point.h
#include <iostream>
-#include "coord.h"
-#include "utils.h"
+#include <2geom/coord.h>
+#include <2geom/utils.h>
namespace Geom {
_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
index fdc1cefe55a585d2205a63fcf925b1a74a587058..2787c4478bb7f3918100f2d2ff821d181a12dd09 100644 (file)
-#include "poly-dk-solve.h"
+#include <2geom/poly-dk-solve.h>
#include <iterator>
/*** implementation of the Durand-Kerner method. seems buggy*/
index f82caf394961e26ee47f5d52575a66e5df829da8..dadd1211094ffd32bb7594c51464aac138d341bd 100644 (file)
-#include "poly.h"
+#include <2geom/poly.h>
#include <complex>
std::vector<std::complex<double> >
index 266e24c2b93c690babeb0f12891fb3f4fbd76f81..bd2ddf305235e358b96aa221a2611ae48420d3a7 100644 (file)
-#include "poly-laguerre-solve.h"
+#include <2geom/poly-laguerre-solve.h>
#include <iterator>
typedef std::complex<double> cdouble;
//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)
}
std::vector<double>
-laguerre_real_interval(Poly const & ply,
- const double lo, const double hi,
- const double tol)
+laguerre_real_interval(Poly const & /*ply*/,
+ const double /*lo*/, const double /*hi*/,
+ const double /*tol*/)
{
/* not implemented*/
assert(false);
index c7bfef24559f93ce59c528a5420aa8a17257b8c7..6f7cb0aee86b0b3ff9dd0156e7998477a1baba51 100644 (file)
-#include "poly.h"
+#include <2geom/poly.h>
#include <complex>
std::vector<std::complex<double> >
diff --git a/src/2geom/poly.cpp b/src/2geom/poly.cpp
index 27f98596b531eb7820ab8091e284cd851ea84837..5f766395779b8661db364916eddcddf0885f5152 100644 (file)
--- a/src/2geom/poly.cpp
+++ b/src/2geom/poly.cpp
-#include "poly.h"
+#include <2geom/poly.h>
Poly Poly::operator*(const Poly& p) const {
Poly result;
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 f653e80a644604e229e2eaa311e15115975dae23..2af24d25fe14c1b96653545e207db44942d5cd54 100644 (file)
--- a/src/2geom/poly.h
+++ b/src/2geom/poly.h
#include <iostream>
#include <algorithm>
#include <complex>
-#include "utils.h"
+#include <2geom/utils.h>
class Poly : public std::vector<double>{
public:
diff --git a/src/2geom/quadtree.cpp b/src/2geom/quadtree.cpp
index bb041edbe1fd7d77359da7c5857a06bf7b8cf4ab..e322a091b8b10b4a27036deed066b0612f6ba226 100644 (file)
--- a/src/2geom/quadtree.cpp
+++ b/src/2geom/quadtree.cpp
-#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 716c56673661cf080bb290e4ff7b5d372ef2b8ca..0770cd5f864bdfa1fe45e747de7357cb5fcbce33 100644 (file)
--- a/src/2geom/quadtree.h
+++ b/src/2geom/quadtree.h
#include <vector>
#include <cassert>
-#include "d2.h"
+#include <2geom/d2.h>
namespace Geom{
diff --git a/src/2geom/rect.h b/src/2geom/rect.h
index c89946606fe86886924c9f54ae778882430dec1e..27a475cf3476d93da06fd4ec52c3824289d4578c 100644 (file)
--- a/src/2geom/rect.h
+++ b/src/2geom/rect.h
#ifndef _2GEOM_RECT
#define _2GEOM_RECT
-#include "matrix.h"
+#include <2geom/matrix.h>
#include <boost/optional/optional.hpp>
namespace Geom {
diff --git a/src/2geom/region.cpp b/src/2geom/region.cpp
index 1bc4a68e79ca77419d0ce7006561fc0085146650..065f3f418b3eb166359708ac6c43c7d7658a28fc 100644 (file)
--- a/src/2geom/region.cpp
+++ b/src/2geom/region.cpp
-#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 561f3184dc1c84d3b577348839760c48f37f7c4b..960a570a4c6622aaa56c1dd21a0730a5b30fd95c 100644 (file)
--- a/src/2geom/region.h
+++ b/src/2geom/region.h
#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 {
index a4decd704503e067f052ec6986c57003f48643fb..79e99519af330e12a308764120f42d575ec02f53 100644 (file)
--- a/src/2geom/sbasis-2d.cpp
+++ b/src/2geom/sbasis-2d.cpp
-#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 921a740b9b9e5c2777109c191244168521eb8a12..1845b5c99305b6bc5923651f13639f41d3f56c29 100644 (file)
--- a/src/2geom/sbasis-2d.h
+++ b/src/2geom/sbasis-2d.h
#include <vector>
#include <cassert>
#include <algorithm>
-#include "d2.h"
-#include "sbasis.h"
+#include <2geom/d2.h>
+#include <2geom/sbasis.h>
#include <iostream>
namespace Geom{
index 1612f46f72cfec8edc1a385d503d15ed45d88fb9..6b3888594a19c691e93461ec466d446535ac4250 100644 (file)
--- a/src/2geom/sbasis-curve.h
+++ b/src/2geom/sbasis-curve.h
#define _2GEOM_SBASIS_CURVE_H_
-#include "curve.h"
+#include <2geom/curve.h>
namespace Geom
index 6558303891278c6ef0ea4e71e420c01cf4d3dfab..887ca9995f5184349ac5b09dd3c9c93ed6de4928 100644 (file)
-#include "sbasis-geometric.h"
-#include "sbasis.h"
-#include "sbasis-math.h"
-//#include "solver.h"
-#include "sbasis-geometric.h"
+#include <2geom/sbasis-geometric.h>
+#include <2geom/sbasis.h>
+#include <2geom/sbasis-math.h>
+//#include <2geom/solver.h>
+#include <2geom/sbasis-geometric.h>
/** Geometric operators on D2<SBasis> (1D->2D).
* Copyright 2007 JF Barraud
index 6e21a6395c51aab5b52212608821feaa72d798fd..7e067d8013c3abe91468956daa9c27f452c2afeb 100644 (file)
#ifndef _SBASIS_GEOMETRIC
#define _SBASIS_GEOMETRIC
-#include "d2.h"
-#include "piecewise.h"
+#include <2geom/d2.h>
+#include <2geom/piecewise.h>
#include <vector>
/** two-dimensional geometric operators.
index b04e15bcca96a81e671c3354c32784e44ce62b64..f5a8ab7a10207768724e4e93cf86b868a5621495 100644 (file)
//TODO: define a truncated compose(sb,sb, order) and extend it to pw<sb>.
//TODO: in all these functions, compute 'order' according to 'tol'.
-#include "sbasis-math.h"
+#include <2geom/sbasis-math.h>
//#define ZERO 1e-3
index 72428bc963d5ae24e5ee6e45d60bca7fcb5316ab..2a51bffb3fdd048ad297ce2667d4c802f0e87843 100644 (file)
--- a/src/2geom/sbasis-math.h
+++ b/src/2geom/sbasis-math.h
#define SEEN_GEOM_SB_CALCULS_H
-#include "sbasis.h"
-#include "piecewise.h"
+#include <2geom/sbasis.h>
+#include <2geom/piecewise.h>
namespace Geom{
//-|x|---------------------------------------------------------------
index ff797920f7d2b376600a41169e7a7c7b2ee68ce6..3fee0afe4e6b1796c98fe393026cbdd0cf507e62 100644 (file)
-#include "sbasis-poly.h"
+#include <2geom/sbasis-poly.h>
namespace Geom{
index 4a8aa1cd8ea8740647385d50815bc67a20eecbf2..09c8e54870f7928c7f576de6902923293996c00c 100644 (file)
--- a/src/2geom/sbasis-poly.h
+++ b/src/2geom/sbasis-poly.h
#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.
index 51a4ec97d998949813c19ee3766318f627c48d69..30f5fb22fdbfd7e853ea1fa39faeb143207a957a 100644 (file)
#include <cmath>
#include <map>
-#include "sbasis.h"
-#include "sbasis-to-bezier.h"
-#include "solver.h"
+#include <2geom/sbasis.h>
+#include <2geom/sbasis-to-bezier.h>
+#include <2geom/solver.h>
using namespace std;
index 9843b78d73c9af44efa1f6641f7a2342eda79370..cca30291d4f1eb4e1a9f238af82998b0b4c746a0 100644 (file)
W_{q,q} = 1 (n even)
*/
-#include "sbasis-to-bezier.h"
-#include "choose.h"
-#include "svg-path.h"
+#include <2geom/sbasis-to-bezier.h>
+#include <2geom/choose.h>
+#include <2geom/svg-path.h>
#include <iostream>
-#include "exception.h"
+#include <2geom/exception.h>
namespace Geom{
index b515ebf8bea141cbc4128a1e6950ee2e621cac0c..046e8299660748f7fe46a86bc7169b62af47cd0f 100644 (file)
#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 d7a3972e1709662586bdf1b792ad73307081368b..377238d9265997e8a52d2f9d0a9e5a9dcafb6f87 100644 (file)
--- a/src/2geom/sbasis.cpp
+++ b/src/2geom/sbasis.cpp
#include <cmath>
-#include "sbasis.h"
-#include "isnan.h"
+#include <2geom/sbasis.h>
+#include <2geom/isnan.h>
namespace Geom{
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]);
}
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 d6f55c8d8c0a8dde929e6bd43f85578f43b957a7..a9939f92a8a2e525d77d65d88acf6d816394e3ff 100644 (file)
--- a/src/2geom/sbasis.h
+++ b/src/2geom/sbasis.h
#include <cassert>
#include <iostream>
-#include "linear.h"
-#include "interval.h"
-#include "utils.h"
-#include "exception.h"
+#include <2geom/linear.h>
+#include <2geom/interval.h>
+#include <2geom/utils.h>
+#include <2geom/exception.h>
namespace Geom {
diff --git a/src/2geom/shape.cpp b/src/2geom/shape.cpp
index a7a40066e9b3a514bfcba85715a48a1acc9b3b56..936dfaf70fdb10ba6656dd66e342b2f57cc07a23 100644 (file)
--- a/src/2geom/shape.cpp
+++ b/src/2geom/shape.cpp
-#include "shape.h"
-#include "utils.h"
-#include "sweep.h"
-#include "ord.h"
+#include <2geom/shape.h>
+#include <2geom/utils.h>
+#include <2geom/sweep.h>
+#include <2geom/ord.h>
#include <iostream>
#include <algorithm>
diff --git a/src/2geom/shape.h b/src/2geom/shape.h
index a55fe9ad3457b95eb8d873355c0b4620af072ee6..147aa9a5064cae21eb95e9b1b73e36676acba596 100644 (file)
--- a/src/2geom/shape.h
+++ b/src/2geom/shape.h
#include <vector>
#include <set>
-#include "region.h"
+#include <2geom/region.h>
//TODO: BBOX optimizations
index 3ec738df1d1be1d65a7d58d4ac33dd0911ab6150..b922e970effd478165526de22624591774d68dff 100644 (file)
-#include "solver.h"
-#include "point.h"
+#include <2geom/solver.h>
+#include <2geom/point.h>
#include <algorithm>
/*** Find the zeros of the bernstein function. The code subdivides until it is happy with the
index fc49dcdb2fcdf565380b16e7263bc0165d190e1f..ad017f596d5258b7adcb50a23586cb65017b880a 100644 (file)
-#include "solver.h"
-#include "point.h"
+#include <2geom/solver.h>
+#include <2geom/point.h>
#include <algorithm>
namespace Geom{
diff --git a/src/2geom/solver.h b/src/2geom/solver.h
index 81ab431a26c882949b8b1ff52677083de2d6de1c..87365ab336553926d8a51e016e4ed37e0c917980 100644 (file)
--- a/src/2geom/solver.h
+++ b/src/2geom/solver.h
#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 9937a6f1517ece1bbfead5d1d58037e50327b36b..097a5120a1d3ee219625acb86cc24066950a1bc7 100644 (file)
--- a/src/2geom/sturm.h
+++ b/src/2geom/sturm.h
#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
--- /dev/null
@@ -0,0 +1,1124 @@
+/*\r
+ * SVG Elliptical Arc Class\r
+ *\r
+ * Copyright 2008 Marco Cecchetti <mrcekets at gmail.com>\r
+ *\r
+ * This library is free software; you can redistribute it and/or\r
+ * modify it either under the terms of the GNU Lesser General Public\r
+ * License version 2.1 as published by the Free Software Foundation\r
+ * (the "LGPL") or, at your option, under the terms of the Mozilla\r
+ * Public License Version 1.1 (the "MPL"). If you do not alter this\r
+ * notice, a recipient may use your version of this file under either\r
+ * the MPL or the LGPL.\r
+ *\r
+ * You should have received a copy of the LGPL along with this library\r
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software\r
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
+ * You should have received a copy of the MPL along with this library\r
+ * in the file COPYING-MPL-1.1\r
+ *\r
+ * The contents of this file are subject to the Mozilla Public License\r
+ * Version 1.1 (the "License"); you may not use this file except in\r
+ * compliance with the License. You may obtain a copy of the License at\r
+ * http://www.mozilla.org/MPL/\r
+ *\r
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY\r
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for\r
+ * the specific language governing rights and limitations.\r
+ */\r
+\r
+\r
+#include <2geom/svg-elliptical-arc.h>\r
+#include <2geom/ellipse.h>\r
+#include <2geom/sbasis-geometric.h>\r
+#include <2geom/bezier-curve.h>\r
+#include <2geom/poly.h>\r
+\r
+#include <cfloat>\r
+#include <limits>\r
+\r
+#include <2geom/numeric/vector.h>\r
+#include <2geom/numeric/fitting-tool.h>\r
+#include <2geom/numeric/fitting-model.h>\r
+\r
+\r
+\r
+namespace Geom\r
+{\r
+\r
+\r
+Rect SVGEllipticalArc::boundsExact() const\r
+{\r
+ if (isDegenerate() && is_svg_compliant())\r
+ return chord().boundsExact();\r
+\r
+ std::vector<double> extremes(4);\r
+ double cosrot = std::cos(rotation_angle());\r
+ double sinrot = std::sin(rotation_angle());\r
+ extremes[0] = std::atan2( -ray(Y) * sinrot, ray(X) * cosrot );\r
+ extremes[1] = extremes[0] + M_PI;\r
+ if ( extremes[0] < 0 ) extremes[0] += 2*M_PI;\r
+ extremes[2] = std::atan2( ray(Y) * cosrot, ray(X) * sinrot );\r
+ extremes[3] = extremes[2] + M_PI;\r
+ if ( extremes[2] < 0 ) extremes[2] += 2*M_PI;\r
+\r
+\r
+ std::vector<double>arc_extremes(4);\r
+ arc_extremes[0] = initialPoint()[X];\r
+ arc_extremes[1] = finalPoint()[X];\r
+ if ( arc_extremes[0] < arc_extremes[1] )\r
+ std::swap(arc_extremes[0], arc_extremes[1]);\r
+ arc_extremes[2] = initialPoint()[Y];\r
+ arc_extremes[3] = finalPoint()[Y];\r
+ if ( arc_extremes[2] < arc_extremes[3] )\r
+ std::swap(arc_extremes[2], arc_extremes[3]);\r
+\r
+\r
+ if ( start_angle() < end_angle() )\r
+ {\r
+ if ( sweep_flag() )\r
+ {\r
+ for ( unsigned int i = 0; i < extremes.size(); ++i )\r
+ {\r
+ if ( start_angle() < extremes[i] && extremes[i] < end_angle() )\r
+ {\r
+ arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];\r
+ }\r
+ }\r
+ }\r
+ else\r
+ {\r
+ for ( unsigned int i = 0; i < extremes.size(); ++i )\r
+ {\r
+ if ( start_angle() > extremes[i] || extremes[i] > end_angle() )\r
+ {\r
+ arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];\r
+ }\r
+ }\r
+ }\r
+ }\r
+ else\r
+ {\r
+ if ( sweep_flag() )\r
+ {\r
+ for ( unsigned int i = 0; i < extremes.size(); ++i )\r
+ {\r
+ if ( start_angle() < extremes[i] || extremes[i] < end_angle() )\r
+ {\r
+ arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];\r
+ }\r
+ }\r
+ }\r
+ else\r
+ {\r
+ for ( unsigned int i = 0; i < extremes.size(); ++i )\r
+ {\r
+ if ( start_angle() > extremes[i] && extremes[i] > end_angle() )\r
+ {\r
+ arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];\r
+ }\r
+ }\r
+ }\r
+ }\r
+\r
+ return Rect( Point(arc_extremes[1], arc_extremes[3]) ,\r
+ Point(arc_extremes[0], arc_extremes[2]) );\r
+}\r
+\r
+\r
+double SVGEllipticalArc::valueAtAngle(Coord t, Dim2 d) const\r
+{\r
+ double sin_rot_angle = std::sin(rotation_angle());\r
+ double cos_rot_angle = std::cos(rotation_angle());\r
+ if ( d == X )\r
+ {\r
+ return ray(X) * cos_rot_angle * std::cos(t)\r
+ - ray(Y) * sin_rot_angle * std::sin(t)\r
+ + center(X);\r
+ }\r
+ else if ( d == Y )\r
+ {\r
+ return ray(X) * sin_rot_angle * std::cos(t)\r
+ + ray(Y) * cos_rot_angle * std::sin(t)\r
+ + center(Y);\r
+ }\r
+ THROW_RANGEERROR("dimension parameter out of range");\r
+}\r
+\r
+\r
+std::vector<double>\r
+SVGEllipticalArc::roots(double v, Dim2 d) const\r
+{\r
+ if ( d > Y )\r
+ {\r
+ THROW_RANGEERROR("dimention out of range");\r
+ }\r
+\r
+ std::vector<double> sol;\r
+\r
+ if (isDegenerate() && is_svg_compliant())\r
+ {\r
+ return chord().roots(v, d);\r
+ }\r
+ else\r
+ {\r
+ if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )\r
+ {\r
+ if ( center(d) == v )\r
+ sol.push_back(0);\r
+ return sol;\r
+ }\r
+\r
+ const char* msg[2][2] =\r
+ {\r
+ { "d == X; ray(X) == 0; "\r
+ "s = (v - center(X)) / ( -ray(Y) * std::sin(rotation_angle()) ); "\r
+ "s should be contained in [-1,1]",\r
+ "d == X; ray(Y) == 0; "\r
+ "s = (v - center(X)) / ( ray(X) * std::cos(rotation_angle()) ); "\r
+ "s should be contained in [-1,1]"\r
+ },\r
+ { "d == Y; ray(X) == 0; "\r
+ "s = (v - center(X)) / ( ray(Y) * std::cos(rotation_angle()) ); "\r
+ "s should be contained in [-1,1]",\r
+ "d == Y; ray(Y) == 0; "\r
+ "s = (v - center(X)) / ( ray(X) * std::sin(rotation_angle()) ); "\r
+ "s should be contained in [-1,1]"\r
+ },\r
+ };\r
+\r
+ for ( unsigned int dim = 0; dim < 2; ++dim )\r
+ {\r
+ if ( are_near(ray(dim), 0) )\r
+ {\r
+ if ( initialPoint()[d] == v && finalPoint()[d] == v )\r
+ {\r
+ THROW_INFINITESOLUTIONS(0);\r
+ }\r
+ if ( (initialPoint()[d] < finalPoint()[d])\r
+ && (initialPoint()[d] > v || finalPoint()[d] < v) )\r
+ {\r
+ return sol;\r
+ }\r
+ if ( (initialPoint()[d] > finalPoint()[d])\r
+ && (finalPoint()[d] > v || initialPoint()[d] < v) )\r
+ {\r
+ return sol;\r
+ }\r
+ double ray_prj;\r
+ switch(d)\r
+ {\r
+ case X:\r
+ switch(dim)\r
+ {\r
+ case X: ray_prj = -ray(Y) * std::sin(rotation_angle());\r
+ break;\r
+ case Y: ray_prj = ray(X) * std::cos(rotation_angle());\r
+ break;\r
+ }\r
+ break;\r
+ case Y:\r
+ switch(dim)\r
+ {\r
+ case X: ray_prj = ray(Y) * std::cos(rotation_angle());\r
+ break;\r
+ case Y: ray_prj = ray(X) * std::sin(rotation_angle());\r
+ break;\r
+ }\r
+ break;\r
+ }\r
+\r
+ double s = (v - center(d)) / ray_prj;\r
+ if ( s < -1 || s > 1 )\r
+ {\r
+ THROW_LOGICALERROR(msg[d][dim]);\r
+ }\r
+ switch(dim)\r
+ {\r
+ case X:\r
+ s = std::asin(s); // return a value in [-PI/2,PI/2]\r
+ if ( logical_xor( sweep_flag(), are_near(start_angle(), M_PI/2) ) )\r
+ {\r
+ if ( s < 0 ) s += 2*M_PI;\r
+ }\r
+ else\r
+ {\r
+ s = M_PI - s;\r
+ if (!(s < 2*M_PI) ) s -= 2*M_PI;\r
+ }\r
+ break;\r
+ case Y:\r
+ s = std::acos(s); // return a value in [0,PI]\r
+ if ( logical_xor( sweep_flag(), are_near(start_angle(), 0) ) )\r
+ {\r
+ s = 2*M_PI - s;\r
+ if ( !(s < 2*M_PI) ) s -= 2*M_PI;\r
+ }\r
+ break;\r
+ }\r
+\r
+ //std::cerr << "s = " << rad_to_deg(s);\r
+ s = map_to_01(s);\r
+ //std::cerr << " -> t: " << s << std::endl;\r
+ if ( !(s < 0 || s > 1) )\r
+ sol.push_back(s);\r
+ return sol;\r
+ }\r
+ }\r
+\r
+ }\r
+\r
+ double rotx, roty;\r
+ switch(d)\r
+ {\r
+ case X:\r
+ rotx = std::cos(rotation_angle());\r
+ roty = -std::sin(rotation_angle());\r
+ break;\r
+ case Y:\r
+ rotx = std::sin(rotation_angle());\r
+ roty = std::cos(rotation_angle());\r
+ break;\r
+ }\r
+ double rxrotx = ray(X) * rotx;\r
+ double c_v = center(d) - v;\r
+\r
+ double a = -rxrotx + c_v;\r
+ double b = ray(Y) * roty;\r
+ double c = rxrotx + c_v;\r
+ //std::cerr << "a = " << a << std::endl;\r
+ //std::cerr << "b = " << b << std::endl;\r
+ //std::cerr << "c = " << c << std::endl;\r
+\r
+ if ( are_near(a,0) )\r
+ {\r
+ sol.push_back(M_PI);\r
+ if ( !are_near(b,0) )\r
+ {\r
+ double s = 2 * std::atan(-c/(2*b));\r
+ if ( s < 0 ) s += 2*M_PI;\r
+ sol.push_back(s);\r
+ }\r
+ }\r
+ else\r
+ {\r
+ double delta = b * b - a * c;\r
+ //std::cerr << "delta = " << delta << std::endl;\r
+ if ( are_near(delta, 0) )\r
+ {\r
+ double s = 2 * std::atan(-b/a);\r
+ if ( s < 0 ) s += 2*M_PI;\r
+ sol.push_back(s);\r
+ }\r
+ else if ( delta > 0 )\r
+ {\r
+ double sq = std::sqrt(delta);\r
+ double s = 2 * std::atan( (-b - sq) / a );\r
+ if ( s < 0 ) s += 2*M_PI;\r
+ sol.push_back(s);\r
+ s = 2 * std::atan( (-b + sq) / a );\r
+ if ( s < 0 ) s += 2*M_PI;\r
+ sol.push_back(s);\r
+ }\r
+ }\r
+\r
+ std::vector<double> arc_sol;\r
+ for (unsigned int i = 0; i < sol.size(); ++i )\r
+ {\r
+ //std::cerr << "s = " << rad_to_deg(sol[i]);\r
+ sol[i] = map_to_01(sol[i]);\r
+ //std::cerr << " -> t: " << sol[i] << std::endl;\r
+ if ( !(sol[i] < 0 || sol[i] > 1) )\r
+ arc_sol.push_back(sol[i]);\r
+ }\r
+ return arc_sol;\r
+}\r
+\r
+\r
+// D(E(t,C),t) = E(t+PI/2,O)\r
+Curve* SVGEllipticalArc::derivative() const\r
+{\r
+ if (isDegenerate() && is_svg_compliant())\r
+ return chord().derivative();\r
+\r
+ SVGEllipticalArc* result = new SVGEllipticalArc(*this);\r
+ result->m_center[X] = result->m_center[Y] = 0;\r
+ result->m_start_angle += M_PI/2;\r
+ if( !( result->m_start_angle < 2*M_PI ) )\r
+ {\r
+ result->m_start_angle -= 2*M_PI;\r
+ }\r
+ result->m_end_angle += M_PI/2;\r
+ if( !( result->m_end_angle < 2*M_PI ) )\r
+ {\r
+ result->m_end_angle -= 2*M_PI;\r
+ }\r
+ result->m_initial_point = result->pointAtAngle( result->start_angle() );\r
+ result->m_final_point = result->pointAtAngle( result->end_angle() );\r
+ return result;\r
+}\r
+\r
+\r
+std::vector<Point>\r
+SVGEllipticalArc::pointAndDerivatives(Coord t, unsigned int n) const\r
+{\r
+ if (isDegenerate() && is_svg_compliant())\r
+ return chord().pointAndDerivatives(t, n);\r
+\r
+ unsigned int nn = n+1; // nn represents the size of the result vector.\r
+ std::vector<Point> result;\r
+ result.reserve(nn);\r
+ double angle = map_unit_interval_on_circular_arc(t, start_angle(),\r
+ end_angle(), sweep_flag());\r
+ SVGEllipticalArc ea(*this);\r
+ ea.m_center = Point(0,0);\r
+ unsigned int m = std::min(nn, 4u);\r
+ for ( unsigned int i = 0; i < m; ++i )\r
+ {\r
+ result.push_back( ea.pointAtAngle(angle) );\r
+ angle += M_PI/2;\r
+ if ( !(angle < 2*M_PI) ) angle -= 2*M_PI;\r
+ }\r
+ m = nn / 4;\r
+ for ( unsigned int i = 1; i < m; ++i )\r
+ {\r
+ for ( unsigned int j = 0; j < 4; ++j )\r
+ result.push_back( result[j] );\r
+ }\r
+ m = nn - 4 * m;\r
+ for ( unsigned int i = 0; i < m; ++i )\r
+ {\r
+ result.push_back( result[i] );\r
+ }\r
+ if ( !result.empty() ) // nn != 0\r
+ result[0] = pointAtAngle(angle);\r
+ return result;\r
+}\r
+\r
+bool SVGEllipticalArc::containsAngle(Coord angle) const\r
+{\r
+ if ( sweep_flag() )\r
+ if ( start_angle() < end_angle() )\r
+ return ( !( angle < start_angle() || angle > end_angle() ) );\r
+ else\r
+ return ( !( angle < start_angle() && angle > end_angle() ) );\r
+ else\r
+ if ( start_angle() > end_angle() )\r
+ return ( !( angle > start_angle() || angle < end_angle() ) );\r
+ else\r
+ return ( !( angle > start_angle() && angle < end_angle() ) );\r
+}\r
+\r
+Curve* SVGEllipticalArc::portion(double f, double t) const\r
+{\r
+ if (f < 0) f = 0;\r
+ if (f > 1) f = 1;\r
+ if (t < 0) t = 0;\r
+ if (t > 1) t = 1;\r
+ if ( are_near(f, t) )\r
+ {\r
+ SVGEllipticalArc* arc = new SVGEllipticalArc();\r
+ arc->m_center = arc->m_initial_point = arc->m_final_point = pointAt(f);\r
+ arc->m_start_angle = arc->m_end_angle = m_start_angle;\r
+ arc->m_rot_angle = m_rot_angle;\r
+ arc->m_sweep = m_sweep;\r
+ arc->m_large_arc = m_large_arc;\r
+ }\r
+ SVGEllipticalArc* arc = new SVGEllipticalArc( *this );\r
+ arc->m_initial_point = pointAt(f);\r
+ arc->m_final_point = pointAt(t);\r
+ double sa = sweep_flag() ? sweep_angle() : -sweep_angle();\r
+ arc->m_start_angle = m_start_angle + sa * f;\r
+ if ( !(arc->m_start_angle < 2*M_PI) )\r
+ arc->m_start_angle -= 2*M_PI;\r
+ if ( arc->m_start_angle < 0 )\r
+ arc->m_start_angle += 2*M_PI;\r
+ arc->m_end_angle = m_start_angle + sa * t;\r
+ if ( !(arc->m_end_angle < 2*M_PI) )\r
+ arc->m_end_angle -= 2*M_PI;\r
+ if ( arc->m_end_angle < 0 )\r
+ arc->m_end_angle += 2*M_PI;\r
+ if ( f > t ) arc->m_sweep = !sweep_flag();\r
+ if ( large_arc_flag() && (arc->sweep_angle() < M_PI) )\r
+ arc->m_large_arc = false;\r
+ return arc;\r
+}\r
+\r
+\r
+std::vector<double> SVGEllipticalArc::\r
+allNearestPoints( Point const& p, double from, double to ) const\r
+{\r
+ std::vector<double> result;\r
+ if (isDegenerate() && is_svg_compliant())\r
+ {\r
+ result.push_back( chord().nearestPoint(p, from, to) );\r
+ return result;\r
+ }\r
+\r
+ if ( from > to ) std::swap(from, to);\r
+ if ( from < 0 || to > 1 )\r
+ {\r
+ THROW_RANGEERROR("[from,to] interval out of range");\r
+ }\r
+\r
+ if ( ( are_near(ray(X), 0) && are_near(ray(Y), 0) ) || are_near(from, to) )\r
+ {\r
+ result.push_back(from);\r
+ return result;\r
+ }\r
+ else if ( are_near(ray(X), 0) || are_near(ray(Y), 0) )\r
+ {\r
+ LineSegment seg(pointAt(from), pointAt(to));\r
+ Point np = seg.pointAt( seg.nearestPoint(p) );\r
+ if ( are_near(ray(Y), 0) )\r
+ {\r
+ if ( are_near(rotation_angle(), M_PI/2)\r
+ || are_near(rotation_angle(), 3*M_PI/2) )\r
+ {\r
+ result = roots(np[Y], Y);\r
+ }\r
+ else\r
+ {\r
+ result = roots(np[X], X);\r
+ }\r
+ }\r
+ else\r
+ {\r
+ if ( are_near(rotation_angle(), M_PI/2)\r
+ || are_near(rotation_angle(), 3*M_PI/2) )\r
+ {\r
+ result = roots(np[X], X);\r
+ }\r
+ else\r
+ {\r
+ result = roots(np[Y], Y);\r
+ }\r
+ }\r
+ return result;\r
+ }\r
+ else if ( are_near(ray(X), ray(Y)) )\r
+ {\r
+ Point r = p - center();\r
+ if ( are_near(r, Point(0,0)) )\r
+ {\r
+ THROW_INFINITESOLUTIONS(0);\r
+ }\r
+ // TODO: implement case r != 0\r
+// Point np = ray(X) * unit_vector(r);\r
+// std::vector<double> solX = roots(np[X],X);\r
+// std::vector<double> solY = roots(np[Y],Y);\r
+// double t;\r
+// if ( are_near(solX[0], solY[0]) || are_near(solX[0], solY[1]))\r
+// {\r
+// t = solX[0];\r
+// }\r
+// else\r
+// {\r
+// t = solX[1];\r
+// }\r
+// if ( !(t < from || t > to) )\r
+// {\r
+// result.push_back(t);\r
+// }\r
+// else\r
+// {\r
+//\r
+// }\r
+ }\r
+\r
+ // solve the equation <D(E(t),t)|E(t)-p> == 0\r
+ // that provides min and max distance points\r
+ // on the ellipse E wrt the point p\r
+ // after the substitutions:\r
+ // cos(t) = (1 - s^2) / (1 + s^2)\r
+ // sin(t) = 2t / (1 + s^2)\r
+ // where s = tan(t/2)\r
+ // we get a 4th degree equation in s\r
+ /*\r
+ * ry s^4 ((-cy + py) Cos[Phi] + (cx - px) Sin[Phi]) +\r
+ * ry ((cy - py) Cos[Phi] + (-cx + px) Sin[Phi]) +\r
+ * 2 s^3 (rx^2 - ry^2 + (-cx + px) rx Cos[Phi] + (-cy + py) rx Sin[Phi]) +\r
+ * 2 s (-rx^2 + ry^2 + (-cx + px) rx Cos[Phi] + (-cy + py) rx Sin[Phi])\r
+ */\r
+\r
+ Point p_c = p - center();\r
+ double rx2_ry2 = (ray(X) - ray(Y)) * (ray(X) + ray(Y));\r
+ double cosrot = std::cos( rotation_angle() );\r
+ double sinrot = std::sin( rotation_angle() );\r
+ double expr1 = ray(X) * (p_c[X] * cosrot + p_c[Y] * sinrot);\r
+ Poly coeff;\r
+ coeff.resize(5);\r
+ coeff[4] = ray(Y) * ( p_c[Y] * cosrot - p_c[X] * sinrot );\r
+ coeff[3] = 2 * ( rx2_ry2 + expr1 );\r
+ coeff[2] = 0;\r
+ coeff[1] = 2 * ( -rx2_ry2 + expr1 );\r
+ coeff[0] = -coeff[4];\r
+\r
+// for ( unsigned int i = 0; i < 5; ++i )\r
+// std::cerr << "c[" << i << "] = " << coeff[i] << std::endl;\r
+\r
+ std::vector<double> real_sol;\r
+ // gsl_poly_complex_solve raises an error\r
+ // if the leading coefficient is zero\r
+ if ( are_near(coeff[4], 0) )\r
+ {\r
+ real_sol.push_back(0);\r
+ if ( !are_near(coeff[3], 0) )\r
+ {\r
+ double sq = -coeff[1] / coeff[3];\r
+ if ( sq > 0 )\r
+ {\r
+ double s = std::sqrt(sq);\r
+ real_sol.push_back(s);\r
+ real_sol.push_back(-s);\r
+ }\r
+ }\r
+ }\r
+ else\r
+ {\r
+ real_sol = solve_reals(coeff);\r
+ }\r
+\r
+ for ( unsigned int i = 0; i < real_sol.size(); ++i )\r
+ {\r
+ real_sol[i] = 2 * std::atan(real_sol[i]);\r
+ if ( real_sol[i] < 0 ) real_sol[i] += 2*M_PI;\r
+ }\r
+ // when s -> Infinity then <D(E)|E-p> -> 0 iff coeff[4] == 0\r
+ // so we add M_PI to the solutions being lim arctan(s) = PI when s->Infinity\r
+ if ( (real_sol.size() % 2) != 0 )\r
+ {\r
+ real_sol.push_back(M_PI);\r
+ }\r
+\r
+ double mindistsq1 = std::numeric_limits<double>::max();\r
+ double mindistsq2 = std::numeric_limits<double>::max();\r
+ double dsq;\r
+ unsigned int mi1, mi2;\r
+ for ( unsigned int i = 0; i < real_sol.size(); ++i )\r
+ {\r
+ dsq = distanceSq(p, pointAtAngle(real_sol[i]));\r
+ if ( mindistsq1 > dsq )\r
+ {\r
+ mindistsq2 = mindistsq1;\r
+ mi2 = mi1;\r
+ mindistsq1 = dsq;\r
+ mi1 = i;\r
+ }\r
+ else if ( mindistsq2 > dsq )\r
+ {\r
+ mindistsq2 = dsq;\r
+ mi2 = i;\r
+ }\r
+ }\r
+\r
+ double t = map_to_01( real_sol[mi1] );\r
+ if ( !(t < from || t > to) )\r
+ {\r
+ result.push_back(t);\r
+ }\r
+\r
+ bool second_sol = false;\r
+ t = map_to_01( real_sol[mi2] );\r
+ if ( real_sol.size() == 4 && !(t < from || t > to) )\r
+ {\r
+ if ( result.empty() || are_near(mindistsq1, mindistsq2) )\r
+ {\r
+ result.push_back(t);\r
+ second_sol = true;\r
+ }\r
+ }\r
+\r
+ // we need to test extreme points too\r
+ double dsq1 = distanceSq(p, pointAt(from));\r
+ double dsq2 = distanceSq(p, pointAt(to));\r
+ if ( second_sol )\r
+ {\r
+ if ( mindistsq2 > dsq1 )\r
+ {\r
+ result.clear();\r
+ result.push_back(from);\r
+ mindistsq2 = dsq1;\r
+ }\r
+ else if ( are_near(mindistsq2, dsq) )\r
+ {\r
+ result.push_back(from);\r
+ }\r
+ if ( mindistsq2 > dsq2 )\r
+ {\r
+ result.clear();\r
+ result.push_back(to);\r
+ }\r
+ else if ( are_near(mindistsq2, dsq2) )\r
+ {\r
+ result.push_back(to);\r
+ }\r
+\r
+ }\r
+ else\r
+ {\r
+ if ( result.empty() )\r
+ {\r
+ if ( are_near(dsq1, dsq2) )\r
+ {\r
+ result.push_back(from);\r
+ result.push_back(to);\r
+ }\r
+ else if ( dsq2 > dsq1 )\r
+ {\r
+ result.push_back(from);\r
+ }\r
+ else\r
+ {\r
+ result.push_back(to);\r
+ }\r
+ }\r
+ }\r
+\r
+ return result;\r
+}\r
+\r
+\r
+/*\r
+ * NOTE: this implementation follows Standard SVG 1.1 implementation guidelines\r
+ * for elliptical arc curves. See Appendix F.6.\r
+ */\r
+void SVGEllipticalArc::calculate_center_and_extreme_angles()\r
+{\r
+ Point d = initialPoint() - finalPoint();\r
+\r
+ if (is_svg_compliant())\r
+ {\r
+ if ( initialPoint() == finalPoint() )\r
+ {\r
+ m_rx = m_ry = m_rot_angle = m_start_angle = m_end_angle = 0;\r
+ m_center = initialPoint();\r
+ m_large_arc = m_sweep = false;\r
+ return;\r
+ }\r
+\r
+ m_rx = std::fabs(m_rx);\r
+ m_ry = std::fabs(m_ry);\r
+\r
+ if ( are_near(ray(X), 0) || are_near(ray(Y), 0) )\r
+ {\r
+ m_rx = L2(d) / 2;\r
+ m_ry = 0;\r
+ m_rot_angle = std::atan2(d[Y], d[X]);\r
+ if (m_rot_angle < 0) m_rot_angle += 2*M_PI;\r
+ m_start_angle = 0;\r
+ m_end_angle = M_PI;\r
+ m_center = middle_point(initialPoint(), finalPoint());\r
+ m_large_arc = false;\r
+ m_sweep = false;\r
+ return;\r
+ }\r
+ }\r
+ else\r
+ {\r
+ if ( are_near(initialPoint(), finalPoint()) )\r
+ {\r
+ if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )\r
+ {\r
+ m_start_angle = m_end_angle = 0;\r
+ m_center = initialPoint();\r
+ return;\r
+ }\r
+ else\r
+ {\r
+ THROW_RANGEERROR("initial and final point are the same");\r
+ }\r
+ }\r
+ if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )\r
+ { // but initialPoint != finalPoint\r
+ THROW_RANGEERROR(\r
+ "there is no ellipse that satisfies the given constraints: "\r
+ "ray(X) == 0 && ray(Y) == 0 but initialPoint != finalPoint"\r
+ );\r
+ }\r
+ if ( are_near(ray(Y), 0) )\r
+ {\r
+ Point v = initialPoint() - finalPoint();\r
+ if ( are_near(L2sq(v), 4*ray(X)*ray(X)) )\r
+ {\r
+ double angle = std::atan2(v[Y], v[X]);\r
+ if (angle < 0) angle += 2*M_PI;\r
+ if ( are_near( angle, rotation_angle() ) )\r
+ {\r
+ m_start_angle = 0;\r
+ m_end_angle = M_PI;\r
+ m_center = v/2 + finalPoint();\r
+ return;\r
+ }\r
+ angle -= M_PI;\r
+ if ( angle < 0 ) angle += 2*M_PI;\r
+ if ( are_near( angle, rotation_angle() ) )\r
+ {\r
+ m_start_angle = M_PI;\r
+ m_end_angle = 0;\r
+ m_center = v/2 + finalPoint();\r
+ return;\r
+ }\r
+ THROW_RANGEERROR(\r
+ "there is no ellipse that satisfies the given constraints: "\r
+ "ray(Y) == 0 "\r
+ "and slope(initialPoint - finalPoint) != rotation_angle "\r
+ "and != rotation_angle + PI"\r
+ );\r
+ }\r
+ if ( L2sq(v) > 4*ray(X)*ray(X) )\r
+ {\r
+ THROW_RANGEERROR(\r
+ "there is no ellipse that satisfies the given constraints: "\r
+ "ray(Y) == 0 and distance(initialPoint, finalPoint) > 2*ray(X)"\r
+ );\r
+ }\r
+ else\r
+ {\r
+ THROW_RANGEERROR(\r
+ "there is infinite ellipses that satisfy the given constraints: "\r
+ "ray(Y) == 0 and distance(initialPoint, finalPoint) < 2*ray(X)"\r
+ );\r
+ }\r
+\r
+ }\r
+\r
+ if ( are_near(ray(X), 0) )\r
+ {\r
+ Point v = initialPoint() - finalPoint();\r
+ if ( are_near(L2sq(v), 4*ray(Y)*ray(Y)) )\r
+ {\r
+ double angle = std::atan2(v[Y], v[X]);\r
+ if (angle < 0) angle += 2*M_PI;\r
+ double rot_angle = rotation_angle() + M_PI/2;\r
+ if ( !(rot_angle < 2*M_PI) ) rot_angle -= 2*M_PI;\r
+ if ( are_near( angle, rot_angle ) )\r
+ {\r
+ m_start_angle = M_PI/2;\r
+ m_end_angle = 3*M_PI/2;\r
+ m_center = v/2 + finalPoint();\r
+ return;\r
+ }\r
+ angle -= M_PI;\r
+ if ( angle < 0 ) angle += 2*M_PI;\r
+ if ( are_near( angle, rot_angle ) )\r
+ {\r
+ m_start_angle = 3*M_PI/2;\r
+ m_end_angle = M_PI/2;\r
+ m_center = v/2 + finalPoint();\r
+ return;\r
+ }\r
+ THROW_RANGEERROR(\r
+ "there is no ellipse that satisfies the given constraints: "\r
+ "ray(X) == 0 "\r
+ "and slope(initialPoint - finalPoint) != rotation_angle + PI/2 "\r
+ "and != rotation_angle + (3/2)*PI"\r
+ );\r
+ }\r
+ if ( L2sq(v) > 4*ray(Y)*ray(Y) )\r
+ {\r
+ THROW_RANGEERROR(\r
+ "there is no ellipse that satisfies the given constraints: "\r
+ "ray(X) == 0 and distance(initialPoint, finalPoint) > 2*ray(Y)"\r
+ );\r
+ }\r
+ else\r
+ {\r
+ THROW_RANGEERROR(\r
+ "there is infinite ellipses that satisfy the given constraints: "\r
+ "ray(X) == 0 and distance(initialPoint, finalPoint) < 2*ray(Y)"\r
+ );\r
+ }\r
+\r
+ }\r
+\r
+ }\r
+\r
+ double sin_rot_angle = std::sin(rotation_angle());\r
+ double cos_rot_angle = std::cos(rotation_angle());\r
+\r
+\r
+ Matrix m( cos_rot_angle, -sin_rot_angle,\r
+ sin_rot_angle, cos_rot_angle,\r
+ 0, 0 );\r
+\r
+ Point p = (d / 2) * m;\r
+ double rx2 = m_rx * m_rx;\r
+ double ry2 = m_ry * m_ry;\r
+ double rxpy = m_rx * p[Y];\r
+ double rypx = m_ry * p[X];\r
+ double rx2py2 = rxpy * rxpy;\r
+ double ry2px2 = rypx * rypx;\r
+ double num = rx2 * ry2;\r
+ double den = rx2py2 + ry2px2;\r
+ assert(den != 0);\r
+ double rad = num / den;\r
+ Point c(0,0);\r
+ if (rad > 1)\r
+ {\r
+ rad -= 1;\r
+ rad = std::sqrt(rad);\r
+\r
+ if (m_large_arc == m_sweep) rad = -rad;\r
+ c = rad * Point(rxpy / m_ry, -rypx / m_rx);\r
+\r
+ m[1] = -m[1];\r
+ m[2] = -m[2];\r
+\r
+ m_center = c * m + middle_point(initialPoint(), finalPoint());\r
+ }\r
+ else if (rad == 1 || is_svg_compliant())\r
+ {\r
+ double lamda = std::sqrt(1 / rad);\r
+ m_rx *= lamda;\r
+ m_ry *= lamda;\r
+ m_center = middle_point(initialPoint(), finalPoint());\r
+ }\r
+ else\r
+ {\r
+ THROW_RANGEERROR(\r
+ "there is no ellipse that satisfies the given constraints"\r
+ );\r
+ }\r
+\r
+ Point sp((p[X] - c[X]) / m_rx, (p[Y] - c[Y]) / m_ry);\r
+ Point ep((-p[X] - c[X]) / m_rx, (-p[Y] - c[Y]) / m_ry);\r
+ Point v(1, 0);\r
+ m_start_angle = angle_between(v, sp);\r
+ double sweep_angle = angle_between(sp, ep);\r
+ if (!m_sweep && sweep_angle > 0) sweep_angle -= 2*M_PI;\r
+ if (m_sweep && sweep_angle < 0) sweep_angle += 2*M_PI;\r
+\r
+ if (m_start_angle < 0) m_start_angle += 2*M_PI;\r
+ m_end_angle = m_start_angle + sweep_angle;\r
+ if (m_end_angle < 0) m_end_angle += 2*M_PI;\r
+ if (m_end_angle >= 2*M_PI) m_end_angle -= 2*M_PI;\r
+}\r
+\r
+\r
+D2<SBasis> SVGEllipticalArc::toSBasis() const\r
+{\r
+\r
+ if (isDegenerate() && is_svg_compliant())\r
+ return chord().toSBasis();\r
+\r
+ D2<SBasis> arc;\r
+ // the interval of parametrization has to be [0,1]\r
+ Coord et = start_angle() + ( sweep_flag() ? sweep_angle() : -sweep_angle() );\r
+ Linear param(start_angle(), et);\r
+ Coord cos_rot_angle = std::cos(rotation_angle());\r
+ Coord sin_rot_angle = std::sin(rotation_angle());\r
+ // order = 4 seems to be enough to get a perfect looking elliptical arc\r
+ // should it be choosen in function of the arc length anyway ?\r
+ // or maybe a user settable parameter: toSBasis(unsigned int order) ?\r
+ SBasis arc_x = ray(X) * cos(param,4);\r
+ SBasis arc_y = ray(Y) * sin(param,4);\r
+ arc[0] = arc_x * cos_rot_angle - arc_y * sin_rot_angle + Linear(center(X),center(X));\r
+ arc[1] = arc_x * sin_rot_angle + arc_y * cos_rot_angle + Linear(center(Y),center(Y));\r
+ return arc;\r
+}\r
+\r
+\r
+Coord SVGEllipticalArc::map_to_02PI(Coord t) const\r
+{\r
+ if ( sweep_flag() )\r
+ {\r
+ Coord angle = start_angle() + sweep_angle() * t;\r
+ if ( !(angle < 2*M_PI) )\r
+ angle -= 2*M_PI;\r
+ return angle;\r
+ }\r
+ else\r
+ {\r
+ Coord angle = start_angle() - sweep_angle() * t;\r
+ if ( angle < 0 ) angle += 2*M_PI;\r
+ return angle;\r
+ }\r
+}\r
+\r
+Coord SVGEllipticalArc::map_to_01(Coord angle) const\r
+{\r
+ return map_circular_arc_on_unit_interval(angle, start_angle(),\r
+ end_angle(), sweep_flag());\r
+}\r
+\r
+\r
+namespace detail\r
+{\r
+\r
+struct ellipse_equation\r
+{\r
+ ellipse_equation(double a, double b, double c, double d, double e, double f)\r
+ : A(a), B(b), C(c), D(d), E(e), F(f)\r
+ {\r
+ }\r
+\r
+ double operator()(double x, double y) const\r
+ {\r
+ // A * x * x + B * x * y + C * y * y + D * x + E * y + F;\r
+ return (A * x + B * y + D) * x + (C * y + E) * y + F;\r
+ }\r
+\r
+ double operator()(Point const& p) const\r
+ {\r
+ return (*this)(p[X], p[Y]);\r
+ }\r
+\r
+ Point normal(double x, double y) const\r
+ {\r
+ Point n( 2 * A * x + B * y + D, 2 * C * y + B * x + E );\r
+ return unit_vector(n);\r
+ }\r
+\r
+ Point normal(Point const& p) const\r
+ {\r
+ return normal(p[X], p[Y]);\r
+ }\r
+\r
+ double A, B, C, D, E, F;\r
+};\r
+\r
+}\r
+\r
+make_elliptical_arc::\r
+make_elliptical_arc( SVGEllipticalArc& _ea,\r
+ curve_type const& _curve,\r
+ unsigned int _total_samples,\r
+ double _tolerance )\r
+ : ea(_ea), curve(_curve),\r
+ dcurve( unitVector(derivative(curve)) ),\r
+ model(), fitter(model, _total_samples),\r
+ tolerance(_tolerance), tol_at_extr(tolerance/2),\r
+ tol_at_center(0.1), angle_tol(0.1),\r
+ initial_point(curve.at0()), final_point(curve.at1()),\r
+ N(_total_samples), last(N-1), partitions(N-1), p(N),\r
+ svg_compliant(true)\r
+{\r
+}\r
+\r
+bool\r
+make_elliptical_arc::\r
+bound_exceeded( unsigned int k, detail::ellipse_equation const & ee,\r
+ double e1x, double e1y, double e2 )\r
+{\r
+ dist_err = std::fabs( ee(p[k]) );\r
+ dist_bound = std::fabs( e1x * p[k][X] + e1y * p[k][Y] + e2 );\r
+ angle_err = std::fabs( dot( dcurve(k/partitions), ee.normal(p[k]) ) );\r
+ //angle_err *= angle_err;\r
+ return ( dist_err > dist_bound || angle_err > angle_tol );\r
+}\r
+\r
+bool\r
+make_elliptical_arc::\r
+check_bound(double A, double B, double C, double D, double E, double F)\r
+{\r
+ // check error magnitude\r
+ detail::ellipse_equation ee(A, B, C, D, E, F);\r
+\r
+ double e1x = (2*A + B) * tol_at_extr;\r
+ double e1y = (B + 2*C) * tol_at_extr;\r
+ double e2 = ((D + E) + (A + B + C) * tol_at_extr) * tol_at_extr;\r
+ if ( bound_exceeded(0, ee, e1x, e1y, e2) )\r
+ {\r
+ print_bound_error(0);\r
+ return false;\r
+ }\r
+ if ( bound_exceeded(0, ee, e1x, e1y, e2) )\r
+ {\r
+ print_bound_error(last);\r
+ return false;\r
+ }\r
+\r
+ e1x = (2*A + B) * tolerance;\r
+ e1y = (B + 2*C) * tolerance;\r
+ e2 = ((D + E) + (A + B + C) * tolerance) * tolerance;\r
+// std::cerr << "e1x = " << e1x << std::endl;\r
+// std::cerr << "e1y = " << e1y << std::endl;\r
+// std::cerr << "e2 = " << e2 << std::endl;\r
+\r
+ for ( unsigned int k = 1; k < last; ++k )\r
+ {\r
+ if ( bound_exceeded(k, ee, e1x, e1y, e2) )\r
+ {\r
+ print_bound_error(k);\r
+ return false;\r
+ }\r
+ }\r
+\r
+ return true;\r
+}\r
+\r
+void make_elliptical_arc::fit()\r
+{\r
+ for (unsigned int k = 0; k < N; ++k)\r
+ {\r
+ p[k] = curve( k / partitions );\r
+ fitter.append(p[k]);\r
+ }\r
+ fitter.update();\r
+\r
+ NL::Vector z(N, 0.0);\r
+ fitter.result(z);\r
+}\r
+\r
+bool make_elliptical_arc::make_elliptiarc()\r
+{\r
+ const NL::Vector & coeff = fitter.result();\r
+ Ellipse e;\r
+ try\r
+ {\r
+ e.set(1, coeff[0], coeff[1], coeff[2], coeff[3], coeff[4]);\r
+ }\r
+ catch(LogicalError exc)\r
+ {\r
+ return false;\r
+ }\r
+\r
+ Point inner_point = curve(0.5);\r
+\r
+ if (svg_compliant_flag())\r
+ {\r
+ ea = e.arc(initial_point, inner_point, final_point);\r
+ }\r
+ else\r
+ {\r
+ try\r
+ {\r
+ ea = e.arc(initial_point, inner_point, final_point, false);\r
+ }\r
+ catch(RangeError exc)\r
+ {\r
+ return false;\r
+ }\r
+ }\r
+\r
+\r
+ if ( !are_near( e.center(),\r
+ ea.center(),\r
+ tol_at_center * std::min(e.ray(X),e.ray(Y))\r
+ )\r
+ )\r
+ {\r
+ return false;\r
+ }\r
+ return true;\r
+}\r
+\r
+\r
+\r
+} // end namespace Geom\r
+\r
+\r
+\r
+\r
+/*\r
+ Local Variables:\r
+ mode:c++\r
+ c-file-style:"stroustrup"\r
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))\r
+ indent-tabs-mode:nil\r
+ fill-column:99\r
+ End:\r
+*/\r
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :\r
+\r
diff --git a/src/2geom/svg-elliptical-arc.h b/src/2geom/svg-elliptical-arc.h
--- /dev/null
@@ -0,0 +1,435 @@
+/*\r
+ * Elliptical Arc - implementation of the svg elliptical arc path element\r
+ *\r
+ * Authors:\r
+ * MenTaLguY <mental@rydia.net>\r
+ * Marco Cecchetti <mrcekets at gmail.com>\r
+ *\r
+ * Copyright 2007-2008 authors\r
+ *\r
+ * This library is free software; you can redistribute it and/or\r
+ * modify it either under the terms of the GNU Lesser General Public\r
+ * License version 2.1 as published by the Free Software Foundation\r
+ * (the "LGPL") or, at your option, under the terms of the Mozilla\r
+ * Public License Version 1.1 (the "MPL"). If you do not alter this\r
+ * notice, a recipient may use your version of this file under either\r
+ * the MPL or the LGPL.\r
+ *\r
+ * You should have received a copy of the LGPL along with this library\r
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software\r
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
+ * You should have received a copy of the MPL along with this library\r
+ * in the file COPYING-MPL-1.1\r
+ *\r
+ * The contents of this file are subject to the Mozilla Public License\r
+ * Version 1.1 (the "License"); you may not use this file except in\r
+ * compliance with the License. You may obtain a copy of the License at\r
+ * http://www.mozilla.org/MPL/\r
+ *\r
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY\r
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for\r
+ * the specific language governing rights and limitations.\r
+ */\r
+\r
+\r
+#ifndef _2GEOM_SVG_ELLIPTICAL_ARC_H_\r
+#define _2GEOM_SVG_ELLIPTICAL_ARC_H_\r
+\r
+\r
+#include <2geom/curve.h>\r
+#include <2geom/angle.h>\r
+#include <2geom/utils.h>\r
+#include <2geom/bezier-curve.h>\r
+#include <2geom/sbasis-curve.h> // for non-native methods\r
+#include <2geom/numeric/vector.h>\r
+#include <2geom/numeric/fitting-tool.h>\r
+#include <2geom/numeric/fitting-model.h>\r
+\r
+\r
+#include <algorithm>\r
+\r
+\r
+\r
+namespace Geom\r
+{\r
+\r
+class SVGEllipticalArc : public Curve\r
+{\r
+ public:\r
+ SVGEllipticalArc(bool _svg_compliant = true)\r
+ : m_initial_point(Point(0,0)), m_final_point(Point(0,0)),\r
+ m_rx(0), m_ry(0), m_rot_angle(0),\r
+ m_large_arc(true), m_sweep(true),\r
+ m_svg_compliant(_svg_compliant)\r
+ {\r
+ m_start_angle = m_end_angle = 0;\r
+ m_center = Point(0,0);\r
+ }\r
+\r
+ SVGEllipticalArc( Point _initial_point, double _rx, double _ry,\r
+ double _rot_angle, bool _large_arc, bool _sweep,\r
+ Point _final_point,\r
+ bool _svg_compliant = true\r
+ )\r
+ : m_initial_point(_initial_point), m_final_point(_final_point),\r
+ m_rx(_rx), m_ry(_ry), m_rot_angle(_rot_angle),\r
+ m_large_arc(_large_arc), m_sweep(_sweep),\r
+ m_svg_compliant(_svg_compliant)\r
+ {\r
+ calculate_center_and_extreme_angles();\r
+ }\r
+\r
+ void set( Point _initial_point, double _rx, double _ry,\r
+ double _rot_angle, bool _large_arc, bool _sweep,\r
+ Point _final_point\r
+ )\r
+ {\r
+ m_initial_point = _initial_point;\r
+ m_final_point = _final_point;\r
+ m_rx = _rx;\r
+ m_ry = _ry;\r
+ m_rot_angle = _rot_angle;\r
+ m_large_arc = _large_arc;\r
+ m_sweep = _sweep;\r
+ calculate_center_and_extreme_angles();\r
+ }\r
+\r
+ Curve* duplicate() const\r
+ {\r
+ return new SVGEllipticalArc(*this);\r
+ }\r
+\r
+ double center(unsigned int i) const\r
+ {\r
+ return m_center[i];\r
+ }\r
+\r
+ Point center() const\r
+ {\r
+ return m_center;\r
+ }\r
+\r
+ Point initialPoint() const\r
+ {\r
+ return m_initial_point;\r
+ }\r
+\r
+ Point finalPoint() const\r
+ {\r
+ return m_final_point;\r
+ }\r
+\r
+ double start_angle() const\r
+ {\r
+ return m_start_angle;\r
+ }\r
+\r
+ double end_angle() const\r
+ {\r
+ return m_end_angle;\r
+ }\r
+\r
+ double ray(unsigned int i) const\r
+ {\r
+ return (i == 0) ? m_rx : m_ry;\r
+ }\r
+\r
+ bool large_arc_flag() const\r
+ {\r
+ return m_large_arc;\r
+ }\r
+\r
+ bool sweep_flag() const\r
+ {\r
+ return m_sweep;\r
+ }\r
+\r
+ double rotation_angle() const\r
+ {\r
+ return m_rot_angle;\r
+ }\r
+\r
+ void setInitial( const Point _point)\r
+ {\r
+ m_initial_point = _point;\r
+ calculate_center_and_extreme_angles();\r
+ }\r
+\r
+ void setFinal( const Point _point)\r
+ {\r
+ m_final_point = _point;\r
+ calculate_center_and_extreme_angles();\r
+ }\r
+\r
+ void setExtremes( const Point& _initial_point, const Point& _final_point )\r
+ {\r
+ m_initial_point = _initial_point;\r
+ m_final_point = _final_point;\r
+ calculate_center_and_extreme_angles();\r
+ }\r
+\r
+ bool isDegenerate() const\r
+ {\r
+ return ( are_near(ray(X), 0) || are_near(ray(Y), 0) );\r
+ }\r
+\r
+ bool is_svg_compliant() const\r
+ {\r
+ return m_svg_compliant;\r
+ }\r
+\r
+ Rect boundsFast() const\r
+ {\r
+ return boundsExact();\r
+ }\r
+\r
+ Rect boundsExact() const;\r
+\r
+ // TODO: native implementation of the following methods\r
+ Rect boundsLocal(Interval i, unsigned int deg) const\r
+ {\r
+ if (isDegenerate() && is_svg_compliant())\r
+ return chord().boundsLocal(i, deg);\r
+ else\r
+ return SBasisCurve(toSBasis()).boundsLocal(i, deg);\r
+ }\r
+\r
+ std::vector<double> roots(double v, Dim2 d) const;\r
+\r
+ std::vector<double>\r
+ allNearestPoints( Point const& p, double from = 0, double to = 1 ) const;\r
+\r
+ double nearestPoint( Point const& p, double from = 0, double to = 1 ) const\r
+ {\r
+ if ( are_near(ray(X), ray(Y)) && are_near(center(), p) )\r
+ {\r
+ return from;\r
+ }\r
+ return allNearestPoints(p, from, to).front();\r
+ }\r
+\r
+ // TODO: native implementation of the following methods\r
+ int winding(Point p) const\r
+ {\r
+ if (isDegenerate() && is_svg_compliant())\r
+ return chord().winding(p);\r
+ else\r
+ return SBasisCurve(toSBasis()).winding(p);\r
+ }\r
+\r
+ Curve *derivative() const;\r
+\r
+ // TODO: native implementation of the following methods\r
+ Curve *transformed(Matrix const &m) const\r
+ {\r
+ return SBasisCurve(toSBasis()).transformed(m);\r
+ }\r
+\r
+ std::vector<Point> pointAndDerivatives(Coord t, unsigned int n) const;\r
+\r
+ D2<SBasis> toSBasis() const;\r
+\r
+ bool containsAngle(Coord angle) const;\r
+\r
+ double valueAtAngle(Coord t, Dim2 d) const;\r
+\r
+ Point pointAtAngle(Coord t) const\r
+ {\r
+ double sin_rot_angle = std::sin(rotation_angle());\r
+ double cos_rot_angle = std::cos(rotation_angle());\r
+ Matrix m( ray(X) * cos_rot_angle, ray(X) * sin_rot_angle,\r
+ -ray(Y) * sin_rot_angle, ray(Y) * cos_rot_angle,\r
+ center(X), center(Y) );\r
+ Point p( std::cos(t), std::sin(t) );\r
+ return p * m;\r
+ }\r
+\r
+ double valueAt(Coord t, Dim2 d) const\r
+ {\r
+ if (isDegenerate() && is_svg_compliant())\r
+ return chord().valueAt(t, d);\r
+\r
+ Coord tt = map_to_02PI(t);\r
+ return valueAtAngle(tt, d);\r
+ }\r
+\r
+ Point pointAt(Coord t) const\r
+ {\r
+ if (isDegenerate() && is_svg_compliant())\r
+ return chord().pointAt(t);\r
+\r
+ Coord tt = map_to_02PI(t);\r
+ return pointAtAngle(tt);\r
+ }\r
+\r
+ std::pair<SVGEllipticalArc, SVGEllipticalArc>\r
+ subdivide(Coord t) const\r
+ {\r
+ SVGEllipticalArc* arc1 = static_cast<SVGEllipticalArc*>(portion(0, t));\r
+ SVGEllipticalArc* arc2 = static_cast<SVGEllipticalArc*>(portion(t, 1));\r
+ assert( arc1 != NULL && arc2 != NULL);\r
+ std::pair<SVGEllipticalArc, SVGEllipticalArc> arc_pair(*arc1, *arc2);\r
+ delete arc1;\r
+ delete arc2;\r
+ return arc_pair;\r
+ }\r
+\r
+ Curve* portion(double f, double t) const;\r
+\r
+ // the arc is the same but traversed in the opposite direction\r
+ Curve* reverse() const\r
+ {\r
+ SVGEllipticalArc* rarc = new SVGEllipticalArc( *this );\r
+ rarc->m_sweep = !m_sweep;\r
+ rarc->m_initial_point = m_final_point;\r
+ rarc->m_final_point = m_initial_point;\r
+ rarc->m_start_angle = m_end_angle;\r
+ rarc->m_end_angle = m_start_angle;\r
+ return rarc;\r
+ }\r
+\r
+ double sweep_angle() const\r
+ {\r
+ Coord d = end_angle() - start_angle();\r
+ if ( !sweep_flag() ) d = -d;\r
+ if ( d < 0 )\r
+ d += 2*M_PI;\r
+ return d;\r
+ }\r
+\r
+ LineSegment chord() const\r
+ {\r
+ return LineSegment(initialPoint(), finalPoint());\r
+ }\r
+\r
+ private:\r
+ Coord map_to_02PI(Coord t) const;\r
+ Coord map_to_01(Coord angle) const;\r
+ void calculate_center_and_extreme_angles();\r
+\r
+ private:\r
+ Point m_initial_point, m_final_point;\r
+ double m_rx, m_ry, m_rot_angle;\r
+ bool m_large_arc, m_sweep;\r
+ double m_start_angle, m_end_angle;\r
+ Point m_center;\r
+ bool m_svg_compliant;\r
+\r
+}; // end class SVGEllipticalArc\r
+\r
+template< class charT >\r
+inline\r
+std::basic_ostream<charT> &\r
+operator<< (std::basic_ostream<charT> & os, const SVGEllipticalArc & ea)\r
+{\r
+ os << "{ cx: " << ea.center(X) << ", cy: " << ea.center(Y)\r
+ << ", rx: " << ea.ray(X) << ", ry: " << ea.ray(Y)\r
+ << ", rot angle: " << decimal_round(rad_to_deg(ea.rotation_angle()),2)\r
+ << ", start angle: " << decimal_round(rad_to_deg(ea.start_angle()),2)\r
+ << ", end angle: " << decimal_round(rad_to_deg(ea.end_angle()),2)\r
+ << " }";\r
+\r
+ return os;\r
+}\r
+\r
+\r
+\r
+\r
+namespace detail\r
+{\r
+ struct ellipse_equation;\r
+}\r
+\r
+\r
+class make_elliptical_arc\r
+{\r
+ public:\r
+ typedef D2<SBasis> curve_type;\r
+\r
+ make_elliptical_arc( SVGEllipticalArc& _ea,\r
+ curve_type const& _curve,\r
+ unsigned int _total_samples,\r
+ double _tolerance );\r
+\r
+ private:\r
+ bool bound_exceeded( unsigned int k, detail::ellipse_equation const & ee,\r
+ double e1x, double e1y, double e2 );\r
+\r
+ bool check_bound(double A, double B, double C, double D, double E, double F);\r
+\r
+ void fit();\r
+\r
+ bool make_elliptiarc();\r
+\r
+ void print_bound_error(unsigned int k)\r
+ {\r
+ std::cerr\r
+ << "tolerance error" << std::endl\r
+ << "at point: " << k << std::endl\r
+ << "error value: "<< dist_err << std::endl\r
+ << "bound: " << dist_bound << std::endl\r
+ << "angle error: " << angle_err\r
+ << " (" << angle_tol << ")" << std::endl;\r
+ }\r
+\r
+ public:\r
+ bool operator()()\r
+ {\r
+ const NL::Vector & coeff = fitter.result();\r
+ fit();\r
+ if ( !check_bound(1, coeff[0], coeff[1], coeff[2], coeff[3], coeff[4]) )\r
+ return false;\r
+ if ( !(make_elliptiarc()) ) return false;\r
+ return true;\r
+ }\r
+\r
+ bool svg_compliant_flag() const\r
+ {\r
+ return svg_compliant;\r
+ }\r
+\r
+ void svg_compliant_on()\r
+ {\r
+ svg_compliant = true;\r
+ }\r
+\r
+ void svg_compliant_off()\r
+ {\r
+ svg_compliant = false;\r
+ }\r
+\r
+ private:\r
+ SVGEllipticalArc& ea;\r
+ const curve_type & curve;\r
+ Piecewise<D2<SBasis> > dcurve;\r
+ NL::LFMEllipse model;\r
+ NL::least_squeares_fitter<NL::LFMEllipse> fitter;\r
+ double tolerance, tol_at_extr, tol_at_center, angle_tol;\r
+ Point initial_point, final_point;\r
+ unsigned int N;\r
+ unsigned int last; // N-1\r
+ double partitions; // N-1\r
+ std::vector<Point> p; // sample points\r
+ double dist_err, dist_bound, angle_err;\r
+ bool svg_compliant;\r
+};\r
+\r
+\r
+} // end namespace Geom\r
+\r
+\r
+\r
+\r
+#endif /* _2GEOM_SVG_ELLIPTICAL_ARC_H_ */\r
+\r
+/*\r
+ Local Variables:\r
+ mode:c++\r
+ c-file-style:"stroustrup"\r
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))\r
+ indent-tabs-mode:nil\r
+ fill-column:99\r
+ End:\r
+*/\r
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :\r
+\r
index fe996cef109f486989fe8002b8bf966d6b7c8393..2a0ee9b1f5222a4784192e503f61fec654ce1bb1 100644 (file)
-#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
*
#include <vector>
#include <glib.h>
-#include "point.h"
-
-#include "svg-path-parser.h"
+#include <2geom/point.h>
+#include <2geom/svg-path-parser.h>
namespace Geom {
};
-#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,
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)
_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;
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);
}
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();
}
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();
}
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();
}
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();
}
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"
}
}
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 ) {
index 98a885361dff03cbeb925abedd855a248aa93821..d24d255517d094c6a4b83f5e00d1f1b9ab198242 100644 (file)
#include <vector>
#include <iterator>
#include <stdexcept>
-#include "exception.h"
-#include "point.h"
-#include "svg-path.h"
+#include <2geom/exception.h>
+#include <2geom/point.h>
+#include <2geom/svg-path.h>
namespace Geom {
diff --git a/src/2geom/svg-path.cpp b/src/2geom/svg-path.cpp
index 819b53aacfc27a7699f2f4762a7e5d250ecb3318..c178b92da43900ffc86f5d906521006ca1573c62 100644 (file)
--- a/src/2geom/svg-path.cpp
+++ b/src/2geom/svg-path.cpp
*
*/
-#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 ad78680deae589d58f3b46d9f3a9a8720080fffa..1ae4b054cba2bb821b248ce163d951b4db056990 100644 (file)
--- a/src/2geom/svg-path.h
+++ b/src/2geom/svg-path.h
#ifndef SEEN_SVG_PATH_H
#define SEEN_SVG_PATH_H
-#include "path.h"
+#include <2geom/path.h>
#include <iterator>
namespace Geom {
diff --git a/src/2geom/sweep.cpp b/src/2geom/sweep.cpp
index 92da524ce076b73f3ff2eb19dd553036d77c57d7..1b881399036103801c60a7cfb02e44d757db34f4 100644 (file)
--- a/src/2geom/sweep.cpp
+++ b/src/2geom/sweep.cpp
-#include "sweep.h"
+#include <2geom/sweep.h>
#include <algorithm>
diff --git a/src/2geom/sweep.h b/src/2geom/sweep.h
index 94ea6855c7f225f48261b94fde87fc4cec11d059..6500ccbc2c880ff3706cf1974a67b7b8ccd88a04 100644 (file)
--- a/src/2geom/sweep.h
+++ b/src/2geom/sweep.h
#define __2GEOM_SWEEP_H__
#include <vector>
-#include "d2.h"
+#include <2geom/d2.h>
namespace Geom {
index 62b340221f5072f37f012448e411f9662128b416..a6426fe811a978e94e5d60524e41d6176ebae9c3 100644 (file)
--- a/src/2geom/transforms.cpp
+++ b/src/2geom/transforms.cpp
-#include "transforms.h"
+#include <2geom/transforms.h>
namespace Geom {
diff --git a/src/2geom/transforms.h b/src/2geom/transforms.h
index fb077a910316256d2313e337185ca070ba1782a7..ca0a81ac7d5d53122b0892506408a04b7b88cae1 100644 (file)
--- a/src/2geom/transforms.h
+++ b/src/2geom/transforms.h
#ifndef SEEN_Geom_TRANSFORMS_H
#define SEEN_Geom_TRANSFORMS_H
-#include "matrix.h"
+#include <2geom/matrix.h>
#include <cmath>
namespace Geom {
diff --git a/src/2geom/utils.h b/src/2geom/utils.h
index 96c63fddd530598cae85818ae53fb35da565797e..ac99b8b007cc4ea2cca51dd8c294e9c8b2e1368a 100644 (file)
--- a/src/2geom/utils.h
+++ b/src/2geom/utils.h
*/
#include <cmath>
+#include <vector>
namespace Geom {
return round( x * multiplier ) / multiplier;
}
+
+void binomial_coefficients(std::vector<size_t>& bc, size_t n);
+
}
#endif
index 282ea82706f44738eaa1d974fde3367564afdfd8..dc98561e1b4e6500a6d759f4a5d8e02f68f7d14a 100644 (file)
Piecewise<D2<SBasis> > n = rot90(derivative(uskeleton));
n = force_continuity(remove_short_cuts(n,.1));
- D2<Piecewise<SBasis> > patternd2 = make_cuts_independant(pwd2_in);
+ D2<Piecewise<SBasis> > patternd2 = make_cuts_independent(pwd2_in);
Piecewise<SBasis> x = vertical_pattern.get_value() ? Piecewise<SBasis>(patternd2[1]) : Piecewise<SBasis>(patternd2[0]);
Piecewise<SBasis> y = vertical_pattern.get_value() ? Piecewise<SBasis>(patternd2[0]) : Piecewise<SBasis>(patternd2[1]);
index 814a381af6ede23ae2c1914369576d9786d7608b..34fb60a31fabf2669fa6c541247be43e288970ac 100644 (file)
startpoint_spacing_variation.resetRandomizer();
endpoint_spacing_variation.resetRandomizer();
- D2<Piecewise<SBasis> > stroke = make_cuts_independant(strokepath.get_pwd2());
+ D2<Piecewise<SBasis> > stroke = make_cuts_independent(strokepath.get_pwd2());
Interval bndsStroke = bounds_exact(stroke[0]);
gdouble scaling = bndsStroke.max() - bndsStroke.min();
Interval bndsStrokeY = bounds_exact(stroke[1]);
for (unsigned int i=0; i < temppath.size(); i++) {
pwd2.concat( temppath[i].toPwSb() );
}
- D2<Piecewise<SBasis> > d2pw = make_cuts_independant(pwd2);
+ D2<Piecewise<SBasis> > d2pw = make_cuts_independent(pwd2);
Interval bndsX = bounds_exact(d2pw[0]);
Interval bndsY = bounds_exact(d2pw[1]);
Point start(bndsX.min(), (bndsY.max()+bndsY.min())/2);
index 9ebe97b9f490c9bd9cd6fef31b498f1f06ecb94a..7d2045d80a95bb478c534221faf535720f300a53 100755 (executable)
n4 = force_continuity(remove_short_cuts(n4,.1));
- D2<Piecewise<SBasis> > patternd2 = make_cuts_independant(pwd2_in);
+ D2<Piecewise<SBasis> > patternd2 = make_cuts_independent(pwd2_in);
Piecewise<SBasis> x = Piecewise<SBasis>(patternd2[0]);
Piecewise<SBasis> y = Piecewise<SBasis>(patternd2[1]);
index 85de609fbf3421a29e525ec4c24147d9710dc132..2e10efe7a0b90bfb2de9ea7f2056dbfd8e063525 100644 (file)
PAPCopyType type = copytype.get_value();
- D2<Piecewise<SBasis> > patternd2 = make_cuts_independant(pattern.get_pwd2());
+ D2<Piecewise<SBasis> > patternd2 = make_cuts_independent(pattern.get_pwd2());
Piecewise<SBasis> x0 = vertical_pattern.get_value() ? Piecewise<SBasis>(patternd2[1]) : Piecewise<SBasis>(patternd2[0]);
Piecewise<SBasis> y0 = vertical_pattern.get_value() ? Piecewise<SBasis>(patternd2[0]) : Piecewise<SBasis>(patternd2[1]);
Interval pattBndsX = bounds_exact(x0);
index 785d7a245ccdbbf473cbc4028500adc2418864c1..9aa26e05d0097e08714aea16ca5cd9bdfcdee308 100644 (file)
// remove this once we have unified coordinate systems
path_a_pw = path_a_pw + Geom::Point(offsetx, -offsety);
- D2<Piecewise<SBasis> > B = make_cuts_independant(path_a_pw);
+ D2<Piecewise<SBasis> > B = make_cuts_independent(path_a_pw);
Piecewise<SBasis> preimage[4];
//Geom::Point orig = Geom::Point(bounds_X.min(), bounds_Y.middle());
index 114dd726306165610eec3397a3334457998927a7..e448682909149d022d4385141fb7ba8fccf54f36 100644 (file)
pwd2.concat( temppath[i].toPwSb() );
}
- D2<Piecewise<SBasis> > d2pw = make_cuts_independant(pwd2);
+ D2<Piecewise<SBasis> > d2pw = make_cuts_independent(pwd2);
Interval bndsX = bounds_exact(d2pw[0]);
Interval bndsY = bounds_exact(d2pw[1]);
Point start(bndsX.min(), (bndsY.max()+bndsY.min())/2);