From 07717cb695cd79993e74908804d04e658a95482f Mon Sep 17 00:00:00 2001 From: jfbarraud Date: Mon, 9 Mar 2009 01:47:39 +0000 Subject: [PATCH] 2geom update --- src/2geom/CMakeLists.txt | 1 + src/2geom/Makefile_insert | 8 +- src/2geom/basic-intersection.cpp | 38 +- src/2geom/basic-intersection.h | 77 +- src/2geom/bezier-clipping.cpp | 1292 +++++++++++++++++++ src/2geom/point.cpp | 12 +- src/2geom/point.h | 6 +- src/2geom/poly-dk-solve.cpp | 67 - src/2geom/poly-dk-solve.h | 25 - src/2geom/poly-laguerre-solve.cpp | 157 --- src/2geom/poly-laguerre-solve.h | 30 - src/2geom/poly.cpp | 9 +- src/2geom/recursive-bezier-intersection.cpp | 472 +++++++ src/2geom/sbasis-geometric.cpp | 12 + src/2geom/sbasis-geometric.h | 3 + src/2geom/sbasis-math.cpp | 5 +- src/2geom/sbasis-to-bezier.cpp | 15 +- src/2geom/sbasis.cpp | 2 +- 18 files changed, 1883 insertions(+), 348 deletions(-) create mode 100644 src/2geom/bezier-clipping.cpp delete mode 100644 src/2geom/poly-dk-solve.cpp delete mode 100644 src/2geom/poly-dk-solve.h delete mode 100644 src/2geom/poly-laguerre-solve.cpp delete mode 100644 src/2geom/poly-laguerre-solve.h create mode 100644 src/2geom/recursive-bezier-intersection.cpp diff --git a/src/2geom/CMakeLists.txt b/src/2geom/CMakeLists.txt index b5461487d..3c8669a9d 100644 --- a/src/2geom/CMakeLists.txt +++ b/src/2geom/CMakeLists.txt @@ -90,6 +90,7 @@ bezier-to-sbasis.h basic-intersection.h basic-intersection.cpp +recursive-bezier-intersection.cpp geom.cpp geom.h diff --git a/src/2geom/Makefile_insert b/src/2geom/Makefile_insert index d03799dd7..901e22f40 100644 --- a/src/2geom/Makefile_insert +++ b/src/2geom/Makefile_insert @@ -7,6 +7,10 @@ 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 \ @@ -26,8 +30,6 @@ 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 \ @@ -82,8 +84,6 @@ 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 \ diff --git a/src/2geom/basic-intersection.cpp b/src/2geom/basic-intersection.cpp index 16b5c0240..c03023e6f 100644 --- a/src/2geom/basic-intersection.cpp +++ b/src/2geom/basic-intersection.cpp @@ -10,6 +10,28 @@ using std::vector; namespace Geom { +//#ifdef USE_RECURSIVE_INTERSECTOR + +// void find_intersections(std::vector > &xs, +// D2 const & A, +// D2 const & B) { +// vector 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 > & xs, +// std::vector const& A, +// std::vector const& B, +// double precision){ +// find_intersections_bezier_recursive(xs, A, B, precision); +// } + +//#else + namespace detail{ namespace bezier_clipping { void portion (std::vector & B, Interval const& I); }; }; @@ -20,12 +42,20 @@ void find_intersections(std::vector > &xs, vector 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 > & xs, + std::vector const& A, + std::vector 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 @@ -79,6 +109,7 @@ find_self_intersections(std::vector > &xs, in = r; } } + for(unsigned i = 0; i < dr.size()-1; i++) { for(unsigned j = i+1; j < dr.size()-1; j++) { std::vector > section; @@ -89,7 +120,8 @@ find_self_intersections(std::vector > &xs, 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])); diff --git a/src/2geom/basic-intersection.h b/src/2geom/basic-intersection.h index 783df370e..a19a10c8c 100644 --- a/src/2geom/basic-intersection.h +++ b/src/2geom/basic-intersection.h @@ -44,38 +44,44 @@ #define USE_RECURSIVE_INTERSECTOR 0 -namespace Geom { - -#ifdef USE_RECURSIVE_INTERSECTOR -void -find_intersections( std::vector > &xs, - D2 const & A, - D2 const & B); -void -find_self_intersections(std::vector > &xs, - D2 const & A); +namespace Geom { -// Bezier form -void -find_intersections( std::vector > &xs, - std::vector const & A, - std::vector const & B, - double precision = 1e-5); +//why not allowing precision to be set here? +void find_intersections(std::vector >& xs, + D2 const & A, + D2 const & B); -// Bezier form -void -find_intersections_bezier_clipping( std::vector > &xs, - std::vector const & A, - std::vector const & B, - double precision = 1e-5); +void find_intersections(std::vector< std::pair > & xs, + std::vector const& A, + std::vector const& B, + double precision = 1e-5); -void -find_self_intersections(std::vector > &xs, - std::vector const & A); +//why not allowing precision to be set here? +void find_self_intersections(std::vector >& xs, + D2 const & A); -#else +//--not implemented +//void find_self_intersections(std::vector >& xs, +// std::vector 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 > & xs, +// std::vector const& A, +// std::vector const& B, +// double precision = 1e-5); +// #else /* * find_intersection * @@ -87,11 +93,14 @@ find_self_intersections(std::vector > &xs, * 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 > & xs, +void find_intersections_bezier_clipping (std::vector< std::pair > & xs, std::vector const& A, std::vector const& B, double precision = 1e-5); -#endif +//#endif + + + /* * find_collinear_normal * @@ -112,18 +121,6 @@ void polish_intersections(std::vector > &xs, D2 const &A, D2 const &B); -void find_intersections(std::vector >& xs, - D2 const & A, - D2 const & B); - -void find_self_intersections(std::vector >& xs, - D2 const & A); - - -void find_self_intersections(std::vector >& xs, - std::vector 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 new file mode 100644 index 000000000..96a06376c --- /dev/null +++ b/src/2geom/bezier-clipping.cpp @@ -0,0 +1,1292 @@ +/* + * Implement the Bezier clipping algorithm for finding + * Bezier curve intersection points and collinear normals + * + * Authors: + * 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/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 +#include +#include +#include +//#include + + + + +#define VERBOSE 0 +#define CHECK 0 + + +namespace Geom { + +namespace detail { namespace bezier_clipping { + +//////////////////////////////////////////////////////////////////////////////// +// for debugging +// + +inline +void print(std::vector 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 & +operator<< (std::basic_ostream & os, const Interval & I) +{ + os << "[" << I.min() << ", " << I.max() << "]"; + return os; +} + +inline +double angle (std::vector 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 & 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(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 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 & D, std::vector 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 & N, std::vector 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 & 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 & 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 & 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 +void clip(Interval & dom, + std::vector const& A, + std::vector const& B); +template +void iterate(std::vector& domsA, + std::vector& domsB, + std::vector const& A, + std::vector 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 & l, + std::vector 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 & l, + std::vector 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 & l, + std::vector const& c, + Point const& p) +{ + if (is_constant(c)) + { + // this should never happen + assert(!is_constant(c)); + } + std::vector 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 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 const& c, + std::vector 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 const& B, + std::vector const& l, + Interval const& bound) +{ + double n = B.size() - 1; // number of sub-intervals + std::vector 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 & 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 (Interval & dom, + std::vector const& A, + std::vector const& B) +{ + std::vector 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 & F, std::vector 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 & D, + std::vector const& B, + std::vector 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 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 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 const& B, + std::vector const& F) +{ + std::vector D; // distance curve control points + distance_control_points(D, B, F); + //print(D, "D"); +// ConvexHull chD(D); +// std::vector& p = chD.boundary; // convex hull vertices + + convex_hull(D); + std::vector & 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 (Interval & dom, + std::vector const& A, + std::vector const& B) +{ + std::vector 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 (std::vector& domsA, + std::vector& domsB, + std::vector const& A, + std::vector 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 pA = A; + std::vector pB = B; + std::vector* C1 = &pA; + std::vector* 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(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 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(domsA, domsB, pC1, pB, + dompC1, dompB, precision); + iterate(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(domsB, domsA, pC1, pA, + dompC1, dompA, precision); + iterate(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 (std::vector& domsA, + std::vector& domsB, + std::vector const& A, + std::vector 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 pA = A; + std::vector pB = B; + std::vector* C1 = &pA; + std::vector* 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(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 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(domsA, domsB, pC1, pB, + dompC1, dompB, precision); + iterate(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(domsB, domsA, pC1, pA, + dompC1, dompA, precision); + iterate(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 +void get_solutions (std::vector< std::pair >& xs, + std::vector const& A, + std::vector const& B, + double precision) +{ + std::pair ci; + std::vector domsA, domsB; + iterate (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 >& xs, + std::vector const& A, + std::vector const& B, + double precision) +{ + using detail::bezier_clipping::get_solutions; + using detail::bezier_clipping::collinear_normal_tag; + get_solutions(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 >& xs, + std::vector const& A, + std::vector const& B, + double precision) +{ + using detail::bezier_clipping::get_solutions; + using detail::bezier_clipping::intersection_point_tag; + get_solutions(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 b8a858f8d..4a0625713 100644 --- a/src/2geom/point.cpp +++ b/src/2geom/point.cpp @@ -150,19 +150,19 @@ Point &Point::operator*=(Matrix const &m) 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 39538c832..af97cbfa5 100644 --- a/src/2geom/point.h +++ b/src/2geom/point.h @@ -229,8 +229,10 @@ Point operator*(Point const &v, Matrix const &m); 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 deleted file mode 100644 index a328b316c..000000000 --- a/src/2geom/poly-dk-solve.cpp +++ /dev/null @@ -1,67 +0,0 @@ -#include <2geom/poly-dk-solve.h> -#include - -/*** implementation of the Durand-Kerner method. seems buggy*/ - -namespace Geom { - -std::complex evalu(Poly const & p, std::complex x) { - std::complex result = 0; - std::complex xx = 1; - - for(unsigned i = 0; i < p.size(); i++) { - result += p[i]*xx; - xx *= x; - } - return result; -} - -std::vector > DK(Poly const & ply, const double tol) { - std::vector > roots; - const int N = ply.degree(); - - std::complex b(0.4, 0.9); - std::complex 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 denom = 1; - std::complex 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 dr = evalu(ply, R)/denom; - error += norm(dr); - roots[r_i] = R - dr; - } - /*std::copy(roots.begin(), roots.end(), std::ostream_iterator >(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 deleted file mode 100644 index 463f854c1..000000000 --- a/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 - -namespace Geom { - -std::vector > -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 deleted file mode 100644 index 5b59690b7..000000000 --- a/src/2geom/poly-laguerre-solve.cpp +++ /dev/null @@ -1,157 +0,0 @@ -#include <2geom/poly-laguerre-solve.h> -#include - - -namespace Geom { - -typedef std::complex 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 -laguerre(Poly p, const double tol) { - std::vector 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 -laguerre_real_interval(Poly const & /*ply*/, - const double /*lo*/, const double /*hi*/, - const double /*tol*/) -{ - /* not implemented*/ - assert(false); - std::vector 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 deleted file mode 100644 index 6836aa49d..000000000 --- a/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 - -namespace Geom { - -std::vector > -laguerre(Poly ply, const double tol=1e-10); - -std::vector -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 d8b379557..30ec1cdb7 100644 --- a/src/2geom/poly.cpp +++ b/src/2geom/poly.cpp @@ -1,9 +1,8 @@ #include <2geom/poly.h> -#define HAVE_GSL -#ifdef HAVE_GSL +//#ifdef HAVE_GSL #include -#endif +//#endif namespace Geom { @@ -39,7 +38,7 @@ void Poly::monicify() { } -#ifdef HAVE_GSL +//#ifdef HAVE_GSL std::vector > solve(Poly const & pp) { Poly p(pp); p.normalize(); @@ -76,7 +75,7 @@ std::vector solve_reals(Poly const & p) { } 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 new file mode 100644 index 000000000..d59e7d9c9 --- /dev/null +++ b/src/2geom/recursive-bezier-intersection.cpp @@ -0,0 +1,472 @@ + + + +#include <2geom/basic-intersection.h> +#include <2geom/sbasis-to-bezier.h> +#include <2geom/exception.h> + + +#include +#include + + +unsigned intersect_steps = 0; + +using std::vector; +namespace Geom { + +class OldBezier { +public: + std::vector 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 > & xs, + OldBezier a, + OldBezier b); + +void +find_intersections_bezier_recursive( std::vector > &xs, + vector const & A, + vector 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= : 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 > ¶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(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 > &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 : diff --git a/src/2geom/sbasis-geometric.cpp b/src/2geom/sbasis-geometric.cpp index 28e9fc964..c37118402 100644 --- a/src/2geom/sbasis-geometric.cpp +++ b/src/2geom/sbasis-geometric.cpp @@ -740,6 +740,18 @@ Geom::cubics_with_prescribed_curvature(Point const &M0, Point const &M1, 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 find_tangents(Point P, D2 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 diff --git a/src/2geom/sbasis-geometric.h b/src/2geom/sbasis-geometric.h index 18c666b11..4f249a7b1 100644 --- a/src/2geom/sbasis-geometric.h +++ b/src/2geom/sbasis-geometric.h @@ -99,6 +99,9 @@ cubics_with_prescribed_curvature(Point const &M0, Point const &M1, int insist_on_speed_signs = 1, double error = 1e-5); + +std::vector find_tangents(Point P, D2 const &A); + }; #endif diff --git a/src/2geom/sbasis-math.cpp b/src/2geom/sbasis-math.cpp index d7045a3a9..7053388dc 100644 --- a/src/2geom/sbasis-math.cpp +++ b/src/2geom/sbasis-math.cpp @@ -35,14 +35,15 @@ //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 #include +//#define ZERO 1e-3 + namespace Geom { -#include <2geom/d2-sbasis.h> //-|x|----------------------------------------------------------------------- /** Return the absolute value of a function pointwise. diff --git a/src/2geom/sbasis-to-bezier.cpp b/src/2geom/sbasis-to-bezier.cpp index dfb24f9d9..ce5bf89bc 100644 --- a/src/2geom/sbasis-to-bezier.cpp +++ b/src/2geom/sbasis-to-bezier.cpp @@ -95,6 +95,7 @@ int sgn(unsigned int j, unsigned int k) 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) { @@ -117,8 +118,8 @@ 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(); @@ -151,15 +152,17 @@ void sbasis_to_bezier (Bezier & bz, SBasis const& sb, size_t sz) \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 & bz, D2 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 const &B, double tol, b pb.lineTo(B.at1()); } else { std::vector 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 b5c1a05a7..49678b19a 100644 --- a/src/2geom/sbasis.cpp +++ b/src/2geom/sbasis.cpp @@ -198,7 +198,7 @@ SBasis shift(SBasis const &a, int sh) { 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; } -- 2.30.2