![[tokkee]](http://tokkee.org/images/avatar.png){kind=avatar}
index f883f8c63da48bbb8eb85b1aae0f2fbc23561e45..2e4eba519c77c18b1e35a4b18d4657f017a9434c 100644 (file)
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
//Looks like it looped, which means everything's flat
return 0;
}
-
+
cont:(void)0;
}
return wind;
} 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;
{
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;
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]|
+ |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.
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];
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);
}
* 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) {
// std::cout << depth << "(" << Al << ", " << Ah << ")\n";
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;
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),
+ 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;
std::vector<double> path_mono_splits(Path const &p) {
std::vector<double> ret;
if(p.empty()) return ret;
-
+
bool pdx=2, pdy=2; //Previous derivative direction
for(unsigned i = 0; i < p.size(); i++) {
std::vector<double> spl = offset_doubles(curve_mono_splits(p[i]), i);
}
/**
- * 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) {
@@ -541,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++) {
@@ -556,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++) {
@@ -567,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);
}
@@ -583,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);