From: johanengelen Date: Mon, 6 Apr 2009 22:29:34 +0000 (+0000) Subject: update 2geom. big commit as it has been a while. (2geom svn rev. 1870, i think) X-Git-Url: https://git.tokkee.org/?a=commitdiff_plain;h=98a704c566ce5a750d76d5fc9675ccc804ac65f5;p=inkscape.git update 2geom. big commit as it has been a while. (2geom svn rev. 1870, i think) i turned some optional compilation stuff *on* per default, to help building inkscape. --- diff --git a/src/2geom/basic-intersection.cpp b/src/2geom/basic-intersection.cpp index 0c8b96ddf..e159839d2 100644 --- a/src/2geom/basic-intersection.cpp +++ b/src/2geom/basic-intersection.cpp @@ -2,10 +2,10 @@ #include <2geom/sbasis-to-bezier.h> #include <2geom/exception.h> - +#ifdef HAVE_GSL #include #include - +#endif using std::vector; namespace Geom { @@ -135,10 +135,8 @@ find_self_intersections(std::vector > &xs, //unique(xs.begin(), xs.end()); } - +#ifdef HAVE_GSL #include - - struct rparams { @@ -161,6 +159,7 @@ intersect_polish_f (const gsl_vector * x, void *params, return GSL_SUCCESS; } +#endif union dbl_64{ long long i64; @@ -178,12 +177,52 @@ static double EpsilonBy(double value, int eps) static void intersect_polish_root (D2 const &A, double &s, D2 const &B, double &t) { +#ifdef HAVE_GSL const gsl_multiroot_fsolver_type *T; gsl_multiroot_fsolver *sol; int status; size_t iter = 0; - +#endif + std::vector as, bs; + as = A.valueAndDerivatives(s, 2); + bs = B.valueAndDerivatives(t, 2); + Point F = as[0] - bs[0]; + double best = dot(F, F); + + for(int i = 0; i < 4; i++) { + + /** + we want to solve + J*(x1 - x0) = f(x0) + + |dA(s)[0] -dB(t)[0]| (X1 - X0) = A(s) - B(t) + |dA(s)[1] -dB(t)[1]| + **/ + + // We're using the standard transformation matricies, which is numerically rather poor. Much better to solve the equation using elimination. + + Matrix jack(as[1][0], as[1][1], + -bs[1][0], -bs[1][1], + 0, 0); + Point soln = (F)*jack.inverse(); + double ns = s - soln[0]; + double nt = t - soln[1]; + + as = A.valueAndDerivatives(ns, 2); + bs = B.valueAndDerivatives(nt, 2); + F = as[0] - bs[0]; + double trial = dot(F, F); + if (trial > best*0.1) {// we have standards, you know + // At this point we could do a line search + break; + } + best = trial; + s = ns; + t = nt; + } + +#ifdef HAVE_GSL const size_t n = 2; struct rparams p = {A, B}; gsl_multiroot_function f = {&intersect_polish_f, n, &p}; @@ -216,7 +255,9 @@ static void intersect_polish_root (D2 const &A, double &s, gsl_multiroot_fsolver_free (sol); gsl_vector_free (x); +#endif + { // This code does a neighbourhood search for minor improvements. double best_v = L1(A(s) - B(t)); //std::cout << "------\n" << best_v << std::endl; @@ -248,6 +289,7 @@ static void intersect_polish_root (D2 const &A, double &s, best_v = trial_v; } } + } } diff --git a/src/2geom/bezier-clipping.cpp b/src/2geom/bezier-clipping.cpp index 96a06376c..c8d99c430 100644 --- a/src/2geom/bezier-clipping.cpp +++ b/src/2geom/bezier-clipping.cpp @@ -1,1292 +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 : +/* + * 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/bezier-curve.h b/src/2geom/bezier-curve.h index acf878ee8..d5259c71f 100644 --- a/src/2geom/bezier-curve.h +++ b/src/2geom/bezier-curve.h @@ -166,15 +166,8 @@ public: return ret; } - Curve *derivative() const { - if(order > 1) - return new BezierCurve(Geom::derivative(inner[X]), Geom::derivative(inner[Y])); - else if (order == 1) { - double dx = inner[X][1] - inner[X][0], dy = inner[Y][1] - inner[Y][0]; - return new BezierCurve<1>(Point(dx,dy),Point(dx,dy)); - } - } - + Curve *derivative() const; + Point pointAt(double t) const { return inner.valueAt(t); } std::vector pointAndDerivatives(Coord t, unsigned n) const { return inner.valueAndDerivatives(t, n); } @@ -193,7 +186,7 @@ protected: }; // BezierCurve<0> is meaningless; specialize it out -template<> class BezierCurve<0> : public BezierCurve<1> { public: BezierCurve(); BezierCurve(Bezier x, Bezier y); }; +template<> class BezierCurve<0> : public BezierCurve<1> { public: BezierCurve();}; typedef BezierCurve<1> LineSegment; typedef BezierCurve<2> QuadraticBezier; @@ -228,6 +221,20 @@ double length(LineSegment const& _segment) return distance(_segment.initialPoint(), _segment.finalPoint()); } +template +inline +Curve *BezierCurve::derivative() const { + return new BezierCurve(Geom::derivative(inner[X]), Geom::derivative(inner[Y])); +} + +template <> +inline +Curve *BezierCurve<1>::derivative() const { + double dx = inner[X][1] - inner[X][0], dy = inner[Y][1] - inner[Y][0]; + return new BezierCurve<1>(Point(dx,dy),Point(dx,dy)); +} + + } // end namespace Geom diff --git a/src/2geom/bezier.h b/src/2geom/bezier.h index 889dde9ed..4ab965f42 100644 --- a/src/2geom/bezier.h +++ b/src/2geom/bezier.h @@ -44,12 +44,44 @@ namespace Geom { inline Coord subdivideArr(Coord t, Coord const *v, Coord *left, Coord *right, unsigned order) { +/* + * Bernstein : + * Evaluate a Bernstein function at a particular parameter value + * Fill in control points for resulting sub-curves. + * + */ + + unsigned N = order+1; + std::valarray vtemp(2*N); + for (unsigned i = 0; i < N; i++) + vtemp[i] = v[i]; + + // Triangle computation + const double omt = (1-t); + if(left) + left[0] = vtemp[0]; + if(right) + right[order] = vtemp[order]; + double *prev_row = &vtemp[0]; + double *row = &vtemp[N]; + for (unsigned i = 1; i < N; i++) { + for (unsigned j = 0; j < N - i; j++) { + row[j] = omt*prev_row[j] + t*prev_row[j+1]; + } + if(left) + left[i] = row[0]; + if(right) + right[order-i] = row[order-i]; + std::swap(prev_row, row); + } + return (prev_row[0]); +/* Coord vtemp[order+1][order+1]; - /* Copy control points */ + // Copy control points std::copy(v, v+order+1, vtemp[0]); - /* Triangle computation */ + // Triangle computation for (unsigned i = 1; i <= order; i++) { for (unsigned j = 0; j <= order - i; j++) { vtemp[i][j] = lerp(t, vtemp[i-1][j], vtemp[i-1][j+1]); @@ -62,7 +94,7 @@ inline Coord subdivideArr(Coord t, Coord const *v, Coord *left, Coord *right, un for (unsigned j = 0; j <= order; j++) right[j] = vtemp[order-j][j]; - return (vtemp[order][0]); + return (vtemp[order][0]);*/ } @@ -162,7 +194,20 @@ public: inline Coord at1() const { return c_[order()]; } inline Coord valueAt(double t) const { - return subdivideArr(t, &c_[0], NULL, NULL, order()); + int n = order(); + double u, bc, tn, tmp; + int i; + u = 1.0 - t; + bc = 1; + tn = 1; + tmp = c_[0]*u; + for(i=1; i&>(c_)[ix]; } + //inline Coord const &operator[](unsigned ix) const { return c_[ix]; } inline void setPoint(unsigned ix, double val) { c_[ix] = val; } /* This is inelegant, as it uses several extra stores. I think there might be a way to @@ -184,14 +230,14 @@ public: std::vector valueAndDerivatives(Coord t, unsigned n_derivs) const { std::vector val_n_der; - Coord d_[order()+1]; + std::valarray d_(order()+1); unsigned nn = n_derivs + 1; // the size of the result vector equals n_derivs+1 ... if(nn > order()) - nn = order()+1; // .. but with a maximum of order() + 1! + nn = order()+1; // .. but with a maximum of order() + 1! for(unsigned i = 0; i < size(); i++) d_[i] = c_[i]; for(unsigned di = 0; di < nn; di++) { - val_n_der.push_back(subdivideArr(t, d_, NULL, NULL, order() - di)); + val_n_der.push_back(subdivideArr(t, &d_[0], NULL, NULL, order() - di)); for(unsigned i = 0; i < order() - di; i++) { d_[i] = (order()-di)*(d_[i+1] - d_[i]); } @@ -202,13 +248,13 @@ public: std::pair subdivide(Coord t) const { Bezier a(Bezier::Order(*this)), b(Bezier::Order(*this)); - subdivideArr(t, &c_[0], &a.c_[0], &b.c_[0], order()); + subdivideArr(t, &const_cast&>(c_)[0], &a.c_[0], &b.c_[0], order()); return std::pair(a, b); } std::vector roots() const { std::vector solutions; - find_bernstein_roots(&c_[0], order(), solutions, 0, 0.0, 1.0); + find_bernstein_roots(&const_cast&>(c_)[0], order(), solutions, 0, 0.0, 1.0); return solutions; } }; @@ -263,17 +309,17 @@ inline Bezier reverse(const Bezier & a) { inline Bezier portion(const Bezier & a, double from, double to) { //TODO: implement better? - Coord res[a.order()+1]; + std::valarray res(a.order() + 1); if(from == 0) { if(to == 1) { return Bezier(a); } - subdivideArr(to, &a.c_[0], res, NULL, a.order()); - return Bezier(res, a.order()); + subdivideArr(to, &const_cast(a).c_[0], &res[0], NULL, a.order()); + return Bezier(&res[0], a.order()); } - subdivideArr(from, &a.c_[0], NULL, res, a.order()); - if(to == 1) return Bezier(res, a.order()); - Coord res2[a.order()+1]; - subdivideArr((to - from)/(1 - from), res, res2, NULL, a.order()); - return Bezier(res2, a.order()); + subdivideArr(from, &const_cast(a).c_[0], NULL, &res[0], a.order()); + if(to == 1) return Bezier(&res[0], a.order()); + std::valarray res2(a.order()+1); + subdivideArr((to - from)/(1 - from), &res[0], &res2[0], NULL, a.order()); + return Bezier(&res2[0], a.order()); } // XXX Todo: how to handle differing orders @@ -309,7 +355,7 @@ inline Bezier integral(const Bezier & a) { } inline OptInterval bounds_fast(Bezier const & b) { - return Interval::fromArray(&b.c_[0], b.size()); + return Interval::fromArray(&const_cast(b).c_[0], b.size()); } //TODO: better bounds exact diff --git a/src/2geom/chebyshev.cpp b/src/2geom/chebyshev.cpp index 447c5183f..4ba3e2325 100644 --- a/src/2geom/chebyshev.cpp +++ b/src/2geom/chebyshev.cpp @@ -1,126 +1,126 @@ -#include <2geom/chebyshev.h> - -#include <2geom/sbasis.h> -#include <2geom/sbasis-poly.h> - -#include -using std::vector; - -#include -#include - -namespace Geom{ - -SBasis cheb(unsigned n) { - static std::vector basis; - if(basis.empty()) { - basis.push_back(Linear(1,1)); - basis.push_back(Linear(0,1)); - } - for(unsigned i = basis.size(); i <= n; i++) { - basis.push_back(Linear(0,2)*basis[i-1] - basis[i-2]); - } - - return basis[n]; -} - -SBasis cheb_series(unsigned n, double* cheb_coeff) { - SBasis r; - for(unsigned i = 0; i < n; i++) { - double cof = cheb_coeff[i]; - //if(i == 0) - //cof /= 2; - r += cheb(i)*cof; - } - - return r; -} - -SBasis clenshaw_series(unsigned m, double* cheb_coeff) { - /** b_n = a_n - b_n-1 = 2*x*b_n + a_n-1 - b_n-k = 2*x*b_{n-k+1} + a_{n-k} - b_{n - k + 2} - b_0 = x*b_1 + a_0 - b_2 - */ - - double a = -1, b = 1; - SBasis d, dd; - SBasis y = (Linear(0, 2) - (a+b)) / (b-a); - SBasis y2 = 2*y; - for(int j = m - 1; j >= 1; j--) { - SBasis sv = d; - d = y2*d - dd + cheb_coeff[j]; - dd = sv; - } - - return y*d - dd + 0.5*cheb_coeff[0]; -} - -SBasis chebyshev_approximant (double (*f)(double,void*), int order, Interval in, void* p) { - gsl_cheb_series *cs = gsl_cheb_alloc (order+2); - - gsl_function F; - - F.function = f; - F.params = p; - - gsl_cheb_init (cs, &F, in[0], in[1]); - - SBasis r = compose(clenshaw_series(order, cs->c), Linear(-1,1)); - - gsl_cheb_free (cs); - return r; -} - -struct wrap { - double (*f)(double,void*); - void* pp; - double fa, fb; - Interval in; -}; - -double f_interp(double x, void* p) { - struct wrap *wr = (struct wrap *)p; - double z = (x - wr->in[0]) / (wr->in[1] - wr->in[0]); - return (wr->f)(x, wr->pp) - ((1 - z)*wr->fa + z*wr->fb); -} - -SBasis chebyshev_approximant_interpolating (double (*f)(double,void*), - int order, Interval in, void* p) { - double fa = f(in[0], p); - double fb = f(in[1], p); - struct wrap wr; - wr.fa = fa; - wr.fb = fb; - wr.in = in; - printf("%f %f\n", fa, fb); - wr.f = f; - wr.pp = p; - return compose(Linear(in[0], in[1]), Linear(fa, fb)) + chebyshev_approximant(f_interp, order, in, &wr) + Linear(fa, fb); -} - -SBasis chebyshev(unsigned n) { - static std::vector basis; - if(basis.empty()) { - basis.push_back(Linear(1,1)); - basis.push_back(Linear(0,1)); - } - for(unsigned i = basis.size(); i <= n; i++) { - basis.push_back(Linear(0,2)*basis[i-1] - basis[i-2]); - } - - return basis[n]; -} - -}; - -/* - 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 : +#include <2geom/chebyshev.h> + +#include <2geom/sbasis.h> +#include <2geom/sbasis-poly.h> + +#include +using std::vector; + +#include +#include + +namespace Geom{ + +SBasis cheb(unsigned n) { + static std::vector basis; + if(basis.empty()) { + basis.push_back(Linear(1,1)); + basis.push_back(Linear(0,1)); + } + for(unsigned i = basis.size(); i <= n; i++) { + basis.push_back(Linear(0,2)*basis[i-1] - basis[i-2]); + } + + return basis[n]; +} + +SBasis cheb_series(unsigned n, double* cheb_coeff) { + SBasis r; + for(unsigned i = 0; i < n; i++) { + double cof = cheb_coeff[i]; + //if(i == 0) + //cof /= 2; + r += cheb(i)*cof; + } + + return r; +} + +SBasis clenshaw_series(unsigned m, double* cheb_coeff) { + /** b_n = a_n + b_n-1 = 2*x*b_n + a_n-1 + b_n-k = 2*x*b_{n-k+1} + a_{n-k} - b_{n - k + 2} + b_0 = x*b_1 + a_0 - b_2 + */ + + double a = -1, b = 1; + SBasis d, dd; + SBasis y = (Linear(0, 2) - (a+b)) / (b-a); + SBasis y2 = 2*y; + for(int j = m - 1; j >= 1; j--) { + SBasis sv = d; + d = y2*d - dd + cheb_coeff[j]; + dd = sv; + } + + return y*d - dd + 0.5*cheb_coeff[0]; +} + +SBasis chebyshev_approximant (double (*f)(double,void*), int order, Interval in, void* p) { + gsl_cheb_series *cs = gsl_cheb_alloc (order+2); + + gsl_function F; + + F.function = f; + F.params = p; + + gsl_cheb_init (cs, &F, in[0], in[1]); + + SBasis r = compose(clenshaw_series(order, cs->c), Linear(-1,1)); + + gsl_cheb_free (cs); + return r; +} + +struct wrap { + double (*f)(double,void*); + void* pp; + double fa, fb; + Interval in; +}; + +double f_interp(double x, void* p) { + struct wrap *wr = (struct wrap *)p; + double z = (x - wr->in[0]) / (wr->in[1] - wr->in[0]); + return (wr->f)(x, wr->pp) - ((1 - z)*wr->fa + z*wr->fb); +} + +SBasis chebyshev_approximant_interpolating (double (*f)(double,void*), + int order, Interval in, void* p) { + double fa = f(in[0], p); + double fb = f(in[1], p); + struct wrap wr; + wr.fa = fa; + wr.fb = fb; + wr.in = in; + printf("%f %f\n", fa, fb); + wr.f = f; + wr.pp = p; + return compose(Linear(in[0], in[1]), Linear(fa, fb)) + chebyshev_approximant(f_interp, order, in, &wr) + Linear(fa, fb); +} + +SBasis chebyshev(unsigned n) { + static std::vector basis; + if(basis.empty()) { + basis.push_back(Linear(1,1)); + basis.push_back(Linear(0,1)); + } + for(unsigned i = basis.size(); i <= n; i++) { + basis.push_back(Linear(0,2)*basis[i-1] - basis[i-2]); + } + + return basis[n]; +} + +}; + +/* + 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/chebyshev.h b/src/2geom/chebyshev.h index 6de9e9cc0..309e09b3e 100644 --- a/src/2geom/chebyshev.h +++ b/src/2geom/chebyshev.h @@ -1,30 +1,30 @@ -#ifndef _CHEBYSHEV -#define _CHEBYSHEV - -#include <2geom/sbasis.h> -#include <2geom/interval.h> - -/*** Conversion between Chebyshev approximation and SBasis. - * - */ - -namespace Geom{ - -SBasis chebyshev_approximant (double (*f)(double,void*), int order, Interval in, void* p=0); -SBasis chebyshev_approximant_interpolating (double (*f)(double,void*), int order, Interval in, void* p=0); -SBasis chebyshev(unsigned n); - -}; - -/* - 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 : - -#endif +#ifndef _CHEBYSHEV +#define _CHEBYSHEV + +#include <2geom/sbasis.h> +#include <2geom/interval.h> + +/*** Conversion between Chebyshev approximation and SBasis. + * + */ + +namespace Geom{ + +SBasis chebyshev_approximant (double (*f)(double,void*), int order, Interval in, void* p=0); +SBasis chebyshev_approximant_interpolating (double (*f)(double,void*), int order, Interval in, void* p=0); +SBasis chebyshev(unsigned n); + +}; + +/* + 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 : + +#endif diff --git a/src/2geom/conjugate_gradient.cpp b/src/2geom/conjugate_gradient.cpp index f5a0f9cd8..9c4ea7776 100644 --- a/src/2geom/conjugate_gradient.cpp +++ b/src/2geom/conjugate_gradient.cpp @@ -46,7 +46,7 @@ matrix_times_vector(valarray const &matrix, /* m * n */ unsigned n = vec.size(); unsigned m = result.size(); assert(m*n == matrix.size()); - const double* mp = &matrix[0]; + const double* mp = &const_cast&>(matrix)[0]; for (unsigned i = 0; i < m; i++) { double res = 0; for (unsigned j = 0; j < n; j++) diff --git a/src/2geom/crossing.h b/src/2geom/crossing.h index 7737b24a9..427848033 100644 --- a/src/2geom/crossing.h +++ b/src/2geom/crossing.h @@ -110,10 +110,15 @@ struct NearF { bool operator()(Crossing a, Crossing b) { return are_near(a, b); struct CrossingOrder { unsigned ix; - CrossingOrder(unsigned i) : ix(i) {} + bool rev; + CrossingOrder(unsigned i, bool r = false) : ix(i), rev(r) {} bool operator()(Crossing a, Crossing b) { - return (ix == a.a ? a.ta : a.tb) < - (ix == b.a ? b.ta : b.tb); + if(rev) + return (ix == a.a ? a.ta : a.tb) < + (ix == b.a ? b.ta : b.tb); + else + return (ix == a.a ? a.ta : a.tb) > + (ix == a.a ? a.ta : a.tb); } }; diff --git a/src/2geom/d2.h b/src/2geom/d2.h index f55f6a596..afa00b40d 100644 --- a/src/2geom/d2.h +++ b/src/2geom/d2.h @@ -397,6 +397,14 @@ D2 integral(D2 const & a) { return D2(integral(a[X]), integral(a[Y])); } +/** A function to print out the Point. It just prints out the coords + on the given output stream */ +template +inline std::ostream &operator<< (std::ostream &out_file, const Geom::D2 &in_d2) { + out_file << "X: " << in_d2[X] << " Y: " << in_d2[Y]; + return out_file; +} + } //end namespace Geom #include <2geom/rect.h> diff --git a/src/2geom/line.h b/src/2geom/line.h index e04c35c29..a7e6a54bb 100644 --- a/src/2geom/line.h +++ b/src/2geom/line.h @@ -79,6 +79,12 @@ class Line return Line(P, P+rot90(n)); } + static Line fromPointDirection(Point o, Point v) { + Line l; + l.m_origin = o; + l.m_versor = v; + return l; + } Line* duplicate() const { diff --git a/src/2geom/numeric/fitting-tool.h b/src/2geom/numeric/fitting-tool.h index 47c588cb2..d6a26bd2d 100644 --- a/src/2geom/numeric/fitting-tool.h +++ b/src/2geom/numeric/fitting-tool.h @@ -463,7 +463,17 @@ class least_squeares_fitter typedef typename base_type::solution_type solution_type; using base_type::m_vector_view; - using base_type::result; + //using base_type::result; // VSC legacy support + solution_type& result( std::vector const& old_sample_values, + std::vector const& new_sample_values ) + { + return base_type::result(old_sample_values, new_sample_values); + } + + solution_type& result() + { + return base_type::result(); + } public: least_squeares_fitter( model_type const& _model, @@ -497,7 +507,18 @@ class least_squeares_fitter typedef typename base_type::solution_type solution_type; using base_type::m_vector_view; - using base_type::result; + //using base_type::result; // VCS legacy support + solution_type& result( std::vector const& old_sample_values, + std::vector const& new_sample_values ) + { + return base_type::result(old_sample_values, new_sample_values); + } + + solution_type& result() + { + return base_type::result(); + } + public: least_squeares_fitter( model_type const& _model, diff --git a/src/2geom/numeric/matrix.h b/src/2geom/numeric/matrix.h index e9b6e6e9e..97db59d56 100644 --- a/src/2geom/numeric/matrix.h +++ b/src/2geom/numeric/matrix.h @@ -232,11 +232,15 @@ class MatrixImpl : public BaseMatrixImpl gsl_matrix_set_identity(m_matrix); } - using base_type::operator(); + using base_type::operator(); // VSC legacy support + const double & operator() (size_t i, size_t j) const + { + return base_type::operator ()(i, j); + } - double & operator() (size_t i, size_t j) + double & operator() (size_t i, size_t j) { - return *gsl_matrix_ptr(m_matrix, i, j); + return *gsl_matrix_ptr(m_matrix, i, j); } using base_type::get_gsl_matrix; diff --git a/src/2geom/path-intersection.cpp b/src/2geom/path-intersection.cpp index bfb1fb96c..b2d5ceabb 100644 --- a/src/2geom/path-intersection.cpp +++ b/src/2geom/path-intersection.cpp @@ -4,8 +4,11 @@ //for path_direction: #include <2geom/sbasis-geometric.h> +#include <2geom/line.h> +#ifdef HAVE_GSL #include #include +#endif namespace Geom { @@ -199,6 +202,7 @@ static double EpsilonOf(double value) } #endif +#ifdef HAVE_GSL struct rparams { Curve const &A; Curve const &B; @@ -219,45 +223,91 @@ intersect_polish_f (const gsl_vector * x, void *params, return GSL_SUCCESS; } +#endif static void intersect_polish_root (Curve const &A, double &s, Curve const &B, double &t) { int status; size_t iter = 0; + std::vector as, bs; + as = A.pointAndDerivatives(s, 2); + bs = B.pointAndDerivatives(t, 2); + Point F = as[0] - bs[0]; + double best = dot(F, F); + + for(int i = 0; i < 4; i++) { + + /** + we want to solve + J*(x1 - x0) = f(x0) + + |dA(s)[0] -dB(t)[0]| (X1 - X0) = A(s) - B(t) + |dA(s)[1] -dB(t)[1]| + **/ + + // We're using the standard transformation matricies, which is numerically rather poor. Much better to solve the equation using elimination. + + Matrix jack(as[1][0], as[1][1], + -bs[1][0], -bs[1][1], + 0, 0); + Point soln = (F)*jack.inverse(); + double ns = s - soln[0]; + double nt = t - soln[1]; + + if (ns<0) ns=0; + else if (ns>1) ns=1; + if (nt<0) nt=0; + else if (nt>1) nt=1; + + as = A.pointAndDerivatives(ns, 2); + bs = B.pointAndDerivatives(nt, 2); + F = as[0] - bs[0]; + double trial = dot(F, F); + if (trial > best*0.1) // we have standards, you know + // At this point we could do a line search + break; + best = trial; + s = ns; + t = nt; + } + +#ifdef HAVE_GSL + if(0) { // the GSL version is more accurate, but taints this with GPL + const size_t n = 2; + struct rparams p = {A, B}; + gsl_multiroot_function f = {&intersect_polish_f, n, &p}; - 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); + 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]); + gsl_vector_set (x, 0, x_init[0]); + gsl_vector_set (x, 1, x_init[1]); - const gsl_multiroot_fsolver_type *T = gsl_multiroot_fsolver_hybrids; - gsl_multiroot_fsolver *sol = gsl_multiroot_fsolver_alloc (T, 2); - gsl_multiroot_fsolver_set (sol, &f, x); + const gsl_multiroot_fsolver_type *T = gsl_multiroot_fsolver_hybrids; + gsl_multiroot_fsolver *sol = gsl_multiroot_fsolver_alloc (T, 2); + gsl_multiroot_fsolver_set (sol, &f, x); - do - { - iter++; - status = gsl_multiroot_fsolver_iterate (sol); + do + { + iter++; + status = gsl_multiroot_fsolver_iterate (sol); - if (status) /* check if solver is stuck */ - break; + if (status) /* check if solver is stuck */ + break; - status = - gsl_multiroot_test_residual (sol->f, 1e-12); - } - while (status == GSL_CONTINUE && iter < 1000); + 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); + s = gsl_vector_get (sol->x, 0); + t = gsl_vector_get (sol->x, 1); - gsl_multiroot_fsolver_free (sol); - gsl_vector_free (x); + gsl_multiroot_fsolver_free (sol); + gsl_vector_free (x); + } +#endif } /** @@ -267,7 +317,7 @@ intersect_polish_root (Curve const &A, double &s, */ void pair_intersect(Curve const & A, double Al, double Ah, Curve const & B, double Bl, double Bh, - Crossings &ret, unsigned depth=0) { + Crossings &ret, unsigned depth = 0) { // std::cout << depth << "(" << Al << ", " << Ah << ")\n"; OptRect Ar = A.boundsLocal(Interval(Al, Ah)); if (!Ar) return; @@ -305,6 +355,13 @@ void pair_intersect(Curve const & A, double Al, double Ah, ret, depth+1); } +Crossings pair_intersect(Curve const & A, Interval const &Ad, + Curve const & B, Interval const &Bd) { + Crossings ret; + pair_intersect(A, Ad.min(), Ad.max(), B, Bd.min(), Bd.max(), ret); + return ret; +} + /** A simple wrapper around pair_intersect */ Crossings SimpleCrosser::crossings(Curve const &a, Curve const &b) { Crossings ret; @@ -312,6 +369,53 @@ Crossings SimpleCrosser::crossings(Curve const &a, Curve const &b) { return ret; } + +//same as below but curves not paths +void mono_intersect(Curve const &A, double Al, double Ah, + Curve const &B, double Bl, double Bh, + Crossings &ret, double tol = 0.1, unsigned depth = 0) { + if( Al >= Ah || Bl >= Bh) return; + //std::cout << " " << depth << "[" << Al << ", " << Ah << "]" << "[" << Bl << ", " << Bh << "]"; + + Point A0 = A.pointAt(Al), A1 = A.pointAt(Ah), + B0 = B.pointAt(Bl), B1 = B.pointAt(Bh); + //inline code that this implies? (without rect/interval construction) + Rect Ar = Rect(A0, A1), Br = Rect(B0, B1); + if(!Ar.intersects(Br) || A0 == A1 || B0 == B1) return; + + if(depth > 12 || (Ar.maxExtent() < tol && Ar.maxExtent() < tol)) { + double tA, tB, c; + if(linear_intersect(A.pointAt(Al), A.pointAt(Ah), + B.pointAt(Bl), B.pointAt(Bh), + tA, tB, c)) { + tA = tA * (Ah - Al) + Al; + tB = tB * (Bh - Bl) + Bl; + intersect_polish_root(A, tA, + B, tB); + if(depth % 2) + ret.push_back(Crossing(tB, tA, c < 0)); + else + ret.push_back(Crossing(tA, tB, c > 0)); + return; + } + } + if(depth > 12) return; + double mid = (Bl + Bh)/2; + mono_intersect(B, Bl, mid, + A, Al, Ah, + ret, tol, depth+1); + mono_intersect(B, mid, Bh, + A, Al, Ah, + ret, tol, depth+1); +} + +Crossings mono_intersect(Curve const & A, Interval const &Ad, + Curve const & B, Interval const &Bd) { + Crossings ret; + mono_intersect(A, Ad.min(), Ad.max(), B, Bd.min(), Bd.max(), ret); + return ret; +} + /** * Takes two paths and time ranges on them, with the invariant that the * paths are monotonic on the range. Splits A when the linear intersection @@ -327,11 +431,10 @@ void mono_pair(Path const &A, double Al, double Ah, Point A0 = A.pointAt(Al), A1 = A.pointAt(Ah), B0 = B.pointAt(Bl), B1 = B.pointAt(Bh); //inline code that this implies? (without rect/interval construction) - if(!Rect(A0, A1).intersects(Rect(B0, B1)) || A0 == A1 || B0 == B1) return; - - //Checks the general linearity of the function - //if((depth > 12) || (A.boundsLocal(Interval(Al, Ah), 1).maxExtent() < 0.1 - // && B.boundsLocal(Interval(Bl, Bh), 1).maxExtent() < 0.1)) { + Rect Ar = Rect(A0, A1), Br = Rect(B0, B1); + if(!Ar.intersects(Br) || A0 == A1 || B0 == B1) return; + + if(depth > 12 || (Ar.maxExtent() < 0.1 && Ar.maxExtent() < 0.1)) { double tA, tB, c; if(linear_intersect(A0, A1, B0, B1, tA, tB, c)) { @@ -343,7 +446,7 @@ void mono_pair(Path const &A, double Al, double Ah, ret.push_back(Crossing(tA, tB, c > 0)); return; } - //} + } if(depth > 12) return; double mid = (Bl + Bh)/2; mono_pair(B, Bl, mid, @@ -356,8 +459,10 @@ void mono_pair(Path const &A, double Al, double Ah, /** This returns the times when the x or y derivative is 0 in the curve. */ std::vector curve_mono_splits(Curve const &d) { + Curve* deriv = d.derivative(); std::vector rs = d.roots(0, X); append(rs, d.roots(0, Y)); + delete deriv; std::sort(rs.begin(), rs.end()); return rs; } @@ -378,17 +483,10 @@ std::vector offset_doubles(std::vector const &x, double offs) { std::vector path_mono_splits(Path const &p) { std::vector ret; if(p.empty()) return ret; - ret.push_back(0); - - Curve* deriv = p[0].derivative(); - append(ret, curve_mono_splits(*deriv)); - delete deriv; bool pdx=2, pdy=2; //Previous derivative direction for(unsigned i = 0; i <= p.size(); i++) { - deriv = p[i].derivative(); - std::vector spl = offset_doubles(curve_mono_splits(*deriv), i); - delete deriv; + std::vector spl = offset_doubles(curve_mono_splits(p[i]), i); bool dx = p[i].initialPoint()[X] > (spl.empty()? p[i].finalPoint()[X] : p.valueAt(spl.front(), X)); bool dy = p[i].initialPoint()[Y] > (spl.empty()? p[i].finalPoint()[Y] : diff --git a/src/2geom/path-intersection.h b/src/2geom/path-intersection.h index 4eef16823..6457b5e43 100644 --- a/src/2geom/path-intersection.h +++ b/src/2geom/path-intersection.h @@ -66,6 +66,11 @@ Crossings curve_sweep(Path const &a, Path const &b) { return ret; } +Crossings pair_intersect(Curve const & A, Interval const &Ad, + Curve const & B, Interval const &Bd); +Crossings mono_intersect(Curve const & A, Interval const &Ad, + Curve const & B, Interval const &Bd); + struct SimpleCrosser : public Crosser { Crossings crossings(Curve const &a, Curve const &b); Crossings crossings(Path const &a, Path const &b) { return curve_sweep(a, b); } diff --git a/src/2geom/path.cpp b/src/2geom/path.cpp index 136e6d481..981c9f044 100644 --- a/src/2geom/path.cpp +++ b/src/2geom/path.cpp @@ -290,7 +290,9 @@ double Path::nearestPoint(Point const &_point, double from, double to, double *d } void Path::appendPortionTo(Path &ret, double from, double to) const { - assert(from >= 0 && to >= 0); + if (!(from >= 0 && to >= 0)) { + THROW_RANGEERROR("from and to must be >=0 in Path::appendPortionTo"); + } if(to == 0) to = size()+0.999999; if(from == to) { return; } double fi, ti; diff --git a/src/2geom/piecewise.h b/src/2geom/piecewise.h index 144d7a9b2..a5be42587 100644 --- a/src/2geom/piecewise.h +++ b/src/2geom/piecewise.h @@ -531,6 +531,19 @@ Piecewise operator*(Piecewise const &a, double b) { return ret; } template +Piecewise operator*(Piecewise const &a, T b) { + boost::function_requires >(); + + if(a.empty()) return Piecewise(); + + Piecewise ret; + ret.segs.reserve(a.size()); + ret.cuts = a.cuts; + for(unsigned i = 0; i < a.size();i++) + ret.push_seg(a[i] * b); + return ret; +} +template Piecewise operator/(Piecewise const &a, double b) { boost::function_requires >(); @@ -753,10 +766,10 @@ Piecewise reverse(Piecewise const &f) { double end = f.cuts.back(); for (unsigned i = 0; i < f.cuts.size(); i++) { double x = f.cuts[f.cuts.size() - 1 - i]; - ret.cuts[i] = end - (x - start); + ret.push_cut(end - (x - start)); } for (unsigned i = 0; i < f.segs.size(); i++) - ret.segs[i] = reverse(f[f.segs.size() - i - 1]); + ret.push_seg(reverse(f[f.segs.size() - i - 1])); return ret; } diff --git a/src/2geom/poly.cpp b/src/2geom/poly.cpp index 30ec1cdb7..d8b379557 100644 --- a/src/2geom/poly.cpp +++ b/src/2geom/poly.cpp @@ -1,8 +1,9 @@ #include <2geom/poly.h> -//#ifdef HAVE_GSL +#define HAVE_GSL +#ifdef HAVE_GSL #include -//#endif +#endif namespace Geom { @@ -38,7 +39,7 @@ void Poly::monicify() { } -//#ifdef HAVE_GSL +#ifdef HAVE_GSL std::vector > solve(Poly const & pp) { Poly p(pp); p.normalize(); @@ -75,7 +76,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/quadtree.cpp b/src/2geom/quadtree.cpp index e322a091b..211590bae 100644 --- a/src/2geom/quadtree.cpp +++ b/src/2geom/quadtree.cpp @@ -43,9 +43,94 @@ Quad* QuadTree::search(double x0, double y0, double x1, double y1) { } return q; } + + +/* +Comments by Vangelis (use with caution :P ) + +Insert Rect (x0, y0), (x1, y1) in the QuadTree Q. + +=================================================================================== +* QuadTree Q has: Quadtree's Quad root R, QuadTree's bounding box B. + +* Each Quad has a Quad::data where we store the id of the Rect that belong to +this Quad. (In reality we'll store a pointer to the shape). + +* Each Quad has 4 Quad children: 0, 1, 2, 3. Each child Quad represents one of the following quarters +of the bounding box B: + ++---------------------+ +| | | +| NW=0 | NE=1 | +| | | +| | | ++---------------------+ +| | | +| SW=2 | SE=3 | +| | | +| | | ++---------------------+ + +Each Quad can further be divided in 4 Quads as above and so on. Below there is an example + + Root + / || \ + / / \ \ + 0 1 2 3 + /\ + / | | \ + 0 1 2 3 + ++---------------------+ +| | 1-0 | 1-1| +| 0 | | | +| |-----|----| +| | 1-2 | 1-3| +| | | | ++---------------------+ +| | | +| | | +| 2 | 3 | +| | | ++---------------------+ + + + +=================================================================================== +Insert Rect (x0, y0), (x1, y1) in the QuadTree Q. Algorithm: +1) check if Rect is bigger than QuadTree's bounding box +2) find in which Quad we should add the Rect: + + + +----------------------------------------------------------------------------------- +How we find in which Quad we should add the Rect R: + +Q = Quadtree's Quad root +B = QuadTree's bounding box B +WHILE (Q) { + IF ( Rect cannot fit in one unique quarter of B ){ + Q = current Quad ; + BREAK; + } + IF ( Rect can fit in the quarter I ) { + IF (Q.children[I] doesn't exist) { + create the Quad Q.children[I]; + } + B = bounding box of the Quad Q.children[I] ; + Q = Q.children[I] ; + CHECK(R, B) ; + } +} +add Rect R to Q ; + + +*/ void QuadTree::insert(double x0, double y0, double x1, double y1, int shape) { // loop until a quad would break the box. + + // empty root => empty QuadTree. Create initial bounding box (0,0), (1,1) if(root == 0) { root = new Quad; @@ -55,14 +140,18 @@ void QuadTree::insert(double x0, double y0, double x1, double y1, int shape) { by1 = 1; } Quad *q = root; - + + //A temp bounding box. Same as root's bounting box (ie of the whole QuadTree) double bxx0 = bx0, bxx1 = bx1; double byy0 = by0, byy1 = by1; + while((bxx0 > x0) || (bxx1 < x1) || (byy0 > y0) || - (byy1 < y1)) { // too small initial size - double + (byy1 < y1)) { + // QuadTree has small size, can't accomodate new rect. Double the size: unsigned i = 0; + if(bxx0 > x0) { bxx0 = 2*bxx0 - bxx1; i += 1; @@ -76,15 +165,22 @@ void QuadTree::insert(double x0, double y0, double x1, double y1, int shape) { byy1 = 2*byy1 - byy0; } q = new Quad; - q->children[i] = root; - root = q; + //check if root is empty (no rects, no quad children) + if( clean_root() ){ + root = q; + } + else{ + q->children[i] = root; + root = q; + } bx0 = bxx0; bx1 = bxx1; by0 = byy0; by1 = byy1; } - + while(q) { + // Find the center of the temp bounding box double cx = (bxx0 + bxx1)/2; double cy = (byy0 + byy1)/2; unsigned i = 0; @@ -92,21 +188,41 @@ void QuadTree::insert(double x0, double y0, double x1, double y1, int shape) { assert(x1 <= bxx1); assert(y0 >= byy0); assert(y1 <= byy1); + if(x0 >= cx) { i += 1; bxx0 = cx; // zoom in a quad } else if(x1 <= cx) { bxx1 = cx; - } else + } else{ + // rect does not fit in one unique quarter (in X axis) of the temp bounding box break; + } if(y0 >= cy) { i += 2; byy0 = cy; } else if(y1 <= cy) { byy1 = cy; - } else + } else{ + // rect does not fit in one unique quarter (in Y axis) of the temp bounding box break; - + } + + // check if rect's bounding box has size 1x1. This means that rect is defined by 2 points + // that are in the same place. + if( ( fabs(bxx0 - bxx1) < 1.0 ) && ( fabs(byy0 - byy1) < 1.0 )){ + bxx0 = floor(bxx0); + bxx1 = floor(bxx1); + byy0 = floor(byy0); + byy1 = floor(byy1); + break; + } + + /* + 1 rect does fit in one unique quarter of the temp bounding box. And we have found which. + 2 temp bounding box = bounding box of this quarter. + 3 "Go in" this quarter (create if doesn't exist) + */ assert(i < 4); Quad *qq = q->children[i]; if(qq == 0) { @@ -129,6 +245,35 @@ void QuadTree::erase(Quad *q, int shape) { return; } +/* +Returns: +false: if root isn't empty +true: if root is empty it cleans root +*/ +bool QuadTree::clean_root() { + assert(root); + + // false if root *has* rects assigned to it. + bool all_clean = root->data.empty(); + + // if root has children we get false + for(unsigned i = 0; i < 4; i++) + { + if(root->children[i]) + { + all_clean = false; + } + } + + if(all_clean) + { + delete root; + root=0; + return true; + } + return false; +} + }; /* diff --git a/src/2geom/quadtree.h b/src/2geom/quadtree.h index 5338698cf..2e114a0a0 100644 --- a/src/2geom/quadtree.h +++ b/src/2geom/quadtree.h @@ -81,6 +81,8 @@ public: Quad* search(Geom::Rect const &r); void insert(Geom::Rect const &r, int shape); void erase(Quad *q, int shape); +private: + bool clean_root(); }; }; diff --git a/src/2geom/sbasis-math.cpp b/src/2geom/sbasis-math.cpp index 7053388dc..f3a984c96 100644 --- a/src/2geom/sbasis-math.cpp +++ b/src/2geom/sbasis-math.cpp @@ -149,7 +149,7 @@ static Piecewise sqrt_internal(SBasis const &f, sqrtf.resize(order+1, Linear(0,0)); sqrtf[0] = Linear(std::sqrt(f[0][0]), std::sqrt(f[0][1])); SBasis r = f - multiply(sqrtf, sqrtf); // remainder - for(unsigned i = 1; int(i) <= order and i roots1(SBasis const & s) { double d = s[0][0] - s[0][1]; if(d != 0) { double r = s[0][0] / d; - if(0 <= r and r <= 1) + if(0 <= r && r <= 1) res.push_back(r); } return res; diff --git a/src/2geom/sbasis.cpp b/src/2geom/sbasis.cpp index 49678b19a..2f7f03bfc 100644 --- a/src/2geom/sbasis.cpp +++ b/src/2geom/sbasis.cpp @@ -375,7 +375,7 @@ SBasis sqrt(SBasis const &a, int k) { c[0] = Linear(std::sqrt(a[0][0]), std::sqrt(a[0][1])); SBasis r = a - multiply(c, c); // remainder - for(unsigned i = 1; i <= (unsigned)k and i::iterator(before), n, l);} - bool operator==(SBasis const&B) { return d == B.d;} + bool operator==(SBasis const&B) const { return d == B.d;} + bool operator!=(SBasis const&B) const { return d != B.d;} operator std::vector() { return d;} @@ -99,6 +100,9 @@ public: explicit SBasis(double a) { push_back(Linear(a,a)); } + explicit SBasis(double a, double b) { + push_back(Linear(a,b)); + } SBasis(SBasis const & a) : d(a.d) {} diff --git a/src/2geom/solve-bezier-one-d.cpp b/src/2geom/solve-bezier-one-d.cpp index 9a7b5633f..a1c0ca557 100644 --- a/src/2geom/solve-bezier-one-d.cpp +++ b/src/2geom/solve-bezier-one-d.cpp @@ -1,6 +1,7 @@ #include <2geom/solver.h> #include <2geom/point.h> #include +#include /*** Find the zeros of the bernstein function. The code subdivides until it is happy with the * linearity of the function. This requires an O(degree^2) subdivision for each step, even when @@ -134,8 +135,7 @@ Bernsteins::find_bernstein_roots(double const *w, /* The control points */ } /* Otherwise, solve recursively after subdividing control polygon */ - double Left[N], /* New left and right */ - Right[N]; /* control polygons */ + std::valarray new_controls(2*N); // New left and right control polygons const double t = 0.5; @@ -150,28 +150,28 @@ Bernsteins::find_bernstein_roots(double const *w, /* The control points */ /* Triangle computation */ const double omt = (1-t); - Left[0] = Vtemp[0]; - Right[degree] = Vtemp[degree]; + new_controls[0] = Vtemp[0]; + new_controls[N+degree] = Vtemp[degree]; double *prev_row = Vtemp; double *row = Vtemp + N; for (unsigned i = 1; i < N; i++) { for (unsigned j = 0; j < N - i; j++) { row[j] = omt*prev_row[j] + t*prev_row[j+1]; } - Left[i] = row[0]; - Right[degree-i] = row[degree-i]; + new_controls[i] = row[0]; + new_controls[N+degree-i] = row[degree-i]; std::swap(prev_row, row); } double mid_t = left_t*(1-t) + right_t*t; - find_bernstein_roots(Left, depth+1, left_t, mid_t); + find_bernstein_roots(&new_controls[0], depth+1, left_t, mid_t); /* Solution is exactly on the subdivision point. */ - if (Right[0] == 0) + if (new_controls[N] == 0) solutions.push_back(mid_t); - find_bernstein_roots(Right, depth+1, mid_t, right_t); + find_bernstein_roots(&new_controls[N], depth+1, mid_t, right_t); } /* diff --git a/src/2geom/sweep.cpp b/src/2geom/sweep.cpp index 227c47e09..7571efe09 100644 --- a/src/2geom/sweep.cpp +++ b/src/2geom/sweep.cpp @@ -8,18 +8,19 @@ namespace Geom { * \brief Make a list of pairs of self intersections in a list of Rects. * * \param rs: vector of Rect. + * \param d: dimension to sweep along * * [(A = rs[i], B = rs[j]) for i,J in enumerate(pairs) for j in J] * then A.left <= B.left */ -std::vector > sweep_bounds(std::vector rs) { +std::vector > sweep_bounds(std::vector rs, Dim2 d) { std::vector events; events.reserve(rs.size()*2); std::vector > pairs(rs.size()); for(unsigned i = 0; i < rs.size(); i++) { - events.push_back(Event(rs[i].left(), i, false)); - events.push_back(Event(rs[i].right(), i, true)); + events.push_back(Event(rs[i][d][0], i, false)); + events.push_back(Event(rs[i][d][1], i, true)); } std::sort(events.begin(), events.end()); @@ -33,7 +34,7 @@ std::vector > sweep_bounds(std::vector rs) { } else { for(unsigned j = 0; j < open.size(); j++) { unsigned jx = open[j]; - if(rs[jx][Y].intersects(rs[ix][Y])) { + if(rs[jx][1-d].intersects(rs[ix][1-d])) { pairs[jx].push_back(ix); } } @@ -48,11 +49,12 @@ std::vector > sweep_bounds(std::vector rs) { * * \param a: vector of Rect. * \param b: vector of Rect. + * \param d: dimension to scan along * * [(A = rs[i], B = rs[j]) for i,J in enumerate(pairs) for j in J] * then A.left <= B.left, A in a, B in b */ -std::vector > sweep_bounds(std::vector a, std::vector b) { +std::vector > sweep_bounds(std::vector a, std::vector b, Dim2 d) { std::vector > pairs(a.size()); if(a.empty() || b.empty()) return pairs; std::vector events[2]; @@ -64,15 +66,17 @@ std::vector > sweep_bounds(std::vector a, std::vecto events[n].reserve(sz*2); for(unsigned i = 0; i < sz; i++) { Rect r = n ? b[i] : a[i]; - events[n].push_back(Event(r.left(), i, false)); - events[n].push_back(Event(r.right(), i, true)); + events[n].push_back(Event(r[d][0], i, false)); + events[n].push_back(Event(r[d][1], i, true)); } std::sort(events[n].begin(), events[n].end()); } std::vector open[2]; bool n = events[1].front() < events[0].front(); - for(unsigned i[] = {0,0}; i[n] < events[n].size();) { + {// As elegant as putting the initialiser in the for was, it upsets some legacy compilers (MS VS C++) + unsigned i[] = {0,0}; + for(; i[n] < events[n].size();) { unsigned ix = events[n][i[n]].ix; bool closing = events[n][i[n]].closing; //std::cout << n << "[" << ix << "] - " << (closing ? "closer" : "opener") << "\n"; @@ -84,7 +88,7 @@ std::vector > sweep_bounds(std::vector a, std::vecto //opening a B, add to all open a for(unsigned j = 0; j < open[0].size(); j++) { unsigned jx = open[0][j]; - if(a[jx][Y].intersects(b[ix][Y])) { + if(a[jx][1-d].intersects(b[ix][1-d])) { pairs[jx].push_back(ix); } } @@ -93,7 +97,7 @@ std::vector > sweep_bounds(std::vector a, std::vecto //opening an A, add all open b for(unsigned j = 0; j < open[1].size(); j++) { unsigned jx = open[1][j]; - if(b[jx][Y].intersects(a[ix][Y])) { + if(b[jx][1-d].intersects(a[ix][1-d])) { pairs[ix].push_back(jx); } } @@ -103,7 +107,7 @@ std::vector > sweep_bounds(std::vector a, std::vecto i[n]++; if(i[n]>=events[n].size()) {break;} n = (events[!n][i[!n]] < events[n][i[n]]) ? !n : n; - } + }} return pairs; } diff --git a/src/2geom/sweep.h b/src/2geom/sweep.h index 0214511ac..9d1643d7a 100644 --- a/src/2geom/sweep.h +++ b/src/2geom/sweep.h @@ -51,10 +51,12 @@ struct Event { if(x > other.x) return false; return closing < other.closing; } - }; -std::vector > sweep_bounds(std::vector); -std::vector > sweep_bounds(std::vector, std::vector); +std::vector > +sweep_bounds(std::vector, Dim2 dim = X); + +std::vector > +sweep_bounds(std::vector, std::vector, Dim2 dim = X); std::vector > fake_cull(unsigned a, unsigned b); diff --git a/src/2geom/utils.cpp b/src/2geom/utils.cpp index 579718553..f0f42ce01 100644 --- a/src/2geom/utils.cpp +++ b/src/2geom/utils.cpp @@ -1,87 +1,87 @@ -/** Various utility functions. - * - * Copyright 2008 Marco Cecchetti - * Copyright 2007 Johan Engelen - * Copyright 2006 Michael G. Sloan - * - * 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/utils.h> - - -namespace Geom -{ - -// return a vector that contains all the binomial coefficients of degree n -void binomial_coefficients(std::vector& bc, size_t n) -{ - size_t s = n+1; - bc.clear(); - bc.resize(s); - bc[0] = 1; - size_t k; - for (size_t i = 1; i < n; ++i) - { - k = i >> 1; - if (i & 1u) - { - bc[k+1] = bc[k] << 1; - } - for (size_t j = k; j > 0; --j) - { - bc[j] += bc[j-1]; - } - } - s >>= 1; - for (size_t i = 0; i < s; ++i) - { - bc[n-i] = bc[i]; - } -} - -} // 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 : +/** Various utility functions. + * + * Copyright 2008 Marco Cecchetti + * Copyright 2007 Johan Engelen + * Copyright 2006 Michael G. Sloan + * + * 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/utils.h> + + +namespace Geom +{ + +// return a vector that contains all the binomial coefficients of degree n +void binomial_coefficients(std::vector& bc, size_t n) +{ + size_t s = n+1; + bc.clear(); + bc.resize(s); + bc[0] = 1; + size_t k; + for (size_t i = 1; i < n; ++i) + { + k = i >> 1; + if (i & 1u) + { + bc[k+1] = bc[k] << 1; + } + for (size_t j = k; j > 0; --j) + { + bc[j] += bc[j-1]; + } + } + s >>= 1; + for (size_t i = 0; i < s; ++i) + { + bc[n-i] = bc[i]; + } +} + +} // 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 :