![[tokkee]](http://tokkee.org/images/avatar.png){kind=avatar}
summary | shortlog | log | commit | commitdiff | tree
raw | patch | inline | side by side (parent: c62b0ca)
raw | patch | inline | side by side (parent: c62b0ca)
| author | johanengelen <johanengelen@users.sourceforge.net> | |
| Sat, 10 May 2008 20:20:11 +0000 (20:20 +0000) | ||
| committer | johanengelen <johanengelen@users.sourceforge.net> | |
| Sat, 10 May 2008 20:20:11 +0000 (20:20 +0000) |
27 files changed:
diff --git a/build.xml b/build.xml
index 5e38e1792bda062b876e9772877fd9cd203712d6..49fa4d76a8f4f6db7b2f3fa664f238cdf68b069b 100644 (file)
--- a/build.xml
+++ b/build.xml
-lgdk_pixbuf-2.0
-lpangocairo-1.0 -lpangoft2-1.0 -lpangowin32-1.0 -lpango-1.0
-lgobject-2.0 -lgmodule-2.0 -lgthread-2.0 -lglib-2.0
+ -lgsl
<!-- PERL -->
-L${gtk}/perl/lib/CORE -lperl58
<!-- PYTHON -->
index f70293a13789abda6107ce0a69cd4fd84aaf2e58..ef82b75251f8ec6a577a7351528b09442440ad80 100644 (file)
2geom/linear.h \\r
2geom/matrix.cpp \\r
2geom/matrix.h \\r
+ 2geom/nearest-point.cpp \\r
+ 2geom/nearest-point.h \\r
+ 2geom/numeric/linear_system.h \\r
+ 2geom/numeric/matrix.h \\r
+ 2geom/numeric/vector.h \\r
2geom/ord.h \\r
2geom/path-intersection.cpp \\r
2geom/path-intersection.h \\r
diff --git a/src/2geom/angle.h b/src/2geom/angle.h
index c6a367d8f75f0045fd2507134bc956cd37f1f973..e95300a174f6dc87575c569af654619b8006a370 100644 (file)
--- a/src/2geom/angle.h
+++ b/src/2geom/angle.h
inline double rad_to_deg(double rad) { return rad*180.0/M_PI;}
+/*
+ * start_angle and angle must belong to [0, 2PI[
+ * and angle must belong to the cirsular arc defined by
+ * start_angle, end_angle and with rotation direction cw
+ */
+inline
+double map_circular_arc_on_unit_interval( double angle, double start_angle, double end_angle, bool cw = true )
+{
+ double d = end_angle - start_angle;
+ double t = angle - start_angle;
+ if ( !cw )
+ {
+ d = -d;
+ t = -t;
+ }
+ if ( d < 0 ) d += 2*M_PI;
+ if ( t < 0 ) t += 2*M_PI;
+ return t / d;
+}
+
+inline
+Coord map_unit_interval_on_circular_arc(Coord t, double start_angle, double end_angle, bool cw = true)
+{
+ double sweep_angle = end_angle - start_angle;
+ if ( !cw ) sweep_angle = -sweep_angle;
+ if ( sweep_angle < 0 ) sweep_angle += 2*M_PI;
+
+ if ( cw )
+ {
+ Coord angle = start_angle + sweep_angle * t;
+ if ( !(angle < 2*M_PI) )
+ angle -= 2*M_PI;
+ return angle;
+ }
+ else
+ {
+ Coord angle = start_angle - sweep_angle * t;
+ if ( angle < 0 ) angle += 2*M_PI;
+ return angle;
+ }
+}
+
+
}
#endif
index 28b3c6f2091ea3ce678e1c05ab8e51087e80bc2e..97c4c6e5ce78d92a4bc2a5dfb123b299c3cc7eae 100644 (file)
std::vector<std::pair<double, double> >
find_self_intersections(OldBezier const &Sb) {
- throwNotImplemented();
+ THROW_NOTIMPLEMENTED();
}
std::vector<std::pair<double, double> >
index 76c90a9150df722fb5054aff63c7ca392dc3379f..59aac895117fe4ab53665636dc655be62a46320d 100644 (file)
@@ -165,8 +165,8 @@ copy_without_nans_or_adjacent_duplicates(Point const src[], unsigned src_len, Po
if ( si == src_len ) {
return 0;
}
- if (!is_nan(src[si][X]) &&
- !is_nan(src[si][Y])) {
+ if (!IS_NAN(src[si][X]) &&
+ !IS_NAN(src[si][Y])) {
dest[0] = Point(src[si]);
++si;
break;
@@ -176,8 +176,8 @@ copy_without_nans_or_adjacent_duplicates(Point const src[], unsigned src_len, Po
for (; si < src_len; ++si) {
Point const src_pt = Point(src[si]);
if ( src_pt != dest[di]
- && !is_nan(src_pt[X])
- && !is_nan(src_pt[Y])) {
+ && !IS_NAN(src_pt[X])
+ && !IS_NAN(src_pt[Y])) {
dest[++di] = src_pt;
}
}
bezier[0] = data[0];
bezier[3] = data[len - 1];
double const dist = distance(bezier[0], bezier[3]) / 3.0;
- if (is_nan(dist)) {
+ if (IS_NAN(dist)) {
/* Numerical problem, fall back to straight line segment. */
bezier[1] = bezier[0];
bezier[2] = bezier[3];
}
}
- if (!is_finite(improved_u)) {
+ if (!IS_FINITE(improved_u)) {
improved_u = u;
} else if ( improved_u < 0.0 ) {
improved_u = 0.0;
double tot_len = u[len - 1];
if(!( tot_len != 0 ))
return;
- if (is_finite(tot_len)) {
+ if (IS_FINITE(tot_len)) {
for (unsigned i = 1; i < len; ++i) {
u[i] /= tot_len;
}
diff --git a/src/2geom/bezier.h b/src/2geom/bezier.h
index 9b7d8fb17f88e28dadb11a066d99993fb5086bd3..289a6772953e69337ccf4df040ce3934c8216061 100644 (file)
--- a/src/2geom/bezier.h
+++ b/src/2geom/bezier.h
}
inline bool isFinite() const {
for(unsigned i = 0; i <= order(); i++) {
- if(!is_finite(c_[i])) return false;
+ if(!IS_FINITE(c_[i])) return false;
}
return true;
}
diff --git a/src/2geom/exception.h b/src/2geom/exception.h
index 88ecfc51b4694fd478595773fe1dcf2eeac5b182..2749ec63a4515645bc39b4bc4a8c0644d4313f58 100644 (file)
--- a/src/2geom/exception.h
+++ b/src/2geom/exception.h
namespace Geom {
-// Base exception class, all 2geom exceptions should be derrived from this one.
+// Base exception class, all 2geom exceptions should be derived from this one.
class Exception : public std::exception {
public:
Exception(const char * message, const char *file, const int line) {
protected:
std::string msgstr;
};
-#define throwException(message) throw(Geom::Exception(message, __FILE__, __LINE__))
+#define THROW_EXCEPTION(message) throw(Geom::Exception(message, __FILE__, __LINE__))
//-----------------------------------------------------------------------
// Two main exception classes: LogicalError and RangeError.
LogicalError(const char * message, const char *file, const int line)
: Exception(message, file, line) {}
};
-#define throwLogicalError(message) throw(LogicalError(message, __FILE__, __LINE__))
+#define THROW_LOGICALERROR(message) throw(LogicalError(message, __FILE__, __LINE__))
class RangeError : public Exception {
public:
RangeError(const char * message, const char *file, const int line)
: Exception(message, file, line) {}
};
-#define throwRangeError(message) throw(RangeError(message, __FILE__, __LINE__))
+#define THROW_RANGEERROR(message) throw(RangeError(message, __FILE__, __LINE__))
//-----------------------------------------------------------------------
// Special case exceptions. Best used with the defines :)
NotImplemented(const char *file, const int line)
: LogicalError("Method not implemented", file, line) {}
};
-#define throwNotImplemented(i) throw(NotImplemented(__FILE__, __LINE__))
+#define THROW_NOTIMPLEMENTED(i) throw(NotImplemented(__FILE__, __LINE__))
class InvariantsViolation : public LogicalError {
public:
InvariantsViolation(const char *file, const int line)
: LogicalError("Invariants violation", file, line) {}
};
-#define throwInvariantsViolation(i) throw(InvariantsViolation(__FILE__, __LINE__))
-#define assert_invariants(e) ((e) ? (void)0 : throwInvariantsViolation())
+#define THROW_INVARIANTSVIOLATION(i) throw(InvariantsViolation(__FILE__, __LINE__))
+#define ASSERT_INVARIANTS(e) ((e) ? (void)0 : THROW_INVARIANTSVIOLATION())
class NotInvertible : public RangeError {
public:
NotInvertible(const char *file, const int line)
: RangeError("Function does not have a unique inverse", file, line) {}
};
-#define throwNotInvertible(i) throw(NotInvertible(__FILE__, __LINE__))
+#define THROW_NOTINVERTIBLE(i) throw(NotInvertible(__FILE__, __LINE__))
class ContinuityError : public RangeError {
public:
ContinuityError(const char *file, const int line)
: RangeError("Non-contiguous path", file, line) {}
};
-#define throwContinuityError(i) throw(ContinuityError(__FILE__, __LINE__))
+#define THROW_CONTINUITYERROR(i) throw(ContinuityError(__FILE__, __LINE__))
struct SVGPathParseError : public std::exception {
char const *what() const throw() { return "parse error"; }
diff --git a/src/2geom/isnan.h b/src/2geom/isnan.h
index 5b068e606da7cb19b21953af4d40abdd79717762..decebc7d27ad41410ecb7baa7a0e6248c28c724e 100644 (file)
--- a/src/2geom/isnan.h
+++ b/src/2geom/isnan.h
*/
#if defined(__isnan)
-# define is_nan(_a) (__isnan(_a))
+# define IS_NAN(_a) (__isnan(_a))
#elif defined(__APPLE__) && __GNUC__ == 3
-# define is_nan(_a) (__isnan(_a)) /* MacOSX/Darwin definition < 10.4 */
+# define IS_NAN(_a) (__isnan(_a)) /* MacOSX/Darwin definition < 10.4 */
#elif defined(WIN32) || defined(_isnan)
-# define is_nan(_a) (_isnan(_a)) /* Win32 definition */
+# define IS_NAN(_a) (_isnan(_a)) /* Win32 definition */
#elif defined(isnan) || defined(__FreeBSD__)
-# define is_nan(_a) (isnan(_a)) /* GNU definition */
+# define IS_NAN(_a) (isnan(_a)) /* GNU definition */
#else
-# define is_nan(_a) (std::isnan(_a))
+# define IS_NAN(_a) (std::isnan(_a))
#endif
/* If the above doesn't work, then try (a != a).
* Also, please report a bug as per http://www.inkscape.org/report_bugs.php,
#if defined(__isfinite)
-# define is_finite(_a) (__isfinite(_a))
+# define IS_FINITE(_a) (__isfinite(_a))
#elif defined(__APPLE__) && __GNUC__ == 3
-# define is_finite(_a) (__isfinite(_a)) /* MacOSX/Darwin definition < 10.4 */
+# define IS_FINITE(_a) (__isfinite(_a)) /* MacOSX/Darwin definition < 10.4 */
#elif defined(isfinite)
-# define is_finite(_a) (isfinite(_a))
+# define IS_FINITE(_a) (isfinite(_a))
#else
-# define is_finite(_a) (std::isfinite(_a))
+# define IS_FINITE(_a) (std::isfinite(_a))
#endif
/* If the above doesn't work, then try (finite(_a) && !isNaN(_a)) or (!isNaN((_a) - (_a))).
* Also, please report a bug as per http://www.inkscape.org/report_bugs.php,
diff --git a/src/2geom/linear.h b/src/2geom/linear.h
index 0778039d33841c0ff434134fc5c683caf724e87d..93cbc1b0413946b4570ba1b6f0fc9e7d400e81ba 100644 (file)
--- a/src/2geom/linear.h
+++ b/src/2geom/linear.h
typedef double output_type;
inline bool isZero() const { return a[0] == 0 && a[1] == 0; }
inline bool isConstant() const { return a[0] == a[1]; }
- inline bool isFinite() const { return is_finite(a[0]) && is_finite(a[1]); }
+ inline bool isFinite() const { return IS_FINITE(a[0]) && IS_FINITE(a[1]); }
inline double at0() const { return a[0]; }
inline double at1() const { return a[1]; }
diff --git a/src/2geom/nearest-point.cpp b/src/2geom/nearest-point.cpp
--- /dev/null
@@ -0,0 +1,278 @@
+/*\r
+ * nearest point routines for D2<SBasis> and Piecewise<D2<SBasis>>\r
+ *\r
+ * Authors:\r
+ * \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
+#include "nearest-point.h"\r
+\r
+namespace Geom\r
+{\r
+\r
+////////////////////////////////////////////////////////////////////////////////\r
+// D2<SBasis> versions\r
+\r
+/*\r
+ * Return the parameter t of a nearest point on the portion of the curve "c", \r
+ * related to the interval [from, to], to the point "p".\r
+ * The needed curve derivative "dc" is passed as parameter.\r
+ * The function return the first nearest point to "p" that is found.\r
+ */\r
+\r
+double nearest_point( Point const& p, \r
+ D2<SBasis> const& c, \r
+ D2<SBasis> const& dc, \r
+ double from, double to )\r
+{\r
+ if ( from > to ) std::swap(from, to);\r
+ if ( from < 0 || to > 1 )\r
+ {\r
+ THROW_RANGEERROR("[from,to] interval out of bounds");\r
+ }\r
+\r
+ SBasis dd = dot(c - p, dc); \r
+ std::vector<double> zeros = Geom::roots(dd);\r
+ \r
+ double closest = from;\r
+ double min_dist_sq = L2sq(c(from) - p);\r
+ double distsq;\r
+ for ( unsigned int i = 0; i < zeros.size(); ++i )\r
+ {\r
+ distsq = L2sq(c(zeros[i]) - p);\r
+ if ( min_dist_sq > L2sq(c(zeros[i]) - p) )\r
+ {\r
+ closest = zeros[i];\r
+ min_dist_sq = distsq;\r
+ }\r
+ }\r
+ if ( min_dist_sq > L2sq( c(to) - p ) )\r
+ closest = to;\r
+ return closest;\r
+\r
+}\r
+\r
+/*\r
+ * Return the parameters t of all the nearest points on the portion of \r
+ * the curve "c", related to the interval [from, to], to the point "p".\r
+ * The needed curve derivative "dc" is passed as parameter.\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
+{\r
+ std::swap(from, to);\r
+ if ( from > to ) std::swap(from, to);\r
+ if ( from < 0 || to > 1 )\r
+ {\r
+ THROW_RANGEERROR("[from,to] interval out of bounds");\r
+ }\r
+\r
+ std::vector<double> result;\r
+ SBasis dd = dot(c - p, Geom::derivative(c));\r
+ \r
+ std::vector<double> zeros = Geom::roots(dd);\r
+ std::vector<double> candidates;\r
+ candidates.push_back(from);\r
+ candidates.insert(candidates.end(), zeros.begin(), zeros.end());\r
+ candidates.push_back(to);\r
+ std::vector<double> distsq;\r
+ distsq.reserve(candidates.size());\r
+ for ( unsigned int i = 0; i < candidates.size(); ++i )\r
+ {\r
+ distsq.push_back( L2sq(c(candidates[i]) - p) );\r
+ }\r
+ unsigned int closest = 0;\r
+ double dsq = distsq[0];\r
+ for ( unsigned int i = 1; i < candidates.size(); ++i )\r
+ {\r
+ if ( dsq > distsq[i] )\r
+ {\r
+ closest = i;\r
+ dsq = distsq[i];\r
+ }\r
+ }\r
+ for ( unsigned int i = 0; i < candidates.size(); ++i )\r
+ {\r
+ if( distsq[closest] == distsq[i] )\r
+ {\r
+ result.push_back(candidates[i]);\r
+ }\r
+ }\r
+ return result;\r
+}\r
+\r
+\r
+////////////////////////////////////////////////////////////////////////////////\r
+// Piecewise< D2<SBasis> > versions\r
+\r
+\r
+double nearest_point( Point const& p, \r
+ Piecewise< D2<SBasis> > const& c,\r
+ double from, double to )\r
+{\r
+ if ( from > to ) std::swap(from, to);\r
+ if ( from < c.cuts[0] || to > c.cuts[c.size()] )\r
+ {\r
+ THROW_RANGEERROR("[from,to] interval out of bounds");\r
+ }\r
+ \r
+ unsigned int si = c.segN(from);\r
+ unsigned int ei = c.segN(to);\r
+ if ( si == ei )\r
+ {\r
+ double nearest=\r
+ nearest_point(p, c[si], c.segT(from, si), c.segT(to, si));\r
+ return c.mapToDomain(nearest, si);\r
+ }\r
+ double t;\r
+ double nearest = nearest_point(p, c[si], c.segT(from, si));\r
+ unsigned int ni = si;\r
+ double dsq;\r
+ double mindistsq = distanceSq(p, c[si](nearest));\r
+ Rect bb;\r
+ for ( unsigned int i = si + 1; i < ei; ++i )\r
+ {\r
+ bb = bounds_fast(c[i]);\r
+ dsq = distanceSq(p, bb);\r
+ if ( mindistsq <= dsq ) continue;\r
+ t = nearest_point(p, c[i]);\r
+ dsq = distanceSq(p, c[i](t));\r
+ if ( mindistsq > dsq )\r
+ {\r
+ nearest = t;\r
+ ni = i;\r
+ mindistsq = dsq;\r
+ }\r
+ }\r
+ bb = bounds_fast(c[ei]);\r
+ dsq = distanceSq(p, bb);\r
+ if ( mindistsq > dsq )\r
+ {\r
+ t = nearest_point(p, c[ei], 0, c.segT(to, ei));\r
+ dsq = distanceSq(p, c[ei](t));\r
+ if ( mindistsq > dsq )\r
+ {\r
+ nearest = t;\r
+ ni = ei;\r
+ }\r
+ }\r
+ return c.mapToDomain(nearest, ni);\r
+}\r
+\r
+std::vector<double> \r
+all_nearest_points( Point const& p, \r
+ Piecewise< D2<SBasis> > const& c, \r
+ double from, double to )\r
+{\r
+ if ( from > to ) std::swap(from, to);\r
+ if ( from < c.cuts[0] || to > c.cuts[c.size()] )\r
+ {\r
+ THROW_RANGEERROR("[from,to] interval out of bounds");\r
+ }\r
+ \r
+ unsigned int si = c.segN(from);\r
+ unsigned int ei = c.segN(to);\r
+ if ( si == ei )\r
+ {\r
+ std::vector<double> all_nearest = \r
+ all_nearest_points(p, c[si], c.segT(from, si), c.segT(to, si));\r
+ for ( unsigned int i = 0; i < all_nearest.size(); ++i )\r
+ {\r
+ all_nearest[i] = c.mapToDomain(all_nearest[i], si);\r
+ }\r
+ return all_nearest;\r
+ }\r
+ std::vector<double> all_t;\r
+ std::vector< std::vector<double> > all_np;\r
+ all_np.push_back( all_nearest_points(p, c[si], c.segT(from, si)) );\r
+ std::vector<unsigned int> ni;\r
+ ni.push_back(si);\r
+ double dsq;\r
+ double mindistsq = distanceSq( p, c[si](all_np.front().front()) );\r
+ Rect bb;\r
+ for ( unsigned int i = si + 1; i < ei; ++i )\r
+ {\r
+ bb = bounds_fast(c[i]);\r
+ dsq = distanceSq(p, bb);\r
+ if ( mindistsq < dsq ) continue;\r
+ all_t = all_nearest_points(p, c[i]);\r
+ dsq = distanceSq( p, c[i](all_t.front()) );\r
+ if ( mindistsq > dsq )\r
+ {\r
+ all_np.clear();\r
+ all_np.push_back(all_t);\r
+ ni.clear();\r
+ ni.push_back(i);\r
+ mindistsq = dsq;\r
+ }\r
+ else if ( mindistsq == dsq )\r
+ {\r
+ all_np.push_back(all_t);\r
+ ni.push_back(i);\r
+ }\r
+ }\r
+ bb = bounds_fast(c[ei]);\r
+ dsq = distanceSq(p, bb);\r
+ if ( mindistsq >= dsq )\r
+ {\r
+ all_t = all_nearest_points(p, c[ei], 0, c.segT(to, ei));\r
+ dsq = distanceSq( p, c[ei](all_t.front()) );\r
+ if ( mindistsq > dsq )\r
+ {\r
+ for ( unsigned int i = 0; i < all_t.size(); ++i )\r
+ {\r
+ all_t[i] = c.mapToDomain(all_t[i], ei);\r
+ }\r
+ return all_t;\r
+ }\r
+ else if ( mindistsq == dsq )\r
+ {\r
+ all_np.push_back(all_t);\r
+ ni.push_back(ei);\r
+ }\r
+ }\r
+ std::vector<double> all_nearest;\r
+ for ( unsigned int i = 0; i < all_np.size(); ++i )\r
+ {\r
+ for ( unsigned int j = 0; j < all_np[i].size(); ++j )\r
+ {\r
+ all_nearest.push_back( c.mapToDomain(all_np[i][j], ni[i]) );\r
+ }\r
+ }\r
+ return all_nearest;\r
+}\r
+\r
+} // end namespace Geom\r
+\r
+\r
diff --git a/src/2geom/nearest-point.h b/src/2geom/nearest-point.h
--- /dev/null
@@ -0,0 +1,122 @@
+/*\r
+ * nearest point routines for D2<SBasis> and Piecewise<D2<SBasis>>\r
+ *\r
+ *\r
+ * Authors:\r
+ * \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 _NEAREST_POINT_H_\r
+#define _NEAREST_POINT_H_\r
+\r
+\r
+#include <vector>\r
+\r
+#include "d2.h"\r
+#include "piecewise.h"\r
+#include "exception.h"\r
+\r
+\r
+\r
+namespace Geom\r
+{\r
+\r
+////////////////////////////////////////////////////////////////////////////////\r
+// D2<SBasis> versions\r
+\r
+/*\r
+ * Return the parameter t of a nearest point on the portion of the curve "c", \r
+ * related to the interval [from, to], to the point "p".\r
+ * The needed curve derivative "dc" is passed as parameter.\r
+ * The function return the first nearest point to "p" that is found.\r
+ */\r
+double nearest_point( Point const& p,\r
+ D2<SBasis> const& c, D2<SBasis> const& dc, \r
+ double from = 0, double to = 1 );\r
+\r
+inline\r
+double nearest_point( Point const& p, \r
+ D2<SBasis> const& c, \r
+ double from = 0, double to = 1 )\r
+{\r
+ return nearest_point(p, c, Geom::derivative(c), from, to);\r
+}\r
+\r
+/*\r
+ * Return the parameters t of all the nearest points on the portion of \r
+ * the curve "c", related to the interval [from, to], to the point "p".\r
+ * The needed curve derivative "dc" is passed as parameter.\r
+ */\r
+std::vector<double> \r
+all_nearest_points( Point const& p, \r
+ D2<SBasis> const& c, D2<SBasis> const& dc, \r
+ double from = 0, double to = 1 );\r
+\r
+inline\r
+std::vector<double> \r
+all_nearest_points( Point const& p, \r
+ D2<SBasis> const& c, \r
+ double from = 0, double to = 1 )\r
+{\r
+ return all_nearest_points(p, c, Geom::derivative(c), from, to);\r
+}\r
+\r
+\r
+////////////////////////////////////////////////////////////////////////////////\r
+// Piecewise< D2<SBasis> > versions\r
+\r
+double nearest_point( Point const& p, \r
+ Piecewise< D2<SBasis> > const& c, \r
+ double from, double to );\r
+\r
+inline\r
+double nearest_point( Point const& p, Piecewise< D2<SBasis> > const& c ) \r
+{\r
+ return nearest_point(p, c, c.cuts[0], c.cuts[c.size()]);\r
+}\r
+\r
+\r
+std::vector<double> \r
+all_nearest_points( Point const& p, \r
+ Piecewise< D2<SBasis> > const& c, \r
+ double from, double to );\r
+\r
+inline\r
+std::vector<double> \r
+all_nearest_points( Point const& p, Piecewise< D2<SBasis> > const& c ) \r
+{\r
+ return all_nearest_points(p, c, c.cuts[0], c.cuts[c.size()]);\r
+}\r
+\r
+} // end namespace Geom\r
+\r
+\r
+\r
+#endif /*_NEAREST_POINT_H_*/\r
diff --git a/src/2geom/numeric/linear_system.h b/src/2geom/numeric/linear_system.h
--- /dev/null
@@ -0,0 +1,89 @@
+#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 "numeric/matrix.h"\r
+#include "numeric/vector.h"\r
+\r
+\r
+namespace Geom { namespace NL {\r
+\r
+\r
+class LinearSystem\r
+{\r
+public:\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
+ Matrix & matrix()\r
+ {\r
+ return m_matrix;\r
+ }\r
+ \r
+ Vector & vector()\r
+ {\r
+ return m_vector;\r
+ }\r
+ \r
+ const Vector & solution() const\r
+ {\r
+ return m_solution;\r
+ }\r
+ \r
+private:\r
+ Matrix & m_matrix;\r
+ Vector & m_vector;\r
+ Vector m_solution;\r
+};\r
+\r
+\r
+} } // end namespaces\r
+\r
+\r
+#endif /*_NL_LINEAR_SYSTEM_H_*/\r
diff --git a/src/2geom/numeric/matrix.h b/src/2geom/numeric/matrix.h
--- /dev/null
@@ -0,0 +1,207 @@
+#ifndef _NL_MATRIX_H_\r
+#define _NL_MATRIX_H_\r
+\r
+\r
+#include <cassert>\r
+#include <utility>\r
+#include <algorithm>\r
+\r
+#include <gsl/gsl_matrix.h>\r
+\r
+\r
+\r
+namespace Geom { namespace NL {\r
+\r
+class Matrix;\r
+void swap(Matrix & m1, Matrix & m2);\r
+\r
+\r
+class Matrix\r
+{\r
+public:\r
+ // the matrix is not inizialized\r
+ Matrix(size_t n1, size_t n2)\r
+ : m_rows(n1), m_columns(n2)\r
+ {\r
+ m_matrix = gsl_matrix_alloc(n1, n2);\r
+ }\r
+ \r
+ Matrix(size_t n1, size_t n2, double x)\r
+ :m_rows(n1), m_columns(n2)\r
+ {\r
+ m_matrix = gsl_matrix_alloc(n1, n2);\r
+ gsl_matrix_set_all(m_matrix, x);\r
+ }\r
+ \r
+ Matrix( Matrix const& _matrix )\r
+ : m_rows(_matrix.rows()), m_columns(_matrix.columns())\r
+ {\r
+ m_matrix = gsl_matrix_alloc(rows(), columns());\r
+ gsl_matrix_memcpy(m_matrix, _matrix.m_matrix);\r
+ }\r
+ \r
+// explicit\r
+// Matrix( gsl_matrix* m, size_t n1, size_t n2)\r
+// : m_rows(n1), m_columns(n2), m_matrix(m)\r
+// { \r
+// }\r
+ \r
+ Matrix & operator=(Matrix const& _matrix)\r
+ {\r
+ assert( rows() == _matrix.rows() && columns() == _matrix.columns() );\r
+ gsl_matrix_memcpy(m_matrix, _matrix.m_matrix);\r
+ \r
+ return *this;\r
+ }\r
+ \r
+ virtual ~Matrix()\r
+ {\r
+ gsl_matrix_free(m_matrix);\r
+ }\r
+ \r
+ void set_all( double x = 0)\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
+ const double & operator() (size_t i, size_t j) const\r
+ {\r
+ return *gsl_matrix_const_ptr(m_matrix, i, j);\r
+ }\r
+ \r
+ double & operator() (size_t i, size_t j)\r
+ {\r
+ return *gsl_matrix_ptr(m_matrix, i, j);\r
+ }\r
+ \r
+ gsl_matrix* get_gsl_matrix()\r
+ {\r
+ return m_matrix;\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
+ Matrix & transpose()\r
+ {\r
+ gsl_matrix_transpose(m_matrix);\r
+ return (*this);\r
+ }\r
+ \r
+ Matrix & scale(double x)\r
+ {\r
+ gsl_matrix_scale(m_matrix, x);\r
+ return (*this);\r
+ }\r
+ \r
+ Matrix & translate(double x)\r
+ {\r
+ gsl_matrix_add_constant(m_matrix, x);\r
+ return (*this);\r
+ }\r
+ \r
+ Matrix & operator+=(Matrix const& _matrix)\r
+ {\r
+ gsl_matrix_add(m_matrix, _matrix.m_matrix);\r
+ return (*this);\r
+ }\r
+ \r
+ Matrix & operator-=(Matrix const& _matrix)\r
+ {\r
+ gsl_matrix_sub(m_matrix, _matrix.m_matrix);\r
+ return (*this);\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
+ friend \r
+ void swap(Matrix & m1, Matrix & m2);\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
+private:\r
+ size_t m_rows, m_columns;\r
+ gsl_matrix* m_matrix;\r
+};\r
+\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
+\r
+} } // end namespaces\r
+\r
+#endif /*_NL_MATRIX_H_*/\r
diff --git a/src/2geom/numeric/vector.h b/src/2geom/numeric/vector.h
--- /dev/null
@@ -0,0 +1,184 @@
+#ifndef _NL_VECTOR_H_\r
+#define _NL_VECTOR_H_\r
+\r
+#include <cassert>\r
+#include <utility>\r
+\r
+#include <gsl/gsl_vector.h>\r
+\r
+\r
+namespace Geom { namespace NL {\r
+\r
+class Vector;\r
+void swap(Vector & v1, Vector & v2);\r
+\r
+\r
+class Vector\r
+{\r
+public:\r
+ Vector(size_t n)\r
+ : m_size(n)\r
+ {\r
+ m_vector = gsl_vector_alloc(n);\r
+ }\r
+ \r
+ Vector(size_t n, double x)\r
+ : m_size(n)\r
+ {\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
+ : m_size(n)\r
+ {\r
+ m_vector = gsl_vector_alloc(n);\r
+ gsl_vector_set_basis(m_vector, i);\r
+ }\r
+ \r
+ Vector(Vector const& _vector)\r
+ : m_size(_vector.size())\r
+ {\r
+ m_vector = gsl_vector_alloc(size());\r
+ gsl_vector_memcpy(m_vector, _vector.m_vector);\r
+ }\r
+ \r
+ virtual ~Vector()\r
+ {\r
+ gsl_vector_free(m_vector);\r
+ }\r
+ \r
+ Vector & operator=(Vector const& _vector)\r
+ {\r
+ assert( size() == _vector.size() );\r
+ gsl_vector_memcpy(m_vector, _vector.m_vector);\r
+ }\r
+ \r
+ void set_all(double x = 0)\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
+ double const& operator[](size_t i) const\r
+ {\r
+ return *gsl_vector_const_ptr(m_vector, i);\r
+ }\r
+ \r
+ double & operator[](size_t i)\r
+ {\r
+ return *gsl_vector_ptr(m_vector, i);\r
+ }\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
+ Vector & scale(double x)\r
+ {\r
+ gsl_vector_scale(m_vector, x);\r
+ return (*this);\r
+ }\r
+ \r
+ Vector & translate(double x)\r
+ {\r
+ gsl_vector_add_constant(m_vector, x);\r
+ return (*this);\r
+ }\r
+ \r
+ Vector & operator+=(Vector const& _vector)\r
+ {\r
+ gsl_vector_add(m_vector, _vector.m_vector);\r
+ return (*this);\r
+ }\r
+ \r
+ Vector & operator-=(Vector const& _vector)\r
+ {\r
+ gsl_vector_sub(m_vector, _vector.m_vector);\r
+ return (*this);\r
+ }\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
+ friend\r
+ void swap(Vector & v1, Vector & v2);\r
+ \r
+ size_t size() const\r
+ {\r
+ return m_size;\r
+ }\r
+ \r
+private:\r
+ size_t m_size;\r
+ gsl_vector* m_vector;\r
+};\r
+\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
+} } // end namespaces\r
+\r
+\r
+#endif /*_NL_VECTOR_H_*/\r
diff --git a/src/2geom/path.cpp b/src/2geom/path.cpp
index 01c1f2b8ee986ad2b6eb7c424efce574acc440d8..fdfa77c791cf36d9a02d98590c53c6bdb5db14ec 100644 (file)
--- a/src/2geom/path.cpp
+++ b/src/2geom/path.cpp
namespace Geom {
+
int CurveHelpers::root_winding(Curve const &c, Point p) {
std::vector<double> ts = c.roots(p[Y], Y);
return wind;
}
+
+template<>
+inline
+double LineSegment::nearestPoint(Point const& p, double from, double to) const
+{
+ if ( from > to ) std::swap(from, to);
+ Point ip = pointAt(from);
+ Point fp = pointAt(to);
+ Point v = fp - ip;
+ double t = dot( p - ip, v ) / L2sq(v);
+ if ( t <= 0 ) return from;
+ else if ( t >= 1 ) return to;
+ else return from + t*(to-from);
+}
+
+
void Path::swap(Path &other) {
std::swap(curves_, other.curves_);
std::swap(closed_, other.closed_);
return ret;
}
+std::vector<double>
+Path::allNearestPoints(Point const& _point, double from, double to) const
+{
+ if ( from > to ) std::swap(from, to);
+ const Path& _path = *this;
+ unsigned int sz = _path.size();
+ if ( _path.closed() ) ++sz;
+ if ( from < 0 || to > sz )
+ {
+ THROW_RANGEERROR("[from,to] interval out of bounds");
+ }
+ double sif, st = modf(from, &sif);
+ double eif, et = modf(to, &eif);
+ unsigned int si = static_cast<unsigned int>(sif);
+ unsigned int ei = static_cast<unsigned int>(eif);
+ if ( si == sz )
+ {
+ --si;
+ st = 1;
+ }
+ if ( ei == sz )
+ {
+ --ei;
+ et = 1;
+ }
+ if ( si == ei )
+ {
+ std::vector<double> all_nearest =
+ _path[si].allNearestPoints(_point, st, et);
+ for ( unsigned int i = 0; i < all_nearest.size(); ++i )
+ {
+ all_nearest[i] = si + all_nearest[i];
+ }
+ return all_nearest;
+ }
+ std::vector<double> all_t;
+ std::vector< std::vector<double> > all_np;
+ all_np.push_back( _path[si].allNearestPoints(_point, st) );
+ std::vector<unsigned int> ni;
+ ni.push_back(si);
+ double dsq;
+ double mindistsq
+ = distanceSq( _point, _path[si].pointAt( all_np.front().front() ) );
+ Rect bb;
+ for ( unsigned int i = si + 1; i < ei; ++i )
+ {
+ bb = _path[i].boundsFast();
+ dsq = distanceSq(_point, bb);
+ if ( mindistsq < dsq ) continue;
+ all_t = _path[i].allNearestPoints(_point);
+ dsq = distanceSq( _point, _path[i].pointAt( all_t.front() ) );
+ if ( mindistsq > dsq )
+ {
+ all_np.clear();
+ all_np.push_back(all_t);
+ ni.clear();
+ ni.push_back(i);
+ mindistsq = dsq;
+ }
+ else if ( mindistsq == dsq )
+ {
+ all_np.push_back(all_t);
+ ni.push_back(i);
+ }
+ }
+ bb = _path[ei].boundsFast();
+ dsq = distanceSq(_point, bb);
+ if ( mindistsq >= dsq )
+ {
+ all_t = _path[ei].allNearestPoints(_point, 0, et);
+ dsq = distanceSq( _point, _path[ei].pointAt( all_t.front() ) );
+ if ( mindistsq > dsq )
+ {
+ for ( unsigned int i = 0; i < all_t.size(); ++i )
+ {
+ all_t[i] = ei + all_t[i];
+ }
+ return all_t;
+ }
+ else if ( mindistsq == dsq )
+ {
+ all_np.push_back(all_t);
+ ni.push_back(ei);
+ }
+ }
+ std::vector<double> all_nearest;
+ for ( unsigned int i = 0; i < all_np.size(); ++i )
+ {
+ for ( unsigned int j = 0; j < all_np[i].size(); ++j )
+ {
+ all_nearest.push_back( ni[i] + all_np[i][j] );
+ }
+ }
+ return all_nearest;
+}
+
+double Path::nearestPoint(Point const& _point, double from, double to) const
+{
+ if ( from > to ) std::swap(from, to);
+ const Path& _path = *this;
+ unsigned int sz = _path.size();
+ if ( _path.closed() ) ++sz;
+ if ( from < 0 || to > sz )
+ {
+ THROW_RANGEERROR("[from,to] interval out of bounds");
+ }
+ double sif, st = modf(from, &sif);
+ double eif, et = modf(to, &eif);
+ unsigned int si = static_cast<unsigned int>(sif);
+ unsigned int ei = static_cast<unsigned int>(eif);
+ if ( si == sz )
+ {
+ --si;
+ st = 1;
+ }
+ if ( ei == sz )
+ {
+ --ei;
+ et = 1;
+ }
+ if ( si == ei )
+ {
+ double nearest =
+ _path[si].nearestPoint(_point, st, et);
+ return si + nearest;
+ }
+ double t;
+ double nearest = _path[si].nearestPoint(_point, st);
+ unsigned int ni = si;
+ double dsq;
+ double mindistsq = distanceSq(_point, _path[si].pointAt(nearest));
+ Rect bb;
+ for ( unsigned int i = si + 1; i < ei; ++i )
+ {
+ bb = _path[i].boundsFast();
+ dsq = distanceSq(_point, bb);
+ if ( mindistsq <= dsq ) continue;
+ t = _path[i].nearestPoint(_point);
+ dsq = distanceSq(_point, _path[i].pointAt(t));
+ if ( mindistsq > dsq )
+ {
+ nearest = t;
+ ni = i;
+ mindistsq = dsq;
+ }
+ }
+ bb = _path[ei].boundsFast();
+ dsq = distanceSq(_point, bb);
+ if ( mindistsq > dsq )
+ {
+ t = _path[ei].nearestPoint(_point, 0, et);
+ dsq = distanceSq(_point, _path[ei].pointAt(t));
+ if ( mindistsq > dsq )
+ {
+ nearest = t;
+ ni = ei;
+ }
+ }
+ return ni + nearest;
+}
+
//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);
void Path::append(Curve const &curve) {
if ( curves_.front() != final_ && !are_near(curve.initialPoint(), (*final_)[0], eps) ) {
- throwContinuityError();
+ THROW_CONTINUITYERROR();
}
do_append(curve.duplicate());
}
if ( curves_.front() != final_ ) {
for ( int i = 0 ; i < 2 ; ++i ) {
if ( !are_near(curve[i][0][0], (*final_)[0][i], eps) ) {
- throwContinuityError();
+ THROW_CONTINUITYERROR();
}
}
}
if ( first != last ) {
if ( first_replaced != curves_.begin() ) {
if ( !are_near( (*first_replaced)->initialPoint(), (*first)->initialPoint(), eps ) ) {
- throwContinuityError();
+ THROW_CONTINUITYERROR();
}
}
if ( last_replaced != (curves_.end()-1) ) {
if ( !are_near( (*(last_replaced-1))->finalPoint(), (*(last-1))->finalPoint(), eps ) ) {
- throwContinuityError();
+ 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 ) ) {
- throwContinuityError();
+ THROW_CONTINUITYERROR();
}
}
}
diff --git a/src/2geom/path.h b/src/2geom/path.h
index 414d69755870dbb6b6eba09dd468ceb2a187c99b..ed15260fc40ba488f95fcf311aeb53120f69e807 100644 (file)
--- a/src/2geom/path.h
+++ b/src/2geom/path.h
#include "bezier.h"
#include "crossing.h"
#include "utils.h"
+#include "nearest-point.h"
namespace Geom {
virtual void setInitial(Point v) = 0;
virtual void setFinal(Point v) = 0;
+
+ virtual
+ double nearestPoint( Point const& p, double from = 0, double to = 1 ) const
+ {
+ return nearest_point(p, toSBasis(), from, to);
+ }
+
+ virtual
+ std::vector<double>
+ allNearestPoints( Point const& p, double from = 0, double to = 1 ) const
+ {
+ return all_nearest_points(p, toSBasis(), from, to);
+ }
virtual Curve *transformed(Matrix const &m) const = 0;
};
class SBasisCurve : public Curve {
+
private:
SBasisCurve();
D2<SBasis> inner;
+
public:
explicit SBasisCurve(D2<SBasis> const &sb) : inner(sb) {}
explicit SBasisCurve(Curve const &other) : inner(other.toSBasis()) {}
Rect boundsLocal(Interval i, unsigned deg) const { return bounds_local(inner, i, deg); }
std::vector<double> roots(double v, Dim2 d) const { return Geom::roots(inner[d] - v); }
+
+ double nearestPoint( Point const& p, double from = 0, double to = 1 ) const
+ {
+ return nearest_point(p, inner, from, to);
+ }
+
+ std::vector<double>
+ allNearestPoints( Point const& p, double from = 0, double to = 1 ) const
+ {
+ return all_nearest_points(p, inner, from, to);
+ }
Curve *portion(double f, double t) const {
return new SBasisCurve(Geom::portion(inner, f, t));
template <unsigned order>
class BezierCurve : public Curve {
+
private:
D2<Bezier > inner;
+
public:
template <unsigned required_degree>
static void assert_degree(BezierCurve<required_degree> const *) {}
roots(double v, Dim2 d) const {
return (inner[d] - v).roots();
}
-
+
+ double nearestPoint( Point const& p, double from = 0, double to = 1 ) const
+ {
+ return Curve::nearestPoint(p, from, to);
+ }
+
void setPoints(std::vector<Point> ps) {
for(unsigned i = 0; i <= order; i++) {
setPoint(i, ps[i]);
typedef BezierCurve<2> QuadraticBezier;
typedef BezierCurve<3> CubicBezier;
+template<>
+double LineSegment::nearestPoint(Point const& p, double from, double to) const;
{
public:
SVGEllipticalArc()
- {}
+ : m_initial_point(Point(0,0)), m_final_point(Point(0,0)),
+ m_rx(0), m_ry(0), m_rot_angle(0),
+ m_large_arc(true), m_sweep(true)
+ {
+ m_start_angle = m_end_angle = 0;
+ m_center = Point(0,0);
+ }
SVGEllipticalArc( Point _initial_point, double _rx, double _ry,
double _rot_angle, bool _large_arc, bool _sweep,
m_rx(_rx), m_ry(_ry), m_rot_angle(_rot_angle),
m_large_arc(_large_arc), m_sweep(_sweep)
{
- //assert( (ray(X) >= 0) && (ray(Y) >= 0) );
- if ( are_near(initialPoint(), finalPoint()) )
- {
- m_start_angle = m_end_angle = 0;
- m_center = initialPoint();
- }
- else
- {
calculate_center_and_extreme_angles();
- }
}
+ void set( Point _initial_point, double _rx, double _ry,
+ double _rot_angle, bool _large_arc, bool _sweep,
+ Point _final_point
+ )
+ {
+ m_initial_point = _initial_point;
+ m_final_point = _final_point;
+ m_rx = _rx;
+ m_ry = _ry;
+ m_rot_angle = _rot_angle;
+ m_large_arc = _large_arc;
+ m_sweep = _sweep;
+ calculate_center_and_extreme_angles();
+ }
+
Curve* duplicate() const
{
return new SVGEllipticalArc(*this);
}
- double center(Dim2 i) const
+ double center(unsigned int i) const
{
return m_center[i];
}
return m_end_angle;
}
- double ray(Dim2 i) const
+ double ray(unsigned int i) const
{
return (i == 0) ? m_rx : m_ry;
}
bool isDegenerate() const
{
- return are_near(initialPoint(), finalPoint());
+ return ( are_near(ray(X), 0) || are_near(ray(Y), 0) );
}
// TODO: native implementation of the following methods
Rect boundsFast() const
{
- return SBasisCurve(toSBasis()).boundsFast();
+ return boundsExact();
}
- Rect boundsExact() const
- {
- return SBasisCurve(toSBasis()).boundsExact();
- }
+ Rect boundsExact() const;
Rect boundsLocal(Interval i, unsigned int deg) const
{
return SBasisCurve(toSBasis()).boundsLocal(i, deg);
}
- std::vector<double> roots(double v, Dim2 d) const
+ std::vector<double> roots(double v, Dim2 d) const;
+
+ std::vector<double>
+ allNearestPoints( Point const& p, double from = 0, double to = 1 ) const;
+
+ double nearestPoint( Point const& p, double from = 0, double to = 1 ) const
{
- return SBasisCurve(toSBasis()).roots(v, d);
+ if ( are_near(ray(X), ray(Y)) && are_near(center(), p) )
+ {
+ return from;
+ }
+ return allNearestPoints(p, from, to).front();
}
int winding(Point p) const
return SBasisCurve(toSBasis()).winding(p);
}
- Curve *derivative() const
- {
- return SBasisCurve(toSBasis()).derivative();
- }
+ Curve *derivative() const;
Curve *transformed(Matrix const &m) const
{
return SBasisCurve(toSBasis()).transformed(m);
}
- std::vector<Point> pointAndDerivatives(Coord t, unsigned n) const
- {
- return SBasisCurve(toSBasis()).pointAndDerivatives(t, n);
- }
+ std::vector<Point> pointAndDerivatives(Coord t, unsigned int n) const;
D2<SBasis> toSBasis() const;
- double valueAt(Coord t, Dim2 d) const;
-
- Point pointAt(Coord t) const
+ bool containsAngle(Coord angle) const;
+
+ double valueAtAngle(Coord t, Dim2 d) const;
+
+ Point pointAtAngle(Coord t) const
{
- Coord tt = from_01_to_02PI(t);
double sin_rot_angle = std::sin(rotation_angle());
double cos_rot_angle = std::cos(rotation_angle());
Matrix m( ray(X) * cos_rot_angle, ray(X) * sin_rot_angle,
-ray(Y) * sin_rot_angle, ray(Y) * cos_rot_angle,
center(X), center(Y) );
- Point p( std::cos(tt), std::sin(tt) );
+ Point p( std::cos(t), std::sin(t) );
return p * m;
}
+
+ double valueAt(Coord t, Dim2 d) const
+ {
+ Coord tt = map_to_02PI(t);
+ return valueAtAngle(tt, d);
+ }
+
+ Point pointAt(Coord t) const
+ {
+ Coord tt = map_to_02PI(t);
+ return pointAtAngle(tt);
+ }
std::pair<SVGEllipticalArc, SVGEllipticalArc>
subdivide(Coord t) const
return rarc;
}
- private:
double sweep_angle() const
{
Coord d = end_angle() - start_angle();
d += 2*M_PI;
return d;
}
-
- Coord from_01_to_02PI(Coord t) const;
-
- void calculate_center_and_extreme_angles() throw(RangeError);
+ private:
+ Coord map_to_02PI(Coord t) const;
+ Coord map_to_01(Coord angle) const;
+ void calculate_center_and_extreme_angles();
private:
Point m_initial_point, m_final_point;
double m_rx, m_ry, m_rot_angle;
return ret;
}
- Point pointAt(double t) const {
- if(empty()) return Point(0,0);
- double i, f = modf(t, &i);
- if(i == size() && f == 0) { i--; }
- assert(i >= 0 && i <= size());
- return (*this)[unsigned(i)].pointAt(f);
+ Point pointAt(double t) const
+ {
+ unsigned int sz = size();
+ if ( closed() ) ++sz;
+ if ( t < 0 || t > sz )
+ {
+ THROW_RANGEERROR("parameter t out of bounds");
+ }
+ if ( empty() ) return Point(0,0);
+ double k, lt = modf(t, &k);
+ unsigned int i = static_cast<unsigned int>(k);
+ if ( i == sz )
+ {
+ --i;
+ lt = 1;
+ }
+ return (*this)[i].pointAt(lt);
}
- double valueAt(double t, Dim2 d) const {
- if(empty()) return 0;
- double i, f = modf(t, &i);
- if(i == size() && f == 0) { i--; }
- assert(i >= 0 && i <= size());
- return (*this)[unsigned(i)].valueAt(f, d);
+ double valueAt(double t, Dim2 d) const
+ {
+ unsigned int sz = size();
+ if ( closed() ) ++sz;
+ if ( t < 0 || t > sz )
+ {
+ THROW_RANGEERROR("parameter t out of bounds");
+ }
+ if ( empty() ) return 0;
+ double k, lt = modf(t, &k);
+ unsigned int i = static_cast<unsigned int>(k);
+ if ( i == sz )
+ {
+ --i;
+ lt = 1;
+ }
+ return (*this)[i].valueAt(lt, d);
}
std::vector<double> roots(double v, Dim2 d) const {
}
return res;
}
-
+
+ std::vector<double>
+ allNearestPoints(Point const& _point, double from, double to) const;
+
+ std::vector<double>
+ allNearestPoints(Point const& _point) const
+ {
+ unsigned int sz = size();
+ if ( closed() ) ++sz;
+ return allNearestPoints(_point, 0, sz);
+ }
+
+
+ double nearestPoint(Point const& _point, double from, double to) const;
+
+ double nearestPoint(Point const& _point) const
+ {
+ unsigned int sz = size();
+ if ( closed() ) ++sz;
+ return nearestPoint(_point, 0, sz);
+ }
+
void appendPortionTo(Path &p, double f, double t) const;
Path portion(double f, double t) const {
}
return ret;
}
-
+
void insert(iterator pos, Curve const &curve) {
Sequence source(1, curve.duplicate());
try {
Rect boundsFast() const { return actualPath().boundsFast; }
Rect boundsExact() const { return actualPath().boundsFast; }
- Rect boundsLocal(Interval i) const { throwNotImplemented(); }
+ Rect boundsLocal(Interval i) const { THROW_NOTIMPLEMENTED(); }
std::vector<double> roots(double v, Dim2 d) const = 0;
diff --git a/src/2geom/piecewise.h b/src/2geom/piecewise.h
index 1d0fd93baf17a9576f9b3b8b4481e974f2e44cba..d548ea0e9f76f98fbd1d37d472663b213528336c 100644 (file)
--- a/src/2geom/piecewise.h
+++ b/src/2geom/piecewise.h
}
//Convenience/implementation hiding function to add cuts.
inline void push_cut(double c) {
- assert_invariants(cuts.empty() || c > cuts.back());
+ ASSERT_INVARIANTS(cuts.empty() || c > cuts.back());
cuts.push_back(c);
}
//Convenience/implementation hiding function to add segments.
//Offsets the piecewise domain
inline void offsetDomain(double o) {
- assert(is_finite(o));
+ assert(IS_FINITE(o));
if(o != 0)
for(unsigned i = 0; i <= size(); i++)
cuts[i] += o;
diff --git a/src/2geom/point.cpp b/src/2geom/point.cpp
index 1e2a3463f006b3627829d74400478a629cd81130..66eaf8ca671f1a06e49ce518be323167beecea46 100644 (file)
--- a/src/2geom/point.cpp
+++ b/src/2geom/point.cpp
void Point::normalize() {
double len = hypot(_pt[0], _pt[1]);
if(len == 0) return;
- if(is_nan(len)) return;
+ if(IS_NAN(len)) return;
static double const inf = 1e400;
if(len != inf) {
*this /= len;
Coord LInfty(Point const &p) {
Coord const a(fabs(p[0]));
Coord const b(fabs(p[1]));
- return ( a < b || is_nan(b)
+ return ( a < b || IS_NAN(b)
? b
: a );
}
diff --git a/src/2geom/rect.h b/src/2geom/rect.h
index 86f9ec1406892eaccb92a9348bde3d58c5a304f6..c89946606fe86886924c9f54ae778882430dec1e 100644 (file)
--- a/src/2geom/rect.h
+++ b/src/2geom/rect.h
return x && y ? boost::optional<Rect>(Rect(*x, *y)) : boost::optional<Rect>();
}
+inline
+double distanceSq( Point const& p, Rect const& rect )
+{
+ double dx = 0, dy = 0;
+ if ( p[X] < rect.left() )
+ {
+ dx = p[X] - rect.left();
+ }
+ else if ( p[X] > rect.right() )
+ {
+ dx = rect.right() - p[X];
+ }
+ if ( p[Y] < rect.top() )
+ {
+ dy = rect.top() - p[Y];
+ }
+ else if ( p[Y] > rect.bottom() )
+ {
+ dy = p[Y] - rect.bottom();
+ }
+ return dx*dx + dy*dy;
}
+inline
+double distance( Point const& p, Rect const& rect )
+{
+ return std::sqrt(distanceSq(p, rect));
+}
+
+
+} // end namespace Geom
+
#endif //_2GEOM_RECT
#endif //_2GEOM_D2
/*
index 497091d1c65f29e7b95b4d4dc9d4e7ffd61ddd6d..9843b78d73c9af44efa1f6641f7a2342eda79370 100644 (file)
// mutating
void
subpath_from_sbasis(Geom::OldPathSetBuilder &pb, D2<SBasis> const &B, double tol, bool initial) {
- assert(B.is_finite());
+ assert(B.IS_FINITE());
if(B.tail_error(2) < tol || B.size() == 2) { // nearly cubic enough
if(B.size() == 1) {
if (initial) {
subpath_from_sbasis_incremental(Geom::OldPathSetBuilder &pb, D2<SBasis> B, double tol, bool initial) {
const unsigned k = 2; // cubic bezier
double te = B.tail_error(k);
- assert(B[0].is_finite());
- assert(B[1].is_finite());
+ assert(B[0].IS_FINITE());
+ assert(B[1].IS_FINITE());
//std::cout << "tol = " << tol << std::endl;
while(1) {
@@ -181,7 +181,7 @@ subpath_from_sbasis_incremental(Geom::OldPathSetBuilder &pb, D2<SBasis> B, doubl
void build_from_sbasis(Geom::PathBuilder &pb, D2<SBasis> const &B, double tol) {
if (!B.isFinite()) {
- throwException("assertion failed: B.isFinite()");
+ THROW_EXCEPTION("assertion failed: B.isFinite()");
}
if(tail_error(B, 2) < tol || sbasis_size(B) == 2) { // nearly cubic enough
if(sbasis_size(B) <= 1) {
diff --git a/src/2geom/sbasis.cpp b/src/2geom/sbasis.cpp
index 7157bc8856313021f405feee0f10371719ff86d2..c6013204374d19a2b69bb506d4e535fab6aeaf28 100644 (file)
--- a/src/2geom/sbasis.cpp
+++ b/src/2geom/sbasis.cpp
return true;
}
+std::vector<double> SBasis::valueAndDerivatives(double t, unsigned n) const {
+ std::vector<double> ret;
+ if(n==1) {
+ ret.push_back(valueAt(t));
+ return ret;
+ }
+ if(n==2) {
+ double der;
+ ret.push_back(valueAndDerivative(t, der));
+ ret.push_back(der);
+ return ret;
+ }
+ SBasis tmp = *this;
+ while(n > 0) {
+ ret.push_back(tmp.valueAt(t));
+ tmp = derivative(tmp);
+ n--;
+ }
+ return ret;
+}
+
+
SBasis operator+(const SBasis& a, const SBasis& b) {
SBasis result;
const unsigned out_size = std::max(a.size(), b.size());
diff --git a/src/2geom/sbasis.h b/src/2geom/sbasis.h
index 82ff989e7b3d01e9b4501926bd02841dcf53bbc9..7f0d88f8cd2211e7443e3bcc6f71ec935a07d9fa 100644 (file)
--- a/src/2geom/sbasis.h
+++ b/src/2geom/sbasis.h
return valueAt(t);
}
- std::vector<double> valueAndDerivatives(double t, unsigned n) const {
- std::vector<double> ret;
- if(n==1) {
- ret.push_back(valueAt(t));
- return ret;
- }
- if(n==2) {
- double der;
- ret.push_back(valueAndDerivative(t, der));
- ret.push_back(der);
- return ret;
- }
- //TODO
- throwNotImplemented();
- }
+ std::vector<double> valueAndDerivatives(double t, unsigned n) const;
SBasis toSBasis() const { return SBasis(*this); }
diff --git a/src/2geom/shape.cpp b/src/2geom/shape.cpp
index 82ae411451e1fd993b4d043281d19c147508a444..a7a40066e9b3a514bfcba85715a48a1acc9b3b56 100644 (file)
--- a/src/2geom/shape.cpp
+++ b/src/2geom/shape.cpp
@@ -191,7 +191,7 @@ Shape shape_boolean_rb(bool rev, Shape const &a, Shape const &b, CrossingSet con
* NOTE: currently doesn't work, as the CrossingSet reversal functions crash
*/
Shape boolop(Shape const &a, Shape const &b, unsigned flags, CrossingSet const &crs) {
- throwNotImplemented();
+ THROW_NOTIMPLEMENTED();
flags &= 15;
if(flags <= BOOLOP_UNION) {
switch(flags) {
index 6583c948f2b16a10b4e53f4eea99717dcffe0240..b3420fba6c9426ad8f6adff23653b1ec7a01a38f 100644 (file)
-/*
- * SVG Elliptical Arc Class
- *
- * Copyright 2008 Marco Cecchetti <mrcekets at gmail.com>
- *
- * This library is free software; you can redistribute it and/or
- * modify it either under the terms of the GNU Lesser General Public
- * License version 2.1 as published by the Free Software Foundation
- * (the "LGPL") or, at your option, under the terms of the Mozilla
- * Public License Version 1.1 (the "MPL"). If you do not alter this
- * notice, a recipient may use your version of this file under either
- * the MPL or the LGPL.
- *
- * You should have received a copy of the LGPL along with this library
- * in the file COPYING-LGPL-2.1; if not, write to the Free Software
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
- * You should have received a copy of the MPL along with this library
- * in the file COPYING-MPL-1.1
- *
- * The contents of this file are subject to the Mozilla Public License
- * Version 1.1 (the "License"); you may not use this file except in
- * compliance with the License. You may obtain a copy of the License at
- * http://www.mozilla.org/MPL/
- *
- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
- * OF ANY KIND, either express or implied. See the LGPL or the MPL for
- * the specific language governing rights and limitations.
- */
-
-
-#include "path.h"
-
-
-namespace Geom
-{
-
-D2<SBasis> SVGEllipticalArc::toSBasis() const
-{
- // the interval of parametrization has to be [0,1]
- Coord et = start_angle() + ( sweep_flag() ? sweep_angle() : -sweep_angle() );
- Linear param(start_angle(), et);
- Coord cos_rot_angle = std::cos(rotation_angle());
- Coord sin_rot_angle = std::sin(rotation_angle());
- // order = 4 seems to be enough to get a perfect looking elliptical arc
- // should it be choosen in function of the arc length anyway ?
- // or maybe a user settable parameter: toSBasis(unsigned int order) ?
- SBasis arc_x = ray(X) * cos(param,4);
- SBasis arc_y = ray(Y) * sin(param,4);
- D2<SBasis> arc;
- arc[0] = arc_x * cos_rot_angle - arc_y * sin_rot_angle + Linear(center(X),center(X));
- arc[1] = arc_x * sin_rot_angle + arc_y * cos_rot_angle + Linear(center(Y),center(Y));
- return arc;
-}
-
-double SVGEllipticalArc::valueAt(Coord t, Dim2 d) const
-{
- Coord tt = from_01_to_02PI(t);
- double sin_rot_angle = std::sin(rotation_angle());
- double cos_rot_angle = std::cos(rotation_angle());
- if ( d == X )
- {
- return ray(X) * cos_rot_angle * std::cos(tt)
- - ray(Y) * sin_rot_angle * std::sin(tt)
- + center(X);
- }
- else
- {
- return ray(X) * sin_rot_angle * std::cos(tt)
- + ray(Y) * cos_rot_angle * std::sin(tt)
- + center(Y);
- }
-}
-
-
-Curve* SVGEllipticalArc::portion(double f, double t) const
-{
- if (f < 0) f = 0;
- if (f > 1) f = 1;
- if (t < 0) t = 0;
- if (t > 1) t = 1;
- SVGEllipticalArc* arc = new SVGEllipticalArc( *this );
- arc->m_initial_point = pointAt(f);
- arc->m_final_point = pointAt(t);
- double sa = sweep_flag() ? sweep_angle() : -sweep_angle();
- arc->m_start_angle = m_start_angle + sa * f;
- if ( !(arc->m_start_angle < 2*M_PI) )
- arc->m_start_angle -= 2*M_PI;
- if ( !(arc->m_start_angle > 0) )
- arc->m_start_angle += 2*M_PI;
- arc->m_end_angle = m_start_angle + sa * t;
- if ( !(arc->m_end_angle < 2*M_PI) )
- arc->m_end_angle -= 2*M_PI;
- if ( !(arc->m_end_angle > 0) )
- arc->m_end_angle += 2*M_PI;
- if ( f > t ) arc->m_sweep = !sweep_flag();
- if ( large_arc_flag() && (arc->sweep_angle() < M_PI) )
- arc->m_large_arc = false;
- return arc;
-}
-
-// NOTE: doesn't work with 360 deg arcs
-void SVGEllipticalArc::calculate_center_and_extreme_angles() throw(RangeError)
-{
- double sin_rot_angle = std::sin(rotation_angle());
- double cos_rot_angle = std::cos(rotation_angle());
-
- Point sp = sweep_flag() ? initialPoint() : finalPoint();
- Point ep = sweep_flag() ? finalPoint() : initialPoint();
-
- Matrix m( ray(X) * cos_rot_angle, ray(X) * sin_rot_angle,
- -ray(Y) * sin_rot_angle, ray(Y) * cos_rot_angle,
- 0, 0 );
- Matrix im = m.inverse();
- Point sol = (ep - sp) * im;
- double half_sum_angle = std::atan2(-sol[X], sol[Y]);
- double half_diff_angle;
- if ( are_near(std::fabs(half_sum_angle), M_PI/2) )
- {
- double anti_sgn_hsa = (half_sum_angle > 0) ? -1 : 1;
- double arg = anti_sgn_hsa * sol[X] / 2;
- // if |arg| is a little bit > 1 acos returns nan
- if ( are_near(arg, 1) )
- half_diff_angle = 0;
- else if ( are_near(arg, -1) )
- half_diff_angle = M_PI;
- else
- {
- if ( !(-1 < arg && arg < 1) )
- throwRangeError("there is no ellipse that satisfies the given "
- "constraints");
- // assert( -1 < arg && arg < 1 );
- // if it fails
- // => there is no ellipse that satisfies the given constraints
- half_diff_angle = std::acos( arg );
- }
-
- half_diff_angle = M_PI/2 - half_diff_angle;
- }
- else
- {
- double arg = sol[Y] / ( 2 * std::cos(half_sum_angle) );
- // if |arg| is a little bit > 1 asin returns nan
- if ( are_near(arg, 1) )
- half_diff_angle = M_PI/2;
- else if ( are_near(arg, -1) )
- half_diff_angle = -M_PI/2;
- else
- {
- if ( !(-1 < arg && arg < 1) )
- throwRangeError("there is no ellipse that satisfies the given "
- "constraints");
- // assert( -1 < arg && arg < 1 );
- // if it fails
- // => there is no ellipse that satisfies the given constraints
- half_diff_angle = std::asin( arg );
- }
- }
-
- if ( ( m_large_arc && half_diff_angle > 0 )
- || (!m_large_arc && half_diff_angle < 0 ) )
- {
- half_diff_angle = -half_diff_angle;
- }
- if ( half_sum_angle < 0 ) half_sum_angle += 2*M_PI;
- if ( half_diff_angle < 0 ) half_diff_angle += M_PI;
-
- m_start_angle = half_sum_angle - half_diff_angle;
- m_end_angle = half_sum_angle + half_diff_angle;
- // 0 <= m_start_angle, m_end_angle < 2PI
- if ( m_start_angle < 0 ) m_start_angle += 2*M_PI;
- if( !(m_end_angle < 2*M_PI) ) m_end_angle -= 2*M_PI;
- sol[0] = std::cos(m_start_angle);
- sol[1] = std::sin(m_start_angle);
- m_center = sp - sol * m;
- if ( !sweep_flag() )
- {
- double angle = m_start_angle;
- m_start_angle = m_end_angle;
- m_end_angle = angle;
- }
-}
-
-Coord SVGEllipticalArc::from_01_to_02PI(Coord t) const
-{
- if ( sweep_flag() )
- {
- Coord angle = start_angle() + sweep_angle() * t;
- if ( !(angle < 2*M_PI) )
- angle -= 2*M_PI;
- return angle;
- }
- else
- {
- Coord angle = start_angle() - sweep_angle() * t;
- if ( angle < 0 ) angle += 2*M_PI;
- return angle;
- }
-}
-
-
-} // end namespace Geom
-
-
-/*
- Local Variables:
- mode:c++
- c-file-style:"stroustrup"
- c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
- indent-tabs-mode:nil
- fill-column:99
- End:
-*/
-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
-
-
+/*\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 "path.h"\r
+#include "angle.h"\r
+\r
+#include <gsl/gsl_poly.h>\r
+\r
+#include <cfloat>\r
+\r
+\r
+\r
+\r
+namespace Geom\r
+{\r
+\r
+\r
+Rect SVGEllipticalArc::boundsExact() const\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
+\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
+ 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
+ \r
+ if ( initialPoint()[d] == v && finalPoint()[d] == v )\r
+ {\r
+ THROW_EXCEPTION("infinite solutions");\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
+ 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
+// return SBasisCurve(toSBasis()).roots(v, d);\r
+}\r
+\r
+// D(E(t,C),t) = E(t+PI/2,O)\r
+Curve* SVGEllipticalArc::derivative() const\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
+ std::vector<Point> result;\r
+ result.reserve(n);\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(n, 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 = n / 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 = n - 4 * m;\r
+ for ( unsigned int i = 0; i < m; ++i )\r
+ {\r
+ result.push_back( result[i] );\r
+ }\r
+ if ( !result.empty() ) // n != 0\r
+ result[0] = pointAtAngle(angle);\r
+ return result;\r
+}\r
+\r
+D2<SBasis> SVGEllipticalArc::toSBasis() const\r
+{\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
+ D2<SBasis> arc;\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
+bool SVGEllipticalArc::containsAngle(Coord angle) const\r
+{\r
+ if ( sweep_flag() )\r
+ if ( start_angle() < end_angle() )\r
+ if ( !( angle < start_angle() || angle > end_angle() ) )\r
+ return true;\r
+ else\r
+ return false;\r
+ else\r
+ if ( !( angle < start_angle() && angle > end_angle() ) )\r
+ return true;\r
+ else\r
+ return false;\r
+ else\r
+ if ( start_angle() > end_angle() )\r
+ if ( !( angle > start_angle() || angle < end_angle() ) )\r
+ return true;\r
+ else\r
+ return false;\r
+ else\r
+ if ( !( angle > start_angle() && angle < end_angle() ) )\r
+ return true;\r
+ else\r
+ return false;\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
+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
+// NOTE: doesn't work with 360 deg arcs\r
+void SVGEllipticalArc::calculate_center_and_extreme_angles()\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
+ double sin_rot_angle = std::sin(rotation_angle());\r
+ double cos_rot_angle = std::cos(rotation_angle());\r
+\r
+ Point sp = sweep_flag() ? initialPoint() : finalPoint();\r
+ Point ep = sweep_flag() ? finalPoint() : initialPoint();\r
+\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
+ 0, 0 );\r
+ Matrix im = m.inverse();\r
+ Point sol = (ep - sp) * im;\r
+ double half_sum_angle = std::atan2(-sol[X], sol[Y]);\r
+ double half_diff_angle;\r
+ if ( are_near(std::fabs(half_sum_angle), M_PI/2) )\r
+ {\r
+ double anti_sgn_hsa = (half_sum_angle > 0) ? -1 : 1;\r
+ double arg = anti_sgn_hsa * sol[X] / 2;\r
+ // if |arg| is a little bit > 1 acos returns nan\r
+ if ( are_near(arg, 1) )\r
+ half_diff_angle = 0;\r
+ else if ( are_near(arg, -1) )\r
+ half_diff_angle = M_PI;\r
+ else\r
+ {\r
+ if ( !(-1 < arg && arg < 1) )\r
+ THROW_RANGEERROR(\r
+ "there is no ellipse that satisfies the given constraints"\r
+ );\r
+ // assert( -1 < arg && arg < 1 );\r
+ // if it fails \r
+ // => there is no ellipse that satisfies the given constraints\r
+ half_diff_angle = std::acos( arg );\r
+ }\r
+\r
+ half_diff_angle = M_PI/2 - half_diff_angle;\r
+ }\r
+ else\r
+ {\r
+ double arg = sol[Y] / ( 2 * std::cos(half_sum_angle) );\r
+ // if |arg| is a little bit > 1 asin returns nan\r
+ if ( are_near(arg, 1) ) \r
+ half_diff_angle = M_PI/2;\r
+ else if ( are_near(arg, -1) )\r
+ half_diff_angle = -M_PI/2;\r
+ else\r
+ {\r
+ if ( !(-1 < arg && arg < 1) )\r
+ THROW_RANGEERROR(\r
+ "there is no ellipse that satisfies the given constraints"\r
+ );\r
+ // assert( -1 < arg && arg < 1 );\r
+ // if it fails \r
+ // => there is no ellipse that satisfies the given constraints\r
+ half_diff_angle = std::asin( arg );\r
+ }\r
+ }\r
+\r
+ if ( ( m_large_arc && half_diff_angle > 0 ) \r
+ || (!m_large_arc && half_diff_angle < 0 ) )\r
+ {\r
+ half_diff_angle = -half_diff_angle;\r
+ }\r
+ if ( half_sum_angle < 0 ) half_sum_angle += 2*M_PI;\r
+ if ( half_diff_angle < 0 ) half_diff_angle += M_PI;\r
+ \r
+ m_start_angle = half_sum_angle - half_diff_angle;\r
+ m_end_angle = half_sum_angle + half_diff_angle;\r
+ // 0 <= m_start_angle, m_end_angle < 2PI\r
+ if ( m_start_angle < 0 ) m_start_angle += 2*M_PI;\r
+ if( !(m_end_angle < 2*M_PI) ) m_end_angle -= 2*M_PI;\r
+ sol[0] = std::cos(m_start_angle);\r
+ sol[1] = std::sin(m_start_angle);\r
+ m_center = sp - sol * m;\r
+ if ( !sweep_flag() )\r
+ {\r
+ double angle = m_start_angle;\r
+ m_start_angle = m_end_angle;\r
+ m_end_angle = angle;\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
+std::vector<double> SVGEllipticalArc::\r
+allNearestPoints( Point const& p, double from, double to ) const\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
+ std::vector<double> result;\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_EXCEPTION("infinite nearest points");\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
+ double coeff[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
+ double sol[8];\r
+ gsl_poly_complex_workspace * w = gsl_poly_complex_workspace_alloc(5);\r
+ gsl_poly_complex_solve(coeff, 5, w, sol );\r
+ gsl_poly_complex_workspace_free(w);\r
+ \r
+ for ( unsigned int i = 0; i < 4; ++i )\r
+ {\r
+ if ( sol[2*i+1] == 0 ) real_sol.push_back(sol[2*i]);\r
+ }\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
+} // 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/svg-path.cpp b/src/2geom/svg-path.cpp
index 60ae957cbeb2d09a32216c5ea96d1f838bb1e65d..c09e5c929e408874cf933d7a8b268653dbfb98d8 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"
namespace Geom {
void output(SVGEllipticalArc const &/*curve*/, SVGPathSink &/*sink*/) {
// FIXME
+ THROW_NOTIMPLEMENTED();
}
template <typename T>
index bdb6bcad3271bdc19fd59c9a6ce62c043c5b90ce..0a27c31f64f30d80e0e882dcb00e13b37a682acc 100644 (file)
--- a/src/2geom/transforms.cpp
+++ b/src/2geom/transforms.cpp
return ret;
}
+Translate pow(Translate const &t, int n) {
+ return Translate(t[0]*n, t[1]*n);
+}
+
+Coord pow(Coord x, long n) // shamelessly lifted from WP
+{
+ Coord result = 1;
+ while ( n ) {
+ if ( n & 1 ) {
+ result = result * x;
+ n = n-1;
+ }
+ x = x*x;
+ n = n/2;
+ }
+ return result;
+}
+Scale pow(Scale const &s, int n) {
+ return Scale(pow(s[0],n), pow(s[1],n));
+
+}
+
+Rotate pow(Rotate x, long n)
+{
+ Rotate result(0,1); // identity
+ while ( n ) {
+ if ( n & 1 ) {
+ result = result * x;
+ n = n-1;
+ }
+ x = x*x;
+ n = n/2;
+ }
+ return result;
+}
+
+Matrix pow(Matrix x, long n)
+{
+ Matrix result;
+ result.setIdentity();
+ while ( n ) {
+ if ( n & 1 ) {
+ result = result * x;
+ n = n-1;
+ }
+ x = x*x;
+ n = n/2;
+ }
+ return result;
+}
+
}
/*
diff --git a/src/2geom/transforms.h b/src/2geom/transforms.h
index 94058bff00f0d067f4035a290bcb23e78185a941..9529d8922f74e14bf12641b448e31934245abeeb 100644 (file)
--- a/src/2geom/transforms.h
+++ b/src/2geom/transforms.h
Matrix operator*(Matrix const &m, Rotate const &r);
Matrix operator*(Matrix const &m1, Matrix const &m2);
+Translate pow(Translate const &t, int n);
+Scale pow(Scale const &t, int n);
+Rotate pow(Rotate t, int n);
+Matrix pow(Matrix t, int n);
+
//TODO: matrix to trans/scale/rotate
} /* namespace Geom */