summaryrefslogtreecommitdiffstats
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/2geom/Makefile_insert5
-rw-r--r--src/2geom/angle.h43
-rw-r--r--src/2geom/basic-intersection.cpp2
-rw-r--r--src/2geom/bezier-utils.cpp14
-rw-r--r--src/2geom/bezier.h2
-rw-r--r--src/2geom/exception.h18
-rw-r--r--src/2geom/isnan.h18
-rw-r--r--src/2geom/linear.h2
-rw-r--r--src/2geom/nearest-point.cpp278
-rw-r--r--src/2geom/nearest-point.h122
-rw-r--r--src/2geom/numeric/linear_system.h89
-rw-r--r--src/2geom/numeric/matrix.h207
-rw-r--r--src/2geom/numeric/vector.h184
-rw-r--r--src/2geom/path.cpp188
-rw-r--r--src/2geom/path.h211
-rw-r--r--src/2geom/piecewise.h4
-rw-r--r--src/2geom/point.cpp4
-rw-r--r--src/2geom/rect.h30
-rw-r--r--src/2geom/sbasis-to-bezier.cpp8
-rw-r--r--src/2geom/sbasis.cpp22
-rw-r--r--src/2geom/sbasis.h16
-rw-r--r--src/2geom/shape.cpp2
-rw-r--r--src/2geom/svg-elliptical-arc.cpp1146
-rw-r--r--src/2geom/svg-path.cpp2
-rw-r--r--src/2geom/transforms.cpp51
-rw-r--r--src/2geom/transforms.h5
26 files changed, 2347 insertions, 326 deletions
diff --git a/src/2geom/Makefile_insert b/src/2geom/Makefile_insert
index f70293a13..ef82b7525 100644
--- a/src/2geom/Makefile_insert
+++ b/src/2geom/Makefile_insert
@@ -34,6 +34,11 @@
2geom/linear.h \
2geom/matrix.cpp \
2geom/matrix.h \
+ 2geom/nearest-point.cpp \
+ 2geom/nearest-point.h \
+ 2geom/numeric/linear_system.h \
+ 2geom/numeric/matrix.h \
+ 2geom/numeric/vector.h \
2geom/ord.h \
2geom/path-intersection.cpp \
2geom/path-intersection.h \
diff --git a/src/2geom/angle.h b/src/2geom/angle.h
index c6a367d8f..e95300a17 100644
--- a/src/2geom/angle.h
+++ b/src/2geom/angle.h
@@ -45,6 +45,49 @@ inline double deg_to_rad(double deg) { return deg*M_PI/180.0;}
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
diff --git a/src/2geom/basic-intersection.cpp b/src/2geom/basic-intersection.cpp
index 28b3c6f20..97c4c6e5c 100644
--- a/src/2geom/basic-intersection.cpp
+++ b/src/2geom/basic-intersection.cpp
@@ -61,7 +61,7 @@ find_intersections( vector<Geom::Point> const & A,
std::vector<std::pair<double, double> >
find_self_intersections(OldBezier const &Sb) {
- throwNotImplemented();
+ THROW_NOTIMPLEMENTED();
}
std::vector<std::pair<double, double> >
diff --git a/src/2geom/bezier-utils.cpp b/src/2geom/bezier-utils.cpp
index 76c90a915..59aac8951 100644
--- a/src/2geom/bezier-utils.cpp
+++ b/src/2geom/bezier-utils.cpp
@@ -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;
}
}
@@ -216,7 +216,7 @@ bezier_fit_cubic_full(Point bezier[], int split_points[],
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];
@@ -619,7 +619,7 @@ NewtonRaphsonRootFind(BezierCurve const Q, Point const &P, double const u)
}
}
- if (!is_finite(improved_u)) {
+ if (!IS_FINITE(improved_u)) {
improved_u = u;
} else if ( improved_u < 0.0 ) {
improved_u = 0.0;
@@ -853,7 +853,7 @@ chord_length_parameterize(Point const d[], double u[], unsigned const len)
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 9b7d8fb17..289a67729 100644
--- a/src/2geom/bezier.h
+++ b/src/2geom/bezier.h
@@ -144,7 +144,7 @@ public:
}
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 88ecfc51b..2749ec63a 100644
--- a/src/2geom/exception.h
+++ b/src/2geom/exception.h
@@ -36,7 +36,7 @@
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) {
@@ -53,7 +53,7 @@ public:
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.
@@ -65,14 +65,14 @@ public:
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 :)
@@ -82,29 +82,29 @@ public:
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 5b068e606..decebc7d2 100644
--- a/src/2geom/isnan.h
+++ b/src/2geom/isnan.h
@@ -27,15 +27,15 @@
*/
#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,
@@ -44,13 +44,13 @@
#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 0778039d3..93cbc1b04 100644
--- a/src/2geom/linear.h
+++ b/src/2geom/linear.h
@@ -88,7 +88,7 @@ public:
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
new file mode 100644
index 000000000..074de1cb3
--- /dev/null
+++ b/src/2geom/nearest-point.cpp
@@ -0,0 +1,278 @@
+/*
+ * nearest point routines for D2<SBasis> and Piecewise<D2<SBasis>>
+ *
+ * Authors:
+ *
+ * Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2007-2008 authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+
+#include "nearest-point.h"
+
+namespace Geom
+{
+
+////////////////////////////////////////////////////////////////////////////////
+// D2<SBasis> versions
+
+/*
+ * Return the parameter t of a nearest point on the portion of the curve "c",
+ * related to the interval [from, to], to the point "p".
+ * The needed curve derivative "dc" is passed as parameter.
+ * The function return the first nearest point to "p" that is found.
+ */
+
+double nearest_point( Point const& p,
+ D2<SBasis> const& c,
+ D2<SBasis> const& dc,
+ double from, double to )
+{
+ if ( from > to ) std::swap(from, to);
+ if ( from < 0 || to > 1 )
+ {
+ THROW_RANGEERROR("[from,to] interval out of bounds");
+ }
+
+ SBasis dd = dot(c - p, dc);
+ std::vector<double> zeros = Geom::roots(dd);
+
+ double closest = from;
+ double min_dist_sq = L2sq(c(from) - p);
+ double distsq;
+ for ( unsigned int i = 0; i < zeros.size(); ++i )
+ {
+ distsq = L2sq(c(zeros[i]) - p);
+ if ( min_dist_sq > L2sq(c(zeros[i]) - p) )
+ {
+ closest = zeros[i];
+ min_dist_sq = distsq;
+ }
+ }
+ if ( min_dist_sq > L2sq( c(to) - p ) )
+ closest = to;
+ return closest;
+
+}
+
+/*
+ * Return the parameters t of all the nearest points on the portion of
+ * the curve "c", related to the interval [from, to], to the point "p".
+ * The needed curve derivative "dc" is passed as parameter.
+ */
+
+std::vector<double>
+all_nearest_points( Point const& p,
+ D2<SBasis> const& c,
+ D2<SBasis> const& dc,
+ double from, double to )
+{
+ std::swap(from, to);
+ if ( from > to ) std::swap(from, to);
+ if ( from < 0 || to > 1 )
+ {
+ THROW_RANGEERROR("[from,to] interval out of bounds");
+ }
+
+ std::vector<double> result;
+ SBasis dd = dot(c - p, Geom::derivative(c));
+
+ std::vector<double> zeros = Geom::roots(dd);
+ std::vector<double> candidates;
+ candidates.push_back(from);
+ candidates.insert(candidates.end(), zeros.begin(), zeros.end());
+ candidates.push_back(to);
+ std::vector<double> distsq;
+ distsq.reserve(candidates.size());
+ for ( unsigned int i = 0; i < candidates.size(); ++i )
+ {
+ distsq.push_back( L2sq(c(candidates[i]) - p) );
+ }
+ unsigned int closest = 0;
+ double dsq = distsq[0];
+ for ( unsigned int i = 1; i < candidates.size(); ++i )
+ {
+ if ( dsq > distsq[i] )
+ {
+ closest = i;
+ dsq = distsq[i];
+ }
+ }
+ for ( unsigned int i = 0; i < candidates.size(); ++i )
+ {
+ if( distsq[closest] == distsq[i] )
+ {
+ result.push_back(candidates[i]);
+ }
+ }
+ return result;
+}
+
+
+////////////////////////////////////////////////////////////////////////////////
+// Piecewise< D2<SBasis> > versions
+
+
+double nearest_point( Point const& p,
+ Piecewise< D2<SBasis> > const& c,
+ double from, double to )
+{
+ if ( from > to ) std::swap(from, to);
+ if ( from < c.cuts[0] || to > c.cuts[c.size()] )
+ {
+ THROW_RANGEERROR("[from,to] interval out of bounds");
+ }
+
+ unsigned int si = c.segN(from);
+ unsigned int ei = c.segN(to);
+ if ( si == ei )
+ {
+ double nearest=
+ nearest_point(p, c[si], c.segT(from, si), c.segT(to, si));
+ return c.mapToDomain(nearest, si);
+ }
+ double t;
+ double nearest = nearest_point(p, c[si], c.segT(from, si));
+ unsigned int ni = si;
+ double dsq;
+ double mindistsq = distanceSq(p, c[si](nearest));
+ Rect bb;
+ for ( unsigned int i = si + 1; i < ei; ++i )
+ {
+ bb = bounds_fast(c[i]);
+ dsq = distanceSq(p, bb);
+ if ( mindistsq <= dsq ) continue;
+ t = nearest_point(p, c[i]);
+ dsq = distanceSq(p, c[i](t));
+ if ( mindistsq > dsq )
+ {
+ nearest = t;
+ ni = i;
+ mindistsq = dsq;
+ }
+ }
+ bb = bounds_fast(c[ei]);
+ dsq = distanceSq(p, bb);
+ if ( mindistsq > dsq )
+ {
+ t = nearest_point(p, c[ei], 0, c.segT(to, ei));
+ dsq = distanceSq(p, c[ei](t));
+ if ( mindistsq > dsq )
+ {
+ nearest = t;
+ ni = ei;
+ }
+ }
+ return c.mapToDomain(nearest, ni);
+}
+
+std::vector<double>
+all_nearest_points( Point const& p,
+ Piecewise< D2<SBasis> > const& c,
+ double from, double to )
+{
+ if ( from > to ) std::swap(from, to);
+ if ( from < c.cuts[0] || to > c.cuts[c.size()] )
+ {
+ THROW_RANGEERROR("[from,to] interval out of bounds");
+ }
+
+ unsigned int si = c.segN(from);
+ unsigned int ei = c.segN(to);
+ if ( si == ei )
+ {
+ std::vector<double> all_nearest =
+ all_nearest_points(p, c[si], c.segT(from, si), c.segT(to, si));
+ for ( unsigned int i = 0; i < all_nearest.size(); ++i )
+ {
+ all_nearest[i] = c.mapToDomain(all_nearest[i], si);
+ }
+ return all_nearest;
+ }
+ std::vector<double> all_t;
+ std::vector< std::vector<double> > all_np;
+ all_np.push_back( all_nearest_points(p, c[si], c.segT(from, si)) );
+ std::vector<unsigned int> ni;
+ ni.push_back(si);
+ double dsq;
+ double mindistsq = distanceSq( p, c[si](all_np.front().front()) );
+ Rect bb;
+ for ( unsigned int i = si + 1; i < ei; ++i )
+ {
+ bb = bounds_fast(c[i]);
+ dsq = distanceSq(p, bb);
+ if ( mindistsq < dsq ) continue;
+ all_t = all_nearest_points(p, c[i]);
+ dsq = distanceSq( p, c[i](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 = bounds_fast(c[ei]);
+ dsq = distanceSq(p, bb);
+ if ( mindistsq >= dsq )
+ {
+ all_t = all_nearest_points(p, c[ei], 0, c.segT(to, ei));
+ dsq = distanceSq( p, c[ei](all_t.front()) );
+ if ( mindistsq > dsq )
+ {
+ for ( unsigned int i = 0; i < all_t.size(); ++i )
+ {
+ all_t[i] = c.mapToDomain(all_t[i], ei);
+ }
+ 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( c.mapToDomain(all_np[i][j], ni[i]) );
+ }
+ }
+ return all_nearest;
+}
+
+} // end namespace Geom
+
+
diff --git a/src/2geom/nearest-point.h b/src/2geom/nearest-point.h
new file mode 100644
index 000000000..43b76c3ab
--- /dev/null
+++ b/src/2geom/nearest-point.h
@@ -0,0 +1,122 @@
+/*
+ * nearest point routines for D2<SBasis> and Piecewise<D2<SBasis>>
+ *
+ *
+ * Authors:
+ *
+ * Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2007-2008 authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+
+#ifndef _NEAREST_POINT_H_
+#define _NEAREST_POINT_H_
+
+
+#include <vector>
+
+#include "d2.h"
+#include "piecewise.h"
+#include "exception.h"
+
+
+
+namespace Geom
+{
+
+////////////////////////////////////////////////////////////////////////////////
+// D2<SBasis> versions
+
+/*
+ * Return the parameter t of a nearest point on the portion of the curve "c",
+ * related to the interval [from, to], to the point "p".
+ * The needed curve derivative "dc" is passed as parameter.
+ * The function return the first nearest point to "p" that is found.
+ */
+double nearest_point( Point const& p,
+ D2<SBasis> const& c, D2<SBasis> const& dc,
+ double from = 0, double to = 1 );
+
+inline
+double nearest_point( Point const& p,
+ D2<SBasis> const& c,
+ double from = 0, double to = 1 )
+{
+ return nearest_point(p, c, Geom::derivative(c), from, to);
+}
+
+/*
+ * Return the parameters t of all the nearest points on the portion of
+ * the curve "c", related to the interval [from, to], to the point "p".
+ * The needed curve derivative "dc" is passed as parameter.
+ */
+std::vector<double>
+all_nearest_points( Point const& p,
+ D2<SBasis> const& c, D2<SBasis> const& dc,
+ double from = 0, double to = 1 );
+
+inline
+std::vector<double>
+all_nearest_points( Point const& p,
+ D2<SBasis> const& c,
+ double from = 0, double to = 1 )
+{
+ return all_nearest_points(p, c, Geom::derivative(c), from, to);
+}
+
+
+////////////////////////////////////////////////////////////////////////////////
+// Piecewise< D2<SBasis> > versions
+
+double nearest_point( Point const& p,
+ Piecewise< D2<SBasis> > const& c,
+ double from, double to );
+
+inline
+double nearest_point( Point const& p, Piecewise< D2<SBasis> > const& c )
+{
+ return nearest_point(p, c, c.cuts[0], c.cuts[c.size()]);
+}
+
+
+std::vector<double>
+all_nearest_points( Point const& p,
+ Piecewise< D2<SBasis> > const& c,
+ double from, double to );
+
+inline
+std::vector<double>
+all_nearest_points( Point const& p, Piecewise< D2<SBasis> > const& c )
+{
+ return all_nearest_points(p, c, c.cuts[0], c.cuts[c.size()]);
+}
+
+} // end namespace Geom
+
+
+
+#endif /*_NEAREST_POINT_H_*/
diff --git a/src/2geom/numeric/linear_system.h b/src/2geom/numeric/linear_system.h
new file mode 100644
index 000000000..2b4963998
--- /dev/null
+++ b/src/2geom/numeric/linear_system.h
@@ -0,0 +1,89 @@
+#ifndef _NL_LINEAR_SYSTEM_H_
+#define _NL_LINEAR_SYSTEM_H_
+
+
+#include <cassert>
+
+#include <gsl/gsl_linalg.h>
+
+#include "numeric/matrix.h"
+#include "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_*/
diff --git a/src/2geom/numeric/matrix.h b/src/2geom/numeric/matrix.h
new file mode 100644
index 000000000..bdfab75ef
--- /dev/null
+++ b/src/2geom/numeric/matrix.h
@@ -0,0 +1,207 @@
+#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_*/
diff --git a/src/2geom/numeric/vector.h b/src/2geom/numeric/vector.h
new file mode 100644
index 000000000..136b2cc3f
--- /dev/null
+++ b/src/2geom/numeric/vector.h
@@ -0,0 +1,184 @@
+#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);
+ }
+
+ 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_*/
diff --git a/src/2geom/path.cpp b/src/2geom/path.cpp
index 01c1f2b8e..fdfa77c79 100644
--- a/src/2geom/path.cpp
+++ b/src/2geom/path.cpp
@@ -33,6 +33,7 @@
namespace Geom {
+
int CurveHelpers::root_winding(Curve const &c, Point p) {
std::vector<double> ts = c.roots(p[Y], Y);
@@ -74,6 +75,22 @@ int CurveHelpers::root_winding(Curve const &c, Point p) {
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_);
@@ -106,6 +123,167 @@ iter inc(iter const &x, unsigned n) {
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);
@@ -145,7 +323,7 @@ const double eps = .1;
void Path::append(Curve const &curve) {
if ( curves_.front() != final_ && !are_near(curve.initialPoint(), (*final_)[0], eps) ) {
- throwContinuityError();
+ THROW_CONTINUITYERROR();
}
do_append(curve.duplicate());
}
@@ -154,7 +332,7 @@ 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) ) {
- throwContinuityError();
+ THROW_CONTINUITYERROR();
}
}
}
@@ -206,17 +384,17 @@ void Path::check_continuity(Sequence::iterator first_replaced,
if ( first != last ) {
if ( first_replaced != curves_.begin() ) {
if ( !are_near( (*first_replaced)->initialPoint(), (*first)->initialPoint(), eps ) ) {
- 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 414d69755..ed15260fc 100644
--- a/src/2geom/path.h
+++ b/src/2geom/path.h
@@ -44,6 +44,7 @@
#include "bezier.h"
#include "crossing.h"
#include "utils.h"
+#include "nearest-point.h"
namespace Geom {
@@ -81,6 +82,19 @@ public:
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;
@@ -91,9 +105,11 @@ public:
};
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()) {}
@@ -116,6 +132,17 @@ public:
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));
@@ -135,8 +162,10 @@ public:
template <unsigned order>
class BezierCurve : public Curve {
+
private:
D2<Bezier > inner;
+
public:
template <unsigned required_degree>
static void assert_degree(BezierCurve<required_degree> const *) {}
@@ -205,7 +234,12 @@ public:
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]);
@@ -269,6 +303,8 @@ typedef BezierCurve<1> LineSegment;
typedef BezierCurve<2> QuadraticBezier;
typedef BezierCurve<3> CubicBezier;
+template<>
+double LineSegment::nearestPoint(Point const& p, double from, double to) const;
@@ -276,7 +312,13 @@ class SVGEllipticalArc : public Curve
{
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,
@@ -286,24 +328,30 @@ class SVGEllipticalArc : public Curve
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];
}
@@ -333,7 +381,7 @@ class SVGEllipticalArc : public Curve
return m_end_angle;
}
- double ray(Dim2 i) const
+ double ray(unsigned int i) const
{
return (i == 0) ? m_rx : m_ry;
}
@@ -374,28 +422,34 @@ class SVGEllipticalArc : public Curve
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
@@ -403,36 +457,43 @@ class SVGEllipticalArc : public Curve
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
@@ -460,7 +521,6 @@ class SVGEllipticalArc : public Curve
return rarc;
}
- private:
double sweep_angle() const
{
Coord d = end_angle() - start_angle();
@@ -469,12 +529,12 @@ class SVGEllipticalArc : public Curve
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;
bool m_large_arc, m_sweep;
@@ -657,20 +717,42 @@ public:
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 {
@@ -682,7 +764,28 @@ public:
}
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 {
@@ -704,7 +807,7 @@ public:
}
return ret;
}
-
+
void insert(iterator pos, Curve const &curve) {
Sequence source(1, curve.duplicate());
try {
@@ -909,7 +1012,7 @@ class PathPortion : public Curve {
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 1d0fd93ba..d548ea0e9 100644
--- a/src/2geom/piecewise.h
+++ b/src/2geom/piecewise.h
@@ -104,7 +104,7 @@ class Piecewise {
}
//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.
@@ -147,7 +147,7 @@ class Piecewise {
//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 1e2a3463f..66eaf8ca6 100644
--- a/src/2geom/point.cpp
+++ b/src/2geom/point.cpp
@@ -18,7 +18,7 @@ namespace Geom {
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;
@@ -71,7 +71,7 @@ Coord L1(Point const &p) {
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 86f9ec140..c89946606 100644
--- a/src/2geom/rect.h
+++ b/src/2geom/rect.h
@@ -153,8 +153,38 @@ inline boost::optional<Rect> intersect(Rect const & a, Rect const & b) {
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
/*
diff --git a/src/2geom/sbasis-to-bezier.cpp b/src/2geom/sbasis-to-bezier.cpp
index 497091d1c..9843b78d7 100644
--- a/src/2geom/sbasis-to-bezier.cpp
+++ b/src/2geom/sbasis-to-bezier.cpp
@@ -112,7 +112,7 @@ D2<Bezier<order> > sbasis_to_bezier(D2<SBasis> const &B) {
// 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) {
@@ -144,8 +144,8 @@ void
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 7157bc885..c60132043 100644
--- a/src/2geom/sbasis.cpp
+++ b/src/2geom/sbasis.cpp
@@ -68,6 +68,28 @@ bool SBasis::isFinite() const {
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 82ff989e7..7f0d88f8c 100644
--- a/src/2geom/sbasis.h
+++ b/src/2geom/sbasis.h
@@ -113,21 +113,7 @@ public:
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 82ae41145..a7a40066e 100644
--- 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) {
diff --git a/src/2geom/svg-elliptical-arc.cpp b/src/2geom/svg-elliptical-arc.cpp
index 6583c948f..b3420fba6 100644
--- a/src/2geom/svg-elliptical-arc.cpp
+++ b/src/2geom/svg-elliptical-arc.cpp
@@ -1,215 +1,931 @@
-/*
- * 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 :
-
-
+/*
+ * 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"
+#include "angle.h"
+
+#include <gsl/gsl_poly.h>
+
+#include <cfloat>
+
+
+
+
+namespace Geom
+{
+
+
+Rect SVGEllipticalArc::boundsExact() const
+{
+ std::vector<double> extremes(4);
+ double cosrot = std::cos(rotation_angle());
+ double sinrot = std::sin(rotation_angle());
+ extremes[0] = std::atan2( -ray(Y) * sinrot, ray(X) * cosrot );
+ extremes[1] = extremes[0] + M_PI;
+ if ( extremes[0] < 0 ) extremes[0] += 2*M_PI;
+ extremes[2] = std::atan2( ray(Y) * cosrot, ray(X) * sinrot );
+ extremes[3] = extremes[2] + M_PI;
+ if ( extremes[2] < 0 ) extremes[2] += 2*M_PI;
+
+
+ std::vector<double>arc_extremes(4);
+ arc_extremes[0] = initialPoint()[X];
+ arc_extremes[1] = finalPoint()[X];
+ if ( arc_extremes[0] < arc_extremes[1] )
+ std::swap(arc_extremes[0], arc_extremes[1]);
+ arc_extremes[2] = initialPoint()[Y];
+ arc_extremes[3] = finalPoint()[Y];
+ if ( arc_extremes[2] < arc_extremes[3] )
+ std::swap(arc_extremes[2], arc_extremes[3]);
+
+
+ if ( start_angle() < end_angle() )
+ {
+ if ( sweep_flag() )
+ {
+ for ( unsigned int i = 0; i < extremes.size(); ++i )
+ {
+ if ( start_angle() < extremes[i] && extremes[i] < end_angle() )
+ {
+ arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];
+ }
+ }
+ }
+ else
+ {
+ for ( unsigned int i = 0; i < extremes.size(); ++i )
+ {
+ if ( start_angle() > extremes[i] || extremes[i] > end_angle() )
+ {
+ arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];
+ }
+ }
+ }
+ }
+ else
+ {
+ if ( sweep_flag() )
+ {
+ for ( unsigned int i = 0; i < extremes.size(); ++i )
+ {
+ if ( start_angle() < extremes[i] || extremes[i] < end_angle() )
+ {
+ arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];
+ }
+ }
+ }
+ else
+ {
+ for ( unsigned int i = 0; i < extremes.size(); ++i )
+ {
+ if ( start_angle() > extremes[i] && extremes[i] > end_angle() )
+ {
+ arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];
+ }
+ }
+ }
+ }
+
+ return Rect( Point(arc_extremes[1], arc_extremes[3]) ,
+ Point(arc_extremes[0], arc_extremes[2]) );
+
+}
+
+
+std::vector<double>
+SVGEllipticalArc::roots(double v, Dim2 d) const
+{
+ if ( d > Y )
+ {
+ THROW_RANGEERROR("dimention out of range");
+ }
+
+ std::vector<double> sol;
+ if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )
+ {
+ if ( center(d) == v )
+ sol.push_back(0);
+ return sol;
+ }
+
+ const char* msg[2][2] =
+ {
+ { "d == X; ray(X) == 0; "
+ "s = (v - center(X)) / ( -ray(Y) * std::sin(rotation_angle()) ); "
+ "s should be contained in [-1,1]",
+ "d == X; ray(Y) == 0; "
+ "s = (v - center(X)) / ( ray(X) * std::cos(rotation_angle()) ); "
+ "s should be contained in [-1,1]"
+ },
+ { "d == Y; ray(X) == 0; "
+ "s = (v - center(X)) / ( ray(Y) * std::cos(rotation_angle()) ); "
+ "s should be contained in [-1,1]",
+ "d == Y; ray(Y) == 0; "
+ "s = (v - center(X)) / ( ray(X) * std::sin(rotation_angle()) ); "
+ "s should be contained in [-1,1]"
+ },
+ };
+
+ for ( unsigned int dim = 0; dim < 2; ++dim )
+ {
+ if ( are_near(ray(dim), 0) )
+ {
+
+ if ( initialPoint()[d] == v && finalPoint()[d] == v )
+ {
+ THROW_EXCEPTION("infinite solutions");
+ }
+ if ( (initialPoint()[d] < finalPoint()[d])
+ && (initialPoint()[d] > v || finalPoint()[d] < v) )
+ {
+ return sol;
+ }
+ if ( (initialPoint()[d] > finalPoint()[d])
+ && (finalPoint()[d] > v || initialPoint()[d] < v) )
+ {
+ return sol;
+ }
+ double ray_prj;
+ switch(d)
+ {
+ case X:
+ switch(dim)
+ {
+ case X: ray_prj = -ray(Y) * std::sin(rotation_angle());
+ break;
+ case Y: ray_prj = ray(X) * std::cos(rotation_angle());
+ break;
+ }
+ break;
+ case Y:
+ switch(dim)
+ {
+ case X: ray_prj = ray(Y) * std::cos(rotation_angle());
+ break;
+ case Y: ray_prj = ray(X) * std::sin(rotation_angle());
+ break;
+ }
+ break;
+ }
+
+ double s = (v - center(d)) / ray_prj;
+ if ( s < -1 || s > 1 )
+ {
+ THROW_LOGICALERROR(msg[d][dim]);
+ }
+ switch(dim)
+ {
+ case X:
+ s = std::asin(s); // return a value in [-PI/2,PI/2]
+ if ( logical_xor( sweep_flag(), are_near(start_angle(), M_PI/2) ) )
+ {
+ if ( s < 0 ) s += 2*M_PI;
+ }
+ else
+ {
+ s = M_PI - s;
+ if (!(s < 2*M_PI) ) s -= 2*M_PI;
+ }
+ break;
+ case Y:
+ s = std::acos(s); // return a value in [0,PI]
+ if ( logical_xor( sweep_flag(), are_near(start_angle(), 0) ) )
+ {
+ s = 2*M_PI - s;
+ if ( !(s < 2*M_PI) ) s -= 2*M_PI;
+ }
+ break;
+ }
+
+ //std::cerr << "s = " << rad_to_deg(s);
+ s = map_to_01(s);
+ //std::cerr << " -> t: " << s << std::endl;
+ if ( !(s < 0 || s > 1) )
+ sol.push_back(s);
+ return sol;
+ }
+ }
+
+ double rotx, roty;
+ switch(d)
+ {
+ case X:
+ rotx = std::cos(rotation_angle());
+ roty = -std::sin(rotation_angle());
+ break;
+ case Y:
+ rotx = std::sin(rotation_angle());
+ roty = std::cos(rotation_angle());
+ break;
+ }
+ double rxrotx = ray(X) * rotx;
+ double c_v = center(d) - v;
+
+ double a = -rxrotx + c_v;
+ double b = ray(Y) * roty;
+ double c = rxrotx + c_v;
+ //std::cerr << "a = " << a << std::endl;
+ //std::cerr << "b = " << b << std::endl;
+ //std::cerr << "c = " << c << std::endl;
+
+ if ( are_near(a,0) )
+ {
+ sol.push_back(M_PI);
+ if ( !are_near(b,0) )
+ {
+ double s = 2 * std::atan(-c/(2*b));
+ if ( s < 0 ) s += 2*M_PI;
+ sol.push_back(s);
+ }
+ }
+ else
+ {
+ double delta = b * b - a * c;
+ //std::cerr << "delta = " << delta << std::endl;
+ if ( are_near(delta, 0) )
+ {
+ double s = 2 * std::atan(-b/a);
+ if ( s < 0 ) s += 2*M_PI;
+ sol.push_back(s);
+ }
+ else if ( delta > 0 )
+ {
+ double sq = std::sqrt(delta);
+ double s = 2 * std::atan( (-b - sq) / a );
+ if ( s < 0 ) s += 2*M_PI;
+ sol.push_back(s);
+ s = 2 * std::atan( (-b + sq) / a );
+ if ( s < 0 ) s += 2*M_PI;
+ sol.push_back(s);
+ }
+ }
+
+ std::vector<double> arc_sol;
+ for (unsigned int i = 0; i < sol.size(); ++i )
+ {
+ //std::cerr << "s = " << rad_to_deg(sol[i]);
+ sol[i] = map_to_01(sol[i]);
+ //std::cerr << " -> t: " << sol[i] << std::endl;
+ if ( !(sol[i] < 0 || sol[i] > 1) )
+ arc_sol.push_back(sol[i]);
+ }
+ return arc_sol;
+
+
+// return SBasisCurve(toSBasis()).roots(v, d);
+}
+
+// D(E(t,C),t) = E(t+PI/2,O)
+Curve* SVGEllipticalArc::derivative() const
+{
+ SVGEllipticalArc* result = new SVGEllipticalArc(*this);
+ result->m_center[X] = result->m_center[Y] = 0;
+ result->m_start_angle += M_PI/2;
+ if( !( result->m_start_angle < 2*M_PI ) )
+ {
+ result->m_start_angle -= 2*M_PI;
+ }
+ result->m_end_angle += M_PI/2;
+ if( !( result->m_end_angle < 2*M_PI ) )
+ {
+ result->m_end_angle -= 2*M_PI;
+ }
+ result->m_initial_point = result->pointAtAngle( result->start_angle() );
+ result->m_final_point = result->pointAtAngle( result->end_angle() );
+ return result;
+
+}
+
+std::vector<Point>
+SVGEllipticalArc::pointAndDerivatives(Coord t, unsigned int n) const
+{
+ std::vector<Point> result;
+ result.reserve(n);
+ double angle = map_unit_interval_on_circular_arc(t, start_angle(),
+ end_angle(), sweep_flag());
+ SVGEllipticalArc ea(*this);
+ ea.m_center = Point(0,0);
+ unsigned int m = std::min(n, 4u);
+ for ( unsigned int i = 0; i < m; ++i )
+ {
+ result.push_back( ea.pointAtAngle(angle) );
+ angle += M_PI/2;
+ if ( !(angle < 2*M_PI) ) angle -= 2*M_PI;
+ }
+ m = n / 4;
+ for ( unsigned int i = 1; i < m; ++i )
+ {
+ for ( unsigned int j = 0; j < 4; ++j )
+ result.push_back( result[j] );
+ }
+ m = n - 4 * m;
+ for ( unsigned int i = 0; i < m; ++i )
+ {
+ result.push_back( result[i] );
+ }
+ if ( !result.empty() ) // n != 0
+ result[0] = pointAtAngle(angle);
+ return result;
+}
+
+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;
+}
+
+
+bool SVGEllipticalArc::containsAngle(Coord angle) const
+{
+ if ( sweep_flag() )
+ if ( start_angle() < end_angle() )
+ if ( !( angle < start_angle() || angle > end_angle() ) )
+ return true;
+ else
+ return false;
+ else
+ if ( !( angle < start_angle() && angle > end_angle() ) )
+ return true;
+ else
+ return false;
+ else
+ if ( start_angle() > end_angle() )
+ if ( !( angle > start_angle() || angle < end_angle() ) )
+ return true;
+ else
+ return false;
+ else
+ if ( !( angle > start_angle() && angle < end_angle() ) )
+ return true;
+ else
+ return false;
+}
+
+
+double SVGEllipticalArc::valueAtAngle(Coord t, Dim2 d) const
+{
+ double sin_rot_angle = std::sin(rotation_angle());
+ double cos_rot_angle = std::cos(rotation_angle());
+ if ( d == X )
+ {
+ return ray(X) * cos_rot_angle * std::cos(t)
+ - ray(Y) * sin_rot_angle * std::sin(t)
+ + center(X);
+ }
+ else if ( d == Y )
+ {
+ return ray(X) * sin_rot_angle * std::cos(t)
+ + ray(Y) * cos_rot_angle * std::sin(t)
+ + center(Y);
+ }
+ THROW_RANGEERROR("dimension parameter out of range");
+}
+
+
+Curve* SVGEllipticalArc::portion(double f, double t) const
+{
+ if (f < 0) f = 0;
+ if (f > 1) f = 1;
+ if (t < 0) t = 0;
+ if (t > 1) t = 1;
+ if ( are_near(f, t) )
+ {
+ SVGEllipticalArc* arc = new SVGEllipticalArc();
+ arc->m_center = arc->m_initial_point = arc->m_final_point = pointAt(f);
+ arc->m_start_angle = arc->m_end_angle = m_start_angle;
+ arc->m_rot_angle = m_rot_angle;
+ arc->m_sweep = m_sweep;
+ arc->m_large_arc = m_large_arc;
+ }
+ SVGEllipticalArc* arc = new SVGEllipticalArc( *this );
+ arc->m_initial_point = pointAt(f);
+ arc->m_final_point = pointAt(t);
+ double sa = sweep_flag() ? sweep_angle() : -sweep_angle();
+ arc->m_start_angle = m_start_angle + sa * f;
+ if ( !(arc->m_start_angle < 2*M_PI) )
+ arc->m_start_angle -= 2*M_PI;
+ if ( arc->m_start_angle < 0 )
+ arc->m_start_angle += 2*M_PI;
+ arc->m_end_angle = m_start_angle + sa * t;
+ if ( !(arc->m_end_angle < 2*M_PI) )
+ arc->m_end_angle -= 2*M_PI;
+ if ( arc->m_end_angle < 0 )
+ arc->m_end_angle += 2*M_PI;
+ if ( f > t ) arc->m_sweep = !sweep_flag();
+ if ( large_arc_flag() && (arc->sweep_angle() < M_PI) )
+ arc->m_large_arc = false;
+ return arc;
+}
+
+// NOTE: doesn't work with 360 deg arcs
+void SVGEllipticalArc::calculate_center_and_extreme_angles()
+{
+ if ( are_near(initialPoint(), finalPoint()) )
+ {
+ if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )
+ {
+ m_start_angle = m_end_angle = 0;
+ m_center = initialPoint();
+ return;
+ }
+ else
+ {
+ THROW_RANGEERROR("initial and final point are the same");
+ }
+ }
+ if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )
+ { // but initialPoint != finalPoint
+ THROW_RANGEERROR(
+ "there is no ellipse that satisfies the given constraints: "
+ "ray(X) == 0 && ray(Y) == 0 but initialPoint != finalPoint"
+ );
+ }
+ if ( are_near(ray(Y), 0) )
+ {
+ Point v = initialPoint() - finalPoint();
+ if ( are_near(L2sq(v), 4*ray(X)*ray(X)) )
+ {
+ double angle = std::atan2(v[Y], v[X]);
+ if (angle < 0) angle += 2*M_PI;
+ if ( are_near( angle, rotation_angle() ) )
+ {
+ m_start_angle = 0;
+ m_end_angle = M_PI;
+ m_center = v/2 + finalPoint();
+ return;
+ }
+ angle -= M_PI;
+ if ( angle < 0 ) angle += 2*M_PI;
+ if ( are_near( angle, rotation_angle() ) )
+ {
+ m_start_angle = M_PI;
+ m_end_angle = 0;
+ m_center = v/2 + finalPoint();
+ return;
+ }
+ THROW_RANGEERROR(
+ "there is no ellipse that satisfies the given constraints: "
+ "ray(Y) == 0 "
+ "and slope(initialPoint - finalPoint) != rotation_angle "
+ "and != rotation_angle + PI"
+ );
+ }
+ if ( L2sq(v) > 4*ray(X)*ray(X) )
+ {
+ THROW_RANGEERROR(
+ "there is no ellipse that satisfies the given constraints: "
+ "ray(Y) == 0 and distance(initialPoint, finalPoint) > 2*ray(X)"
+ );
+ }
+ else
+ {
+ THROW_RANGEERROR(
+ "there is infinite ellipses that satisfy the given constraints: "
+ "ray(Y) == 0 and distance(initialPoint, finalPoint) < 2*ray(X)"
+ );
+ }
+
+ }
+
+ if ( are_near(ray(X), 0) )
+ {
+ Point v = initialPoint() - finalPoint();
+ if ( are_near(L2sq(v), 4*ray(Y)*ray(Y)) )
+ {
+ double angle = std::atan2(v[Y], v[X]);
+ if (angle < 0) angle += 2*M_PI;
+ double rot_angle = rotation_angle() + M_PI/2;
+ if ( !(rot_angle < 2*M_PI) ) rot_angle -= 2*M_PI;
+ if ( are_near( angle, rot_angle ) )
+ {
+ m_start_angle = M_PI/2;
+ m_end_angle = 3*M_PI/2;
+ m_center = v/2 + finalPoint();
+ return;
+ }
+ angle -= M_PI;
+ if ( angle < 0 ) angle += 2*M_PI;
+ if ( are_near( angle, rot_angle ) )
+ {
+ m_start_angle = 3*M_PI/2;
+ m_end_angle = M_PI/2;
+ m_center = v/2 + finalPoint();
+ return;
+ }
+ THROW_RANGEERROR(
+ "there is no ellipse that satisfies the given constraints: "
+ "ray(X) == 0 "
+ "and slope(initialPoint - finalPoint) != rotation_angle + PI/2 "
+ "and != rotation_angle + (3/2)*PI"
+ );
+ }
+ if ( L2sq(v) > 4*ray(Y)*ray(Y) )
+ {
+ THROW_RANGEERROR(
+ "there is no ellipse that satisfies the given constraints: "
+ "ray(X) == 0 and distance(initialPoint, finalPoint) > 2*ray(Y)"
+ );
+ }
+ else
+ {
+ THROW_RANGEERROR(
+ "there is infinite ellipses that satisfy the given constraints: "
+ "ray(X) == 0 and distance(initialPoint, finalPoint) < 2*ray(Y)"
+ );
+ }
+
+ }
+
+ double sin_rot_angle = std::sin(rotation_angle());
+ double cos_rot_angle = std::cos(rotation_angle());
+
+ 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) )
+ THROW_RANGEERROR(
+ "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) )
+ THROW_RANGEERROR(
+ "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::map_to_02PI(Coord t) const
+{
+ if ( sweep_flag() )
+ {
+ Coord angle = start_angle() + sweep_angle() * t;
+ if ( !(angle < 2*M_PI) )
+ angle -= 2*M_PI;
+ return angle;
+ }
+ else
+ {
+ Coord angle = start_angle() - sweep_angle() * t;
+ if ( angle < 0 ) angle += 2*M_PI;
+ return angle;
+ }
+}
+
+Coord SVGEllipticalArc::map_to_01(Coord angle) const
+{
+ return map_circular_arc_on_unit_interval(angle, start_angle(),
+ end_angle(), sweep_flag());
+}
+
+
+std::vector<double> SVGEllipticalArc::
+allNearestPoints( Point const& p, double from, double to ) const
+{
+ if ( from > to ) std::swap(from, to);
+ if ( from < 0 || to > 1 )
+ {
+ THROW_RANGEERROR("[from,to] interval out of range");
+ }
+ std::vector<double> result;
+ if ( ( are_near(ray(X), 0) && are_near(ray(Y), 0) ) || are_near(from, to) )
+ {
+ result.push_back(from);
+ return result;
+ }
+ else if ( are_near(ray(X), 0) || are_near(ray(Y), 0) )
+ {
+ LineSegment seg(pointAt(from), pointAt(to));
+ Point np = seg.pointAt( seg.nearestPoint(p) );
+ if ( are_near(ray(Y), 0) )
+ {
+ if ( are_near(rotation_angle(), M_PI/2)
+ || are_near(rotation_angle(), 3*M_PI/2) )
+ {
+ result = roots(np[Y], Y);
+ }
+ else
+ {
+ result = roots(np[X], X);
+ }
+ }
+ else
+ {
+ if ( are_near(rotation_angle(), M_PI/2)
+ || are_near(rotation_angle(), 3*M_PI/2) )
+ {
+ result = roots(np[X], X);
+ }
+ else
+ {
+ result = roots(np[Y], Y);
+ }
+ }
+ return result;
+ }
+ else if ( are_near(ray(X), ray(Y)) )
+ {
+ Point r = p - center();
+ if ( are_near(r, Point(0,0)) )
+ {
+ THROW_EXCEPTION("infinite nearest points");
+ }
+ // TODO: implement case r != 0
+// Point np = ray(X) * unit_vector(r);
+// std::vector<double> solX = roots(np[X],X);
+// std::vector<double> solY = roots(np[Y],Y);
+// double t;
+// if ( are_near(solX[0], solY[0]) || are_near(solX[0], solY[1]))
+// {
+// t = solX[0];
+// }
+// else
+// {
+// t = solX[1];
+// }
+// if ( !(t < from || t > to) )
+// {
+// result.push_back(t);
+// }
+// else
+// {
+//
+// }
+ }
+
+ // solve the equation <D(E(t),t)|E(t)-p> == 0
+ // that provides min and max distance points
+ // on the ellipse E wrt the point p
+ // after the substitutions:
+ // cos(t) = (1 - s^2) / (1 + s^2)
+ // sin(t) = 2t / (1 + s^2)
+ // where s = tan(t/2)
+ // we get a 4th degree equation in s
+ /*
+ * ry s^4 ((-cy + py) Cos[Phi] + (cx - px) Sin[Phi]) +
+ * ry ((cy - py) Cos[Phi] + (-cx + px) Sin[Phi]) +
+ * 2 s^3 (rx^2 - ry^2 + (-cx + px) rx Cos[Phi] + (-cy + py) rx Sin[Phi]) +
+ * 2 s (-rx^2 + ry^2 + (-cx + px) rx Cos[Phi] + (-cy + py) rx Sin[Phi])
+ */
+
+ Point p_c = p - center();
+ double rx2_ry2 = (ray(X) - ray(Y)) * (ray(X) + ray(Y));
+ double cosrot = std::cos( rotation_angle() );
+ double sinrot = std::sin( rotation_angle() );
+ double expr1 = ray(X) * (p_c[X] * cosrot + p_c[Y] * sinrot);
+ double coeff[5];
+ coeff[4] = ray(Y) * ( p_c[Y] * cosrot - p_c[X] * sinrot );
+ coeff[3] = 2 * ( rx2_ry2 + expr1 );
+ coeff[2] = 0;
+ coeff[1] = 2 * ( -rx2_ry2 + expr1 );
+ coeff[0] = -coeff[4];
+
+// for ( unsigned int i = 0; i < 5; ++i )
+// std::cerr << "c[" << i << "] = " << coeff[i] << std::endl;
+
+ std::vector<double> real_sol;
+ // gsl_poly_complex_solve raises an error
+ // if the leading coefficient is zero
+ if ( are_near(coeff[4], 0) )
+ {
+ real_sol.push_back(0);
+ if ( !are_near(coeff[3], 0) )
+ {
+ double sq = -coeff[1] / coeff[3];
+ if ( sq > 0 )
+ {
+ double s = std::sqrt(sq);
+ real_sol.push_back(s);
+ real_sol.push_back(-s);
+ }
+ }
+ }
+ else
+ {
+ double sol[8];
+ gsl_poly_complex_workspace * w = gsl_poly_complex_workspace_alloc(5);
+ gsl_poly_complex_solve(coeff, 5, w, sol );
+ gsl_poly_complex_workspace_free(w);
+
+ for ( unsigned int i = 0; i < 4; ++i )
+ {
+ if ( sol[2*i+1] == 0 ) real_sol.push_back(sol[2*i]);
+ }
+ }
+
+ for ( unsigned int i = 0; i < real_sol.size(); ++i )
+ {
+ real_sol[i] = 2 * std::atan(real_sol[i]);
+ if ( real_sol[i] < 0 ) real_sol[i] += 2*M_PI;
+ }
+ // when s -> Infinity then <D(E)|E-p> -> 0 iff coeff[4] == 0
+ // so we add M_PI to the solutions being lim arctan(s) = PI when s->Infinity
+ if ( (real_sol.size() % 2) != 0 )
+ {
+ real_sol.push_back(M_PI);
+ }
+
+ double mindistsq1 = std::numeric_limits<double>::max();
+ double mindistsq2 = std::numeric_limits<double>::max();
+ double dsq;
+ unsigned int mi1, mi2;
+ for ( unsigned int i = 0; i < real_sol.size(); ++i )
+ {
+ dsq = distanceSq(p, pointAtAngle(real_sol[i]));
+ if ( mindistsq1 > dsq )
+ {
+ mindistsq2 = mindistsq1;
+ mi2 = mi1;
+ mindistsq1 = dsq;
+ mi1 = i;
+ }
+ else if ( mindistsq2 > dsq )
+ {
+ mindistsq2 = dsq;
+ mi2 = i;
+ }
+ }
+
+ double t = map_to_01( real_sol[mi1] );
+ if ( !(t < from || t > to) )
+ {
+ result.push_back(t);
+ }
+
+ bool second_sol = false;
+ t = map_to_01( real_sol[mi2] );
+ if ( real_sol.size() == 4 && !(t < from || t > to) )
+ {
+ if ( result.empty() || are_near(mindistsq1, mindistsq2) )
+ {
+ result.push_back(t);
+ second_sol = true;
+ }
+ }
+
+ // we need to test extreme points too
+ double dsq1 = distanceSq(p, pointAt(from));
+ double dsq2 = distanceSq(p, pointAt(to));
+ if ( second_sol )
+ {
+ if ( mindistsq2 > dsq1 )
+ {
+ result.clear();
+ result.push_back(from);
+ mindistsq2 = dsq1;
+ }
+ else if ( are_near(mindistsq2, dsq) )
+ {
+ result.push_back(from);
+ }
+ if ( mindistsq2 > dsq2 )
+ {
+ result.clear();
+ result.push_back(to);
+ }
+ else if ( are_near(mindistsq2, dsq2) )
+ {
+ result.push_back(to);
+ }
+
+ }
+ else
+ {
+ if ( result.empty() )
+ {
+ if ( are_near(dsq1, dsq2) )
+ {
+ result.push_back(from);
+ result.push_back(to);
+ }
+ else if ( dsq2 > dsq1 )
+ {
+ result.push_back(from);
+ }
+ else
+ {
+ result.push_back(to);
+ }
+ }
+ }
+
+ return result;
+}
+
+
+} // end namespace Geom
+
+
+/*
+ Local Variables:
+ mode:c++
+ c-file-style:"stroustrup"
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+ indent-tabs-mode:nil
+ fill-column:99
+ End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
+
+
diff --git a/src/2geom/svg-path.cpp b/src/2geom/svg-path.cpp
index 60ae957cb..c09e5c929 100644
--- a/src/2geom/svg-path.cpp
+++ b/src/2geom/svg-path.cpp
@@ -30,6 +30,7 @@
#include "sbasis-to-bezier.h"
#include "svg-path.h"
+#include "exception.h"
namespace Geom {
@@ -52,6 +53,7 @@ void output(QuadraticBezier const &curve, SVGPathSink &sink) {
void output(SVGEllipticalArc const &/*curve*/, SVGPathSink &/*sink*/) {
// FIXME
+ THROW_NOTIMPLEMENTED();
}
template <typename T>
diff --git a/src/2geom/transforms.cpp b/src/2geom/transforms.cpp
index bdb6bcad3..0a27c31f6 100644
--- a/src/2geom/transforms.cpp
+++ b/src/2geom/transforms.cpp
@@ -45,6 +45,57 @@ Matrix operator*(Matrix const &m, Scale const &s) {
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 94058bff0..9529d8922 100644
--- a/src/2geom/transforms.h
+++ b/src/2geom/transforms.h
@@ -112,6 +112,11 @@ Matrix operator*(Matrix const &m, Scale const &s);
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 */
generated by cgit v1.2.3 (git 2.25.1) at 2026-05-19 02:55:44 +0000