summaryrefslogtreecommitdiffstats
path: root/src/2geom/path-intersection.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'src/2geom/path-intersection.cpp')
-rw-r--r--src/2geom/path-intersection.cpp178
1 files changed, 138 insertions, 40 deletions
diff --git a/src/2geom/path-intersection.cpp b/src/2geom/path-intersection.cpp
index bfb1fb96c..b2d5ceabb 100644
--- a/src/2geom/path-intersection.cpp
+++ b/src/2geom/path-intersection.cpp
@@ -4,8 +4,11 @@
//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 {
@@ -199,6 +202,7 @@ static double EpsilonOf(double value)
}
#endif
+#ifdef HAVE_GSL
struct rparams {
Curve const &A;
Curve const &B;
@@ -219,45 +223,91 @@ intersect_polish_f (const gsl_vector * x, void *params,
return GSL_SUCCESS;
}
+#endif
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]|
+ **/
+
+ // 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;
+ 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};
- 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);
+ 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]);
+ 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);
+ 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);
+ do
+ {
+ iter++;
+ status = gsl_multiroot_fsolver_iterate (sol);
- if (status) /* check if solver is stuck */
- break;
+ if (status) /* check if solver is stuck */
+ break;
- status =
- gsl_multiroot_test_residual (sol->f, 1e-12);
- }
- while (status == GSL_CONTINUE && iter < 1000);
+ 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);
+ s = gsl_vector_get (sol->x, 0);
+ t = gsl_vector_get (sol->x, 1);
- gsl_multiroot_fsolver_free (sol);
- gsl_vector_free (x);
+ gsl_multiroot_fsolver_free (sol);
+ gsl_vector_free (x);
+ }
+#endif
}
/**
@@ -267,7 +317,7 @@ intersect_polish_root (Curve const &A, double &s,
*/
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;
@@ -305,6 +355,13 @@ void pair_intersect(Curve const & A, double Al, double Ah,
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;
@@ -312,6 +369,53 @@ Crossings SimpleCrosser::crossings(Curve const &a, Curve const &b) {
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
@@ -327,11 +431,10 @@ void mono_pair(Path const &A, double Al, double Ah,
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)) {
@@ -343,7 +446,7 @@ void mono_pair(Path const &A, double Al, double Ah,
ret.push_back(Crossing(tA, tB, c > 0));
return;
}
- //}
+ }
if(depth > 12) return;
double mid = (Bl + Bh)/2;
mono_pair(B, Bl, mid,
@@ -356,8 +459,10 @@ void mono_pair(Path const &A, double Al, double Ah,
/** 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;
}
@@ -378,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;
+ 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] :
generated by cgit v1.2.3 (git 2.25.1) at 2026-05-19 02:12:26 +0000