![[tokkee]](http://tokkee.org/images/avatar.png){kind=avatar}
summary | shortlog | log | commit | commitdiff | tree
raw | patch | inline | side by side (parent: 884edd7)
raw | patch | inline | side by side (parent: 884edd7)
| author | johanengelen <johanengelen@users.sourceforge.net> | |
| Mon, 6 Apr 2009 22:29:34 +0000 (22:29 +0000) | ||
| committer | johanengelen <johanengelen@users.sourceforge.net> | |
| Mon, 6 Apr 2009 22:29:34 +0000 (22:29 +0000) |
i turned some optional compilation stuff *on* per default, to help building inkscape.
27 files changed:
index 0c8b96ddff3084f28449044dfb253889172f8aa0..e159839d208626c77a30de4d1c93f9c4bb66ea33 100644 (file)
#include <2geom/sbasis-to-bezier.h>
#include <2geom/exception.h>
-
+#ifdef HAVE_GSL
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multiroots.h>
-
+#endif
using std::vector;
namespace Geom {
//unique(xs.begin(), xs.end());
}
-
+#ifdef HAVE_GSL
#include <gsl/gsl_multiroots.h>
-
-
struct rparams
{
return GSL_SUCCESS;
}
+#endif
union dbl_64{
long long i64;
static void intersect_polish_root (D2<SBasis> const &A, double &s,
D2<SBasis> 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<Point> 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};
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;
best_v = trial_v;
}
}
+ }
}
index 96a06376cb43a1b9a2299f2bf27e08290fc0d75d..c8d99c4308fbcdb22c7918627d774205b0144ab3 100644 (file)
-/*
- * Implement the Bezier clipping algorithm for finding
- * Bezier curve intersection points and collinear normals
- *
- * Authors:
- * Marco Cecchetti <mrcekets at gmail.com>
- *
- * Copyright 2008 authors
- *
- * This library is free software; you can redistribute it and/or
- * modify it either under the terms of the GNU Lesser General Public
- * License version 2.1 as published by the Free Software Foundation
- * (the "LGPL") or, at your option, under the terms of the Mozilla
- * Public License Version 1.1 (the "MPL"). If you do not alter this
- * notice, a recipient may use your version of this file under either
- * the MPL or the LGPL.
- *
- * You should have received a copy of the LGPL along with this library
- * in the file COPYING-LGPL-2.1; if not, write to the Free Software
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
- * You should have received a copy of the MPL along with this library
- * in the file COPYING-MPL-1.1
- *
- * The contents of this file are subject to the Mozilla Public License
- * Version 1.1 (the "License"); you may not use this file except in
- * compliance with the License. You may obtain a copy of the License at
- * http://www.mozilla.org/MPL/
- *
- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
- * OF ANY KIND, either express or implied. See the LGPL or the MPL for
- * the specific language governing rights and limitations.
- */
-
-
-
-
-#include <2geom/basic-intersection.h>
-#include <2geom/choose.h>
-#include <2geom/point.h>
-#include <2geom/interval.h>
-#include <2geom/bezier.h>
-//#include <2geom/convex-cover.h>
-#include <2geom/numeric/matrix.h>
-
-#include <cassert>
-#include <vector>
-#include <algorithm>
-#include <utility>
-//#include <iomanip>
-
-
-
-
-#define VERBOSE 0
-#define CHECK 0
-
-
-namespace Geom {
-
-namespace detail { namespace bezier_clipping {
-
-////////////////////////////////////////////////////////////////////////////////
-// for debugging
-//
-
-inline
-void print(std::vector<Point> const& cp, const char* msg = "")
-{
- std::cerr << msg << std::endl;
- for (size_t i = 0; i < cp.size(); ++i)
- std::cerr << i << " : " << cp[i] << std::endl;
-}
-
-template< class charT >
-inline
-std::basic_ostream<charT> &
-operator<< (std::basic_ostream<charT> & os, const Interval & I)
-{
- os << "[" << I.min() << ", " << I.max() << "]";
- return os;
-}
-
-inline
-double angle (std::vector<Point> const& A)
-{
- size_t n = A.size() -1;
- double a = std::atan2(A[n][Y] - A[0][Y], A[n][X] - A[0][X]);
- return (180 * a / M_PI);
-}
-
-inline
-size_t get_precision(Interval const& I)
-{
- double d = I.extent();
- double e = 0.1, p = 10;
- int n = 0;
- while (n < 16 && d < e)
- {
- p *= 10;
- e = 1/p;
- ++n;
- }
- return n;
-}
-
-inline
-void range_assertion(int k, int m, int n, const char* msg)
-{
- if ( k < m || k > n)
- {
- std::cerr << "range assertion failed: \n"
- << msg << std::endl
- << "value: " << k
- << " range: " << m << ", " << n << std::endl;
- assert (k >= m && k <= n);
- }
-}
-
-
-////////////////////////////////////////////////////////////////////////////////
-// convex hull
-
-/*
- * return true in case the oriented polyline p0, p1, p2 is a right turn
- */
-inline
-bool is_a_right_turn (Point const& p0, Point const& p1, Point const& p2)
-{
- if (p1 == p2) return false;
- Point q1 = p1 - p0;
- Point q2 = p2 - p0;
- if (q1 == -q2) return false;
- return (cross (q1, q2) < 0);
-}
-
-/*
- * return true if p < q wrt the lexicographyc order induced by the coordinates
- */
-struct lex_less
-{
- bool operator() (Point const& p, Point const& q)
- {
- return ((p[X] < q[X]) || (p[X] == q[X] && p[Y] < q[Y]));
- }
-};
-
-/*
- * return true if p > q wrt the lexicographyc order induced by the coordinates
- */
-struct lex_greater
-{
- bool operator() (Point const& p, Point const& q)
- {
- return ((p[X] > q[X]) || (p[X] == q[X] && p[Y] > q[Y]));
- }
-};
-
-/*
- * Compute the convex hull of a set of points.
- * The implementation is based on the Andrew's scan algorithm
- * note: in the Bezier clipping for collinear normals it seems
- * to be more stable wrt the Graham's scan algorithm and in general
- * a bit quikier
- */
-void convex_hull (std::vector<Point> & P)
-{
- size_t n = P.size();
- if (n < 2) return;
- std::sort(P.begin(), P.end(), lex_less());
- if (n < 4) return;
- // upper hull
- size_t u = 2;
- for (size_t i = 2; i < n; ++i)
- {
- while (u > 1 && !is_a_right_turn(P[u-2], P[u-1], P[i]))
- {
- --u;
- }
- std::swap(P[u], P[i]);
- ++u;
- }
- std::sort(P.begin() + u, P.end(), lex_greater());
- std::rotate(P.begin(), P.begin() + 1, P.end());
- // lower hull
- size_t l = u;
- size_t k = u - 1;
- for (size_t i = l; i < n; ++i)
- {
- while (l > k && !is_a_right_turn(P[l-2], P[l-1], P[i]))
- {
- --l;
- }
- std::swap(P[l], P[i]);
- ++l;
- }
- P.resize(l);
-}
-
-
-////////////////////////////////////////////////////////////////////////////////
-// numerical routines
-
-/*
- * Compute the binomial coefficient (n, k)
- */
-inline
-double binomial(unsigned int n, unsigned int k)
-{
- return choose<double>(n, k);
-}
-
-/*
- * Compute the determinant of the 2x2 matrix with column the point P1, P2
- */
-inline
-double det(Point const& P1, Point const& P2)
-{
- return P1[X]*P2[Y] - P1[Y]*P2[X];
-}
-
-/*
- * Solve the linear system [P1,P2] * P = Q
- * in case there isn't exactly one solution the routine returns false
- */
-inline
-bool solve(Point & P, Point const& P1, Point const& P2, Point const& Q)
-{
- double d = det(P1, P2);
- if (d == 0) return false;
- d = 1 / d;
- P[X] = det(Q, P2) * d;
- P[Y] = det(P1, Q) * d;
- return true;
-}
-
-////////////////////////////////////////////////////////////////////////////////
-// interval routines
-
-/*
- * Map the sub-interval I in [0,1] into the interval J and assign it to J
- */
-inline
-void map_to(Interval & J, Interval const& I)
-{
- double length = J.extent();
- J[1] = I.max() * length + J[0];
- J[0] = I.min() * length + J[0];
-}
-
-/*
- * The interval [1,0] is used to represent the empty interval, this routine
- * is just an helper function for creating such an interval
- */
-inline
-Interval make_empty_interval()
-{
- Interval I(0);
- I[0] = 1;
- return I;
-}
-
-
-////////////////////////////////////////////////////////////////////////////////
-// bezier curve routines
-
-/*
- * Return true if all the Bezier curve control points are near,
- * false otherwise
- */
-inline
-bool is_constant(std::vector<Point> const& A, double precision = EPSILON)
-{
- for (unsigned int i = 1; i < A.size(); ++i)
- {
- if(!are_near(A[i], A[0], precision))
- return false;
- }
- return true;
-}
-
-/*
- * Compute the hodograph of the bezier curve B and return it in D
- */
-inline
-void derivative(std::vector<Point> & D, std::vector<Point> const& B)
-{
- D.clear();
- size_t sz = B.size();
- if (sz == 0) return;
- if (sz == 1)
- {
- D.resize(1, Point(0,0));
- return;
- }
- size_t n = sz-1;
- D.reserve(n);
- for (size_t i = 0; i < n; ++i)
- {
- D.push_back(n*(B[i+1] - B[i]));
- }
-}
-
-/*
- * Compute the hodograph of the Bezier curve B rotated of 90 degree
- * and return it in D; we have N(t) orthogonal to B(t) for any t
- */
-inline
-void normal(std::vector<Point> & N, std::vector<Point> const& B)
-{
- derivative(N,B);
- for (size_t i = 0; i < N.size(); ++i)
- {
- N[i] = rot90(N[i]);
- }
-}
-
-/*
- * Compute the portion of the Bezier curve "B" wrt the interval [0,t]
- */
-inline
-void left_portion(Coord t, std::vector<Point> & B)
-{
- size_t n = B.size();
- for (size_t i = 1; i < n; ++i)
- {
- for (size_t j = n-1; j > i-1 ; --j)
- {
- B[j] = lerp(t, B[j-1], B[j]);
- }
- }
-}
-
-/*
- * Compute the portion of the Bezier curve "B" wrt the interval [t,1]
- */
-inline
-void right_portion(Coord t, std::vector<Point> & B)
-{
- size_t n = B.size();
- for (size_t i = 1; i < n; ++i)
- {
- for (size_t j = 0; j < n-i; ++j)
- {
- B[j] = lerp(t, B[j], B[j+1]);
- }
- }
-}
-
-/*
- * Compute the portion of the Bezier curve "B" wrt the interval "I"
- */
-inline
-void portion (std::vector<Point> & B , Interval const& I)
-{
- if (I.min() == 0)
- {
- if (I.max() == 1) return;
- left_portion(I.max(), B);
- return;
- }
- right_portion(I.min(), B);
- if (I.max() == 1) return;
- double t = I.extent() / (1 - I.min());
- left_portion(t, B);
-}
-
-
-////////////////////////////////////////////////////////////////////////////////
-// tags
-
-struct intersection_point_tag;
-struct collinear_normal_tag;
-template <typename Tag>
-void clip(Interval & dom,
- std::vector<Point> const& A,
- std::vector<Point> const& B);
-template <typename Tag>
-void iterate(std::vector<Interval>& domsA,
- std::vector<Interval>& domsB,
- std::vector<Point> const& A,
- std::vector<Point> const& B,
- Interval const& domA,
- Interval const& domB,
- double precision );
-
-
-////////////////////////////////////////////////////////////////////////////////
-// intersection
-
-/*
- * Make up an orientation line using the control points c[i] and c[j]
- * the line is returned in the output parameter "l" in the form of a 3 element
- * vector : l[0] * x + l[1] * y + l[2] == 0; the line is normalized.
- */
-inline
-void orientation_line (std::vector<double> & l,
- std::vector<Point> const& c,
- size_t i, size_t j)
-{
- l[0] = c[j][Y] - c[i][Y];
- l[1] = c[i][X] - c[j][X];
- l[2] = cross(c[i], c[j]);
- double length = std::sqrt(l[0] * l[0] + l[1] * l[1]);
- assert (length != 0);
- l[0] /= length;
- l[1] /= length;
- l[2] /= length;
-}
-
-/*
- * Pick up an orientation line for the Bezier curve "c" and return it in
- * the output parameter "l"
- */
-inline
-void pick_orientation_line (std::vector<double> & l,
- std::vector<Point> const& c)
-{
- size_t i = c.size();
- while (--i > 0 && are_near(c[0], c[i]))
- {}
- if (i == 0)
- {
- // this should never happen because when a new curve portion is created
- // we check that it is not constant;
- // however this requires that the precision used in the is_constant
- // routine has to be the same used here in the are_near test
- assert(i != 0);
- }
- orientation_line(l, c, 0, i);
- //std::cerr << "i = " << i << std::endl;
-}
-
-/*
- * Make up an orientation line for constant bezier curve;
- * the orientation line is made up orthogonal to the other curve base line;
- * the line is returned in the output parameter "l" in the form of a 3 element
- * vector : l[0] * x + l[1] * y + l[2] == 0; the line is normalized.
- */
-inline
-void orthogonal_orientation_line (std::vector<double> & l,
- std::vector<Point> const& c,
- Point const& p)
-{
- if (is_constant(c))
- {
- // this should never happen
- assert(!is_constant(c));
- }
- std::vector<Point> ol(2);
- ol[0] = p;
- ol[1] = (c.back() - c.front()).cw() + p;
- orientation_line(l, ol, 0, 1);
-}
-
-/*
- * Compute the signed distance of the point "P" from the normalized line l
- */
-inline
-double distance (Point const& P, std::vector<double> const& l)
-{
- return l[X] * P[X] + l[Y] * P[Y] + l[2];
-}
-
-/*
- * Compute the min and max distance of the control points of the Bezier
- * curve "c" from the normalized orientation line "l".
- * This bounds are returned through the output Interval parameter"bound".
- */
-inline
-void fat_line_bounds (Interval& bound,
- std::vector<Point> const& c,
- std::vector<double> const& l)
-{
- bound[0] = 0;
- bound[1] = 0;
- double d;
- for (size_t i = 0; i < c.size(); ++i)
- {
- d = distance(c[i], l);
- if (bound[0] > d) bound[0] = d;
- if (bound[1] < d) bound[1] = d;
- }
-}
-
-/*
- * return the x component of the intersection point between the line
- * passing through points p1, p2 and the line Y = "y"
- */
-inline
-double intersect (Point const& p1, Point const& p2, double y)
-{
- // we are sure that p2[Y] != p1[Y] because this routine is called
- // only when the lower or the upper bound is crossed
- double dy = (p2[Y] - p1[Y]);
- double s = (y - p1[Y]) / dy;
- return (p2[X]-p1[X])*s + p1[X];
-}
-
-/*
- * Clip the Bezier curve "B" wrt the fat line defined by the orientation
- * line "l" and the interval range "bound", the new parameter interval for
- * the clipped curve is returned through the output parameter "dom"
- */
-void clip_interval (Interval& dom,
- std::vector<Point> const& B,
- std::vector<double> const& l,
- Interval const& bound)
-{
- double n = B.size() - 1; // number of sub-intervals
- std::vector<Point> D; // distance curve control points
- D.reserve (B.size());
- double d;
- for (size_t i = 0; i < B.size(); ++i)
- {
- d = distance (B[i], l);
- D.push_back (Point(i/n, d));
- }
- //print(D);
-
- convex_hull(D);
- std::vector<Point> & p = D;
- //print(p);
-
- bool plower, phigher;
- bool clower, chigher;
- double t, tmin = 1, tmax = 0;
-// std::cerr << "bound : " << bound << std::endl;
-
- plower = (p[0][Y] < bound.min());
- phigher = (p[0][Y] > bound.max());
- if (!(plower || phigher)) // inside the fat line
- {
- if (tmin > p[0][X]) tmin = p[0][X];
- if (tmax < p[0][X]) tmax = p[0][X];
-// std::cerr << "0 : inside " << p[0]
-// << " : tmin = " << tmin << ", tmax = " << tmax << std::endl;
- }
-
- for (size_t i = 1; i < p.size(); ++i)
- {
- clower = (p[i][Y] < bound.min());
- chigher = (p[i][Y] > bound.max());
- if (!(clower || chigher)) // inside the fat line
- {
- if (tmin > p[i][X]) tmin = p[i][X];
- if (tmax < p[i][X]) tmax = p[i][X];
-// std::cerr << i << " : inside " << p[i]
-// << " : tmin = " << tmin << ", tmax = " << tmax
-// << std::endl;
- }
- if (clower != plower) // cross the lower bound
- {
- t = intersect(p[i-1], p[i], bound.min());
- if (tmin > t) tmin = t;
- if (tmax < t) tmax = t;
- plower = clower;
-// std::cerr << i << " : lower " << p[i]
-// << " : tmin = " << tmin << ", tmax = " << tmax
-// << std::endl;
- }
- if (chigher != phigher) // cross the upper bound
- {
- t = intersect(p[i-1], p[i], bound.max());
- if (tmin > t) tmin = t;
- if (tmax < t) tmax = t;
- phigher = chigher;
-// std::cerr << i << " : higher " << p[i]
-// << " : tmin = " << tmin << ", tmax = " << tmax
-// << std::endl;
- }
- }
-
- // we have to test the closing segment for intersection
- size_t last = p.size() - 1;
- clower = (p[0][Y] < bound.min());
- chigher = (p[0][Y] > bound.max());
- if (clower != plower) // cross the lower bound
- {
- t = intersect(p[last], p[0], bound.min());
- if (tmin > t) tmin = t;
- if (tmax < t) tmax = t;
-// std::cerr << "0 : lower " << p[0]
-// << " : tmin = " << tmin << ", tmax = " << tmax << std::endl;
- }
- if (chigher != phigher) // cross the upper bound
- {
- t = intersect(p[last], p[0], bound.max());
- if (tmin > t) tmin = t;
- if (tmax < t) tmax = t;
-// std::cerr << "0 : higher " << p[0]
-// << " : tmin = " << tmin << ", tmax = " << tmax << std::endl;
- }
-
- dom[0] = tmin;
- dom[1] = tmax;
-}
-
-/*
- * Clip the Bezier curve "B" wrt the Bezier curve "A" for individuating
- * intersection points the new parameter interval for the clipped curve
- * is returned through the output parameter "dom"
- */
-template <>
-inline
-void clip<intersection_point_tag> (Interval & dom,
- std::vector<Point> const& A,
- std::vector<Point> const& B)
-{
- std::vector<double> bl(3);
- Interval bound;
- if (is_constant(A))
- {
- Point M = middle_point(A.front(), A.back());
- orthogonal_orientation_line(bl, B, M);
- }
- else
- {
- pick_orientation_line(bl, A);
- }
- fat_line_bounds(bound, A, bl);
- clip_interval(dom, B, bl, bound);
-}
-
-
-///////////////////////////////////////////////////////////////////////////////
-// collinear normal
-
-/*
- * Compute a closed focus for the Bezier curve B and return it in F
- * A focus is any curve through which all lines perpendicular to B(t) pass.
- */
-inline
-void make_focus (std::vector<Point> & F, std::vector<Point> const& B)
-{
- assert (B.size() > 2);
- size_t n = B.size() - 1;
- normal(F, B);
- Point c(1, 1);
-#if VERBOSE
- if (!solve(c, F[0], -F[n-1], B[n]-B[0]))
- {
- std::cerr << "make_focus: unable to make up a closed focus" << std::endl;
- }
-#else
- solve(c, F[0], -F[n-1], B[n]-B[0]);
-#endif
-// std::cerr << "c = " << c << std::endl;
-
-
- // B(t) + c(t) * N(t)
- double n_inv = 1 / (double)(n);
- Point c0ni;
- F.push_back(c[1] * F[n-1]);
- F[n] += B[n];
- for (size_t i = n-1; i > 0; --i)
- {
- F[i] *= -c[0];
- c0ni = F[i];
- F[i] += (c[1] * F[i-1]);
- F[i] *= (i * n_inv);
- F[i] -= c0ni;
- F[i] += B[i];
- }
- F[0] *= c[0];
- F[0] += B[0];
-}
-
-/*
- * Compute the projection on the plane (t, d) of the control points
- * (t, u, D(t,u)) where D(t,u) = <(B(t) - F(u)), B'(t)> with 0 <= t, u <= 1
- * B is a Bezier curve and F is a focus of another Bezier curve.
- * See Sederberg, Nishita, 1990 - Curve intersection using Bezier clipping.
- */
-void distance_control_points (std::vector<Point> & D,
- std::vector<Point> const& B,
- std::vector<Point> const& F)
-{
- assert (B.size() > 1);
- assert (F.size() > 0);
- const size_t n = B.size() - 1;
- const size_t m = F.size() - 1;
- const size_t r = 2 * n - 1;
- const double r_inv = 1 / (double)(r);
- D.clear();
- D.reserve (B.size() * F.size());
-
- std::vector<Point> dB;
- dB.reserve(n);
- for (size_t k = 0; k < n; ++k)
- {
- dB.push_back (B[k+1] - B[k]);
- }
- NL::Matrix dBB(n,B.size());
- for (size_t i = 0; i < n; ++i)
- for (size_t j = 0; j < B.size(); ++j)
- dBB(i,j) = dot (dB[i], B[j]);
- NL::Matrix dBF(n, F.size());
- for (size_t i = 0; i < n; ++i)
- for (size_t j = 0; j < F.size(); ++j)
- dBF(i,j) = dot (dB[i], F[j]);
-
- size_t k0, kn, l;
- double bc, bri;
- Point dij;
- std::vector<double> d(F.size());
- for (size_t i = 0; i <= r; ++i)
- {
- for (size_t j = 0; j <= m; ++j)
- {
- d[j] = 0;
- }
- k0 = std::max(i, n) - n;
- kn = std::min(i, n-1);
- bri = n / binomial(r,i);
- for (size_t k = k0; k <= kn; ++k)
- {
- //if (k > i || (i-k) > n) continue;
- l = i - k;
-#if CHECK
- assert (l <= n);
-#endif
- bc = bri * binomial(n,l) * binomial(n-1, k);
- for (size_t j = 0; j <= m; ++j)
- {
- //d[j] += bc * dot(dB[k], B[l] - F[j]);
- d[j] += bc * (dBB(k,l) - dBF(k,j));
- }
- }
- double dmin, dmax;
- dmin = dmax = d[m];
- for (size_t j = 0; j < m; ++j)
- {
- if (dmin > d[j]) dmin = d[j];
- if (dmax < d[j]) dmax = d[j];
- }
- dij[0] = i * r_inv;
- dij[1] = dmin;
- D.push_back (dij);
- dij[1] = dmax;
- D.push_back (dij);
- }
-}
-
-/*
- * Clip the Bezier curve "B" wrt the focus "F"; the new parameter interval for
- * the clipped curve is returned through the output parameter "dom"
- */
-void clip_interval (Interval& dom,
- std::vector<Point> const& B,
- std::vector<Point> const& F)
-{
- std::vector<Point> D; // distance curve control points
- distance_control_points(D, B, F);
- //print(D, "D");
-// ConvexHull chD(D);
-// std::vector<Point>& p = chD.boundary; // convex hull vertices
-
- convex_hull(D);
- std::vector<Point> & p = D;
- //print(p, "CH(D)");
-
- bool plower, clower;
- double t, tmin = 1, tmax = 0;
-
- plower = (p[0][Y] < 0);
- if (p[0][Y] == 0) // on the x axis
- {
- if (tmin > p[0][X]) tmin = p[0][X];
- if (tmax < p[0][X]) tmax = p[0][X];
-// std::cerr << "0 : on x axis " << p[0]
-// << " : tmin = " << tmin << ", tmax = " << tmax << std::endl;
- }
-
- for (size_t i = 1; i < p.size(); ++i)
- {
- clower = (p[i][Y] < 0);
- if (p[i][Y] == 0) // on x axis
- {
- if (tmin > p[i][X]) tmin = p[i][X];
- if (tmax < p[i][X]) tmax = p[i][X];
-// std::cerr << i << " : on x axis " << p[i]
-// << " : tmin = " << tmin << ", tmax = " << tmax
-// << std::endl;
- }
- else if (clower != plower) // cross the x axis
- {
- t = intersect(p[i-1], p[i], 0);
- if (tmin > t) tmin = t;
- if (tmax < t) tmax = t;
- plower = clower;
-// std::cerr << i << " : lower " << p[i]
-// << " : tmin = " << tmin << ", tmax = " << tmax
-// << std::endl;
- }
- }
-
- // we have to test the closing segment for intersection
- size_t last = p.size() - 1;
- clower = (p[0][Y] < 0);
- if (clower != plower) // cross the x axis
- {
- t = intersect(p[last], p[0], 0);
- if (tmin > t) tmin = t;
- if (tmax < t) tmax = t;
-// std::cerr << "0 : lower " << p[0]
-// << " : tmin = " << tmin << ", tmax = " << tmax << std::endl;
- }
- dom[0] = tmin;
- dom[1] = tmax;
-}
-
-/*
- * Clip the Bezier curve "B" wrt the Bezier curve "A" for individuating
- * points which have collinear normals; the new parameter interval
- * for the clipped curve is returned through the output parameter "dom"
- */
-template <>
-inline
-void clip<collinear_normal_tag> (Interval & dom,
- std::vector<Point> const& A,
- std::vector<Point> const& B)
-{
- std::vector<Point> F;
- make_focus(F, A);
- clip_interval(dom, B, F);
-}
-
-
-
-const double MAX_PRECISION = 1e-8;
-const double MIN_CLIPPED_SIZE_THRESHOLD = 0.8;
-const Interval UNIT_INTERVAL(0,1);
-const Interval EMPTY_INTERVAL = make_empty_interval();
-const Interval H1_INTERVAL(0, 0.5);
-const Interval H2_INTERVAL(0.5 + MAX_PRECISION, 1.0);
-
-/*
- * iterate
- *
- * input:
- * A, B: control point sets of two bezier curves
- * domA, domB: real parameter intervals of the two curves
- * precision: required computational precision of the returned parameter ranges
- * output:
- * domsA, domsB: sets of parameter intervals
- *
- * The parameter intervals are computed by using a Bezier clipping algorithm,
- * in case the clipping doesn't shrink the initial interval more than 20%,
- * a subdivision step is performed.
- * If during the computation both curves collapse to a single point
- * the routine exits indipendently by the precision reached in the computation
- * of the curve intervals.
- */
-template <>
-void iterate<intersection_point_tag> (std::vector<Interval>& domsA,
- std::vector<Interval>& domsB,
- std::vector<Point> const& A,
- std::vector<Point> const& B,
- Interval const& domA,
- Interval const& domB,
- double precision )
-{
- // in order to limit recursion
- static size_t counter = 0;
- if (domA.extent() == 1 && domB.extent() == 1) counter = 0;
- if (++counter > 100) return;
-#if VERBOSE
- std::cerr << std::fixed << std::setprecision(16);
- std::cerr << ">> curve subdision performed <<" << std::endl;
- std::cerr << "dom(A) : " << domA << std::endl;
- std::cerr << "dom(B) : " << domB << std::endl;
-// std::cerr << "angle(A) : " << angle(A) << std::endl;
-// std::cerr << "angle(B) : " << angle(B) << std::endl;
-#endif
-
- if (precision < MAX_PRECISION)
- precision = MAX_PRECISION;
-
- std::vector<Point> pA = A;
- std::vector<Point> pB = B;
- std::vector<Point>* C1 = &pA;
- std::vector<Point>* C2 = &pB;
-
- Interval dompA = domA;
- Interval dompB = domB;
- Interval* dom1 = &dompA;
- Interval* dom2 = &dompB;
-
- Interval dom;
-
- if ( is_constant(A) && is_constant(B) ){
- Point M1 = middle_point(C1->front(), C1->back());
- Point M2 = middle_point(C2->front(), C2->back());
- if (are_near(M1,M2)){
- domsA.push_back(domA);
- domsB.push_back(domB);
- }
- return;
- }
-
- size_t iter = 0;
- while (++iter < 100
- && (dompA.extent() >= precision || dompB.extent() >= precision))
- {
-#if VERBOSE
- std::cerr << "iter: " << iter << std::endl;
-#endif
- clip<intersection_point_tag>(dom, *C1, *C2);
-
- // [1,0] is utilized to represent an empty interval
- if (dom == EMPTY_INTERVAL)
- {
-#if VERBOSE
- std::cerr << "dom: empty" << std::endl;
-#endif
- return;
- }
-#if VERBOSE
- std::cerr << "dom : " << dom << std::endl;
-#endif
- // all other cases where dom[0] > dom[1] are invalid
- if (dom.min() > dom.max())
- {
- assert(dom.min() < dom.max());
- }
-
- map_to(*dom2, dom);
-
- portion(*C2, dom);
- if (is_constant(*C2) && is_constant(*C1))
- {
- Point M1 = middle_point(C1->front(), C1->back());
- Point M2 = middle_point(C2->front(), C2->back());
-#if VERBOSE
- std::cerr << "both curves are constant: \n"
- << "M1: " << M1 << "\n"
- << "M2: " << M2 << std::endl;
- print(*C2, "C2");
- print(*C1, "C1");
-#endif
- if (are_near(M1,M2))
- break; // append the new interval
- else
- return; // exit without appending any new interval
- }
-
-
- // if we have clipped less than 20% than we need to subdive the curve
- // with the largest domain into two sub-curves
- if ( dom.extent() > MIN_CLIPPED_SIZE_THRESHOLD)
- {
-#if VERBOSE
- std::cerr << "clipped less than 20% : " << dom.extent() << std::endl;
- std::cerr << "angle(pA) : " << angle(pA) << std::endl;
- std::cerr << "angle(pB) : " << angle(pB) << std::endl;
-#endif
- std::vector<Point> pC1, pC2;
- Interval dompC1, dompC2;
- if (dompA.extent() > dompB.extent())
- {
- pC1 = pC2 = pA;
- portion(pC1, H1_INTERVAL);
- portion(pC2, H2_INTERVAL);
- dompC1 = dompC2 = dompA;
- map_to(dompC1, H1_INTERVAL);
- map_to(dompC2, H2_INTERVAL);
- iterate<intersection_point_tag>(domsA, domsB, pC1, pB,
- dompC1, dompB, precision);
- iterate<intersection_point_tag>(domsA, domsB, pC2, pB,
- dompC2, dompB, precision);
- }
- else
- {
- pC1 = pC2 = pB;
- portion(pC1, H1_INTERVAL);
- portion(pC2, H2_INTERVAL);
- dompC1 = dompC2 = dompB;
- map_to(dompC1, H1_INTERVAL);
- map_to(dompC2, H2_INTERVAL);
- iterate<intersection_point_tag>(domsB, domsA, pC1, pA,
- dompC1, dompA, precision);
- iterate<intersection_point_tag>(domsB, domsA, pC2, pA,
- dompC2, dompA, precision);
- }
- return;
- }
-
- std::swap(C1, C2);
- std::swap(dom1, dom2);
-#if VERBOSE
- std::cerr << "dom(pA) : " << dompA << std::endl;
- std::cerr << "dom(pB) : " << dompB << std::endl;
-#endif
- }
- domsA.push_back(dompA);
- domsB.push_back(dompB);
-}
-
-
-/*
- * iterate
- *
- * input:
- * A, B: control point sets of two bezier curves
- * domA, domB: real parameter intervals of the two curves
- * precision: required computational precision of the returned parameter ranges
- * output:
- * domsA, domsB: sets of parameter intervals
- *
- * The parameter intervals are computed by using a Bezier clipping algorithm,
- * in case the clipping doesn't shrink the initial interval more than 20%,
- * a subdivision step is performed.
- * If during the computation one of the two curve interval length becomes less
- * than MAX_PRECISION the routine exits indipendently by the precision reached
- * in the computation of the other curve interval.
- */
-template <>
-void iterate<collinear_normal_tag> (std::vector<Interval>& domsA,
- std::vector<Interval>& domsB,
- std::vector<Point> const& A,
- std::vector<Point> const& B,
- Interval const& domA,
- Interval const& domB,
- double precision)
-{
- // in order to limit recursion
- static size_t counter = 0;
- if (domA.extent() == 1 && domB.extent() == 1) counter = 0;
- if (++counter > 100) return;
-#if VERBOSE
- std::cerr << std::fixed << std::setprecision(16);
- std::cerr << ">> curve subdision performed <<" << std::endl;
- std::cerr << "dom(A) : " << domA << std::endl;
- std::cerr << "dom(B) : " << domB << std::endl;
-// std::cerr << "angle(A) : " << angle(A) << std::endl;
-// std::cerr << "angle(B) : " << angle(B) << std::endl;
-#endif
-
- if (precision < MAX_PRECISION)
- precision = MAX_PRECISION;
-
- std::vector<Point> pA = A;
- std::vector<Point> pB = B;
- std::vector<Point>* C1 = &pA;
- std::vector<Point>* C2 = &pB;
-
- Interval dompA = domA;
- Interval dompB = domB;
- Interval* dom1 = &dompA;
- Interval* dom2 = &dompB;
-
- Interval dom;
-
- size_t iter = 0;
- while (++iter < 100
- && (dompA.extent() >= precision || dompB.extent() >= precision))
- {
-#if VERBOSE
- std::cerr << "iter: " << iter << std::endl;
-#endif
- clip<collinear_normal_tag>(dom, *C1, *C2);
-
- // [1,0] is utilized to represent an empty interval
- if (dom == EMPTY_INTERVAL)
- {
-#if VERBOSE
- std::cerr << "dom: empty" << std::endl;
-#endif
- return;
- }
-#if VERBOSE
- std::cerr << "dom : " << dom << std::endl;
-#endif
- // all other cases where dom[0] > dom[1] are invalid
- if (dom.min() > dom.max())
- {
- assert(dom.min() < dom.max());
- }
-
- map_to(*dom2, dom);
-
- // it's better to stop before losing computational precision
- if (iter > 1 && (dom2->extent() <= MAX_PRECISION))
- {
-#if VERBOSE
- std::cerr << "beyond max precision limit" << std::endl;
-#endif
- break;
- }
-
- portion(*C2, dom);
- if (iter > 1 && is_constant(*C2))
- {
-#if VERBOSE
- std::cerr << "new curve portion pC1 is constant" << std::endl;
-#endif
- break;
- }
-
-
- // if we have clipped less than 20% than we need to subdive the curve
- // with the largest domain into two sub-curves
- if ( dom.extent() > MIN_CLIPPED_SIZE_THRESHOLD)
- {
-#if VERBOSE
- std::cerr << "clipped less than 20% : " << dom.extent() << std::endl;
- std::cerr << "angle(pA) : " << angle(pA) << std::endl;
- std::cerr << "angle(pB) : " << angle(pB) << std::endl;
-#endif
- std::vector<Point> pC1, pC2;
- Interval dompC1, dompC2;
- if (dompA.extent() > dompB.extent())
- {
- if ((dompA.extent() / 2) < MAX_PRECISION)
- {
- break;
- }
- pC1 = pC2 = pA;
- portion(pC1, H1_INTERVAL);
- if (false && is_constant(pC1))
- {
-#if VERBOSE
- std::cerr << "new curve portion pC1 is constant" << std::endl;
-#endif
- break;
- }
- portion(pC2, H2_INTERVAL);
- if (is_constant(pC2))
- {
-#if VERBOSE
- std::cerr << "new curve portion pC2 is constant" << std::endl;
-#endif
- break;
- }
- dompC1 = dompC2 = dompA;
- map_to(dompC1, H1_INTERVAL);
- map_to(dompC2, H2_INTERVAL);
- iterate<collinear_normal_tag>(domsA, domsB, pC1, pB,
- dompC1, dompB, precision);
- iterate<collinear_normal_tag>(domsA, domsB, pC2, pB,
- dompC2, dompB, precision);
- }
- else
- {
- if ((dompB.extent() / 2) < MAX_PRECISION)
- {
- break;
- }
- pC1 = pC2 = pB;
- portion(pC1, H1_INTERVAL);
- if (is_constant(pC1))
- {
-#if VERBOSE
- std::cerr << "new curve portion pC1 is constant" << std::endl;
-#endif
- break;
- }
- portion(pC2, H2_INTERVAL);
- if (is_constant(pC2))
- {
-#if VERBOSE
- std::cerr << "new curve portion pC2 is constant" << std::endl;
-#endif
- break;
- }
- dompC1 = dompC2 = dompB;
- map_to(dompC1, H1_INTERVAL);
- map_to(dompC2, H2_INTERVAL);
- iterate<collinear_normal_tag>(domsB, domsA, pC1, pA,
- dompC1, dompA, precision);
- iterate<collinear_normal_tag>(domsB, domsA, pC2, pA,
- dompC2, dompA, precision);
- }
- return;
- }
-
- std::swap(C1, C2);
- std::swap(dom1, dom2);
-#if VERBOSE
- std::cerr << "dom(pA) : " << dompA << std::endl;
- std::cerr << "dom(pB) : " << dompB << std::endl;
-#endif
- }
- domsA.push_back(dompA);
- domsB.push_back(dompB);
-}
-
-
-/*
- * get_solutions
- *
- * input: A, B - set of control points of two Bezier curve
- * input: precision - required precision of computation
- * input: clip - the routine used for clipping
- * output: xs - set of pairs of parameter values
- * at which the clipping algorithm converges
- *
- * This routine is based on the Bezier Clipping Algorithm,
- * see: Sederberg - Computer Aided Geometric Design
- */
-template <typename Tag>
-void get_solutions (std::vector< std::pair<double, double> >& xs,
- std::vector<Point> const& A,
- std::vector<Point> const& B,
- double precision)
-{
- std::pair<double, double> ci;
- std::vector<Interval> domsA, domsB;
- iterate<Tag> (domsA, domsB, A, B, UNIT_INTERVAL, UNIT_INTERVAL, precision);
- if (domsA.size() != domsB.size())
- {
- assert (domsA.size() == domsB.size());
- }
- xs.clear();
- xs.reserve(domsA.size());
- for (size_t i = 0; i < domsA.size(); ++i)
- {
-#if VERBOSE
- std::cerr << i << " : domA : " << domsA[i] << std::endl;
- std::cerr << "extent A: " << domsA[i].extent() << " ";
- std::cerr << "precision A: " << get_precision(domsA[i]) << std::endl;
- std::cerr << i << " : domB : " << domsB[i] << std::endl;
- std::cerr << "extent B: " << domsB[i].extent() << " ";
- std::cerr << "precision B: " << get_precision(domsB[i]) << std::endl;
-#endif
- ci.first = domsA[i].middle();
- ci.second = domsB[i].middle();
- xs.push_back(ci);
- }
-}
-
-} /* end namespace bezier_clipping */ } /* end namespace detail */
-
-
-/*
- * find_collinear_normal
- *
- * input: A, B - set of control points of two Bezier curve
- * input: precision - required precision of computation
- * output: xs - set of pairs of parameter values
- * at which there are collinear normals
- *
- * This routine is based on the Bezier Clipping Algorithm,
- * see: Sederberg, Nishita, 1990 - Curve intersection using Bezier clipping
- */
-void find_collinear_normal (std::vector< std::pair<double, double> >& xs,
- std::vector<Point> const& A,
- std::vector<Point> const& B,
- double precision)
-{
- using detail::bezier_clipping::get_solutions;
- using detail::bezier_clipping::collinear_normal_tag;
- get_solutions<collinear_normal_tag>(xs, A, B, precision);
-}
-
-
-/*
- * find_intersections_bezier_clipping
- *
- * input: A, B - set of control points of two Bezier curve
- * input: precision - required precision of computation
- * output: xs - set of pairs of parameter values
- * at which crossing happens
- *
- * This routine is based on the Bezier Clipping Algorithm,
- * see: Sederberg, Nishita, 1990 - Curve intersection using Bezier clipping
- */
-void find_intersections_bezier_clipping (std::vector< std::pair<double, double> >& xs,
- std::vector<Point> const& A,
- std::vector<Point> const& B,
- double precision)
-{
- using detail::bezier_clipping::get_solutions;
- using detail::bezier_clipping::intersection_point_tag;
- get_solutions<intersection_point_tag>(xs, A, B, precision);
-}
-
-} // end namespace Geom
-
-
-
-
-/*
- Local Variables:
- mode:c++
- c-file-style:"stroustrup"
- c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
- indent-tabs-mode:nil
- fill-column:99
- End:
-*/
-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
+/*\r
+ * Implement the Bezier clipping algorithm for finding\r
+ * Bezier curve intersection points and collinear normals\r
+ *\r
+ * Authors:\r
+ * Marco Cecchetti <mrcekets at gmail.com>\r
+ *\r
+ * Copyright 2008 authors\r
+ *\r
+ * This library is free software; you can redistribute it and/or\r
+ * modify it either under the terms of the GNU Lesser General Public\r
+ * License version 2.1 as published by the Free Software Foundation\r
+ * (the "LGPL") or, at your option, under the terms of the Mozilla\r
+ * Public License Version 1.1 (the "MPL"). If you do not alter this\r
+ * notice, a recipient may use your version of this file under either\r
+ * the MPL or the LGPL.\r
+ *\r
+ * You should have received a copy of the LGPL along with this library\r
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software\r
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
+ * You should have received a copy of the MPL along with this library\r
+ * in the file COPYING-MPL-1.1\r
+ *\r
+ * The contents of this file are subject to the Mozilla Public License\r
+ * Version 1.1 (the "License"); you may not use this file except in\r
+ * compliance with the License. You may obtain a copy of the License at\r
+ * http://www.mozilla.org/MPL/\r
+ *\r
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY\r
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for\r
+ * the specific language governing rights and limitations.\r
+ */\r
+\r
+\r
+\r
+\r
+#include <2geom/basic-intersection.h>\r
+#include <2geom/choose.h>\r
+#include <2geom/point.h>\r
+#include <2geom/interval.h>\r
+#include <2geom/bezier.h>\r
+//#include <2geom/convex-cover.h>\r
+#include <2geom/numeric/matrix.h>\r
+\r
+#include <cassert>\r
+#include <vector>\r
+#include <algorithm>\r
+#include <utility>\r
+//#include <iomanip>\r
+\r
+\r
+\r
+\r
+#define VERBOSE 0\r
+#define CHECK 0\r
+\r
+\r
+namespace Geom {\r
+\r
+namespace detail { namespace bezier_clipping {\r
+\r
+////////////////////////////////////////////////////////////////////////////////\r
+// for debugging\r
+//\r
+\r
+inline\r
+void print(std::vector<Point> const& cp, const char* msg = "")\r
+{\r
+ std::cerr << msg << std::endl;\r
+ for (size_t i = 0; i < cp.size(); ++i)\r
+ std::cerr << i << " : " << cp[i] << std::endl;\r
+}\r
+\r
+template< class charT >\r
+inline\r
+std::basic_ostream<charT> &\r
+operator<< (std::basic_ostream<charT> & os, const Interval & I)\r
+{\r
+ os << "[" << I.min() << ", " << I.max() << "]";\r
+ return os;\r
+}\r
+\r
+inline\r
+double angle (std::vector<Point> const& A)\r
+{\r
+ size_t n = A.size() -1;\r
+ double a = std::atan2(A[n][Y] - A[0][Y], A[n][X] - A[0][X]);\r
+ return (180 * a / M_PI);\r
+}\r
+\r
+inline\r
+size_t get_precision(Interval const& I)\r
+{\r
+ double d = I.extent();\r
+ double e = 0.1, p = 10;\r
+ int n = 0;\r
+ while (n < 16 && d < e)\r
+ {\r
+ p *= 10;\r
+ e = 1/p;\r
+ ++n;\r
+ }\r
+ return n;\r
+}\r
+\r
+inline\r
+void range_assertion(int k, int m, int n, const char* msg)\r
+{\r
+ if ( k < m || k > n)\r
+ {\r
+ std::cerr << "range assertion failed: \n"\r
+ << msg << std::endl\r
+ << "value: " << k\r
+ << " range: " << m << ", " << n << std::endl;\r
+ assert (k >= m && k <= n);\r
+ }\r
+}\r
+\r
+\r
+////////////////////////////////////////////////////////////////////////////////\r
+// convex hull\r
+\r
+/*\r
+ * return true in case the oriented polyline p0, p1, p2 is a right turn\r
+ */\r
+inline\r
+bool is_a_right_turn (Point const& p0, Point const& p1, Point const& p2)\r
+{\r
+ if (p1 == p2) return false;\r
+ Point q1 = p1 - p0;\r
+ Point q2 = p2 - p0;\r
+ if (q1 == -q2) return false;\r
+ return (cross (q1, q2) < 0);\r
+}\r
+\r
+/*\r
+ * return true if p < q wrt the lexicographyc order induced by the coordinates\r
+ */\r
+struct lex_less\r
+{\r
+ bool operator() (Point const& p, Point const& q)\r
+ {\r
+ return ((p[X] < q[X]) || (p[X] == q[X] && p[Y] < q[Y]));\r
+ }\r
+};\r
+\r
+/*\r
+ * return true if p > q wrt the lexicographyc order induced by the coordinates\r
+ */\r
+struct lex_greater\r
+{\r
+ bool operator() (Point const& p, Point const& q)\r
+ {\r
+ return ((p[X] > q[X]) || (p[X] == q[X] && p[Y] > q[Y]));\r
+ }\r
+};\r
+\r
+/*\r
+ * Compute the convex hull of a set of points.\r
+ * The implementation is based on the Andrew's scan algorithm\r
+ * note: in the Bezier clipping for collinear normals it seems\r
+ * to be more stable wrt the Graham's scan algorithm and in general\r
+ * a bit quikier\r
+ */\r
+void convex_hull (std::vector<Point> & P)\r
+{\r
+ size_t n = P.size();\r
+ if (n < 2) return;\r
+ std::sort(P.begin(), P.end(), lex_less());\r
+ if (n < 4) return;\r
+ // upper hull\r
+ size_t u = 2;\r
+ for (size_t i = 2; i < n; ++i)\r
+ {\r
+ while (u > 1 && !is_a_right_turn(P[u-2], P[u-1], P[i]))\r
+ {\r
+ --u;\r
+ }\r
+ std::swap(P[u], P[i]);\r
+ ++u;\r
+ }\r
+ std::sort(P.begin() + u, P.end(), lex_greater());\r
+ std::rotate(P.begin(), P.begin() + 1, P.end());\r
+ // lower hull\r
+ size_t l = u;\r
+ size_t k = u - 1;\r
+ for (size_t i = l; i < n; ++i)\r
+ {\r
+ while (l > k && !is_a_right_turn(P[l-2], P[l-1], P[i]))\r
+ {\r
+ --l;\r
+ }\r
+ std::swap(P[l], P[i]);\r
+ ++l;\r
+ }\r
+ P.resize(l);\r
+}\r
+\r
+\r
+////////////////////////////////////////////////////////////////////////////////\r
+// numerical routines\r
+\r
+/*\r
+ * Compute the binomial coefficient (n, k)\r
+ */\r
+inline\r
+double binomial(unsigned int n, unsigned int k)\r
+{\r
+ return choose<double>(n, k);\r
+}\r
+\r
+/*\r
+ * Compute the determinant of the 2x2 matrix with column the point P1, P2\r
+ */\r
+inline\r
+double det(Point const& P1, Point const& P2)\r
+{\r
+ return P1[X]*P2[Y] - P1[Y]*P2[X];\r
+}\r
+\r
+/*\r
+ * Solve the linear system [P1,P2] * P = Q\r
+ * in case there isn't exactly one solution the routine returns false\r
+ */\r
+inline\r
+bool solve(Point & P, Point const& P1, Point const& P2, Point const& Q)\r
+{\r
+ double d = det(P1, P2);\r
+ if (d == 0) return false;\r
+ d = 1 / d;\r
+ P[X] = det(Q, P2) * d;\r
+ P[Y] = det(P1, Q) * d;\r
+ return true;\r
+}\r
+\r
+////////////////////////////////////////////////////////////////////////////////\r
+// interval routines\r
+\r
+/*\r
+ * Map the sub-interval I in [0,1] into the interval J and assign it to J\r
+ */\r
+inline\r
+void map_to(Interval & J, Interval const& I)\r
+{\r
+ double length = J.extent();\r
+ J[1] = I.max() * length + J[0];\r
+ J[0] = I.min() * length + J[0];\r
+}\r
+\r
+/*\r
+ * The interval [1,0] is used to represent the empty interval, this routine\r
+ * is just an helper function for creating such an interval\r
+ */\r
+inline\r
+Interval make_empty_interval()\r
+{\r
+ Interval I(0);\r
+ I[0] = 1;\r
+ return I;\r
+}\r
+\r
+\r
+////////////////////////////////////////////////////////////////////////////////\r
+// bezier curve routines\r
+\r
+/*\r
+ * Return true if all the Bezier curve control points are near,\r
+ * false otherwise\r
+ */\r
+inline\r
+bool is_constant(std::vector<Point> const& A, double precision = EPSILON)\r
+{\r
+ for (unsigned int i = 1; i < A.size(); ++i)\r
+ {\r
+ if(!are_near(A[i], A[0], precision))\r
+ return false;\r
+ }\r
+ return true;\r
+}\r
+\r
+/*\r
+ * Compute the hodograph of the bezier curve B and return it in D\r
+ */\r
+inline\r
+void derivative(std::vector<Point> & D, std::vector<Point> const& B)\r
+{\r
+ D.clear();\r
+ size_t sz = B.size();\r
+ if (sz == 0) return;\r
+ if (sz == 1)\r
+ {\r
+ D.resize(1, Point(0,0));\r
+ return;\r
+ }\r
+ size_t n = sz-1;\r
+ D.reserve(n);\r
+ for (size_t i = 0; i < n; ++i)\r
+ {\r
+ D.push_back(n*(B[i+1] - B[i]));\r
+ }\r
+}\r
+\r
+/*\r
+ * Compute the hodograph of the Bezier curve B rotated of 90 degree\r
+ * and return it in D; we have N(t) orthogonal to B(t) for any t\r
+ */\r
+inline\r
+void normal(std::vector<Point> & N, std::vector<Point> const& B)\r
+{\r
+ derivative(N,B);\r
+ for (size_t i = 0; i < N.size(); ++i)\r
+ {\r
+ N[i] = rot90(N[i]);\r
+ }\r
+}\r
+\r
+/*\r
+ * Compute the portion of the Bezier curve "B" wrt the interval [0,t]\r
+ */\r
+inline\r
+void left_portion(Coord t, std::vector<Point> & B)\r
+{\r
+ size_t n = B.size();\r
+ for (size_t i = 1; i < n; ++i)\r
+ {\r
+ for (size_t j = n-1; j > i-1 ; --j)\r
+ {\r
+ B[j] = lerp(t, B[j-1], B[j]);\r
+ }\r
+ }\r
+}\r
+\r
+/*\r
+ * Compute the portion of the Bezier curve "B" wrt the interval [t,1]\r
+ */\r
+inline\r
+void right_portion(Coord t, std::vector<Point> & B)\r
+{\r
+ size_t n = B.size();\r
+ for (size_t i = 1; i < n; ++i)\r
+ {\r
+ for (size_t j = 0; j < n-i; ++j)\r
+ {\r
+ B[j] = lerp(t, B[j], B[j+1]);\r
+ }\r
+ }\r
+}\r
+\r
+/*\r
+ * Compute the portion of the Bezier curve "B" wrt the interval "I"\r
+ */\r
+inline\r
+void portion (std::vector<Point> & B , Interval const& I)\r
+{\r
+ if (I.min() == 0)\r
+ {\r
+ if (I.max() == 1) return;\r
+ left_portion(I.max(), B);\r
+ return;\r
+ }\r
+ right_portion(I.min(), B);\r
+ if (I.max() == 1) return;\r
+ double t = I.extent() / (1 - I.min());\r
+ left_portion(t, B);\r
+}\r
+\r
+\r
+////////////////////////////////////////////////////////////////////////////////\r
+// tags\r
+\r
+struct intersection_point_tag;\r
+struct collinear_normal_tag;\r
+template <typename Tag>\r
+void clip(Interval & dom,\r
+ std::vector<Point> const& A,\r
+ std::vector<Point> const& B);\r
+template <typename Tag>\r
+void iterate(std::vector<Interval>& domsA,\r
+ std::vector<Interval>& domsB,\r
+ std::vector<Point> const& A,\r
+ std::vector<Point> const& B,\r
+ Interval const& domA,\r
+ Interval const& domB,\r
+ double precision );\r
+\r
+\r
+////////////////////////////////////////////////////////////////////////////////\r
+// intersection\r
+\r
+/*\r
+ * Make up an orientation line using the control points c[i] and c[j]\r
+ * the line is returned in the output parameter "l" in the form of a 3 element\r
+ * vector : l[0] * x + l[1] * y + l[2] == 0; the line is normalized.\r
+ */\r
+inline\r
+void orientation_line (std::vector<double> & l,\r
+ std::vector<Point> const& c,\r
+ size_t i, size_t j)\r
+{\r
+ l[0] = c[j][Y] - c[i][Y];\r
+ l[1] = c[i][X] - c[j][X];\r
+ l[2] = cross(c[i], c[j]);\r
+ double length = std::sqrt(l[0] * l[0] + l[1] * l[1]);\r
+ assert (length != 0);\r
+ l[0] /= length;\r
+ l[1] /= length;\r
+ l[2] /= length;\r
+}\r
+\r
+/*\r
+ * Pick up an orientation line for the Bezier curve "c" and return it in\r
+ * the output parameter "l"\r
+ */\r
+inline\r
+void pick_orientation_line (std::vector<double> & l,\r
+ std::vector<Point> const& c)\r
+{\r
+ size_t i = c.size();\r
+ while (--i > 0 && are_near(c[0], c[i]))\r
+ {}\r
+ if (i == 0)\r
+ {\r
+ // this should never happen because when a new curve portion is created\r
+ // we check that it is not constant;\r
+ // however this requires that the precision used in the is_constant\r
+ // routine has to be the same used here in the are_near test\r
+ assert(i != 0);\r
+ }\r
+ orientation_line(l, c, 0, i);\r
+ //std::cerr << "i = " << i << std::endl;\r
+}\r
+\r
+/*\r
+ * Make up an orientation line for constant bezier curve;\r
+ * the orientation line is made up orthogonal to the other curve base line;\r
+ * the line is returned in the output parameter "l" in the form of a 3 element\r
+ * vector : l[0] * x + l[1] * y + l[2] == 0; the line is normalized.\r
+ */\r
+inline\r
+void orthogonal_orientation_line (std::vector<double> & l,\r
+ std::vector<Point> const& c,\r
+ Point const& p)\r
+{\r
+ if (is_constant(c))\r
+ {\r
+ // this should never happen\r
+ assert(!is_constant(c));\r
+ }\r
+ std::vector<Point> ol(2);\r
+ ol[0] = p;\r
+ ol[1] = (c.back() - c.front()).cw() + p;\r
+ orientation_line(l, ol, 0, 1);\r
+}\r
+\r
+/*\r
+ * Compute the signed distance of the point "P" from the normalized line l\r
+ */\r
+inline\r
+double distance (Point const& P, std::vector<double> const& l)\r
+{\r
+ return l[X] * P[X] + l[Y] * P[Y] + l[2];\r
+}\r
+\r
+/*\r
+ * Compute the min and max distance of the control points of the Bezier\r
+ * curve "c" from the normalized orientation line "l".\r
+ * This bounds are returned through the output Interval parameter"bound".\r
+ */\r
+inline\r
+void fat_line_bounds (Interval& bound,\r
+ std::vector<Point> const& c,\r
+ std::vector<double> const& l)\r
+{\r
+ bound[0] = 0;\r
+ bound[1] = 0;\r
+ double d;\r
+ for (size_t i = 0; i < c.size(); ++i)\r
+ {\r
+ d = distance(c[i], l);\r
+ if (bound[0] > d) bound[0] = d;\r
+ if (bound[1] < d) bound[1] = d;\r
+ }\r
+}\r
+\r
+/*\r
+ * return the x component of the intersection point between the line\r
+ * passing through points p1, p2 and the line Y = "y"\r
+ */\r
+inline\r
+double intersect (Point const& p1, Point const& p2, double y)\r
+{\r
+ // we are sure that p2[Y] != p1[Y] because this routine is called\r
+ // only when the lower or the upper bound is crossed\r
+ double dy = (p2[Y] - p1[Y]);\r
+ double s = (y - p1[Y]) / dy;\r
+ return (p2[X]-p1[X])*s + p1[X];\r
+}\r
+\r
+/*\r
+ * Clip the Bezier curve "B" wrt the fat line defined by the orientation\r
+ * line "l" and the interval range "bound", the new parameter interval for\r
+ * the clipped curve is returned through the output parameter "dom"\r
+ */\r
+void clip_interval (Interval& dom,\r
+ std::vector<Point> const& B,\r
+ std::vector<double> const& l,\r
+ Interval const& bound)\r
+{\r
+ double n = B.size() - 1; // number of sub-intervals\r
+ std::vector<Point> D; // distance curve control points\r
+ D.reserve (B.size());\r
+ double d;\r
+ for (size_t i = 0; i < B.size(); ++i)\r
+ {\r
+ d = distance (B[i], l);\r
+ D.push_back (Point(i/n, d));\r
+ }\r
+ //print(D);\r
+\r
+ convex_hull(D);\r
+ std::vector<Point> & p = D;\r
+ //print(p);\r
+\r
+ bool plower, phigher;\r
+ bool clower, chigher;\r
+ double t, tmin = 1, tmax = 0;\r
+// std::cerr << "bound : " << bound << std::endl;\r
+\r
+ plower = (p[0][Y] < bound.min());\r
+ phigher = (p[0][Y] > bound.max());\r
+ if (!(plower || phigher)) // inside the fat line\r
+ {\r
+ if (tmin > p[0][X]) tmin = p[0][X];\r
+ if (tmax < p[0][X]) tmax = p[0][X];\r
+// std::cerr << "0 : inside " << p[0]\r
+// << " : tmin = " << tmin << ", tmax = " << tmax << std::endl;\r
+ }\r
+\r
+ for (size_t i = 1; i < p.size(); ++i)\r
+ {\r
+ clower = (p[i][Y] < bound.min());\r
+ chigher = (p[i][Y] > bound.max());\r
+ if (!(clower || chigher)) // inside the fat line\r
+ {\r
+ if (tmin > p[i][X]) tmin = p[i][X];\r
+ if (tmax < p[i][X]) tmax = p[i][X];\r
+// std::cerr << i << " : inside " << p[i]\r
+// << " : tmin = " << tmin << ", tmax = " << tmax\r
+// << std::endl;\r
+ }\r
+ if (clower != plower) // cross the lower bound\r
+ {\r
+ t = intersect(p[i-1], p[i], bound.min());\r
+ if (tmin > t) tmin = t;\r
+ if (tmax < t) tmax = t;\r
+ plower = clower;\r
+// std::cerr << i << " : lower " << p[i]\r
+// << " : tmin = " << tmin << ", tmax = " << tmax\r
+// << std::endl;\r
+ }\r
+ if (chigher != phigher) // cross the upper bound\r
+ {\r
+ t = intersect(p[i-1], p[i], bound.max());\r
+ if (tmin > t) tmin = t;\r
+ if (tmax < t) tmax = t;\r
+ phigher = chigher;\r
+// std::cerr << i << " : higher " << p[i]\r
+// << " : tmin = " << tmin << ", tmax = " << tmax\r
+// << std::endl;\r
+ }\r
+ }\r
+\r
+ // we have to test the closing segment for intersection\r
+ size_t last = p.size() - 1;\r
+ clower = (p[0][Y] < bound.min());\r
+ chigher = (p[0][Y] > bound.max());\r
+ if (clower != plower) // cross the lower bound\r
+ {\r
+ t = intersect(p[last], p[0], bound.min());\r
+ if (tmin > t) tmin = t;\r
+ if (tmax < t) tmax = t;\r
+// std::cerr << "0 : lower " << p[0]\r
+// << " : tmin = " << tmin << ", tmax = " << tmax << std::endl;\r
+ }\r
+ if (chigher != phigher) // cross the upper bound\r
+ {\r
+ t = intersect(p[last], p[0], bound.max());\r
+ if (tmin > t) tmin = t;\r
+ if (tmax < t) tmax = t;\r
+// std::cerr << "0 : higher " << p[0]\r
+// << " : tmin = " << tmin << ", tmax = " << tmax << std::endl;\r
+ }\r
+\r
+ dom[0] = tmin;\r
+ dom[1] = tmax;\r
+}\r
+\r
+/*\r
+ * Clip the Bezier curve "B" wrt the Bezier curve "A" for individuating\r
+ * intersection points the new parameter interval for the clipped curve\r
+ * is returned through the output parameter "dom"\r
+ */\r
+template <>\r
+inline\r
+void clip<intersection_point_tag> (Interval & dom,\r
+ std::vector<Point> const& A,\r
+ std::vector<Point> const& B)\r
+{\r
+ std::vector<double> bl(3);\r
+ Interval bound;\r
+ if (is_constant(A))\r
+ {\r
+ Point M = middle_point(A.front(), A.back());\r
+ orthogonal_orientation_line(bl, B, M);\r
+ }\r
+ else\r
+ {\r
+ pick_orientation_line(bl, A);\r
+ }\r
+ fat_line_bounds(bound, A, bl);\r
+ clip_interval(dom, B, bl, bound);\r
+}\r
+\r
+\r
+///////////////////////////////////////////////////////////////////////////////\r
+// collinear normal\r
+\r
+/*\r
+ * Compute a closed focus for the Bezier curve B and return it in F\r
+ * A focus is any curve through which all lines perpendicular to B(t) pass.\r
+ */\r
+inline\r
+void make_focus (std::vector<Point> & F, std::vector<Point> const& B)\r
+{\r
+ assert (B.size() > 2);\r
+ size_t n = B.size() - 1;\r
+ normal(F, B);\r
+ Point c(1, 1);\r
+#if VERBOSE\r
+ if (!solve(c, F[0], -F[n-1], B[n]-B[0]))\r
+ {\r
+ std::cerr << "make_focus: unable to make up a closed focus" << std::endl;\r
+ }\r
+#else\r
+ solve(c, F[0], -F[n-1], B[n]-B[0]);\r
+#endif\r
+// std::cerr << "c = " << c << std::endl;\r
+\r
+\r
+ // B(t) + c(t) * N(t)\r
+ double n_inv = 1 / (double)(n);\r
+ Point c0ni;\r
+ F.push_back(c[1] * F[n-1]);\r
+ F[n] += B[n];\r
+ for (size_t i = n-1; i > 0; --i)\r
+ {\r
+ F[i] *= -c[0];\r
+ c0ni = F[i];\r
+ F[i] += (c[1] * F[i-1]);\r
+ F[i] *= (i * n_inv);\r
+ F[i] -= c0ni;\r
+ F[i] += B[i];\r
+ }\r
+ F[0] *= c[0];\r
+ F[0] += B[0];\r
+}\r
+\r
+/*\r
+ * Compute the projection on the plane (t, d) of the control points\r
+ * (t, u, D(t,u)) where D(t,u) = <(B(t) - F(u)), B'(t)> with 0 <= t, u <= 1\r
+ * B is a Bezier curve and F is a focus of another Bezier curve.\r
+ * See Sederberg, Nishita, 1990 - Curve intersection using Bezier clipping.\r
+ */\r
+void distance_control_points (std::vector<Point> & D,\r
+ std::vector<Point> const& B,\r
+ std::vector<Point> const& F)\r
+{\r
+ assert (B.size() > 1);\r
+ assert (F.size() > 0);\r
+ const size_t n = B.size() - 1;\r
+ const size_t m = F.size() - 1;\r
+ const size_t r = 2 * n - 1;\r
+ const double r_inv = 1 / (double)(r);\r
+ D.clear();\r
+ D.reserve (B.size() * F.size());\r
+\r
+ std::vector<Point> dB;\r
+ dB.reserve(n);\r
+ for (size_t k = 0; k < n; ++k)\r
+ {\r
+ dB.push_back (B[k+1] - B[k]);\r
+ }\r
+ NL::Matrix dBB(n,B.size());\r
+ for (size_t i = 0; i < n; ++i)\r
+ for (size_t j = 0; j < B.size(); ++j)\r
+ dBB(i,j) = dot (dB[i], B[j]);\r
+ NL::Matrix dBF(n, F.size());\r
+ for (size_t i = 0; i < n; ++i)\r
+ for (size_t j = 0; j < F.size(); ++j)\r
+ dBF(i,j) = dot (dB[i], F[j]);\r
+\r
+ size_t k0, kn, l;\r
+ double bc, bri;\r
+ Point dij;\r
+ std::vector<double> d(F.size());\r
+ for (size_t i = 0; i <= r; ++i)\r
+ {\r
+ for (size_t j = 0; j <= m; ++j)\r
+ {\r
+ d[j] = 0;\r
+ }\r
+ k0 = std::max(i, n) - n;\r
+ kn = std::min(i, n-1);\r
+ bri = n / binomial(r,i);\r
+ for (size_t k = k0; k <= kn; ++k)\r
+ {\r
+ //if (k > i || (i-k) > n) continue;\r
+ l = i - k;\r
+#if CHECK\r
+ assert (l <= n);\r
+#endif\r
+ bc = bri * binomial(n,l) * binomial(n-1, k);\r
+ for (size_t j = 0; j <= m; ++j)\r
+ {\r
+ //d[j] += bc * dot(dB[k], B[l] - F[j]);\r
+ d[j] += bc * (dBB(k,l) - dBF(k,j));\r
+ }\r
+ }\r
+ double dmin, dmax;\r
+ dmin = dmax = d[m];\r
+ for (size_t j = 0; j < m; ++j)\r
+ {\r
+ if (dmin > d[j]) dmin = d[j];\r
+ if (dmax < d[j]) dmax = d[j];\r
+ }\r
+ dij[0] = i * r_inv;\r
+ dij[1] = dmin;\r
+ D.push_back (dij);\r
+ dij[1] = dmax;\r
+ D.push_back (dij);\r
+ }\r
+}\r
+\r
+/*\r
+ * Clip the Bezier curve "B" wrt the focus "F"; the new parameter interval for\r
+ * the clipped curve is returned through the output parameter "dom"\r
+ */\r
+void clip_interval (Interval& dom,\r
+ std::vector<Point> const& B,\r
+ std::vector<Point> const& F)\r
+{\r
+ std::vector<Point> D; // distance curve control points\r
+ distance_control_points(D, B, F);\r
+ //print(D, "D");\r
+// ConvexHull chD(D);\r
+// std::vector<Point>& p = chD.boundary; // convex hull vertices\r
+\r
+ convex_hull(D);\r
+ std::vector<Point> & p = D;\r
+ //print(p, "CH(D)");\r
+\r
+ bool plower, clower;\r
+ double t, tmin = 1, tmax = 0;\r
+\r
+ plower = (p[0][Y] < 0);\r
+ if (p[0][Y] == 0) // on the x axis\r
+ {\r
+ if (tmin > p[0][X]) tmin = p[0][X];\r
+ if (tmax < p[0][X]) tmax = p[0][X];\r
+// std::cerr << "0 : on x axis " << p[0]\r
+// << " : tmin = " << tmin << ", tmax = " << tmax << std::endl;\r
+ }\r
+\r
+ for (size_t i = 1; i < p.size(); ++i)\r
+ {\r
+ clower = (p[i][Y] < 0);\r
+ if (p[i][Y] == 0) // on x axis\r
+ {\r
+ if (tmin > p[i][X]) tmin = p[i][X];\r
+ if (tmax < p[i][X]) tmax = p[i][X];\r
+// std::cerr << i << " : on x axis " << p[i]\r
+// << " : tmin = " << tmin << ", tmax = " << tmax\r
+// << std::endl;\r
+ }\r
+ else if (clower != plower) // cross the x axis\r
+ {\r
+ t = intersect(p[i-1], p[i], 0);\r
+ if (tmin > t) tmin = t;\r
+ if (tmax < t) tmax = t;\r
+ plower = clower;\r
+// std::cerr << i << " : lower " << p[i]\r
+// << " : tmin = " << tmin << ", tmax = " << tmax\r
+// << std::endl;\r
+ }\r
+ }\r
+\r
+ // we have to test the closing segment for intersection\r
+ size_t last = p.size() - 1;\r
+ clower = (p[0][Y] < 0);\r
+ if (clower != plower) // cross the x axis\r
+ {\r
+ t = intersect(p[last], p[0], 0);\r
+ if (tmin > t) tmin = t;\r
+ if (tmax < t) tmax = t;\r
+// std::cerr << "0 : lower " << p[0]\r
+// << " : tmin = " << tmin << ", tmax = " << tmax << std::endl;\r
+ }\r
+ dom[0] = tmin;\r
+ dom[1] = tmax;\r
+}\r
+\r
+/*\r
+ * Clip the Bezier curve "B" wrt the Bezier curve "A" for individuating\r
+ * points which have collinear normals; the new parameter interval\r
+ * for the clipped curve is returned through the output parameter "dom"\r
+ */\r
+template <>\r
+inline\r
+void clip<collinear_normal_tag> (Interval & dom,\r
+ std::vector<Point> const& A,\r
+ std::vector<Point> const& B)\r
+{\r
+ std::vector<Point> F;\r
+ make_focus(F, A);\r
+ clip_interval(dom, B, F);\r
+}\r
+\r
+\r
+\r
+const double MAX_PRECISION = 1e-8;\r
+const double MIN_CLIPPED_SIZE_THRESHOLD = 0.8;\r
+const Interval UNIT_INTERVAL(0,1);\r
+const Interval EMPTY_INTERVAL = make_empty_interval();\r
+const Interval H1_INTERVAL(0, 0.5);\r
+const Interval H2_INTERVAL(0.5 + MAX_PRECISION, 1.0);\r
+\r
+/*\r
+ * iterate\r
+ *\r
+ * input:\r
+ * A, B: control point sets of two bezier curves\r
+ * domA, domB: real parameter intervals of the two curves\r
+ * precision: required computational precision of the returned parameter ranges\r
+ * output:\r
+ * domsA, domsB: sets of parameter intervals\r
+ *\r
+ * The parameter intervals are computed by using a Bezier clipping algorithm,\r
+ * in case the clipping doesn't shrink the initial interval more than 20%,\r
+ * a subdivision step is performed.\r
+ * If during the computation both curves collapse to a single point\r
+ * the routine exits indipendently by the precision reached in the computation\r
+ * of the curve intervals.\r
+ */\r
+template <>\r
+void iterate<intersection_point_tag> (std::vector<Interval>& domsA,\r
+ std::vector<Interval>& domsB,\r
+ std::vector<Point> const& A,\r
+ std::vector<Point> const& B,\r
+ Interval const& domA,\r
+ Interval const& domB,\r
+ double precision )\r
+{\r
+ // in order to limit recursion\r
+ static size_t counter = 0;\r
+ if (domA.extent() == 1 && domB.extent() == 1) counter = 0;\r
+ if (++counter > 100) return;\r
+#if VERBOSE\r
+ std::cerr << std::fixed << std::setprecision(16);\r
+ std::cerr << ">> curve subdision performed <<" << std::endl;\r
+ std::cerr << "dom(A) : " << domA << std::endl;\r
+ std::cerr << "dom(B) : " << domB << std::endl;\r
+// std::cerr << "angle(A) : " << angle(A) << std::endl;\r
+// std::cerr << "angle(B) : " << angle(B) << std::endl;\r
+#endif\r
+\r
+ if (precision < MAX_PRECISION)\r
+ precision = MAX_PRECISION;\r
+\r
+ std::vector<Point> pA = A;\r
+ std::vector<Point> pB = B;\r
+ std::vector<Point>* C1 = &pA;\r
+ std::vector<Point>* C2 = &pB;\r
+\r
+ Interval dompA = domA;\r
+ Interval dompB = domB;\r
+ Interval* dom1 = &dompA;\r
+ Interval* dom2 = &dompB;\r
+\r
+ Interval dom;\r
+\r
+ if ( is_constant(A) && is_constant(B) ){\r
+ Point M1 = middle_point(C1->front(), C1->back());\r
+ Point M2 = middle_point(C2->front(), C2->back());\r
+ if (are_near(M1,M2)){\r
+ domsA.push_back(domA);\r
+ domsB.push_back(domB);\r
+ }\r
+ return;\r
+ }\r
+\r
+ size_t iter = 0;\r
+ while (++iter < 100\r
+ && (dompA.extent() >= precision || dompB.extent() >= precision))\r
+ {\r
+#if VERBOSE\r
+ std::cerr << "iter: " << iter << std::endl;\r
+#endif\r
+ clip<intersection_point_tag>(dom, *C1, *C2);\r
+\r
+ // [1,0] is utilized to represent an empty interval\r
+ if (dom == EMPTY_INTERVAL)\r
+ {\r
+#if VERBOSE\r
+ std::cerr << "dom: empty" << std::endl;\r
+#endif\r
+ return;\r
+ }\r
+#if VERBOSE\r
+ std::cerr << "dom : " << dom << std::endl;\r
+#endif\r
+ // all other cases where dom[0] > dom[1] are invalid\r
+ if (dom.min() > dom.max())\r
+ {\r
+ assert(dom.min() < dom.max());\r
+ }\r
+\r
+ map_to(*dom2, dom);\r
+\r
+ portion(*C2, dom);\r
+ if (is_constant(*C2) && is_constant(*C1))\r
+ {\r
+ Point M1 = middle_point(C1->front(), C1->back());\r
+ Point M2 = middle_point(C2->front(), C2->back());\r
+#if VERBOSE\r
+ std::cerr << "both curves are constant: \n"\r
+ << "M1: " << M1 << "\n"\r
+ << "M2: " << M2 << std::endl;\r
+ print(*C2, "C2");\r
+ print(*C1, "C1");\r
+#endif\r
+ if (are_near(M1,M2))\r
+ break; // append the new interval\r
+ else\r
+ return; // exit without appending any new interval\r
+ }\r
+\r
+\r
+ // if we have clipped less than 20% than we need to subdive the curve\r
+ // with the largest domain into two sub-curves\r
+ if ( dom.extent() > MIN_CLIPPED_SIZE_THRESHOLD)\r
+ {\r
+#if VERBOSE\r
+ std::cerr << "clipped less than 20% : " << dom.extent() << std::endl;\r
+ std::cerr << "angle(pA) : " << angle(pA) << std::endl;\r
+ std::cerr << "angle(pB) : " << angle(pB) << std::endl;\r
+#endif\r
+ std::vector<Point> pC1, pC2;\r
+ Interval dompC1, dompC2;\r
+ if (dompA.extent() > dompB.extent())\r
+ {\r
+ pC1 = pC2 = pA;\r
+ portion(pC1, H1_INTERVAL);\r
+ portion(pC2, H2_INTERVAL);\r
+ dompC1 = dompC2 = dompA;\r
+ map_to(dompC1, H1_INTERVAL);\r
+ map_to(dompC2, H2_INTERVAL);\r
+ iterate<intersection_point_tag>(domsA, domsB, pC1, pB,\r
+ dompC1, dompB, precision);\r
+ iterate<intersection_point_tag>(domsA, domsB, pC2, pB,\r
+ dompC2, dompB, precision);\r
+ }\r
+ else\r
+ {\r
+ pC1 = pC2 = pB;\r
+ portion(pC1, H1_INTERVAL);\r
+ portion(pC2, H2_INTERVAL);\r
+ dompC1 = dompC2 = dompB;\r
+ map_to(dompC1, H1_INTERVAL);\r
+ map_to(dompC2, H2_INTERVAL);\r
+ iterate<intersection_point_tag>(domsB, domsA, pC1, pA,\r
+ dompC1, dompA, precision);\r
+ iterate<intersection_point_tag>(domsB, domsA, pC2, pA,\r
+ dompC2, dompA, precision);\r
+ }\r
+ return;\r
+ }\r
+\r
+ std::swap(C1, C2);\r
+ std::swap(dom1, dom2);\r
+#if VERBOSE\r
+ std::cerr << "dom(pA) : " << dompA << std::endl;\r
+ std::cerr << "dom(pB) : " << dompB << std::endl;\r
+#endif\r
+ }\r
+ domsA.push_back(dompA);\r
+ domsB.push_back(dompB);\r
+}\r
+\r
+\r
+/*\r
+ * iterate\r
+ *\r
+ * input:\r
+ * A, B: control point sets of two bezier curves\r
+ * domA, domB: real parameter intervals of the two curves\r
+ * precision: required computational precision of the returned parameter ranges\r
+ * output:\r
+ * domsA, domsB: sets of parameter intervals\r
+ *\r
+ * The parameter intervals are computed by using a Bezier clipping algorithm,\r
+ * in case the clipping doesn't shrink the initial interval more than 20%,\r
+ * a subdivision step is performed.\r
+ * If during the computation one of the two curve interval length becomes less\r
+ * than MAX_PRECISION the routine exits indipendently by the precision reached\r
+ * in the computation of the other curve interval.\r
+ */\r
+template <>\r
+void iterate<collinear_normal_tag> (std::vector<Interval>& domsA,\r
+ std::vector<Interval>& domsB,\r
+ std::vector<Point> const& A,\r
+ std::vector<Point> const& B,\r
+ Interval const& domA,\r
+ Interval const& domB,\r
+ double precision)\r
+{\r
+ // in order to limit recursion\r
+ static size_t counter = 0;\r
+ if (domA.extent() == 1 && domB.extent() == 1) counter = 0;\r
+ if (++counter > 100) return;\r
+#if VERBOSE\r
+ std::cerr << std::fixed << std::setprecision(16);\r
+ std::cerr << ">> curve subdision performed <<" << std::endl;\r
+ std::cerr << "dom(A) : " << domA << std::endl;\r
+ std::cerr << "dom(B) : " << domB << std::endl;\r
+// std::cerr << "angle(A) : " << angle(A) << std::endl;\r
+// std::cerr << "angle(B) : " << angle(B) << std::endl;\r
+#endif\r
+\r
+ if (precision < MAX_PRECISION)\r
+ precision = MAX_PRECISION;\r
+\r
+ std::vector<Point> pA = A;\r
+ std::vector<Point> pB = B;\r
+ std::vector<Point>* C1 = &pA;\r
+ std::vector<Point>* C2 = &pB;\r
+\r
+ Interval dompA = domA;\r
+ Interval dompB = domB;\r
+ Interval* dom1 = &dompA;\r
+ Interval* dom2 = &dompB;\r
+\r
+ Interval dom;\r
+\r
+ size_t iter = 0;\r
+ while (++iter < 100\r
+ && (dompA.extent() >= precision || dompB.extent() >= precision))\r
+ {\r
+#if VERBOSE\r
+ std::cerr << "iter: " << iter << std::endl;\r
+#endif\r
+ clip<collinear_normal_tag>(dom, *C1, *C2);\r
+\r
+ // [1,0] is utilized to represent an empty interval\r
+ if (dom == EMPTY_INTERVAL)\r
+ {\r
+#if VERBOSE\r
+ std::cerr << "dom: empty" << std::endl;\r
+#endif\r
+ return;\r
+ }\r
+#if VERBOSE\r
+ std::cerr << "dom : " << dom << std::endl;\r
+#endif\r
+ // all other cases where dom[0] > dom[1] are invalid\r
+ if (dom.min() > dom.max())\r
+ {\r
+ assert(dom.min() < dom.max());\r
+ }\r
+\r
+ map_to(*dom2, dom);\r
+\r
+ // it's better to stop before losing computational precision\r
+ if (iter > 1 && (dom2->extent() <= MAX_PRECISION))\r
+ {\r
+#if VERBOSE\r
+ std::cerr << "beyond max precision limit" << std::endl;\r
+#endif\r
+ break;\r
+ }\r
+\r
+ portion(*C2, dom);\r
+ if (iter > 1 && is_constant(*C2))\r
+ {\r
+#if VERBOSE\r
+ std::cerr << "new curve portion pC1 is constant" << std::endl;\r
+#endif\r
+ break;\r
+ }\r
+\r
+\r
+ // if we have clipped less than 20% than we need to subdive the curve\r
+ // with the largest domain into two sub-curves\r
+ if ( dom.extent() > MIN_CLIPPED_SIZE_THRESHOLD)\r
+ {\r
+#if VERBOSE\r
+ std::cerr << "clipped less than 20% : " << dom.extent() << std::endl;\r
+ std::cerr << "angle(pA) : " << angle(pA) << std::endl;\r
+ std::cerr << "angle(pB) : " << angle(pB) << std::endl;\r
+#endif\r
+ std::vector<Point> pC1, pC2;\r
+ Interval dompC1, dompC2;\r
+ if (dompA.extent() > dompB.extent())\r
+ {\r
+ if ((dompA.extent() / 2) < MAX_PRECISION)\r
+ {\r
+ break;\r
+ }\r
+ pC1 = pC2 = pA;\r
+ portion(pC1, H1_INTERVAL);\r
+ if (false && is_constant(pC1))\r
+ {\r
+#if VERBOSE\r
+ std::cerr << "new curve portion pC1 is constant" << std::endl;\r
+#endif\r
+ break;\r
+ }\r
+ portion(pC2, H2_INTERVAL);\r
+ if (is_constant(pC2))\r
+ {\r
+#if VERBOSE\r
+ std::cerr << "new curve portion pC2 is constant" << std::endl;\r
+#endif\r
+ break;\r
+ }\r
+ dompC1 = dompC2 = dompA;\r
+ map_to(dompC1, H1_INTERVAL);\r
+ map_to(dompC2, H2_INTERVAL);\r
+ iterate<collinear_normal_tag>(domsA, domsB, pC1, pB,\r
+ dompC1, dompB, precision);\r
+ iterate<collinear_normal_tag>(domsA, domsB, pC2, pB,\r
+ dompC2, dompB, precision);\r
+ }\r
+ else\r
+ {\r
+ if ((dompB.extent() / 2) < MAX_PRECISION)\r
+ {\r
+ break;\r
+ }\r
+ pC1 = pC2 = pB;\r
+ portion(pC1, H1_INTERVAL);\r
+ if (is_constant(pC1))\r
+ {\r
+#if VERBOSE\r
+ std::cerr << "new curve portion pC1 is constant" << std::endl;\r
+#endif\r
+ break;\r
+ }\r
+ portion(pC2, H2_INTERVAL);\r
+ if (is_constant(pC2))\r
+ {\r
+#if VERBOSE\r
+ std::cerr << "new curve portion pC2 is constant" << std::endl;\r
+#endif\r
+ break;\r
+ }\r
+ dompC1 = dompC2 = dompB;\r
+ map_to(dompC1, H1_INTERVAL);\r
+ map_to(dompC2, H2_INTERVAL);\r
+ iterate<collinear_normal_tag>(domsB, domsA, pC1, pA,\r
+ dompC1, dompA, precision);\r
+ iterate<collinear_normal_tag>(domsB, domsA, pC2, pA,\r
+ dompC2, dompA, precision);\r
+ }\r
+ return;\r
+ }\r
+\r
+ std::swap(C1, C2);\r
+ std::swap(dom1, dom2);\r
+#if VERBOSE\r
+ std::cerr << "dom(pA) : " << dompA << std::endl;\r
+ std::cerr << "dom(pB) : " << dompB << std::endl;\r
+#endif\r
+ }\r
+ domsA.push_back(dompA);\r
+ domsB.push_back(dompB);\r
+}\r
+\r
+\r
+/*\r
+ * get_solutions\r
+ *\r
+ * input: A, B - set of control points of two Bezier curve\r
+ * input: precision - required precision of computation\r
+ * input: clip - the routine used for clipping\r
+ * output: xs - set of pairs of parameter values\r
+ * at which the clipping algorithm converges\r
+ *\r
+ * This routine is based on the Bezier Clipping Algorithm,\r
+ * see: Sederberg - Computer Aided Geometric Design\r
+ */\r
+template <typename Tag>\r
+void get_solutions (std::vector< std::pair<double, double> >& xs,\r
+ std::vector<Point> const& A,\r
+ std::vector<Point> const& B,\r
+ double precision)\r
+{\r
+ std::pair<double, double> ci;\r
+ std::vector<Interval> domsA, domsB;\r
+ iterate<Tag> (domsA, domsB, A, B, UNIT_INTERVAL, UNIT_INTERVAL, precision);\r
+ if (domsA.size() != domsB.size())\r
+ {\r
+ assert (domsA.size() == domsB.size());\r
+ }\r
+ xs.clear();\r
+ xs.reserve(domsA.size());\r
+ for (size_t i = 0; i < domsA.size(); ++i)\r
+ {\r
+#if VERBOSE\r
+ std::cerr << i << " : domA : " << domsA[i] << std::endl;\r
+ std::cerr << "extent A: " << domsA[i].extent() << " ";\r
+ std::cerr << "precision A: " << get_precision(domsA[i]) << std::endl;\r
+ std::cerr << i << " : domB : " << domsB[i] << std::endl;\r
+ std::cerr << "extent B: " << domsB[i].extent() << " ";\r
+ std::cerr << "precision B: " << get_precision(domsB[i]) << std::endl;\r
+#endif\r
+ ci.first = domsA[i].middle();\r
+ ci.second = domsB[i].middle();\r
+ xs.push_back(ci);\r
+ }\r
+}\r
+\r
+} /* end namespace bezier_clipping */ } /* end namespace detail */\r
+\r
+\r
+/*\r
+ * find_collinear_normal\r
+ *\r
+ * input: A, B - set of control points of two Bezier curve\r
+ * input: precision - required precision of computation\r
+ * output: xs - set of pairs of parameter values\r
+ * at which there are collinear normals\r
+ *\r
+ * This routine is based on the Bezier Clipping Algorithm,\r
+ * see: Sederberg, Nishita, 1990 - Curve intersection using Bezier clipping\r
+ */\r
+void find_collinear_normal (std::vector< std::pair<double, double> >& xs,\r
+ std::vector<Point> const& A,\r
+ std::vector<Point> const& B,\r
+ double precision)\r
+{\r
+ using detail::bezier_clipping::get_solutions;\r
+ using detail::bezier_clipping::collinear_normal_tag;\r
+ get_solutions<collinear_normal_tag>(xs, A, B, precision);\r
+}\r
+\r
+\r
+/*\r
+ * find_intersections_bezier_clipping\r
+ *\r
+ * input: A, B - set of control points of two Bezier curve\r
+ * input: precision - required precision of computation\r
+ * output: xs - set of pairs of parameter values\r
+ * at which crossing happens\r
+ *\r
+ * This routine is based on the Bezier Clipping Algorithm,\r
+ * see: Sederberg, Nishita, 1990 - Curve intersection using Bezier clipping\r
+ */\r
+void find_intersections_bezier_clipping (std::vector< std::pair<double, double> >& xs,\r
+ std::vector<Point> const& A,\r
+ std::vector<Point> const& B,\r
+ double precision)\r
+{\r
+ using detail::bezier_clipping::get_solutions;\r
+ using detail::bezier_clipping::intersection_point_tag;\r
+ get_solutions<intersection_point_tag>(xs, A, B, precision);\r
+}\r
+\r
+} // end namespace Geom\r
+\r
+\r
+\r
+\r
+/*\r
+ Local Variables:\r
+ mode:c++\r
+ c-file-style:"stroustrup"\r
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))\r
+ indent-tabs-mode:nil\r
+ fill-column:99\r
+ End:\r
+*/\r
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :\r
index acf878ee89fddb768b3ee45928ea5bcf0dbd6724..d5259c71f61ade37848c799444b3bf263c54e24b 100644 (file)
--- a/src/2geom/bezier-curve.h
+++ b/src/2geom/bezier-curve.h
return ret;
}
- Curve *derivative() const {
- if(order > 1)
- return new BezierCurve<order-1>(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<Point> pointAndDerivatives(Coord t, unsigned n) const { return inner.valueAndDerivatives(t, n); }
};
// 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;
return distance(_segment.initialPoint(), _segment.finalPoint());
}
+template <unsigned order>
+inline
+Curve *BezierCurve<order>::derivative() const {
+ return new BezierCurve<order-1>(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 889dde9edb43a825e8d1c9df2c3f740ce86ef26e..4ab965f42ecd5d2b036bf25d2f87af7f3ef12291 100644 (file)
--- a/src/2geom/bezier.h
+++ b/src/2geom/bezier.h
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<Coord> 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]);
for (unsigned j = 0; j <= order; j++)
right[j] = vtemp[order-j][j];
- return (vtemp[order][0]);
+ return (vtemp[order][0]);*/
}
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<n; i++){
+ tn = tn*t;
+ bc = bc*(n-i+1)/i;
+ tmp = (tmp + tn*bc*c_[i])*u;
+ }
+ return (tmp + tn*t*c_[n]);
+ //return subdivideArr(t, &c_[0], NULL, NULL, order());
}
inline Coord operator()(double t) const { return valueAt(t); }
//Only mutator
inline Coord &operator[](unsigned ix) { return c_[ix]; }
- inline Coord const &operator[](unsigned ix) const { return c_[ix]; }
+ inline Coord const &operator[](unsigned ix) const { return const_cast<std::valarray<Coord>&>(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
std::vector<Coord> valueAndDerivatives(Coord t, unsigned n_derivs) const {
std::vector<Coord> val_n_der;
- Coord d_[order()+1];
+ std::valarray<Coord> 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]);
}
std::pair<Bezier, Bezier > 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<std::valarray<Coord>&>(c_)[0], &a.c_[0], &b.c_[0], order());
return std::pair<Bezier, Bezier >(a, b);
}
std::vector<double> roots() const {
std::vector<double> solutions;
- find_bernstein_roots(&c_[0], order(), solutions, 0, 0.0, 1.0);
+ find_bernstein_roots(&const_cast<std::valarray<Coord>&>(c_)[0], order(), solutions, 0, 0.0, 1.0);
return solutions;
}
};
inline Bezier portion(const Bezier & a, double from, double to) {
//TODO: implement better?
- Coord res[a.order()+1];
+ std::valarray<Coord> 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<Bezier&>(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<Bezier&>(a).c_[0], NULL, &res[0], a.order());
+ if(to == 1) return Bezier(&res[0], a.order());
+ std::valarray<Coord> 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
}
inline OptInterval bounds_fast(Bezier const & b) {
- return Interval::fromArray(&b.c_[0], b.size());
+ return Interval::fromArray(&const_cast<Bezier&>(b).c_[0], b.size());
}
//TODO: better bounds exact
index 447c5183f920aeab0913118741f7b20b5dc84d77..4ba3e2325a3b2c1b362da8059e1fe88e4f3c1e03 100644 (file)
--- a/src/2geom/chebyshev.cpp
+++ b/src/2geom/chebyshev.cpp
-#include <2geom/chebyshev.h>
-
-#include <2geom/sbasis.h>
-#include <2geom/sbasis-poly.h>
-
-#include <vector>
-using std::vector;
-
-#include <gsl/gsl_math.h>
-#include <gsl/gsl_chebyshev.h>
-
-namespace Geom{
-
-SBasis cheb(unsigned n) {
- static std::vector<SBasis> 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<SBasis> 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>\r
+\r
+#include <2geom/sbasis.h>\r
+#include <2geom/sbasis-poly.h>\r
+\r
+#include <vector>\r
+using std::vector;\r
+\r
+#include <gsl/gsl_math.h>\r
+#include <gsl/gsl_chebyshev.h>\r
+\r
+namespace Geom{\r
+\r
+SBasis cheb(unsigned n) {\r
+ static std::vector<SBasis> basis;\r
+ if(basis.empty()) {\r
+ basis.push_back(Linear(1,1));\r
+ basis.push_back(Linear(0,1));\r
+ }\r
+ for(unsigned i = basis.size(); i <= n; i++) {\r
+ basis.push_back(Linear(0,2)*basis[i-1] - basis[i-2]);\r
+ }\r
+ \r
+ return basis[n];\r
+}\r
+\r
+SBasis cheb_series(unsigned n, double* cheb_coeff) {\r
+ SBasis r;\r
+ for(unsigned i = 0; i < n; i++) {\r
+ double cof = cheb_coeff[i];\r
+ //if(i == 0)\r
+ //cof /= 2;\r
+ r += cheb(i)*cof;\r
+ }\r
+ \r
+ return r;\r
+}\r
+\r
+SBasis clenshaw_series(unsigned m, double* cheb_coeff) {\r
+ /** b_n = a_n\r
+ b_n-1 = 2*x*b_n + a_n-1\r
+ b_n-k = 2*x*b_{n-k+1} + a_{n-k} - b_{n - k + 2}\r
+ b_0 = x*b_1 + a_0 - b_2\r
+ */\r
+ \r
+ double a = -1, b = 1;\r
+ SBasis d, dd;\r
+ SBasis y = (Linear(0, 2) - (a+b)) / (b-a);\r
+ SBasis y2 = 2*y;\r
+ for(int j = m - 1; j >= 1; j--) {\r
+ SBasis sv = d;\r
+ d = y2*d - dd + cheb_coeff[j];\r
+ dd = sv;\r
+ }\r
+ \r
+ return y*d - dd + 0.5*cheb_coeff[0];\r
+}\r
+\r
+SBasis chebyshev_approximant (double (*f)(double,void*), int order, Interval in, void* p) {\r
+ gsl_cheb_series *cs = gsl_cheb_alloc (order+2);\r
+\r
+ gsl_function F;\r
+\r
+ F.function = f;\r
+ F.params = p;\r
+\r
+ gsl_cheb_init (cs, &F, in[0], in[1]);\r
+ \r
+ SBasis r = compose(clenshaw_series(order, cs->c), Linear(-1,1));\r
+ \r
+ gsl_cheb_free (cs);\r
+ return r;\r
+}\r
+\r
+struct wrap {\r
+ double (*f)(double,void*);\r
+ void* pp;\r
+ double fa, fb;\r
+ Interval in;\r
+};\r
+\r
+double f_interp(double x, void* p) {\r
+ struct wrap *wr = (struct wrap *)p;\r
+ double z = (x - wr->in[0]) / (wr->in[1] - wr->in[0]);\r
+ return (wr->f)(x, wr->pp) - ((1 - z)*wr->fa + z*wr->fb);\r
+}\r
+\r
+SBasis chebyshev_approximant_interpolating (double (*f)(double,void*), \r
+ int order, Interval in, void* p) {\r
+ double fa = f(in[0], p);\r
+ double fb = f(in[1], p);\r
+ struct wrap wr;\r
+ wr.fa = fa;\r
+ wr.fb = fb;\r
+ wr.in = in;\r
+ printf("%f %f\n", fa, fb);\r
+ wr.f = f;\r
+ wr.pp = p;\r
+ return compose(Linear(in[0], in[1]), Linear(fa, fb)) + chebyshev_approximant(f_interp, order, in, &wr) + Linear(fa, fb);\r
+}\r
+\r
+SBasis chebyshev(unsigned n) {\r
+ static std::vector<SBasis> basis;\r
+ if(basis.empty()) {\r
+ basis.push_back(Linear(1,1));\r
+ basis.push_back(Linear(0,1));\r
+ }\r
+ for(unsigned i = basis.size(); i <= n; i++) {\r
+ basis.push_back(Linear(0,2)*basis[i-1] - basis[i-2]);\r
+ }\r
+ \r
+ return basis[n];\r
+}\r
+\r
+};\r
+\r
+/*\r
+ Local Variables:\r
+ mode:c++\r
+ c-file-style:"stroustrup"\r
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))\r
+ indent-tabs-mode:nil\r
+ fill-column:99\r
+ End:\r
+*/\r
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :\r
diff --git a/src/2geom/chebyshev.h b/src/2geom/chebyshev.h
index 6de9e9cc0ed1f742e3182ee2dd3939ec56f60c9f..309e09b3e788533c187d8cea4f436540f493e02e 100644 (file)
--- a/src/2geom/chebyshev.h
+++ b/src/2geom/chebyshev.h
-#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\r
+#define _CHEBYSHEV\r
+\r
+#include <2geom/sbasis.h>\r
+#include <2geom/interval.h>\r
+\r
+/*** Conversion between Chebyshev approximation and SBasis.\r
+ * \r
+ */\r
+\r
+namespace Geom{\r
+\r
+SBasis chebyshev_approximant (double (*f)(double,void*), int order, Interval in, void* p=0);\r
+SBasis chebyshev_approximant_interpolating (double (*f)(double,void*), int order, Interval in, void* p=0);\r
+SBasis chebyshev(unsigned n);\r
+\r
+};\r
+\r
+/*\r
+ Local Variables:\r
+ mode:c++\r
+ c-file-style:"stroustrup"\r
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))\r
+ indent-tabs-mode:nil\r
+ fill-column:99\r
+ End:\r
+*/\r
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :\r
+\r
+#endif\r
index f5a0f9cd82dffcd00b0c3bebb002a280ab8e1c7b..9c4ea7776d232f12a6e195ba6b35debf0e4a11af 100644 (file)
unsigned n = vec.size();
unsigned m = result.size();
assert(m*n == matrix.size());
- const double* mp = &matrix[0];
+ const double* mp = &const_cast<valarray<double>&>(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 7737b24a9685a8b7d3d128586cbdedbf9c449bdf..427848033dc22150733afaad60d43bd5e8785f26 100644 (file)
--- 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 f55f6a5969fe770de62134c9f8e904b66296f0e6..afa00b40d2e3b5f00af3268682ba3181a21c756f 100644 (file)
--- a/src/2geom/d2.h
+++ b/src/2geom/d2.h
return D2<T>(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 <typename T>
+inline std::ostream &operator<< (std::ostream &out_file, const Geom::D2<T> &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 e04c35c29d0df80fd22bd7b6ca9f185b967dc31f..a7e6a54bb0a2ce4b5c32caf0405892e1976b347b 100644 (file)
--- a/src/2geom/line.h
+++ b/src/2geom/line.h
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
{
index 47c588cb22fcec488fdb91a336dbecc7a973eb41..d6a26bd2d26a1eeff1fcd69a76403d9b629476fd 100644 (file)
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<Point> const& old_sample_values,
+ std::vector<Point> 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<ModelT, double, true>( model_type const& _model,
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<Point> const& old_sample_values,
+ std::vector<Point> 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<ModelT, Point, true>( model_type const& _model,
index e9b6e6e9e13b53c8f4a9962a38b20559653d97ca..97db59d56a4c18bdd336f03927ef681d714679f9 100644 (file)
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;
index bfb1fb96c508f1820fbc3742aa87f5608598a0bf..b2d5ceabb7a79d63f275f6805f6377e69c1c2c77 100644 (file)
//for path_direction:
#include <2geom/sbasis-geometric.h>
+#include <2geom/line.h>
+#ifdef HAVE_GSL
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multiroots.h>
+#endif
namespace Geom {
}
#endif
+#ifdef HAVE_GSL
struct rparams {
Curve const &A;
Curve const &B;
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<Point> 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
}
/**
*/
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;
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;
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
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)) {
ret.push_back(Crossing(tA, tB, c > 0));
return;
}
- //}
+ }
if(depth > 12) return;
double mid = (Bl + Bh)/2;
mono_pair(B, Bl, mid,
/** This returns the times when the x or y derivative is 0 in the curve. */
std::vector<double> curve_mono_splits(Curve const &d) {
+ Curve* deriv = d.derivative();
std::vector<double> 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<double> offset_doubles(std::vector<double> const &x, double offs) {
std::vector<double> path_mono_splits(Path const &p) {
std::vector<double> 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<double> spl = offset_doubles(curve_mono_splits(*deriv), i);
- delete deriv;
+ std::vector<double> 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] :
index 4eef168237a59c00a6ce04726aa41c56de164e9e..6457b5e431d2ab7f2bcfe0c4354594521f0421f8 100644 (file)
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<Path> {
Crossings crossings(Curve const &a, Curve const &b);
Crossings crossings(Path const &a, Path const &b) { return curve_sweep<SimpleCrosser>(a, b); }
diff --git a/src/2geom/path.cpp b/src/2geom/path.cpp
index 136e6d4818fee5048c70065a9aeecb9f2a67dd05..981c9f044717747981ad3c7a85a3d42d5495c7cc 100644 (file)
--- 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 144d7a9b22d098b2c14e76cb6194dd9c8b1eeb5b..a5be42587c0be602eb02c28677769e14943bfed7 100644 (file)
--- a/src/2geom/piecewise.h
+++ b/src/2geom/piecewise.h
return ret;
}
template<typename T>
+Piecewise<T> operator*(Piecewise<T> const &a, T b) {
+ boost::function_requires<ScalableConcept<T> >();
+
+ if(a.empty()) return Piecewise<T>();
+
+ Piecewise<T> 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<typename T>
Piecewise<T> operator/(Piecewise<T> const &a, double b) {
boost::function_requires<ScalableConcept<T> >();
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 30ec1cdb74686ed1f4197350aba3188363b3804d..d8b379557f154b812dd151dd62d3a50df3ac73e3 100644 (file)
--- a/src/2geom/poly.cpp
+++ b/src/2geom/poly.cpp
#include <2geom/poly.h>
-//#ifdef HAVE_GSL
+#define HAVE_GSL
+#ifdef HAVE_GSL
#include <gsl/gsl_poly.h>
-//#endif
+#endif
namespace Geom {
}
-//#ifdef HAVE_GSL
+#ifdef HAVE_GSL
std::vector<std::complex<double> > solve(Poly const & pp) {
Poly p(pp);
p.normalize();
}
return real_roots;
}
-//#endif
+#endif
double polish_root(Poly const & p, double guess, double tol) {
Poly dp = derivative(p);
diff --git a/src/2geom/quadtree.cpp b/src/2geom/quadtree.cpp
index e322a091b8b10b4a27036deed066b0612f6ba226..211590baeff75dcc69de5e28f01deff416525ed5 100644 (file)
--- a/src/2geom/quadtree.cpp
+++ b/src/2geom/quadtree.cpp
}
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;
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;
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;
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) {
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 5338698cff0efc6e8eb0071f602685d3107227c2..2e114a0a08fb101a3234d60cdcc12370e57598d0 100644 (file)
--- a/src/2geom/quadtree.h
+++ b/src/2geom/quadtree.h
Quad* search(Geom::Rect const &r);
void insert(Geom::Rect const &r, int shape);
void erase(Quad *q, int shape);
+private:
+ bool clean_root();
};
};
index 7053388dce013ce1635bbc1571e3e9e266247948..f3a984c967636cc26435adeb72cc10fbe506ddca 100644 (file)
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<r.size(); i++) {
+ for(unsigned i = 1; int(i) <= order && i<r.size(); i++) {
Linear ci(r[i][0]/(2*sqrtf[0][0]), r[i][1]/(2*sqrtf[0][1]));
SBasis cisi = shift(ci, i);
r -= multiply(shift((sqrtf*2 + cisi), i), SBasis(ci));
index 5249053fa5955e3e0c1b9c0a6a6aa94968046a9c..37e07cbe8a69f96d8344f985d8fa6dd18d725fc3 100644 (file)
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 49678b19ab56e0a21878bdb3505deae115e030d1..2f7f03bfcdf81df230521b27940862f452f9150c 100644 (file)
--- a/src/2geom/sbasis.cpp
+++ b/src/2geom/sbasis.cpp
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<r.size(); i++) {
+ for(unsigned i = 1; i <= (unsigned)k && i<r.size(); i++) {
Linear ci(r[i][0]/(2*c[0][0]), r[i][1]/(2*c[0][1]));
SBasis cisi = shift(ci, i);
r -= multiply(shift((c*2 + cisi), i), SBasis(ci));
diff --git a/src/2geom/sbasis.h b/src/2geom/sbasis.h
index ae045b2229aad25f051630459940fd8cc4a0f25a..a32823f133a6d838bf65a226d81d514b2d2321da 100644 (file)
--- a/src/2geom/sbasis.h
+++ b/src/2geom/sbasis.h
//void insert(Linear* aa, Linear* bb, Linear* cc} { d.insert(aa, bb, cc);}
Linear& at(unsigned i) { return d.at(i);}
//void insert(Linear* before, int& n, Linear const &l) { d.insert(std::vector<Linear>::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<Linear>() { return d;}
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)
{}
index 9a7b5633fff5d37fba9ecc8e9115ec28e8b03e5a..a1c0ca5578853b662d5219a333c4a61bcdedeb9e 100644 (file)
#include <2geom/solver.h>
#include <2geom/point.h>
#include <algorithm>
+#include <valarray>
/*** 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
}
/* Otherwise, solve recursively after subdividing control polygon */
- double Left[N], /* New left and right */
- Right[N]; /* control polygons */
+ std::valarray<double> new_controls(2*N); // New left and right control polygons
const double t = 0.5;
/* 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 227c47e098f41da0a621cc78ed0ac175ba9773c4..7571efe094b671c8080bcd64927d102cb43564db 100644 (file)
--- a/src/2geom/sweep.cpp
+++ b/src/2geom/sweep.cpp
* \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<std::vector<unsigned> > sweep_bounds(std::vector<Rect> rs) {
+std::vector<std::vector<unsigned> > sweep_bounds(std::vector<Rect> rs, Dim2 d) {
std::vector<Event> events; events.reserve(rs.size()*2);
std::vector<std::vector<unsigned> > 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());
} 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);
}
}
*
* \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<std::vector<unsigned> > sweep_bounds(std::vector<Rect> a, std::vector<Rect> b) {
+std::vector<std::vector<unsigned> > sweep_bounds(std::vector<Rect> a, std::vector<Rect> b, Dim2 d) {
std::vector<std::vector<unsigned> > pairs(a.size());
if(a.empty() || b.empty()) return pairs;
std::vector<Event> events[2];
@@ -64,15 +66,17 @@ std::vector<std::vector<unsigned> > sweep_bounds(std::vector<Rect> 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<unsigned> 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";
//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);
}
}
//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<std::vector<unsigned> > sweep_bounds(std::vector<Rect> 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 0214511ac328cdf85ec5370d95f7a46dbe7ca6ba..9d1643d7af83edd5e0970cc8300b36525d8a0a21 100644 (file)
--- a/src/2geom/sweep.h
+++ b/src/2geom/sweep.h
if(x > other.x) return false;
return closing < other.closing;
}
-
};
-std::vector<std::vector<unsigned> > sweep_bounds(std::vector<Rect>);
-std::vector<std::vector<unsigned> > sweep_bounds(std::vector<Rect>, std::vector<Rect>);
+std::vector<std::vector<unsigned> >
+sweep_bounds(std::vector<Rect>, Dim2 dim = X);
+
+std::vector<std::vector<unsigned> >
+sweep_bounds(std::vector<Rect>, std::vector<Rect>, Dim2 dim = X);
std::vector<std::vector<unsigned> > fake_cull(unsigned a, unsigned b);
diff --git a/src/2geom/utils.cpp b/src/2geom/utils.cpp
index 579718553503dd13f0509e839468dff1cadbe44f..f0f42ce0179e8ac33ed7a5e085207e0298e394e3 100644 (file)
--- a/src/2geom/utils.cpp
+++ b/src/2geom/utils.cpp
-/** Various utility functions.
- *
- * Copyright 2008 Marco Cecchetti <mrcekets at gmail.com>
- * Copyright 2007 Johan Engelen <goejendaagh@zonnet.nl>
- * Copyright 2006 Michael G. Sloan <mgsloan@gmail.com>
- *
- * 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<size_t>& 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.\r
+ *\r
+ * Copyright 2008 Marco Cecchetti <mrcekets at gmail.com>\r
+ * Copyright 2007 Johan Engelen <goejendaagh@zonnet.nl>\r
+ * Copyright 2006 Michael G. Sloan <mgsloan@gmail.com>\r
+ *\r
+ * This library is free software; you can redistribute it and/or\r
+ * modify it either under the terms of the GNU Lesser General Public\r
+ * License version 2.1 as published by the Free Software Foundation\r
+ * (the "LGPL") or, at your option, under the terms of the Mozilla\r
+ * Public License Version 1.1 (the "MPL"). If you do not alter this\r
+ * notice, a recipient may use your version of this file under either\r
+ * the MPL or the LGPL.\r
+ *\r
+ * You should have received a copy of the LGPL along with this library\r
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software\r
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
+ * You should have received a copy of the MPL along with this library\r
+ * in the file COPYING-MPL-1.1\r
+ *\r
+ * The contents of this file are subject to the Mozilla Public License\r
+ * Version 1.1 (the "License"); you may not use this file except in\r
+ * compliance with the License. You may obtain a copy of the License at\r
+ * http://www.mozilla.org/MPL/\r
+ *\r
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY\r
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for\r
+ * the specific language governing rights and limitations.\r
+ *\r
+ */\r
+\r
+\r
+#include <2geom/utils.h>\r
+\r
+\r
+namespace Geom \r
+{\r
+\r
+// return a vector that contains all the binomial coefficients of degree n \r
+void binomial_coefficients(std::vector<size_t>& bc, size_t n)\r
+{\r
+ size_t s = n+1;\r
+ bc.clear();\r
+ bc.resize(s);\r
+ bc[0] = 1;\r
+ size_t k;\r
+ for (size_t i = 1; i < n; ++i)\r
+ {\r
+ k = i >> 1;\r
+ if (i & 1u)\r
+ {\r
+ bc[k+1] = bc[k] << 1;\r
+ }\r
+ for (size_t j = k; j > 0; --j)\r
+ {\r
+ bc[j] += bc[j-1];\r
+ }\r
+ }\r
+ s >>= 1;\r
+ for (size_t i = 0; i < s; ++i)\r
+ {\r
+ bc[n-i] = bc[i];\r
+ }\r
+}\r
+\r
+} // end namespace Geom\r
+\r
+\r
+\r
+\r
+\r
+\r
+\r
+\r
+\r
+\r
+\r
+/*\r
+ Local Variables:\r
+ mode:c++\r
+ c-file-style:"stroustrup"\r
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))\r
+ indent-tabs-mode:nil\r
+ fill-column:99\r
+ End:\r
+*/\r
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :\r