dbb33756bd786e9432e18ec7be93f8c416e1b492Jon A. Cruz#include <2geom/path-intersection.h>

dbb33756bd786e9432e18ec7be93f8c416e1b492Jon A. Cruz

dbb33756bd786e9432e18ec7be93f8c416e1b492Jon A. Cruz#include <2geom/ord.h>

e9b6af083e34e2397a8ddbe9781920733d09d151Ted Gould

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm#include <2geom/sbasis-geometric.h>

23d859f2ce09c04ed802cb4912cc9c50f512f0a2bgk#include <2geom/line.h>

23d859f2ce09c04ed802cb4912cc9c50f512f0a2bgk#ifdef HAVE_GSL

9dc68827cbd515262ecb8d5ae8547d9e82c72e00Jon A. Cruz#include <gsl/gsl_vector.h>

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm#include <gsl/gsl_multiroots.h>

e9b6af083e34e2397a8ddbe9781920733d09d151Ted Gould#endif

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrmnamespace Geom {

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm

035e109c2dce9f6a9552f75d09b1573311d02546tweenkint winding(Path const &path, Point p) {

035e109c2dce9f6a9552f75d09b1573311d02546tweenk //start on a segment which is not a horizontal line with y = p[y]

035e109c2dce9f6a9552f75d09b1573311d02546tweenk Path::const_iterator start;

035e109c2dce9f6a9552f75d09b1573311d02546tweenk for(Path::const_iterator iter = path.begin(); ; ++iter) {

47badd0035ae8c9135c51444f3770b17ae504ddaAlex Valavanis if(iter == path.end_closed()) { return 0; }

035e109c2dce9f6a9552f75d09b1573311d02546tweenk if(iter->initialPoint()[Y]!=p[Y]) { start = iter; break; }

035e109c2dce9f6a9552f75d09b1573311d02546tweenk if(iter->finalPoint()[Y]!=p[Y]) { start = iter; break; }

035e109c2dce9f6a9552f75d09b1573311d02546tweenk if(iter->boundsFast()->height()!=0.){ start = iter; break; }

035e109c2dce9f6a9552f75d09b1573311d02546tweenk }

035e109c2dce9f6a9552f75d09b1573311d02546tweenk int wind = 0;

035e109c2dce9f6a9552f75d09b1573311d02546tweenk unsigned cnt = 0;

edcba9f1706559e93aa06a7173daa6cd6516acb5Johan B. C. Engelen bool starting = true;

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm for (Path::const_iterator iter = start; iter != start || starting

23d859f2ce09c04ed802cb4912cc9c50f512f0a2bgk ; ++iter, iter = (iter == path.end_closed()) ? path.begin() : iter )

035e109c2dce9f6a9552f75d09b1573311d02546tweenk {

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm cnt++;

ae22ad7adc4a7a418e71f5dbab8c1f0f7f464562johanengelen if(cnt > path.size()) return wind; //some bug makes this required

edcba9f1706559e93aa06a7173daa6cd6516acb5Johan B. C. Engelen starting = false;

edcba9f1706559e93aa06a7173daa6cd6516acb5Johan B. C. Engelen Rect bounds = *(iter->boundsFast());

035e109c2dce9f6a9552f75d09b1573311d02546tweenk Coord x = p[X], y = p[Y];

035e109c2dce9f6a9552f75d09b1573311d02546tweenk

035e109c2dce9f6a9552f75d09b1573311d02546tweenk if(x > bounds.right() || !bounds[Y].contains(y)) continue; //ray doesn't intersect box

035e109c2dce9f6a9552f75d09b1573311d02546tweenk

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm Point final = iter->finalPoint();

0d00bc9f32167e81375a4be524572b27e2894ee4John Smith Point initial = iter->initialPoint();

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm Cmp final_to_ray = cmp(final[Y], y);

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm Cmp initial_to_ray = cmp(initial[Y], y);

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm // if y is included, these will have opposite values, giving order.

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm Cmp c = cmp(final_to_ray, initial_to_ray);

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm if(x < bounds.left()) {

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm // ray goes through bbox

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm // winding delta determined by position of endpoints

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm if(final_to_ray != EQUAL_TO) {

738092bcf0d040b2431137e191dfd7cf3ce3afadJohan Engelen wind += int(c); // GT = counter-clockwise = 1; LT = clockwise = -1; EQ = not-included = 0

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm //std::cout << int(c) << " ";

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm goto cont;

738092bcf0d040b2431137e191dfd7cf3ce3afadJohan Engelen }

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm } else {

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm //inside bbox, use custom per-curve winding thingie

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm int delt = iter->winding(p);

a95be1234ba4df33d6d074589edaa56f0d546069buliabyak wind += delt;

a95be1234ba4df33d6d074589edaa56f0d546069buliabyak //std::cout << "n" << delt << " ";

a95be1234ba4df33d6d074589edaa56f0d546069buliabyak }

a95be1234ba4df33d6d074589edaa56f0d546069buliabyak //Handling the special case of an endpoint on the ray:

a95be1234ba4df33d6d074589edaa56f0d546069buliabyak if(final[Y] == y) {

a95be1234ba4df33d6d074589edaa56f0d546069buliabyak //Traverse segments until it breaks away from y

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm //99.9% of the time this will happen the first go

0d00bc9f32167e81375a4be524572b27e2894ee4John Smith Path::const_iterator next = iter;

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm next++;

0d00bc9f32167e81375a4be524572b27e2894ee4John Smith for(; ; next++) {

0d00bc9f32167e81375a4be524572b27e2894ee4John Smith if(next == path.end_closed()) next = path.begin();

0d00bc9f32167e81375a4be524572b27e2894ee4John Smith Rect bnds = *(next->boundsFast());

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm //TODO: X considerations

963f23115db07f460bdd862b957f8bd9dba88b9bgustav_b if(bnds.height() > 0) {

852b4f6c7a572bc2ccbd96e80c4063a38f77153bjohanengelen //It has diverged

e9b6af083e34e2397a8ddbe9781920733d09d151Ted Gould if(bnds.contains(p)) {

d67024fee7c560721468212767fb6aa221d5b542John Smith const double fudge = 0.01;

f0ed14f45951d21de3ff2c7c131f6aafa6e30c17buliabyak if(cmp(y, next->valueAt(fudge, Y)) == initial_to_ray) {

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm wind += int(c);

0d00bc9f32167e81375a4be524572b27e2894ee4John Smith //std::cout << "!!!!!" << int(c) << " ";

0d00bc9f32167e81375a4be524572b27e2894ee4John Smith }

c67c19fc5d1f6d97cc795b5a53998434b431c641John Smith iter = next; // No increment, as the rest of the thing hasn't been counted.

c67c19fc5d1f6d97cc795b5a53998434b431c641John Smith } else {

c67c19fc5d1f6d97cc795b5a53998434b431c641John Smith Coord ny = next->initialPoint()[Y];

c67c19fc5d1f6d97cc795b5a53998434b431c641John Smith if(cmp(y, ny) == initial_to_ray) {

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen //Is a continuation through the ray, so counts windingwise

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen wind += int(c);

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm //std::cout << "!!!!!" << int(c) << " ";

963f23115db07f460bdd862b957f8bd9dba88b9bgustav_b }

963f23115db07f460bdd862b957f8bd9dba88b9bgustav_b iter = ++next;

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm }

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen goto cont;

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen }

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen if(next==start) return wind;

0d00bc9f32167e81375a4be524572b27e2894ee4John Smith }

0d00bc9f32167e81375a4be524572b27e2894ee4John Smith //Looks like it looped, which means everything's flat

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen return 0;

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm }

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm cont:(void)0;

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen }

0d00bc9f32167e81375a4be524572b27e2894ee4John Smith return wind;

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm}

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm

02db7fad736ff0812658ef9c5f82ac2e64ffdefcAlex Valavanisbool path_direction(Path const &p) {

02db7fad736ff0812658ef9c5f82ac2e64ffdefcAlex Valavanis if(p.empty()) return false;

02db7fad736ff0812658ef9c5f82ac2e64ffdefcAlex Valavanis //could probably be more efficient, but this is a quick job

c67c19fc5d1f6d97cc795b5a53998434b431c641John Smith double y = p.initialPoint()[Y];

c67c19fc5d1f6d97cc795b5a53998434b431c641John Smith double x = p.initialPoint()[X];

0d00bc9f32167e81375a4be524572b27e2894ee4John Smith Cmp res = cmp(p[0].finalPoint()[Y], y);

260ccdd5a6e2a22e13ce13a71bd292da1be3e1adAlex Valavanis goto doh;

02db7fad736ff0812658ef9c5f82ac2e64ffdefcAlex Valavanis for(unsigned i = 1; i < p.size(); i++) {

02db7fad736ff0812658ef9c5f82ac2e64ffdefcAlex Valavanis Cmp final_to_ray = cmp(p[i].finalPoint()[Y], y);

02db7fad736ff0812658ef9c5f82ac2e64ffdefcAlex Valavanis Cmp initial_to_ray = cmp(p[i].initialPoint()[Y], y);

02db7fad736ff0812658ef9c5f82ac2e64ffdefcAlex Valavanis // if y is included, these will have opposite values, giving order.

02db7fad736ff0812658ef9c5f82ac2e64ffdefcAlex Valavanis Cmp c = cmp(final_to_ray, initial_to_ray);

02db7fad736ff0812658ef9c5f82ac2e64ffdefcAlex Valavanis if(c != EQUAL_TO) {

c67c19fc5d1f6d97cc795b5a53998434b431c641John Smith std::vector<double> rs = p[i].roots(y, Y);

c67c19fc5d1f6d97cc795b5a53998434b431c641John Smith for(unsigned j = 0; j < rs.size(); j++) {

0d00bc9f32167e81375a4be524572b27e2894ee4John Smith double nx = p[i].valueAt(rs[j], X);

260ccdd5a6e2a22e13ce13a71bd292da1be3e1adAlex Valavanis if(nx > x) {

02db7fad736ff0812658ef9c5f82ac2e64ffdefcAlex Valavanis x = nx;

02db7fad736ff0812658ef9c5f82ac2e64ffdefcAlex Valavanis res = c;

02db7fad736ff0812658ef9c5f82ac2e64ffdefcAlex Valavanis }

02db7fad736ff0812658ef9c5f82ac2e64ffdefcAlex Valavanis }

02db7fad736ff0812658ef9c5f82ac2e64ffdefcAlex Valavanis } else if(final_to_ray == EQUAL_TO) goto doh;

02db7fad736ff0812658ef9c5f82ac2e64ffdefcAlex Valavanis }

c67c19fc5d1f6d97cc795b5a53998434b431c641John Smith return res < 0;

c67c19fc5d1f6d97cc795b5a53998434b431c641John Smith

0d00bc9f32167e81375a4be524572b27e2894ee4John Smith doh:

260ccdd5a6e2a22e13ce13a71bd292da1be3e1adAlex Valavanis //Otherwise fallback on area

02db7fad736ff0812658ef9c5f82ac2e64ffdefcAlex Valavanis

02db7fad736ff0812658ef9c5f82ac2e64ffdefcAlex Valavanis Piecewise<D2<SBasis> > pw = p.toPwSb();

02db7fad736ff0812658ef9c5f82ac2e64ffdefcAlex Valavanis double area;

02db7fad736ff0812658ef9c5f82ac2e64ffdefcAlex Valavanis Point centre;

02db7fad736ff0812658ef9c5f82ac2e64ffdefcAlex Valavanis Geom::centroid(pw, centre, area);

02db7fad736ff0812658ef9c5f82ac2e64ffdefcAlex Valavanis return area > 0;

c67c19fc5d1f6d97cc795b5a53998434b431c641John Smith}

0d00bc9f32167e81375a4be524572b27e2894ee4John Smith

0d00bc9f32167e81375a4be524572b27e2894ee4John Smith

8e7d13ca30e4b9671afc47f03dc11affd5507077Alex Valavanistemplate<typename T>

8e7d13ca30e4b9671afc47f03dc11affd5507077Alex Valavanisvoid append(T &a, T const &b) {

8e7d13ca30e4b9671afc47f03dc11affd5507077Alex Valavanis a.insert(a.end(), b.begin(), b.end());

0d00bc9f32167e81375a4be524572b27e2894ee4John Smith}

8e7d13ca30e4b9671afc47f03dc11affd5507077Alex Valavanis

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelenlinear_intersect(Point A0, Point A1, Point B0, Point B1,

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen double &tA, double &tB, double &det) {

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen // kramers rule as cross products

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen Point Ad = A1 - A0,

e018daf5696bee513e8574a6730e2ce79d2a3a23johanengelen Bd = B1 - B0,

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen d = B0 - A0;

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen det = cross(Ad, Bd);

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen if( 1.0 + det == 1.0 )

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen return false;

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen else

b3722066cc416b6582d9a9a3ffc8ba2e868e102cjohanengelen {

c42ca006cc2593efa9a0ea19f2680ea9d5792654Markus Engel double detinv = 1.0 / det;

11d3fae21261f21488c39b12b716385b5a91c847Johan Engelen tA = cross(d, Bd) * detinv;

b3722066cc416b6582d9a9a3ffc8ba2e868e102cjohanengelen tB = cross(d, Ad) * detinv;

607ec3f31779845d307f157ff34472da27b8bdbcjohanengelen return tA >= 0. && tA <= 1. && tB >= 0. && tB <= 1.;

b3722066cc416b6582d9a9a3ffc8ba2e868e102cjohanengelen }

b3722066cc416b6582d9a9a3ffc8ba2e868e102cjohanengelen}

b3722066cc416b6582d9a9a3ffc8ba2e868e102cjohanengelen

b3722066cc416b6582d9a9a3ffc8ba2e868e102cjohanengelen

b3722066cc416b6582d9a9a3ffc8ba2e868e102cjohanengelen#if 0

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelentypedef union dbl_64{

b3722066cc416b6582d9a9a3ffc8ba2e868e102cjohanengelen long long i64;

b3722066cc416b6582d9a9a3ffc8ba2e868e102cjohanengelen double d64;

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen};

0d00bc9f32167e81375a4be524572b27e2894ee4John Smith

0d00bc9f32167e81375a4be524572b27e2894ee4John Smithstatic double EpsilonOf(double value)

b4998608f5fbde14c744b6ab8020664300e11f80jucablues{

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm dbl_64 s;

0d00bc9f32167e81375a4be524572b27e2894ee4John Smith s.d64 = value;

0d00bc9f32167e81375a4be524572b27e2894ee4John Smith if(s.i64 == 0)

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm {

0d00bc9f32167e81375a4be524572b27e2894ee4John Smith s.i64++;

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm return s.d64 - value;

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen }

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen else if(s.i64-- < 0)

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen return s.d64 - value;

d67024fee7c560721468212767fb6aa221d5b542John Smith else

d67024fee7c560721468212767fb6aa221d5b542John Smith return value - s.d64;

d67024fee7c560721468212767fb6aa221d5b542John Smith}

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm#endif

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm#ifdef HAVE_GSL

852b4f6c7a572bc2ccbd96e80c4063a38f77153bjohanengelenstruct rparams {

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm Curve const &A;

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm Curve const &B;

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm};

c0cd5511d3b975ebe07d019c1f5528108725e438johanengelen

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrmstatic int

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrmintersect_polish_f (const gsl_vector * x, void *params,

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm gsl_vector * f)

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm{

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm const double x0 = gsl_vector_get (x, 0);

a95be1234ba4df33d6d074589edaa56f0d546069buliabyak const double x1 = gsl_vector_get (x, 1);

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm Geom::Point dx = ((struct rparams *) params)->A(x0) -

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm ((struct rparams *) params)->B(x1);

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm

90a3966dd44e306d23febc15ebd65cde07d7a4ddTed Gould gsl_vector_set (f, 0, dx[0]);

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm gsl_vector_set (f, 1, dx[1]);

c5526a2c3001be486990d816757dd5ac028b3c3fjohanengelen

c5526a2c3001be486990d816757dd5ac028b3c3fjohanengelen return GSL_SUCCESS;

c5526a2c3001be486990d816757dd5ac028b3c3fjohanengelen}

c5526a2c3001be486990d816757dd5ac028b3c3fjohanengelen#endif

c5526a2c3001be486990d816757dd5ac028b3c3fjohanengelen

fdad0bc3aa765bdc61fe39e8c4da03f717525dccjohanengelenstatic void

0d00bc9f32167e81375a4be524572b27e2894ee4John Smithintersect_polish_root (Curve const &A, double &s,

0d00bc9f32167e81375a4be524572b27e2894ee4John Smith Curve const &B, double &t) {

5cf332777b4c27336d64c273ac63bce3ee27a53dAlex Valavanis int status;

c5526a2c3001be486990d816757dd5ac028b3c3fjohanengelen size_t iter = 0;

c5526a2c3001be486990d816757dd5ac028b3c3fjohanengelen std::vector<Point> as, bs;

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm as = A.pointAndDerivatives(s, 2);

c5526a2c3001be486990d816757dd5ac028b3c3fjohanengelen bs = B.pointAndDerivatives(t, 2);

c5526a2c3001be486990d816757dd5ac028b3c3fjohanengelen Point F = as[0] - bs[0];

0d00bc9f32167e81375a4be524572b27e2894ee4John Smith double best = dot(F, F);

0d00bc9f32167e81375a4be524572b27e2894ee4John Smith

c5526a2c3001be486990d816757dd5ac028b3c3fjohanengelen for(int i = 0; i < 4; i++) {

c5526a2c3001be486990d816757dd5ac028b3c3fjohanengelen

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm /**

90a3966dd44e306d23febc15ebd65cde07d7a4ddTed Gould

f0ed14f45951d21de3ff2c7c131f6aafa6e30c17buliabyak // We're using the standard transformation matricies, which is numerically rather poor. Much better to solve the equation using elimination.

f0ed14f45951d21de3ff2c7c131f6aafa6e30c17buliabyak

f0ed14f45951d21de3ff2c7c131f6aafa6e30c17buliabyak Matrix jack(as[1][0], as[1][1],

f0ed14f45951d21de3ff2c7c131f6aafa6e30c17buliabyak -bs[1][0], -bs[1][1],

f0ed14f45951d21de3ff2c7c131f6aafa6e30c17buliabyak 0, 0);

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm Point soln = (F)*jack.inverse();

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm double ns = s - soln[0];

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm double nt = t - soln[1];

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm if (ns<0) ns=0;

c0cd5511d3b975ebe07d019c1f5528108725e438johanengelen else if (ns>1) ns=1;

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm if (nt<0) nt=0;

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm else if (nt>1) nt=1;

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm

0d00bc9f32167e81375a4be524572b27e2894ee4John Smith as = A.pointAndDerivatives(ns, 2);

0d00bc9f32167e81375a4be524572b27e2894ee4John Smith bs = B.pointAndDerivatives(nt, 2);

0d00bc9f32167e81375a4be524572b27e2894ee4John Smith F = as[0] - bs[0];

0d00bc9f32167e81375a4be524572b27e2894ee4John Smith double trial = dot(F, F);

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm if (trial > best*0.1) // we have standards, you know

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm // At this point we could do a line search

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm break;

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm best = trial;

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm s = ns;

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm t = nt;

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm }

0d00bc9f32167e81375a4be524572b27e2894ee4John Smith

0d00bc9f32167e81375a4be524572b27e2894ee4John Smith#ifdef HAVE_GSL

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm if(0) { // the GSL version is more accurate, but taints this with GPL

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen const size_t n = 2;

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen struct rparams p = {A, B};

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen gsl_multiroot_function f = {&intersect_polish_f, n, &p};

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm double x_init[2] = {s, t};

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen gsl_vector *x = gsl_vector_alloc (n);

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm gsl_vector_set (x, 0, x_init[0]);

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm gsl_vector_set (x, 1, x_init[1]);

f0ed14f45951d21de3ff2c7c131f6aafa6e30c17buliabyak

f0ed14f45951d21de3ff2c7c131f6aafa6e30c17buliabyak const gsl_multiroot_fsolver_type *T = gsl_multiroot_fsolver_hybrids;

f0ed14f45951d21de3ff2c7c131f6aafa6e30c17buliabyak gsl_multiroot_fsolver *sol = gsl_multiroot_fsolver_alloc (T, 2);

f0ed14f45951d21de3ff2c7c131f6aafa6e30c17buliabyak gsl_multiroot_fsolver_set (sol, &f, x);

c730b7fa118ad88d7cc6047dcc78e7098cdf79caKrzysztof Kosiński

f0ed14f45951d21de3ff2c7c131f6aafa6e30c17buliabyak do

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen {

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen iter++;

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen status = gsl_multiroot_fsolver_iterate (sol);

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm if (status) /* check if solver is stuck */

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm break;

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz status =

852b4f6c7a572bc2ccbd96e80c4063a38f77153bjohanengelen gsl_multiroot_test_residual (sol->f, 1e-12);

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen }

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen while (status == GSL_CONTINUE && iter < 1000);

852b4f6c7a572bc2ccbd96e80c4063a38f77153bjohanengelen

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm s = gsl_vector_get (sol->x, 0);

8261c3f0a43f2bbfa6f25edbae152278c40d9d97Johan B. C. Engelen t = gsl_vector_get (sol->x, 1);

f4ce66b44581cea80a736cf20244c0539a7525daJohan B. C. Engelen

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen gsl_multiroot_fsolver_free (sol);

90a3966dd44e306d23febc15ebd65cde07d7a4ddTed Gould gsl_vector_free (x);

c730b7fa118ad88d7cc6047dcc78e7098cdf79caKrzysztof Kosiński }

f0ed14f45951d21de3ff2c7c131f6aafa6e30c17buliabyak#endif

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm}

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm

0d00bc9f32167e81375a4be524572b27e2894ee4John Smithvoid pair_intersect(Curve const & A, double Al, double Ah,

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm Curve const & B, double Bl, double Bh,

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz Crossings &ret, unsigned depth = 0) {

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz // std::cout << depth << "(" << Al << ", " << Ah << ")\n";

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz OptRect Ar = A.boundsLocal(Interval(Al, Ah));

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz if (!Ar) return;

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz OptRect Br = B.boundsLocal(Interval(Bl, Bh));

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz if (!Br) return;

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz if(! Ar->intersects(*Br)) return;

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz //Checks the general linearity of the function

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz if((depth > 12)) { // || (A.boundsLocal(Interval(Al, Ah), 1).maxExtent() < 0.1

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz //&& B.boundsLocal(Interval(Bl, Bh), 1).maxExtent() < 0.1)) {

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz double tA, tB, c;

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz if(linear_intersect(A.pointAt(Al), A.pointAt(Ah),

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz B.pointAt(Bl), B.pointAt(Bh),

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz tA, tB, c)) {

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz tA = tA * (Ah - Al) + Al;

edcba9f1706559e93aa06a7173daa6cd6516acb5Johan B. C. Engelen tB = tB * (Bh - Bl) + Bl;

edcba9f1706559e93aa06a7173daa6cd6516acb5Johan B. C. Engelen intersect_polish_root(A, tA,

edcba9f1706559e93aa06a7173daa6cd6516acb5Johan B. C. Engelen B, tB);

edcba9f1706559e93aa06a7173daa6cd6516acb5Johan B. C. Engelen if(depth % 2)

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm ret.push_back(Crossing(tB, tA, c < 0));

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm else

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm ret.push_back(Crossing(tA, tB, c > 0));

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm return;

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm }

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm }

0d00bc9f32167e81375a4be524572b27e2894ee4John Smith if(depth > 12) return;

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm double mid = (Bl + Bh)/2;

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm pair_intersect(B, Bl, mid,

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm A, Al, Ah,

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm ret, depth+1);

852b4f6c7a572bc2ccbd96e80c4063a38f77153bjohanengelen pair_intersect(B, mid, Bh,

852b4f6c7a572bc2ccbd96e80c4063a38f77153bjohanengelen A, Al, Ah,

852b4f6c7a572bc2ccbd96e80c4063a38f77153bjohanengelen ret, depth+1);

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen}

852b4f6c7a572bc2ccbd96e80c4063a38f77153bjohanengelen

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelenCrossings pair_intersect(Curve const & A, Interval const &Ad,

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen Curve const & B, Interval const &Bd) {

852b4f6c7a572bc2ccbd96e80c4063a38f77153bjohanengelen Crossings ret;

32143ea5edb30f496040c1f24538e38d9453fc06Johan B. C. Engelen pair_intersect(A, Ad.min(), Ad.max(), B, Bd.min(), Bd.max(), ret);

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen return ret;

84d6d1f7365e49f2936df9df890ce49d2c000ce2Kris}

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen

cd392fe863cb06ba22f2b54bdb96ef3a6ee6bc5bJohan EngelenCrossings SimpleCrosser::crossings(Curve const &a, Curve const &b) {

cd392fe863cb06ba22f2b54bdb96ef3a6ee6bc5bJohan Engelen Crossings ret;

cd392fe863cb06ba22f2b54bdb96ef3a6ee6bc5bJohan Engelen pair_intersect(a, 0, 1, b, 0, 1, ret);

90a3966dd44e306d23febc15ebd65cde07d7a4ddTed Gould return ret;

90a3966dd44e306d23febc15ebd65cde07d7a4ddTed Gould}

90a3966dd44e306d23febc15ebd65cde07d7a4ddTed Gould

90a3966dd44e306d23febc15ebd65cde07d7a4ddTed Gould

90a3966dd44e306d23febc15ebd65cde07d7a4ddTed Gouldvoid mono_intersect(Curve const &A, double Al, double Ah,

90a3966dd44e306d23febc15ebd65cde07d7a4ddTed Gould Curve const &B, double Bl, double Bh,

aea6c63514922bfc46ced140fa2877576aa21203johanengelen Crossings &ret, double tol = 0.1, unsigned depth = 0) {

90a3966dd44e306d23febc15ebd65cde07d7a4ddTed Gould if( Al >= Ah || Bl >= Bh) return;

90a3966dd44e306d23febc15ebd65cde07d7a4ddTed Gould //std::cout << " " << depth << "[" << Al << ", " << Ah << "]" << "[" << Bl << ", " << Bh << "]";

90a3966dd44e306d23febc15ebd65cde07d7a4ddTed Gould

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen Point A0 = A.pointAt(Al), A1 = A.pointAt(Ah),

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen B0 = B.pointAt(Bl), B1 = B.pointAt(Bh);

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen //inline code that this implies? (without rect/interval construction)

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen Rect Ar = Rect(A0, A1), Br = Rect(B0, B1);

852b4f6c7a572bc2ccbd96e80c4063a38f77153bjohanengelen if(!Ar.intersects(Br) || A0 == A1 || B0 == B1) return;

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm if(depth > 12 || (Ar.maxExtent() < tol && Ar.maxExtent() < tol)) {

1e7c268648bcbae15fc13b8c94dee677b401d9b3gustav_b double tA, tB, c;

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm if(linear_intersect(A.pointAt(Al), A.pointAt(Ah),

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm B.pointAt(Bl), B.pointAt(Bh),

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm tA, tB, c)) {

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm tA = tA * (Ah - Al) + Al;

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm tB = tB * (Bh - Bl) + Bl;

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm intersect_polish_root(A, tA,

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm B, tB);

a95be1234ba4df33d6d074589edaa56f0d546069buliabyak if(depth % 2)

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm ret.push_back(Crossing(tB, tA, c < 0));

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm else

738092bcf0d040b2431137e191dfd7cf3ce3afadJohan Engelen ret.push_back(Crossing(tA, tB, c > 0));

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm return;

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm }

769a6887551cf7ff7bce4b48d3ac303cbea69507Liam P. White }

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm if(depth > 12) return;

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm double mid = (Bl + Bh)/2;

a95be1234ba4df33d6d074589edaa56f0d546069buliabyak mono_intersect(B, Bl, mid,

a95be1234ba4df33d6d074589edaa56f0d546069buliabyak A, Al, Ah,

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm ret, tol, depth+1);

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm mono_intersect(B, mid, Bh,

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm A, Al, Ah,

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm ret, tol, depth+1);

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm}

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrmCrossings mono_intersect(Curve const & A, Interval const &Ad,

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm Curve const & B, Interval const &Bd) {

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm Crossings ret;

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm mono_intersect(A, Ad.min(), Ad.max(), B, Bd.min(), Bd.max(), ret);

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm return ret;

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm}

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrmvoid mono_pair(Path const &A, double Al, double Ah,

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm Path const &B, double Bl, double Bh,

edcba9f1706559e93aa06a7173daa6cd6516acb5Johan B. C. Engelen Crossings &ret, double /*tol*/, unsigned depth = 0) {

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz if( Al >= Ah || Bl >= Bh) return;

edcba9f1706559e93aa06a7173daa6cd6516acb5Johan B. C. Engelen std::cout << " " << depth << "[" << Al << ", " << Ah << "]" << "[" << Bl << ", " << Bh << "]";

edcba9f1706559e93aa06a7173daa6cd6516acb5Johan B. C. Engelen

edcba9f1706559e93aa06a7173daa6cd6516acb5Johan B. C. Engelen Point A0 = A.pointAt(Al), A1 = A.pointAt(Ah),

edcba9f1706559e93aa06a7173daa6cd6516acb5Johan B. C. Engelen B0 = B.pointAt(Bl), B1 = B.pointAt(Bh);

edcba9f1706559e93aa06a7173daa6cd6516acb5Johan B. C. Engelen //inline code that this implies? (without rect/interval construction)

edcba9f1706559e93aa06a7173daa6cd6516acb5Johan B. C. Engelen Rect Ar = Rect(A0, A1), Br = Rect(B0, B1);

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm if(!Ar.intersects(Br) || A0 == A1 || B0 == B1) return;

edcba9f1706559e93aa06a7173daa6cd6516acb5Johan B. C. Engelen

ae22ad7adc4a7a418e71f5dbab8c1f0f7f464562johanengelen if(depth > 12 || (Ar.maxExtent() < 0.1 && Ar.maxExtent() < 0.1)) {

edcba9f1706559e93aa06a7173daa6cd6516acb5Johan B. C. Engelen double tA, tB, c;

edcba9f1706559e93aa06a7173daa6cd6516acb5Johan B. C. Engelen if(linear_intersect(A0, A1, B0, B1,

edcba9f1706559e93aa06a7173daa6cd6516acb5Johan B. C. Engelen tA, tB, c)) {

edcba9f1706559e93aa06a7173daa6cd6516acb5Johan B. C. Engelen tA = tA * (Ah - Al) + Al;

3ae3dd4eb0191cc986a035c790b8b97e6c6e4ee4johanengelen tB = tB * (Bh - Bl) + Bl;

edcba9f1706559e93aa06a7173daa6cd6516acb5Johan B. C. Engelen if(depth % 2)

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz ret.push_back(Crossing(tB, tA, c < 0));

e25076a0c8caf01a05b34c96485471798f61b6cbEric Greveson else

edcba9f1706559e93aa06a7173daa6cd6516acb5Johan B. C. Engelen ret.push_back(Crossing(tA, tB, c > 0));

edcba9f1706559e93aa06a7173daa6cd6516acb5Johan B. C. Engelen return;

48c5a523ec28475a6f2b5e1f1a087e44d56f684fbuliabyak }

edcba9f1706559e93aa06a7173daa6cd6516acb5Johan B. C. Engelen }

edcba9f1706559e93aa06a7173daa6cd6516acb5Johan B. C. Engelen if(depth > 12) return;

edcba9f1706559e93aa06a7173daa6cd6516acb5Johan B. C. Engelen double mid = (Bl + Bh)/2;

edcba9f1706559e93aa06a7173daa6cd6516acb5Johan B. C. Engelen mono_pair(B, Bl, mid,

edcba9f1706559e93aa06a7173daa6cd6516acb5Johan B. C. Engelen A, Al, Ah,

edcba9f1706559e93aa06a7173daa6cd6516acb5Johan B. C. Engelen ret, depth+1);

edcba9f1706559e93aa06a7173daa6cd6516acb5Johan B. C. Engelen mono_pair(B, mid, Bh,

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz A, Al, Ah,

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz ret, depth+1);

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz}

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruzstd::vector<double> curve_mono_splits(Curve const &d) {

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz Curve* deriv = d.derivative();

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz std::vector<double> rs = d.roots(0, X);

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz append(rs, d.roots(0, Y));

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz delete deriv;

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz std::sort(rs.begin(), rs.end());

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz return rs;

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz}

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruzstd::vector<double> offset_doubles(std::vector<double> const &x, double offs) {

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz std::vector<double> ret;

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz for(unsigned i = 0; i < x.size(); i++) {

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz ret.push_back(x[i] + offs);

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz }

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz return ret;

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz}

dd0ae0a28fda34d3805f7fe6deece97c4192910aJohan B. C. Engelen

dd0ae0a28fda34d3805f7fe6deece97c4192910aJohan B. C. Engelenstd::vector<double> path_mono_splits(Path const &p) {

dd0ae0a28fda34d3805f7fe6deece97c4192910aJohan B. C. Engelen std::vector<double> ret;

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz if(p.empty()) return ret;

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz bool pdx=2, pdy=2; //Previous derivative direction

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz for(unsigned i = 0; i < p.size(); i++) {

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz std::vector<double> spl = offset_doubles(curve_mono_splits(p[i]), i);

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz bool dx = p[i].initialPoint()[X] > (spl.empty()? p[i].finalPoint()[X] :

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz p.valueAt(spl.front(), X));

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz bool dy = p[i].initialPoint()[Y] > (spl.empty()? p[i].finalPoint()[Y] :

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz p.valueAt(spl.front(), Y));

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz //The direction changed, include the split time

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz if(dx != pdx || dy != pdy) {

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz ret.push_back(i);

edcba9f1706559e93aa06a7173daa6cd6516acb5Johan B. C. Engelen pdx = dx; pdy = dy;

edcba9f1706559e93aa06a7173daa6cd6516acb5Johan B. C. Engelen }

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm append(ret, spl);

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm }

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm return ret;

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm}

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrmstd::vector<std::vector<double> > paths_mono_splits(std::vector<Path> const &ps) {

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz std::vector<std::vector<double> > ret;

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz for(unsigned i = 0; i < ps.size(); i++)

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz ret.push_back(path_mono_splits(ps[i]));

852b4f6c7a572bc2ccbd96e80c4063a38f77153bjohanengelen return ret;

9dc68827cbd515262ecb8d5ae8547d9e82c72e00Jon A. Cruz}

9dc68827cbd515262ecb8d5ae8547d9e82c72e00Jon A. Cruz

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelenstd::vector<std::vector<Rect> > split_bounds(std::vector<Path> const &p, std::vector<std::vector<double> > splits) {

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen std::vector<std::vector<Rect> > ret;

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen for(unsigned i = 0; i < p.size(); i++) {

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen std::vector<Rect> res;

852b4f6c7a572bc2ccbd96e80c4063a38f77153bjohanengelen for(unsigned j = 1; j < splits[i].size(); j++)

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen res.push_back(Rect(p[i].pointAt(splits[i][j-1]), p[i].pointAt(splits[i][j])));

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen ret.push_back(res);

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz }

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz return ret;

32143ea5edb30f496040c1f24538e38d9453fc06Johan B. C. Engelen}

852b4f6c7a572bc2ccbd96e80c4063a38f77153bjohanengelen

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelenCrossingSet MonoCrosser::crossings(std::vector<Path> const &a, std::vector<Path> const &b) {

852b4f6c7a572bc2ccbd96e80c4063a38f77153bjohanengelen if(b.empty()) return CrossingSet(a.size(), Crossings());

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen CrossingSet results(a.size() + b.size(), Crossings());

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen if(a.empty()) return results;

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz

f4d474ff2d58b3d32dacd5feed0c164e8df4936cJon A. Cruz std::vector<std::vector<double> > splits_a = paths_mono_splits(a), splits_b = paths_mono_splits(b);

32143ea5edb30f496040c1f24538e38d9453fc06Johan B. C. Engelen std::vector<std::vector<Rect> > bounds_a = split_bounds(a, splits_a), bounds_b = split_bounds(b, splits_b);

852b4f6c7a572bc2ccbd96e80c4063a38f77153bjohanengelen

9dc68827cbd515262ecb8d5ae8547d9e82c72e00Jon A. Cruz std::vector<Rect> bounds_a_union, bounds_b_union;

9dc68827cbd515262ecb8d5ae8547d9e82c72e00Jon A. Cruz for(unsigned i = 0; i < bounds_a.size(); i++) bounds_a_union.push_back(union_list(bounds_a[i]));

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen for(unsigned i = 0; i < bounds_b.size(); i++) bounds_b_union.push_back(union_list(bounds_b[i]));

32143ea5edb30f496040c1f24538e38d9453fc06Johan B. C. Engelen

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm std::vector<std::vector<unsigned> > cull = sweep_bounds(bounds_a_union, bounds_b_union);

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm Crossings n;

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm for(unsigned i = 0; i < cull.size(); i++) {

fd4b29a5cdef220804dfed85fec8acb5daceec5fpjrm for(unsigned jx = 0; jx < cull[i].size(); jx++) {

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen unsigned j = cull[i][jx];

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen unsigned jc = j + a.size();

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen Crossings res;

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen //Sweep of the monotonic portions

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen std::vector<std::vector<unsigned> > cull2 = sweep_bounds(bounds_a[i], bounds_b[j]);

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen for(unsigned k = 0; k < cull2.size(); k++) {

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen for(unsigned lx = 0; lx < cull2[k].size(); lx++) {

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen unsigned l = cull2[k][lx];

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen mono_pair(a[i], splits_a[i][k-1], splits_a[i][k],

90a3966dd44e306d23febc15ebd65cde07d7a4ddTed Gould b[j], splits_b[j][l-1], splits_b[j][l],

90a3966dd44e306d23febc15ebd65cde07d7a4ddTed Gould res, .1);

f4ce66b44581cea80a736cf20244c0539a7525daJohan B. C. Engelen }

90a3966dd44e306d23febc15ebd65cde07d7a4ddTed Gould }

90a3966dd44e306d23febc15ebd65cde07d7a4ddTed Gould

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen for(unsigned k = 0; k < res.size(); k++) { res[k].a = i; res[k].b = jc; }

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen

0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen merge_crossings(results[i], res, i);

607ec3f31779845d307f157ff34472da27b8bdbcjohanengelen merge_crossings(results[i], res, jc);

607ec3f31779845d307f157ff34472da27b8bdbcjohanengelen }

607ec3f31779845d307f157ff34472da27b8bdbcjohanengelen }

607ec3f31779845d307f157ff34472da27b8bdbcjohanengelen

607ec3f31779845d307f157ff34472da27b8bdbcjohanengelen return results;

607ec3f31779845d307f157ff34472da27b8bdbcjohanengelen}

607ec3f31779845d307f157ff34472da27b8bdbcjohanengelen

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;

}

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(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<Path> 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(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;

}

}