#include "path-intersection.h" #include "ord.h" //for path_direction: #include "sbasis-geometric.h" namespace Geom { /* This function computes the winding of the path, given a reference point. * Positive values correspond to counter-clockwise in the mathematical coordinate system, * and clockwise in screen coordinates. This particular implementation casts a ray in * the positive x direction. It iterates the path, checking for intersection with the * bounding boxes. If an intersection is found, the initial/final Y value of the curve is * used to derive a delta on the winding value. If the point is within the bounding box, * the curve specific winding function is called. */ int winding(Path const &path, Point p) { //start on a segment which is not a horizontal line with y = p[y] Path::const_iterator start; for(Path::const_iterator iter = path.begin(); ; ++iter) { if(iter == path.end_closed()) { return 0; } if(iter->initialPoint()[Y]!=p[Y]) { start = iter; break; } if(iter->finalPoint()[Y]!=p[Y]) { start = iter; break; } if(iter->boundsFast().height()!=0.){ start = iter; break; } } int wind = 0; unsigned cnt = 0; bool starting = true; for (Path::const_iterator iter = start; iter != start || starting ; ++iter, iter = (iter == path.end_closed()) ? path.begin() : iter ) { cnt++; if(cnt > path.size()) return wind; //some bug makes this required starting = false; Rect bounds = iter->boundsFast(); Coord x = p[X], y = p[Y]; if(x > bounds.right() || !bounds[Y].contains(y)) continue; //ray doesn't intersect box Point final = iter->finalPoint(); Point initial = iter->initialPoint(); Cmp final_to_ray = cmp(final[Y], y); Cmp initial_to_ray = cmp(initial[Y], y); // if y is included, these will have opposite values, giving order. Cmp c = cmp(final_to_ray, initial_to_ray); if(x < bounds.left()) { // ray goes through bbox // winding delta determined by position of endpoints if(final_to_ray != EQUAL_TO) { wind += int(c); // GT = counter-clockwise = 1; LT = clockwise = -1; EQ = not-included = 0 //std::cout << int(c) << " "; goto cont; } } else { //inside bbox, use custom per-curve winding thingie int delt = iter->winding(p); wind += delt; //std::cout << "n" << delt << " "; } //Handling the special case of an endpoint on the ray: if(final[Y] == y) { //Traverse segments until it breaks away from y //99.9% of the time this will happen the first go Path::const_iterator next = iter; next++; for(; ; next++) { if(next == path.end_closed()) next = path.begin(); Rect bnds = next->boundsFast(); //TODO: X considerations if(bnds.height() > 0) { //It has diverged if(bnds.contains(p)) { const double fudge = 0.01; if(cmp(y, next->valueAt(fudge, Y)) == initial_to_ray) { wind += int(c); // std::cout << "!!!!!" << int(c) << " "; } iter = next; // No increment, as the rest of the thing hasn't been counted. } else { Coord ny = next->initialPoint()[Y]; if(cmp(y, ny) == initial_to_ray) { //Is a continuation through the ray, so counts windingwise wind += int(c); // std::cout << "!!!!!" << int(c) << " "; } iter = ++next; } goto cont; } if(next==start) return wind; } //Looks like it looped, which means everything's flat return 0; } cont:(void)0; } return wind; } /* 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; //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); goto doh; 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 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 > 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 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 A0, Point A1, Point B0, Point B1, double &tA, double &tB, double &det) { // kramers rule as cross products Point Ad = A1 - A0, Bd = B1 - B0, d = B0 - A0; det = cross(Ad, Bd); if( 1.0 + det == 1.0 ) return false; else { double detinv = 1.0 / det; tA = cross(d, Bd) * detinv; tB = cross(d, Ad) * detinv; return tA >= 0. && tA <= 1. && tB >= 0. && tB <= 1.; } } /* 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"; Rect Ar = A.boundsLocal(Interval(Al, Ah)); if(Ar.isEmpty()) return; Rect Br = B.boundsLocal(Interval(Bl, Bh)); if(Br.isEmpty()) 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; 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); } // 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; } /* 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) 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_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 curve_mono_splits(Curve const &d) { std::vector rs = d.roots(0, X); append(rs, d.roots(0, Y)); std::sort(rs.begin(), rs.end()); return rs; } // Convenience function to add a value to each entry in a vector of doubles. std::vector offset_doubles(std::vector const &x, double offs) { std::vector ret; for(unsigned i = 0; i < x.size(); i++) { ret.push_back(x[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 path_mono_splits(Path const &p) { std::vector ret; if(p.empty()) return ret; ret.push_back(0); Curve* deriv = p[0].derivative(); append(ret, curve_mono_splits(*deriv)); delete deriv; int pdx=2, pdy=2; //Previous derivative direction for(unsigned i = 0; i <= p.size(); i++) { deriv = p[i].derivative(); std::vector spl = offset_doubles(curve_mono_splits(*deriv), i); delete deriv; int dx = p[i].initialPoint()[X] > (spl.empty()? p[i].finalPoint()[X] : p.valueAt(spl.front(), X)); int dy = p[i].initialPoint()[Y] > (spl.empty()? p[i].finalPoint()[Y] : p.valueAt(spl.front(), Y)); //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 > paths_mono_splits(std::vector const &ps) { std::vector > ret; for(unsigned i = 0; i < ps.size(); i++) ret.push_back(path_mono_splits(ps[i])); 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 > split_bounds(std::vector const &p, std::vector > splits) { std::vector > ret; for(unsigned i = 0; i < p.size(); i++) { std::vector res; for(unsigned j = 1; j < splits[i].size(); j++) res.push_back(Rect(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(std::vector const &a, std::vector const &b) { if(b.empty()) return CrossingSet(a.size(), Crossings()); CrossingSet results(a.size() + b.size(), Crossings()); if(a.empty()) return results; std::vector > splits_a = paths_mono_splits(a), splits_b = paths_mono_splits(b); std::vector > bounds_a = split_bounds(a, splits_a), bounds_b = split_bounds(b, splits_b); std::vector bounds_a_union, bounds_b_union; for(unsigned i = 0; i < bounds_a.size(); i++) bounds_a_union.push_back(union_list(bounds_a[i])); for(unsigned i = 0; i < bounds_b.size(); i++) bounds_b_union.push_back(union_list(bounds_b[i])); std::vector > 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 > 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(unsigned k = 0; k < res.size(); k++) { res[k].a = i; res[k].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(std::vector const &p) { CrossingSet results(p.size(), Crossings()); if(p.empty()) return results; std::vector > splits = paths_mono_splits(p); std::vector > prs = split_bounds(p, splits); std::vector rs; for(unsigned i = 0; i < prs.size(); i++) rs.push_back(union_list(prs[i])); std::vector > 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 > 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 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 > curves_mono_splits(Path const &p) { std::vector > ret; for(unsigned i = 0; i <= p.size(); i++) { std::vector spl; spl.push_back(0); append(spl, curve_mono_splits(p[i])); spl.push_back(1); ret.push_back(spl); } } std::vector > curves_split_bounds(Path const &p, std::vector > splits) { std::vector > ret; for(unsigned i = 0; i < splits.size(); i++) { std::vector 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 > cull = sweep_bounds(bounds(p)); std::vector > spl = curves_mono_splits(p); std::vector > 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 > 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 > 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(unsigned k = 0; k < res.size(); k++) { if(res[k].ta != 0 && res[k].ta != 1 && res[k].tb != 0 && res[k].tb != 1) { res2.push_back(res[k]); } } res = res2; //} offset_crossings(res, i, j); append(ret, res); } } return ret; } void flip_crossings(Crossings &crs) { for(unsigned i = 0; i < crs.size(); i++) crs[i] = Crossing(crs[i].tb, crs[i].ta, crs[i].b, crs[i].a, !crs[i].dir); } CrossingSet crossings_among(std::vector const &p) { CrossingSet results(p.size(), Crossings()); if(p.empty()) return results; SimpleCrosser cc; std::vector > cull = sweep_bounds(bounds(p)); for(unsigned i = 0; i < cull.size(); i++) { Crossings res = self_crossings(p[i]); for(unsigned k = 0; k < res.size(); k++) { res[k].a = res[k].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(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; } }