summaryrefslogtreecommitdiffstats
path: root/src/2geom/basic-intersection.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'src/2geom/basic-intersection.cpp')
-rw-r--r--src/2geom/basic-intersection.cpp177
1 files changed, 126 insertions, 51 deletions
diff --git a/src/2geom/basic-intersection.cpp b/src/2geom/basic-intersection.cpp
index 48b990445..334c23fea 100644
--- a/src/2geom/basic-intersection.cpp
+++ b/src/2geom/basic-intersection.cpp
@@ -1,8 +1,14 @@
+
+
+
#include <2geom/basic-intersection.h>
+#include <2geom/sbasis-to-bezier.h>
#include <2geom/exception.h>
+
+
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multiroots.h>
-
+
unsigned intersect_steps = 0;
@@ -16,10 +22,10 @@ public:
}
void split(double t, OldBezier &a, OldBezier &b) const;
Point operator()(double t) const;
-
+
~OldBezier() {}
- void bounds(double &minax, double &maxax,
+ void bounds(double &minax, double &maxax,
double &minay, double &maxay) {
// Compute bounding box for a
minax = p[0][X]; // These are the most likely to be extremal
@@ -45,17 +51,17 @@ public:
}
}
-
+
};
static std::vector<std::pair<double, double> >
find_intersections( OldBezier a, OldBezier b);
-static std::vector<std::pair<double, double> >
+static std::vector<std::pair<double, double> >
find_self_intersections(OldBezier const &Sb, D2<SBasis> const & A);
std::vector<std::pair<double, double> >
-find_intersections( vector<Geom::Point> const & A,
+find_intersections( vector<Geom::Point> const & A,
vector<Geom::Point> const & B) {
OldBezier a, b;
a.p = A;
@@ -63,23 +69,23 @@ find_intersections( vector<Geom::Point> const & A,
return find_intersections(a,b);
}
-std::vector<std::pair<double, double> >
+std::vector<std::pair<double, double> >
find_self_intersections(OldBezier const &/*Sb*/) {
THROW_NOTIMPLEMENTED();
}
-std::vector<std::pair<double, double> >
+std::vector<std::pair<double, double> >
find_self_intersections(D2<SBasis> const & A) {
OldBezier Sb;
- Sb.p = sbasis_to_bezier(A);
+ sbasis_to_bezier(Sb.p, A);
return find_self_intersections(Sb, A);
}
-static std::vector<std::pair<double, double> >
+static std::vector<std::pair<double, double> >
find_self_intersections(OldBezier const &Sb, D2<SBasis> const & A) {
-
+
vector<double> dr = roots(derivative(A[X]));
{
vector<double> dyr = roots(derivative(A[Y]));
@@ -90,9 +96,9 @@ find_self_intersections(OldBezier const &Sb, D2<SBasis> const & A) {
// We want to be sure that we have no empty segments
sort(dr.begin(), dr.end());
unique(dr.begin(), dr.end());
-
+
std::vector<std::pair<double, double> > all_si;
-
+
vector<OldBezier> pieces;
{
OldBezier in = Sb, l, r;
@@ -104,7 +110,7 @@ find_self_intersections(OldBezier const &Sb, D2<SBasis> const & A) {
}
for(unsigned i = 0; i < dr.size()-1; i++) {
for(unsigned j = i+1; j < dr.size()-1; j++) {
- std::vector<std::pair<double, double> > section =
+ std::vector<std::pair<double, double> > section =
find_intersections( pieces[i], pieces[j]);
for(unsigned k = 0; k < section.size(); k++) {
double l = section[k].first;
@@ -118,17 +124,17 @@ find_self_intersections(OldBezier const &Sb, D2<SBasis> const & A) {
}
}
}
-
+
// Because i is in order, all_si should be roughly already in order?
//sort(all_si.begin(), all_si.end());
//unique(all_si.begin(), all_si.end());
-
+
return all_si;
}
/* The value of 1.0 / (1L<<14) is enough for most applications */
const double INV_EPS = (1L<<14);
-
+
/*
* split the curve at the midpoint, returning an array with the two parts
* Temporary storage is minimized by using part of the storage for the result
@@ -142,12 +148,12 @@ void OldBezier::split(double t, OldBezier &left, OldBezier &right) const {
std::copy(p.begin(), p.end(), Vtemp[0]);
/* Triangle computation */
- for (unsigned i = 1; i < sz; i++) {
+ for (unsigned i = 1; i < sz; i++) {
for (unsigned j = 0; j < sz - i; j++) {
Vtemp[i][j] = lerp(t, Vtemp[i-1][j], Vtemp[i-1][j+1]);
}
}
-
+
left.p.resize(sz);
right.p.resize(sz);
for (unsigned j = 0; j < sz; j++)
@@ -156,6 +162,7 @@ void OldBezier::split(double t, OldBezier &left, OldBezier &right) const {
right.p[j] = Vtemp[sz-1-j][j];
}
+#if 0
/*
* split the curve at the midpoint, returning an array with the two parts
* Temporary storage is minimized by using part of the storage for the result
@@ -169,13 +176,35 @@ Point OldBezier::operator()(double t) const {
std::copy(p.begin(), p.end(), Vtemp[0]);
/* Triangle computation */
- for (unsigned i = 1; i < sz; i++) {
+ for (unsigned i = 1; i < sz; i++) {
for (unsigned j = 0; j < sz - i; j++) {
Vtemp[i][j] = lerp(t, Vtemp[i-1][j], Vtemp[i-1][j+1]);
}
}
return Vtemp[sz-1][0];
}
+#endif
+
+// suggested by Sederberg.
+Point OldBezier::operator()(double t) const {
+ int n = p.size()-1;
+ double u, bc, tn, tmp;
+ int i;
+ Point r;
+ for(int dim = 0; dim < 2; dim++) {
+ u = 1.0 - t;
+ bc = 1;
+ tn = 1;
+ tmp = p[0][dim]*u;
+ for(i=1; i<n; i++){
+ tn = tn*t;
+ bc = bc*(n-i+1)/i;
+ tmp = (tmp + tn*bc*p[i][dim])*u;
+ }
+ r[dim] = (tmp + tn*t*p[n][dim]);
+ }
+ return r;
+}
/*
@@ -183,11 +212,11 @@ Point OldBezier::operator()(double t) const {
* Several observations:
* First, it is cheaper to compute the bounding box of the second curve
* and test its bounding box for interference than to use a more direct
- * approach of comparing all control points of the second curve with
+ * approach of comparing all control points of the second curve with
* the various edges of the bounding box of the first curve to test
* for interference.
* Second, after a few subdivisions it is highly probable that two corners
- * of the bounding box of a given Bezier curve are the first and last
+ * of the bounding box of a given Bezier curve are the first and last
* control point. Once this happens once, it happens for all subsequent
* subcurves. It might be worth putting in a test and then short-circuit
* code for further subdivision levels.
@@ -209,39 +238,39 @@ bool intersect_BB( OldBezier a, OldBezier b ) {
return not( ( minax > maxbx ) || ( minay > maxby )
|| ( minbx > maxax ) || ( minby > maxay ) );
}
-
-/*
- * Recursively intersect two curves keeping track of their real parameters
+
+/*
+ * Recursively intersect two curves keeping track of their real parameters
* and depths of intersection.
* The results are returned in a 2-D array of doubles indicating the parameters
* for which intersections are found. The parameters are in the order the
* intersections were found, which is probably not in sorted order.
- * When an intersection is found, the parameter value for each of the two
+ * When an intersection is found, the parameter value for each of the two
* is stored in the index elements array, and the index is incremented.
- *
+ *
* If either of the curves has subdivisions left before it is straight
* (depth > 0)
* that curve (possibly both) is (are) subdivided at its (their) midpoint(s).
* the depth(s) is (are) decremented, and the parameter value(s) corresponding
* to the midpoints(s) is (are) computed.
- * Then each of the subcurves of one curve is intersected with each of the
+ * Then each of the subcurves of one curve is intersected with each of the
* subcurves of the other curve, first by testing the bounding boxes for
* interference. If there is any bounding box interference, the corresponding
* subcurves are recursively intersected.
- *
+ *
* If neither curve has subdivisions left, the line segments from the first
- * to last control point of each segment are intersected. (Actually the
+ * to last control point of each segment are intersected. (Actually the
* only the parameter value corresponding to the intersection point is found).
*
* The apriori flatness test is probably more efficient than testing at each
* level of recursion, although a test after three or four levels would
- * probably be worthwhile, since many curves become flat faster than their
+ * probably be worthwhile, since many curves become flat faster than their
* asymptotic rate for the first few levels of recursion.
*
* The bounding box test fails much more frequently than it succeeds, providing
* substantial pruning of the search space.
*
- * Each (sub)curve is subdivided only once, hence it is not possible that for
+ * Each (sub)curve is subdivided only once, hence it is not possible that for
* one final line intersection test the subdivision was at one level, while
* for another final line intersection test the subdivision (of the same curve)
* was at another. Since the line segments share endpoints, the intersection
@@ -335,7 +364,7 @@ void recursively_intersect( OldBezier a, double t0, double t1, int deptha,
}
inline double log4( double x ) { return log(x)/log(4.); }
-
+
/*
* Wang's theorem is used to estimate the level of subdivision required,
* but only if the bounding boxes interfere at the top level.
@@ -354,7 +383,7 @@ unsigned wangs_theorem(OldBezier a) {
double la2 = Lmax( ( a.p[3] - a.p[2] ) - (a.p[2] - a.p[1]) );
double l0 = std::max(la1, la2);
unsigned ra;
- if( l0 * 0.75 * M_SQRT2 + 1.0 == 1.0 )
+ if( l0 * 0.75 * M_SQRT2 + 1.0 == 1.0 )
ra = 0;
else
ra = (unsigned)ceil( log4( M_SQRT2 * 6.0 / 8.0 * INV_EPS * l0 ) );
@@ -367,63 +396,109 @@ struct rparams
OldBezier &A;
OldBezier &B;
};
-
+
static int
intersect_polish_f (const gsl_vector * x, void *params,
gsl_vector * f)
{
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;
}
+typedef union dbl_64{
+ long long i64;
+ double d64;
+};
+
+static double EpsilonBy(double value, int eps)
+{
+ dbl_64 s;
+ s.d64 = value;
+ s.i64 += eps;
+ return s.d64;
+}
+
+
static void intersect_polish_root (OldBezier &A, double &s,
OldBezier &B, double &t) {
const gsl_multiroot_fsolver_type *T;
gsl_multiroot_fsolver *sol;
-
+
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]);
-
+
T = gsl_multiroot_fsolver_hybrids;
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 */
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);
+
+ // This code does a neighbourhood search for minor improvements.
+ double best_v = L1(A(s) - B(t));
+ //std::cout << "------\n" << best_v << std::endl;
+ Point best(s,t);
+ while (true) {
+ Point trial = best;
+ double trial_v = best_v;
+ for(int nsi = -1; nsi < 2; nsi++) {
+ for(int nti = -1; nti < 2; nti++) {
+ Point n(EpsilonBy(best[0], nsi),
+ EpsilonBy(best[1], nti));
+ double c = L1(A(n[0]) - B(n[1]));
+ //std::cout << c << "; ";
+ if (c < trial_v) {
+ trial = n;
+ trial_v = c;
+ }
+ }
+ }
+ if(trial == best) {
+ //std::cout << "\n" << s << " -> " << s - best[0] << std::endl;
+ //std::cout << t << " -> " << t - best[1] << std::endl;
+ //std::cout << best_v << std::endl;
+ s = best[0];
+ t = best[1];
+ return;
+ } else {
+ best = trial;
+ best_v = trial_v;
+ }
+ }
}
@@ -432,8 +507,8 @@ std::vector<std::pair<double, double> > find_intersections( OldBezier a, OldBezi
std::vector<std::pair<double, double> > parameters;
if( intersect_BB( a, b ) )
{
- recursively_intersect( a, 0., 1., wangs_theorem(a),
- b, 0., 1., wangs_theorem(b),
+ recursively_intersect( a, 0., 1., wangs_theorem(a),
+ b, 0., 1., wangs_theorem(b),
parameters);
}
for(unsigned i = 0; i < parameters.size(); i++)
generated by cgit v1.2.3 (git 2.25.1) at 2026-05-19 16:24:18 +0000