src/2geom/path-intersection.cpp

	path-intersection.cpp revision 76906e4a110b558e62eab96cd9549c2d8197123b
8001ba81cb851b38d86650a2fef5817facffb763johanengelen#include <2geom/path-intersection.h>
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen
8001ba81cb851b38d86650a2fef5817facffb763johanengelen#include <2geom/ord.h>
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen//for path_direction:
8001ba81cb851b38d86650a2fef5817facffb763johanengelen#include <2geom/sbasis-geometric.h>
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen#include <2geom/line.h>
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen#ifdef HAVE_GSL
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen#include <gsl/gsl_vector.h>
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen#include <gsl/gsl_multiroots.h>
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen#endif
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelennamespace Geom {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould/**
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould * This function computes the winding of the path, given a reference point.
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen * Positive values correspond to counter-clockwise in the mathematical coordinate system,
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen * and clockwise in screen coordinates.  This particular implementation casts a ray in
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen * the positive x direction.  It iterates the path, checking for intersection with the
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen * bounding boxes.  If an intersection is found, the initial/final Y value of the curve is
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen * used to derive a delta on the winding value.  If the point is within the bounding box,
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen * the curve specific winding function is called.
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen */
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelenint winding(Path const &path, Point p) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen  //start on a segment which is not a horizontal line with y = p[y]
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen  Path::const_iterator start;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen  for(Path::const_iterator iter = path.begin(); ; ++iter) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    if(iter == path.end_closed()) { return 0; }
6c3e745a94ef6b25a4ef9f018d350a7535aa45afTed Gould    if(iter->initialPoint()[Y]!=p[Y])   { start = iter; break; }
6c3e745a94ef6b25a4ef9f018d350a7535aa45afTed Gould    if(iter->finalPoint()[Y]!=p[Y])     { start = iter; break; }
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen    if(iter->boundsFast().height()!=0.){ start = iter; break; }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen  }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen  int wind = 0;
e6bdf746e2d9e775704a475a29cc1bb167ec271cjohanengelen  unsigned cnt = 0;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen  bool starting = true;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen  for (Path::const_iterator iter = start; iter != start || starting
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen       ; ++iter, iter = (iter == path.end_closed()) ? path.begin() : iter )
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen  {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    cnt++;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    if(cnt > path.size()) return wind;  //some bug makes this required
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    starting = false;
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen    Rect bounds = (iter->boundsFast());
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    Coord x = p[X], y = p[Y];
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    if(x > bounds.right() || !bounds[Y].contains(y)) continue; //ray doesn't intersect box
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    Point final = iter->finalPoint();
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    Point initial = iter->initialPoint();
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    Cmp final_to_ray = cmp(final[Y], y);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    Cmp initial_to_ray = cmp(initial[Y], y);
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    // if y is included, these will have opposite values, giving order.
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen    Cmp c = cmp(final_to_ray, initial_to_ray);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    if(x < bounds.left()) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        // ray goes through bbox
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        // winding delta determined by position of endpoints
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        if(final_to_ray != EQUAL_TO) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            wind += int(c); // GT = counter-clockwise = 1; LT = clockwise = -1; EQ = not-included = 0
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            //std::cout << int(c) << " ";
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            goto cont;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    } else {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        //inside bbox, use custom per-curve winding thingie
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        int delt = iter->winding(p);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        wind += delt;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        //std::cout << "n" << delt << " ";
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    //Handling the special case of an endpoint on the ray:
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    if(final[Y] == y) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        //Traverse segments until it breaks away from y
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        //99.9% of the time this will happen the first go
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        Path::const_iterator next = iter;
ea8dd7683dd12883474f6cf9b5f424f8ed831166Kris        ++next;
ea8dd7683dd12883474f6cf9b5f424f8ed831166Kris        for(; ; ++next) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            if(next == path.end_closed()) next = path.begin();
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen            Rect bnds = (next->boundsFast());
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            //TODO: X considerations
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            if(bnds.height() > 0) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                //It has diverged
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                if(bnds.contains(p)) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                    const double fudge = 0.01;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                    if(cmp(y, next->valueAt(fudge, Y)) == initial_to_ray) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                        wind += int(c);
c3a8ad9235ff81909bd472707550aef5b91daf7bjohanengelen                        //std::cout << "!!!!!" << int(c) << " ";
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                    }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                    iter = next; // No increment, as the rest of the thing hasn't been counted.
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                } else {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                    Coord ny = next->initialPoint()[Y];
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                    if(cmp(y, ny) == initial_to_ray) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                        //Is a continuation through the ray, so counts windingwise
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                        wind += int(c);
c3a8ad9235ff81909bd472707550aef5b91daf7bjohanengelen                        //std::cout << "!!!!!" << int(c) << " ";
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                    }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                    iter = ++next;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                goto cont;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            if(next==start) return wind;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        //Looks like it looped, which means everything's flat
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        return 0;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    }
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    cont:(void)0;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen  }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen  return wind;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen}
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould/**
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould * This function should only be applied to simple paths (regions), as otherwise
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen * a boolean winding direction is undefined.  It returns true for fill, false for
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen * hole.  Defaults to using the sign of area when it reaches funny cases.
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen */
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelenbool path_direction(Path const &p) {
ed0798e33e14e86b60f8cd254d7131f9d83ba8c0johanengelen    if(p.empty()) return false;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    //could probably be more efficient, but this is a quick job
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    double y = p.initialPoint()[Y];
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    double x = p.initialPoint()[X];
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    Cmp res = cmp(p[0].finalPoint()[Y], y);
76906e4a110b558e62eab96cd9549c2d8197123bKris    /*goto doh;
a9e97816d32eb4bb1a9b34ad9633c7e9749b5c14Johan Engelen    for(unsigned i = 1; i < p.size(); i++) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        Cmp final_to_ray = cmp(p[i].finalPoint()[Y], y);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        Cmp initial_to_ray = cmp(p[i].initialPoint()[Y], y);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        // if y is included, these will have opposite values, giving order.
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        Cmp c = cmp(final_to_ray, initial_to_ray);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        if(c != EQUAL_TO) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            std::vector<double> rs = p[i].roots(y, Y);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            for(unsigned j = 0; j < rs.size(); j++) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                double nx = p[i].valueAt(rs[j], X);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                if(nx > x) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                    x = nx;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                    res = c;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        } else if(final_to_ray == EQUAL_TO) goto doh;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    return res < 0;
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
76906e4a110b558e62eab96cd9549c2d8197123bKris    doh:*/
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        //Otherwise fallback on area
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        Piecewise<D2<SBasis> > pw = p.toPwSb();
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        double area;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        Point centre;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        Geom::centroid(pw, centre, area);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        return area > 0;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen}
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen//pair intersect code based on njh's pair-intersect
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould/** A little sugar for appending a list to another */
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelentemplate<typename T>
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelenvoid append(T &a, T const &b) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    a.insert(a.end(), b.begin(), b.end());
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen}
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould/**
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould * Finds the intersection between the lines defined by A0 & A1, and B0 & B1.
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen * Returns through the last 3 parameters, returning the t-values on the lines
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen * and the cross-product of the deltas (a useful byproduct).  The return value
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen * indicates if the time values are within their proper range on the line segments.
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen */
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelenbool
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelenlinear_intersect(Point A0, Point A1, Point B0, Point B1,
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                 double &tA, double &tB, double &det) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    // kramers rule as cross products
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    Point Ad = A1 - A0,
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen          Bd = B1 - B0,
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen           d = B0 - A0;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    det = cross(Ad, Bd);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    if( 1.0 + det == 1.0 )
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        return false;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    else
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        double detinv = 1.0 / det;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        tA = cross(d, Bd) * detinv;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        tB = cross(d, Ad) * detinv;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        return tA >= 0. && tA <= 1. && tB >= 0. && tB <= 1.;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen}
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen
c3a8ad9235ff81909bd472707550aef5b91daf7bjohanengelen#if 0
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelentypedef union dbl_64{
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen    long long i64;
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen    double d64;
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen};
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelenstatic double EpsilonOf(double value)
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen{
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen    dbl_64 s;
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen    s.d64 = value;
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen    if(s.i64 == 0)
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen    {
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen        s.i64++;
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen        return s.d64 - value;
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen    }
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen    else if(s.i64-- < 0)
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen        return s.d64 - value;
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen    else
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen        return value - s.d64;
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen}
c3a8ad9235ff81909bd472707550aef5b91daf7bjohanengelen#endif
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen#ifdef HAVE_GSL
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelenstruct rparams {
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen    Curve const &A;
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen    Curve const &B;
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen};
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelenstatic int
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelenintersect_polish_f (const gsl_vector * x, void *params,
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen                    gsl_vector * f)
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen{
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen    const double x0 = gsl_vector_get (x, 0);
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen    const double x1 = gsl_vector_get (x, 1);
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen    Geom::Point dx = ((struct rparams *) params)->A(x0) -
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen        ((struct rparams *) params)->B(x1);
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen    gsl_vector_set (f, 0, dx[0]);
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen    gsl_vector_set (f, 1, dx[1]);
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen    return GSL_SUCCESS;
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen}
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen#endif
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelenstatic void
d6519bf53baba32bd74436ad9c85f1fa2c6b6ae9Krzysztof Kosińskiintersect_polish_root (Curve const &A, double &s, Curve const &B, double &t)
d6519bf53baba32bd74436ad9c85f1fa2c6b6ae9Krzysztof Kosiński{
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen    std::vector<Point> as, bs;
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen    as = A.pointAndDerivatives(s, 2);
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen    bs = B.pointAndDerivatives(t, 2);
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen    Point F = as[0] - bs[0];
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen    double best = dot(F, F);
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen    for(int i = 0; i < 4; i++) {
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        /**
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen           we want to solve
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen           J*(x1 - x0) = f(x0)
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen           |dA(s)[0]  -dB(t)[0]|  (X1 - X0) = A(s) - B(t)
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen           |dA(s)[1]  -dB(t)[1]|
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        **/
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        // We're using the standard transformation matricies, which is numerically rather poor.  Much better to solve the equation using elimination.
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen        Affine jack(as[1][0], as[1][1],
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen                    -bs[1][0], -bs[1][1],
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen                    0, 0);
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        Point soln = (F)*jack.inverse();
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        double ns = s - soln[0];
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        double nt = t - soln[1];
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        if (ns<0) ns=0;
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        else if (ns>1) ns=1;
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        if (nt<0) nt=0;
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        else if (nt>1) nt=1;
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        as = A.pointAndDerivatives(ns, 2);
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        bs = B.pointAndDerivatives(nt, 2);
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        F = as[0] - bs[0];
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        double trial = dot(F, F);
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        if (trial > best*0.1) // we have standards, you know
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen            // At this point we could do a line search
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen            break;
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        best = trial;
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        s = ns;
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        t = nt;
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen    }
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen#ifdef HAVE_GSL
d6519bf53baba32bd74436ad9c85f1fa2c6b6ae9Krzysztof Kosiński    int status;
d6519bf53baba32bd74436ad9c85f1fa2c6b6ae9Krzysztof Kosiński    size_t iter = 0;
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen    if(0) { // the GSL version is more accurate, but taints this with GPL
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        const size_t n = 2;
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        struct rparams p = {A, B};
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        gsl_multiroot_function f = {&intersect_polish_f, n, &p};
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        double x_init[2] = {s, t};
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        gsl_vector *x = gsl_vector_alloc (n);
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        gsl_vector_set (x, 0, x_init[0]);
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        gsl_vector_set (x, 1, x_init[1]);
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        const gsl_multiroot_fsolver_type *T = gsl_multiroot_fsolver_hybrids;
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        gsl_multiroot_fsolver *sol = gsl_multiroot_fsolver_alloc (T, 2);
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        gsl_multiroot_fsolver_set (sol, &f, x);
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        do
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        {
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen            iter++;
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen            status = gsl_multiroot_fsolver_iterate (sol);
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen            if (status)   /* check if solver is stuck */
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen                break;
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen            status =
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen                gsl_multiroot_test_residual (sol->f, 1e-12);
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        }
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        while (status == GSL_CONTINUE && iter < 1000);
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        s = gsl_vector_get (sol->x, 0);
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        t = gsl_vector_get (sol->x, 1);
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        gsl_multiroot_fsolver_free (sol);
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        gsl_vector_free (x);
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen    }
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen#endif
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen}
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould/**
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould * This uses the local bounds functions of curves to generically intersect two.
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen * It passes in the curves, time intervals, and keeps track of depth, while
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen * returning the results through the Crossings parameter.
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen */
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelenvoid pair_intersect(Curve const & A, double Al, double Ah,
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                    Curve const & B, double Bl, double Bh,
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen                    Crossings &ret,  unsigned depth = 0) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen   // std::cout << depth << "(" << Al << ", " << Ah << ")\n";
6c3e745a94ef6b25a4ef9f018d350a7535aa45afTed Gould    OptRect Ar = A.boundsLocal(Interval(Al, Ah));
6c3e745a94ef6b25a4ef9f018d350a7535aa45afTed Gould    if (!Ar) return;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen
6c3e745a94ef6b25a4ef9f018d350a7535aa45afTed Gould    OptRect Br = B.boundsLocal(Interval(Bl, Bh));
6c3e745a94ef6b25a4ef9f018d350a7535aa45afTed Gould    if (!Br) return;
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
6c3e745a94ef6b25a4ef9f018d350a7535aa45afTed Gould    if(! Ar->intersects(*Br)) return;
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    //Checks the general linearity of the function
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen    if((depth > 12)) { // || (A.boundsLocal(Interval(Al, Ah), 1).maxExtent() < 0.1
ed0798e33e14e86b60f8cd254d7131f9d83ba8c0johanengelen                    //&&  B.boundsLocal(Interval(Bl, Bh), 1).maxExtent() < 0.1)) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        double tA, tB, c;
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen        if(linear_intersect(A.pointAt(Al), A.pointAt(Ah),
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen                            B.pointAt(Bl), B.pointAt(Bh),
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                            tA, tB, c)) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            tA = tA * (Ah - Al) + Al;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            tB = tB * (Bh - Bl) + Bl;
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen            intersect_polish_root(A, tA,
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen                                  B, tB);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            if(depth % 2)
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                ret.push_back(Crossing(tB, tA, c < 0));
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            else
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                ret.push_back(Crossing(tA, tB, c > 0));
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            return;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    if(depth > 12) return;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    double mid = (Bl + Bh)/2;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    pair_intersect(B, Bl, mid,
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                    A, Al, Ah,
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                    ret, depth+1);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    pair_intersect(B, mid, Bh,
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                    A, Al, Ah,
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                    ret, depth+1);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen}
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelenCrossings pair_intersect(Curve const & A, Interval const &Ad,
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen                         Curve const & B, Interval const &Bd) {
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen    Crossings ret;
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen    pair_intersect(A, Ad.min(), Ad.max(), B, Bd.min(), Bd.max(), ret);
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen    return ret;
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen}
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould/** A simple wrapper around pair_intersect */
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelenCrossings SimpleCrosser::crossings(Curve const &a, Curve const &b) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    Crossings ret;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    pair_intersect(a, 0, 1, b, 0, 1, ret);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    return ret;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen}
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen//same as below but curves not paths
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelenvoid mono_intersect(Curve const &A, double Al, double Ah,
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen                    Curve const &B, double Bl, double Bh,
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen                    Crossings &ret, double tol = 0.1, unsigned depth = 0) {
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen    if( Al >= Ah || Bl >= Bh) return;
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen    //std::cout << " " << depth << "[" << Al << ", " << Ah << "]" << "[" << Bl << ", " << Bh << "]";
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen    Point A0 = A.pointAt(Al), A1 = A.pointAt(Ah),
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen          B0 = B.pointAt(Bl), B1 = B.pointAt(Bh);
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen    //inline code that this implies? (without rect/interval construction)
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen    Rect Ar = Rect(A0, A1), Br = Rect(B0, B1);
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen    if(!Ar.intersects(Br) || A0 == A1 || B0 == B1) return;
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen    if(depth > 12 || (Ar.maxExtent() < tol && Ar.maxExtent() < tol)) {
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        double tA, tB, c;
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen        if(linear_intersect(A.pointAt(Al), A.pointAt(Ah),
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen                            B.pointAt(Bl), B.pointAt(Bh),
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen                            tA, tB, c)) {
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen            tA = tA * (Ah - Al) + Al;
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen            tB = tB * (Bh - Bl) + Bl;
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen            intersect_polish_root(A, tA,
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen                                  B, tB);
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen            if(depth % 2)
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen                ret.push_back(Crossing(tB, tA, c < 0));
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen            else
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen                ret.push_back(Crossing(tA, tB, c > 0));
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen            return;
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        }
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen    }
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen    if(depth > 12) return;
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen    double mid = (Bl + Bh)/2;
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen    mono_intersect(B, Bl, mid,
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen              A, Al, Ah,
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen              ret, tol, depth+1);
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen    mono_intersect(B, mid, Bh,
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen              A, Al, Ah,
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen              ret, tol, depth+1);
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen}
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelenCrossings mono_intersect(Curve const & A, Interval const &Ad,
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen                         Curve const & B, Interval const &Bd) {
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen    Crossings ret;
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen    mono_intersect(A, Ad.min(), Ad.max(), B, Bd.min(), Bd.max(), ret);
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen    return ret;
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen}
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould/**
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould * Takes two paths and time ranges on them, with the invariant that the
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen * paths are monotonic on the range.  Splits A when the linear intersection
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen * doesn't exist or is inaccurate.  Uses the fact that it is monotonic to
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen * do very fast local bounds.
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen */
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelenvoid mono_pair(Path const &A, double Al, double Ah,
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen               Path const &B, double Bl, double Bh,
8001ba81cb851b38d86650a2fef5817facffb763johanengelen               Crossings &ret, double /*tol*/, unsigned depth = 0) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    if( Al >= Ah || Bl >= Bh) return;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    std::cout << " " << depth << "[" << Al << ", " << Ah << "]" << "[" << Bl << ", " << Bh << "]";
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    Point A0 = A.pointAt(Al), A1 = A.pointAt(Ah),
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen          B0 = B.pointAt(Bl), B1 = B.pointAt(Bh);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    //inline code that this implies? (without rect/interval construction)
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen    Rect Ar = Rect(A0, A1), Br = Rect(B0, B1);
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen    if(!Ar.intersects(Br) || A0 == A1 || B0 == B1) return;
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen    if(depth > 12 || (Ar.maxExtent() < 0.1 && Ar.maxExtent() < 0.1)) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        double tA, tB, c;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        if(linear_intersect(A0, A1, B0, B1,
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                            tA, tB, c)) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            tA = tA * (Ah - Al) + Al;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            tB = tB * (Bh - Bl) + Bl;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            if(depth % 2)
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                ret.push_back(Crossing(tB, tA, c < 0));
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            else
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                ret.push_back(Crossing(tA, tB, c > 0));
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            return;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        }
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen    }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    if(depth > 12) return;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    double mid = (Bl + Bh)/2;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    mono_pair(B, Bl, mid,
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen              A, Al, Ah,
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen              ret, depth+1);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    mono_pair(B, mid, Bh,
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen              A, Al, Ah,
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen              ret, depth+1);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen}
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould/** This returns the times when the x or y derivative is 0 in the curve. */
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelenstd::vector<double> curve_mono_splits(Curve const &d) {
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen    Curve* deriv = d.derivative();
07bda0b13ae048815f53f21ad1edbe3cc1b7e4e8Johan Engelen    std::vector<double> rs = deriv->roots(0, X);
07bda0b13ae048815f53f21ad1edbe3cc1b7e4e8Johan Engelen    append(rs, deriv->roots(0, Y));
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen    delete deriv;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    std::sort(rs.begin(), rs.end());
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    return rs;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen}
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould/** Convenience function to add a value to each entry in a vector of doubles. */
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelenstd::vector<double> offset_doubles(std::vector<double> const &x, double offs) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    std::vector<double> ret;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    for(unsigned i = 0; i < x.size(); i++) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        ret.push_back(x[i] + offs);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    return ret;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen}
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould/**
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould * Finds all the monotonic splits for a path.  Only includes the split between
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen * curves if they switch derivative directions at that point.
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen */
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelenstd::vector<double> path_mono_splits(Path const &p) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    std::vector<double> ret;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    if(p.empty()) return ret;
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    bool pdx=2, pdy=2;  //Previous derivative direction
a9e97816d32eb4bb1a9b34ad9633c7e9749b5c14Johan Engelen    for(unsigned i = 0; i < p.size(); i++) {
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen        std::vector<double> spl = offset_doubles(curve_mono_splits(p[i]), i);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        bool dx = p[i].initialPoint()[X] > (spl.empty()? p[i].finalPoint()[X] :
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                                                         p.valueAt(spl.front(), X));
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        bool dy = p[i].initialPoint()[Y] > (spl.empty()? p[i].finalPoint()[Y] :
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                                                         p.valueAt(spl.front(), Y));
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        //The direction changed, include the split time
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        if(dx != pdx || dy != pdy) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            ret.push_back(i);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            pdx = dx; pdy = dy;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        append(ret, spl);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    return ret;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen}
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould/**
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen * Applies path_mono_splits to multiple paths, and returns the results such that
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen * time-set i corresponds to Path i.
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen */
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelenstd::vector<std::vector<double> > paths_mono_splits(std::vector<Path> const &ps) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    std::vector<std::vector<double> > ret;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    for(unsigned i = 0; i < ps.size(); i++)
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        ret.push_back(path_mono_splits(ps[i]));
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    return ret;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen}
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould/**
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould * Processes the bounds for a list of paths and a list of splits on them, yielding a list of rects for each.
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen * Each entry i corresponds to path i of the input.  The number of rects in each entry is guaranteed to be the
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen * number of splits for that path, subtracted by one.
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen */
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelenstd::vector<std::vector<Rect> > split_bounds(std::vector<Path> const &p, std::vector<std::vector<double> > splits) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    std::vector<std::vector<Rect> > ret;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    for(unsigned i = 0; i < p.size(); i++) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        std::vector<Rect> res;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        for(unsigned j = 1; j < splits[i].size(); j++)
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            res.push_back(Rect(p[i].pointAt(splits[i][j-1]), p[i].pointAt(splits[i][j])));
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        ret.push_back(res);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    return ret;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen}
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould/**
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould * This is the main routine of "MonoCrosser", and implements a monotonic strategy on multiple curves.
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen * Finds crossings between two sets of paths, yielding a CrossingSet.  [0, a.size()) of the return correspond
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen * to the sorted crossings of a with paths of b.  The rest of the return, [a.size(), a.size() + b.size()],
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen * corresponds to the sorted crossings of b with paths of a.
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen *
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen * This function does two sweeps, one on the bounds of each path, and after that cull, one on the curves within.
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen * This leads to a certain amount of code complexity, however, most of that is factored into the above functions
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen */
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelenCrossingSet MonoCrosser::crossings(std::vector<Path> const &a, std::vector<Path> const &b) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    if(b.empty()) return CrossingSet(a.size(), Crossings());
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    CrossingSet results(a.size() + b.size(), Crossings());
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    if(a.empty()) return results;
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    std::vector<std::vector<double> > splits_a = paths_mono_splits(a), splits_b = paths_mono_splits(b);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    std::vector<std::vector<Rect> > bounds_a = split_bounds(a, splits_a), bounds_b = split_bounds(b, splits_b);
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen    std::vector<Rect> bounds_a_union, bounds_b_union;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    for(unsigned i = 0; i < bounds_a.size(); i++) bounds_a_union.push_back(union_list(bounds_a[i]));
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    for(unsigned i = 0; i < bounds_b.size(); i++) bounds_b_union.push_back(union_list(bounds_b[i]));
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    std::vector<std::vector<unsigned> > cull = sweep_bounds(bounds_a_union, bounds_b_union);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    Crossings n;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    for(unsigned i = 0; i < cull.size(); i++) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        for(unsigned jx = 0; jx < cull[i].size(); jx++) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            unsigned j = cull[i][jx];
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            unsigned jc = j + a.size();
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            Crossings res;
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            //Sweep of the monotonic portions
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            std::vector<std::vector<unsigned> > cull2 = sweep_bounds(bounds_a[i], bounds_b[j]);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            for(unsigned k = 0; k < cull2.size(); k++) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                for(unsigned lx = 0; lx < cull2[k].size(); lx++) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                    unsigned l = cull2[k][lx];
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                    mono_pair(a[i], splits_a[i][k-1], splits_a[i][k],
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                              b[j], splits_b[j][l-1], splits_b[j][l],
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                              res, .1);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            }
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            for(unsigned k = 0; k < res.size(); k++) { res[k].a = i; res[k].b = jc; }
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            merge_crossings(results[i], res, i);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            merge_crossings(results[i], res, jc);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    return results;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen}
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen/* This function is similar codewise to the MonoCrosser, the main difference is that it deals with
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen * only one set of paths and includes self intersection
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelenCrossingSet crossings_among(std::vector<Path> const &p) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    CrossingSet results(p.size(), Crossings());
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    if(p.empty()) return results;
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    std::vector<std::vector<double> > splits = paths_mono_splits(p);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    std::vector<std::vector<Rect> > prs = split_bounds(p, splits);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    std::vector<Rect> rs;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    for(unsigned i = 0; i < prs.size(); i++) rs.push_back(union_list(prs[i]));
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    std::vector<std::vector<unsigned> > cull = sweep_bounds(rs);
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    //we actually want to do the self-intersections, so add em in:
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    for(unsigned i = 0; i < cull.size(); i++) cull[i].push_back(i);
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    for(unsigned i = 0; i < cull.size(); i++) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        for(unsigned jx = 0; jx < cull[i].size(); jx++) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            unsigned j = cull[i][jx];
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            Crossings res;
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            //Sweep of the monotonic portions
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            std::vector<std::vector<unsigned> > cull2 = sweep_bounds(prs[i], prs[j]);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            for(unsigned k = 0; k < cull2.size(); k++) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                for(unsigned lx = 0; lx < cull2[k].size(); lx++) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                    unsigned l = cull2[k][lx];
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                    mono_pair(p[i], splits[i][k-1], splits[i][k],
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                              p[j], splits[j][l-1], splits[j][l],
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                              res, .1);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            }
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            for(unsigned k = 0; k < res.size(); k++) { res[k].a = i; res[k].b = j; }
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            merge_crossings(results[i], res, i);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            merge_crossings(results[j], res, j);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    }
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    return results;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen}
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen*/
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen
ed0798e33e14e86b60f8cd254d7131f9d83ba8c0johanengelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelenCrossings curve_self_crossings(Curve const &a) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    Crossings res;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    std::vector<double> spl;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    spl.push_back(0);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    append(spl, curve_mono_splits(a));
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    spl.push_back(1);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    for(unsigned i = 1; i < spl.size(); i++)
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        for(unsigned j = i+1; j < spl.size(); j++)
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            pair_intersect(a, spl[i-1], spl[i], a, spl[j-1], spl[j], res);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    return res;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen}
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen/*
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelenvoid mono_curve_intersect(Curve const & A, double Al, double Ah,
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                          Curve const & B, double Bl, double Bh,
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                          Crossings &ret,  unsigned depth=0) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen   // std::cout << depth << "(" << Al << ", " << Ah << ")\n";
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    Point A0 = A.pointAt(Al), A1 = A.pointAt(Ah),
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen          B0 = B.pointAt(Bl), B1 = B.pointAt(Bh);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    //inline code that this implies? (without rect/interval construction)
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    if(!Rect(A0, A1).intersects(Rect(B0, B1)) || A0 == A1 || B0 == B1) return;
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    //Checks the general linearity of the function
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen    if((depth > 12) || (A.boundsLocal(Interval(Al, Ah), 1).maxExtent() < 0.1
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                    &&  B.boundsLocal(Interval(Bl, Bh), 1).maxExtent() < 0.1)) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        double tA, tB, c;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        if(linear_intersect(A0, A1, B0, B1, tA, tB, c)) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            tA = tA * (Ah - Al) + Al;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            tB = tB * (Bh - Bl) + Bl;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            if(depth % 2)
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                ret.push_back(Crossing(tB, tA, c < 0));
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            else
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                ret.push_back(Crossing(tA, tB, c > 0));
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            return;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    if(depth > 12) return;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    double mid = (Bl + Bh)/2;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    mono_curve_intersect(B, Bl, mid,
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                         A, Al, Ah,
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                         ret, depth+1);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    mono_curve_intersect(B, mid, Bh,
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                         A, Al, Ah,
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                         ret, depth+1);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen}
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelenstd::vector<std::vector<double> > curves_mono_splits(Path const &p) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    std::vector<std::vector<double> > ret;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    for(unsigned i = 0; i <= p.size(); i++) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        std::vector<double> spl;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        spl.push_back(0);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        append(spl, curve_mono_splits(p[i]));
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        spl.push_back(1);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        ret.push_back(spl);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen}
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelenstd::vector<std::vector<Rect> > curves_split_bounds(Path const &p, std::vector<std::vector<double> > splits) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    std::vector<std::vector<Rect> > ret;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    for(unsigned i = 0; i < splits.size(); i++) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        std::vector<Rect> res;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        for(unsigned j = 1; j < splits[i].size(); j++)
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            res.push_back(Rect(p.pointAt(splits[i][j-1]+i), p.pointAt(splits[i][j]+i)));
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        ret.push_back(res);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    return ret;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen}
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelenCrossings path_self_crossings(Path const &p) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    Crossings ret;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    std::vector<std::vector<unsigned> > cull = sweep_bounds(bounds(p));
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    std::vector<std::vector<double> > spl = curves_mono_splits(p);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    std::vector<std::vector<Rect> > bnds = curves_split_bounds(p, spl);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    for(unsigned i = 0; i < cull.size(); i++) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        Crossings res;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        for(unsigned k = 1; k < spl[i].size(); k++)
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            for(unsigned l = k+1; l < spl[i].size(); l++)
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                mono_curve_intersect(p[i], spl[i][k-1], spl[i][k], p[i], spl[i][l-1], spl[i][l], res);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        offset_crossings(res, i, i);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        append(ret, res);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        for(unsigned jx = 0; jx < cull[i].size(); jx++) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            unsigned j = cull[i][jx];
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            res.clear();
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            std::vector<std::vector<unsigned> > cull2 = sweep_bounds(bnds[i], bnds[j]);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            for(unsigned k = 0; k < cull2.size(); k++) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                for(unsigned lx = 0; lx < cull2[k].size(); lx++) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                    unsigned l = cull2[k][lx];
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                    mono_curve_intersect(p[i], spl[i][k-1], spl[i][k], p[j], spl[j][l-1], spl[j][l], res);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            }
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            //if(fabs(int(i)-j) == 1 || fabs(int(i)-j) == p.size()-1) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                Crossings res2;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                for(unsigned k = 0; k < res.size(); k++) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                    if(res[k].ta != 0 && res[k].ta != 1 && res[k].tb != 0 && res[k].tb != 1) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                        res.push_back(res[k]);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                    }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                res = res2;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            //}
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            offset_crossings(res, i, j);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            append(ret, res);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    return ret;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen}
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen*/
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen
ed0798e33e14e86b60f8cd254d7131f9d83ba8c0johanengelenCrossings self_crossings(Path const &p) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    Crossings ret;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    std::vector<std::vector<unsigned> > cull = sweep_bounds(bounds(p));
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    for(unsigned i = 0; i < cull.size(); i++) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        Crossings res = curve_self_crossings(p[i]);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        offset_crossings(res, i, i);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        append(ret, res);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        for(unsigned jx = 0; jx < cull[i].size(); jx++) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            unsigned j = cull[i][jx];
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            res.clear();
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            pair_intersect(p[i], 0, 1, p[j], 0, 1, res);
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            //if(fabs(int(i)-j) == 1 || fabs(int(i)-j) == p.size()-1) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                Crossings res2;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                for(unsigned k = 0; k < res.size(); k++) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                    if(res[k].ta != 0 && res[k].ta != 1 && res[k].tb != 0 && res[k].tb != 1) {
ed0798e33e14e86b60f8cd254d7131f9d83ba8c0johanengelen                        res2.push_back(res[k]);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                    }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen                res = res2;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            //}
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            offset_crossings(res, i, j);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            append(ret, res);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    return ret;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen}
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen
ed0798e33e14e86b60f8cd254d7131f9d83ba8c0johanengelenvoid flip_crossings(Crossings &crs) {
ed0798e33e14e86b60f8cd254d7131f9d83ba8c0johanengelen    for(unsigned i = 0; i < crs.size(); i++)
ed0798e33e14e86b60f8cd254d7131f9d83ba8c0johanengelen        crs[i] = Crossing(crs[i].tb, crs[i].ta, crs[i].b, crs[i].a, !crs[i].dir);
ed0798e33e14e86b60f8cd254d7131f9d83ba8c0johanengelen}
ed0798e33e14e86b60f8cd254d7131f9d83ba8c0johanengelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelenCrossingSet crossings_among(std::vector<Path> const &p) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    CrossingSet results(p.size(), Crossings());
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    if(p.empty()) return results;
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    SimpleCrosser cc;
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    std::vector<std::vector<unsigned> > cull = sweep_bounds(bounds(p));
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    for(unsigned i = 0; i < cull.size(); i++) {
ed0798e33e14e86b60f8cd254d7131f9d83ba8c0johanengelen        Crossings res = self_crossings(p[i]);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        for(unsigned k = 0; k < res.size(); k++) { res[k].a = res[k].b = i; }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        merge_crossings(results[i], res, i);
ed0798e33e14e86b60f8cd254d7131f9d83ba8c0johanengelen        flip_crossings(res);
ed0798e33e14e86b60f8cd254d7131f9d83ba8c0johanengelen        merge_crossings(results[i], res, i);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        for(unsigned jx = 0; jx < cull[i].size(); jx++) {
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            unsigned j = cull[i][jx];
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            Crossings res = cc.crossings(p[i], p[j]);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            for(unsigned k = 0; k < res.size(); k++) { res[k].a = i; res[k].b = j; }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            merge_crossings(results[i], res, i);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen            merge_crossings(results[j], res, j);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        }
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen    }
ed0798e33e14e86b60f8cd254d7131f9d83ba8c0johanengelen    return results;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen}
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen}
63267518b4ce196caab66ef8cbdcfc0921206b3djohanengelen
63267518b4ce196caab66ef8cbdcfc0921206b3djohanengelen/*
63267518b4ce196caab66ef8cbdcfc0921206b3djohanengelen  Local Variables:
63267518b4ce196caab66ef8cbdcfc0921206b3djohanengelen  mode:c++
63267518b4ce196caab66ef8cbdcfc0921206b3djohanengelen  c-file-style:"stroustrup"
63267518b4ce196caab66ef8cbdcfc0921206b3djohanengelen  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
63267518b4ce196caab66ef8cbdcfc0921206b3djohanengelen  indent-tabs-mode:nil
63267518b4ce196caab66ef8cbdcfc0921206b3djohanengelen  fill-column:99
63267518b4ce196caab66ef8cbdcfc0921206b3djohanengelen  End:
63267518b4ce196caab66ef8cbdcfc0921206b3djohanengelen*/
a4030d5ca449e7e384bc699cd249ee704faaeab0Chris Morgan// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :