diff options
Diffstat (limited to 'src/2geom/path-intersection.cpp')
| -rw-r--r-- | src/2geom/path-intersection.cpp | 728 |
1 files changed, 728 insertions, 0 deletions
diff --git a/src/2geom/path-intersection.cpp b/src/2geom/path-intersection.cpp new file mode 100644 index 0000000..280d7ba --- /dev/null +++ b/src/2geom/path-intersection.cpp @@ -0,0 +1,728 @@ +#include <2geom/path-intersection.h> + +#include <2geom/ord.h> + +//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 { + +/// Compute winding number of the path at the specified point +int winding(Path const &path, Point const &p) { + return path.winding(p); +} + +/** + * This function should only be applied to simple paths (regions), as otherwise + * a boolean winding direction is undefined. It returns true for fill, false for + * hole. Defaults to using the sign of area when it reaches funny cases. + */ +bool path_direction(Path const &p) { + if(p.empty()) return false; + + /*goto doh; + //could probably be more efficient, but this is a quick job + double y = p.initialPoint()[Y]; + double x = p.initialPoint()[X]; + Cmp res = cmp(p[0].finalPoint()[Y], y); + for(unsigned i = 1; i < p.size(); i++) { + Cmp final_to_ray = cmp(p[i].finalPoint()[Y], y); + Cmp initial_to_ray = cmp(p[i].initialPoint()[Y], y); + // if y is included, these will have opposite values, giving order. + Cmp c = cmp(final_to_ray, initial_to_ray); + if(c != EQUAL_TO) { + std::vector<double> rs = p[i].roots(y, Y); + for(unsigned j = 0; j < rs.size(); j++) { + double nx = p[i].valueAt(rs[j], X); + if(nx > x) { + x = nx; + res = c; + } + } + } else if(final_to_ray == EQUAL_TO) goto doh; + } + return res < 0; + + doh:*/ + //Otherwise fallback on area + + Piecewise<D2<SBasis> > pw = p.toPwSb(); + double area; + Point centre; + Geom::centroid(pw, centre, area); + return area > 0; +} + +//pair intersect code based on njh's pair-intersect + +/** A little sugar for appending a list to another */ +template<typename T> +void append(T &a, T const &b) { + a.insert(a.end(), b.begin(), b.end()); +} + +/** + * Finds the intersection between the lines defined by A0 & A1, and B0 & B1. + * Returns through the last 3 parameters, returning the t-values on the lines + * and the cross-product of the deltas (a useful byproduct). The return value + * indicates if the time values are within their proper range on the line segments. + */ +bool +linear_intersect(Point const &A0, Point const &A1, Point const &B0, Point const &B1, + double &tA, double &tB, double &det) { + bool both_lines_non_zero = (!are_near(A0, A1)) && (!are_near(B0, B1)); + + // Cramer's rule as cross products + Point Ad = A1 - A0, + Bd = B1 - B0, + d = B0 - A0; + det = cross(Ad, Bd); + + double det_rel = det; // Calculate the determinant of the normalized vectors + if (both_lines_non_zero) { + det_rel /= Ad.length(); + det_rel /= Bd.length(); + } + + if( fabs(det_rel) < 1e-12 ) { // If the cross product is NEARLY zero, + // Then one of the linesegments might have length zero + if (both_lines_non_zero) { + // If that's not the case, then we must have either: + // - parallel lines, with no intersections, or + // - coincident lines, with an infinite number of intersections + // Either is quite useless, so we'll just bail out + return false; + } // Else, one of the linesegments is zero, and we might still be able to calculate a single intersection point + } // Else we haven't bailed out, and we'll try to calculate the intersections + + double detinv = 1.0 / det; + tA = cross(d, Bd) * detinv; + tB = cross(d, Ad) * detinv; + return (tA >= 0.) && (tA <= 1.) && (tB >= 0.) && (tB <= 1.); +} + + +#if 0 +typedef union dbl_64{ + long long i64; + double d64; +}; + +static double EpsilonOf(double value) +{ + dbl_64 s; + s.d64 = value; + if(s.i64 == 0) + { + s.i64++; + return s.d64 - value; + } + else if(s.i64-- < 0) + return s.d64 - value; + else + return value - s.d64; +} +#endif + +#ifdef HAVE_GSL +struct rparams { + Curve const &A; + Curve const &B; +}; + +static int +intersect_polish_f (const gsl_vector * x, void *params, + gsl_vector * f) +{ + const double x0 = gsl_vector_get (x, 0); + const double x1 = gsl_vector_get (x, 1); + + Geom::Point dx = ((struct rparams *) params)->A(x0) - + ((struct rparams *) params)->B(x1); + + gsl_vector_set (f, 0, dx[0]); + gsl_vector_set (f, 1, dx[1]); + + return GSL_SUCCESS; +} +#endif + +static void +intersect_polish_root (Curve const &A, double &s, Curve const &B, double &t) +{ + 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. + + Affine 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 + int status; + size_t iter = 0; + const size_t n = 2; + struct rparams p = {A, B}; + gsl_multiroot_function f = {&intersect_polish_f, n, &p}; + + double x_init[2] = {s, t}; + gsl_vector *x = gsl_vector_alloc (n); + + gsl_vector_set (x, 0, x_init[0]); + gsl_vector_set (x, 1, x_init[1]); + + 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); + + if (status) /* check if solver is stuck */ + break; + + status = + gsl_multiroot_test_residual (sol->f, 1e-12); + } + while (status == GSL_CONTINUE && iter < 1000); + + s = gsl_vector_get (sol->x, 0); + t = gsl_vector_get (sol->x, 1); + + gsl_multiroot_fsolver_free (sol); + gsl_vector_free (x); + } +#endif +} + +/** + * This uses the local bounds functions of curves to generically intersect two. + * It passes in the curves, time intervals, and keeps track of depth, while + * returning the results through the Crossings parameter. + */ +void pair_intersect(Curve const & A, double Al, double Ah, + Curve const & B, double Bl, double Bh, + Crossings &ret, unsigned depth = 0) { + // std::cout << depth << "(" << Al << ", " << Ah << ")\n"; + OptRect Ar = A.boundsLocal(Interval(Al, Ah)); + if (!Ar) return; + + OptRect Br = B.boundsLocal(Interval(Bl, Bh)); + if (!Br) return; + + if(! Ar->intersects(*Br)) 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)) { + 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; + pair_intersect(B, Bl, mid, + A, Al, Ah, + ret, depth+1); + pair_intersect(B, mid, Bh, + A, Al, Ah, + ret, depth+1); +} + +Crossings pair_intersect(Curve const & A, Interval const &Ad, + Curve const & B, Interval const &Bd) { + Crossings ret; + pair_intersect(A, Ad.min(), Ad.max(), B, Bd.min(), Bd.max(), ret); + return ret; +} + +/** A simple wrapper around pair_intersect */ +Crossings SimpleCrosser::crossings(Curve const &a, Curve const &b) { + Crossings ret; + pair_intersect(a, 0, 1, b, 0, 1, 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 + * doesn't exist or is inaccurate. Uses the fact that it is monotonic to + * do very fast local bounds. + */ +void mono_pair(Path const &A, double Al, double Ah, + Path const &B, double Bl, double Bh, + Crossings &ret, double /*tol*/, 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() < 0.1 && Ar.maxExtent() < 0.1)) { + double tA, tB, c; + if(linear_intersect(A0, A1, B0, B1, + tA, tB, c)) { + tA = tA * (Ah - Al) + Al; + tB = tB * (Bh - Bl) + Bl; + 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_pair(B, Bl, mid, + A, Al, Ah, + ret, depth+1); + mono_pair(B, mid, Bh, + A, Al, Ah, + ret, depth+1); +} + +/** 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 = deriv->roots(0, X); + append(rs, deriv->roots(0, Y)); + delete deriv; + std::sort(rs.begin(), rs.end()); + return rs; +} + +/** Convenience function to add a value to each entry in a vector of doubles. */ +std::vector<double> offset_doubles(std::vector<double> const &x, double offs) { + std::vector<double> ret; + for(double i : x) { + ret.push_back(i + offs); + } + return ret; +} + +/** + * Finds all the monotonic splits for a path. Only includes the split between + * curves if they switch derivative directions at that point. + */ +std::vector<double> path_mono_splits(Path const &p) { + std::vector<double> ret; + if(p.empty()) return ret; + + int pdx = 2, pdy = 2; // Previous derivative direction + for(unsigned i = 0; i < p.size(); i++) { + std::vector<double> spl = offset_doubles(curve_mono_splits(p[i]), i); + int dx = p[i].initialPoint()[X] > (spl.empty() ? p[i].finalPoint()[X] : p.valueAt(spl.front(), X)) ? 1 : 0; + int dy = p[i].initialPoint()[Y] > (spl.empty() ? p[i].finalPoint()[Y] : p.valueAt(spl.front(), Y)) ? 1 : 0; + //The direction changed, include the split time + if(dx != pdx || dy != pdy) { + ret.push_back(i); + pdx = dx; pdy = dy; + } + append(ret, spl); + } + return ret; +} + +/** + * Applies path_mono_splits to multiple paths, and returns the results such that + * time-set i corresponds to Path i. + */ +std::vector<std::vector<double> > paths_mono_splits(PathVector const &ps) { + std::vector<std::vector<double> > ret; + for(const auto & p : ps) + ret.push_back(path_mono_splits(p)); + return ret; +} + +/** + * Processes the bounds for a list of paths and a list of splits on them, yielding a list of rects for each. + * Each entry i corresponds to path i of the input. The number of rects in each entry is guaranteed to be the + * number of splits for that path, subtracted by one. + */ +std::vector<std::vector<Rect> > split_bounds(PathVector const &p, std::vector<std::vector<double> > splits) { + std::vector<std::vector<Rect> > ret; + for(unsigned i = 0; i < p.size(); i++) { + std::vector<Rect> res; + for(unsigned j = 1; j < splits[i].size(); j++) + res.emplace_back(p[i].pointAt(splits[i][j-1]), p[i].pointAt(splits[i][j])); + ret.push_back(res); + } + return ret; +} + +/** + * This is the main routine of "MonoCrosser", and implements a monotonic strategy on multiple curves. + * Finds crossings between two sets of paths, yielding a CrossingSet. [0, a.size()) of the return correspond + * to the sorted crossings of a with paths of b. The rest of the return, [a.size(), a.size() + b.size()], + * corresponds to the sorted crossings of b with paths of a. + * + * This function does two sweeps, one on the bounds of each path, and after that cull, one on the curves within. + * This leads to a certain amount of code complexity, however, most of that is factored into the above functions + */ +CrossingSet MonoCrosser::crossings(PathVector const &a, PathVector const &b) { + if(b.empty()) return CrossingSet(a.size(), Crossings()); + CrossingSet results(a.size() + b.size(), Crossings()); + if(a.empty()) return results; + + std::vector<std::vector<double> > splits_a = paths_mono_splits(a), splits_b = paths_mono_splits(b); + std::vector<std::vector<Rect> > bounds_a = split_bounds(a, splits_a), bounds_b = split_bounds(b, splits_b); + + std::vector<Rect> bounds_a_union, bounds_b_union; + for(auto & i : bounds_a) bounds_a_union.push_back(union_list(i)); + for(auto & i : bounds_b) bounds_b_union.push_back(union_list(i)); + + std::vector<std::vector<unsigned> > cull = sweep_bounds(bounds_a_union, bounds_b_union); + Crossings n; + for(unsigned i = 0; i < cull.size(); i++) { + for(unsigned jx = 0; jx < cull[i].size(); jx++) { + unsigned j = cull[i][jx]; + unsigned jc = j + a.size(); + Crossings res; + + //Sweep of the monotonic portions + std::vector<std::vector<unsigned> > cull2 = sweep_bounds(bounds_a[i], bounds_b[j]); + for(unsigned k = 0; k < cull2.size(); k++) { + for(unsigned lx = 0; lx < cull2[k].size(); lx++) { + unsigned l = cull2[k][lx]; + mono_pair(a[i], splits_a[i][k-1], splits_a[i][k], + b[j], splits_b[j][l-1], splits_b[j][l], + res, .1); + } + } + + for(auto & re : res) { re.a = i; re.b = jc; } + + merge_crossings(results[i], res, i); + merge_crossings(results[i], res, jc); + } + } + + return results; +} + +/* This function is similar codewise to the MonoCrosser, the main difference is that it deals with + * only one set of paths and includes self intersection +CrossingSet crossings_among(PathVector const &p) { + CrossingSet results(p.size(), Crossings()); + if(p.empty()) return results; + + std::vector<std::vector<double> > splits = paths_mono_splits(p); + std::vector<std::vector<Rect> > prs = split_bounds(p, splits); + std::vector<Rect> rs; + for(unsigned i = 0; i < prs.size(); i++) rs.push_back(union_list(prs[i])); + + std::vector<std::vector<unsigned> > cull = sweep_bounds(rs); + + //we actually want to do the self-intersections, so add em in: + for(unsigned i = 0; i < cull.size(); i++) cull[i].push_back(i); + + for(unsigned i = 0; i < cull.size(); i++) { + for(unsigned jx = 0; jx < cull[i].size(); jx++) { + unsigned j = cull[i][jx]; + Crossings res; + + //Sweep of the monotonic portions + std::vector<std::vector<unsigned> > cull2 = sweep_bounds(prs[i], prs[j]); + for(unsigned k = 0; k < cull2.size(); k++) { + for(unsigned lx = 0; lx < cull2[k].size(); lx++) { + unsigned l = cull2[k][lx]; + mono_pair(p[i], splits[i][k-1], splits[i][k], + p[j], splits[j][l-1], splits[j][l], + res, .1); + } + } + + for(unsigned k = 0; k < res.size(); k++) { res[k].a = i; res[k].b = j; } + + merge_crossings(results[i], res, i); + merge_crossings(results[j], res, j); + } + } + + return results; +} +*/ + + +Crossings curve_self_crossings(Curve const &a) { + Crossings res; + std::vector<double> spl; + spl.push_back(0); + append(spl, curve_mono_splits(a)); + spl.push_back(1); + for(unsigned i = 1; i < spl.size(); i++) + for(unsigned j = i+1; j < spl.size(); j++) + pair_intersect(a, spl[i-1], spl[i], a, spl[j-1], spl[j], res); + return res; +} + +/* +void mono_curve_intersect(Curve const & A, double Al, double Ah, + Curve const & B, double Bl, double Bh, + Crossings &ret, unsigned depth=0) { + // std::cout << depth << "(" << Al << ", " << Ah << ")\n"; + 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)) { + double tA, tB, c; + if(linear_intersect(A0, A1, B0, B1, tA, tB, c)) { + tA = tA * (Ah - Al) + Al; + tB = tB * (Bh - Bl) + Bl; + 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_curve_intersect(B, Bl, mid, + A, Al, Ah, + ret, depth+1); + mono_curve_intersect(B, mid, Bh, + A, Al, Ah, + ret, depth+1); +} + +std::vector<std::vector<double> > curves_mono_splits(Path const &p) { + std::vector<std::vector<double> > ret; + for(unsigned i = 0; i <= p.size(); i++) { + std::vector<double> spl; + spl.push_back(0); + append(spl, curve_mono_splits(p[i])); + spl.push_back(1); + ret.push_back(spl); + } +} + +std::vector<std::vector<Rect> > curves_split_bounds(Path const &p, std::vector<std::vector<double> > splits) { + std::vector<std::vector<Rect> > ret; + for(unsigned i = 0; i < splits.size(); i++) { + std::vector<Rect> res; + for(unsigned j = 1; j < splits[i].size(); j++) + res.push_back(Rect(p.pointAt(splits[i][j-1]+i), p.pointAt(splits[i][j]+i))); + ret.push_back(res); + } + return ret; +} + +Crossings path_self_crossings(Path const &p) { + Crossings ret; + std::vector<std::vector<unsigned> > cull = sweep_bounds(bounds(p)); + std::vector<std::vector<double> > spl = curves_mono_splits(p); + std::vector<std::vector<Rect> > bnds = curves_split_bounds(p, spl); + for(unsigned i = 0; i < cull.size(); i++) { + Crossings res; + for(unsigned k = 1; k < spl[i].size(); k++) + for(unsigned l = k+1; l < spl[i].size(); l++) + mono_curve_intersect(p[i], spl[i][k-1], spl[i][k], p[i], spl[i][l-1], spl[i][l], res); + offset_crossings(res, i, i); + append(ret, res); + for(unsigned jx = 0; jx < cull[i].size(); jx++) { + unsigned j = cull[i][jx]; + res.clear(); + + std::vector<std::vector<unsigned> > cull2 = sweep_bounds(bnds[i], bnds[j]); + for(unsigned k = 0; k < cull2.size(); k++) { + for(unsigned lx = 0; lx < cull2[k].size(); lx++) { + unsigned l = cull2[k][lx]; + mono_curve_intersect(p[i], spl[i][k-1], spl[i][k], p[j], spl[j][l-1], spl[j][l], res); + } + } + + //if(fabs(int(i)-j) == 1 || fabs(int(i)-j) == p.size()-1) { + Crossings res2; + for(unsigned k = 0; k < res.size(); k++) { + if(res[k].ta != 0 && res[k].ta != 1 && res[k].tb != 0 && res[k].tb != 1) { + res.push_back(res[k]); + } + } + res = res2; + //} + offset_crossings(res, i, j); + append(ret, res); + } + } + return ret; +} +*/ + +Crossings self_crossings(Path const &p) { + Crossings ret; + std::vector<std::vector<unsigned> > cull = sweep_bounds(bounds(p)); + for(unsigned i = 0; i < cull.size(); i++) { + Crossings res = curve_self_crossings(p[i]); + offset_crossings(res, i, i); + append(ret, res); + for(unsigned jx = 0; jx < cull[i].size(); jx++) { + unsigned j = cull[i][jx]; + res.clear(); + pair_intersect(p[i], 0, 1, p[j], 0, 1, res); + + //if(fabs(int(i)-j) == 1 || fabs(int(i)-j) == p.size()-1) { + Crossings res2; + for(auto & re : res) { + if(re.ta != 0 && re.ta != 1 && re.tb != 0 && re.tb != 1) { + res2.push_back(re); + } + } + res = res2; + //} + offset_crossings(res, i, j); + append(ret, res); + } + } + return ret; +} + +void flip_crossings(Crossings &crs) { + for(auto & cr : crs) + cr = Crossing(cr.tb, cr.ta, cr.b, cr.a, !cr.dir); +} + +CrossingSet crossings_among(PathVector const &p) { + CrossingSet results(p.size(), Crossings()); + if(p.empty()) return results; + + SimpleCrosser cc; + + std::vector<std::vector<unsigned> > cull = sweep_bounds(bounds(p)); + for(unsigned i = 0; i < cull.size(); i++) { + Crossings res = self_crossings(p[i]); + for(auto & re : res) { re.a = re.b = i; } + merge_crossings(results[i], res, i); + flip_crossings(res); + merge_crossings(results[i], res, i); + for(unsigned jx = 0; jx < cull[i].size(); jx++) { + unsigned j = cull[i][jx]; + + Crossings res = cc.crossings(p[i], p[j]); + for(auto & re : res) { re.a = i; re.b = j; } + merge_crossings(results[i], res, i); + merge_crossings(results[j], res, j); + } + } + return results; +} + +} + +/* + Local Variables: + mode:c++ + c-file-style:"stroustrup" + c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) + indent-tabs-mode:nil + fill-column:99 + End: +*/ +// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 : |