![[tokkee]](http://tokkee.org/images/avatar.png){kind=avatar}
index 5a4e7602056e1234948c82690f558d1d1e17a874..2e4eba519c77c18b1e35a4b18d4657f017a9434c 100644 (file)
//for path_direction:
#include <2geom/sbasis-geometric.h>
+#include <2geom/line.h>
+#ifdef HAVE_GSL
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multiroots.h>
+#endif
namespace Geom {
starting = false;
Rect bounds = *(iter->boundsFast());
Coord x = p[X], y = p[Y];
-
+
if(x > bounds.right() || !bounds[Y].contains(y)) continue; //ray doesn't intersect box
-
+
Point final = iter->finalPoint();
Point initial = iter->initialPoint();
Cmp final_to_ray = cmp(final[Y], y);
Cmp initial_to_ray = cmp(initial[Y], y);
-
+
// if y is included, these will have opposite values, giving order.
- Cmp c = cmp(final_to_ray, initial_to_ray);
+ Cmp c = cmp(final_to_ray, initial_to_ray);
if(x < bounds.left()) {
// ray goes through bbox
// winding delta determined by position of endpoints
const double fudge = 0.01;
if(cmp(y, next->valueAt(fudge, Y)) == initial_to_ray) {
wind += int(c);
- std::cout << "!!!!!" << int(c) << " ";
+ //std::cout << "!!!!!" << int(c) << " ";
}
iter = next; // No increment, as the rest of the thing hasn't been counted.
} else {
if(cmp(y, ny) == initial_to_ray) {
//Is a continuation through the ray, so counts windingwise
wind += int(c);
- std::cout << "!!!!!" << int(c) << " ";
+ //std::cout << "!!!!!" << int(c) << " ";
}
iter = ++next;
}
//Looks like it looped, which means everything's flat
return 0;
}
-
+
cont:(void)0;
}
return wind;
double x = p.initialPoint()[X];
Cmp res = cmp(p[0].finalPoint()[Y], y);
goto doh;
- for(unsigned i = 1; i <= p.size(); i++) {
+ for(unsigned i = 1; i < p.size(); i++) {
Cmp final_to_ray = cmp(p[i].finalPoint()[Y], y);
Cmp initial_to_ray = cmp(p[i].initialPoint()[Y], y);
// if y is included, these will have opposite values, giving order.
} else if(final_to_ray == EQUAL_TO) goto doh;
}
return res < 0;
-
+
doh:
//Otherwise fallback on area
-
+
Piecewise<D2<SBasis> > pw = p.toPwSb();
double area;
Point centre;
}
+#if 0
typedef union dbl_64{
long long i64;
double d64;
else
return value - s.d64;
}
+#endif
+#ifdef HAVE_GSL
struct rparams {
Curve const &A;
Curve const &B;
{
const double x0 = gsl_vector_get (x, 0);
const double x1 = gsl_vector_get (x, 1);
-
- Geom::Point dx = ((struct rparams *) params)->A(x0) -
+
+ Geom::Point dx = ((struct rparams *) params)->A(x0) -
((struct rparams *) params)->B(x1);
-
+
gsl_vector_set (f, 0, dx[0]);
gsl_vector_set (f, 1, dx[1]);
-
+
return GSL_SUCCESS;
}
+#endif
-static void
+static void
intersect_polish_root (Curve const &A, double &s,
Curve const &B, double &t) {
- int status;
- size_t iter = 0;
-
- const size_t n = 2;
- struct rparams p = {A, B};
- gsl_multiroot_function f = {&intersect_polish_f, n, &p};
-
- double x_init[2] = {s, t};
- gsl_vector *x = gsl_vector_alloc (n);
-
- gsl_vector_set (x, 0, x_init[0]);
- gsl_vector_set (x, 1, x_init[1]);
-
- const gsl_multiroot_fsolver_type *T = gsl_multiroot_fsolver_hybrids;
- gsl_multiroot_fsolver *sol = gsl_multiroot_fsolver_alloc (T, 2);
- gsl_multiroot_fsolver_set (sol, &f, x);
-
- do
- {
- iter++;
- status = gsl_multiroot_fsolver_iterate (sol);
-
- if (status) /* check if solver is stuck */
+ std::vector<Point> as, bs;
+ as = A.pointAndDerivatives(s, 2);
+ bs = B.pointAndDerivatives(t, 2);
+ Point F = as[0] - bs[0];
+ double best = dot(F, F);
+
+ for(int i = 0; i < 4; i++) {
+
+ /**
+ we want to solve
+ J*(x1 - x0) = f(x0)
+
+ |dA(s)[0] -dB(t)[0]| (X1 - X0) = A(s) - B(t)
+ |dA(s)[1] -dB(t)[1]|
+ **/
+
+ // We're using the standard transformation matricies, which is numerically rather poor. Much better to solve the equation using elimination.
+
+ Matrix jack(as[1][0], as[1][1],
+ -bs[1][0], -bs[1][1],
+ 0, 0);
+ Point soln = (F)*jack.inverse();
+ double ns = s - soln[0];
+ double nt = t - soln[1];
+
+ if (ns<0) ns=0;
+ else if (ns>1) ns=1;
+ if (nt<0) nt=0;
+ else if (nt>1) nt=1;
+
+ as = A.pointAndDerivatives(ns, 2);
+ bs = B.pointAndDerivatives(nt, 2);
+ F = as[0] - bs[0];
+ double trial = dot(F, F);
+ if (trial > best*0.1) // we have standards, you know
+ // At this point we could do a line search
break;
-
- status =
- gsl_multiroot_test_residual (sol->f, 1e-12);
+ best = trial;
+ s = ns;
+ t = nt;
+ }
+
+#ifdef HAVE_GSL
+ if(0) { // the GSL version is more accurate, but taints this with GPL
+ const size_t n = 2;
+ struct rparams p = {A, B};
+ gsl_multiroot_function f = {&intersect_polish_f, n, &p};
+
+ double x_init[2] = {s, t};
+ gsl_vector *x = gsl_vector_alloc (n);
+
+ gsl_vector_set (x, 0, x_init[0]);
+ gsl_vector_set (x, 1, x_init[1]);
+
+ const gsl_multiroot_fsolver_type *T = gsl_multiroot_fsolver_hybrids;
+ gsl_multiroot_fsolver *sol = gsl_multiroot_fsolver_alloc (T, 2);
+ gsl_multiroot_fsolver_set (sol, &f, x);
+
+ int status = 0;
+ size_t iter = 0;
+ do
+ {
+ iter++;
+ status = gsl_multiroot_fsolver_iterate (sol);
+
+ if (status) /* check if solver is stuck */
+ break;
+
+ status =
+ gsl_multiroot_test_residual (sol->f, 1e-12);
+ }
+ while (status == GSL_CONTINUE && iter < 1000);
+
+ s = gsl_vector_get (sol->x, 0);
+ t = gsl_vector_get (sol->x, 1);
+
+ gsl_multiroot_fsolver_free (sol);
+ gsl_vector_free (x);
}
- while (status == GSL_CONTINUE && iter < 1000);
-
- s = gsl_vector_get (sol->x, 0);
- t = gsl_vector_get (sol->x, 1);
-
- gsl_multiroot_fsolver_free (sol);
- gsl_vector_free (x);
+#endif
}
/**
* It passes in the curves, time intervals, and keeps track of depth, while
* returning the results through the Crossings parameter.
*/
-void pair_intersect(Curve const & A, double Al, double Ah,
+void pair_intersect(Curve const & A, double Al, double Ah,
Curve const & B, double Bl, double Bh,
- Crossings &ret, unsigned depth=0) {
+ Crossings &ret, unsigned depth = 0) {
// std::cout << depth << "(" << Al << ", " << Ah << ")\n";
OptRect Ar = A.boundsLocal(Interval(Al, Ah));
if (!Ar) return;
OptRect Br = B.boundsLocal(Interval(Bl, Bh));
if (!Br) return;
-
+
if(! Ar->intersects(*Br)) return;
-
+
//Checks the general linearity of the function
- if((depth > 12)) { // || (A.boundsLocal(Interval(Al, Ah), 1).maxExtent() < 0.1
+ if((depth > 12)) { // || (A.boundsLocal(Interval(Al, Ah), 1).maxExtent() < 0.1
//&& B.boundsLocal(Interval(Bl, Bh), 1).maxExtent() < 0.1)) {
double tA, tB, c;
- if(linear_intersect(A.pointAt(Al), A.pointAt(Ah),
- B.pointAt(Bl), B.pointAt(Bh),
+ if(linear_intersect(A.pointAt(Al), A.pointAt(Ah),
+ B.pointAt(Bl), B.pointAt(Bh),
tA, tB, c)) {
tA = tA * (Ah - Al) + Al;
tB = tB * (Bh - Bl) + Bl;
ret, depth+1);
}
+Crossings pair_intersect(Curve const & A, Interval const &Ad,
+ Curve const & B, Interval const &Bd) {
+ Crossings ret;
+ pair_intersect(A, Ad.min(), Ad.max(), B, Bd.min(), Bd.max(), ret);
+ return ret;
+}
+
/** A simple wrapper around pair_intersect */
Crossings SimpleCrosser::crossings(Curve const &a, Curve const &b) {
Crossings ret;
return ret;
}
+
+//same as below but curves not paths
+void mono_intersect(Curve const &A, double Al, double Ah,
+ Curve const &B, double Bl, double Bh,
+ Crossings &ret, double tol = 0.1, unsigned depth = 0) {
+ if( Al >= Ah || Bl >= Bh) return;
+ //std::cout << " " << depth << "[" << Al << ", " << Ah << "]" << "[" << Bl << ", " << Bh << "]";
+
+ Point A0 = A.pointAt(Al), A1 = A.pointAt(Ah),
+ B0 = B.pointAt(Bl), B1 = B.pointAt(Bh);
+ //inline code that this implies? (without rect/interval construction)
+ Rect Ar = Rect(A0, A1), Br = Rect(B0, B1);
+ if(!Ar.intersects(Br) || A0 == A1 || B0 == B1) return;
+
+ if(depth > 12 || (Ar.maxExtent() < tol && Ar.maxExtent() < tol)) {
+ double tA, tB, c;
+ if(linear_intersect(A.pointAt(Al), A.pointAt(Ah),
+ B.pointAt(Bl), B.pointAt(Bh),
+ tA, tB, c)) {
+ tA = tA * (Ah - Al) + Al;
+ tB = tB * (Bh - Bl) + Bl;
+ intersect_polish_root(A, tA,
+ B, tB);
+ if(depth % 2)
+ ret.push_back(Crossing(tB, tA, c < 0));
+ else
+ ret.push_back(Crossing(tA, tB, c > 0));
+ return;
+ }
+ }
+ if(depth > 12) return;
+ double mid = (Bl + Bh)/2;
+ mono_intersect(B, Bl, mid,
+ A, Al, Ah,
+ ret, tol, depth+1);
+ mono_intersect(B, mid, Bh,
+ A, Al, Ah,
+ ret, tol, depth+1);
+}
+
+Crossings mono_intersect(Curve const & A, Interval const &Ad,
+ Curve const & B, Interval const &Bd) {
+ Crossings ret;
+ mono_intersect(A, Ad.min(), Ad.max(), B, Bd.min(), Bd.max(), ret);
+ return ret;
+}
+
/**
* Takes two paths and time ranges on them, with the invariant that the
* paths are monotonic on the range. Splits A when the linear intersection
Point A0 = A.pointAt(Al), A1 = A.pointAt(Ah),
B0 = B.pointAt(Bl), B1 = B.pointAt(Bh);
//inline code that this implies? (without rect/interval construction)
- if(!Rect(A0, A1).intersects(Rect(B0, B1)) || A0 == A1 || B0 == B1) return;
-
- //Checks the general linearity of the function
- //if((depth > 12) || (A.boundsLocal(Interval(Al, Ah), 1).maxExtent() < 0.1
- // && B.boundsLocal(Interval(Bl, Bh), 1).maxExtent() < 0.1)) {
+ Rect Ar = Rect(A0, A1), Br = Rect(B0, B1);
+ if(!Ar.intersects(Br) || A0 == A1 || B0 == B1) return;
+
+ if(depth > 12 || (Ar.maxExtent() < 0.1 && Ar.maxExtent() < 0.1)) {
double tA, tB, c;
if(linear_intersect(A0, A1, B0, B1,
tA, tB, c)) {
ret.push_back(Crossing(tA, tB, c > 0));
return;
}
- //}
+ }
if(depth > 12) return;
double mid = (Bl + Bh)/2;
mono_pair(B, Bl, mid,
/** This returns the times when the x or y derivative is 0 in the curve. */
std::vector<double> curve_mono_splits(Curve const &d) {
+ Curve* deriv = d.derivative();
std::vector<double> rs = d.roots(0, X);
append(rs, d.roots(0, Y));
+ delete deriv;
std::sort(rs.begin(), rs.end());
return rs;
}
@@ -376,17 +483,10 @@ std::vector<double> offset_doubles(std::vector<double> const &x, double offs) {
std::vector<double> path_mono_splits(Path const &p) {
std::vector<double> ret;
if(p.empty()) return ret;
- ret.push_back(0);
-
- Curve* deriv = p[0].derivative();
- append(ret, curve_mono_splits(*deriv));
- delete deriv;
-
+
bool pdx=2, pdy=2; //Previous derivative direction
- for(unsigned i = 0; i <= p.size(); i++) {
- deriv = p[i].derivative();
- std::vector<double> spl = offset_doubles(curve_mono_splits(*deriv), i);
- delete deriv;
+ for(unsigned i = 0; i < p.size(); i++) {
+ std::vector<double> spl = offset_doubles(curve_mono_splits(p[i]), i);
bool dx = p[i].initialPoint()[X] > (spl.empty()? p[i].finalPoint()[X] :
p.valueAt(spl.front(), X));
bool dy = p[i].initialPoint()[Y] > (spl.empty()? p[i].finalPoint()[Y] :
}
/**
- * Applies path_mono_splits to multiple paths, and returns the results such that
+ * Applies path_mono_splits to multiple paths, and returns the results such that
* time-set i corresponds to Path i.
*/
std::vector<std::vector<double> > paths_mono_splits(std::vector<Path> const &ps) {
@@ -441,14 +541,14 @@ CrossingSet MonoCrosser::crossings(std::vector<Path> const &a, std::vector<Path>
if(b.empty()) return CrossingSet(a.size(), Crossings());
CrossingSet results(a.size() + b.size(), Crossings());
if(a.empty()) return results;
-
+
std::vector<std::vector<double> > splits_a = paths_mono_splits(a), splits_b = paths_mono_splits(b);
std::vector<std::vector<Rect> > bounds_a = split_bounds(a, splits_a), bounds_b = split_bounds(b, splits_b);
-
- std::vector<Rect> bounds_a_union, bounds_b_union;
+
+ std::vector<Rect> bounds_a_union, bounds_b_union;
for(unsigned i = 0; i < bounds_a.size(); i++) bounds_a_union.push_back(union_list(bounds_a[i]));
for(unsigned i = 0; i < bounds_b.size(); i++) bounds_b_union.push_back(union_list(bounds_b[i]));
-
+
std::vector<std::vector<unsigned> > cull = sweep_bounds(bounds_a_union, bounds_b_union);
Crossings n;
for(unsigned i = 0; i < cull.size(); i++) {
@@ -456,7 +556,7 @@ CrossingSet MonoCrosser::crossings(std::vector<Path> const &a, std::vector<Path>
unsigned j = cull[i][jx];
unsigned jc = j + a.size();
Crossings res;
-
+
//Sweep of the monotonic portions
std::vector<std::vector<unsigned> > cull2 = sweep_bounds(bounds_a[i], bounds_b[j]);
for(unsigned k = 0; k < cull2.size(); k++) {
@@ -467,9 +567,9 @@ CrossingSet MonoCrosser::crossings(std::vector<Path> const &a, std::vector<Path>
res, .1);
}
}
-
+
for(unsigned k = 0; k < res.size(); k++) { res[k].a = i; res[k].b = jc; }
-
+
merge_crossings(results[i], res, i);
merge_crossings(results[i], res, jc);
}
@@ -483,22 +583,22 @@ CrossingSet MonoCrosser::crossings(std::vector<Path> const &a, std::vector<Path>
CrossingSet crossings_among(std::vector<Path> const &p) {
CrossingSet results(p.size(), Crossings());
if(p.empty()) return results;
-
+
std::vector<std::vector<double> > splits = paths_mono_splits(p);
std::vector<std::vector<Rect> > prs = split_bounds(p, splits);
std::vector<Rect> rs;
for(unsigned i = 0; i < prs.size(); i++) rs.push_back(union_list(prs[i]));
-
+
std::vector<std::vector<unsigned> > cull = sweep_bounds(rs);
-
+
//we actually want to do the self-intersections, so add em in:
for(unsigned i = 0; i < cull.size(); i++) cull[i].push_back(i);
-
+
for(unsigned i = 0; i < cull.size(); i++) {
for(unsigned jx = 0; jx < cull[i].size(); jx++) {
unsigned j = cull[i][jx];
Crossings res;
-
+
//Sweep of the monotonic portions
std::vector<std::vector<unsigned> > cull2 = sweep_bounds(prs[i], prs[j]);
for(unsigned k = 0; k < cull2.size(); k++) {
res, .1);
}
}
-
+
for(unsigned k = 0; k < res.size(); k++) { res[k].a = i; res[k].b = j; }
-
+
merge_crossings(results[i], res, i);
merge_crossings(results[j], res, j);
}
}
-
+
return results;
}
*/
}
/*
-void mono_curve_intersect(Curve const & A, double Al, double Ah,
+void mono_curve_intersect(Curve const & A, double Al, double Ah,
Curve const & B, double Bl, double Bh,
Crossings &ret, unsigned depth=0) {
// std::cout << depth << "(" << Al << ", " << Ah << ")\n";
B0 = B.pointAt(Bl), B1 = B.pointAt(Bh);
//inline code that this implies? (without rect/interval construction)
if(!Rect(A0, A1).intersects(Rect(B0, B1)) || A0 == A1 || B0 == B1) return;
-
+
//Checks the general linearity of the function
- if((depth > 12) || (A.boundsLocal(Interval(Al, Ah), 1).maxExtent() < 0.1
+ if((depth > 12) || (A.boundsLocal(Interval(Al, Ah), 1).maxExtent() < 0.1
&& B.boundsLocal(Interval(Bl, Bh), 1).maxExtent() < 0.1)) {
double tA, tB, c;
if(linear_intersect(A0, A1, B0, B1, tA, tB, c)) {
for(unsigned jx = 0; jx < cull[i].size(); jx++) {
unsigned j = cull[i][jx];
res.clear();
-
+
std::vector<std::vector<unsigned> > cull2 = sweep_bounds(bnds[i], bnds[j]);
for(unsigned k = 0; k < cull2.size(); k++) {
for(unsigned lx = 0; lx < cull2[k].size(); lx++) {
mono_curve_intersect(p[i], spl[i][k-1], spl[i][k], p[j], spl[j][l-1], spl[j][l], res);
}
}
-
+
//if(fabs(int(i)-j) == 1 || fabs(int(i)-j) == p.size()-1) {
Crossings res2;
for(unsigned k = 0; k < res.size(); k++) {
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++) {
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]);
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);