![[tokkee]](http://tokkee.org/images/avatar.png){kind=avatar}
summary | shortlog | log | commit | commitdiff | tree
raw | patch | inline | side by side (parent: 22e5d80)
raw | patch | inline | side by side (parent: 22e5d80)
| author | jfbarraud <jfbarraud@users.sourceforge.net> | |
| Mon, 9 Mar 2009 01:47:39 +0000 (01:47 +0000) | ||
| committer | jfbarraud <jfbarraud@users.sourceforge.net> | |
| Mon, 9 Mar 2009 01:47:39 +0000 (01:47 +0000) |
18 files changed:
index b5461487dd659dcaf0b32c866bf6b3019646f704..3c8669a9d8e6bdcced2cc2d91a2c7fd996ca2546 100644 (file)
--- a/src/2geom/CMakeLists.txt
+++ b/src/2geom/CMakeLists.txt
basic-intersection.h
basic-intersection.cpp
+recursive-bezier-intersection.cpp
geom.cpp
geom.h
index d03799dd7ff75c1e3e9f533e41ed688f00a96424..901e22f40b8033efc7cb561bd61480472ba38bfc 100644 (file)
2geom_lib2geom_a_SOURCES = \
2geom/basic-intersection.cpp \
+ 2geom/bezier-clipping.cpp \
+ 2geom/chebyshev.cpp \
+ 2geom/chebyshev.h \
+ 2geom/utils.cpp \
2geom/bezier-utils.cpp \
2geom/circle-circle.cpp \
2geom/circle.cpp \
2geom/pathvector.cpp \
2geom/piecewise.cpp \
2geom/point.cpp \
- 2geom/poly-dk-solve.cpp \
- 2geom/poly-laguerre-solve.cpp \
2geom/poly.cpp \
2geom/quadtree.cpp \
2geom/region.cpp \
2geom/point-l.h \
2geom/point-ops.h \
2geom/point.h \
- 2geom/poly-dk-solve.h \
- 2geom/poly-laguerre-solve.h \
2geom/poly.h \
2geom/quadtree.h \
2geom/rect.h \
index 16b5c0240f59019caa250848a1a6ea2df6816991..c03023e6fecddbb8bb952e521f1664ff7ae69526 100644 (file)
using std::vector;
namespace Geom {
+//#ifdef USE_RECURSIVE_INTERSECTOR
+
+// void find_intersections(std::vector<std::pair<double, double> > &xs,
+// D2<SBasis> const & A,
+// D2<SBasis> const & B) {
+// vector<Point> BezA, BezB;
+// sbasis_to_bezier(BezA, A);
+// sbasis_to_bezier(BezB, B);
+
+// xs.clear();
+
+// find_intersections_bezier_recursive(xs, BezA, BezB);
+// }
+// void find_intersections(std::vector< std::pair<double, double> > & xs,
+// std::vector<Point> const& A,
+// std::vector<Point> const& B,
+// double precision){
+// find_intersections_bezier_recursive(xs, A, B, precision);
+// }
+
+//#else
+
namespace detail{ namespace bezier_clipping {
void portion (std::vector<Point> & B, Interval const& I);
}; };
vector<Point> BezA, BezB;
sbasis_to_bezier(BezA, A);
sbasis_to_bezier(BezB, B);
-
+
xs.clear();
- find_intersections(xs, BezA, BezB);
+ find_intersections_bezier_clipping(xs, BezA, BezB);
+}
+void find_intersections(std::vector< std::pair<double, double> > & xs,
+ std::vector<Point> const& A,
+ std::vector<Point> const& B,
+ double precision){
+ find_intersections_bezier_clipping(xs, A, B, precision);
}
+//#endif
+
/*
* split the curve at the midpoint, returning an array with the two parts
* Temporary storage is minimized by using part of the storage for the result
in = r;
}
}
+
for(unsigned i = 0; i < dr.size()-1; i++) {
for(unsigned j = i+1; j < dr.size()-1; j++) {
std::vector<std::pair<double, double> > section;
double r = section[k].second;
// XXX: This condition will prune out false positives, but it might create some false negatives. Todo: Confirm it is correct.
if(j == i+1)
- if((l == 1) && (r == 0))
+ //if((l == 1) && (r == 0))
+ if( ( l > 1-1e-4 ) && (r < 1e-4) )//FIXME: what precision should be used here???
continue;
xs.push_back(std::make_pair((1-l)*dr[i] + l*dr[i+1],
(1-r)*dr[j] + r*dr[j+1]));
index 783df370ec006458f102aeabc61ce631f67cd3f5..a19a10c8c59cf2a2ce20f20d51299da1e69b2d7f 100644 (file)
#define USE_RECURSIVE_INTERSECTOR 0
-namespace Geom {
-
-#ifdef USE_RECURSIVE_INTERSECTOR
-void
-find_intersections( std::vector<std::pair<double, double> > &xs,
- D2<SBasis> const & A,
- D2<SBasis> const & B);
-void
-find_self_intersections(std::vector<std::pair<double, double> > &xs,
- D2<SBasis> const & A);
+namespace Geom {
-// Bezier form
-void
-find_intersections( std::vector<std::pair<double, double> > &xs,
- std::vector<Point> const & A,
- std::vector<Point> const & B,
- double precision = 1e-5);
+//why not allowing precision to be set here?
+void find_intersections(std::vector<std::pair<double, double> >& xs,
+ D2<SBasis> const & A,
+ D2<SBasis> const & B);
-// Bezier form
-void
-find_intersections_bezier_clipping( std::vector<std::pair<double, double> > &xs,
- std::vector<Point> const & A,
- std::vector<Point> const & B,
- double precision = 1e-5);
+void find_intersections(std::vector< std::pair<double, double> > & xs,
+ std::vector<Point> const& A,
+ std::vector<Point> const& B,
+ double precision = 1e-5);
-void
-find_self_intersections(std::vector<std::pair<double, double> > &xs,
- std::vector<Point> const & A);
+//why not allowing precision to be set here?
+void find_self_intersections(std::vector<std::pair<double, double> >& xs,
+ D2<SBasis> const & A);
-#else
+//--not implemented
+//void find_self_intersections(std::vector<std::pair<double, double> >& xs,
+// std::vector<Point> const & A);
+
+
+//TODO: this should be moved to .cpp, shouldn't it?
+// #ifdef USE_RECURSIVE_INTERSECTOR
+// /*
+// * find_intersection
+// *
+// * input: A, B - set of control points of two Bezier curve
+// * input: precision - required precision of computation
+// * output: xs - set of pairs of parameter values
+// * at which crossing happens
+// */
+// void find_intersections_bezier_recursive (std::vector< std::pair<double, double> > & xs,
+// std::vector<Point> const& A,
+// std::vector<Point> const& B,
+// double precision = 1e-5);
+// #else
/*
* find_intersection
*
* This routine is based on the Bezier Clipping Algorithm,
* see: Sederberg, Nishita, 1990 - Curve intersection using Bezier clipping
*/
-void find_intersections_clipping (std::vector< std::pair<double, double> > & xs,
+void find_intersections_bezier_clipping (std::vector< std::pair<double, double> > & xs,
std::vector<Point> const& A,
std::vector<Point> const& B,
double precision = 1e-5);
-#endif
+//#endif
+
+
+
/*
* find_collinear_normal
*
D2<SBasis> const &A,
D2<SBasis> const &B);
-void find_intersections(std::vector<std::pair<double, double> >& xs,
- D2<SBasis> const & A,
- D2<SBasis> const & B);
-
-void find_self_intersections(std::vector<std::pair<double, double> >& xs,
- D2<SBasis> const & A);
-
-
-void find_self_intersections(std::vector<std::pair<double, double> >& xs,
- std::vector<Point> const & A);
-
-
/**
* Compute the Hausdorf distance from A to B only.
diff --git a/src/2geom/bezier-clipping.cpp b/src/2geom/bezier-clipping.cpp
--- /dev/null
@@ -0,0 +1,1292 @@
+/*
+ * Implement the Bezier clipping algorithm for finding
+ * Bezier curve intersection points and collinear normals
+ *
+ * Authors:
+ * Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2008 authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+
+
+
+#include <2geom/basic-intersection.h>
+#include <2geom/choose.h>
+#include <2geom/point.h>
+#include <2geom/interval.h>
+#include <2geom/bezier.h>
+//#include <2geom/convex-cover.h>
+#include <2geom/numeric/matrix.h>
+
+#include <cassert>
+#include <vector>
+#include <algorithm>
+#include <utility>
+//#include <iomanip>
+
+
+
+
+#define VERBOSE 0
+#define CHECK 0
+
+
+namespace Geom {
+
+namespace detail { namespace bezier_clipping {
+
+////////////////////////////////////////////////////////////////////////////////
+// for debugging
+//
+
+inline
+void print(std::vector<Point> const& cp, const char* msg = "")
+{
+ std::cerr << msg << std::endl;
+ for (size_t i = 0; i < cp.size(); ++i)
+ std::cerr << i << " : " << cp[i] << std::endl;
+}
+
+template< class charT >
+inline
+std::basic_ostream<charT> &
+operator<< (std::basic_ostream<charT> & os, const Interval & I)
+{
+ os << "[" << I.min() << ", " << I.max() << "]";
+ return os;
+}
+
+inline
+double angle (std::vector<Point> const& A)
+{
+ size_t n = A.size() -1;
+ double a = std::atan2(A[n][Y] - A[0][Y], A[n][X] - A[0][X]);
+ return (180 * a / M_PI);
+}
+
+inline
+size_t get_precision(Interval const& I)
+{
+ double d = I.extent();
+ double e = 0.1, p = 10;
+ int n = 0;
+ while (n < 16 && d < e)
+ {
+ p *= 10;
+ e = 1/p;
+ ++n;
+ }
+ return n;
+}
+
+inline
+void range_assertion(int k, int m, int n, const char* msg)
+{
+ if ( k < m || k > n)
+ {
+ std::cerr << "range assertion failed: \n"
+ << msg << std::endl
+ << "value: " << k
+ << " range: " << m << ", " << n << std::endl;
+ assert (k >= m && k <= n);
+ }
+}
+
+
+////////////////////////////////////////////////////////////////////////////////
+// convex hull
+
+/*
+ * return true in case the oriented polyline p0, p1, p2 is a right turn
+ */
+inline
+bool is_a_right_turn (Point const& p0, Point const& p1, Point const& p2)
+{
+ if (p1 == p2) return false;
+ Point q1 = p1 - p0;
+ Point q2 = p2 - p0;
+ if (q1 == -q2) return false;
+ return (cross (q1, q2) < 0);
+}
+
+/*
+ * return true if p < q wrt the lexicographyc order induced by the coordinates
+ */
+struct lex_less
+{
+ bool operator() (Point const& p, Point const& q)
+ {
+ return ((p[X] < q[X]) || (p[X] == q[X] && p[Y] < q[Y]));
+ }
+};
+
+/*
+ * return true if p > q wrt the lexicographyc order induced by the coordinates
+ */
+struct lex_greater
+{
+ bool operator() (Point const& p, Point const& q)
+ {
+ return ((p[X] > q[X]) || (p[X] == q[X] && p[Y] > q[Y]));
+ }
+};
+
+/*
+ * Compute the convex hull of a set of points.
+ * The implementation is based on the Andrew's scan algorithm
+ * note: in the Bezier clipping for collinear normals it seems
+ * to be more stable wrt the Graham's scan algorithm and in general
+ * a bit quikier
+ */
+void convex_hull (std::vector<Point> & P)
+{
+ size_t n = P.size();
+ if (n < 2) return;
+ std::sort(P.begin(), P.end(), lex_less());
+ if (n < 4) return;
+ // upper hull
+ size_t u = 2;
+ for (size_t i = 2; i < n; ++i)
+ {
+ while (u > 1 && !is_a_right_turn(P[u-2], P[u-1], P[i]))
+ {
+ --u;
+ }
+ std::swap(P[u], P[i]);
+ ++u;
+ }
+ std::sort(P.begin() + u, P.end(), lex_greater());
+ std::rotate(P.begin(), P.begin() + 1, P.end());
+ // lower hull
+ size_t l = u;
+ size_t k = u - 1;
+ for (size_t i = l; i < n; ++i)
+ {
+ while (l > k && !is_a_right_turn(P[l-2], P[l-1], P[i]))
+ {
+ --l;
+ }
+ std::swap(P[l], P[i]);
+ ++l;
+ }
+ P.resize(l);
+}
+
+
+////////////////////////////////////////////////////////////////////////////////
+// numerical routines
+
+/*
+ * Compute the binomial coefficient (n, k)
+ */
+inline
+double binomial(unsigned int n, unsigned int k)
+{
+ return choose<double>(n, k);
+}
+
+/*
+ * Compute the determinant of the 2x2 matrix with column the point P1, P2
+ */
+inline
+double det(Point const& P1, Point const& P2)
+{
+ return P1[X]*P2[Y] - P1[Y]*P2[X];
+}
+
+/*
+ * Solve the linear system [P1,P2] * P = Q
+ * in case there isn't exactly one solution the routine returns false
+ */
+inline
+bool solve(Point & P, Point const& P1, Point const& P2, Point const& Q)
+{
+ double d = det(P1, P2);
+ if (d == 0) return false;
+ d = 1 / d;
+ P[X] = det(Q, P2) * d;
+ P[Y] = det(P1, Q) * d;
+ return true;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+// interval routines
+
+/*
+ * Map the sub-interval I in [0,1] into the interval J and assign it to J
+ */
+inline
+void map_to(Interval & J, Interval const& I)
+{
+ double length = J.extent();
+ J[1] = I.max() * length + J[0];
+ J[0] = I.min() * length + J[0];
+}
+
+/*
+ * The interval [1,0] is used to represent the empty interval, this routine
+ * is just an helper function for creating such an interval
+ */
+inline
+Interval make_empty_interval()
+{
+ Interval I(0);
+ I[0] = 1;
+ return I;
+}
+
+
+////////////////////////////////////////////////////////////////////////////////
+// bezier curve routines
+
+/*
+ * Return true if all the Bezier curve control points are near,
+ * false otherwise
+ */
+inline
+bool is_constant(std::vector<Point> const& A, double precision = EPSILON)
+{
+ for (unsigned int i = 1; i < A.size(); ++i)
+ {
+ if(!are_near(A[i], A[0], precision))
+ return false;
+ }
+ return true;
+}
+
+/*
+ * Compute the hodograph of the bezier curve B and return it in D
+ */
+inline
+void derivative(std::vector<Point> & D, std::vector<Point> const& B)
+{
+ D.clear();
+ size_t sz = B.size();
+ if (sz == 0) return;
+ if (sz == 1)
+ {
+ D.resize(1, Point(0,0));
+ return;
+ }
+ size_t n = sz-1;
+ D.reserve(n);
+ for (size_t i = 0; i < n; ++i)
+ {
+ D.push_back(n*(B[i+1] - B[i]));
+ }
+}
+
+/*
+ * Compute the hodograph of the Bezier curve B rotated of 90 degree
+ * and return it in D; we have N(t) orthogonal to B(t) for any t
+ */
+inline
+void normal(std::vector<Point> & N, std::vector<Point> const& B)
+{
+ derivative(N,B);
+ for (size_t i = 0; i < N.size(); ++i)
+ {
+ N[i] = rot90(N[i]);
+ }
+}
+
+/*
+ * Compute the portion of the Bezier curve "B" wrt the interval [0,t]
+ */
+inline
+void left_portion(Coord t, std::vector<Point> & B)
+{
+ size_t n = B.size();
+ for (size_t i = 1; i < n; ++i)
+ {
+ for (size_t j = n-1; j > i-1 ; --j)
+ {
+ B[j] = lerp(t, B[j-1], B[j]);
+ }
+ }
+}
+
+/*
+ * Compute the portion of the Bezier curve "B" wrt the interval [t,1]
+ */
+inline
+void right_portion(Coord t, std::vector<Point> & B)
+{
+ size_t n = B.size();
+ for (size_t i = 1; i < n; ++i)
+ {
+ for (size_t j = 0; j < n-i; ++j)
+ {
+ B[j] = lerp(t, B[j], B[j+1]);
+ }
+ }
+}
+
+/*
+ * Compute the portion of the Bezier curve "B" wrt the interval "I"
+ */
+inline
+void portion (std::vector<Point> & B , Interval const& I)
+{
+ if (I.min() == 0)
+ {
+ if (I.max() == 1) return;
+ left_portion(I.max(), B);
+ return;
+ }
+ right_portion(I.min(), B);
+ if (I.max() == 1) return;
+ double t = I.extent() / (1 - I.min());
+ left_portion(t, B);
+}
+
+
+////////////////////////////////////////////////////////////////////////////////
+// tags
+
+struct intersection_point_tag;
+struct collinear_normal_tag;
+template <typename Tag>
+void clip(Interval & dom,
+ std::vector<Point> const& A,
+ std::vector<Point> const& B);
+template <typename Tag>
+void iterate(std::vector<Interval>& domsA,
+ std::vector<Interval>& domsB,
+ std::vector<Point> const& A,
+ std::vector<Point> const& B,
+ Interval const& domA,
+ Interval const& domB,
+ double precision );
+
+
+////////////////////////////////////////////////////////////////////////////////
+// intersection
+
+/*
+ * Make up an orientation line using the control points c[i] and c[j]
+ * the line is returned in the output parameter "l" in the form of a 3 element
+ * vector : l[0] * x + l[1] * y + l[2] == 0; the line is normalized.
+ */
+inline
+void orientation_line (std::vector<double> & l,
+ std::vector<Point> const& c,
+ size_t i, size_t j)
+{
+ l[0] = c[j][Y] - c[i][Y];
+ l[1] = c[i][X] - c[j][X];
+ l[2] = cross(c[i], c[j]);
+ double length = std::sqrt(l[0] * l[0] + l[1] * l[1]);
+ assert (length != 0);
+ l[0] /= length;
+ l[1] /= length;
+ l[2] /= length;
+}
+
+/*
+ * Pick up an orientation line for the Bezier curve "c" and return it in
+ * the output parameter "l"
+ */
+inline
+void pick_orientation_line (std::vector<double> & l,
+ std::vector<Point> const& c)
+{
+ size_t i = c.size();
+ while (--i > 0 && are_near(c[0], c[i]))
+ {}
+ if (i == 0)
+ {
+ // this should never happen because when a new curve portion is created
+ // we check that it is not constant;
+ // however this requires that the precision used in the is_constant
+ // routine has to be the same used here in the are_near test
+ assert(i != 0);
+ }
+ orientation_line(l, c, 0, i);
+ //std::cerr << "i = " << i << std::endl;
+}
+
+/*
+ * Make up an orientation line for constant bezier curve;
+ * the orientation line is made up orthogonal to the other curve base line;
+ * the line is returned in the output parameter "l" in the form of a 3 element
+ * vector : l[0] * x + l[1] * y + l[2] == 0; the line is normalized.
+ */
+inline
+void orthogonal_orientation_line (std::vector<double> & l,
+ std::vector<Point> const& c,
+ Point const& p)
+{
+ if (is_constant(c))
+ {
+ // this should never happen
+ assert(!is_constant(c));
+ }
+ std::vector<Point> ol(2);
+ ol[0] = p;
+ ol[1] = (c.back() - c.front()).cw() + p;
+ orientation_line(l, ol, 0, 1);
+}
+
+/*
+ * Compute the signed distance of the point "P" from the normalized line l
+ */
+inline
+double distance (Point const& P, std::vector<double> const& l)
+{
+ return l[X] * P[X] + l[Y] * P[Y] + l[2];
+}
+
+/*
+ * Compute the min and max distance of the control points of the Bezier
+ * curve "c" from the normalized orientation line "l".
+ * This bounds are returned through the output Interval parameter"bound".
+ */
+inline
+void fat_line_bounds (Interval& bound,
+ std::vector<Point> const& c,
+ std::vector<double> const& l)
+{
+ bound[0] = 0;
+ bound[1] = 0;
+ double d;
+ for (size_t i = 0; i < c.size(); ++i)
+ {
+ d = distance(c[i], l);
+ if (bound[0] > d) bound[0] = d;
+ if (bound[1] < d) bound[1] = d;
+ }
+}
+
+/*
+ * return the x component of the intersection point between the line
+ * passing through points p1, p2 and the line Y = "y"
+ */
+inline
+double intersect (Point const& p1, Point const& p2, double y)
+{
+ // we are sure that p2[Y] != p1[Y] because this routine is called
+ // only when the lower or the upper bound is crossed
+ double dy = (p2[Y] - p1[Y]);
+ double s = (y - p1[Y]) / dy;
+ return (p2[X]-p1[X])*s + p1[X];
+}
+
+/*
+ * Clip the Bezier curve "B" wrt the fat line defined by the orientation
+ * line "l" and the interval range "bound", the new parameter interval for
+ * the clipped curve is returned through the output parameter "dom"
+ */
+void clip_interval (Interval& dom,
+ std::vector<Point> const& B,
+ std::vector<double> const& l,
+ Interval const& bound)
+{
+ double n = B.size() - 1; // number of sub-intervals
+ std::vector<Point> D; // distance curve control points
+ D.reserve (B.size());
+ double d;
+ for (size_t i = 0; i < B.size(); ++i)
+ {
+ d = distance (B[i], l);
+ D.push_back (Point(i/n, d));
+ }
+ //print(D);
+
+ convex_hull(D);
+ std::vector<Point> & p = D;
+ //print(p);
+
+ bool plower, phigher;
+ bool clower, chigher;
+ double t, tmin = 1, tmax = 0;
+// std::cerr << "bound : " << bound << std::endl;
+
+ plower = (p[0][Y] < bound.min());
+ phigher = (p[0][Y] > bound.max());
+ if (!(plower || phigher)) // inside the fat line
+ {
+ if (tmin > p[0][X]) tmin = p[0][X];
+ if (tmax < p[0][X]) tmax = p[0][X];
+// std::cerr << "0 : inside " << p[0]
+// << " : tmin = " << tmin << ", tmax = " << tmax << std::endl;
+ }
+
+ for (size_t i = 1; i < p.size(); ++i)
+ {
+ clower = (p[i][Y] < bound.min());
+ chigher = (p[i][Y] > bound.max());
+ if (!(clower || chigher)) // inside the fat line
+ {
+ if (tmin > p[i][X]) tmin = p[i][X];
+ if (tmax < p[i][X]) tmax = p[i][X];
+// std::cerr << i << " : inside " << p[i]
+// << " : tmin = " << tmin << ", tmax = " << tmax
+// << std::endl;
+ }
+ if (clower != plower) // cross the lower bound
+ {
+ t = intersect(p[i-1], p[i], bound.min());
+ if (tmin > t) tmin = t;
+ if (tmax < t) tmax = t;
+ plower = clower;
+// std::cerr << i << " : lower " << p[i]
+// << " : tmin = " << tmin << ", tmax = " << tmax
+// << std::endl;
+ }
+ if (chigher != phigher) // cross the upper bound
+ {
+ t = intersect(p[i-1], p[i], bound.max());
+ if (tmin > t) tmin = t;
+ if (tmax < t) tmax = t;
+ phigher = chigher;
+// std::cerr << i << " : higher " << p[i]
+// << " : tmin = " << tmin << ", tmax = " << tmax
+// << std::endl;
+ }
+ }
+
+ // we have to test the closing segment for intersection
+ size_t last = p.size() - 1;
+ clower = (p[0][Y] < bound.min());
+ chigher = (p[0][Y] > bound.max());
+ if (clower != plower) // cross the lower bound
+ {
+ t = intersect(p[last], p[0], bound.min());
+ if (tmin > t) tmin = t;
+ if (tmax < t) tmax = t;
+// std::cerr << "0 : lower " << p[0]
+// << " : tmin = " << tmin << ", tmax = " << tmax << std::endl;
+ }
+ if (chigher != phigher) // cross the upper bound
+ {
+ t = intersect(p[last], p[0], bound.max());
+ if (tmin > t) tmin = t;
+ if (tmax < t) tmax = t;
+// std::cerr << "0 : higher " << p[0]
+// << " : tmin = " << tmin << ", tmax = " << tmax << std::endl;
+ }
+
+ dom[0] = tmin;
+ dom[1] = tmax;
+}
+
+/*
+ * Clip the Bezier curve "B" wrt the Bezier curve "A" for individuating
+ * intersection points the new parameter interval for the clipped curve
+ * is returned through the output parameter "dom"
+ */
+template <>
+inline
+void clip<intersection_point_tag> (Interval & dom,
+ std::vector<Point> const& A,
+ std::vector<Point> const& B)
+{
+ std::vector<double> bl(3);
+ Interval bound;
+ if (is_constant(A))
+ {
+ Point M = middle_point(A.front(), A.back());
+ orthogonal_orientation_line(bl, B, M);
+ }
+ else
+ {
+ pick_orientation_line(bl, A);
+ }
+ fat_line_bounds(bound, A, bl);
+ clip_interval(dom, B, bl, bound);
+}
+
+
+///////////////////////////////////////////////////////////////////////////////
+// collinear normal
+
+/*
+ * Compute a closed focus for the Bezier curve B and return it in F
+ * A focus is any curve through which all lines perpendicular to B(t) pass.
+ */
+inline
+void make_focus (std::vector<Point> & F, std::vector<Point> const& B)
+{
+ assert (B.size() > 2);
+ size_t n = B.size() - 1;
+ normal(F, B);
+ Point c(1, 1);
+#if VERBOSE
+ if (!solve(c, F[0], -F[n-1], B[n]-B[0]))
+ {
+ std::cerr << "make_focus: unable to make up a closed focus" << std::endl;
+ }
+#else
+ solve(c, F[0], -F[n-1], B[n]-B[0]);
+#endif
+// std::cerr << "c = " << c << std::endl;
+
+
+ // B(t) + c(t) * N(t)
+ double n_inv = 1 / (double)(n);
+ Point c0ni;
+ F.push_back(c[1] * F[n-1]);
+ F[n] += B[n];
+ for (size_t i = n-1; i > 0; --i)
+ {
+ F[i] *= -c[0];
+ c0ni = F[i];
+ F[i] += (c[1] * F[i-1]);
+ F[i] *= (i * n_inv);
+ F[i] -= c0ni;
+ F[i] += B[i];
+ }
+ F[0] *= c[0];
+ F[0] += B[0];
+}
+
+/*
+ * Compute the projection on the plane (t, d) of the control points
+ * (t, u, D(t,u)) where D(t,u) = <(B(t) - F(u)), B'(t)> with 0 <= t, u <= 1
+ * B is a Bezier curve and F is a focus of another Bezier curve.
+ * See Sederberg, Nishita, 1990 - Curve intersection using Bezier clipping.
+ */
+void distance_control_points (std::vector<Point> & D,
+ std::vector<Point> const& B,
+ std::vector<Point> const& F)
+{
+ assert (B.size() > 1);
+ assert (F.size() > 0);
+ const size_t n = B.size() - 1;
+ const size_t m = F.size() - 1;
+ const size_t r = 2 * n - 1;
+ const double r_inv = 1 / (double)(r);
+ D.clear();
+ D.reserve (B.size() * F.size());
+
+ std::vector<Point> dB;
+ dB.reserve(n);
+ for (size_t k = 0; k < n; ++k)
+ {
+ dB.push_back (B[k+1] - B[k]);
+ }
+ NL::Matrix dBB(n,B.size());
+ for (size_t i = 0; i < n; ++i)
+ for (size_t j = 0; j < B.size(); ++j)
+ dBB(i,j) = dot (dB[i], B[j]);
+ NL::Matrix dBF(n, F.size());
+ for (size_t i = 0; i < n; ++i)
+ for (size_t j = 0; j < F.size(); ++j)
+ dBF(i,j) = dot (dB[i], F[j]);
+
+ size_t k0, kn, l;
+ double bc, bri;
+ Point dij;
+ std::vector<double> d(F.size());
+ for (size_t i = 0; i <= r; ++i)
+ {
+ for (size_t j = 0; j <= m; ++j)
+ {
+ d[j] = 0;
+ }
+ k0 = std::max(i, n) - n;
+ kn = std::min(i, n-1);
+ bri = n / binomial(r,i);
+ for (size_t k = k0; k <= kn; ++k)
+ {
+ //if (k > i || (i-k) > n) continue;
+ l = i - k;
+#if CHECK
+ assert (l <= n);
+#endif
+ bc = bri * binomial(n,l) * binomial(n-1, k);
+ for (size_t j = 0; j <= m; ++j)
+ {
+ //d[j] += bc * dot(dB[k], B[l] - F[j]);
+ d[j] += bc * (dBB(k,l) - dBF(k,j));
+ }
+ }
+ double dmin, dmax;
+ dmin = dmax = d[m];
+ for (size_t j = 0; j < m; ++j)
+ {
+ if (dmin > d[j]) dmin = d[j];
+ if (dmax < d[j]) dmax = d[j];
+ }
+ dij[0] = i * r_inv;
+ dij[1] = dmin;
+ D.push_back (dij);
+ dij[1] = dmax;
+ D.push_back (dij);
+ }
+}
+
+/*
+ * Clip the Bezier curve "B" wrt the focus "F"; the new parameter interval for
+ * the clipped curve is returned through the output parameter "dom"
+ */
+void clip_interval (Interval& dom,
+ std::vector<Point> const& B,
+ std::vector<Point> const& F)
+{
+ std::vector<Point> D; // distance curve control points
+ distance_control_points(D, B, F);
+ //print(D, "D");
+// ConvexHull chD(D);
+// std::vector<Point>& p = chD.boundary; // convex hull vertices
+
+ convex_hull(D);
+ std::vector<Point> & p = D;
+ //print(p, "CH(D)");
+
+ bool plower, clower;
+ double t, tmin = 1, tmax = 0;
+
+ plower = (p[0][Y] < 0);
+ if (p[0][Y] == 0) // on the x axis
+ {
+ if (tmin > p[0][X]) tmin = p[0][X];
+ if (tmax < p[0][X]) tmax = p[0][X];
+// std::cerr << "0 : on x axis " << p[0]
+// << " : tmin = " << tmin << ", tmax = " << tmax << std::endl;
+ }
+
+ for (size_t i = 1; i < p.size(); ++i)
+ {
+ clower = (p[i][Y] < 0);
+ if (p[i][Y] == 0) // on x axis
+ {
+ if (tmin > p[i][X]) tmin = p[i][X];
+ if (tmax < p[i][X]) tmax = p[i][X];
+// std::cerr << i << " : on x axis " << p[i]
+// << " : tmin = " << tmin << ", tmax = " << tmax
+// << std::endl;
+ }
+ else if (clower != plower) // cross the x axis
+ {
+ t = intersect(p[i-1], p[i], 0);
+ if (tmin > t) tmin = t;
+ if (tmax < t) tmax = t;
+ plower = clower;
+// std::cerr << i << " : lower " << p[i]
+// << " : tmin = " << tmin << ", tmax = " << tmax
+// << std::endl;
+ }
+ }
+
+ // we have to test the closing segment for intersection
+ size_t last = p.size() - 1;
+ clower = (p[0][Y] < 0);
+ if (clower != plower) // cross the x axis
+ {
+ t = intersect(p[last], p[0], 0);
+ if (tmin > t) tmin = t;
+ if (tmax < t) tmax = t;
+// std::cerr << "0 : lower " << p[0]
+// << " : tmin = " << tmin << ", tmax = " << tmax << std::endl;
+ }
+ dom[0] = tmin;
+ dom[1] = tmax;
+}
+
+/*
+ * Clip the Bezier curve "B" wrt the Bezier curve "A" for individuating
+ * points which have collinear normals; the new parameter interval
+ * for the clipped curve is returned through the output parameter "dom"
+ */
+template <>
+inline
+void clip<collinear_normal_tag> (Interval & dom,
+ std::vector<Point> const& A,
+ std::vector<Point> const& B)
+{
+ std::vector<Point> F;
+ make_focus(F, A);
+ clip_interval(dom, B, F);
+}
+
+
+
+const double MAX_PRECISION = 1e-8;
+const double MIN_CLIPPED_SIZE_THRESHOLD = 0.8;
+const Interval UNIT_INTERVAL(0,1);
+const Interval EMPTY_INTERVAL = make_empty_interval();
+const Interval H1_INTERVAL(0, 0.5);
+const Interval H2_INTERVAL(0.5 + MAX_PRECISION, 1.0);
+
+/*
+ * iterate
+ *
+ * input:
+ * A, B: control point sets of two bezier curves
+ * domA, domB: real parameter intervals of the two curves
+ * precision: required computational precision of the returned parameter ranges
+ * output:
+ * domsA, domsB: sets of parameter intervals
+ *
+ * The parameter intervals are computed by using a Bezier clipping algorithm,
+ * in case the clipping doesn't shrink the initial interval more than 20%,
+ * a subdivision step is performed.
+ * If during the computation both curves collapse to a single point
+ * the routine exits indipendently by the precision reached in the computation
+ * of the curve intervals.
+ */
+template <>
+void iterate<intersection_point_tag> (std::vector<Interval>& domsA,
+ std::vector<Interval>& domsB,
+ std::vector<Point> const& A,
+ std::vector<Point> const& B,
+ Interval const& domA,
+ Interval const& domB,
+ double precision )
+{
+ // in order to limit recursion
+ static size_t counter = 0;
+ if (domA.extent() == 1 && domB.extent() == 1) counter = 0;
+ if (++counter > 100) return;
+#if VERBOSE
+ std::cerr << std::fixed << std::setprecision(16);
+ std::cerr << ">> curve subdision performed <<" << std::endl;
+ std::cerr << "dom(A) : " << domA << std::endl;
+ std::cerr << "dom(B) : " << domB << std::endl;
+// std::cerr << "angle(A) : " << angle(A) << std::endl;
+// std::cerr << "angle(B) : " << angle(B) << std::endl;
+#endif
+
+ if (precision < MAX_PRECISION)
+ precision = MAX_PRECISION;
+
+ std::vector<Point> pA = A;
+ std::vector<Point> pB = B;
+ std::vector<Point>* C1 = &pA;
+ std::vector<Point>* C2 = &pB;
+
+ Interval dompA = domA;
+ Interval dompB = domB;
+ Interval* dom1 = &dompA;
+ Interval* dom2 = &dompB;
+
+ Interval dom;
+
+ if ( is_constant(A) && is_constant(B) ){
+ Point M1 = middle_point(C1->front(), C1->back());
+ Point M2 = middle_point(C2->front(), C2->back());
+ if (are_near(M1,M2)){
+ domsA.push_back(domA);
+ domsB.push_back(domB);
+ }
+ return;
+ }
+
+ size_t iter = 0;
+ while (++iter < 100
+ && (dompA.extent() >= precision || dompB.extent() >= precision))
+ {
+#if VERBOSE
+ std::cerr << "iter: " << iter << std::endl;
+#endif
+ clip<intersection_point_tag>(dom, *C1, *C2);
+
+ // [1,0] is utilized to represent an empty interval
+ if (dom == EMPTY_INTERVAL)
+ {
+#if VERBOSE
+ std::cerr << "dom: empty" << std::endl;
+#endif
+ return;
+ }
+#if VERBOSE
+ std::cerr << "dom : " << dom << std::endl;
+#endif
+ // all other cases where dom[0] > dom[1] are invalid
+ if (dom.min() > dom.max())
+ {
+ assert(dom.min() < dom.max());
+ }
+
+ map_to(*dom2, dom);
+
+ portion(*C2, dom);
+ if (is_constant(*C2) && is_constant(*C1))
+ {
+ Point M1 = middle_point(C1->front(), C1->back());
+ Point M2 = middle_point(C2->front(), C2->back());
+#if VERBOSE
+ std::cerr << "both curves are constant: \n"
+ << "M1: " << M1 << "\n"
+ << "M2: " << M2 << std::endl;
+ print(*C2, "C2");
+ print(*C1, "C1");
+#endif
+ if (are_near(M1,M2))
+ break; // append the new interval
+ else
+ return; // exit without appending any new interval
+ }
+
+
+ // if we have clipped less than 20% than we need to subdive the curve
+ // with the largest domain into two sub-curves
+ if ( dom.extent() > MIN_CLIPPED_SIZE_THRESHOLD)
+ {
+#if VERBOSE
+ std::cerr << "clipped less than 20% : " << dom.extent() << std::endl;
+ std::cerr << "angle(pA) : " << angle(pA) << std::endl;
+ std::cerr << "angle(pB) : " << angle(pB) << std::endl;
+#endif
+ std::vector<Point> pC1, pC2;
+ Interval dompC1, dompC2;
+ if (dompA.extent() > dompB.extent())
+ {
+ pC1 = pC2 = pA;
+ portion(pC1, H1_INTERVAL);
+ portion(pC2, H2_INTERVAL);
+ dompC1 = dompC2 = dompA;
+ map_to(dompC1, H1_INTERVAL);
+ map_to(dompC2, H2_INTERVAL);
+ iterate<intersection_point_tag>(domsA, domsB, pC1, pB,
+ dompC1, dompB, precision);
+ iterate<intersection_point_tag>(domsA, domsB, pC2, pB,
+ dompC2, dompB, precision);
+ }
+ else
+ {
+ pC1 = pC2 = pB;
+ portion(pC1, H1_INTERVAL);
+ portion(pC2, H2_INTERVAL);
+ dompC1 = dompC2 = dompB;
+ map_to(dompC1, H1_INTERVAL);
+ map_to(dompC2, H2_INTERVAL);
+ iterate<intersection_point_tag>(domsB, domsA, pC1, pA,
+ dompC1, dompA, precision);
+ iterate<intersection_point_tag>(domsB, domsA, pC2, pA,
+ dompC2, dompA, precision);
+ }
+ return;
+ }
+
+ std::swap(C1, C2);
+ std::swap(dom1, dom2);
+#if VERBOSE
+ std::cerr << "dom(pA) : " << dompA << std::endl;
+ std::cerr << "dom(pB) : " << dompB << std::endl;
+#endif
+ }
+ domsA.push_back(dompA);
+ domsB.push_back(dompB);
+}
+
+
+/*
+ * iterate
+ *
+ * input:
+ * A, B: control point sets of two bezier curves
+ * domA, domB: real parameter intervals of the two curves
+ * precision: required computational precision of the returned parameter ranges
+ * output:
+ * domsA, domsB: sets of parameter intervals
+ *
+ * The parameter intervals are computed by using a Bezier clipping algorithm,
+ * in case the clipping doesn't shrink the initial interval more than 20%,
+ * a subdivision step is performed.
+ * If during the computation one of the two curve interval length becomes less
+ * than MAX_PRECISION the routine exits indipendently by the precision reached
+ * in the computation of the other curve interval.
+ */
+template <>
+void iterate<collinear_normal_tag> (std::vector<Interval>& domsA,
+ std::vector<Interval>& domsB,
+ std::vector<Point> const& A,
+ std::vector<Point> const& B,
+ Interval const& domA,
+ Interval const& domB,
+ double precision)
+{
+ // in order to limit recursion
+ static size_t counter = 0;
+ if (domA.extent() == 1 && domB.extent() == 1) counter = 0;
+ if (++counter > 100) return;
+#if VERBOSE
+ std::cerr << std::fixed << std::setprecision(16);
+ std::cerr << ">> curve subdision performed <<" << std::endl;
+ std::cerr << "dom(A) : " << domA << std::endl;
+ std::cerr << "dom(B) : " << domB << std::endl;
+// std::cerr << "angle(A) : " << angle(A) << std::endl;
+// std::cerr << "angle(B) : " << angle(B) << std::endl;
+#endif
+
+ if (precision < MAX_PRECISION)
+ precision = MAX_PRECISION;
+
+ std::vector<Point> pA = A;
+ std::vector<Point> pB = B;
+ std::vector<Point>* C1 = &pA;
+ std::vector<Point>* C2 = &pB;
+
+ Interval dompA = domA;
+ Interval dompB = domB;
+ Interval* dom1 = &dompA;
+ Interval* dom2 = &dompB;
+
+ Interval dom;
+
+ size_t iter = 0;
+ while (++iter < 100
+ && (dompA.extent() >= precision || dompB.extent() >= precision))
+ {
+#if VERBOSE
+ std::cerr << "iter: " << iter << std::endl;
+#endif
+ clip<collinear_normal_tag>(dom, *C1, *C2);
+
+ // [1,0] is utilized to represent an empty interval
+ if (dom == EMPTY_INTERVAL)
+ {
+#if VERBOSE
+ std::cerr << "dom: empty" << std::endl;
+#endif
+ return;
+ }
+#if VERBOSE
+ std::cerr << "dom : " << dom << std::endl;
+#endif
+ // all other cases where dom[0] > dom[1] are invalid
+ if (dom.min() > dom.max())
+ {
+ assert(dom.min() < dom.max());
+ }
+
+ map_to(*dom2, dom);
+
+ // it's better to stop before losing computational precision
+ if (iter > 1 && (dom2->extent() <= MAX_PRECISION))
+ {
+#if VERBOSE
+ std::cerr << "beyond max precision limit" << std::endl;
+#endif
+ break;
+ }
+
+ portion(*C2, dom);
+ if (iter > 1 && is_constant(*C2))
+ {
+#if VERBOSE
+ std::cerr << "new curve portion pC1 is constant" << std::endl;
+#endif
+ break;
+ }
+
+
+ // if we have clipped less than 20% than we need to subdive the curve
+ // with the largest domain into two sub-curves
+ if ( dom.extent() > MIN_CLIPPED_SIZE_THRESHOLD)
+ {
+#if VERBOSE
+ std::cerr << "clipped less than 20% : " << dom.extent() << std::endl;
+ std::cerr << "angle(pA) : " << angle(pA) << std::endl;
+ std::cerr << "angle(pB) : " << angle(pB) << std::endl;
+#endif
+ std::vector<Point> pC1, pC2;
+ Interval dompC1, dompC2;
+ if (dompA.extent() > dompB.extent())
+ {
+ if ((dompA.extent() / 2) < MAX_PRECISION)
+ {
+ break;
+ }
+ pC1 = pC2 = pA;
+ portion(pC1, H1_INTERVAL);
+ if (false && is_constant(pC1))
+ {
+#if VERBOSE
+ std::cerr << "new curve portion pC1 is constant" << std::endl;
+#endif
+ break;
+ }
+ portion(pC2, H2_INTERVAL);
+ if (is_constant(pC2))
+ {
+#if VERBOSE
+ std::cerr << "new curve portion pC2 is constant" << std::endl;
+#endif
+ break;
+ }
+ dompC1 = dompC2 = dompA;
+ map_to(dompC1, H1_INTERVAL);
+ map_to(dompC2, H2_INTERVAL);
+ iterate<collinear_normal_tag>(domsA, domsB, pC1, pB,
+ dompC1, dompB, precision);
+ iterate<collinear_normal_tag>(domsA, domsB, pC2, pB,
+ dompC2, dompB, precision);
+ }
+ else
+ {
+ if ((dompB.extent() / 2) < MAX_PRECISION)
+ {
+ break;
+ }
+ pC1 = pC2 = pB;
+ portion(pC1, H1_INTERVAL);
+ if (is_constant(pC1))
+ {
+#if VERBOSE
+ std::cerr << "new curve portion pC1 is constant" << std::endl;
+#endif
+ break;
+ }
+ portion(pC2, H2_INTERVAL);
+ if (is_constant(pC2))
+ {
+#if VERBOSE
+ std::cerr << "new curve portion pC2 is constant" << std::endl;
+#endif
+ break;
+ }
+ dompC1 = dompC2 = dompB;
+ map_to(dompC1, H1_INTERVAL);
+ map_to(dompC2, H2_INTERVAL);
+ iterate<collinear_normal_tag>(domsB, domsA, pC1, pA,
+ dompC1, dompA, precision);
+ iterate<collinear_normal_tag>(domsB, domsA, pC2, pA,
+ dompC2, dompA, precision);
+ }
+ return;
+ }
+
+ std::swap(C1, C2);
+ std::swap(dom1, dom2);
+#if VERBOSE
+ std::cerr << "dom(pA) : " << dompA << std::endl;
+ std::cerr << "dom(pB) : " << dompB << std::endl;
+#endif
+ }
+ domsA.push_back(dompA);
+ domsB.push_back(dompB);
+}
+
+
+/*
+ * get_solutions
+ *
+ * input: A, B - set of control points of two Bezier curve
+ * input: precision - required precision of computation
+ * input: clip - the routine used for clipping
+ * output: xs - set of pairs of parameter values
+ * at which the clipping algorithm converges
+ *
+ * This routine is based on the Bezier Clipping Algorithm,
+ * see: Sederberg - Computer Aided Geometric Design
+ */
+template <typename Tag>
+void get_solutions (std::vector< std::pair<double, double> >& xs,
+ std::vector<Point> const& A,
+ std::vector<Point> const& B,
+ double precision)
+{
+ std::pair<double, double> ci;
+ std::vector<Interval> domsA, domsB;
+ iterate<Tag> (domsA, domsB, A, B, UNIT_INTERVAL, UNIT_INTERVAL, precision);
+ if (domsA.size() != domsB.size())
+ {
+ assert (domsA.size() == domsB.size());
+ }
+ xs.clear();
+ xs.reserve(domsA.size());
+ for (size_t i = 0; i < domsA.size(); ++i)
+ {
+#if VERBOSE
+ std::cerr << i << " : domA : " << domsA[i] << std::endl;
+ std::cerr << "extent A: " << domsA[i].extent() << " ";
+ std::cerr << "precision A: " << get_precision(domsA[i]) << std::endl;
+ std::cerr << i << " : domB : " << domsB[i] << std::endl;
+ std::cerr << "extent B: " << domsB[i].extent() << " ";
+ std::cerr << "precision B: " << get_precision(domsB[i]) << std::endl;
+#endif
+ ci.first = domsA[i].middle();
+ ci.second = domsB[i].middle();
+ xs.push_back(ci);
+ }
+}
+
+} /* end namespace bezier_clipping */ } /* end namespace detail */
+
+
+/*
+ * find_collinear_normal
+ *
+ * input: A, B - set of control points of two Bezier curve
+ * input: precision - required precision of computation
+ * output: xs - set of pairs of parameter values
+ * at which there are collinear normals
+ *
+ * This routine is based on the Bezier Clipping Algorithm,
+ * see: Sederberg, Nishita, 1990 - Curve intersection using Bezier clipping
+ */
+void find_collinear_normal (std::vector< std::pair<double, double> >& xs,
+ std::vector<Point> const& A,
+ std::vector<Point> const& B,
+ double precision)
+{
+ using detail::bezier_clipping::get_solutions;
+ using detail::bezier_clipping::collinear_normal_tag;
+ get_solutions<collinear_normal_tag>(xs, A, B, precision);
+}
+
+
+/*
+ * find_intersections_bezier_clipping
+ *
+ * input: A, B - set of control points of two Bezier curve
+ * input: precision - required precision of computation
+ * output: xs - set of pairs of parameter values
+ * at which crossing happens
+ *
+ * This routine is based on the Bezier Clipping Algorithm,
+ * see: Sederberg, Nishita, 1990 - Curve intersection using Bezier clipping
+ */
+void find_intersections_bezier_clipping (std::vector< std::pair<double, double> >& xs,
+ std::vector<Point> const& A,
+ std::vector<Point> const& B,
+ double precision)
+{
+ using detail::bezier_clipping::get_solutions;
+ using detail::bezier_clipping::intersection_point_tag;
+ get_solutions<intersection_point_tag>(xs, A, B, precision);
+}
+
+} // 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/point.cpp b/src/2geom/point.cpp
index b8a858f8d48d6f9910441a93100cef69a9a8a743..4a0625713385337d9e211cbbf627fb49f4d3a45b 100644 (file)
--- a/src/2geom/point.cpp
+++ b/src/2geom/point.cpp
return *this;
}
-Point constrain_angle(Point const &ref, Point const &pt, unsigned int n, Point const &dir)
+Point constrain_angle(Point const &A, Point const &B, unsigned int n, Point const &dir)
{
- // for special cases we could perhaps use faster routines
+ // for special cases we could perhaps use explicit testing (which might be faster)
if (n == 0.0) {
- return pt;
+ return B;
}
- Point diff(pt - ref);
+ Point diff(B - A);
double angle = -angle_between(diff, dir);
double k = round(angle * (double)n / (2.0*M_PI));
- return ref + dir * Rotate(k * 2.0 * M_PI / (double)n) * L2(diff);
+ return A + dir * Rotate(k * 2.0 * M_PI / (double)n) * L2(diff);
}
-} //Namespace Geom
+} //namespace Geom
/*
Local Variables:
diff --git a/src/2geom/point.h b/src/2geom/point.h
index 39538c832686a60b7469fd709ef147d085f558f7..af97cbfa56d04d1864675f738ab2e4112fa4ed58 100644 (file)
--- a/src/2geom/point.h
+++ b/src/2geom/point.h
Point operator/(Point const &p, Matrix const &m);
-/** Constrains the angle between a and b to a multiple of pi/n with respect to dir. */
-Point constrain_angle(Point const &ref, Point const &pt, unsigned int n = 4, Geom::Point const &dir = Geom::Point(1,0));
+/** Constrains the angle (with respect to dir) of the line
+ * joining A and B to a multiple of pi/n.
+ */
+Point constrain_angle(Point const &A, Point const &B, unsigned int n = 4, Geom::Point const &dir = Geom::Point(1,0));
} /* namespace Geom */
diff --git a/src/2geom/poly-dk-solve.cpp b/src/2geom/poly-dk-solve.cpp
+++ /dev/null
@@ -1,67 +0,0 @@
-#include <2geom/poly-dk-solve.h>
-#include <iterator>
-
-/*** implementation of the Durand-Kerner method. seems buggy*/
-
-namespace Geom {
-
-std::complex<double> evalu(Poly const & p, std::complex<double> x) {
- std::complex<double> result = 0;
- std::complex<double> xx = 1;
-
- for(unsigned i = 0; i < p.size(); i++) {
- result += p[i]*xx;
- xx *= x;
- }
- return result;
-}
-
-std::vector<std::complex<double> > DK(Poly const & ply, const double tol) {
- std::vector<std::complex<double> > roots;
- const int N = ply.degree();
-
- std::complex<double> b(0.4, 0.9);
- std::complex<double> p = 1;
- for(int i = 0; i < N; i++) {
- roots.push_back(p);
- p *= b;
- }
- assert(roots.size() == ply.degree());
-
- double error = 0;
- int i;
- for( i = 0; i < 30; i++) {
- error = 0;
- for(int r_i = 0; r_i < N; r_i++) {
- std::complex<double> denom = 1;
- std::complex<double> R = roots[r_i];
- for(int d_i = 0; d_i < N; d_i++) {
- if(r_i != d_i)
- denom *= R-roots[d_i];
- }
- assert(norm(denom) != 0);
- std::complex<double> dr = evalu(ply, R)/denom;
- error += norm(dr);
- roots[r_i] = R - dr;
- }
- /*std::copy(roots.begin(), roots.end(), std::ostream_iterator<std::complex<double> >(std::cout, ",\t"));
- std::cout << std::endl;*/
- if(error < tol)
- break;
- }
- //std::cout << error << ", " << i<< std::endl;
- return roots;
-}
-
-} // 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/poly-dk-solve.h b/src/2geom/poly-dk-solve.h
+++ /dev/null
@@ -1,25 +0,0 @@
-#ifndef LIB2GEOM_SEEN_POLY_DK_SOLVE_H
-#define LIB2GEOM_SEEN_POLY_DK_SOLVE_H
-
-#include <2geom/poly.h>
-#include <complex>
-
-namespace Geom {
-
-std::vector<std::complex<double> >
-DK(Poly const & ply, const double tol=1e-10);
-
-} // namespace Geom
-
-#endif // LIB2GEOM_SEEN_POLY_DK_SOLVE_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:encoding=utf-8:textwidth=99 :
diff --git a/src/2geom/poly-laguerre-solve.cpp b/src/2geom/poly-laguerre-solve.cpp
+++ /dev/null
@@ -1,157 +0,0 @@
-#include <2geom/poly-laguerre-solve.h>
-#include <iterator>
-
-
-namespace Geom {
-
-typedef std::complex<double> cdouble;
-
-cdouble laguerre_internal_complex(Poly const & p,
- double x0,
- double tol,
- bool & quad_root) {
- cdouble a = 2*tol;
- cdouble xk = x0;
- double n = p.degree();
- quad_root = false;
- const unsigned shuffle_rate = 10;
- //static double shuffle[] = {0, 0.5, 0.25, 0.75, 0.125, 0.375, 0.625, 0.875, 1.0};
- unsigned shuffle_counter = 0;
- while(std::norm(a) > (tol*tol)) {
- //std::cout << "xk = " << xk << std::endl;
- cdouble b = p.back();
- cdouble d = 0, f = 0;
- double err = abs(b);
- double abx = abs(xk);
- for(int j = p.size()-2; j >= 0; j--) {
- f = xk*f + d;
- d = xk*d + b;
- b = xk*b + p[j];
- err = abs(b) + abx*err;
- }
-
- err *= 1e-7; // magic epsilon for convergence, should be computed from tol
-
- cdouble px = b;
- if(abs(b) < err)
- return xk;
- //if(std::norm(px) < tol*tol)
- // return xk;
- cdouble G = d / px;
- cdouble H = G*G - f / px;
-
- //std::cout << "G = " << G << "H = " << H;
- cdouble radicand = (n - 1)*(n*H-G*G);
- //assert(radicand.real() > 0);
- if(radicand.real() < 0)
- quad_root = true;
- //std::cout << "radicand = " << radicand << std::endl;
- if(G.real() < 0) // here we try to maximise the denominator avoiding cancellation
- a = - sqrt(radicand);
- else
- a = sqrt(radicand);
- //std::cout << "a = " << a << std::endl;
- a = n / (a + G);
- //std::cout << "a = " << a << std::endl;
- //if(shuffle_counter % shuffle_rate == 0)
- //a *= shuffle[shuffle_counter / shuffle_rate];
- xk -= a;
- shuffle_counter++;
- if(shuffle_counter >= 90)
- break;
- }
- //std::cout << "xk = " << xk << std::endl;
- return xk;
-}
-
-double laguerre_internal(Poly const & p,
- Poly const & pp,
- Poly const & ppp,
- double x0,
- double tol,
- bool & quad_root) {
- double a = 2*tol;
- double xk = x0;
- double n = p.degree();
- quad_root = false;
- while(a*a > (tol*tol)) {
- //std::cout << "xk = " << xk << std::endl;
- double px = p(xk);
- if(px*px < tol*tol)
- return xk;
- double G = pp(xk) / px;
- double H = G*G - ppp(xk) / px;
-
- //std::cout << "G = " << G << "H = " << H;
- double radicand = (n - 1)*(n*H-G*G);
- assert(radicand > 0);
- //std::cout << "radicand = " << radicand << std::endl;
- if(G < 0) // here we try to maximise the denominator avoiding cancellation
- a = - sqrt(radicand);
- else
- a = sqrt(radicand);
- //std::cout << "a = " << a << std::endl;
- a = n / (a + G);
- //std::cout << "a = " << a << std::endl;
- xk -= a;
- }
- //std::cout << "xk = " << xk << std::endl;
- return xk;
-}
-
-
-std::vector<cdouble >
-laguerre(Poly p, const double tol) {
- std::vector<cdouble > solutions;
- //std::cout << "p = " << p << " = ";
- while(p.size() > 1)
- {
- double x0 = 0;
- bool quad_root = false;
- cdouble sol = laguerre_internal_complex(p, x0, tol, quad_root);
- //if(abs(sol) > 1) break;
- Poly dvs;
- if(quad_root) {
- dvs.push_back((sol*conj(sol)).real());
- dvs.push_back(-(sol + conj(sol)).real());
- dvs.push_back(1.0);
- //std::cout << "(" << dvs << ")";
- //solutions.push_back(sol);
- //solutions.push_back(conj(sol));
- } else {
- //std::cout << sol << std::endl;
- dvs.push_back(-sol.real());
- dvs.push_back(1.0);
- solutions.push_back(sol);
- //std::cout << "(" << dvs << ")";
- }
- Poly r;
- p = divide(p, dvs, r);
- //std::cout << r << std::endl;
- }
- return solutions;
-}
-
-std::vector<double>
-laguerre_real_interval(Poly const & /*ply*/,
- const double /*lo*/, const double /*hi*/,
- const double /*tol*/)
-{
- /* not implemented*/
- assert(false);
- std::vector<double> result;
- return result;
-}
-
-} // 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/poly-laguerre-solve.h b/src/2geom/poly-laguerre-solve.h
+++ /dev/null
@@ -1,30 +0,0 @@
-#ifndef LIB2GEOM_SEEN_POLY_LAGUERRE_SOLVE_H
-#define LIB2GEOM_SEEN_POLY_LAGUERRE_SOLVE_H
-
-#include <2geom/poly.h>
-#include <complex>
-
-namespace Geom {
-
-std::vector<std::complex<double> >
-laguerre(Poly ply, const double tol=1e-10);
-
-std::vector<double>
-laguerre_real_interval(Poly ply,
- const double lo, const double hi,
- const double tol=1e-10);
-
-} // namespace Geom
-
-#endif // LIB2GEOM_SEEN_POLY_LAGUERRE_SOLVE_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:encoding=utf-8:textwidth=99 :
diff --git a/src/2geom/poly.cpp b/src/2geom/poly.cpp
index d8b379557f154b812dd151dd62d3a50df3ac73e3..30ec1cdb74686ed1f4197350aba3188363b3804d 100644 (file)
--- a/src/2geom/poly.cpp
+++ b/src/2geom/poly.cpp
#include <2geom/poly.h>
-#define HAVE_GSL
-#ifdef HAVE_GSL
+//#ifdef HAVE_GSL
#include <gsl/gsl_poly.h>
-#endif
+//#endif
namespace Geom {
}
-#ifdef HAVE_GSL
+//#ifdef HAVE_GSL
std::vector<std::complex<double> > solve(Poly const & pp) {
Poly p(pp);
p.normalize();
}
return real_roots;
}
-#endif
+//#endif
double polish_root(Poly const & p, double guess, double tol) {
Poly dp = derivative(p);
diff --git a/src/2geom/recursive-bezier-intersection.cpp b/src/2geom/recursive-bezier-intersection.cpp
--- /dev/null
@@ -0,0 +1,472 @@
+
+
+
+#include <2geom/basic-intersection.h>
+#include <2geom/sbasis-to-bezier.h>
+#include <2geom/exception.h>
+
+
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_multiroots.h>
+
+
+unsigned intersect_steps = 0;
+
+using std::vector;
+namespace Geom {
+
+class OldBezier {
+public:
+ std::vector<Geom::Point> p;
+ OldBezier() {
+ }
+ void split(double t, OldBezier &a, OldBezier &b) const;
+ Point operator()(double t) const;
+
+ ~OldBezier() {}
+
+ void bounds(double &minax, double &maxax,
+ double &minay, double &maxay) {
+ // Compute bounding box for a
+ minax = p[0][X]; // These are the most likely to be extremal
+ maxax = p.back()[X];
+ if( minax > maxax )
+ std::swap(minax, maxax);
+ for(unsigned i = 1; i < p.size()-1; i++) {
+ if( p[i][X] < minax )
+ minax = p[i][X];
+ else if( p[i][X] > maxax )
+ maxax = p[i][X];
+ }
+
+ minay = p[0][Y]; // These are the most likely to be extremal
+ maxay = p.back()[Y];
+ if( minay > maxay )
+ std::swap(minay, maxay);
+ for(unsigned i = 1; i < p.size()-1; i++) {
+ if( p[i][Y] < minay )
+ minay = p[i][Y];
+ else if( p[i][Y] > maxay )
+ maxay = p[i][Y];
+ }
+
+ }
+
+};
+
+static void
+find_intersections_bezier_recursive(std::vector<std::pair<double, double> > & xs,
+ OldBezier a,
+ OldBezier b);
+
+void
+find_intersections_bezier_recursive( std::vector<std::pair<double, double> > &xs,
+ vector<Geom::Point> const & A,
+ vector<Geom::Point> const & B,
+ double precision) {
+ OldBezier a, b;
+ a.p = A;
+ b.p = B;
+ return find_intersections_bezier_recursive(xs, a,b);
+}
+
+
+/* The value of 1.0 / (1L<<14) is enough for most applications */
+const double INV_EPS = (1L<<14);
+
+/*
+ * split the curve at the midpoint, returning an array with the two parts
+ * Temporary storage is minimized by using part of the storage for the result
+ * to hold an intermediate value until it is no longer needed.
+ */
+void OldBezier::split(double t, OldBezier &left, OldBezier &right) const {
+ const unsigned sz = p.size();
+ Geom::Point Vtemp[sz][sz];
+
+ /* Copy control points */
+ std::copy(p.begin(), p.end(), Vtemp[0]);
+
+ /* Triangle computation */
+ for (unsigned i = 1; i < sz; i++) {
+ for (unsigned j = 0; j < sz - i; j++) {
+ Vtemp[i][j] = lerp(t, Vtemp[i-1][j], Vtemp[i-1][j+1]);
+ }
+ }
+
+ left.p.resize(sz);
+ right.p.resize(sz);
+ for (unsigned j = 0; j < sz; j++)
+ left.p[j] = Vtemp[j][0];
+ for (unsigned j = 0; j < sz; j++)
+ right.p[j] = Vtemp[sz-1-j][j];
+}
+
+#if 0
+/*
+ * split the curve at the midpoint, returning an array with the two parts
+ * Temporary storage is minimized by using part of the storage for the result
+ * to hold an intermediate value until it is no longer needed.
+ */
+Point OldBezier::operator()(double t) const {
+ const unsigned sz = p.size();
+ Geom::Point Vtemp[sz][sz];
+
+ /* Copy control points */
+ std::copy(p.begin(), p.end(), Vtemp[0]);
+
+ /* Triangle computation */
+ for (unsigned i = 1; i < sz; i++) {
+ for (unsigned j = 0; j < sz - i; j++) {
+ Vtemp[i][j] = lerp(t, Vtemp[i-1][j], Vtemp[i-1][j+1]);
+ }
+ }
+ return Vtemp[sz-1][0];
+}
+#endif
+
+// suggested by Sederberg.
+Point OldBezier::operator()(double t) const {
+ int n = p.size()-1;
+ double u, bc, tn, tmp;
+ int i;
+ Point r;
+ for(int dim = 0; dim < 2; dim++) {
+ u = 1.0 - t;
+ bc = 1;
+ tn = 1;
+ tmp = p[0][dim]*u;
+ for(i=1; i<n; i++){
+ tn = tn*t;
+ bc = bc*(n-i+1)/i;
+ tmp = (tmp + tn*bc*p[i][dim])*u;
+ }
+ r[dim] = (tmp + tn*t*p[n][dim]);
+ }
+ return r;
+}
+
+
+/*
+ * Test the bounding boxes of two OldBezier curves for interference.
+ * Several observations:
+ * First, it is cheaper to compute the bounding box of the second curve
+ * and test its bounding box for interference than to use a more direct
+ * approach of comparing all control points of the second curve with
+ * the various edges of the bounding box of the first curve to test
+ * for interference.
+ * Second, after a few subdivisions it is highly probable that two corners
+ * of the bounding box of a given Bezier curve are the first and last
+ * control point. Once this happens once, it happens for all subsequent
+ * subcurves. It might be worth putting in a test and then short-circuit
+ * code for further subdivision levels.
+ * Third, in the final comparison (the interference test) the comparisons
+ * should both permit equality. We want to find intersections even if they
+ * occur at the ends of segments.
+ * Finally, there are tighter bounding boxes that can be derived. It isn't
+ * clear whether the higher probability of rejection (and hence fewer
+ * subdivisions and tests) is worth the extra work.
+ */
+
+bool intersect_BB( OldBezier a, OldBezier b ) {
+ double minax, maxax, minay, maxay;
+ a.bounds(minax, maxax, minay, maxay);
+ double minbx, maxbx, minby, maxby;
+ b.bounds(minbx, maxbx, minby, maxby);
+ // Test bounding box of b against bounding box of a
+ // Not >= : need boundary case
+ return not( ( minax > maxbx ) || ( minay > maxby )
+ || ( minbx > maxax ) || ( minby > maxay ) );
+}
+
+/*
+ * Recursively intersect two curves keeping track of their real parameters
+ * and depths of intersection.
+ * The results are returned in a 2-D array of doubles indicating the parameters
+ * for which intersections are found. The parameters are in the order the
+ * intersections were found, which is probably not in sorted order.
+ * When an intersection is found, the parameter value for each of the two
+ * is stored in the index elements array, and the index is incremented.
+ *
+ * If either of the curves has subdivisions left before it is straight
+ * (depth > 0)
+ * that curve (possibly both) is (are) subdivided at its (their) midpoint(s).
+ * the depth(s) is (are) decremented, and the parameter value(s) corresponding
+ * to the midpoints(s) is (are) computed.
+ * Then each of the subcurves of one curve is intersected with each of the
+ * subcurves of the other curve, first by testing the bounding boxes for
+ * interference. If there is any bounding box interference, the corresponding
+ * subcurves are recursively intersected.
+ *
+ * If neither curve has subdivisions left, the line segments from the first
+ * to last control point of each segment are intersected. (Actually the
+ * only the parameter value corresponding to the intersection point is found).
+ *
+ * The apriori flatness test is probably more efficient than testing at each
+ * level of recursion, although a test after three or four levels would
+ * probably be worthwhile, since many curves become flat faster than their
+ * asymptotic rate for the first few levels of recursion.
+ *
+ * The bounding box test fails much more frequently than it succeeds, providing
+ * substantial pruning of the search space.
+ *
+ * Each (sub)curve is subdivided only once, hence it is not possible that for
+ * one final line intersection test the subdivision was at one level, while
+ * for another final line intersection test the subdivision (of the same curve)
+ * was at another. Since the line segments share endpoints, the intersection
+ * is robust: a near-tangential intersection will yield zero or two
+ * intersections.
+ */
+void recursively_intersect( OldBezier a, double t0, double t1, int deptha,
+ OldBezier b, double u0, double u1, int depthb,
+ std::vector<std::pair<double, double> > ¶meters)
+{
+ intersect_steps ++;
+ //std::cout << deptha << std::endl;
+ if( deptha > 0 )
+ {
+ OldBezier A[2];
+ a.split(0.5, A[0], A[1]);
+ double tmid = (t0+t1)*0.5;
+ deptha--;
+ if( depthb > 0 )
+ {
+ OldBezier B[2];
+ b.split(0.5, B[0], B[1]);
+ double umid = (u0+u1)*0.5;
+ depthb--;
+ if( intersect_BB( A[0], B[0] ) )
+ recursively_intersect( A[0], t0, tmid, deptha,
+ B[0], u0, umid, depthb,
+ parameters );
+ if( intersect_BB( A[1], B[0] ) )
+ recursively_intersect( A[1], tmid, t1, deptha,
+ B[0], u0, umid, depthb,
+ parameters );
+ if( intersect_BB( A[0], B[1] ) )
+ recursively_intersect( A[0], t0, tmid, deptha,
+ B[1], umid, u1, depthb,
+ parameters );
+ if( intersect_BB( A[1], B[1] ) )
+ recursively_intersect( A[1], tmid, t1, deptha,
+ B[1], umid, u1, depthb,
+ parameters );
+ }
+ else
+ {
+ if( intersect_BB( A[0], b ) )
+ recursively_intersect( A[0], t0, tmid, deptha,
+ b, u0, u1, depthb,
+ parameters );
+ if( intersect_BB( A[1], b ) )
+ recursively_intersect( A[1], tmid, t1, deptha,
+ b, u0, u1, depthb,
+ parameters );
+ }
+ }
+ else
+ if( depthb > 0 )
+ {
+ OldBezier B[2];
+ b.split(0.5, B[0], B[1]);
+ double umid = (u0 + u1)*0.5;
+ depthb--;
+ if( intersect_BB( a, B[0] ) )
+ recursively_intersect( a, t0, t1, deptha,
+ B[0], u0, umid, depthb,
+ parameters );
+ if( intersect_BB( a, B[1] ) )
+ recursively_intersect( a, t0, t1, deptha,
+ B[0], umid, u1, depthb,
+ parameters );
+ }
+ else // Both segments are fully subdivided; now do line segments
+ {
+ double xlk = a.p.back()[X] - a.p[0][X];
+ double ylk = a.p.back()[Y] - a.p[0][Y];
+ double xnm = b.p.back()[X] - b.p[0][X];
+ double ynm = b.p.back()[Y] - b.p[0][Y];
+ double xmk = b.p[0][X] - a.p[0][X];
+ double ymk = b.p[0][Y] - a.p[0][Y];
+ double det = xnm * ylk - ynm * xlk;
+ if( 1.0 + det == 1.0 )
+ return;
+ else
+ {
+ double detinv = 1.0 / det;
+ double s = ( xnm * ymk - ynm *xmk ) * detinv;
+ double t = ( xlk * ymk - ylk * xmk ) * detinv;
+ if( ( s < 0.0 ) || ( s > 1.0 ) || ( t < 0.0 ) || ( t > 1.0 ) )
+ return;
+ parameters.push_back(std::pair<double, double>(t0 + s * ( t1 - t0 ),
+ u0 + t * ( u1 - u0 )));
+ }
+ }
+}
+
+inline double log4( double x ) { return log(x)/log(4.); }
+
+/*
+ * Wang's theorem is used to estimate the level of subdivision required,
+ * but only if the bounding boxes interfere at the top level.
+ * Assuming there is a possible intersection, recursively_intersect is
+ * used to find all the parameters corresponding to intersection points.
+ * these are then sorted and returned in an array.
+ */
+
+double Lmax(Point p) {
+ return std::max(fabs(p[X]), fabs(p[Y]));
+}
+
+unsigned wangs_theorem(OldBezier a) {
+ return 6; // seems a good approximation!
+ double la1 = Lmax( ( a.p[2] - a.p[1] ) - (a.p[1] - a.p[0]) );
+ double la2 = Lmax( ( a.p[3] - a.p[2] ) - (a.p[2] - a.p[1]) );
+ double l0 = std::max(la1, la2);
+ unsigned ra;
+ if( l0 * 0.75 * M_SQRT2 + 1.0 == 1.0 )
+ ra = 0;
+ else
+ ra = (unsigned)ceil( log4( M_SQRT2 * 6.0 / 8.0 * INV_EPS * l0 ) );
+ //std::cout << ra << std::endl;
+ return ra;
+}
+
+struct rparams
+{
+ OldBezier &A;
+ OldBezier &B;
+};
+
+static int
+intersect_polish_f (const gsl_vector * x, void *params,
+ gsl_vector * f)
+{
+ const double x0 = gsl_vector_get (x, 0);
+ const double x1 = gsl_vector_get (x, 1);
+
+ Geom::Point dx = ((struct rparams *) params)->A(x0) -
+ ((struct rparams *) params)->B(x1);
+
+ gsl_vector_set (f, 0, dx[0]);
+ gsl_vector_set (f, 1, dx[1]);
+
+ return GSL_SUCCESS;
+}
+
+union dbl_64{
+ long long i64;
+ double d64;
+};
+
+static double EpsilonBy(double value, int eps)
+{
+ dbl_64 s;
+ s.d64 = value;
+ s.i64 += eps;
+ return s.d64;
+}
+
+
+static void intersect_polish_root (OldBezier &A, double &s,
+ OldBezier &B, double &t) {
+ const gsl_multiroot_fsolver_type *T;
+ gsl_multiroot_fsolver *sol;
+
+ int status;
+ size_t iter = 0;
+
+ const size_t n = 2;
+ struct rparams p = {A, B};
+ gsl_multiroot_function f = {&intersect_polish_f, n, &p};
+
+ double x_init[2] = {s, t};
+ gsl_vector *x = gsl_vector_alloc (n);
+
+ gsl_vector_set (x, 0, x_init[0]);
+ gsl_vector_set (x, 1, x_init[1]);
+
+ T = gsl_multiroot_fsolver_hybrids;
+ sol = gsl_multiroot_fsolver_alloc (T, 2);
+ gsl_multiroot_fsolver_set (sol, &f, x);
+
+ do
+ {
+ iter++;
+ status = gsl_multiroot_fsolver_iterate (sol);
+
+ if (status) /* check if solver is stuck */
+ break;
+
+ status =
+ gsl_multiroot_test_residual (sol->f, 1e-12);
+ }
+ while (status == GSL_CONTINUE && iter < 1000);
+
+ s = gsl_vector_get (sol->x, 0);
+ t = gsl_vector_get (sol->x, 1);
+
+ gsl_multiroot_fsolver_free (sol);
+ gsl_vector_free (x);
+
+ // This code does a neighbourhood search for minor improvements.
+ double best_v = L1(A(s) - B(t));
+ //std::cout << "------\n" << best_v << std::endl;
+ Point best(s,t);
+ while (true) {
+ Point trial = best;
+ double trial_v = best_v;
+ for(int nsi = -1; nsi < 2; nsi++) {
+ for(int nti = -1; nti < 2; nti++) {
+ Point n(EpsilonBy(best[0], nsi),
+ EpsilonBy(best[1], nti));
+ double c = L1(A(n[0]) - B(n[1]));
+ //std::cout << c << "; ";
+ if (c < trial_v) {
+ trial = n;
+ trial_v = c;
+ }
+ }
+ }
+ if(trial == best) {
+ //std::cout << "\n" << s << " -> " << s - best[0] << std::endl;
+ //std::cout << t << " -> " << t - best[1] << std::endl;
+ //std::cout << best_v << std::endl;
+ s = best[0];
+ t = best[1];
+ return;
+ } else {
+ best = trial;
+ best_v = trial_v;
+ }
+ }
+}
+
+
+void find_intersections_bezier_recursive( std::vector<std::pair<double, double> > &xs,
+ OldBezier a, OldBezier b)
+{
+ if( intersect_BB( a, b ) )
+ {
+ recursively_intersect( a, 0., 1., wangs_theorem(a),
+ b, 0., 1., wangs_theorem(b),
+ xs);
+ }
+ /*for(unsigned i = 0; i < xs.size(); i++)
+ intersect_polish_root(a, xs[i].first,
+ b, xs[i].second);*/
+ std::sort(xs.begin(), xs.end());
+}
+
+
+};
+
+/*
+ Local Variables:
+ mode:c++
+ c-file-style:"stroustrup"
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+ indent-tabs-mode:nil
+ fill-column:99
+ End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
index 28e9fc964006df620bfe2c2b9353787d05f200c1..c3711840270f31d0d0c2c40a9d0122ca496a6bfc 100644 (file)
return cubics_fitting_curvature(M0,M1,dM0,dM1,d2M0xdM0,d2M1xdM1,insist_on_speed_signs,epsilon);
}
+
+/**
+* \brief returns all the parameter values of A whose tangent passes through P.
+*/
+std::vector<double> find_tangents(Point P, D2<SBasis> const &A) {
+ SBasis crs (cross(A - P, derivative(A)));
+ crs = shift(crs*Linear(-1, 0)*Linear(-1, 0), -2); // We know that there is a double root at t=0 so we divide out t^2
+// JFB points out that this is equivalent to (t-1)^2 followed by a divide by s^2 (shift)
+ return roots(crs);
+}
+
+
//}; // namespace
index 18c666b11141b4cd2513fdc12ff8d11ce4ce9955..4f249a7b1888060e43f70072ac58d7e2817e0443 100644 (file)
int insist_on_speed_signs = 1,
double error = 1e-5);
+
+std::vector<double> find_tangents(Point P, D2<SBasis> const &A);
+
};
#endif
index d7045a3a9b879adfa20f54ca259edece3a3640fd..7053388dce013ce1635bbc1571e3e9e266247948 100644 (file)
//TODO: in all these functions, compute 'order' according to 'tol'.
#include <2geom/sbasis-math.h>
-//#define ZERO 1e-3
+#include <2geom/d2-sbasis.h>
#include <stdio.h>
#include <math.h>
+//#define ZERO 1e-3
+
namespace Geom {
-#include <2geom/d2-sbasis.h>
//-|x|-----------------------------------------------------------------------
/** Return the absolute value of a function pointwise.
index dfb24f9d90275a423d0fd773b2e039bb3540e62a..ce5bf89bce3999a9e7a2e4e21e0ad8a7bc28bdb3 100644 (file)
if the degree is even q is the order in the symmetrical power basis,
if the degree is odd q is the order + 1
n is always the polynomial degree, i. e. the Bezier order
+ sz is the number of bezier handles.
*/
void sbasis_to_bezier (Bezier & bz, SBasis const& sb, size_t sz)
{
}
else
{
- q = (sz > sb.size()) ? sb.size() : sz;
- n = 2*sz-1;
+ q = (sz > 2*sb.size()-1) ? sb.size() : (sz+1)/2;
+ n = sz-1;
even = false;
}
bz.clear();
\param p the D2 Symmetric basis polynomial
\returns the D2 Bernstein basis polynomial
- if the degree is even q is the order in the symmetrical power basis,
- if the degree is odd q is the order + 1
- n is always the polynomial degree, i. e. the Bezier order
+ sz is always the polynomial degree, i. e. the Bezier order
*/
void sbasis_to_bezier (std::vector<Point> & bz, D2<SBasis> const& sb, size_t sz)
{
Bezier bzx, bzy;
+ if(sz == 0) {
+ sz = std::max(sb[X].size(), sb[Y].size())*2;
+ }
sbasis_to_bezier(bzx, sb[X], sz);
sbasis_to_bezier(bzy, sb[Y], sz);
+ assert(bzx.size() == bzy.size());
size_t n = (bzx.size() >= bzy.size()) ? bzx.size() : bzy.size();
bz.resize(n, Point(0,0));
@@ -345,7 +348,7 @@ void build_from_sbasis(Geom::PathBuilder &pb, D2<SBasis> const &B, double tol, b
pb.lineTo(B.at1());
} else {
std::vector<Geom::Point> bez;
- sbasis_to_bezier(bez, B, 2);
+ sbasis_to_bezier(bez, B, 4);
pb.curveTo(bez[1], bez[2], bez[3]);
}
} else {
diff --git a/src/2geom/sbasis.cpp b/src/2geom/sbasis.cpp
index b5c1a05a7bec28c6c6dfeb8f052466d154a65745..49678b19ab56e0a21878bdb3505deae115e030d1 100644 (file)
--- a/src/2geom/sbasis.cpp
+++ b/src/2geom/sbasis.cpp
for(int i = 0; i < sh; i++)
c[i] = Linear(0,0);
- for(size_t i = m, j = 0; i < n; i++, j++)
+ for(size_t i = m, j = std::max(0,-sh); i < n; i++, j++)
c[i] = a[j];
return c;
}