diff --git a/src/2geom/CMakeLists.txt b/src/2geom/CMakeLists.txt index 62c78dfc..d5e3fa14 100755 --- a/src/2geom/CMakeLists.txt +++ b/src/2geom/CMakeLists.txt @@ -92,6 +92,9 @@ math-utils.h nearest-time.cpp nearest-time.h +furthest-time.cpp +furthest-time.h + numeric/matrix.cpp ord.h diff --git a/src/2geom/bezier-curve.cpp b/src/2geom/bezier-curve.cpp index 86626304..8088463b 100644 --- a/src/2geom/bezier-curve.cpp +++ b/src/2geom/bezier-curve.cpp @@ -35,6 +35,7 @@ #include <2geom/path-sink.h> #include <2geom/basic-intersection.h> #include <2geom/nearest-time.h> +#include <2geom/furthest-time.h> namespace Geom { @@ -222,6 +223,11 @@ Coord BezierCurve::nearestTime(Point const &p, Coord from, Coord to) const return nearest_time(p, inner, from, to); } +Coord BezierCurve::furthestTime(Point const &p, Coord from, Coord to) const +{ + return furthest_time(p, inner, from, to); +} + void BezierCurve::feed(PathSink &sink, bool moveto_initial) const { if (size() > 4) { @@ -294,6 +300,18 @@ Coord BezierCurveN<1>::nearestTime(Point const& p, Coord from, Coord to) const else return from + t*(to-from); } +template<> +Coord BezierCurveN<1>::furthestTime(Point const& p, Coord from, Coord to) const +{ + using std::swap; + + if ( from > to ) swap(from, to); + Point ip = pointAt(from); + Point fp = pointAt(to); + if ( Geom::distance(fp,p) >= Geom::distance(ip,p) ) return to; + else return from; +} + template <> std::vector BezierCurveN<1>::intersect(Curve const &other, Coord eps) const { diff --git a/src/2geom/bezier-curve.h b/src/2geom/bezier-curve.h index 9416ba78..2d84ae53 100644 --- a/src/2geom/bezier-curve.h +++ b/src/2geom/bezier-curve.h @@ -158,6 +158,8 @@ class BezierCurve : public Curve { return (inner[d] - v).roots(); } virtual Coord nearestTime(Point const &p, Coord from = 0, Coord to = 1) const; + virtual Coord furthestTime(Point const& p, Coord from = 0, Coord to = 1) const; + virtual Coord length(Coord tolerance) const; virtual std::vector intersect(Curve const &other, Coord eps = EPSILON) const; virtual Point pointAt(Coord t) const { return inner.pointAt(t); } @@ -275,6 +277,9 @@ class BezierCurveN virtual Coord nearestTime(Point const &p, Coord from = 0, Coord to = 1) const { return BezierCurve::nearestTime(p, from, to); } + virtual Coord furthestTime(Point const &p, Coord from = 0, Coord to = 1) const { + return BezierCurve::furthestTime(p, from, to); + } virtual std::vector intersect(Curve const &other, Coord eps = EPSILON) const { // call super. this is implemented only to allow specializations return BezierCurve::intersect(other, eps); @@ -318,6 +323,7 @@ template <> inline bool BezierCurveN<1>::isDegenerate() const { template <> inline bool BezierCurveN<1>::isLineSegment() const { return true; } template <> Curve *BezierCurveN<1>::derivative() const; template <> Coord BezierCurveN<1>::nearestTime(Point const &, Coord, Coord) const; +template <> Coord BezierCurveN<1>::furthestTime(Point const& p, Coord from, Coord to) const; template <> std::vector BezierCurveN<1>::intersect(Curve const &, Coord) const; template <> int BezierCurveN<1>::winding(Point const &) const; template <> void BezierCurveN<1>::feed(PathSink &sink, bool moveto_initial) const; diff --git a/src/2geom/circle.cpp b/src/2geom/circle.cpp index 934a8d3a..1592f1c0 100644 --- a/src/2geom/circle.cpp +++ b/src/2geom/circle.cpp @@ -126,6 +126,12 @@ Coord Circle::nearestTime(Point const &p) const { return timeAt(p); } +Coord Circle::furthestTime(Point const &p) const { + Coord angle = std::fmod(timeAt(p) + rad_from_deg(180), 2*M_PI); + if (angle < 0) angle += 2*M_PI; + return angle; +} + bool Circle::contains(Rect const &r) const { for (unsigned i = 0; i < 4; ++i) { diff --git a/src/2geom/circle.h b/src/2geom/circle.h index a4d5f209..33899d1b 100644 --- a/src/2geom/circle.h +++ b/src/2geom/circle.h @@ -90,6 +90,7 @@ class Circle Coord valueAt(Coord t, Dim2 d) const; Coord timeAt(Point const &p) const; Coord nearestTime(Point const &p) const; + Coord furthestTime(Point const &p) const; bool contains(Point const &p) const { return distance(p, _center) <= _radius; } bool contains(Rect const &other) const; diff --git a/src/2geom/conicsec.cpp b/src/2geom/conicsec.cpp index 2d396ba3..b17c51c0 100644 --- a/src/2geom/conicsec.cpp +++ b/src/2geom/conicsec.cpp @@ -1470,6 +1470,14 @@ std::vector xAx::allNearestTimes (const Point &P) const return points; } +/* + * Maybe we also need do this?. + * std::vector xAx::allFurthestTimes (const Point &P) const + * + * P: the point to compute the furthest one + */ + + bool clip (std::vector & rq, const xAx & cs, const Rect & R) diff --git a/src/2geom/curve.cpp b/src/2geom/curve.cpp index 387a9180..02c5552a 100644 --- a/src/2geom/curve.cpp +++ b/src/2geom/curve.cpp @@ -34,6 +34,7 @@ #include <2geom/curve.h> #include <2geom/exception.h> #include <2geom/nearest-time.h> +#include <2geom/furthest-time.h> #include <2geom/sbasis-geometric.h> #include <2geom/sbasis-to-bezier.h> #include <2geom/ord.h> @@ -49,11 +50,23 @@ Coord Curve::nearestTime(Point const& p, Coord a, Coord b) const return nearest_time(p, toSBasis(), a, b); } +Coord Curve::furthestTime(Point const& p, Coord a, Coord b) const +{ + return furthest_time(p, toSBasis(), a, b); +} + std::vector Curve::allNearestTimes(Point const& p, Coord from, Coord to) const { return all_nearest_times(p, toSBasis(), from, to); } +std::vector Curve::allFurthestTimes(Point const& p, Coord from, Coord to) const +{ + return all_furthest_times(p, toSBasis(), from, to); +} + + + Coord Curve::length(Coord tolerance) const { return ::Geom::length(toSBasis(), tolerance); diff --git a/src/2geom/curve.h b/src/2geom/curve.h index 41db9ca8..7198ef4a 100644 --- a/src/2geom/curve.h +++ b/src/2geom/curve.h @@ -265,6 +265,33 @@ class Curve return allNearestTimes(p, i.min(), i.max()); } + /** @brief Compute a time value at which the curve comes furthest to a specified point. + * The first value with the smallest distance is returned if there are multiple such points. + * @param p Query point + * @param a Minimum time value to consider + * @param b Maximum time value to consider; \f$a < b\f$ + * @return \f$q \in [a, b]: ||\mathbf{C}(q) - \mathbf{p}|| = + \inf(\{r \in \mathbb{R} : ||\mathbf{C}(r) - \mathbf{p}||\})\f$ */ + virtual Coord furthestTime( Point const& p, Coord a = 0, Coord b = 1 ) const; + + /** @brief A version that takes an Interval. */ + Coord furthestTime(Point const &p, Interval const &i) const { + return furthestTime(p, i.min(), i.max()); + } + + /** @brief Compute a time value at which the curve comes furthest to a specified point. + * @param p Query point + * @param a Minimum time value to consider + * @param b Maximum time value to consider; \f$a < b\f$ + * @return Vector of points closest and equally far away from the query point */ + virtual std::vector allFurthestTimes( Point const& p, Coord from = 0, + Coord to = 1 ) const; + + /** @brief A version that takes an Interval. */ + std::vector allFurthestTimes(Point const &p, Interval const &i) { + return allFurthestTimes(p, i.min(), i.max()); + } + /** @brief Compute the arc length of this curve. * For a curve \f$\mathbf{C}(t) = (C_x(t), C_y(t))\f$, arc length is defined for 2D curves as * \f[ \ell = \int_{0}^{1} \sqrt { [C_x'(t)]^2 + [C_y'(t)]^2 }\, \text{d}t \f] @@ -346,6 +373,11 @@ Coord nearest_time(Point const& p, Curve const& c) { return c.nearestTime(p); } +inline +Coord furthest_time(Point const& p, Curve const& c) { + return c.furthestTime(p); +} + // for make benefit glorious library of Boost Pointer Container inline Curve *new_clone(Curve const &c) { diff --git a/src/2geom/elliptical-arc.cpp b/src/2geom/elliptical-arc.cpp index ec62b4be..c30f1e71 100644 --- a/src/2geom/elliptical-arc.cpp +++ b/src/2geom/elliptical-arc.cpp @@ -561,8 +561,234 @@ std::vector EllipticalArc::allNearestTimes( Point const& p, double from, return result; } -#endif +//TODO: Fix it currently give nearest points +std::vector EllipticalArc::allFurthestTimes( Point const& p, double from, double to ) const +{ + std::vector result; + + if ( from > to ) std::swap(from, to); + if ( from < 0 || to > 1 ) + { + THROW_RANGEERROR("[from,to] interval out of range"); + } + + if ( ( are_near(ray(X), 0) && are_near(ray(Y), 0) ) || are_near(from, to) ) + { + result.push_back(to); + return result; + } + else if ( are_near(ray(X), 0) || are_near(ray(Y), 0) ) + { + LineSegment seg(pointAt(from), pointAt(to)); + Point fp = seg.pointAt( seg.furthestTime(p) ); + if ( are_near(ray(Y), 0) ) + { + if ( are_near(rotationAngle(), M_PI/2) + || are_near(rotationAngle(), 3*M_PI/2) ) + { + result = roots(fp[Y], Y); + } + else + { + result = roots(fp[X], X); + } + } + else + { + if ( are_near(rotationAngle(), M_PI/2) + || are_near(rotationAngle(), 3*M_PI/2) ) + { + result = roots(fp[X], X); + } + else + { + result = roots(fp[Y], Y); + } + } + return result; + } + else if ( are_near(ray(X), ray(Y)) ) + { + Point r = p - center(); + if ( are_near(r, Point(0,0)) ) + { + THROW_INFINITESOLUTIONS(0); + } + // TODO: implement case r != 0 +// Point np = ray(X) * unit_vector(r); +// std::vector solX = roots(np[X],X); +// std::vector solY = roots(np[Y],Y); +// double t; +// if ( are_near(solX[0], solY[0]) || are_near(solX[0], solY[1])) +// { +// t = solX[0]; +// } +// else +// { +// t = solX[1]; +// } +// if ( !(t < from || t > to) ) +// { +// result.push_back(t); +// } +// else +// { +// +// } + } + + // solve the equation == 0 + // that provides min and max distance points + // on the ellipse E wrt the point p + // after the substitutions: + // cos(t) = (1 - s^2) / (1 + s^2) + // sin(t) = 2t / (1 + s^2) + // where s = tan(t/2) + // we get a 4th degree equation in s + /* + * ry s^4 ((-cy + py) Cos[Phi] + (cx - px) Sin[Phi]) + + * ry ((cy - py) Cos[Phi] + (-cx + px) Sin[Phi]) + + * 2 s^3 (rx^2 - ry^2 + (-cx + px) rx Cos[Phi] + (-cy + py) rx Sin[Phi]) + + * 2 s (-rx^2 + ry^2 + (-cx + px) rx Cos[Phi] + (-cy + py) rx Sin[Phi]) + */ + + Point p_c = p - center(); + double rx2_ry2 = (ray(X) - ray(Y)) * (ray(X) + ray(Y)); + double sinrot, cosrot; + sincos(rotationAngle(), sinrot, cosrot); + double expr1 = ray(X) * (p_c[X] * cosrot + p_c[Y] * sinrot); + Poly coeff; + coeff.resize(5); + coeff[4] = ray(Y) * ( p_c[Y] * cosrot - p_c[X] * sinrot ); + coeff[3] = 2 * ( rx2_ry2 + expr1 ); + coeff[2] = 0; + coeff[1] = 2 * ( -rx2_ry2 + expr1 ); + coeff[0] = -coeff[4]; +// for ( unsigned int i = 0; i < 5; ++i ) +// std::cerr << "c[" << i << "] = " << coeff[i] << std::endl; + + std::vector real_sol; + // gsl_poly_complex_solve raises an error + // if the leading coefficient is zero + if ( are_near(coeff[4], 0) ) + { + real_sol.push_back(0); + if ( !are_near(coeff[3], 0) ) + { + double sq = -coeff[1] / coeff[3]; + if ( sq > 0 ) + { + double s = std::sqrt(sq); + real_sol.push_back(s); + real_sol.push_back(-s); + } + } + } + else + { + real_sol = solve_reals(coeff); + } + + for ( unsigned int i = 0; i < real_sol.size(); ++i ) + { + real_sol[i] = 2 * std::atan(real_sol[i]); + if ( real_sol[i] < 0 ) real_sol[i] += 2*M_PI; + } + // when s -> Infinity then -> 0 iff coeff[4] == 0 + // so we add M_PI to the solutions being lim arctan(s) = PI when s->Infinity + if ( (real_sol.size() % 2) != 0 ) + { + real_sol.push_back(M_PI); + } + + double maxdistsq1 = 0; + double maxdistsq2 = 0; + double dsq = 0; + unsigned int mi1 = 0, mi2 = 0; + for ( unsigned int i = 0; i < real_sol.size(); ++i ) + { + dsq = distanceSq(p, pointAtAngle(real_sol[i])); + if ( maxdistsq1 < dsq ) + { + maxdistsq2 = maxdistsq1; + mi2 = mi1; + maxdistsq1 = dsq; + mi1 = i; + } + else if ( maxdistsq2 < dsq ) + { + maxdistsq2 = dsq; + mi2 = i; + } + } + + double t = timeAtAngle(real_sol[mi1]); + if ( !(t < from || t > to) ) + { + result.push_back(t); + } + + bool second_sol = false; + t = timeAtAngle(real_sol[mi2]); + if ( real_sol.size() == 4 && !(t < from || t > to) ) + { + if ( result.empty() || are_near(maxdistsq1, maxdistsq2) ) + { + 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 ( maxdistsq2 < dsq1 ) + { + result.clear(); + result.push_back(from); + maxdistsq2 = dsq1; + } + else if ( are_near(maxdistsq2, dsq1) ) + { + result.push_back(from); + } + if ( maxdistsq2 < dsq2 ) + { + result.clear(); + result.push_back(to); + } + else if ( are_near(maxdistsq2, 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(to); + } + else + { + result.push_back(from); + } + } + } + + return result; +} +#endif void EllipticalArc::_filterIntersections(std::vector &xs, bool is_first) const { @@ -570,13 +796,13 @@ void EllipticalArc::_filterIntersections(std::vector &xs, boo std::vector::reverse_iterator i = xs.rbegin(), last = xs.rend(); while (i != last) { Coord &t = is_first ? i->first : i->second; - assert(are_near(_ellipse.pointAt(t), i->point(), 1e-5)); + assert(are_near(_ellipse.pointAt(t), i->point(), 1e-6)); t = timeAtAngle(t); if (!unit.contains(t)) { xs.erase((++i).base()); continue; } else { - assert(are_near(pointAt(t), i->point(), 1e-5)); + assert(are_near(pointAt(t), i->point(), 1e-6)); ++i; } } diff --git a/src/2geom/elliptical-arc.h b/src/2geom/elliptical-arc.h index 96a92eb8..48436ac9 100644 --- a/src/2geom/elliptical-arc.h +++ b/src/2geom/elliptical-arc.h @@ -283,6 +283,10 @@ class EllipticalArc : public Curve } return allNearestTimes(p, from, to).front(); } + virtual std::vector allFurthestTimes( Point const& p, double from = 0, double to = 1 ) const; + virtual double furthestTime( Point const& p, double from = 0, double to = 1 ) const { + return allFurthestTimes(p, from, to).front(); + } #endif virtual std::vector intersect(Curve const &other, Coord eps=EPSILON) const; virtual int degreesOfFreedom() const { return 7; } diff --git a/src/2geom/furthest-time.cpp b/src/2geom/furthest-time.cpp new file mode 100644 index 00000000..b39c542f --- /dev/null +++ b/src/2geom/furthest-time.cpp @@ -0,0 +1,323 @@ +/** @file + * @brief Furthest time routines for D2 and Piecewise> + *//* + * Authors: + * Marco Cecchetti + * Jabier Arraiza + * + * 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 <2geom/furthest-time.h> +#include + +namespace Geom +{ + +Coord furthest_time(Point const &p, D2 const &input, Coord from, Coord to) +{ + Interval domain(from, to); + bool partial = false; + + if (domain.min() < 0 || domain.max() > 1) { + THROW_RANGEERROR("[from,to] interval out of bounds"); + } + + if (input.isConstant(0)) return from; + + D2 bez; + if (domain.min() != 0 || domain.max() != 1) { + bez = portion(input, domain) - p; + partial = true; + } else { + bez = input - p; + } + + // find extrema of the function x(t)^2 + y(t)^2 + // use the fact that (f^2)' = 2 f f' + // this reduces the order of the distance function by 1 + D2 deriv = derivative(bez); + std::vector ts = (multiply(bez[X], deriv[X]) + multiply(bez[Y], deriv[Y])).roots(); + + Coord t = -1, maxd = 0; + for (unsigned i = 0; i < ts.size(); ++i) { + Coord droot = L2sq(bez.valueAt(ts[i])); + if (droot > maxd) { + maxd = droot; + t = ts[i]; + } + } + + // also check endpoints + Coord dinitial = L2sq(bez.at0()); + Coord dfinal = L2sq(bez.at1()); + + if (dinitial > maxd) { + maxd = dinitial; + t = 0; + } + if (dfinal > maxd) { + //maxd = dfinal; + t = 1; + } + + if (partial) { + t = domain.valueAt(t); + } + return t; +} + +//////////////////////////////////////////////////////////////////////////////// +// D2 versions + +/* + * Return the parameter t of the furthest time value 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 furthest time value to "p" that is found. + */ + +double furthest_time(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"); + } + if (c.isConstant()) return from; + SBasis dd = dot(c - p, dc); + //std::cout << dd << std::endl; + std::vector zeros = Geom::roots(dd); + + double closest = from; + double max_dist_sq = L2sq(c(from) - p); + for ( size_t i = 0; i < zeros.size(); ++i ) + { + double distsq = L2sq(c(zeros[i]) - p); + if ( max_dist_sq < L2sq(c(zeros[i]) - p) ) + { + closest = zeros[i]; + max_dist_sq = distsq; + } + } + if ( max_dist_sq < L2sq( c(to) - p ) ) + closest = to; + return closest; + +} + +/* + * Return the parameters t of all the furthest 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_furthest_times(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"); + } + + std::vector result; + if (c.isConstant()) { + result.push_back(from); + return result; + } + SBasis dd = dot(c - p, dc); + + 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 i = 0; i < candidates.size(); ++i) { + distsq.push_back(L2sq(c(candidates[i]) - p)); + } + unsigned closest = 0; + double dsq = distsq[0]; + for (unsigned i = 1; i < candidates.size(); ++i) { + if (dsq < distsq[i]) { + closest = i; + dsq = distsq[i]; + } + } + for (unsigned i = 0; i < candidates.size(); ++i) { + if (distsq[closest] == distsq[i]) { + result.push_back(candidates[i]); + } + } + return result; +} + + +//////////////////////////////////////////////////////////////////////////////// +// Piecewise< D2 > versions + + +double furthest_time(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 si = c.segN(from); + unsigned ei = c.segN(to); + if (si == ei) { + double furthest = + furthest_time(p, c[si], c.segT(from, si), c.segT(to, si)); + return c.mapToDomain(furthest, si); + } + + double t; + double furthest = furthest_time(p, c[si], c.segT(from, si)); + unsigned int ni = si; + double dsq; + double maxdistsq = distanceSq(p, c[si](furthest)); + Rect bb; + for (unsigned i = si + 1; i < ei; ++i) { + bb = *bounds_fast(c[i]); + dsq = distanceSq(p, bb); + if ( maxdistsq > dsq ) continue; + + t = furthest_time(p, c[i]); + dsq = distanceSq(p, c[i](t)); + if (maxdistsq < dsq) { + furthest = t; + ni = i; + maxdistsq = dsq; + } + } + bb = *bounds_fast(c[ei]); + dsq = distanceSq(p, bb); + if (maxdistsq < dsq) { + t = furthest_time(p, c[ei], 0, c.segT(to, ei)); + dsq = distanceSq(p, c[ei](t)); + if (maxdistsq < dsq) { + furthest = t; + ni = ei; + } + } + return c.mapToDomain(furthest, ni); +} + +std::vector +all_furthest_times(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 si = c.segN(from); + unsigned ei = c.segN(to); + if ( si == ei ) + { + std::vector all_furthest = + all_furthest_times(p, c[si], c.segT(from, si), c.segT(to, si)); + for ( unsigned int i = 0; i < all_furthest.size(); ++i ) + { + all_furthest[i] = c.mapToDomain(all_furthest[i], si); + } + return all_furthest; + } + std::vector all_t; + std::vector< std::vector > all_np; + all_np.push_back( all_furthest_times(p, c[si], c.segT(from, si)) ); + std::vector ni; + ni.push_back(si); + double dsq; + double maxdistsq = distanceSq( p, c[si](all_np.front().front()) ); + Rect bb; + + for (unsigned i = si + 1; i < ei; ++i) { + bb = *bounds_fast(c[i]); + dsq = distanceSq(p, bb); + if ( maxdistsq > dsq ) continue; + all_t = all_furthest_times(p, c[i]); + dsq = distanceSq( p, c[i](all_t.front()) ); + if ( maxdistsq < dsq ) + { + all_np.clear(); + all_np.push_back(all_t); + ni.clear(); + ni.push_back(i); + maxdistsq = dsq; + } + else if ( maxdistsq == dsq ) + { + all_np.push_back(all_t); + ni.push_back(i); + } + } + bb = *bounds_fast(c[ei]); + dsq = distanceSq(p, bb); + if (maxdistsq < dsq) { + all_t = all_furthest_times(p, c[ei], 0, c.segT(to, ei)); + dsq = distanceSq( p, c[ei](all_t.front()) ); + if (maxdistsq < 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 (maxdistsq == dsq) { + all_np.push_back(all_t); + ni.push_back(ei); + } + } + std::vector all_furthest; + for (unsigned i = 0; i < all_np.size(); ++i) { + for (unsigned int j = 0; j < all_np[i].size(); ++j) { + all_furthest.push_back( c.mapToDomain(all_np[i][j], ni[i]) ); + } + } + all_furthest.erase(std::unique(all_furthest.begin(), all_furthest.end()), + all_furthest.end()); + return all_furthest; +} + +} // end namespace Geom + + diff --git a/src/2geom/furthest-time.h b/src/2geom/furthest-time.h new file mode 100644 index 00000000..e17a2797 --- /dev/null +++ b/src/2geom/furthest-time.h @@ -0,0 +1,141 @@ +/** @file + * @brief Furthest time routines for D2 and Piecewise> + *//* + * Authors: + * Marco Cecchetti + * Jabier Arraiza + * + * 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 LIB2GEOM_SEEN_FURTHEST_TIME_H +#define LIB2GEOM_SEEN_FURTHEST_TIME_H + + +#include + +#include <2geom/d2.h> +#include <2geom/piecewise.h> +#include <2geom/exception.h> +#include <2geom/bezier.h> + + +namespace Geom +{ + +/* + * Given a line L specified by a point A and direction vector v, + * return the point on L furthest to p. Note that the returned value + * is with respect to the _normalized_ direction of v! + */ +inline double furthest_time(Point const, Point const, Point const) +{ + return infinity(); +} + +Coord furthest_time(Point const &p, D2 const &bez, Coord from = 0, Coord to = 1); + +//////////////////////////////////////////////////////////////////////////////// +// D2 versions + +/* + * Return the parameter t of a furthest point on the portion of the curve "c", + * related to the interval [from, to], to the point "p". + * The needed curve derivative "deriv" is passed as parameter. + * The function return the first furthest point to "p" that is found. + */ +double furthest_time(Point const &p, + D2 const &c, D2 const &deriv, + double from = 0, double to = 1); + +inline +double furthest_time(Point const &p, + D2 const &c, + double from = 0, double to = 1 ) +{ + return furthest_time(p, c, Geom::derivative(c), from, to); +} + +/* + * Return the parameters t of all the furthest times 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_furthest_times(Point const& p, + D2 const& c, D2 const& dc, + double from = 0, double to = 1 ); + +inline +std::vector +all_furthest_times(Point const &p, + D2 const &c, + double from = 0, double to = 1) +{ + return all_furthest_times(p, c, Geom::derivative(c), from, to); +} + + +//////////////////////////////////////////////////////////////////////////////// +// Piecewise< D2 > versions + +double furthest_time(Point const &p, + Piecewise< D2 > const &c, + double from, double to); + +inline +double furthest_time(Point const& p, Piecewise< D2 > const &c) +{ + return furthest_time(p, c, c.cuts[0], c.cuts[c.size()]); +} + + +std::vector +all_furthest_times(Point const &p, + Piecewise< D2 > const &c, + double from, double to); + +inline +std::vector +all_furthest_times( Point const& p, Piecewise< D2 > const& c ) +{ + return all_furthest_times(p, c, c.cuts[0], c.cuts[c.size()]); +} + +} // end namespace Geom + +#endif // LIB2GEOM_SEEN_FURTHEST_TIME_H +/* + Local Variables: + mode:c++ + c-file-style:"stroustrup" + c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) + indent-tabs-mode:nil + fill-column:99 + End: +*/ +// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 : diff --git a/src/2geom/interval.h b/src/2geom/interval.h index fd446a1c..c5bb4a30 100644 --- a/src/2geom/interval.h +++ b/src/2geom/interval.h @@ -113,6 +113,11 @@ class Interval if (v >= max()) return 1; return timeAt(v); } + /// Find furthest time in [0,1] that maps to the given value. */ + Coord furthestTime(Coord v) { + if (v <= 0.5) return 1; + return 0; + } /// @} /// @name Test coordinates and other intervals for inclusion. diff --git a/src/2geom/path.cpp b/src/2geom/path.cpp index 2e2bce9f..00932c9e 100644 --- a/src/2geom/path.cpp +++ b/src/2geom/path.cpp @@ -747,6 +747,135 @@ PathTime Path::nearestTime(Point const &p, Coord *dist) const return ret; } +std::vector Path::furthestTimePerCurve(Point const &p) const +{ + // return a single furthest time for each curve in this path + std::vector np; + for (const_iterator it = begin(); it != end_default(); ++it) { + np.push_back(it->furthestTime(p)); + } + return np; +} + +std::vector Path::allFurthestTimes(Point const &_point, double from, double to) const +{ + // TODO from and to are not used anywhere. + // rewrite this to simplify. + using std::swap; + + if (from > to) + 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_furthest = _path[si].allFurthestTimes(_point, st, et); + for (unsigned int i = 0; i < all_furthest.size(); ++i) { + all_furthest[i] = si + all_furthest[i]; + } + return all_furthest; + } + std::vector all_t; + std::vector > all_fp; + all_fp.push_back(_path[si].allFurthestTimes(_point, st)); + std::vector ni; + ni.push_back(si); + double dsq; + double maxdistsq = distanceSq(_point, _path[si].pointAt(all_fp.front().front())); + Rect bb(Geom::Point(0, 0), Geom::Point(0, 0)); + for (unsigned int i = si + 1; i < ei; ++i) { + bb = (_path[i].boundsFast()); + dsq = distanceSq(_point, bb); + if (maxdistsq > dsq) + continue; + all_t = _path[i].allFurthestTimes(_point); + dsq = distanceSq(_point, _path[i].pointAt(all_t.front())); + if (maxdistsq < dsq) { + all_fp.clear(); + all_fp.push_back(all_t); + ni.clear(); + ni.push_back(i); + maxdistsq = dsq; + } else if (maxdistsq == dsq) { + all_fp.push_back(all_t); + ni.push_back(i); + } + } + bb = (_path[ei].boundsFast()); + dsq = distanceSq(_point, bb); + if (maxdistsq < dsq) { + all_t = _path[ei].allFurthestTimes(_point, 0, et); + dsq = distanceSq(_point, _path[ei].pointAt(all_t.front())); + if (maxdistsq < dsq) { + for (unsigned int i = 0; i < all_t.size(); ++i) { + all_t[i] = ei + all_t[i]; + } + return all_t; + } else if (maxdistsq == dsq) { + all_fp.push_back(all_t); + ni.push_back(ei); + } + } + std::vector all_furthest; + for (unsigned int i = 0; i < all_fp.size(); ++i) { + for (unsigned int j = 0; j < all_fp[i].size(); ++j) { + all_furthest.push_back(ni[i] + all_fp[i][j]); + } + } + all_furthest.erase(std::unique(all_furthest.begin(), all_furthest.end()), all_furthest.end()); + return all_furthest; +} + +PathTime Path::furthestTime(Point const &p, Coord *dist) const +{ + Coord maxdist = 0; + PathTime ret; + + if (_data->curves.size() == 1) { + // naked moveto + ret.curve_index = 0; + ret.t = 0; + if (dist) { + *dist = distance(_closing_seg->initialPoint(), p); + } + return ret; + } + + for (size_type i = 0; i < size_default(); ++i) { + Curve const &c = at(i); + if (distance(p, c.boundsFast()) < maxdist) continue; + + Coord t = c.furthestTime(p); + Coord d = distance(c.pointAt(t), p); + if (d > maxdist) { + maxdist = d; + ret.curve_index = i; + ret.t = t; + } + } + if (dist) { + *dist = maxdist; + } + + return ret; +} + std::vector Path::nodes() const { std::vector result; diff --git a/src/2geom/path.h b/src/2geom/path.h index 3552ec43..1c1806bc 100644 --- a/src/2geom/path.h +++ b/src/2geom/path.h @@ -573,10 +573,15 @@ class Path std::vector allNearestTimes(Point const &p) const { return allNearestTimes(p, 0, size_default()); } + std::vector allFurthestTimes(Point const &p, Coord from, Coord to) const; + std::vector allFurthestTimes(Point const &p) const { + return allFurthestTimes(p, 0, size_default()); + } PathTime nearestTime(Point const &p, Coord *dist = NULL) const; + PathTime furthestTime(Point const &p, Coord *dist = NULL) const; std::vector nearestTimePerCurve(Point const &p) const; - + std::vector furthestTimePerCurve(Point const &p) const; std::vector nodes() const; void appendPortionTo(Path &p, Coord f, Coord t) const; @@ -853,6 +858,11 @@ inline Coord nearest_time(Point const &p, Path const &c) { return pt.curve_index + pt.t; } +inline Coord furthest_time(Point const &p, Path const &c) { + PathTime pt = c.furthestTime(p); + return pt.curve_index + pt.t; +} + bool are_near(Path const &a, Path const &b, Coord precision = EPSILON); std::ostream &operator<<(std::ostream &out, Path const &path); diff --git a/src/2geom/pathvector.cpp b/src/2geom/pathvector.cpp index 28ab2623..0524fafa 100644 --- a/src/2geom/pathvector.cpp +++ b/src/2geom/pathvector.cpp @@ -280,6 +280,50 @@ std::vector PathVector::allNearestTimes(Point const &p, Coord *d return retval; } + +boost::optional PathVector::furthestTime(Point const &p, Coord *dist) const +{ + boost::optional retval; + + Coord maxdist = 0; + for (size_type i = 0; i < size(); ++i) { + Coord d; + PathTime pos = (*this)[i].furthestTime(p, &d); + if (d > maxdist) { + maxdist = d; + retval = PathVectorTime(i, pos.curve_index, pos.t); + } + } + + if (dist) { + *dist = maxdist; + } + return retval; +} + +std::vector PathVector::allFurthestTimes(Point const &p, Coord *dist) const +{ + std::vector retval; + + Coord maxdist = 0; + for (size_type i = 0; i < size(); ++i) { + Coord d; + PathTime pos = (*this)[i].furthestTime(p, &d); + if (d > maxdist) { + maxdist = d; + retval.clear(); + } + if (d > maxdist) { + retval.push_back(PathVectorTime(i, pos.curve_index, pos.t)); + } + } + + if (dist) { + *dist = maxdist; + } + return retval; +} + std::vector PathVector::nodes() const { std::vector result; diff --git a/src/2geom/pathvector.h b/src/2geom/pathvector.h index 0cbe5bf5..8fc37726 100644 --- a/src/2geom/pathvector.h +++ b/src/2geom/pathvector.h @@ -272,6 +272,9 @@ class PathVector boost::optional nearestTime(Point const &p, Coord *dist = NULL) const; std::vector allNearestTimes(Point const &p, Coord *dist = NULL) const; + boost::optional furthestTime(Point const &p, Coord *dist = NULL) const; + std::vector allFurthestTimes(Point const &p, Coord *dist = NULL) const; + std::vector nodes() const; private: diff --git a/src/2geom/ray.h b/src/2geom/ray.h index 4e60fd81..d6cb515c 100644 --- a/src/2geom/ray.h +++ b/src/2geom/ray.h @@ -105,6 +105,9 @@ class Ray { if (t < 0) t = 0; return t; } + Coord furthestTime(Point const) const { + return infinity(); + } Ray reverse() const { Ray result; result.setOrigin(_origin); diff --git a/src/2geom/sbasis-curve.h b/src/2geom/sbasis-curve.h index cfc4ee9a..a20b06c6 100644 --- a/src/2geom/sbasis-curve.h +++ b/src/2geom/sbasis-curve.h @@ -39,6 +39,7 @@ #include <2geom/curve.h> #include <2geom/exception.h> #include <2geom/nearest-time.h> +#include <2geom/furthest-time.h> #include <2geom/sbasis-geometric.h> #include <2geom/transforms.h> @@ -115,6 +116,14 @@ class SBasisCurve : public Curve { { return all_nearest_times(p, inner, from, to); } + virtual Coord furthestTime( Point const& p, Coord from = 0, Coord to = 1 ) const { + return furthest_time(p, inner, from, to); + } + virtual std::vector allFurthestTimes( Point const& p, Coord from = 0, + Coord to = 1 ) const + { + return all_furthest_times(p, inner, from, to); + } virtual Coord length(Coord tolerance) const { return ::Geom::length(inner, tolerance); } virtual Curve *portion(Coord f, Coord t) const { return new SBasisCurve(Geom::portion(inner, f, t)); diff --git a/src/toys/CMakeLists.txt b/src/toys/CMakeLists.txt index 3090b2f0..7f13577f 100644 --- a/src/toys/CMakeLists.txt +++ b/src/toys/CMakeLists.txt @@ -34,6 +34,7 @@ convole curvature-curve curvature-test curve-curve-distance +curve-curve-furthest-time curve-curve-nearest-time curve-intersection-by-bezier-clipping curve-intersection-by-implicitization @@ -78,6 +79,7 @@ path-effects pencil pencil-2 plane3d +point-curve-furthest-time point-curve-nearest-time portion-test precise-flat diff --git a/src/toys/curve-curve-furthest-time.cpp b/src/toys/curve-curve-furthest-time.cpp new file mode 100644 index 00000000..17a47882 --- /dev/null +++ b/src/toys/curve-curve-furthest-time.cpp @@ -0,0 +1,608 @@ +/* + * Furthest Points Toy 3 + * + * Authors: + * Nathan Hurst + * Marco Cecchetti + * + * Copyright 2008 authors + * + * This library is free software; you can redistribute it and/or + * modify it either under the terms of the GNU Lesser General Public + * License version 2.1 as published by the Free Software Foundation + * (the "LGPL") or, at your option, under the terms of the Mozilla + * Public License Version 1.1 (the "MPL"). If you do not alter this + * notice, a recipient may use your version of this file under either + * the MPL or the LGPL. + * + * You should have received a copy of the LGPL along with this library + * in the file COPYING-LGPL-2.1; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * You should have received a copy of the MPL along with this library + * in the file COPYING-MPL-1.1 + * + * The contents of this file are subject to the Mozilla Public License + * Version 1.1 (the "License"); you may not use this file except in + * compliance with the License. You may obtain a copy of the License at + * http://www.mozilla.org/MPL/ + * + * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY + * OF ANY KIND, either express or implied. See the LGPL or the MPL for + * the specific language governing rights and limitations. + */ + +#include <2geom/d2.h> +#include <2geom/sbasis.h> +#include <2geom/path.h> +#include <2geom/bezier-to-sbasis.h> +#include <2geom/sbasis-geometric.h> +#include <2geom/piecewise.h> +#include <2geom/path-intersection.h> + +#include +#include + +#include + + +using namespace Geom; + + +class fp_finder +{ +public: + fp_finder(cairo_t* _cr, D2 const& _c1, D2 const& _c2) + : cr(_cr), cc1(_c1), cc2(_c2), c1(_c1), c2(_c2) + { + + dc1 = derivative(_c1); + dc2 = derivative(_c2); + cd1 = dot(_c1,dc1); + cd2 = dot(_c2,dc2); + dsq = 10e30; + + Piecewise< D2 > uv1 = unitVector(dc1, EPSILON); + Piecewise< D2 > uv2 = unitVector(dc2, EPSILON); + + dcn1 = dot(Piecewise< D2 >(dc1), uv1); + dcn2 = dot(Piecewise< D2 >(dc2), uv2); + + r_dcn1 = cross(derivative(uv1), uv1); + r_dcn2 = cross(derivative(uv2), uv2); + + k1 = Geom::divide(r_dcn1, dcn1, EPSILON, 3); + k2 = Geom::divide(r_dcn2, dcn2, EPSILON, 3); + + + n1 = divide(rot90(uv1), k1, EPSILON, 3); + n2 = divide(rot90(uv2), k2, EPSILON, 3); + + std::vector cuts1, cuts2; + + // add cuts at points where the curvature is discontinuos + for ( unsigned int i = 1; i < k1.size(); ++i ) + { + if( !are_near(k1[i-1].at1(), k1[i].at0()) ) + { + cuts1.push_back(k1.cuts[i]); + } + } + for ( unsigned int i = 1; i < k2.size(); ++i ) + { + if( !are_near(k2[i-1].at1(), k2[i].at0()) ) + { + cuts2.push_back(k2.cuts[i]); + } + } + + c1 = partition(c1, cuts1); + c2 = partition(c2, cuts2); + +// std::cerr << "# k1 discontinuitis" << std::endl; +// for( unsigned int i = 0; i < cuts1.size(); ++i ) +// { +// std::cerr << "[" << i << "]= " << cuts1[i] << std::endl; +// } +// std::cerr << "# k2 discontinuitis" << std::endl; +// for( unsigned int i = 0; i < cuts2.size(); ++i ) +// { +// std::cerr << "[" << i << "]= " << cuts2[i] << std::endl; +// } + + // add cuts at points were the curvature is zero + std::vector k1_roots = roots(k1); + std::vector k2_roots = roots(k2); + std::sort(k1_roots.begin(), k1_roots.end()); + std::sort(k2_roots.begin(), k2_roots.end()); + c1 = partition(c1, k1_roots); + c2 = partition(c2, k2_roots); + +// std::cerr << "# k1 zeros" << std::endl; +// for( unsigned int i = 0; i < k1_roots.size(); ++i ) +// { +// std::cerr << "[" << i << "]= " << k1_roots[i] << std::endl; +// } +// std::cerr << "# k2 zeros" << std::endl; +// for( unsigned int i = 0; i < k2_roots.size(); ++i ) +// { +// std::cerr << "[" << i << "]= " << k2_roots[i] << std::endl; +// } + + + cairo_set_line_width(cr, 0.2); +// cairo_set_source_rgba(cr, 0.0, 0.5, 0.0, 1.0); +// for( unsigned int i = 1; i < c1.size(); ++i ) +// { +// draw_circ(cr, c1[i].at0() ); +// } +// for( unsigned int i = 1; i < c2.size(); ++i ) +// { +// draw_circ(cr, c2[i].at0() ); +// } + + + // add cuts at furthest points to the other curve cuts points + cuts1.clear(); + cuts1.reserve(c1.size()+1); + for ( unsigned int i = 0; i < c1.size(); ++i ) + { + cuts1.push_back( furthest_time(c1[i].at0(), _c2, dc2, cd2) ); + } + cuts1.push_back( furthest_time(c1[c1.size()-1].at1(), _c2, dc2, cd2) ); + +// for ( unsigned int i = 0; i < c1.size(); ++i ) +// { +// cairo_move_to( cr, c1[i].at0() ); +// cairo_line_to(cr, c2(cuts1[i]) ); +// } +// cairo_move_to( cr, c1[c1.size()-1].at1() ); +// cairo_line_to(cr, c2(cuts1[c1.size()])); + + std::sort(cuts1.begin(), cuts1.end()); + + cuts2.clear(); + cuts2.reserve(c2.size()+1); + for ( unsigned int i = 0; i < c2.size(); ++i ) + { + cuts2.push_back( furthest_time(c2[i].at0(), _c1, dc1, cd1) ); + } + cuts2.push_back( furthest_time(c2[c2.size()-1].at1(), _c1, dc1, cd1) ); + +// for ( unsigned int i = 0; i < c2.size(); ++i ) +// { +// cairo_move_to( cr, c2[i].at0() ); +// cairo_line_to(cr, c1(cuts2[i]) ); +// } +// cairo_move_to( cr, c2[c2.size()-1].at1() ); +// cairo_line_to(cr, c1(cuts2[c2.size()])); +// cairo_stroke(cr); + + std::sort(cuts2.begin(), cuts2.end()); + + c1 = partition(c1, cuts2); + c2 = partition(c2, cuts1); + + + // copy curve to preserve cuts status + Piecewise< D2 > pwc1 = c1; + n1 = partition(n1, pwc1.cuts); + pwc1 = partition(pwc1, n1.cuts); + r_dcn1 = partition(r_dcn1, n1.cuts); + Piecewise< D2 > pwc2 = c2; + n2 = partition(n2, pwc2.cuts); + pwc2 = partition(pwc2, n2.cuts); + + assert( pwc1.size() == n1.size() ); + assert( pwc2.size() == n2.size() ); + assert( r_dcn1.size() == n1.size() ); + + // add cuts at curvature max and min points + Piecewise dk1 = derivative(k1); + Piecewise dk2 = derivative(k2); + std::vector dk1_roots = roots(dk1); + std::vector dk2_roots = roots(dk2); + std::sort(dk1_roots.begin(), dk1_roots.end()); + std::sort(dk2_roots.begin(), dk2_roots.end()); + + c1 = partition(c1, dk1_roots); + c2 = partition(c2, dk2_roots); + +// std::cerr << "# k1 min/max" << std::endl; +// for( unsigned int i = 0; i < dk1_roots.size(); ++i ) +// { +// std::cerr << "[" << i << "]= " << dk1_roots[i] << std::endl; +// } +// std::cerr << "# k2 min/max" << std::endl; +// for( unsigned int i = 0; i < dk2_roots.size(); ++i ) +// { +// std::cerr << "[" << i << "]= " << dk2_roots[i] << std::endl; +// } + +// cairo_set_source_rgba(cr, 0.0, 0.0, 0.6, 1.0); +// for( unsigned int i = 0; i < dk1_roots.size(); ++i ) +// { +// draw_handle(cr, c1(dk1_roots[i])); +// } +// for( unsigned int i = 0; i < dk2_roots.size(); ++i ) +// { +// draw_handle(cr, c2(dk2_roots[i])); +// } + + + // add cuts at furthest points to max and min curvature points + // of the other curve + cuts1.clear(); + cuts1.reserve(dk2_roots.size()); + for ( unsigned int i = 0; i < dk2_roots.size(); ++i ) + { + cuts1.push_back(furthest_time(_c2(dk2_roots[i]), _c1, dc1, cd1)); + } + +// for( unsigned int i = 0; i < dk2_roots.size(); ++i ) +// { +// cairo_move_to(cr, c2(dk2_roots[i])); +// cairo_line_to(cr, c1(cuts1[i])); +// } +// cairo_stroke(cr); + + std::sort(cuts1.begin(), cuts1.end()); + c1 = partition(c1, cuts1); + + + // swap normal vector direction and fill the skip list + skip_list.clear(); + skip_list.resize(c1.size(), false); + double fpt; + Point p, nv; + unsigned int si; + for ( unsigned int i = 0; i < pwc1.size(); ++i ) + { + p = pwc1[i](0.5); + nv = n1[i](0.5); + fpt = furthest_time(p, _c2, dc2, cd2); + if( dot( _c2(fpt) - p, nv ) > 0 ) + { + if ( dot( nv, n2(fpt) ) > 0 ) + { + n1[i] = -n1[i]; + r_dcn1[i] = -r_dcn1[i]; + } + else + { + si = c1.segN( n1.mapToDomain(0.5, i) ); + skip_list[si] = true; + } + } + } + + + for ( unsigned int i = 0; i < pwc2.size(); ++i ) + { + p = pwc2[i](0.5); + nv = n2[i](0.5); + fpt = furthest_time(p, _c1, dc1, cd1); + if( dot( _c1(fpt) - p, nv ) > 0 ) + { + if ( dot( nv, n1(fpt) ) > 0 ) + { + n2[i] = -n2[i]; + } + } + } + + + evl1 = c1 + n1; + evl2 = c2 + n2; + +// cairo_set_source_rgba(cr, 0.3, 0.3, 0.3, 1.0); +// for ( unsigned int i = 0; i < c1.size(); ++i ) +// { +// double t = c1.mapToDomain(0.5, i); +// cairo_move_to(cr, c1(t)); +// cairo_line_to(cr, c1(t) + 30*unit_vector(n1(t))); +// } +// +// for ( unsigned int i = 0; i < c2.size(); ++i ) +// { +// double t = c2.mapToDomain(0.5, i); +// cairo_move_to(cr, c2(t)); +// cairo_line_to(cr, c2(t) + 30*unit_vector(n2(t))); +// } +// cairo_stroke(cr); + + std::cerr << "# skip list: "; + for( unsigned int i = 0; i < c1.cuts.size(); ++i ) + { + if ( skip_list[i] ) + std::cerr << i << " "; + } + std::cerr << std::endl; + + cairo_set_line_width(cr, 0.4); + cairo_set_source_rgba(cr, 0.6, 0.0, 0.0, 1.0); + for( unsigned int i = 0; i < c1.size(); ++i ) + { + if ( skip_list[i] ) + { + cairo_move_to(cr, c1[i].at0()); + cairo_line_to(cr, c1[i].at1()); + } + } + cairo_stroke(cr); + + cairo_set_source_rgba(cr, 0.2, 0.2, 0.2, 1.0); + for( unsigned int i = 1; i < c1.size(); ++i ) + { + draw_circ(cr, c1[i].at0() ); + } + cairo_stroke(cr); + + std::cerr << "# c1 cuts: " << std::endl; + for( unsigned int i = 0; i < c1.cuts.size(); ++i ) + { + std::cerr << "c1.cuts[" << i << "]= " << c1.cuts[i] << std::endl; + } + + } + + void operator() () + { + furthest_times_impl(); + d = sqrt(dsq); + } + + Point firstPoint() const + { + return p1; + } + + Point secondPoint() const + { + return p2; + } + + double firstValue() const + { + return t1; + } + + double secondValue() const + { + return t2; + } + + double distance() const + { + return d; + } + +private: + void furthest_times_impl() + { + double t; + for ( unsigned int i = 0; i < c1.size(); ++i ) + { + if ( skip_list[i] ) continue; + std::cerr << i << " "; + t = c1.mapToDomain(0.5, i); + std::pair fpc = loc_furthest_times(t, c1.cuts[i], c1.cuts[i+1]); + if ( fpc.second != -1 && dsq > L2sq(c1(fpc.first) - c2(fpc.second)) ) + { + t1 = fpc.first; + t2 = fpc.second; + p1 = c1(t1); + p2 = c2(t2); + dsq = L2sq(p1 - p2); + } + } + } + + std::pair + loc_furthest_times( double t, double from = 0, double to = 1 ) + { + std::cerr << "[" << from << "," << to << "] t: " << t << std::endl; + unsigned int iter = 0, iter1 = 0, iter2 = 0; + std::pair fp(-1,-1); + std::pair fpf(from, -1); + std::pair fpt(to, -1); + double ct = t; + double pt = -1; + double s = furthest_time(c1(t), cc2, dc2, cd2); + cairo_set_source_rgba(cr, 1/(t+1), t*t, t, 1.0); + cairo_move_to(cr, c1(t)); + while( !are_near(ct, pt) && iter < 1000 ) + { + pt = ct; + double angle = angle_between( n1(ct), evl2(s) - evl1(ct) ); + assert( !IS_NAN(angle) ); + angle = (angle > 0) ? angle - M_PI : angle + M_PI; + if ( std::fabs(angle) < M_PI/12 ) + { + ++iter2; +// cairo_move_to(cr, c1(ct)); +// cairo_line_to(cr, evl1(ct)); +// cairo_line_to(cr, evl2(s)); + //std::cerr << "s: " << s << std::endl; + //std::cerr << "t: " << ct << std::endl; + + ct = ct + angle / r_dcn1(ct); + s = furthest_time(c1(ct), cc2, dc2, cd2); +// angle = angle_between( n2(s), evl1(ct) - evl2(s) ); +// assert( !IS_NAN(angle) ); +// angle = (angle > 0) ? angle - M_PI : angle + M_PI; +// s = s + angle / (dcn2(s) * k2(s)); + } + else + { + ++iter1; + ct = furthest_time(c2(s), cc1, dc1, cd1, from, to); + s = furthest_time(c1(ct), cc2, dc2, cd2); + } + iter = iter1 + iter2; + //std::cerr << "s: " << s << std::endl; + //std::cerr << "t: " << ct << std::endl; + //cairo_line_to(cr, c2(s)); + //cairo_line_to(cr, c1(ct)); + //std::cerr << "d(pt, ct) = " << std::fabs(ct - pt) << std::endl; + if ( ct < from ) + { + std::cerr << "break left" << std::endl; + fp = fpf; + break; + } + if ( ct > to ) + { + std::cerr << "break right" << std::endl; + fp =fpt; + break; + } + } + //std::cerr << "\n \n"; + std::cerr << "iterations: " << iter1 << " + " << iter2 << " = "<< iter << std::endl; + assert(iter < 3000); + //cairo_move_to(cr, c1(ct)); + //cairo_line_to(cr, c2(s)); + cairo_stroke(cr); + fp.first = ct; + fp.second = s; + return fp; + } + + double furthest_time( Point const& p, D2 const&c, D2 const& dc, SBasis const& cd, double from = 0, double to = 1 ) + { + D2 sbc = c - p; + SBasis dd = cd - dotp(p, dc); + std::vector zeros = roots(dd); + double closest = from; + double distsq = L2sq(sbc(from)); + for ( unsigned int i = 0; i < zeros.size(); ++i ) + { + if ( distsq > L2sq(sbc(zeros[i])) ) + { + closest = zeros[i]; + distsq = L2sq(sbc(closest)); + } + } + if ( distsq > L2sq(sbc(to)) ) + closest = to; + return closest; + } + + SBasis dotp(Point const& p, D2 const& c) + { + SBasis d; + d.resize(c[X].size()); + for ( unsigned int i = 0; i < c[0].size(); ++i ) + { + for( unsigned int j = 0; j < 2; ++j ) + d[i][j] = p[X] * c[X][i][j] + p[Y] * c[Y][i][j]; + } + return d; + } + + Piecewise< D2 > + divide( Piecewise< D2 > const& a, Piecewise const& b, double tol, unsigned int k, double zero=1.e-3) + { + D2< Piecewise > aa = make_cuts_independent(a); + D2< Piecewise > q(Geom::divide(aa[0], b, tol, k, zero), Geom::divide(aa[1], b, tol, k, zero)); + return sectionize(q); + } + + struct are_near_ + { + bool operator() (double x, double y, double eps = Geom::EPSILON ) + { + return are_near(x, y, eps); + } + }; + +private: + cairo_t* cr; + D2 const& cc1, cc2; + Piecewise< D2 > c1, c2; + D2 dc1, dc2; + SBasis cd1, cd2; + Piecewise< D2 > n1, n2, evl1, evl2; + Piecewise k1, k2, dcn1, dcn2, r_dcn1, r_dcn2; + double t1, t2, d, dsq; + Point p1, p2; + std::vector skip_list; +}; + + + + +class FurthestPoints : public Toy +{ + private: + void draw( cairo_t *cr, std::ostringstream *notify, + int width, int height, bool save, std::ostringstream *timer_stream) + { + cairo_set_line_width (cr, 0.3); + D2 A = pshA.asBezier(); + cairo_d2_sb(cr, A); + D2 B = pshB.asBezier(); + cairo_d2_sb(cr, B); + cairo_stroke(cr); + + fp_finder fp(cr, A, B); + Path AP, BP; + AP.append(A); BP.append(B); + Crossings ip_list = curve_sweep(AP, BP); + if( ip_list.empty() ) + { + fp(); + cairo_set_line_width (cr, 0.4); + cairo_set_source_rgba(cr, 0.7, 0.0, 0.7, 1.0); + cairo_move_to(cr, fp.firstPoint()); + cairo_line_to(cr, fp.secondPoint()); + cairo_stroke(cr); + //std::cerr << "fp: (" << fp.firstValue() << "," << fp.secondValue() << ")" << std::endl; + } + Toy::draw(cr, notify, width, height, save,timer_stream); + } + + public: + FurthestPoints(unsigned int _A_bez_ord, unsigned int _B_bez_ord) + : A_bez_ord(_A_bez_ord), B_bez_ord(_B_bez_ord) + { + handles.push_back(&pshA); + handles.push_back(&pshB); + for ( unsigned int i = 0; i < A_bez_ord; ++i ) + pshA.push_back(Geom::Point(uniform()*400, uniform()*400)); + for ( unsigned int i = 0; i < B_bez_ord; ++i ) + pshB.push_back(Geom::Point(uniform()*400, uniform()*400)); + + } + + private: + PointSetHandle pshA, pshB; + unsigned int A_bez_ord; + unsigned int B_bez_ord; +}; + + +int main(int argc, char **argv) +{ + unsigned int A_bez_ord=8; + unsigned int B_bez_ord=5; + if(argc > 2) + sscanf(argv[2], "%d", &B_bez_ord); + if(argc > 1) + sscanf(argv[1], "%d", &A_bez_ord); + + init( argc, argv, new FurthestPoints(A_bez_ord, B_bez_ord)); + return 0; +} + + +/* + 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:fileencoding=utf-8:textwidth=99 : diff --git a/src/toys/elliptical-arc-toy.cpp b/src/toys/elliptical-arc-toy.cpp index 64cb6989..277e3b13 100644 --- a/src/toys/elliptical-arc-toy.cpp +++ b/src/toys/elliptical-arc-toy.cpp @@ -62,6 +62,7 @@ class EllipticalArcToy: public Toy TEST_PORTION, TEST_REVERSE, TEST_NEAREST_POINTS, + TEST_FURTHEST_POINTS, TEST_DERIVATIVE, TEST_ROOTS, TEST_BOUNDS, @@ -390,7 +391,28 @@ class EllipticalArcToy: public Toy } cairo_stroke(cr); } + + void init_fp() + { + init_common(); + fph.pos = Point(10,10); + handles.push_back(&fph); + } + void draw_fp(cairo_t *cr, std::ostringstream *notify, + int width, int height, bool save, std::ostringstream */*timer_stream*/) + { + draw_common(cr, notify, width, height, save); + if ( no_solution || point_overlap ) return; + + std::vector times = ea.allFurthestTimes( fph.pos ); + for ( unsigned int i = 0; i < times.size(); ++i ) + { + cairo_move_to(cr,fph.pos); + cairo_line_to( cr, ea.pointAt(times[i]) ); + } + cairo_stroke(cr); + } void init_derivative() { @@ -700,22 +722,26 @@ class EllipticalArcToy: public Toy draw_f = &EllipticalArcToy::draw_np; break; case 'G': + init_fp(); + draw_f = &EllipticalArcToy::draw_fp; + break; + case 'H': init_derivative(); draw_f = &EllipticalArcToy::draw_derivative; break; - case 'H': + case 'I': init_roots(); draw_f = &EllipticalArcToy::draw_roots; break; - case 'I': + case 'J': init_bounds(); draw_f = &EllipticalArcToy::draw_bounds; break; - case 'J': + case 'K': init_fitting(); draw_f = &EllipticalArcToy::draw_fitting; break; - case 'K': + case 'L': init_transform(); draw_f = &EllipticalArcToy::draw_transform; break; @@ -851,7 +877,7 @@ class EllipticalArcToy: public Toy bool no_solution, point_overlap; bool degenerate; PointHandle initial_point, final_point; - PointHandle nph, ph; + PointHandle nph, ph, fph; std::vector toggles; std::vector sliders; EllipticalArc ea; @@ -871,6 +897,7 @@ const char* EllipticalArcToy::menu_items[] = "portion, pointAt", "reverse, valueAt", "nearest points", + "furthest points", "derivative", "roots", "bounding box", diff --git a/src/toys/furthest-times.cpp b/src/toys/furthest-times.cpp new file mode 100644 index 00000000..f0ff1601 --- /dev/null +++ b/src/toys/furthest-times.cpp @@ -0,0 +1,257 @@ +/* + * Furthest Points Toy high based on Nearest Points Toy + * + * Authors: + * Jabier Arraiza + * + * Copyright 2008 authors + * + * This library is free software; you can redistribute it and/or + * modify it either under the terms of the GNU Lesser General Public + * License version 2.1 as published by the Free Software Foundation + * (the "LGPL") or, at your option, under the terms of the Mozilla + * Public License Version 1.1 (the "MPL"). If you do not alter this + * notice, a recipient may use your version of this file under either + * the MPL or the LGPL. + * + * You should have received a copy of the LGPL along with this library + * in the file COPYING-LGPL-2.1; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * You should have received a copy of the MPL along with this library + * in the file COPYING-MPL-1.1 + * + * The contents of this file are subject to the Mozilla Public License + * Version 1.1 (the "License"); you may not use this file except in + * compliance with the License. You may obtain a copy of the License at + * http://www.mozilla.org/MPL/ + * + * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY + * OF ANY KIND, either express or implied. See the LGPL or the MPL for + * the specific language governing rights and limitations. + */ + +#include <2geom/d2.h> +#include <2geom/sbasis.h> +#include <2geom/path.h> +#include <2geom/bezier-to-sbasis.h> + +#include +#include + + +using namespace Geom; + + +class fp_finder +{ +public: + fp_finder(cairo_t* _cr, D2 const& _c1, D2 const& _c2) + : cr(_cr), c1(_c1), c2(_c2) + { + dc1 = derivative(c1); + dc2 = derivative(c2); + cd1 = dot(c1,dc1); + cd2 = dot(c2,dc2); + dsq = 10e30; + } + + void operator() () + { + furthest_times_impl(0.5, 0, 1); + d = sqrt(dsq); + } + + Point firstPoint() const + { + return p1; + } + + Point secondPoint() const + { + return p2; + } + + double firstValue() const + { + return t1; + } + + double secondValue() const + { + return t2; + } + + double distance() const + { + return d; + } + +private: + void furthest_times_impl( double t, double from = 0, double to = 1 ) + { + //std::cerr << "[" << from << "," << to << "] t: " << t << std::endl; + + double st = t, et = t; + std::pair fpc = loc_furthest_times(t, from, to); + //std::cerr << "(" << fpc.first << "," << fpc.second << ")" << std::endl; + if ( fpc.second != -1 && dsq > L2sq(c1(fpc.first) - c2(fpc.second)) ) + { + t1 = fpc.first; + t2 = fpc.second; + p1 = c1(t1); + p2 = c2(t2); + dsq = L2sq(p1 - p2); + } + if (fpc.first < t) + st = fpc.first; + else + et = fpc.first; + //std::cerr << "[" << st << "," << et << "]" << std::endl; + double d1 = std::fabs(st - from); + double d2 = std::fabs(to - et); + if ( d1 > EPSILON ) + furthest_times_impl(from + d1/2, from, st); + if ( d2 > EPSILON ) + furthest_times_impl(et + d2/2, et, to); + } + + std::pair + loc_furthest_times( double t, double from = 0, double to = 1 ) + { + unsigned int i = 0; + std::pair fp(-1,-1); + std::pair fpf(from, -1); + std::pair fpt(to, -1); + double ct = t; + double pt = -1; + double s; + //cairo_set_source_rgba(cr, 1/(t+1), t*t, t, 1.0); + //cairo_move_to(cr, c1(t)); + while( !are_near(ct, pt) ) + { + ++i; + pt = ct; + s = furthest_time(c1(ct), c2, dc2, cd2); + //std::cerr << "s: " << s << std::endl; + //cairo_line_to(cr, c2(s)); + ct = furthest_time(c2(s), c1, dc1, cd1, from, to); + //std::cerr << "t: " << t1 << std::endl; + //cairo_line_to(cr, c1(ct)); + if ( ct < from ) return fpf; + if ( ct > to ) return fpt; + } + //std::cerr << "\n \n"; + //std::cerr << "iterations: " << i << std::endl; + cairo_stroke(cr); + fp.first = ct; + fp.second = s; + return fp; + } + + double furthest_time( Point const& p, D2 const&c, D2 const& dc, SBasis const& cd, double from = 0, double to = 1 ) + { + D2 sbc = c - p; + SBasis dd = cd - dotp(p, dc); + std::vector zeros = roots(dd); + double closest = from; + double distsq = L2sq(sbc(from)); + for ( unsigned int i = 0; i < zeros.size(); ++i ) + { + if ( distsq > L2sq(sbc(zeros[i])) ) + { + closest = zeros[i]; + distsq = L2sq(sbc(closest)); + } + } + if ( distsq > L2sq(sbc(to)) ) + closest = to; + return closest; + } + + SBasis dotp(Point const& p, D2 const& c) + { + SBasis d; + d.resize(c[X].size()); + for ( unsigned int i = 0; i < c[0].size(); ++i ) + { + for( unsigned int j = 0; j < 2; ++j ) + d[i][j] = p[X] * c[X][i][j] + p[Y] * c[Y][i][j]; + } + return d; + } + +private: + static const Coord EPSILON = 10e-3; + cairo_t* cr; + D2 const& c1, c2; + D2 dc1, dc2; + SBasis cd1, cd2; + double t1, t2, d, dsq; + Point p1, p2; +}; + + + + +class FurthestPoints : public Toy +{ + private: + void draw( cairo_t *cr, std::ostringstream *notify, + int width, int height, bool save, std::ostringstream *timer_stream) + { + cairo_set_line_width (cr, 0.2); + D2 A = handles_to_sbasis(handles.begin(), A_bez_ord-1); + cairo_d2_sb(cr, A); + D2 B = handles_to_sbasis(handles.begin()+A_bez_ord, B_bez_ord-1); + cairo_d2_sb(cr, B); + + + fp_finder fp(cr, A, B); + fp(); + cairo_move_to(cr, fp.firstPoint()); + cairo_line_to(cr, fp.secondPoint()); + cairo_stroke(cr); + //std::cerr << "fp: (" << fp.firstValue() << "," << fp.secondValue() << ")" << std::endl; + + Toy::draw(cr, notify, width, height, save,timer_stream); + } + + public: + FurthestPoints(unsigned int _A_bez_ord, unsigned int _B_bez_ord) + : A_bez_ord(_A_bez_ord), B_bez_ord(_B_bez_ord) + { + unsigned int total_handles = A_bez_ord + B_bez_ord; + for ( unsigned int i = 0; i < total_handles; ++i ) + handles.push_back(Geom::Point(uniform()*400, uniform()*400)); + } + + private: + unsigned int A_bez_ord; + unsigned int B_bez_ord; +}; + + +int main(int argc, char **argv) +{ + unsigned int A_bez_ord=8; + unsigned int B_bez_ord=5; + if(argc > 2) + sscanf(argv[2], "%d", &B_bez_ord); + if(argc > 1) + sscanf(argv[1], "%d", &A_bez_ord); + + init( argc, argv, new FurthestPoints(A_bez_ord, B_bez_ord)); + return 0; +} + + +/* + 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:fileencoding=utf-8:textwidth=99 : diff --git a/src/toys/furthest-times2.cpp b/src/toys/furthest-times2.cpp new file mode 100644 index 00000000..9c42d556 --- /dev/null +++ b/src/toys/furthest-times2.cpp @@ -0,0 +1,312 @@ +/* + * Furthest Points Toy 2 high based on Nearest Points Toy 2 + * + * Authors: + * Jabier Arraiza + * + * Copyright 2008 authors + * + * This library is free software; you can redistribute it and/or + * modify it either under the terms of the GNU Lesser General Public + * License version 2.1 as published by the Free Software Foundation + * (the "LGPL") or, at your option, under the terms of the Mozilla + * Public License Version 1.1 (the "MPL"). If you do not alter this + * notice, a recipient may use your version of this file under either + * the MPL or the LGPL. + * + * You should have received a copy of the LGPL along with this library + * in the file COPYING-LGPL-2.1; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * You should have received a copy of the MPL along with this library + * in the file COPYING-MPL-1.1 + * + * The contents of this file are subject to the Mozilla Public License + * Version 1.1 (the "License"); you may not use this file except in + * compliance with the License. You may obtain a copy of the License at + * http://www.mozilla.org/MPL/ + * + * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY + * OF ANY KIND, either express or implied. See the LGPL or the MPL for + * the specific language governing rights and limitations. + */ + + +#include <2geom/d2.h> +#include <2geom/sbasis.h> +#include <2geom/path.h> +#include <2geom/bezier-to-sbasis.h> +#include <2geom/sbasis-geometric.h> + +#include +#include + +#include + + +using namespace Geom; + + +class fp_finder +{ +public: + fp_finder(cairo_t* _cr, D2 const& _c1, D2 const& _c2) + : cr(_cr), c1(_c1), c2(_c2) + { + dc1 = derivative(c1); + dc2 = derivative(c2); + cd1 = dot(c1,dc1); + cd2 = dot(c2,dc2); + dsq = 10e30; + + Piecewise k1 = curvature(c1, EPSILON); + Piecewise k2 = curvature(c2, EPSILON); + Piecewise dk1 = derivative(k1); + Piecewise dk2 = derivative(k2); + std::vector k1_roots = roots(k1); + std::vector k2_roots = roots(k2); + std::vector dk1_roots = roots(dk1); + std::vector dk2_roots = roots(dk2); + tlist.clear(); + tlist.resize(k1_roots.size() + k2_roots.size() + dk1_roots.size() + dk2_roots.size() + 4); + tlist.push_back(0); + tlist.insert(tlist.end(), dk1_roots.begin(), dk1_roots.end()); + tlist.insert(tlist.end(), k1_roots.begin(), k1_roots.end()); +// std::cerr << "dk1 roots: "; +// for ( unsigned int i = 0; i < dk1_roots.size(); ++i ) +// { +// std::cerr << dk1_roots[i] << " "; +// } +// std::cerr << "\n"; + for ( unsigned int i = 0; i < k2_roots.size(); ++i ) + { + tlist.push_back(furthest_time(c2(k2_roots[i]), c1, dc1, cd1)); + } + + for ( unsigned int i = 0; i < dk2_roots.size(); ++i ) + { + tlist.push_back(furthest_time(c2(dk2_roots[i]), c1, dc1, cd1)); + } + tlist.push_back(furthest_time(c2(0), c1, dc1, cd1)); + tlist.push_back(furthest_time(c2(1), c1, dc1, cd1)); + tlist.push_back(1); + std::sort(tlist.begin(), tlist.end()); + std::vector::iterator pos + = std::unique(tlist.begin(), tlist.end(), are_near_() ); + if (pos != tlist.end()) + { + tlist.erase(pos, tlist.end()); + } + for( unsigned int i = 0; i < tlist.size(); ++i ) + { + draw_circ(cr, c1(tlist[i]) ); + } + } + + void operator() () + { + //furthest_times_impl( tlist.size() / 2, 0, tlist.size() - 1 ); + furthest_times_impl(); + d = sqrt(dsq); + } + + Point firstPoint() const + { + return p1; + } + + Point secondPoint() const + { + return p2; + } + + double firstValue() const + { + return t1; + } + + double secondValue() const + { + return t2; + } + + double distance() const + { + return d; + } + +private: + void furthest_times_impl() + { + double t; + double from = tlist[0]; + double to; + for ( unsigned int i = 1; i < tlist.size(); ++i ) + { + to = tlist[i]; + t = from + (to - from) / 2 ; + std::pair fpc = loc_furthest_times(t, from, to); + if ( fpc.second != -1 && dsq > L2sq(c1(fpc.first) - c2(fpc.second)) ) + { + t1 = fpc.first; + t2 = fpc.second; + p1 = c1(t1); + p2 = c2(t2); + dsq = L2sq(p1 - p2); + } + from = tlist[i]; + } + } + + std::pair + loc_furthest_times( double t, double from = 0, double to = 1 ) + { + //std::cerr << "[" << from << "," << to << "] t: " << t << std::endl; + unsigned int i = 0; + std::pair fp(-1,-1); + std::pair fpf(from, -1); + std::pair fpt(to, -1); + double ct = t; + double pt = -1; + double s; + //cairo_set_source_rgba(cr, 1/(t+1), t*t, t, 1.0); + cairo_move_to(cr, c1(t)); + while( !are_near(ct, pt) ) + { + ++i; + pt = ct; + s = furthest_time(c1(ct), c2, dc2, cd2); + //std::cerr << "s: " << s << std::endl; + //cairo_line_to(cr, c2(s)); + ct = furthest_time(c2(s), c1, dc1, cd1, from, to); + //std::cerr << "t: " << ct << std::endl; + //cairo_line_to(cr, c1(ct)); + //std::cerr << "d(pt, ct) = " << std::fabs(ct - pt) << std::endl; + if ( ct < from ) return fpf; + if ( ct > to ) return fpt; + } + //std::cerr << "\n \n"; + std::cerr << "iterations: " << i << std::endl; + //cairo_move_to(cr, c1(ct)); + //cairo_line_to(cr, c2(s)); + cairo_stroke(cr); + fp.first = ct; + fp.second = s; + return fp; + } + + double furthest_time( Point const& p, D2 const&c, D2 const& dc, SBasis const& cd, double from = 0, double to = 1 ) + { + D2 sbc = c - p; + SBasis dd = cd - dotp(p, dc); + std::vector zeros = roots(dd); + double closest = from; + double distsq = L2sq(sbc(from)); + for ( unsigned int i = 0; i < zeros.size(); ++i ) + { + if ( distsq > L2sq(sbc(zeros[i])) ) + { + closest = zeros[i]; + distsq = L2sq(sbc(closest)); + } + } + if ( distsq > L2sq(sbc(to)) ) + closest = to; + return closest; + } + + SBasis dotp(Point const& p, D2 const& c) + { + SBasis d; + d.resize(c[X].size()); + for ( unsigned int i = 0; i < c[0].size(); ++i ) + { + for( unsigned int j = 0; j < 2; ++j ) + d[i][j] = p[X] * c[X][i][j] + p[Y] * c[Y][i][j]; + } + return d; + } + + struct are_near_ + { + bool operator() (double x, double y, double eps = Geom::EPSILON ) + { + return are_near(x, y, eps); + } + }; + +private: + static const Coord EPSILON = 10e-5; + cairo_t* cr; + D2 const& c1, c2; + D2 dc1, dc2; + SBasis cd1, cd2; + double t1, t2, d, dsq; + Point p1, p2; + std::vector tlist; +}; + + + + +class FurthestPoints : public Toy +{ + private: + void draw( cairo_t *cr, std::ostringstream *notify, + int width, int height, bool save, std::ostringstream *timer_stream) + { + cairo_set_line_width (cr, 0.2); + D2 A = handles_to_sbasis(handles.begin(), A_bez_ord-1); + cairo_d2_sb(cr, A); + D2 B = handles_to_sbasis(handles.begin()+A_bez_ord, B_bez_ord-1); + cairo_d2_sb(cr, B); + + + fp_finder fp(cr, A, B); + fp(); + cairo_move_to(cr, fp.firstPoint()); + cairo_line_to(cr, fp.secondPoint()); + cairo_stroke(cr); + //std::cerr << "fp: (" << fp.firstValue() << "," << fp.secondValue() << ")" << std::endl; + + Toy::draw(cr, notify, width, height, save,timer_stream); + } + + public: + FurthestPoints(unsigned int _A_bez_ord, unsigned int _B_bez_ord) + : A_bez_ord(_A_bez_ord), B_bez_ord(_B_bez_ord) + { + unsigned int total_handles = A_bez_ord + B_bez_ord; + for ( unsigned int i = 0; i < total_handles; ++i ) + handles.push_back(Geom::Point(uniform()*400, uniform()*400)); + } + + private: + unsigned int A_bez_ord; + unsigned int B_bez_ord; +}; + + +int main(int argc, char **argv) +{ + unsigned int A_bez_ord=8; + unsigned int B_bez_ord=5; + if(argc > 2) + sscanf(argv[2], "%d", &B_bez_ord); + if(argc > 1) + sscanf(argv[1], "%d", &A_bez_ord); + + init( argc, argv, new FurthestPoints(A_bez_ord, B_bez_ord)); + return 0; +} + + +/* + 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:fileencoding=utf-8:textwidth=99 : diff --git a/src/toys/point-curve-furthest-time.cpp b/src/toys/point-curve-furthest-time.cpp new file mode 100644 index 00000000..d699c77b --- /dev/null +++ b/src/toys/point-curve-furthest-time.cpp @@ -0,0 +1,397 @@ +/* + * point-curve furthest point routines testing high based in point-curve-.nearest-point toy + * + * Authors: + * Jabier Arraiza + * + * Copyright 2008 authors + * + * This library is free software; you can redistribute it and/or + * modify it either under the terms of the GNU Lesser General Public + * License version 2.1 as published by the Free Software Foundation + * (the "LGPL") or, at your option, under the terms of the Mozilla + * Public License Version 1.1 (the "MPL"). If you do not alter this + * notice, a recipient may use your version of this file under either + * the MPL or the LGPL. + * + * You should have received a copy of the LGPL along with this library + * in the file COPYING-LGPL-2.1; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * You should have received a copy of the MPL along with this library + * in the file COPYING-MPL-1.1 + * + * The contents of this file are subject to the Mozilla Public License + * Version 1.1 (the "License"); you may not use this file except in + * compliance with the License. You may obtain a copy of the License at + * http://www.mozilla.org/MPL/ + * + * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY + * OF ANY KIND, either express or implied. See the LGPL or the MPL for + * the specific language governing rights and limitations. + */ + +#include <2geom/d2.h> +#include <2geom/sbasis.h> +#include <2geom/path.h> +#include <2geom/curves.h> +#include <2geom/angle.h> +#include <2geom/bezier-to-sbasis.h> +#include <2geom/sbasis-geometric.h> +#include <2geom/piecewise.h> + +#include <2geom/svg-path-parser.h> + +#include +#include + +#include <2geom/transforms.h> +#include <2geom/pathvector.h> + + +#include + +using namespace Geom; + +std::ostream& +operator<< (std::ostream &out, PathVectorTime const &pvp) +{ + return out << pvp.path_index << "." << pvp.curve_index << "." << pvp.t; +} + +class FurthestPoints : public Toy +{ + enum menu_item_t + { + FIRST_ITEM = 1, + LINE_SEGMENT = FIRST_ITEM, + ELLIPTICAL_ARC, + SBASIS_CURVE, + PIECEWISE, + PATH, + PATH_SVGD, + TOTAL_ITEMS + }; + + static const char* menu_items[TOTAL_ITEMS]; + +private: + PathVector paths_b; + void draw( cairo_t *cr, std::ostringstream *notify, + int width, int height, bool save, std::ostringstream *timer_stream) + { + + Point p = ph.pos; + Point fp = p; + std::vector fps; + + cairo_set_line_width (cr, 0.3); + cairo_set_source_rgb(cr, 0,0,0); + switch ( choice ) + { + case '1': + { + LineSegment seg(psh.pts[0], psh.pts[1]); + cairo_move_to(cr, psh.pts[0]); + cairo_curve(cr, seg); + double t = seg.furthestTime(p); + fp = seg.pointAt(t); + if ( toggles[0].on ) + { + fps.push_back(fp); + } + break; + } + case '2': + { + EllipticalArc earc; + bool earc_constraints_satisfied = true; + try + { + earc.set(psh.pts[0], 200, 150, 0, true, true, psh.pts[1]); + } + catch( RangeError e ) + { + earc_constraints_satisfied = false; + } + if ( earc_constraints_satisfied ) + { + cairo_d2_sb(cr, earc.toSBasis()); + if ( toggles[0].on ) + { + std::vector t = earc.allFurthestTimes(p); + for ( unsigned int i = 0; i < t.size(); ++i ) + fps.push_back(earc.pointAt(t[i])); + } + else + { + double t = earc.furthestTime(p); + fp = earc.pointAt(t); + } + } + break; + } + case '3': + { + D2 A = psh.asBezier(); + cairo_d2_sb(cr, A); + if ( toggles[0].on ) + { + std::vector t = Geom::all_furthest_times(p, A); + for ( unsigned int i = 0; i < t.size(); ++i ) + fps.push_back(A(t[i])); + } + else + { + double t = furthest_time(p, A); + fp = A(t); + } + break; + } + case '4': + { + D2 A = handles_to_sbasis(psh.pts.begin(), 3); + D2 B = handles_to_sbasis(psh.pts.begin() + 3, 3); + D2 C = handles_to_sbasis(psh.pts.begin() + 6, 3); + D2 D = handles_to_sbasis(psh.pts.begin() + 9, 3); + cairo_d2_sb(cr, A); + cairo_d2_sb(cr, B); + cairo_d2_sb(cr, C); + cairo_d2_sb(cr, D); + Piecewise< D2 > pwc; + pwc.push_cut(0); + pwc.push_seg(A); + pwc.push_cut(0.25); + pwc.push_seg(B); + pwc.push_cut(0.50); + pwc.push_seg(C); + pwc.push_cut(0.75); + pwc.push_seg(D); + pwc.push_cut(1); + if ( toggles[0].on ) + { + std::vector t = Geom::all_furthest_times(p, pwc); + for ( unsigned int i = 0; i < t.size(); ++i ) + fps.push_back(pwc(t[i])); + } + else + { + double t = Geom::furthest_time(p, pwc); + fp = pwc(t); + } + break; + } + case '5': + { + closed_toggle = true; + BezierCurveN<3> A(psh.pts[0], psh.pts[1], psh.pts[2], psh.pts[3]); + BezierCurveN<2> B(psh.pts[3], psh.pts[4], psh.pts[5]); + BezierCurveN<3> C(psh.pts[5], psh.pts[6], psh.pts[7], psh.pts[8]); + Path path; + path.append(A); + path.append(B); + path.append(C); + EllipticalArc D; + bool earc_constraints_satisfied = true; + try + { + D.set(psh.pts[8], 160, 80, 0, true, true, psh.pts[9]); + } + catch( RangeError e ) + { + earc_constraints_satisfied = false; + } + if ( earc_constraints_satisfied ) path.append(D); + if ( toggles[1].on ) path.close(true); + + cairo_path(cr, path); + + if ( toggles[0].on ) + { + std::vector t = path.allFurthestTimes(p); + for ( unsigned int i = 0; i < t.size(); ++i ) + fps.push_back(path.pointAt(t[i])); + } + else + { + PathTime pt = path.furthestTime(p); + fp = path.pointAt(pt); + } + break; + } + case '6': + { + closed_toggle = true; + PathVector pathv = paths_b*Translate(psh.pts[0]-paths_b[0][0].initialPoint()); + //std::cout << pathv.size() << std::endl; + OptRect optRect = bounds_fast(pathv); + + cairo_rectangle(cr, *optRect); + cairo_stroke(cr); + + cairo_path(cr, pathv); + + if ( toggles[0].on ) + { + std::vector t = pathv.allFurthestTimes(p); + for ( unsigned int i = 0; i < t.size(); ++i ) + fps.push_back(pathv.pointAt(t[i])); + } + else + { + //boost::optional + double s = 0, e = 1; + draw_cross(cr, pathv[0].pointAt(s)); + draw_cross(cr, pathv[0].pointAt(e)); + double t = pathv[0][0].furthestTime(p, 0, 1); + if(t) { + *notify << p+psh.pts[0] << std::endl; + *notify << t << std::endl; + fp = pathv[0].pointAt(t); + } + } + break; + } + default: + { + *notify << std::endl; + for (int i = FIRST_ITEM; i < TOTAL_ITEMS; ++i) + { + *notify << " " << i << " - " << menu_items[i] << std::endl; + } + Toy::draw(cr, notify, width, height, save,timer_stream); + return; + } + } + + if ( toggles[0].on ) + { + for ( unsigned int i = 0; i < fps.size(); ++i ) + { + cairo_move_to(cr, p); + cairo_line_to(cr, fps[i]); + } + } + else + { + cairo_move_to(cr, p); + cairo_line_to(cr, fp); + } + cairo_stroke(cr); + + toggles[0].bounds = Rect( Point(10, height - 50), Point(10, height - 50) + Point(80,25) ); + toggles[0].draw(cr); + if ( closed_toggle ) + { + toggles[1].bounds = Rect( Point(100, height - 50), Point(100, height - 50) + Point(80,25) ); + toggles[1].draw(cr); + closed_toggle = false; + } + + Toy::draw(cr, notify, width, height, save,timer_stream); + } + + void key_hit(GdkEventKey *e) + { + choice = e->keyval; + switch ( choice ) + { + case '1': + total_handles = 2; + break; + case '2': + total_handles = 2; + break; + case '3': + total_handles = 6; + break; + case '4': + total_handles = 13; + break; + case '5': + total_handles = 10; + break; + case '6': + total_handles = 1; + break; + default: + total_handles = 0; + } + psh.pts.clear(); + psh.push_back( 262.6037,35.824151); + psh.pts[0] += Point(300,300); + psh.push_back(0,0); + psh.pts[1] += psh.pts[0]; + psh.push_back(-92.64892,-187.405851); + psh.pts[2] += psh.pts[0]; + psh.push_back(30,-149.999981); + psh.pts[3] += psh.pts[0]; + for ( unsigned int i = 0; i < total_handles; ++i ) + { + psh.push_back(uniform()*400, uniform()*400); + } + ph.pos = Point(uniform()*400, uniform()*400); + redraw(); + } + + void mouse_pressed(GdkEventButton* e) + { + toggle_events(toggles, e); + Toy::mouse_pressed(e); + } + +public: + void first_time(int argc, char** argv) { + const char *path_b_name = "svgd/star.svgd"; + if(argc > 1) + path_b_name = argv[1]; + paths_b = read_svgd(path_b_name); + } + FurthestPoints() + : total_handles(0), choice('0'), closed_toggle(false) + { + handles.push_back(&psh); + handles.push_back(&ph); + ph.pos = Point(20,20); + toggles.push_back( Toggle("ALL FP", false) ); + toggles.push_back( Toggle("CLOSED", false) ); + } + +private: + PointSetHandle psh; + PointHandle ph; + std::vector toggles; + unsigned int total_handles; + char choice; + bool closed_toggle; +}; + +const char* FurthestPoints::menu_items[] = +{ + "", + "LineSegment", + "EllipticalArc", + "SBasisCurve", + "Piecewise", + "Path", + "Path from SVGD" +}; + + + +int main(int argc, char **argv) +{ + init( argc, argv, new FurthestPoints() ); + return 0; +} + + +/* + 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:fileencoding=utf-8:textwidth=99 : diff --git a/src/toys/point-curve-nearest-time.cpp b/src/toys/point-curve-nearest-time.cpp index 8e36dd33..261258cf 100644 --- a/src/toys/point-curve-nearest-time.cpp +++ b/src/toys/point-curve-nearest-time.cpp @@ -341,7 +341,7 @@ class NearestPoints : public Toy public: void first_time(int argc, char** argv) { - const char *path_b_name="star.svgd"; + const char *path_b_name = "svgd/star.svgd"; if(argc > 1) path_b_name = argv[1]; paths_b = read_svgd(path_b_name);