From: johanengelen Date: Sat, 10 May 2008 20:20:11 +0000 (+0000) Subject: update to latest 2geom. this adds gsl dependency, doesn't seem to make inskape execut... X-Git-Url: https://git.tokkee.org/?a=commitdiff_plain;h=d8fa3c4faade9a5a8e7f79450544b1925e1ade41;p=inkscape.git update to latest 2geom. this adds gsl dependency, doesn't seem to make inskape executable bigger --- diff --git a/build.xml b/build.xml index 5e38e1792..49fa4d76a 100644 --- a/build.xml +++ b/build.xml @@ -404,6 +404,7 @@ -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 -L${gtk}/perl/lib/CORE -lperl58 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 const & A, std::vector > find_self_intersections(OldBezier const &Sb) { - throwNotImplemented(); + THROW_NOTIMPLEMENTED(); } std::vector > 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 and Piecewise> + * + * Authors: + * + * Marco Cecchetti + * + * Copyright 2007-2008 authors + * + * This library is free software; you can redistribute it and/or + * modify it either under the terms of the GNU Lesser General Public + * License version 2.1 as published by the Free Software Foundation + * (the "LGPL") or, at your option, under the terms of the Mozilla + * Public License Version 1.1 (the "MPL"). If you do not alter this + * notice, a recipient may use your version of this file under either + * the MPL or the LGPL. + * + * You should have received a copy of the LGPL along with this library + * in the file COPYING-LGPL-2.1; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * You should have received a copy of the MPL along with this library + * in the file COPYING-MPL-1.1 + * + * The contents of this file are subject to the Mozilla Public License + * Version 1.1 (the "License"); you may not use this file except in + * compliance with the License. You may obtain a copy of the License at + * http://www.mozilla.org/MPL/ + * + * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY + * OF ANY KIND, either express or implied. See the LGPL or the MPL for + * the specific language governing rights and limitations. + */ + + +#include "nearest-point.h" + +namespace Geom +{ + +//////////////////////////////////////////////////////////////////////////////// +// D2 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 const& c, + D2 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 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 +all_nearest_points( Point const& p, + D2 const& c, + D2 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 result; + SBasis dd = dot(c - p, Geom::derivative(c)); + + std::vector zeros = Geom::roots(dd); + std::vector candidates; + candidates.push_back(from); + candidates.insert(candidates.end(), zeros.begin(), zeros.end()); + candidates.push_back(to); + std::vector 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 > versions + + +double nearest_point( Point const& p, + Piecewise< D2 > 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 +all_nearest_points( Point const& p, + Piecewise< D2 > 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 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 all_t; + std::vector< std::vector > all_np; + all_np.push_back( all_nearest_points(p, c[si], c.segT(from, si)) ); + std::vector 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 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 and Piecewise> + * + * + * Authors: + * + * Marco Cecchetti + * + * Copyright 2007-2008 authors + * + * This library is free software; you can redistribute it and/or + * modify it either under the terms of the GNU Lesser General Public + * License version 2.1 as published by the Free Software Foundation + * (the "LGPL") or, at your option, under the terms of the Mozilla + * Public License Version 1.1 (the "MPL"). If you do not alter this + * notice, a recipient may use your version of this file under either + * the MPL or the LGPL. + * + * You should have received a copy of the LGPL along with this library + * in the file COPYING-LGPL-2.1; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * You should have received a copy of the MPL along with this library + * in the file COPYING-MPL-1.1 + * + * The contents of this file are subject to the Mozilla Public License + * Version 1.1 (the "License"); you may not use this file except in + * compliance with the License. You may obtain a copy of the License at + * http://www.mozilla.org/MPL/ + * + * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY + * OF ANY KIND, either express or implied. See the LGPL or the MPL for + * the specific language governing rights and limitations. + */ + + +#ifndef _NEAREST_POINT_H_ +#define _NEAREST_POINT_H_ + + +#include + +#include "d2.h" +#include "piecewise.h" +#include "exception.h" + + + +namespace Geom +{ + +//////////////////////////////////////////////////////////////////////////////// +// D2 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 const& c, D2 const& dc, + double from = 0, double to = 1 ); + +inline +double nearest_point( Point const& p, + D2 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 +all_nearest_points( Point const& p, + D2 const& c, D2 const& dc, + double from = 0, double to = 1 ); + +inline +std::vector +all_nearest_points( Point const& p, + D2 const& c, + double from = 0, double to = 1 ) +{ + return all_nearest_points(p, c, Geom::derivative(c), from, to); +} + + +//////////////////////////////////////////////////////////////////////////////// +// Piecewise< D2 > versions + +double nearest_point( Point const& p, + Piecewise< D2 > const& c, + double from, double to ); + +inline +double nearest_point( Point const& p, Piecewise< D2 > const& c ) +{ + return nearest_point(p, c, c.cuts[0], c.cuts[c.size()]); +} + + +std::vector +all_nearest_points( Point const& p, + Piecewise< D2 > const& c, + double from, double to ); + +inline +std::vector +all_nearest_points( Point const& p, Piecewise< D2 > 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 + +#include + +#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 +#include +#include + +#include + + + +namespace Geom { namespace NL { + +class Matrix; +void swap(Matrix & m1, Matrix & m2); + + +class Matrix +{ +public: + // the matrix is not inizialized + Matrix(size_t n1, size_t n2) + : m_rows(n1), m_columns(n2) + { + m_matrix = gsl_matrix_alloc(n1, n2); + } + + Matrix(size_t n1, size_t n2, double x) + :m_rows(n1), m_columns(n2) + { + m_matrix = gsl_matrix_alloc(n1, n2); + gsl_matrix_set_all(m_matrix, x); + } + + Matrix( Matrix const& _matrix ) + : m_rows(_matrix.rows()), m_columns(_matrix.columns()) + { + m_matrix = gsl_matrix_alloc(rows(), columns()); + gsl_matrix_memcpy(m_matrix, _matrix.m_matrix); + } + +// explicit +// Matrix( gsl_matrix* m, size_t n1, size_t n2) +// : m_rows(n1), m_columns(n2), m_matrix(m) +// { +// } + + Matrix & operator=(Matrix const& _matrix) + { + assert( rows() == _matrix.rows() && columns() == _matrix.columns() ); + gsl_matrix_memcpy(m_matrix, _matrix.m_matrix); + + return *this; + } + + virtual ~Matrix() + { + gsl_matrix_free(m_matrix); + } + + void set_all( double x = 0) + { + gsl_matrix_set_all(m_matrix, x); + } + + void set_identity() + { + gsl_matrix_set_identity(m_matrix); + } + + const double & operator() (size_t i, size_t j) const + { + return *gsl_matrix_const_ptr(m_matrix, i, j); + } + + double & operator() (size_t i, size_t j) + { + return *gsl_matrix_ptr(m_matrix, i, j); + } + + gsl_matrix* get_gsl_matrix() + { + return m_matrix; + } + + void swap_rows(size_t i, size_t j) + { + gsl_matrix_swap_rows(m_matrix, i, j); + } + + void swap_columns(size_t i, size_t j) + { + gsl_matrix_swap_columns(m_matrix, i, j); + } + + Matrix & transpose() + { + gsl_matrix_transpose(m_matrix); + return (*this); + } + + Matrix & scale(double x) + { + gsl_matrix_scale(m_matrix, x); + return (*this); + } + + Matrix & translate(double x) + { + gsl_matrix_add_constant(m_matrix, x); + return (*this); + } + + Matrix & operator+=(Matrix const& _matrix) + { + gsl_matrix_add(m_matrix, _matrix.m_matrix); + return (*this); + } + + Matrix & operator-=(Matrix const& _matrix) + { + gsl_matrix_sub(m_matrix, _matrix.m_matrix); + return (*this); + } + + bool is_zero() const + { + return gsl_matrix_isnull(m_matrix); + } + + bool is_positive() const + { + return gsl_matrix_ispos(m_matrix); + } + + bool is_negative() const + { + return gsl_matrix_isneg(m_matrix); + } + + bool is_non_negative() const + { + for ( unsigned int i = 0; i < rows(); ++i ) + { + for ( unsigned int j = 0; j < columns(); ++j ) + { + if ( (*this)(i,j) < 0 ) return false; + } + } + return true; + } + + double max() const + { + return gsl_matrix_max(m_matrix); + } + + double min() const + { + return gsl_matrix_min(m_matrix); + } + + std::pair + max_index() const + { + std::pair indices; + gsl_matrix_max_index(m_matrix, &(indices.first), &(indices.second)); + return indices; + } + + std::pair + min_index() const + { + std::pair indices; + gsl_matrix_min_index(m_matrix, &(indices.first), &(indices.second)); + return indices; + } + + friend + void swap(Matrix & m1, Matrix & m2); + + size_t rows() const + { + return m_rows; + } + + size_t columns() const + { + return m_columns; + } + +private: + size_t m_rows, m_columns; + gsl_matrix* m_matrix; +}; + +void swap(Matrix & m1, Matrix & m2) +{ + assert( m1.rows() == m2.rows() && m1.columns() == m2.columns() ); + std::swap(m1.m_matrix, m2.m_matrix); +} + + +} } // end namespaces + +#endif /*_NL_MATRIX_H_*/ 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 +#include + +#include + + +namespace Geom { namespace NL { + +class Vector; +void swap(Vector & v1, Vector & v2); + + +class Vector +{ +public: + Vector(size_t n) + : m_size(n) + { + m_vector = gsl_vector_alloc(n); + } + + Vector(size_t n, double x) + : m_size(n) + { + m_vector = gsl_vector_alloc(n); + gsl_vector_set_all(m_vector, x); + } + + // create a vector with n elements all set to zero + // but the i-th that is set to 1 + Vector(size_t n, size_t i) + : m_size(n) + { + m_vector = gsl_vector_alloc(n); + gsl_vector_set_basis(m_vector, i); + } + + Vector(Vector const& _vector) + : m_size(_vector.size()) + { + m_vector = gsl_vector_alloc(size()); + gsl_vector_memcpy(m_vector, _vector.m_vector); + } + + virtual ~Vector() + { + gsl_vector_free(m_vector); + } + + Vector & operator=(Vector const& _vector) + { + assert( size() == _vector.size() ); + gsl_vector_memcpy(m_vector, _vector.m_vector); + } + + 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 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 +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(sif); + unsigned int ei = static_cast(eif); + if ( si == sz ) + { + --si; + st = 1; + } + if ( ei == sz ) + { + --ei; + et = 1; + } + if ( si == ei ) + { + std::vector 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 all_t; + std::vector< std::vector > all_np; + all_np.push_back( _path[si].allNearestPoints(_point, st) ); + std::vector 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 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(sif); + unsigned int ei = static_cast(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 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 + 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 inner; + public: explicit SBasisCurve(D2 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 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 + 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 class BezierCurve : public Curve { + private: D2 inner; + public: template static void assert_degree(BezierCurve 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 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 roots(double v, Dim2 d) const + std::vector roots(double v, Dim2 d) const; + + std::vector + allNearestPoints( Point const& p, double from = 0, double to = 1 ) const; + + double nearestPoint( Point const& p, double from = 0, double to = 1 ) const { - 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 pointAndDerivatives(Coord t, unsigned n) const - { - return SBasisCurve(toSBasis()).pointAndDerivatives(t, n); - } + std::vector pointAndDerivatives(Coord t, unsigned int n) const; D2 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 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,11 +529,11 @@ 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; @@ -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(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(k); + if ( i == sz ) + { + --i; + lt = 1; + } + return (*this)[i].valueAt(lt, d); } std::vector roots(double v, Dim2 d) const { @@ -682,7 +764,28 @@ public: } return res; } - + + std::vector + allNearestPoints(Point const& _point, double from, double to) const; + + std::vector + 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 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 intersect(Rect const & a, Rect const & b) { return x && y ? boost::optional(Rect(*x, *y)) : boost::optional(); } +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 > sbasis_to_bezier(D2 const &B) { // mutating void subpath_from_sbasis(Geom::OldPathSetBuilder &pb, D2 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 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 B, doubl void build_from_sbasis(Geom::PathBuilder &pb, D2 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 SBasis::valueAndDerivatives(double t, unsigned n) const { + std::vector 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 valueAndDerivatives(double t, unsigned n) const { - std::vector 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 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 - * - * 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 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 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 + * + * 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 + +#include + + + + +namespace Geom +{ + + +Rect SVGEllipticalArc::boundsExact() const +{ + std::vector extremes(4); + double cosrot = std::cos(rotation_angle()); + double sinrot = std::sin(rotation_angle()); + extremes[0] = std::atan2( -ray(Y) * sinrot, ray(X) * cosrot ); + extremes[1] = extremes[0] + M_PI; + if ( extremes[0] < 0 ) extremes[0] += 2*M_PI; + extremes[2] = std::atan2( ray(Y) * cosrot, ray(X) * sinrot ); + extremes[3] = extremes[2] + M_PI; + if ( extremes[2] < 0 ) extremes[2] += 2*M_PI; + + + std::vectorarc_extremes(4); + arc_extremes[0] = initialPoint()[X]; + arc_extremes[1] = finalPoint()[X]; + if ( arc_extremes[0] < arc_extremes[1] ) + std::swap(arc_extremes[0], arc_extremes[1]); + arc_extremes[2] = initialPoint()[Y]; + arc_extremes[3] = finalPoint()[Y]; + if ( arc_extremes[2] < arc_extremes[3] ) + std::swap(arc_extremes[2], arc_extremes[3]); + + + if ( start_angle() < end_angle() ) + { + if ( sweep_flag() ) + { + for ( unsigned int i = 0; i < extremes.size(); ++i ) + { + if ( start_angle() < extremes[i] && extremes[i] < end_angle() ) + { + arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1]; + } + } + } + else + { + for ( unsigned int i = 0; i < extremes.size(); ++i ) + { + if ( start_angle() > extremes[i] || extremes[i] > end_angle() ) + { + arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1]; + } + } + } + } + else + { + if ( sweep_flag() ) + { + for ( unsigned int i = 0; i < extremes.size(); ++i ) + { + if ( start_angle() < extremes[i] || extremes[i] < end_angle() ) + { + arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1]; + } + } + } + else + { + for ( unsigned int i = 0; i < extremes.size(); ++i ) + { + if ( start_angle() > extremes[i] && extremes[i] > end_angle() ) + { + arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1]; + } + } + } + } + + return Rect( Point(arc_extremes[1], arc_extremes[3]) , + Point(arc_extremes[0], arc_extremes[2]) ); + +} + + +std::vector +SVGEllipticalArc::roots(double v, Dim2 d) const +{ + if ( d > Y ) + { + THROW_RANGEERROR("dimention out of range"); + } + + std::vector 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 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 +SVGEllipticalArc::pointAndDerivatives(Coord t, unsigned int n) const +{ + std::vector 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 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 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 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 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 solX = roots(np[X],X); +// std::vector solY = roots(np[Y],Y); +// double t; +// if ( are_near(solX[0], solY[0]) || are_near(solX[0], solY[1])) +// { +// t = solX[0]; +// } +// else +// { +// t = solX[1]; +// } +// if ( !(t < from || t > to) ) +// { +// result.push_back(t); +// } +// else +// { +// +// } + } + + // solve the equation == 0 + // that provides min and max distance points + // on the ellipse E wrt the point p + // after the substitutions: + // cos(t) = (1 - s^2) / (1 + s^2) + // sin(t) = 2t / (1 + s^2) + // where s = tan(t/2) + // we get a 4th degree equation in s + /* + * ry s^4 ((-cy + py) Cos[Phi] + (cx - px) Sin[Phi]) + + * ry ((cy - py) Cos[Phi] + (-cx + px) Sin[Phi]) + + * 2 s^3 (rx^2 - ry^2 + (-cx + px) rx Cos[Phi] + (-cy + py) rx Sin[Phi]) + + * 2 s (-rx^2 + ry^2 + (-cx + px) rx Cos[Phi] + (-cy + py) rx Sin[Phi]) + */ + + Point p_c = p - center(); + double rx2_ry2 = (ray(X) - ray(Y)) * (ray(X) + ray(Y)); + double cosrot = std::cos( rotation_angle() ); + double sinrot = std::sin( rotation_angle() ); + double expr1 = ray(X) * (p_c[X] * cosrot + p_c[Y] * sinrot); + 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 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 -> 0 iff coeff[4] == 0 + // so we add M_PI to the solutions being lim arctan(s) = PI when s->Infinity + if ( (real_sol.size() % 2) != 0 ) + { + real_sol.push_back(M_PI); + } + + double mindistsq1 = std::numeric_limits::max(); + double mindistsq2 = std::numeric_limits::max(); + double dsq; + unsigned int mi1, mi2; + for ( unsigned int i = 0; i < real_sol.size(); ++i ) + { + dsq = distanceSq(p, pointAtAngle(real_sol[i])); + if ( mindistsq1 > dsq ) + { + mindistsq2 = mindistsq1; + mi2 = mi1; + mindistsq1 = dsq; + mi1 = i; + } + else if ( mindistsq2 > dsq ) + { + mindistsq2 = dsq; + mi2 = i; + } + } + + double t = map_to_01( real_sol[mi1] ); + if ( !(t < from || t > to) ) + { + result.push_back(t); + } + + bool second_sol = false; + t = map_to_01( real_sol[mi2] ); + if ( real_sol.size() == 4 && !(t < from || t > to) ) + { + if ( result.empty() || are_near(mindistsq1, mindistsq2) ) + { + result.push_back(t); + second_sol = true; + } + } + + // we need to test extreme points too + double dsq1 = distanceSq(p, pointAt(from)); + double dsq2 = distanceSq(p, pointAt(to)); + if ( second_sol ) + { + if ( mindistsq2 > dsq1 ) + { + result.clear(); + result.push_back(from); + mindistsq2 = dsq1; + } + else if ( are_near(mindistsq2, dsq) ) + { + result.push_back(from); + } + if ( mindistsq2 > dsq2 ) + { + result.clear(); + result.push_back(to); + } + else if ( are_near(mindistsq2, dsq2) ) + { + result.push_back(to); + } + + } + else + { + if ( result.empty() ) + { + if ( are_near(dsq1, dsq2) ) + { + result.push_back(from); + result.push_back(to); + } + else if ( dsq2 > dsq1 ) + { + result.push_back(from); + } + else + { + result.push_back(to); + } + } + } + + return result; +} + + +} // 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 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 */