summaryrefslogtreecommitdiffstats
path: root/src/2geom/orphan-code
diff options
context:
space:
mode:
authorDaniel Baumann <daniel.baumann@progress-linux.org>2024-04-13 11:57:42 +0000
committerDaniel Baumann <daniel.baumann@progress-linux.org>2024-04-13 11:57:42 +0000
commit61f3ab8f23f4c924d455757bf3e65f8487521b5a (patch)
tree885599a36a308f422af98616bc733a0494fe149a /src/2geom/orphan-code
parentInitial commit. (diff)
downloadlib2geom-upstream.tar.xz
lib2geom-upstream.zip
Adding upstream version 1.3.upstream/1.3upstream
Signed-off-by: Daniel Baumann <daniel.baumann@progress-linux.org>
Diffstat (limited to '')
-rw-r--r--src/2geom/orphan-code/arc-length.cpp292
-rw-r--r--src/2geom/orphan-code/chebyshev.cpp126
-rw-r--r--src/2geom/orphan-code/intersection-by-bezier-clipping.cpp560
-rw-r--r--src/2geom/orphan-code/intersection-by-smashing.cpp349
-rw-r--r--src/2geom/orphan-code/nearestpoint.cpp405
-rw-r--r--src/2geom/orphan-code/redblack-toy.cpp327
-rw-r--r--src/2geom/orphan-code/redblacktree.cpp575
-rw-r--r--src/2geom/orphan-code/rtree.cpp1350
8 files changed, 3984 insertions, 0 deletions
diff --git a/src/2geom/orphan-code/arc-length.cpp b/src/2geom/orphan-code/arc-length.cpp
new file mode 100644
index 0000000..3f72862
--- /dev/null
+++ b/src/2geom/orphan-code/arc-length.cpp
@@ -0,0 +1,292 @@
+/*
+ * arc-length.cpp
+ *
+ * Copyright 2006 Nathan Hurst <njh@mail.csse.monash.edu.au>
+ * Copyright 2006 Michael G. Sloan <mgsloan@gmail.com>
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ *
+ */
+
+#include <2geom/arc-length.h>
+#include <2geom/bezier-utils.h>
+#include <2geom/polynomial.h>
+using namespace Geom;
+
+/** Calculates the length of a cubic element through subdivision.
+ * The 'tol' parameter is the maximum error allowed. This is used to subdivide the curve where necessary.
+ */
+double cubic_length_subdividing(Path::Elem const & e, double tol) {
+ Point v[3];
+ for(int i = 0; i < 3; i++)
+ v[i] = e[i+1] - e[0];
+ Point orth = v[2]; // unit normal to path line
+ rot90(orth);
+ orth.normalize();
+ double err = fabs(dot(orth, v[1])) + fabs(dot(orth, v[0]));
+ if(err < tol) {
+ return distance(e.first(), e.last()); // approximately a line
+ } else {
+ Point mid[3];
+ double result;
+ for(int i = 0; i < 3; i++)
+ mid[i] = lerp(0.5, e[i], e[i+1]);
+ Point midmid[2];
+ for(int i = 0; i < 2; i++)
+ midmid[i] = lerp(0.5, mid[i], mid[i+1]);
+ Point midmidmid = lerp(0.5, midmid[0], midmid[1]);
+ {
+ Point curve[4] = {e[0], mid[0], midmid[0], midmidmid};
+ Path::Elem e0(cubicto, std::vector<Point>::const_iterator(curve), std::vector<Point>::const_iterator(curve) + 4);
+ result = cubic_length_subdividing(e0, tol);
+ } {
+ Point curve[4] = {midmidmid, midmid[1], mid[2], e[3]};
+ Path::Elem e1(cubicto, std::vector<Point>::const_iterator(curve), std::vector<Point>::const_iterator(curve) + 4);
+ return result + cubic_length_subdividing(e1, tol);
+ }
+ }
+}
+
+/** Calculates the length of a path through iteration and subsequent subdivision.
+ * Currently handles cubic curves and lines.
+ * The 'tol' parameter is the maximum error allowed. This is used to subdivide the curve where necessary.
+ */
+double arc_length_subdividing(Path const & p, double tol) {
+ double result = 0;
+
+ for(Path::const_iterator iter(p.begin()), end(p.end()); iter != end; ++iter) {
+ if(dynamic_cast<LineTo *>(iter.cmd()))
+ result += distance((*iter).first(), (*iter).last());
+ else if(dynamic_cast<CubicTo *>(iter.cmd()))
+ result += cubic_length_subdividing(*iter, tol);
+ else
+ ;
+ }
+
+ return result;
+}
+
+
+#ifdef HAVE_GSL
+#include <gsl/gsl_integration.h>
+static double poly_length_integrating(double t, void* param) {
+ Poly* pc = (Poly*)param;
+ return hypot(pc[0].eval(t), pc[1].eval(t));
+}
+
+/** Calculates the length of a path Element through gsl integration.
+ \param pe the Element.
+ \param t the parametric input 0 to 1 which specifies the amount of the curve to use.
+ \param tol the maximum error allowed.
+ \param result variable to be incremented with the length of the path
+ \param abs_error variable to be incremented with the estimated error
+*/
+void arc_length_integrating(Path::Elem pe, double t, double tol, double &result, double &abs_error) {
+ if(dynamic_cast<LineTo *>(iter.cmd()))
+ result += distance(pe.first(), pe.last()) * t;
+ else if(dynamic_cast<QuadTo *>(iter.cmd()) ||
+ dynamic_cast<CubicTo *>(iter.cmd())) {
+ Poly B[2] = {get_parametric_poly(pe, X), get_parametric_poly(pe, Y)};
+ for(int i = 0; i < 2; i++)
+ B[i] = derivative(B[i]);
+
+ gsl_function F;
+ gsl_integration_workspace * w
+ = gsl_integration_workspace_alloc (20);
+ F.function = &poly_length_integrating;
+ F.params = (void*)B;
+ double quad_result, err;
+ /* We could probably use the non adaptive code here if we removed any cusps first. */
+ int returncode =
+ gsl_integration_qag (&F, 0, t, 0, tol, 20,
+ GSL_INTEG_GAUSS21, w, &quad_result, &err);
+
+ abs_error += err;
+ result += quad_result;
+ } else
+ return;
+}
+
+/** Calculates the length of a Path through gsl integration. The parameter 'tol' is the maximum error allowed. */
+double arc_length_integrating(Path const & p, double tol) {
+ double result = 0, abserr = 0;
+
+ for(Path::const_iterator iter(p.begin()), end(p.end()); iter != end; ++iter) {
+ arc_length_integrating(*iter, 1.0, tol, result, abserr);
+ }
+ //printf("got %g with err %g\n", result, abserr);
+
+ return result;
+}
+
+/** Calculates the arc length to a specific location on the path. The parameter 'tol' is the maximum error allowed. */
+double arc_length_integrating(Path const & p, Path::Location const & pl, double tol) {
+ double result = 0, abserr = 0;
+ ptrdiff_t offset = pl.it - p.begin();
+
+ assert(offset >= 0);
+ assert(offset < p.size());
+
+ for(Path::const_iterator iter(p.begin()), end(p.end());
+ (iter != pl.it); ++iter) {
+ arc_length_integrating(*iter, 1.0, tol, result, abserr);
+ }
+ arc_length_integrating(*pl.it, pl.t, tol, result, abserr);
+
+ return result;
+}
+
+/* We use a somewhat surprising result for this that s'(t) = |p'(t)|
+ Thus, we can use a derivative based root finder.
+*/
+
+#include <stdio.h>
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_roots.h>
+
+struct arc_length_params
+{
+ Path::Elem pe;
+ double s,tol, result, abs_error;
+ double left, right;
+};
+
+double
+arc_length (double t, void *params)
+{
+ struct arc_length_params *p
+ = (struct arc_length_params *) params;
+
+ double result = 0, abs_error = 0;
+ if(t < 0) t = 0;
+ if(t > 1) t = 1;
+ if(!((t >= 0) && (t <= 1))) {
+ printf("t = %g\n", t);
+ }
+ assert((t >= 0) && (t <= 1));
+ arc_length_integrating(p->pe, t, p->tol, result, abs_error);
+ return result - p->s ;
+}
+
+double
+arc_length_deriv (double t, void *params)
+{
+ struct arc_length_params *p
+ = (struct arc_length_params *) params;
+
+ Point pos, tgt, acc;
+ p->pe.point_tangent_acc_at(t, pos, tgt, acc);
+ return L2(tgt);
+}
+
+void
+arc_length_fdf (double t, void *params,
+ double *y, double *dy)
+{
+ *y = arc_length(t, params);
+ *dy = arc_length_deriv(t, params);
+}
+
+double polish_brent(double t, arc_length_params &alp) {
+ int status;
+ int iter = 0, max_iter = 10;
+ const gsl_root_fsolver_type *T;
+ gsl_root_fsolver *solver;
+ double x_lo = 0.0, x_hi = 1.0;
+ gsl_function F;
+
+ F.function = &arc_length;
+ F.params = &alp;
+
+ T = gsl_root_fsolver_brent;
+ solver = gsl_root_fsolver_alloc (T);
+ gsl_root_fsolver_set (solver, &F, x_lo, x_hi);
+
+ do
+ {
+ iter++;
+ status = gsl_root_fsolver_iterate (solver);
+ t = gsl_root_fsolver_root (solver);
+ x_lo = gsl_root_fsolver_x_lower (solver);
+ x_hi = gsl_root_fsolver_x_upper (solver);
+ status = gsl_root_test_interval (x_lo, x_hi,
+ 0, alp.tol);
+
+ //if (status == GSL_SUCCESS)
+ // printf ("Converged:\n");
+
+ }
+ while (status == GSL_CONTINUE && iter < max_iter);
+ return t;
+}
+
+double polish (double t, arc_length_params &alp) {
+ int status;
+ int iter = 0, max_iter = 5;
+ const gsl_root_fdfsolver_type *T;
+ gsl_root_fdfsolver *solver;
+ double t0;
+ gsl_function_fdf FDF;
+
+ FDF.f = &arc_length;
+ FDF.df = &arc_length_deriv;
+ FDF.fdf = &arc_length_fdf;
+ FDF.params = &alp;
+
+ T = gsl_root_fdfsolver_newton;
+ solver = gsl_root_fdfsolver_alloc (T);
+ gsl_root_fdfsolver_set (solver, &FDF, t);
+
+ do
+ {
+ iter++;
+ status = gsl_root_fdfsolver_iterate (solver);
+ t0 = t;
+ t = gsl_root_fdfsolver_root (solver);
+ status = gsl_root_test_delta (t, t0, 0, alp.tol);
+
+ if (status == GSL_SUCCESS)
+ ;//printf ("Converged:\n");
+
+ printf ("%5d %10.7f %+10.7f\n",
+ iter, t, t - t0);
+ }
+ while (status == GSL_CONTINUE && iter < max_iter);
+ return t;
+}
+
+
+#endif
+
+/*
+ Local Variables:
+ mode:c++
+ c-file-style:"stroustrup"
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+ indent-tabs-mode:nil
+ fill-column:99
+ End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
diff --git a/src/2geom/orphan-code/chebyshev.cpp b/src/2geom/orphan-code/chebyshev.cpp
new file mode 100644
index 0000000..c886daf
--- /dev/null
+++ b/src/2geom/orphan-code/chebyshev.cpp
@@ -0,0 +1,126 @@
+#include <2geom/chebyshev.h>
+
+#include <2geom/sbasis.h>
+#include <2geom/sbasis-poly.h>
+
+#include <vector>
+using std::vector;
+
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_chebyshev.h>
+
+namespace Geom{
+
+SBasis cheb(unsigned n) {
+ static std::vector<SBasis> basis;
+ if(basis.empty()) {
+ basis.push_back(Linear(1,1));
+ basis.push_back(Linear(0,1));
+ }
+ for(unsigned i = basis.size(); i <= n; i++) {
+ basis.push_back(Linear(0,2)*basis[i-1] - basis[i-2]);
+ }
+
+ return basis[n];
+}
+
+SBasis cheb_series(unsigned n, double* cheb_coeff) {
+ SBasis r;
+ for(unsigned i = 0; i < n; i++) {
+ double cof = cheb_coeff[i];
+ //if(i == 0)
+ //cof /= 2;
+ r += cheb(i)*cof;
+ }
+
+ return r;
+}
+
+SBasis clenshaw_series(unsigned m, double* cheb_coeff) {
+ /** b_n = a_n
+ b_n-1 = 2*x*b_n + a_n-1
+ b_n-k = 2*x*b_{n-k+1} + a_{n-k} - b_{n - k + 2}
+ b_0 = x*b_1 + a_0 - b_2
+ */
+
+ double a = -1, b = 1;
+ SBasis d, dd;
+ SBasis y = (Linear(0, 2) - (a+b)) / (b-a);
+ SBasis y2 = 2*y;
+ for(int j = m - 1; j >= 1; j--) {
+ SBasis sv = d;
+ d = y2*d - dd + cheb_coeff[j];
+ dd = sv;
+ }
+
+ return y*d - dd + 0.5*cheb_coeff[0];
+}
+
+SBasis chebyshev_approximant (double (*f)(double,void*), int order, Interval in, void* p) {
+ gsl_cheb_series *cs = gsl_cheb_alloc (order+2);
+
+ gsl_function F;
+
+ F.function = f;
+ F.params = p;
+
+ gsl_cheb_init (cs, &F, in[0], in[1]);
+
+ SBasis r = compose(clenshaw_series(order, cs->c), Linear(-1,1));
+
+ gsl_cheb_free (cs);
+ return r;
+}
+
+struct wrap {
+ double (*f)(double,void*);
+ void* pp;
+ double fa, fb;
+ Interval in;
+};
+
+double f_interp(double x, void* p) {
+ struct wrap *wr = (struct wrap *)p;
+ double z = (x - wr->in[0]) / (wr->in[1] - wr->in[0]);
+ return (wr->f)(x, wr->pp) - ((1 - z)*wr->fa + z*wr->fb);
+}
+
+SBasis chebyshev_approximant_interpolating (double (*f)(double,void*),
+ int order, Interval in, void* p) {
+ double fa = f(in[0], p);
+ double fb = f(in[1], p);
+ struct wrap wr;
+ wr.fa = fa;
+ wr.fb = fb;
+ wr.in = in;
+ printf("%f %f\n", fa, fb);
+ wr.f = f;
+ wr.pp = p;
+ return compose(Linear(in[0], in[1]), Linear(fa, fb)) + chebyshev_approximant(f_interp, order, in, &wr) + Linear(fa, fb);
+}
+
+SBasis chebyshev(unsigned n) {
+ static std::vector<SBasis> basis;
+ if(basis.empty()) {
+ basis.push_back(Linear(1,1));
+ basis.push_back(Linear(0,1));
+ }
+ for(unsigned i = basis.size(); i <= n; i++) {
+ basis.push_back(Linear(0,2)*basis[i-1] - basis[i-2]);
+ }
+
+ return basis[n];
+}
+
+};
+
+/*
+ Local Variables:
+ mode:c++
+ c-file-style:"stroustrup"
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+ indent-tabs-mode:nil
+ fill-column:99
+ End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
diff --git a/src/2geom/orphan-code/intersection-by-bezier-clipping.cpp b/src/2geom/orphan-code/intersection-by-bezier-clipping.cpp
new file mode 100644
index 0000000..c55f623
--- /dev/null
+++ b/src/2geom/orphan-code/intersection-by-bezier-clipping.cpp
@@ -0,0 +1,560 @@
+
+/*
+ * Find intersecions between two Bezier curves.
+ * The intersection points are found by using Bezier clipping.
+ *
+ * Authors:
+ * Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2008 authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+
+
+
+
+#include <2geom/basic-intersection.h>
+#include <2geom/bezier.h>
+#include <2geom/interval.h>
+#include <2geom/convex-hull.h>
+
+
+#include <vector>
+#include <utility>
+#include <iomanip>
+
+
+namespace Geom {
+
+namespace detail { namespace bezier_clipping {
+
+
+////////////////////////////////////////////////////////////////////////////////
+// for debugging
+//
+
+inline
+void print(std::vector<Point> const& cp)
+{
+ for (size_t i = 0; i < cp.size(); ++i)
+ std::cerr << i << " : " << cp[i] << std::endl;
+}
+
+template< class charT >
+inline
+std::basic_ostream<charT> &
+operator<< (std::basic_ostream<charT> & os, const Interval & I)
+{
+ os << "[" << I.min() << ", " << I.max() << "]";
+ return os;
+}
+
+inline
+double angle (std::vector<Point> const& A)
+{
+ size_t n = A.size() -1;
+ double a = std::atan2(A[n][Y] - A[0][Y], A[n][X] - A[0][X]);
+ return (180 * a / M_PI);
+}
+
+inline
+size_t get_precision(Interval const& I)
+{
+ double d = I.extent();
+ double e = 1, p = 1;
+ size_t n = 0;
+ while (n < 16 && (are_near(d, 0, e)))
+ {
+ p *= 10;
+ e = 1 /p;
+ ++n;
+ }
+ return n;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+
+
+/*
+ * return true if all the Bezier curve control points are near,
+ * false otherwise
+ */
+inline
+bool is_constant(std::vector<Point> const& A, double precision = EPSILON)
+{
+ for (unsigned int i = 1; i < A.size(); ++i)
+ {
+ if(!are_near(A[i], A[0], precision))
+ return false;
+ }
+ return true;
+}
+
+/*
+ * Make up an orientation line using the control points c[i] and c[j]
+ * the line is returned in the output parameter "l" in the form of a 3 element
+ * vector : l[0] * x + l[1] * y + l[2] == 0; the line is normalized.
+ */
+inline
+void orientation_line (std::vector<double> & l,
+ std::vector<Point> const& c,
+ size_t i, size_t j)
+{
+ l[0] = c[j][Y] - c[i][Y];
+ l[1] = c[i][X] - c[j][X];
+ l[2] = cross(c[i], c[j]);
+ double length = std::sqrt(l[0] * l[0] + l[1] * l[1]);
+ assert (length != 0);
+ l[0] /= length;
+ l[1] /= length;
+ l[2] /= length;
+}
+
+/*
+ * Pick up an orientation line for the Bezier curve "c" and return it in
+ * the output parameter "l"
+ */
+inline
+void pick_orientation_line (std::vector<double> & l,
+ std::vector<Point> const& c)
+{
+ size_t i = c.size();
+ while (--i > 0 && are_near(c[0], c[i]))
+ {}
+ if (i == 0)
+ {
+ // this should never happen because when a new curve portion is created
+ // we check that it is not constant;
+ // however this requires that the precision used in the is_constant
+ // routine has to be the same used here in the are_near test
+ assert(i != 0);
+ }
+ orientation_line(l, c, 0, i);
+ //std::cerr << "i = " << i << std::endl;
+}
+
+/*
+ * Compute the signed distance of the point "P" from the normalized line l
+ */
+inline
+double distance (Point const& P, std::vector<double> const& l)
+{
+ return l[X] * P[X] + l[Y] * P[Y] + l[2];
+}
+
+/*
+ * Compute the min and max distance of the control points of the Bezier
+ * curve "c" from the normalized orientation line "l".
+ * This bounds are returned through the output Interval parameter"bound".
+ */
+inline
+void fat_line_bounds (Interval& bound,
+ std::vector<Point> const& c,
+ std::vector<double> const& l)
+{
+ bound[0] = 0;
+ bound[1] = 0;
+ double d;
+ for (size_t i = 0; i < c.size(); ++i)
+ {
+ d = distance(c[i], l);
+ if (bound[0] > d) bound[0] = d;
+ if (bound[1] < d) bound[1] = d;
+ }
+}
+
+/*
+ * return the x component of the intersection point between the line
+ * passing through points p1, p2 and the line Y = "y"
+ */
+inline
+double intersect (Point const& p1, Point const& p2, double y)
+{
+ // we are sure that p2[Y] != p1[Y] because this routine is called
+ // only when the lower or the upper bound is crossed
+ double dy = (p2[Y] - p1[Y]);
+ double s = (y - p1[Y]) / dy;
+ return (p2[X]-p1[X])*s + p1[X];
+}
+
+/*
+ * Clip the Bezier curve "B" wrt the fat line defined by the orientation
+ * line "l" and the interval range "bound", the new parameter interval for
+ * the clipped curve is returned through the output parameter "dom"
+ */
+void clip (Interval& dom,
+ std::vector<Point> const& B,
+ std::vector<double> const& l,
+ Interval const& bound)
+{
+ double n = B.size() - 1; // number of sub-intervals
+ std::vector<Point> D; // distance curve control points
+ D.reserve (B.size());
+ double d;
+ for (size_t i = 0; i < B.size(); ++i)
+ {
+ d = distance (B[i], l);
+ D.push_back (Point(i/n, d));
+ }
+ //print(D);
+ ConvexHull chD(D);
+ std::vector<Point> & p = chD.boundary; // convex hull vertices
+
+ //print(p);
+
+ bool plower, phigher;
+ bool clower, chigher;
+ double t, tmin = 1, tmax = 0;
+ //std::cerr << "bound : " << bound << std::endl;
+
+ plower = (p[0][Y] < bound.min());
+ phigher = (p[0][Y] > bound.max());
+ if (!(plower || phigher)) // inside the fat line
+ {
+ if (tmin > p[0][X]) tmin = p[0][X];
+ if (tmax < p[0][X]) tmax = p[0][X];
+ //std::cerr << "0 : inside " << p[0]
+ // << " : tmin = " << tmin << ", tmax = " << tmax << std::endl;
+ }
+
+ for (size_t i = 1; i < p.size(); ++i)
+ {
+ clower = (p[i][Y] < bound.min());
+ chigher = (p[i][Y] > bound.max());
+ if (!(clower || chigher)) // inside the fat line
+ {
+ if (tmin > p[i][X]) tmin = p[i][X];
+ if (tmax < p[i][X]) tmax = p[i][X];
+ //std::cerr << i << " : inside " << p[i]
+ // << " : tmin = " << tmin << ", tmax = " << tmax
+ // << std::endl;
+ }
+ if (clower != plower) // cross the lower bound
+ {
+ t = intersect(p[i-1], p[i], bound.min());
+ if (tmin > t) tmin = t;
+ if (tmax < t) tmax = t;
+ plower = clower;
+ //std::cerr << i << " : lower " << p[i]
+ // << " : tmin = " << tmin << ", tmax = " << tmax
+ // << std::endl;
+ }
+ if (chigher != phigher) // cross the upper bound
+ {
+ t = intersect(p[i-1], p[i], bound.max());
+ if (tmin > t) tmin = t;
+ if (tmax < t) tmax = t;
+ phigher = chigher;
+ //std::cerr << i << " : higher " << p[i]
+ // << " : tmin = " << tmin << ", tmax = " << tmax
+ // << std::endl;
+ }
+ }
+
+ // we have to test the closing segment for intersection
+ size_t last = p.size() - 1;
+ clower = (p[0][Y] < bound.min());
+ chigher = (p[0][Y] > bound.max());
+ if (clower != plower) // cross the lower bound
+ {
+ t = intersect(p[last], p[0], bound.min());
+ if (tmin > t) tmin = t;
+ if (tmax < t) tmax = t;
+ //std::cerr << "0 : lower " << p[0]
+ // << " : tmin = " << tmin << ", tmax = " << tmax << std::endl;
+ }
+ if (chigher != phigher) // cross the upper bound
+ {
+ t = intersect(p[last], p[0], bound.max());
+ if (tmin > t) tmin = t;
+ if (tmax < t) tmax = t;
+ //std::cerr << "0 : higher " << p[0]
+ // << " : tmin = " << tmin << ", tmax = " << tmax << std::endl;
+ }
+
+ dom[0] = tmin;
+ dom[1] = tmax;
+}
+
+/*
+ * Compute the portion of the Bezier curve "B" wrt the interval "I"
+ */
+void portion (std::vector<Point> & B, Interval const& I)
+{
+ Bezier::Order bo(B.size()-1);
+ Bezier Bx(bo), By(bo);
+ for (size_t i = 0; i < B.size(); ++i)
+ {
+ Bx[i] = B[i][X];
+ By[i] = B[i][Y];
+ }
+ Bx = portion(Bx, I.min(), I.max());
+ By = portion(By, I.min(), I.max());
+ assert (Bx.size() == By.size());
+ B.resize(Bx.size());
+ for (size_t i = 0; i < Bx.size(); ++i)
+ {
+ B[i][X] = Bx[i];
+ B[i][Y] = By[i];
+ }
+}
+
+/*
+ * Map the sub-interval I in [0,1] into the interval J and assign it to J
+ */
+inline
+void map_to(Interval & J, Interval const& I)
+{
+ double length = J.extent();
+ J[1] = I.max() * length + J[0];
+ J[0] = I.min() * length + J[0];
+}
+
+/*
+ * The interval [1,0] is used to represent the empty interval, this routine
+ * is just an helper function for creating such an interval
+ */
+inline
+Interval make_empty_interval()
+{
+ Interval I(0);
+ I[0] = 1;
+ return I;
+}
+
+
+
+
+const double MAX_PRECISION = 1e-8;
+const double MIN_CLIPPED_SIZE_THRESHOLD = 0.8;
+const Interval UNIT_INTERVAL(0,1);
+const Interval EMPTY_INTERVAL = make_empty_interval();
+const Interval H1_INTERVAL(0, 0.5);
+const Interval H2_INTERVAL(0.5 + MAX_PRECISION, 1.0);
+
+/*
+ * intersection
+ *
+ * input:
+ * A, B: control point sets of two bezier curves
+ * domA, domB: real parameter intervals of the two curves
+ * precision: required computational precision of the returned parameter ranges
+ * output:
+ * domsA, domsB: sets of parameter intervals describing an intersection point
+ *
+ * The parameter intervals are computed by using a Bezier clipping algorithm,
+ * in case the clipping doesn't shrink the initial interval more than 20%,
+ * a subdivision step is performed.
+ * If during the computation one of the two curve interval length becomes less
+ * than MAX_PRECISION the routine exits independently by the precision reached
+ * in the computation of the other curve interval.
+ */
+void intersection (std::vector<Interval>& domsA,
+ std::vector<Interval>& domsB,
+ std::vector<Point> const& A,
+ std::vector<Point> const& B,
+ Interval const& domA,
+ Interval const& domB,
+ double precision)
+{
+// std::cerr << ">> curve subdision performed <<" << std::endl;
+// std::cerr << "dom(A) : " << domA << std::endl;
+// std::cerr << "dom(B) : " << domB << std::endl;
+// std::cerr << "angle(A) : " << angle(A) << std::endl;
+// std::cerr << "angle(B) : " << angle(B) << std::endl;
+
+
+ if (precision < MAX_PRECISION)
+ precision = MAX_PRECISION;
+
+ std::vector<Point> pA = A;
+ std::vector<Point> pB = B;
+ std::vector<Point>* C1 = &pA;
+ std::vector<Point>* C2 = &pB;
+
+ Interval dompA = domA;
+ Interval dompB = domB;
+ Interval* dom1 = &dompA;
+ Interval* dom2 = &dompB;
+
+ std::vector<double> bl(3);
+ Interval bound, dom;
+
+
+ size_t iter = 0;
+ while (++iter < 100
+ && (dompA.extent() >= precision || dompB.extent() >= precision))
+ {
+// std::cerr << "iter: " << iter << std::endl;
+
+ pick_orientation_line(bl, *C1);
+ fat_line_bounds(bound, *C1, bl);
+ clip(dom, *C2, bl, bound);
+
+ // [1,0] is utilized to represent an empty interval
+ if (dom == EMPTY_INTERVAL)
+ {
+// std::cerr << "dom: empty" << std::endl;
+ return;
+ }
+// std::cerr << "dom : " << dom << std::endl;
+
+ // all other cases where dom[0] > dom[1] are invalid
+ if (dom.min() > dom.max())
+ {
+ assert(dom.min() < dom.max());
+ }
+
+ map_to(*dom2, dom);
+
+ // it's better to stop before losing computational precision
+ if (dom2->extent() <= MAX_PRECISION)
+ {
+// std::cerr << "beyond max precision limit" << std::endl;
+ break;
+ }
+
+ portion(*C2, dom);
+ if (is_constant(*C2))
+ {
+// std::cerr << "new curve portion is constant" << std::endl;
+ break;
+ }
+ // if we have clipped less than 20% than we need to subdive the curve
+ // with the largest domain into two sub-curves
+ if (dom.extent() > MIN_CLIPPED_SIZE_THRESHOLD)
+ {
+// std::cerr << "clipped less than 20% : " << dom.extent() << std::endl;
+// std::cerr << "angle(pA) : " << angle(pA) << std::endl;
+// std::cerr << "angle(pB) : " << angle(pB) << std::endl;
+
+ std::vector<Point> pC1, pC2;
+ Interval dompC1, dompC2;
+ if (dompA.extent() > dompB.extent())
+ {
+ if ((dompA.extent() / 2) < MAX_PRECISION)
+ {
+ break;
+ }
+ pC1 = pC2 = pA;
+ portion(pC1, H1_INTERVAL);
+ portion(pC2, H2_INTERVAL);
+ dompC1 = dompC2 = dompA;
+ map_to(dompC1, H1_INTERVAL);
+ map_to(dompC2, H2_INTERVAL);
+ intersection(domsA, domsB, pC1, pB, dompC1, dompB, precision);
+ intersection(domsA, domsB, pC2, pB, dompC2, dompB, precision);
+ }
+ else
+ {
+ if ((dompB.extent() / 2) < MAX_PRECISION)
+ {
+ break;
+ }
+ pC1 = pC2 = pB;
+ portion(pC1, H1_INTERVAL);
+ portion(pC2, H2_INTERVAL);
+ dompC1 = dompC2 = dompB;
+ map_to(dompC1, H1_INTERVAL);
+ map_to(dompC2, H2_INTERVAL);
+ intersection(domsB, domsA, pC1, pA, dompC1, dompA, precision);
+ intersection(domsB, domsA, pC2, pA, dompC2, dompA, precision);
+ }
+ return;
+ }
+
+ using std::swap;
+ swap(C1, C2);
+ swap(dom1, dom2);
+// std::cerr << "dom(pA) : " << dompA << std::endl;
+// std::cerr << "dom(pB) : " << dompB << std::endl;
+ }
+ domsA.push_back(dompA);
+ domsB.push_back(dompB);
+}
+
+} /* end namespace bezier_clipping */ } /* end namespace detail */
+
+
+/*
+ * find_intersection
+ *
+ * input: A, B - set of control points of two Bezier curve
+ * input: precision - required precision of computation
+ * output: xs - set of pairs of parameter values
+ * at which crossing happens
+ *
+ * This routine is based on the Bezier Clipping Algorithm,
+ * see: Sederberg - Computer Aided Geometric Design
+ */
+void find_intersections (std::vector< std::pair<double, double> > & xs,
+ std::vector<Point> const& A,
+ std::vector<Point> const& B,
+ double precision)
+{
+ std::cout << "find_intersections: intersection-by-clipping.cpp version\n";
+// std::cerr << std::fixed << std::setprecision(16);
+
+ using detail::bezier_clipping::get_precision;
+ using detail::bezier_clipping::operator<<;
+ using detail::bezier_clipping::intersection;
+ using detail::bezier_clipping::UNIT_INTERVAL;
+
+ std::pair<double, double> ci;
+ std::vector<Interval> domsA, domsB;
+ intersection (domsA, domsB, A, B, UNIT_INTERVAL, UNIT_INTERVAL, precision);
+ if (domsA.size() != domsB.size())
+ {
+ assert (domsA.size() == domsB.size());
+ }
+ xs.clear();
+ xs.reserve(domsA.size());
+ for (size_t i = 0; i < domsA.size(); ++i)
+ {
+// std::cerr << i << " : domA : " << domsA[i] << std::endl;
+// std::cerr << "precision A: " << get_precision(domsA[i]) << std::endl;
+// std::cerr << i << " : domB : " << domsB[i] << std::endl;
+// std::cerr << "precision B: " << get_precision(domsB[i]) << std::endl;
+
+ ci.first = domsA[i].middle();
+ ci.second = domsB[i].middle();
+ xs.push_back(ci);
+ }
+}
+
+} // end namespace Geom
+
+
+/*
+ Local Variables:
+ mode:c++
+ c-file-style:"stroustrup"
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+ indent-tabs-mode:nil
+ fill-column:99
+ End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
diff --git a/src/2geom/orphan-code/intersection-by-smashing.cpp b/src/2geom/orphan-code/intersection-by-smashing.cpp
new file mode 100644
index 0000000..02e44b1
--- /dev/null
+++ b/src/2geom/orphan-code/intersection-by-smashing.cpp
@@ -0,0 +1,349 @@
+#include <2geom/d2.h>
+#include <2geom/sbasis.h>
+#include <2geom/sbasis-geometric.h>
+#include <2geom/orphan-code/intersection-by-smashing.h>
+
+#include <cstdlib>
+#include <cstdio>
+#include <vector>
+#include <algorithm>
+
+namespace Geom {
+using namespace Geom;
+
+/*
+ * Computes the top and bottom boundaries of the L_\infty neighborhood
+ * of a curve. The curve is supposed to be a graph over the x-axis.
+ */
+static
+void computeLinfinityNeighborhood( D2<SBasis > const &f, double tol, D2<Piecewise<SBasis> > &topside, D2<Piecewise<SBasis> > &botside ){
+ double signx = ( f[X].at0() > f[X].at1() )? -1 : 1;
+ double signy = ( f[Y].at0() > f[Y].at1() )? -1 : 1;
+
+ Piecewise<D2<SBasis> > top, bot;
+ top = Piecewise<D2<SBasis> > (f);
+ top.cuts.insert( top.cuts.end(), 2);
+ top.segs.insert( top.segs.end(), D2<SBasis>(SBasis(Linear( f[X].at1(), f[X].at1()+2*tol*signx)),
+ SBasis(Linear( f[Y].at1() )) ));
+ bot = Piecewise<D2<SBasis> >(f);
+ bot.cuts.insert( bot.cuts.begin(), - 1 );
+ bot.segs.insert( bot.segs.begin(), D2<SBasis>(SBasis(Linear( f[X].at0()-2*tol*signx, f[X].at0())),
+ SBasis(Linear( f[Y].at0() )) ));
+ top += Point(-tol*signx, tol);
+ bot += Point( tol*signx, -tol);
+
+ if ( signy < 0 ){
+ std::swap( top, bot );
+ top += Point( 0, 2*tol);
+ bot += Point( 0, -2*tol);
+ }
+ topside = make_cuts_independent(top);
+ botside = make_cuts_independent(bot);
+}
+
+
+/*
+ * Compute top and bottom boundaries of the L^infty nbhd of the graph of a *monotonic* function f.
+ * if f is increasing, it is given by [f(t-tol)-tol, f(t+tol)+tol].
+ * if not, it is [f(t+tol)-tol, f(t-tol)+tol].
+ */
+static
+void computeLinfinityNeighborhood( Piecewise<SBasis> const &f, double tol, Piecewise<SBasis> &top, Piecewise<SBasis> &bot){
+ top = f + tol;
+ top.offsetDomain( - tol );
+ top.cuts.insert( top.cuts.end(), f.domain().max() + tol);
+ top.segs.insert( top.segs.end(), SBasis(Linear( f.lastValue() + tol )) );
+
+ bot = f - tol;
+ bot.offsetDomain( tol );
+ bot.cuts.insert( bot.cuts.begin(), f.domain().min() - tol);
+ bot.segs.insert( bot.segs.begin(), SBasis(Linear( f.firstValue() - tol )) );
+
+ if (f.firstValue() > f.lastValue()) {
+ std::swap(top, bot);
+ top += 2 * tol;
+ bot -= 2 * tol;
+ }
+}
+
+/*
+ * Returns the intervals over which the curve keeps its slope
+ * in one of the 8 sectors delimited by x=0, y=0, y=x, y=-x.
+ */
+std::vector<Interval> monotonicSplit(D2<SBasis> const &p){
+ std::vector<Interval> result;
+
+ D2<SBasis> v = derivative(p);
+
+ std::vector<double> someroots;
+ std::vector<double> cuts (2,0.);
+ cuts[1] = 1.;
+
+ someroots = roots(v[X]);
+ cuts.insert( cuts.end(), someroots.begin(), someroots.end() );
+
+ someroots = roots(v[Y]);
+ cuts.insert( cuts.end(), someroots.begin(), someroots.end() );
+
+ //we could split in the middle to avoid computing roots again...
+ someroots = roots(v[X]-v[Y]);
+ cuts.insert( cuts.end(), someroots.begin(), someroots.end() );
+
+ someroots = roots(v[X]+v[Y]);
+ cuts.insert( cuts.end(), someroots.begin(), someroots.end() );
+
+ sort(cuts.begin(),cuts.end());
+ unique(cuts.begin(), cuts.end() );
+
+ for (unsigned i=1; i<cuts.size(); i++){
+ result.push_back( Interval( cuts[i-1], cuts[i] ) );
+ }
+ return result;
+}
+
+//std::vector<Interval> level_set( D2<SBasis> const &f, Rect region){
+// std::vector<Interval> x_in_reg = level_set( f[X], region[X] );
+// std::vector<Interval> y_in_reg = level_set( f[Y], region[Y] );
+// std::vector<Interval> result = intersect ( x_in_reg, y_in_reg );
+// return result;
+//}
+
+/*TODO: remove this!!!
+ * the minimum would be to move it to piecewise.h but this would be stupid.
+ * The best would be to let 'compose' be aware of extension modes (constant, linear, polynomial..)
+ * (I think the extension modes (at start and end) should be properties of the pwsb).
+ */
+static
+void prolongateByConstants( Piecewise<SBasis> &f, double paddle_width ){
+ if ( f.size() == 0 ) return; //do we have a covention about the domain of empty pwsb?
+ f.cuts.insert( f.cuts.begin(), f.cuts.front() - paddle_width );
+ f.segs.insert( f.segs.begin(), SBasis( f.segs.front().at0() ) );
+ f.cuts.insert( f.cuts.end(), f.cuts.back() + paddle_width );
+ f.segs.insert( f.segs.end(), SBasis( f.segs.back().at1() ) );
+}
+
+static
+bool compareIntersectionsTimesX( SmashIntersection const &inter1, SmashIntersection const &inter2 ){
+ return inter1.times[X].min() < inter2.times[Y].min();
+}
+/*Fuse contiguous intersection domains
+ *
+ */
+static
+void cleanup_and_fuse( std::vector<SmashIntersection> &inters ){
+ std::sort( inters.begin(), inters.end(), compareIntersectionsTimesX);
+ for (unsigned i=0; i < inters.size(); i++ ){
+ for (unsigned j=i+1; j < inters.size() && inters[i].times[X].intersects( inters[j].times[X]) ; j++ ){
+ if (inters[i].times[Y].intersects( inters[j].times[Y] ) ){
+ inters[i].times.unionWith(inters[j].times);
+ inters[i].bbox.unionWith(inters[j].bbox);
+ inters.erase( inters.begin() + j );
+ }
+ }
+ }
+}
+
+/* Computes the intersection of two sets given as (ordered) union intervals.
+ */
+static
+std::vector<Interval> intersect( std::vector<Interval> const &a, std::vector<Interval> const &b){
+ std::vector<Interval> result;
+ //TODO: use order to optimize this!
+ for (auto i : a){
+ for (auto j : b){
+ OptInterval c( i );
+ c &= j;
+ if ( c ) {
+ result.push_back( *c );
+ }
+ }
+ }
+ return result;
+}
+
+/* Returns the intervals over which the curves are in the
+ * tol-neighborhood one of the other for the L_\infty metric.
+ * WARNING: each curve is supposed to be a graph over x or y axis
+ * (but not necessarily the same axis for both) and the smaller
+ * the slope the better (typically <=45°).
+ */
+std::vector<SmashIntersection> monotonic_smash_intersect( D2<SBasis> const &a, D2<SBasis> const &b, double tol){
+ using std::swap;
+
+ // a and b or X and Y may have to be exchanged, so make local copies.
+ D2<SBasis> aa = a;
+ D2<SBasis> bb = b;
+ bool swapresult = false;
+ bool swapcoord = false;//debug only!
+
+ //if the (enlarged) bounding boxes don't intersect, stop.
+ OptRect abounds = bounds_fast( a );
+ OptRect bbounds = bounds_fast( b );
+ if ( !abounds || !bbounds ) return std::vector<SmashIntersection>();
+ abounds->expandBy(tol);
+ if ( !(abounds->intersects(*bbounds))){
+ return std::vector<SmashIntersection>();
+ }
+
+ //Choose the best curve to be re-parametrized by x or y values.
+ OptRect dabounds = bounds_exact(derivative(a));
+ OptRect dbbounds = bounds_exact(derivative(b));
+ if ( dbbounds->min().length() > dabounds->min().length() ){
+ aa=b;
+ bb=a;
+ swap( dabounds, dbbounds );
+ swapresult = true;
+ }
+
+ //Choose the best coordinate to use as new parameter
+ double dxmin = std::min( std::abs((*dabounds)[X].max()), std::abs((*dabounds)[X].min()) );
+ double dymin = std::min( std::abs((*dabounds)[Y].max()), std::abs((*dabounds)[Y].min()) );
+ if ( (*dabounds)[X].max()*(*dabounds)[X].min() < 0 ) dxmin=0;
+ if ( (*dabounds)[Y].max()*(*dabounds)[Y].min() < 0 ) dymin=0;
+ assert (dxmin>=0 && dymin>=0);
+
+ if (dxmin < dymin) {
+ aa = D2<SBasis>( aa[Y], aa[X] );
+ bb = D2<SBasis>( bb[Y], bb[X] );
+ swapcoord = true;
+ }
+
+ //re-parametrize aa by the value of x.
+ Interval x_range_strict( aa[X].at0(), aa[X].at1() );
+ Piecewise<SBasis> y_of_x = pw_compose_inverse(aa[Y],aa[X], 2, 1e-5);
+
+ //Compute top and bottom boundaries of the L^infty nbhd of aa.
+ Piecewise<SBasis> top_ay, bot_ay;
+ computeLinfinityNeighborhood( y_of_x, tol, top_ay, bot_ay);
+
+ Interval ax_range = top_ay.domain();//i.e. aa[X] domain ewpanded by tol.
+ std::vector<Interval> bx_in_ax_range = level_set(bb[X], ax_range );
+
+ // find times when bb is in the neighborhood of aa.
+ std::vector<Interval> tbs;
+ for (auto & i : bx_in_ax_range){
+ D2<Piecewise<SBasis> > bb_in;
+ bb_in[X] = Piecewise<SBasis> ( portion( bb[X], i ) );
+ bb_in[Y] = Piecewise<SBasis> ( portion( bb[Y], i) );
+ bb_in[X].setDomain( i );
+ bb_in[Y].setDomain( i );
+
+ Piecewise<SBasis> h;
+ Interval level;
+ h = bb_in[Y] - compose( top_ay, bb_in[X] );
+ level = Interval( -infinity(), 0 );
+ std::vector<Interval> rts_lo = level_set( h, level);
+ h = bb_in[Y] - compose( bot_ay, bb_in[X] );
+ level = Interval( 0, infinity());
+ std::vector<Interval> rts_hi = level_set( h, level);
+
+ std::vector<Interval> rts = intersect( rts_lo, rts_hi );
+ tbs.insert(tbs.end(), rts.begin(), rts.end() );
+ }
+
+ std::vector<SmashIntersection> result(tbs.size(), SmashIntersection());
+
+ /* for each solution I, find times when aa is in the neighborhood of bb(I).
+ * (Note: the preimage of bb[X](I) by aa[X], enlarged by tol, is a good approximation of this:
+ * it would give points in the 2*tol neighborhood of bb (if the slope of aa is never more than 1).
+ * + faster computation.
+ * - implies little jumps depending on the subdivision of the input curve into monotonic pieces
+ * and on the choice of preferred axis. If noticeable, these jumps would feel random to the user :-(
+ */
+ for (unsigned j=0; j<tbs.size(); j++){
+ result[j].times[Y] = tbs[j];
+ std::vector<Interval> tas;
+ //TODO: replace this by some option in the "compose(pw,pw)" method!
+ Piecewise<SBasis> fat_y_of_x = y_of_x;
+ prolongateByConstants( fat_y_of_x, 100*(1+tol) );
+
+ D2<Piecewise<SBasis> > top_b, bot_b;
+ D2<SBasis> bbj = portion( bb, tbs[j] );
+ computeLinfinityNeighborhood( bbj, tol, top_b, bot_b );
+
+ Piecewise<SBasis> h;
+ Interval level;
+ h = top_b[Y] - compose( fat_y_of_x, top_b[X] );
+ level = Interval( +infinity(), 0 );
+ std::vector<Interval> rts_top = level_set( h, level);
+ for (auto & idx : rts_top){
+ idx = Interval( top_b[X].valueAt( idx.min() ),
+ top_b[X].valueAt( idx.max() ) );
+ }
+ assert( rts_top.size() == 1 );
+
+ h = bot_b[Y] - compose( fat_y_of_x, bot_b[X] );
+ level = Interval( 0, -infinity());
+ std::vector<Interval> rts_bot = level_set( h, level);
+ for (auto & idx : rts_bot){
+ idx = Interval( bot_b[X].valueAt( idx.min() ),
+ bot_b[X].valueAt( idx.max() ) );
+ }
+ assert( rts_bot.size() == 1 );
+ rts_top = intersect( rts_top, rts_bot );
+ assert (rts_top.size() == 1);
+ Interval x_dom = rts_top[0];
+
+ if ( x_dom.max() <= x_range_strict.min() ){
+ tas.push_back( Interval ( ( aa[X].at0() < aa[X].at1() ) ? 0 : 1 ) );
+ }else if ( x_dom.min() >= x_range_strict.max() ){
+ tas.push_back( Interval ( ( aa[X].at0() < aa[X].at1() ) ? 1 : 0 ) );
+ }else{
+ tas = level_set(aa[X], x_dom );
+ }
+ assert( tas.size()==1 );
+ result[j].times[X] = tas.front();
+
+ result[j].bbox = Rect( bbj.at0(), bbj.at1() );
+ Interval y_dom( aa[Y](result[j].times[X].min()), aa[Y](result[j].times[X].max()) );
+ result[j].bbox.unionWith( Rect( x_dom, y_dom ) );
+ }
+
+ if (swapresult) {
+ for (auto & i : result){
+ swap( i.times[X], i.times[Y]);
+ }
+ }
+ if (swapcoord) {
+ for (auto & i : result){
+ swap( i.bbox[X], i.bbox[Y] );
+ }
+ }
+
+ //TODO: cleanup result? fuse contiguous intersections?
+ return result;
+}
+
+std::vector<SmashIntersection> smash_intersect( D2<SBasis> const &a, D2<SBasis> const &b, double tol){
+ std::vector<SmashIntersection> result;
+
+ std::vector<Interval> acuts = monotonicSplit(a);
+ std::vector<Interval> bcuts = monotonicSplit(b);
+ for (auto & acut : acuts){
+ D2<SBasis> ai = portion( a, acut);
+ for (auto & bcut : bcuts){
+ D2<SBasis> bj = portion( b, bcut);
+ std::vector<SmashIntersection> ai_cap_bj = monotonic_smash_intersect( ai, bj, tol );
+ for (auto & k : ai_cap_bj){
+ k.times[X] = k.times[X] * acut.extent() + acut.min();
+ k.times[Y] = k.times[Y] * bcut.extent() + bcut.min();
+ }
+ result.insert( result.end(), ai_cap_bj.begin(), ai_cap_bj.end() );
+ }
+ }
+ cleanup_and_fuse( result );
+ return result;
+}
+
+}
+
+/*
+ Local Variables:
+ mode:c++
+ c-file-style:"stroustrup"
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+ indent-tabs-mode:nil
+ fill-column:99
+ End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
diff --git a/src/2geom/orphan-code/nearestpoint.cpp b/src/2geom/orphan-code/nearestpoint.cpp
new file mode 100644
index 0000000..870ed09
--- /dev/null
+++ b/src/2geom/orphan-code/nearestpoint.cpp
@@ -0,0 +1,405 @@
+/*
+** vim: ts=4 sw=4 et tw=0 wm=0
+**
+** RCS Information:
+** $Author: mjw $
+** $Revision: 1 $
+** $Date: 2006-03-28 15:59:38 +1100 (Tue, 28 Mar 2006) $
+**
+** Solving the Nearest Point-on-Curve Problem and
+** A Bezier Curve-Based Root-Finder
+** by Philip J. Schneider
+** from "Graphics Gems", Academic Press, 1990
+** modified by mwybrow, njh
+*/
+
+/* point_on_curve.c */
+
+static double SquaredLength(const Geom::Point a)
+{
+ return dot(a, a);
+}
+
+
+/*
+ * Forward declarations
+ */
+static int FindRoots(Geom::Point *w, int degree, double *t, int depth);
+static Geom::Point *ConvertToBezierForm( Geom::Point P, Geom::Point *V);
+static double ComputeXIntercept( Geom::Point *V, int degree);
+static int ControlPolygonFlatEnough( Geom::Point *V, int degree);
+static int CrossingCount(Geom::Point *V, int degree);
+static Geom::Point Bez(Geom::Point *V, int degree, double t, Geom::Point *Left,
+ Geom::Point *Right);
+
+int MAXDEPTH = 64; /* Maximum depth for recursion */
+
+#define EPSILON (ldexp(1.0,-MAXDEPTH-1)) /*Flatness control value */
+#define DEGREE 3 /* Cubic Bezier curve */
+#define W_DEGREE 5 /* Degree of eqn to find roots of */
+
+
+/*
+ * NearestPointOnCurve :
+ * Compute the parameter value of the point on a Bezier
+ * curve segment closest to some arbtitrary, user-input point.
+ * Return the point on the curve at that parameter value.
+ *
+ Geom::Point P; The user-supplied point
+ Geom::Point *V; Control points of cubic Bezier
+*/
+double NearestPointOnCurve(Geom::Point P, Geom::Point *V)
+{
+ double t_candidate[W_DEGREE]; /* Possible roots */
+
+ /* Convert problem to 5th-degree Bezier form */
+ Geom::Point *w = ConvertToBezierForm(P, V);
+
+ /* Find all possible roots of 5th-degree equation */
+ int n_solutions = FindRoots(w, W_DEGREE, t_candidate, 0);
+ std::free((char *)w);
+
+ /* Check distance to end of the curve, where t = 1 */
+ double dist = SquaredLength(P - V[DEGREE]);
+ double t = 1.0;
+
+ /* Find distances for candidate points */
+ for (int i = 0; i < n_solutions; i++) {
+ Geom::Point p = Bez(V, DEGREE, t_candidate[i], NULL, NULL);
+ double new_dist = SquaredLength(P - p);
+ if (new_dist < dist) {
+ dist = new_dist;
+ t = t_candidate[i];
+ }
+ }
+
+ /* Return the parameter value t */
+ return t;
+}
+
+
+/*
+ * ConvertToBezierForm :
+ * Given a point and a Bezier curve, generate a 5th-degree
+ * Bezier-format equation whose solution finds the point on the
+ * curve nearest the user-defined point.
+ */
+static Geom::Point *ConvertToBezierForm(
+ Geom::Point P, /* The point to find t for */
+ Geom::Point *V) /* The control points */
+{
+ Geom::Point c[DEGREE+1]; /* V(i)'s - P */
+ Geom::Point d[DEGREE]; /* V(i+1) - V(i) */
+ Geom::Point *w; /* Ctl pts of 5th-degree curve */
+ double cdTable[3][4]; /* Dot product of c, d */
+ static double z[3][4] = { /* Precomputed "z" for cubics */
+ {1.0, 0.6, 0.3, 0.1},
+ {0.4, 0.6, 0.6, 0.4},
+ {0.1, 0.3, 0.6, 1.0},
+ };
+
+
+ /*Determine the c's -- these are vectors created by subtracting*/
+ /* point P from each of the control points */
+ for (int i = 0; i <= DEGREE; i++) {
+ c[i] = V[i] - P;
+ }
+ /* Determine the d's -- these are vectors created by subtracting*/
+ /* each control point from the next */
+ for (int i = 0; i <= DEGREE - 1; i++) {
+ d[i] = 3.0*(V[i+1] - V[i]);
+ }
+
+ /* Create the c,d table -- this is a table of dot products of the */
+ /* c's and d's */
+ for (int row = 0; row <= DEGREE - 1; row++) {
+ for (int column = 0; column <= DEGREE; column++) {
+ cdTable[row][column] = dot(d[row], c[column]);
+ }
+ }
+
+ /* Now, apply the z's to the dot products, on the skew diagonal*/
+ /* Also, set up the x-values, making these "points" */
+ w = (Geom::Point *)malloc((unsigned)(W_DEGREE+1) * sizeof(Geom::Point));
+ for (int i = 0; i <= W_DEGREE; i++) {
+ w[i][Geom::Y] = 0.0;
+ w[i][Geom::X] = (double)(i) / W_DEGREE;
+ }
+
+ const int n = DEGREE;
+ const int m = DEGREE-1;
+ for (int k = 0; k <= n + m; k++) {
+ const int lb = std::max(0, k - m);
+ const int ub = std::min(k, n);
+ for (int i = lb; i <= ub; i++) {
+ int j = k - i;
+ w[i+j][Geom::Y] += cdTable[j][i] * z[j][i];
+ }
+ }
+
+ return w;
+}
+
+
+/*
+ * FindRoots :
+ * Given a 5th-degree equation in Bernstein-Bezier form, find
+ * all of the roots in the interval [0, 1]. Return the number
+ * of roots found.
+ */
+static int FindRoots(
+ Geom::Point *w, /* The control points */
+ int degree, /* The degree of the polynomial */
+ double *t, /* RETURN candidate t-values */
+ int depth) /* The depth of the recursion */
+{
+ int i;
+ Geom::Point Left[W_DEGREE+1], /* New left and right */
+ Right[W_DEGREE+1]; /* control polygons */
+ int left_count, /* Solution count from */
+ right_count; /* children */
+ double left_t[W_DEGREE+1], /* Solutions from kids */
+ right_t[W_DEGREE+1];
+
+ switch (CrossingCount(w, degree)) {
+ case 0 : { /* No solutions here */
+ return 0;
+ break;
+ }
+ case 1 : { /* Unique solution */
+ /* Stop recursion when the tree is deep enough */
+ /* if deep enough, return 1 solution at midpoint */
+ if (depth >= MAXDEPTH) {
+ t[0] = (w[0][Geom::X] + w[W_DEGREE][Geom::X]) / 2.0;
+ return 1;
+ }
+ if (ControlPolygonFlatEnough(w, degree)) {
+ t[0] = ComputeXIntercept(w, degree);
+ return 1;
+ }
+ break;
+ }
+ }
+
+ /* Otherwise, solve recursively after */
+ /* subdividing control polygon */
+ Bez(w, degree, 0.5, Left, Right);
+ left_count = FindRoots(Left, degree, left_t, depth+1);
+ right_count = FindRoots(Right, degree, right_t, depth+1);
+
+
+ /* Gather solutions together */
+ for (i = 0; i < left_count; i++) {
+ t[i] = left_t[i];
+ }
+ for (i = 0; i < right_count; i++) {
+ t[i+left_count] = right_t[i];
+ }
+
+ /* Send back total number of solutions */
+ return (left_count+right_count);
+}
+
+
+/*
+ * CrossingCount :
+ * Count the number of times a Bezier control polygon
+ * crosses the 0-axis. This number is >= the number of roots.
+ *
+ */
+static int CrossingCount(
+ Geom::Point *V, /* Control pts of Bezier curve */
+ int degree) /* Degree of Bezier curve */
+{
+ int n_crossings = 0; /* Number of zero-crossings */
+ int old_sign; /* Sign of coefficients */
+
+ old_sign = Geom::sgn(V[0][Geom::Y]);
+ for (int i = 1; i <= degree; i++) {
+ int sign = Geom::sgn(V[i][Geom::Y]);
+ if (sign != old_sign)
+ n_crossings++;
+ old_sign = sign;
+ }
+ return n_crossings;
+}
+
+
+
+/*
+ * ControlPolygonFlatEnough :
+ * Check if the control polygon of a Bezier curve is flat enough
+ * for recursive subdivision to bottom out.
+ *
+ */
+static int ControlPolygonFlatEnough(
+ Geom::Point *V, /* Control points */
+ int degree) /* Degree of polynomial */
+{
+ int i; /* Index variable */
+ double *distance; /* Distances from pts to line */
+ double max_distance_above; /* maximum of these */
+ double max_distance_below;
+ double error; /* Precision of root */
+ //Geom::Point t; /* Vector from V[0] to V[degree]*/
+ double intercept_1,
+ intercept_2,
+ left_intercept,
+ right_intercept;
+ double a, b, c; /* Coefficients of implicit */
+ /* eqn for line from V[0]-V[deg]*/
+
+ /* Find the perpendicular distance */
+ /* from each interior control point to */
+ /* line connecting V[0] and V[degree] */
+ distance = (double *)malloc((unsigned)(degree + 1) * sizeof(double));
+ {
+ double abSquared;
+
+ /* Derive the implicit equation for line connecting first */
+ /* and last control points */
+ a = V[0][Geom::Y] - V[degree][Geom::Y];
+ b = V[degree][Geom::X] - V[0][Geom::X];
+ c = V[0][Geom::X] * V[degree][Geom::Y] - V[degree][Geom::X] * V[0][Geom::Y];
+
+ abSquared = (a * a) + (b * b);
+
+ for (i = 1; i < degree; i++) {
+ /* Compute distance from each of the points to that line */
+ distance[i] = a * V[i][Geom::X] + b * V[i][Geom::Y] + c;
+ if (distance[i] > 0.0) {
+ distance[i] = (distance[i] * distance[i]) / abSquared;
+ }
+ if (distance[i] < 0.0) {
+ distance[i] = -((distance[i] * distance[i]) / abSquared);
+ }
+ }
+ }
+
+
+ /* Find the largest distance */
+ max_distance_above = 0.0;
+ max_distance_below = 0.0;
+ for (i = 1; i < degree; i++) {
+ if (distance[i] < 0.0) {
+ max_distance_below = std::min(max_distance_below, distance[i]);
+ };
+ if (distance[i] > 0.0) {
+ max_distance_above = std::max(max_distance_above, distance[i]);
+ }
+ }
+ free((char *)distance);
+
+ {
+ double det;
+ double a1, b1, c1, a2, b2, c2;
+
+ /* Implicit equation for zero line */
+ a1 = 0.0;
+ b1 = 1.0;
+ c1 = 0.0;
+
+ /* Implicit equation for "above" line */
+ a2 = a;
+ b2 = b;
+ c2 = c + max_distance_above;
+
+ det = a1 * b2 - a2 * b1;
+
+ intercept_1 = (b1 * c2 - b2 * c1) / det;
+
+ /* Implicit equation for "below" line */
+ a2 = a;
+ b2 = b;
+ c2 = c + max_distance_below;
+
+ det = a1 * b2 - a2 * b1;
+
+ intercept_2 = (b1 * c2 - b2 * c1) / det;
+ }
+
+ /* Compute intercepts of bounding box */
+ left_intercept = std::min(intercept_1, intercept_2);
+ right_intercept = std::max(intercept_1, intercept_2);
+
+ error = 0.5 * (right_intercept-left_intercept);
+ if (error < EPSILON) {
+ return 1;
+ }
+ else {
+ return 0;
+ }
+}
+
+
+
+/*
+ * ComputeXIntercept :
+ * Compute intersection of chord from first control point to last
+ * with 0-axis.
+ *
+ */
+static double ComputeXIntercept(
+ Geom::Point *V, /* Control points */
+ int degree) /* Degree of curve */
+{
+ const Geom::Point A = V[degree] - V[0];
+
+ return (A[Geom::X]*V[0][Geom::Y] - A[Geom::Y]*V[0][Geom::X]) / -A[Geom::Y];
+}
+
+
+/*
+ * Bez :
+ * Evaluate a Bezier curve at a particular parameter value
+ * Fill in control points for resulting sub-curves if "Left" and
+ * "Right" are non-null.
+ *
+ */
+static Geom::Point Bez(
+ Geom::Point *V, /* Control pts */
+ int degree, /* Degree of bezier curve */
+ double t, /* Parameter value */
+ Geom::Point *Left, /* RETURN left half ctl pts */
+ Geom::Point *Right) /* RETURN right half ctl pts */
+{
+ Geom::Point Vtemp[W_DEGREE+1][W_DEGREE+1];
+
+
+ /* Copy control points */
+ for (int j =0; j <= degree; j++) {
+ Vtemp[0][j] = V[j];
+ }
+
+ /* Triangle computation */
+ for (int i = 1; i <= degree; i++) {
+ for (int j =0 ; j <= degree - i; j++) {
+ Vtemp[i][j] =
+ (1.0 - t) * Vtemp[i-1][j] + t * Vtemp[i-1][j+1];
+ }
+ }
+
+ if (Left != NULL) {
+ for (int j = 0; j <= degree; j++) {
+ Left[j] = Vtemp[j][0];
+ }
+ }
+ if (Right != NULL) {
+ for (int j = 0; j <= degree; j++) {
+ Right[j] = Vtemp[degree-j][j];
+ }
+ }
+
+ return (Vtemp[degree][0]);
+}
+
+/*
+ Local Variables:
+ mode:c++
+ c-file-style:"stroustrup"
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+ indent-tabs-mode:nil
+ fill-column:99
+ End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
diff --git a/src/2geom/orphan-code/redblack-toy.cpp b/src/2geom/orphan-code/redblack-toy.cpp
new file mode 100644
index 0000000..01ffa7d
--- /dev/null
+++ b/src/2geom/orphan-code/redblack-toy.cpp
@@ -0,0 +1,327 @@
+/*
+ * Copyright 2009 Evangelos Katsikaros <vkatsikaros at yahoo dot gr>
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+/*
+ initial toy for redblack trees
+*/
+
+#include <2geom/toys/path-cairo.h>
+#include <2geom/toys/toy-framework-2.h>
+
+#include <2geom/orphan-code/redblacktree.h>
+#include <2geom/orphan-code/redblacktree.cpp>
+
+#include <time.h>
+using std::vector;
+using namespace Geom;
+using namespace std;
+
+
+class RedBlackToy: public Toy
+{
+ PointSetHandle handle_set;
+ Geom::Point starting_point; // during click and drag: start point of click
+ Geom::Point ending_point; // during click and drag: end point of click (release)
+ Geom::Point highlight_point; // not used
+
+ Geom::RedBlackTree rbtree_x;
+ RedBlack* search_result;
+ RedBlack temp_deleted_node;
+
+ // colors we are going to use for different purposes
+ colour color_rect, color_rect_guide; // black(a=0.6), black
+ colour color_select_area, color_select_area_guide; // red(a=0.6), red
+
+ int alter_existing_rect;
+ int add_new_rect;
+
+ Rect rect_chosen; // the rectangle of the search area
+ Rect dummy_draw; // the "helper" rectangle that is shown during the click and drag (before the mouse release)
+ int mode; // insert/alter, search, delete modes
+
+ // printing of the tree
+ int help_counter; // the "x" of the label of each node
+ static const int label_size = 15 ; // size the label of each node
+
+ // used for the keys that switch between modes
+ enum menu_item_t
+ {
+ INSERT = 0,
+ DELETE,
+ SEARCH,
+ TOTAL_ITEMS // this one must be the last item
+ };
+ static const char* menu_items[TOTAL_ITEMS];
+ static const char keys[TOTAL_ITEMS];
+
+
+
+ void draw(cairo_t *cr, std::ostringstream *notify, int width, int height, bool save, std::ostringstream *timer_stream) {
+ cairo_set_line_width( cr, 1 );
+
+ // draw the rects that we have in the handles
+ for( unsigned i=0; i<handle_set.pts.size(); i=i+2 ){
+ Rect r1( handle_set.pts[i], handle_set.pts[i+1] );
+ cairo_rectangle( cr, r1 );
+ }
+ cairo_set_source_rgba( cr, color_rect);
+ cairo_stroke( cr );
+
+ // draw a rect if we click & drag (so that we know what we are going to create)
+ if(add_new_rect){
+ dummy_draw = Rect( starting_point, ending_point );
+ cairo_rectangle( cr, dummy_draw );
+ if( mode == 0){
+ cairo_set_source_rgba( cr, color_rect_guide);
+ }
+ else if( mode == 1){
+ cairo_set_source_rgba( cr, color_select_area_guide );
+ }
+ cairo_stroke( cr );
+ }
+
+ // draw a rect for the search area
+ cairo_rectangle( cr, rect_chosen );
+ cairo_set_source_rgba( cr, color_select_area);
+ cairo_stroke( cr );
+
+ Toy::draw( cr, notify, width, height, save,timer_stream );
+ draw_tree_in_toy( cr ,rbtree_x.root, 0);
+ help_counter=0;
+ }
+
+ void mouse_moved(GdkEventMotion* e){
+ if( !( alter_existing_rect && mode == 1 ) ){
+ Toy::mouse_moved(e);
+ }
+
+ if(add_new_rect){
+ ending_point = Point(e->x, e->y);
+ }
+ }
+
+ void mouse_pressed(GdkEventButton* e) {
+ Toy::mouse_pressed(e);
+ if(e->button == 1){ // left mouse button
+ if( mode == 0 ){ // mode: insert / alter
+ if(!selected) {
+ starting_point = Point(e->x, e->y);
+ ending_point = starting_point;
+ add_new_rect = 1;
+ }
+ else
+ {
+ // TODO find the selected rect
+ // ideas : from Handle *selected ???
+ //std::cout <<find_selected_rect(selected) << std::endl ;
+ alter_existing_rect = 1;
+ }
+ }
+ else if( mode == 1 ){ // mode: search
+ if(!selected) {
+ starting_point = Point(e->x, e->y);
+ ending_point = starting_point;
+ add_new_rect = 1;
+ }
+ else{
+ alter_existing_rect = 1;
+ }
+ }
+ else if( mode == 2) { // mode: delete
+ }
+ }
+ else if(e->button == 2){ //middle button
+ }
+ else if(e->button == 3){ //right button
+ }
+ }
+
+ virtual void mouse_released(GdkEventButton* e) {
+ Toy::mouse_released(e);
+ if( e->button == 1 ) { //left mouse button
+ if( mode == 0) { // mode: insert / alter
+ if( add_new_rect ){
+ ending_point = Point(e->x, e->y);
+ handle_set.push_back(starting_point);
+ handle_set.push_back(ending_point);
+ insert_in_tree_the_last_rect();
+ add_new_rect = 0;
+ }
+ else if( alter_existing_rect ){
+ //TODO update rect (and tree)
+ // delete selected rect
+ // insert altered
+ alter_existing_rect = 0;
+ }
+ }
+ else if( mode == 1 ){ // mode: search
+ if( add_new_rect ){
+ ending_point = Point(e->x, e->y);
+ rect_chosen = Rect(starting_point, ending_point);
+
+ // search in the X axis
+ Coord a = rect_chosen[0].min();
+ Coord b = rect_chosen[0].max();
+ search_result = rbtree_x.search( Interval( a, b ) );
+ if(search_result){
+ std::cout << "Found: (" << search_result->data << ": " << search_result->key()
+ << ", " << search_result->high() << " : " << search_result->subtree_max << ") "
+ << std::endl;
+ }
+ else{
+ std::cout << "Nothing found..."<< std::endl;
+ }
+ add_new_rect = 0;
+ }
+ else if(alter_existing_rect){
+ // do nothing
+ alter_existing_rect = 0;
+ }
+ }
+ else if( mode == 2) { // mode: delete
+
+ }
+ }
+ else if(e->button == 2){ //middle button
+ }
+ else if(e->button == 3){ //right button
+
+ }
+ }
+
+
+ void key_hit(GdkEventKey *e)
+ {
+ char choice = std::toupper(e->keyval);
+ switch ( choice )
+ {
+ case 'A':
+ mode = 0;
+ break;
+ case 'B':
+ mode = 1;
+ break;
+ case 'C':
+ mode = 2;
+ break;
+ }
+ //redraw();
+ }
+
+ void insert_in_tree_the_last_rect(){
+ unsigned i = handle_set.pts.size() - 2;
+ Rect r1(handle_set.pts[i], handle_set.pts[i+1]);
+ // insert in X axis tree
+ rbtree_x.insert(r1, i, 0);
+ rbtree_x.print_tree();
+ };
+
+ void draw_tree_in_toy(cairo_t* cr, Geom::RedBlack* n, int depth = 0) {
+ if(n){
+ if(n->left){
+ draw_tree_in_toy(cr, n->left, depth+1);
+ }
+ help_counter += 1;
+ //drawthisnode(cr, x*10, depth*10);
+ if(n->isRed){
+ cairo_set_source_rgba (cr, color_select_area_guide);
+ }
+ else{
+ cairo_set_source_rgba (cr, color_rect_guide);
+ }
+
+ cairo_stroke(cr);
+
+ Geom::Point text_point = Point( help_counter*15, depth*15 );
+ char label[4];
+ sprintf( label,"%d",n->data ); // instead of std::itoa(depth, label, 10);
+
+ draw_text(cr, text_point, label);
+ ////////////////////////////////////////////////////////////////
+ if(n->right){
+ draw_tree_in_toy(cr, n->right, depth+1);
+ }
+ }
+ };
+
+/*
+ int find_selected_rect(PointHandle * selected){
+
+ for( unsigned i=0; i<handle_set.pts.size(); i=i+2 ){
+ if( handle_set.pts[i] == selected || handle_set.pts[i+1] == selected ){
+ return i;
+ }
+ }
+
+ return -1;
+ };
+*/
+
+
+public:
+ RedBlackToy(): color_rect(0, 0, 0, 0.6), color_rect_guide(0, 0, 0, 1),
+ color_select_area(1, 0, 0, 0.6 ), color_select_area_guide(1, 0, 0, 1 ),
+ alter_existing_rect(0), add_new_rect(0), mode(0), help_counter(0)
+ {
+ if(handles.empty()) {
+ handles.push_back(&handle_set);
+ }
+ Rect rect_chosen();
+ Rect dummy_draw();
+ }
+
+
+};
+
+
+
+int main(int argc, char **argv) {
+ std::cout << "---------------------------------------------------------"<< std::endl;
+ std::cout << "Let's play with the Red Black Tree! ONLY Insert works now!!!"<< std::endl;
+ std::cout << " Key A: insert/alter mode "<< std::endl;
+ std::cout << " * Left click and drag on white area: create a rectangle"<< std::endl;
+ std::cout << " *NOT READY: Left click and drag on handler: alter a rectangle"<< std::endl;
+ std::cout << " Key B: search mode "<< std::endl;
+ std::cout << " * Left click and drag on white area: \"search\" for nodes that intersect red area"<< std::endl;
+ std::cout << " NOT READY: Key C: delete mode "<< std::endl;
+ std::cout << " * Left click on handler: delete for a rectangle"<< std::endl;
+ std::cout << "---------------------------------------------------------"<< std::endl;
+ init(argc, argv, new RedBlackToy);
+ return 0;
+}
+
+const char* RedBlackToy::menu_items[] =
+{
+ "Insert / Alter Rectangle",
+ "Search Rectangle",
+ "Delete Reactangle"
+};
+
+const char RedBlackToy::keys[] =
+{
+ 'A', 'B', 'C'
+};
diff --git a/src/2geom/orphan-code/redblacktree.cpp b/src/2geom/orphan-code/redblacktree.cpp
new file mode 100644
index 0000000..bf9a728
--- /dev/null
+++ b/src/2geom/orphan-code/redblacktree.cpp
@@ -0,0 +1,575 @@
+#include <2geom/orphan-code/redblacktree.h>
+//#include <algorithm>
+
+
+#define _REDBLACK_PRINT(x) std::cout << x << std::endl;
+//comment the following if you want output during RedBlack Tree operations
+//#define _REDBLACK_PRINT(x) ;
+
+
+namespace Geom{
+
+
+
+RedBlack* RedBlackTree::search(Rect const &r, int dimension){
+ return search( Interval( r[dimension].min(), r[dimension].max() ) );
+ // TODO get rid of dimension
+ // TODO put 2 trees (X, Y axis in one lump)
+}
+
+/*
+INTERVAL-SEARCH(T, i)
+1 x <- root[T]
+2 while x != nil[T] and i does not overlap int[x]
+3 do if left[x] != nil[T] and max[left[x]] >= low[i]
+4 then x <- left[x]
+5 else x <- right[x]
+6 return x
+
+Two intervals i,x overlap in the 4 following cases:
+ 1) |--------| i
+ |---| x
+
+ 2) |-----| i
+ |----------| x
+
+ 3) |------| i
+ |------| x
+
+ 4) |----| i
+ |----| x
+
+And do not overlap when (the one os left or right of the other)
+ 1) |--------| i
+ |---| x
+
+ 2) |-----| i
+ |----------| x
+
+
+*/
+RedBlack* RedBlackTree::search(Interval i){
+ _REDBLACK_PRINT( "==============================================================" << std::endl << "ENTER: search(Interval i) : (" << i.min() << ", " << i.max() << ")" )
+ RedBlack *x;
+ x = root;
+
+ while( x!=0 &&
+ ( i.max() < x->interval.min() ||
+ i.min() > x->interval.max() )
+ ){
+ _REDBLACK_PRINT( "(" << x->data << ": " << x->key() << ", " << x->high() << " : " << x->subtree_max << ") "
+ << " i do not overlap with x")
+
+ if(x->left != 0 && (x->left)->subtree_max >= i.min() ){
+ x = x->left;
+ }
+ else{
+ x = x->right;
+ }
+ }
+ _REDBLACK_PRINT( "RETURN: search" << std::endl )
+ return x;
+}
+
+
+
+void RedBlackTree::insert(Rect const &r, int shape, int dimension) {
+ _REDBLACK_PRINT( "==============================================================" << std::endl << "ENTER: insert(Rect, int, dimension): " << " dimension:" << dimension << " shape:" << shape )
+ insert(r[dimension].min(), r[dimension].max(), shape);
+ _REDBLACK_PRINT( "RETURN: insert(Rect, int, dimension)")
+}
+
+// source: book pp 251
+void RedBlackTree::insert(Coord dimension_min, Coord dimension_max, int shape) {
+ _REDBLACK_PRINT( std::endl << "ENTER: insert(Coord, Coord, int): " << dimension_min << ", " << dimension_max << " , shape: " << shape )
+ // x is the new node we insert
+ RedBlack *x = new RedBlack();
+ x->interval = Interval( dimension_min, dimension_max );
+ x->data = shape;
+ x->isRed = true;
+
+ _REDBLACK_PRINT( " x: " << x << " KEY: " << x->key() << " high: " << x->high() )
+
+ tree_insert(x);
+
+ print_tree();
+
+ _REDBLACK_PRINT( " Begin coloring" )
+ // we now do the coloring of the tree.
+ _REDBLACK_PRINT( " while( x!= root && (x->parent)->isRed )" )
+ while( x!= root && (x->parent)->isRed ){
+ _REDBLACK_PRINT( " ((x->parent)->parent)->left:" << ((x->parent)->parent)->left << " ((x->parent)->parent)->right:" << ((x->parent)->parent)->right )
+
+ if( x->parent == ((x->parent)->parent)->left ){
+ _REDBLACK_PRINT( " Left:" )
+ RedBlack *y = new RedBlack;
+ y = ((x->parent)->parent)->right;
+ if( y == 0 ){
+ /*
+ This 1st brach is not in the book, but is needed. We must check y->isRed but it is
+ undefined, so we get segfault. But 0 (undefined) means that y is a leaf, so it is
+ black by definition. So, do the same as in the else part.
+ */
+ _REDBLACK_PRINT( " y==0" )
+ if( x == (x->parent)->right ){
+ x = x->parent;
+ left_rotate(x);
+ }
+ (x->parent)->isRed = false;
+ ((x->parent)->parent)->isRed = true;
+ right_rotate((x->parent)->parent);
+ }
+ else if( y->isRed ){
+ _REDBLACK_PRINT( " y->isRed" )
+ (x->parent)->isRed = false;
+ y->isRed = false;
+ ((x->parent)->parent)->isRed = true;
+ x = (x->parent)->parent;
+ }
+ else{
+ _REDBLACK_PRINT( " !( y->isRed)" )
+ if( x == (x->parent)->right ){
+ x = x->parent;
+ left_rotate(x);
+ }
+ (x->parent)->isRed = false;
+ ((x->parent)->parent)->isRed = true;
+ right_rotate((x->parent)->parent);
+ }
+ }
+ else{ // this branch is the same with the above if clause with "right", "left" exchanged
+ _REDBLACK_PRINT( " Right:" )
+ RedBlack *y = new RedBlack;
+ y = ((x->parent)->parent)->left;
+ if( y == 0 ){
+ /*
+ This 1st brach is not in the book, but is needed. We must check y->isRed but it is
+ undefined, so we get segfault. But 0 (undefined) means that y is a leaf, so it is
+ black by definition. So, do the same as in the else part.
+ */
+ _REDBLACK_PRINT( " y==0" )
+ if( x == (x->parent)->left ){
+ x = x->parent;
+ right_rotate(x);
+ }
+ (x->parent)->isRed = false;
+ ((x->parent)->parent)->isRed = true;
+ left_rotate((x->parent)->parent);
+ }
+ else if( y->isRed ){
+ _REDBLACK_PRINT( " y->isRed" )
+ (x->parent)->isRed = false;
+ y->isRed = false;
+ ((x->parent)->parent)->isRed = true;
+ x = (x->parent)->parent;
+ }
+ else{
+ _REDBLACK_PRINT( " !( y->isRed)" )
+ if( x == (x->parent)->left ){
+ x = x->parent;
+ right_rotate(x);
+ }
+ (x->parent)->isRed = false;
+ ((x->parent)->parent)->isRed = true;
+ left_rotate((x->parent)->parent);
+ }
+ }
+ }
+ root->isRed = false;
+
+ // update the max value with a slow/stupid yet certain way, walking all the tree :P
+ // TODO find better way
+ _REDBLACK_PRINT( " Update max" )
+
+ update_max(root);
+
+ _REDBLACK_PRINT( "RETURN: insert(Coord, Coord, int)" << std::endl)
+}
+
+// from book p. 266)
+void RedBlackTree::left_rotate(RedBlack* x){
+ // x->right != 0 (assumption book page 266)
+ // ??? hm problem ???
+ _REDBLACK_PRINT( "ENTER: left_rotate" )
+ RedBlack* y = new RedBlack;
+ y = x->right;
+ x->right = y->left;
+
+ if( y->left != 0){
+ (y->left)->parent = x;
+ }
+
+ y->parent = x->parent;
+
+ if( x->parent == 0){
+ root = y;
+ }
+ else{
+ if( x == (x->parent)->left ){
+ (x->parent)->left = y;
+ }
+ else{
+ (x->parent)->right = y;
+ }
+ }
+ y->left = x;
+ x->parent = y;
+ _REDBLACK_PRINT( "RETURN: left_rotate" << std::endl )
+}
+
+// from book p. 266: right_rotate is inverse of left_rotate
+// same to left_rotate with "right", "left" exchanged
+void RedBlackTree::right_rotate(RedBlack* x){
+ // x->right != 0 (assumption book page 266)
+ // ??? hm problem ??
+ _REDBLACK_PRINT( "ENTER: right_rotate" )
+ RedBlack* y = new RedBlack;
+
+ _REDBLACK_PRINT( "x->left: " << x->left )
+ y = x->left;
+ x->left = y->right;
+
+ if( y->right != 0){
+ (y->right)->parent = x;
+ }
+
+ y->parent = x->parent;
+
+ if( x->parent == 0){
+ root = y;
+ }
+ else{
+ if( x == (x->parent)->left ){
+ (x->parent)->left = y;
+ }
+ else{
+ (x->parent)->right = y;
+ }
+ }
+ y->right = x;
+ x->parent = y;
+ _REDBLACK_PRINT( "RETURN: right_rotate" << std::endl )
+}
+
+// insertion in binary search tree: book page 251
+// then the redblack insert performs the coloring
+void RedBlackTree::tree_insert(RedBlack* z){
+ _REDBLACK_PRINT( "ENTER: tree_insert(RedBlack* z)" )
+ RedBlack* y = 0; // y <- nil
+
+ RedBlack* x = root;
+
+ _REDBLACK_PRINT( " while x!=0 " )
+ while( x != 0 ){
+ y = x;
+// _REDBLACK_PRINT( " x:" << x << " y:" << y << " z:" << z )
+ _REDBLACK_PRINT( " z->key: " << z->key() << " y->key: " << y->key() << " compare")
+ if( z->key() < x->key() ){
+ _REDBLACK_PRINT( " z smaller: go left" )
+ x = x->left;
+ }
+ else{
+ _REDBLACK_PRINT( " z bigger: go right" )
+ x = x->right;
+ }
+ }
+
+ _REDBLACK_PRINT( " z->parent = y" )
+ z->parent = y;
+
+ if( y == 0 ){
+ _REDBLACK_PRINT( " set z root (empty tree)" )
+ root = z;
+ }
+ else{
+ _REDBLACK_PRINT( " z->key: " << z->key() << " y->key: " << y->key() << " compare")
+ if( z->key() < y->key() ){
+ _REDBLACK_PRINT( " z->key() smaller: y->left = z; " )
+ y->left = z;
+ }
+ else{
+ _REDBLACK_PRINT( " z->key() bigger: y->right = z " )
+ y->right = z;
+ }
+ }
+ _REDBLACK_PRINT( "RETURN: tree_insert(RedBlack* z)" << std::endl )
+}
+
+
+/*
+RB-DELETE(T, z)
+ 1 if left[z] = nil[T] or right[z] = nil[T]
+ 2 then y <- z
+ 3 else y <- TREE-SUCCESSOR(z)
+ 4 if left[y] != nil[T]
+ 5 then x <- left[y]
+ 6 else x <- right[y]
+ 7 p[x] <- p[y]
+ 8 if p[y] = nil[T]
+ 9 then root[T] <- x
+10 else if y = left[p[y]]
+11 then left[p[y]] <- x
+12 else right[p[y]] <- x
+13 if y != z
+14 then key[z] <- key[y]
+15 copy y's satellite data into z
+16 if color[y] = BLACK
+17 then RB-DELETE-FIXUP(T, x)
+18 return y
+*/
+RedBlack* RedBlackTree::erase(RedBlack* z){
+ _REDBLACK_PRINT( "==============================================================" << std::endl << "ENTER: earse(z)" )
+ RedBlack* x = new RedBlack();
+ RedBlack* y = new RedBlack();
+ if( z->left == 0 || z->right == 0 ){
+ y = z;
+ }
+ else{
+ y = tree_successor(z);
+ }
+
+ if( y->left != 0 ){
+ x = y->left;
+ }
+ else{
+ x = y->right;
+ }
+
+ x->parent = y->parent;
+
+ if( y->parent == 0){
+ root = x;
+ }
+ else {
+ if( y == (y->parent)->left ){
+ (y->parent)->left = x;
+ }
+ else{
+ (y->parent)->right = x;
+ }
+ }
+
+ if( y != z){
+ z->interval = y->interval ; // key[z] <- key[y] TODO check this
+ //copy y's satellite data into z
+ z->data = y->data;
+ z->isRed = y->isRed;
+
+ z->left = y->left;
+ z->right = y->right;
+ z->parent = y->parent;
+ }
+
+ if( y->isRed == false){
+ erase_fixup(x);
+ }
+
+ _REDBLACK_PRINT( "Update max" )
+ update_max(root);
+
+ _REDBLACK_PRINT( "RETURN: erase" )
+ return y;
+}
+
+/*
+RB-DELETE-FIXUP(T, x)
+ 1 while x != root[T] and color[x] = BLACK
+ 2 do if x = left[p[x]]
+ 3 then w <- right[p[x]]
+ 4 if color[w] = RED
+ 5 then color[w] <- BLACK Case 1
+ 6 color[p[x]] <- RED Case 1
+ 7 LEFT-ROTATE(T, p[x]) Case 1
+ 8 w <- right[p[x]]
+ 9 if color[left[w]] = BLACK and color[right[w]] = BLACK
+10 then color[w] <- RED Case 2
+11 x p[x] Case 2
+12 else if color[right[w]] = BLACK
+13 then color[left[w]] <- BLACK Case 3
+14 color[w] <- RED Case 3
+15 RIGHT-ROTATE(T, w) Case 3
+16 w <- right[p[x]] Case 3
+17 color[w] <- color[p[x]] Case 4
+18 color[p[x]] <- BLACK Case 4
+19 color[right[w]] <- BLACK Case 4
+20 LEFT-ROTATE(T, p[x]) Case 4
+21 x <- root[T] Case 4
+22 else (same as then clause with "right" and "left" exchanged)
+23 color[x] <- BLACK
+*/
+void RedBlackTree::erase_fixup(RedBlack* x){
+ RedBlack* w = 0;
+ while( x != root && x->isRed == false ){
+ if( x == (x->parent)->left ){
+ w = (x->parent)->right;
+ if(w->isRed == true){
+ w->isRed = false;
+ (w->parent)->isRed = true;
+ left_rotate(x->parent);
+ w = (x->parent)->right;
+ }
+ if( (w->left)->isRed == false && (w->right)->isRed == false ){
+ w->isRed = true;
+ x = x->parent; // TODO understand why this happens ???
+ }
+ else{
+ if( (w->right)->isRed == false ){
+ (w->left)->isRed = false;
+ right_rotate(w);
+ w = (x->parent)->right;
+ }
+ else{ // TODO ??? is this correct ???
+ w->isRed = (x->parent)->isRed;
+ (x->parent)->isRed = false;
+ (w->right)->isRed = false;
+ left_rotate(x->parent);
+ x = root; // TODO ??? is this correct ???
+ }
+ }
+ }
+ else{ // same as then clause with "right" and "left" exchanged
+ w = (x->parent)->left;
+ if(w->isRed == true){
+ w->isRed = false;
+ (w->parent)->isRed = true;
+ right_rotate(x->parent);
+ w = (x->parent)->left;
+ }
+ if( (w->right)->isRed == false && (w->left)->isRed == false ){
+ w->isRed = true;
+ x = x->parent; // ??? is this correct ???
+ }
+ else{
+ if( (w->left)->isRed == false ){
+ (w->right)->isRed = false;
+ left_rotate(w);
+ w = (x->parent)->left;
+ }
+ else{ // TODO ??? is this correct ???
+ w->isRed = (x->parent)->isRed;
+ (x->parent)->isRed = false;
+ (w->left)->isRed = false;
+ right_rotate(x->parent);
+ x = root; // TODO ??? is this correct ???
+ }
+ }
+ }
+ }
+ x->isRed = false;
+}
+
+
+void RedBlackTree::print_tree(){
+ std::cout << "Print RedBlackTree status:" << std::endl;
+ inorder_tree_walk(root);
+}
+
+
+void RedBlackTree::inorder_tree_walk(RedBlack* x){
+ int oops =0;
+ if( x != 0 ){
+ inorder_tree_walk(x->left);
+ std::cout<< "(" << x->data << ": " << x->key() << ", " << x->high() << " : " << x->subtree_max << ") " ;
+
+ if( x->left != 0 ){
+ std::cout<< "L:(" << (x->left)->data << ", " << (x->left)->key() << ") " ;
+ if( x->key() < (x->left)->key()){
+ std::cout<<" !!! ";
+ oops = 1;
+ }
+ }
+ else{
+ std::cout<< "L:0 " ;
+ }
+
+ if( x->right != 0 ){
+ std::cout<< "R:(" << (x->right)->data << ", "<< (x->right)->key() << ") " ;
+ if( x->key() > (x->right)->key() ){
+ std::cout<<" !!! ";
+ oops = 1;
+ }
+ }
+ else{
+ std::cout<< "R:0 " ;
+ }
+
+ if(oops){
+ std::cout<<" ....... !!! Problem " << oops ;
+ }
+ std::cout << std::endl;
+ inorder_tree_walk(x->right);
+ }
+}
+
+// not an norder walk of the tree
+void RedBlackTree::update_max(RedBlack* x){
+ Coord max_left, max_right;
+ if( x != 0 ){
+ update_max(x->left);
+ update_max(x->right);
+
+ // check for child
+ // if child is Nil then max = DBL_MIN
+ // could there be problems when comparing for max between two DBL_MIN ???
+ if( x->left == 0 ){
+ max_left = DBL_MIN ;
+ }
+ else{
+ max_left = (x->left)->subtree_max;
+ }
+
+ if( x->right == 0 ){
+ max_right = DBL_MIN ;
+ }
+ else{
+ max_right = (x->right)->subtree_max;
+ }
+
+ //find max of: x->high(), max_left, max_right
+ Coord temp_max;
+ temp_max = std::max( x->high(), max_left );
+ temp_max = std::max( temp_max, max_right );
+ x->subtree_max = temp_max;
+
+ }
+}
+
+
+RedBlack* RedBlackTree::tree_minimum(RedBlack* x){
+ _REDBLACK_PRINT( "==============================================================" << std::endl << "ENTER: tree_minimum" )
+ while( x->left <- 0 ) {
+ x->left = x;
+ }
+ _REDBLACK_PRINT( "RETURN: tree_minimum" << std::endl )
+ return x;
+}
+
+RedBlack* RedBlackTree::tree_successor(RedBlack* x){
+ _REDBLACK_PRINT( "==============================================================" << std::endl << "ENTER: tree_successor" )
+ if( x->right <- 0 ){
+ _REDBLACK_PRINT( "RETURN: tree_successor" << std::endl )
+ return tree_minimum(x);
+ }
+ RedBlack* y = x->parent;
+ _REDBLACK_PRINT( "y->parent: y->parent" )
+ while( y <- 0 && x == y->right ){
+ x = y;
+ y = y->parent;
+ }
+ _REDBLACK_PRINT( "RETURN: tree_successor" << std::endl )
+ return y;
+}
+
+
+};
+
+/*
+ Local Variables:
+ mode:c++
+ c-file-style:"stroustrup"
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+ indent-tabs-mode:nil
+ fill-column:99
+ End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
diff --git a/src/2geom/orphan-code/rtree.cpp b/src/2geom/orphan-code/rtree.cpp
new file mode 100644
index 0000000..4264292
--- /dev/null
+++ b/src/2geom/orphan-code/rtree.cpp
@@ -0,0 +1,1350 @@
+#include <2geom/orphan-code/rtree.h>
+#include <limits>
+
+/*
+Based on source (BibTex):
+@inproceedings{DBLP:conf/sigmod/Guttman84,
+ author = {Antonin Guttman},
+ title = {R-Trees: A Dynamic Index Structure for Spatial Searching},
+ booktitle = {SIGMOD Conference},
+ year = {1984},
+ pages = {47-57},
+ ee = {http://doi.acm.org/10.1145/602259.602266, db/conf/sigmod/Guttman84.html},
+}
+*/
+
+/*
+#define _RTREE_PRINT(x) std::cout << x << std::endl;
+#define _RTREE_PRINT_TREE( x, y ) print_tree( x, y );
+#define _RTREE_PRINT_TREE_INS( x, y, z ) print_tree( x, y, z );
+*/
+//comment the following if you want output during RTree operations
+
+
+#define _RTREE_PRINT(x) ;
+#define _RTREE_PRINT_TREE( x, y ) ;
+#define _RTREE_PRINT_TREE_INS( x, y, z ) ;
+
+
+
+/*
+TODO 1
+some if(non_leaf)
+ else // leaf
+could be eliminated when function starts from a leaf
+do leaf action
+then repeat function for non-leafs only
+candidates:
+- adjust_tree
+- condense_tree
+
+TODO 2
+generalize in a different way the splitting techniques
+
+*/
+
+
+namespace Geom{
+
+/*=============================================================================
+ insert
+===============================================================================
+insert a new index entry E into the R-tree:
+
+I1) find position of new record:
+ choose_node will find a leaf node L (position) in which to place r
+I2) add record to leaf node:
+ if L has room for another entry install E
+ else split_node will obtain L and LL containing E and all the old entries of L
+ from the available splitting strategies we chose quadtratic-cost algorithm (just to begin
+ with sth)
+ // TODO implement more of them
+I3) propagate changes upward:
+ Invoke adjust_tree on L, also passing LL if a split was performed.
+I4) grow tree taller:
+ if a node spilt propagation, cuased the root to split
+ create new root whose children are the 2 resulting nodes
+*/
+
+void RTree::insert( Rect const &r, unsigned shape ){
+ _RTREE_PRINT("\n=====================================");
+ _RTREE_PRINT("insert");
+ RTreeRecord_Leaf* leaf_record= new RTreeRecord_Leaf( r, shape );
+ insert( *leaf_record );
+}
+
+
+
+void RTree::insert( const RTreeRecord_Leaf &leaf_record,
+ const bool &insert_high /* false */,
+ const unsigned &stop_height /* 0 */,
+ const RTreeRecord_NonLeaf &nonleaf_record /* 0 */
+ )
+{
+ _RTREE_PRINT("\n--------------");
+ _RTREE_PRINT("insert private. element:" << leaf_record.data << " insert high:" << insert_high << " stop height:" << stop_height );
+ RTreeNode *position = 0;
+
+ // if tree is unused create the root Node, not described in source, stupid me :P
+ if(root == 0){
+ root = new RTreeNode();
+ }
+
+ _RTREE_PRINT("I1"); // I1
+ if( insert_high == false ){ // choose leaf node
+ position = choose_node( leaf_record.bounding_box );
+ }
+ else { // choose nonleaf node
+ position = choose_node( nonleaf_record.bounding_box, insert_high, stop_height );
+ }
+ _RTREE_PRINT("leaf node chosen: " );
+ _RTREE_PRINT_TREE( position , 0 );
+ std::pair< RTreeNode*, RTreeNode* > node_division;
+
+ bool split_performed = false;
+
+ if( position->children_nodes.size() > 0 ){ // non-leaf node: position
+ // we must reach here only to insert high non leaf node, not insert leaf node
+ assert( insert_high == true );
+
+ // put new element in node temporarily. Later on, if we need to, we will split the node.
+ position->children_nodes.push_back( nonleaf_record );
+ if( position->children_nodes.size() <= max_records ){
+ _RTREE_PRINT("I2 nonleaf: no split: " << position->children_nodes.size() ); // I2
+ }
+ else{
+ _RTREE_PRINT("I2 nonleaf: split: " << position->children_nodes.size() ); // I2
+ node_division = split_node( position );
+ split_performed = true;
+ }
+
+ }
+ else { // leaf node: position:
+ // we must reach here only to insert leaf node, not insert high non leaf node
+ assert( insert_high == false );
+
+
+ // put new element in node temporarily. Later on, if we need to, we will split the node.
+ position->children_leaves.push_back( leaf_record );
+ if( position->children_leaves.size() <= max_records ){
+ _RTREE_PRINT("I2 leaf: no split: " << position->children_leaves.size() ); // I2
+ }
+ else{
+ _RTREE_PRINT("I2 leaf: split: " << position->children_leaves.size() << " max_records:" << max_records); // I2
+ node_division = split_node( position );
+ split_performed = true;
+
+ _RTREE_PRINT(" group A");
+ _RTREE_PRINT_TREE( node_division.first , 3 );
+ _RTREE_PRINT(" group B");
+ _RTREE_PRINT_TREE( node_division.second , 3 );
+
+ }
+
+ }
+
+ _RTREE_PRINT("I3"); // I3
+ bool root_split_performed = adjust_tree( position, node_division, split_performed );
+ _RTREE_PRINT("root split: " << root_split_performed);
+
+
+// _RTREE_PRINT("TREE:");
+// print_tree( root , 2 );
+
+ _RTREE_PRINT("I4"); // I4
+ if( root_split_performed ){
+ std::pair<RTreeNode*, RTreeNode*> root_division;
+ root_division = quadratic_split( root ); // AT5
+
+ Rect first_record_bounding_box;
+ Rect second_record_bounding_box;
+
+ RTreeRecord_NonLeaf first_new_record = create_nonleaf_record_from_rtreenode( first_record_bounding_box, root_division.first );
+ RTreeRecord_NonLeaf second_new_record = create_nonleaf_record_from_rtreenode( second_record_bounding_box, root_division.second );
+ _RTREE_PRINT(" 1st:");
+ _RTREE_PRINT_TREE( first_new_record.data, 5 );
+ _RTREE_PRINT(" 2nd:");
+ _RTREE_PRINT_TREE( second_new_record.data, 5 );
+
+ // *new* root is by definition non-leaf. Install the new records there
+ RTreeNode* new_root = new RTreeNode();
+ new_root->children_nodes.push_back( first_new_record );
+ new_root->children_nodes.push_back( second_new_record );
+
+ delete root;
+
+ root = new_root;
+ tree_height++; // increse tree height
+
+ _RTREE_PRINT_TREE( root, 5 );
+ sanity_check( root, 0 );
+ }
+ _RTREE_PRINT("done");
+
+ /*
+ the node_division.second is saved on the tree
+ the node_division.first was copied in existing tree in node
+ so we don't need this anymore
+ */
+ delete node_division.first;
+}
+
+/* I1 =========================================================================
+
+original: choose_node will find a leaf node L in which to place r
+changed to choose_node will find a node L in which to place r
+the node L is:
+non-leaf: if flag is set. the height of the node is insert_at_height
+leaf: if flag is NOT set
+
+1) Initialize: set N to be the root node
+2) Leaf Check:
+ insert_height = false
+ if N is leaf return N
+ insert_height = true
+3) Choose subtree: If N not leaf OR not we are not in the proper height then
+ let F be an entry in N whose rect Fi needs least enlargement to include r
+ ties resolved with rect of smallest area
+4) descend until a leaf is reached OR proper height is reached: set N to the child node pointed to by F and goto 2.
+*/
+
+// TODO keep stack with visited nodes
+
+RTreeNode* RTree::choose_node( const Rect &r, const bool &insert_high /* false */, const unsigned &stop_height /* 0 */) const {
+
+ _RTREE_PRINT(" CL1");// CL1
+ RTreeNode *pos = root;
+
+ double min_enlargement;
+ double current_enlargement;
+ int node_min_enlargement;
+ unsigned current_height = 0; // zero is the root
+
+ _RTREE_PRINT(" CL2 current_height:" << current_height << " stop_height:" << stop_height << " insert_high:" << insert_high);
+ // CL2 Leaf check && Height check
+ while( ( insert_high ? true : pos->children_nodes.size() != 0 )
+ && ( insert_high ? current_height < stop_height : true ) )
+ /* Leaf check, during insert leaf */
+ /* node height check, during insert non-leaf */
+ {
+ _RTREE_PRINT(" CL3 current_height:" << current_height << " stop_height:" << stop_height ); // CL3
+ min_enlargement = std::numeric_limits<double>::max();
+ current_enlargement = 0;
+ node_min_enlargement = 0;
+
+ for(unsigned i=0; i< pos->children_nodes.size(); i++){
+ current_enlargement = find_enlargement( pos->children_nodes[i].bounding_box, r );
+
+ // TODO tie not solved!
+ if( current_enlargement < min_enlargement ){
+ node_min_enlargement = i;
+ min_enlargement = current_enlargement;
+ }
+ }
+ _RTREE_PRINT(" CL4"); // CL4
+ // descend to the node with the min_enlargement
+ pos = pos->children_nodes[node_min_enlargement].data;
+ current_height++; // increase current visiting height
+ }
+
+ return pos;
+}
+
+
+/*
+find_enlargement:
+
+enlargement that "a" needs in order to include "b"
+b is the new rect we want to insert.
+a is the rect of the node we try to see if b should go in.
+*/
+double RTree::find_enlargement( const Rect &a, const Rect &b ) const{
+
+
+ Rect union_rect(a);
+ union_rect.unionWith(b);
+
+ OptRect a_intersection_b = intersect( a, b );
+
+ // a, b do not intersect
+ if( a_intersection_b.empty() ){
+ _RTREE_PRINT(" find_enlargement (no intersect): " << union_rect.area() - a.area() - b.area() );
+ return union_rect.area() - a.area() - b.area();
+ }
+
+ // a, b intersect
+
+ // a contains b
+ if( a.contains( b ) ){
+ _RTREE_PRINT(" find_enlargement (intersect: a cont b): " << a.area() - b.area() );
+ //return a.area() - b.area();
+ return b.area() - a.area(); // enlargement is negative in this case.
+ }
+
+ // b contains a
+ if( b.contains( a ) ){
+ _RTREE_PRINT(" find_enlargement (intersect: b cont a): " << a.area() - b.area() );
+ return b.area() - a.area();
+ }
+
+ // a partially cover b
+ _RTREE_PRINT(" find_enlargement (intersect: a partial cover b): " << union_rect.area() - a.area() - b.area() - a_intersection_b->area() );
+ return union_rect.area()
+ - ( a.area() - a_intersection_b->area() )
+ - ( b.area() - a_intersection_b->area() );
+}
+
+
+/* I2 =========================================================================
+use one split strategy
+*/
+
+std::pair<RTreeNode*, RTreeNode*> RTree::split_node( RTreeNode *s ){
+/*
+ if( split_strategy == LINEAR_COST ){
+ linear_cost_split( ............. );
+ }
+*/
+ return quadratic_split( s ); // else QUADRATIC_SPIT
+}
+
+
+/*-----------------------------------------------------------------------------
+ Quadratic Split
+
+QS1) Pick first entry for each group:
+ Appy pick_seeds to choose 2 entries to be the first elements of the groups. Assign each one of
+ them to one group
+QS2) check if done:
+ a) if all entries have been assinged stop
+ b) if one group has so few entries that all the rest must be assignmed to it, in order for it to
+ have the min number , assign them and stop
+QS3) select entry and assign:
+ Inkvoke pick_next() to choose the next entry to assign.
+ *[in pick_next] Add it to the group whose covering rectangle will have to be enlrarged least to
+ accommodate it. Resolve ties by adding the entry to the group with the smaller are, then to the
+ one with fewer entries, then to either of the two.
+ goto 2.
+*/
+std::pair<RTreeNode*, RTreeNode*> RTree::quadratic_split( RTreeNode *s ) {
+
+ // s is the original leaf node or non-leaf node
+ RTreeNode* group_a = new RTreeNode(); // a is the 1st group
+ RTreeNode* group_b = new RTreeNode(); // b is the 2nd group
+
+
+ _RTREE_PRINT(" QS1"); // QS1
+ std::pair<unsigned, unsigned> initial_seeds;
+ initial_seeds = pick_seeds(s);
+
+ // if non-leaf node: s
+ if( s->children_nodes.size() > 0 ){
+ _RTREE_PRINT(" non-leaf node");
+ // each element is true if the node has been assinged to either "a" or "b"
+ std::vector<bool> assigned_v( s->children_nodes.size() );
+ std::fill( assigned_v.begin(), assigned_v.end(), false );
+
+ group_a->children_nodes.push_back( s->children_nodes[initial_seeds.first] );
+ assert(initial_seeds.first < assigned_v.size());
+ assigned_v[ initial_seeds.first ] = true;
+
+ group_b->children_nodes.push_back( s->children_nodes[initial_seeds.second] );
+ assert(initial_seeds.second < assigned_v.size());
+ assigned_v[ initial_seeds.second ] = true;
+
+ _RTREE_PRINT(" QS2"); // QS2
+ unsigned num_of_not_assigned = s->children_nodes.size() - 2;
+ // so far we have assinged 2 out of all
+
+ while( num_of_not_assigned ){// QS2 a
+ _RTREE_PRINT(" QS2 b, num_of_not_assigned:" << num_of_not_assigned); // QS2 b
+ /*
+ we are on NON leaf node so children of split groups must be nodes
+
+ Check each group to see if one group has so few entries that all the rest must
+ be assignmed to it, in order for it to have the min number.
+ */
+ if( group_a->children_nodes.size() + num_of_not_assigned <= min_records ){
+ // add the non-assigned to group_a
+ for(unsigned i = 0; i < assigned_v.size(); i++){
+ if(assigned_v[i] == false){
+ group_a->children_nodes.push_back( s->children_nodes[i] );
+ assigned_v[i] = true;
+ }
+ }
+ break;
+ }
+
+ if( group_b->children_nodes.size() + num_of_not_assigned <= min_records ){
+ // add the non-assigned to group_b
+ for( unsigned i = 0; i < assigned_v.size(); i++ ){
+ if( assigned_v[i] == false ){
+ group_b->children_nodes.push_back( s->children_nodes[i] );
+ assigned_v[i] = true;
+ }
+ }
+ break;
+ }
+
+ _RTREE_PRINT(" QS3"); // QS3
+ std::pair<unsigned, enum_add_to_group> next_element;
+ next_element = pick_next( group_a, group_b, s, assigned_v );
+ if( next_element.second == ADD_TO_GROUP_A ){
+ group_a->children_nodes.push_back( s->children_nodes[ next_element.first ] );
+ }
+ else{
+ group_b->children_nodes.push_back( s->children_nodes[ next_element.first ] );
+ }
+
+ num_of_not_assigned--;
+ }
+ }
+ // else leaf node: s
+ else{
+ _RTREE_PRINT(" leaf node");
+ // each element is true if the node has been assinged to either "a" or "b"
+ std::vector<bool> assigned_v( s->children_leaves.size() );
+ std::fill( assigned_v.begin(), assigned_v.end(), false );
+
+ // assign 1st seed to group a
+ group_a->children_leaves.push_back( s->children_leaves[initial_seeds.first] );
+ assert(initial_seeds.first < assigned_v.size());
+ assigned_v[ initial_seeds.first ] = true;
+
+ // assign 2nd seed to group b
+ group_b->children_leaves.push_back( s->children_leaves[initial_seeds.second] );
+ assert(initial_seeds.second < assigned_v.size());
+ assigned_v[ initial_seeds.second ] = true;
+
+ _RTREE_PRINT(" QS2"); // QS2
+ unsigned num_of_not_assigned = s->children_leaves.size() - 2;
+ // so far we have assinged 2 out of all
+
+ while( num_of_not_assigned ){// QS2 a
+ _RTREE_PRINT(" QS2 b, num_of_not_assigned:" << num_of_not_assigned); // QS2 b
+ /*
+ we are on leaf node so children of split groups must be leaves
+
+ Check each group to see if one group has so few entries that all the rest must
+ be assignmed to it, in order for it to have the min number.
+ */
+ if( group_a->children_leaves.size() + num_of_not_assigned <= min_records ){
+ _RTREE_PRINT(" add the non-assigned to group_a");
+ // add the non-assigned to group_a
+ for( unsigned i = 0; i < assigned_v.size(); i++ ){
+ if( assigned_v[i] == false ){
+ group_a->children_leaves.push_back( s->children_leaves[i] );
+ assigned_v[i] = true;
+ }
+ }
+ break;
+ }
+
+ if( group_b->children_leaves.size() + num_of_not_assigned <= min_records ){
+ _RTREE_PRINT(" add the non-assigned to group_b");
+ // add the non-assigned to group_b
+ for( unsigned i = 0; i < assigned_v.size(); i++ ){
+ if( assigned_v[i] == false ){
+ group_b->children_leaves.push_back( s->children_leaves[i] );
+ assigned_v[i] = true;
+ }
+ }
+ break;
+ }
+
+ _RTREE_PRINT(" QS3"); // QS3
+ std::pair<unsigned, enum_add_to_group> next_element;
+ next_element = pick_next(group_a, group_b, s, assigned_v);
+ if( next_element.second == ADD_TO_GROUP_A ){
+ group_a->children_leaves.push_back( s->children_leaves[ next_element.first ] );
+ }
+ else{
+ group_b->children_leaves.push_back( s->children_leaves[ next_element.first ] );
+ }
+
+ num_of_not_assigned--;
+ }
+ }
+ assert( initial_seeds.first != initial_seeds.second );
+ return std::make_pair( group_a, group_b );
+}
+
+/*
+PS1) caclulate ineffeciency of grouping entries together:
+ Foreach pair of entries E1 (i), E2 (j) compose rectangle J (i_union_j) including E1, E2.
+ Calculate d = area(i_union_j) - area(i) - area(j)
+PS2) choose the most wastefull pair:
+ Choose pair with largest d
+*/
+
+std::pair<unsigned, unsigned> RTree::pick_seeds( RTreeNode *s ) const{
+ double current_d = 0;
+ double max_d = std::numeric_limits<double>::min();
+ unsigned seed_a = 0;
+ unsigned seed_b = 1;
+ _RTREE_PRINT(" pick_seeds");
+
+ // if non leaf node: s
+ if( s->children_nodes.size() > 0 ){
+ _RTREE_PRINT(" non leaf");
+ _RTREE_PRINT(" PS1"); // PS1
+ for( unsigned a = 0; a < s->children_nodes.size(); a++ ){
+ // with j=i we check only the upper (diagonal) half
+ // with j=i+1 we also avoid checking for b==a (we don't need it)
+ for( unsigned b = a+1; b < s->children_nodes.size(); b++ ){
+ _RTREE_PRINT(" PS2 " << a << " - " << b ); // PS2
+ current_d = find_waste_area( s->children_nodes[a].bounding_box, s->children_nodes[b].bounding_box );
+
+ if( current_d > max_d ){
+ max_d = current_d;
+ seed_a = a;
+ seed_b = b;
+ }
+ }
+ }
+ }
+ // else leaf node: s
+ else{
+ _RTREE_PRINT(" leaf node");
+ _RTREE_PRINT(" PS1"); // PS1
+ for( unsigned a = 0; a < s->children_leaves.size(); a++ ){
+ // with j=i we check only the upper (diagonal) half
+ // with j=i+1 we also avoid checking for j==i (we don't need this one)
+ for( unsigned b = a+1; b < s->children_leaves.size(); b++ ){
+ _RTREE_PRINT(" PS2 " << s->children_leaves[a].data << ":" << s->children_leaves[a].bounding_box.area()
+ << " - " << s->children_leaves[b].data << ":" << s->children_leaves[b].bounding_box.area() ); // PS2
+ current_d = find_waste_area( s->children_leaves[a].bounding_box, s->children_leaves[b].bounding_box );
+
+ if( current_d > max_d ){
+ max_d = current_d;
+ seed_a = a;
+ seed_b = b;
+ }
+ }
+ }
+ }
+ _RTREE_PRINT(" seed_a: " << seed_a << " seed_b: " << seed_b );
+ return std::make_pair( seed_a, seed_b );
+}
+
+/*
+find_waste_area (used in pick_seeds step 1)
+
+for a pair A, B compose a rect union_rect that includes a and b
+calculate area of union_rect - area of a - area b
+*/
+double RTree::find_waste_area( const Rect &a, const Rect &b ) const{
+ Rect union_rect(a);
+ union_rect.unionWith(b);
+
+ return union_rect.area() - a.area() - b.area();
+}
+
+/*
+pick_next:
+select one remaining entry for classification in a group
+
+PN1) Determine cost of putting each entry in each group:
+ Foreach entry E not yet in a group, calculate
+ d1= area increase required in the cover rect of Group 1 to include E
+ d2= area increase required in the cover rect of Group 2 to include E
+PN2) Find entry with greatest preference for each group:
+ Choose any entry with the maximum difference between d1 and d2
+
+*/
+
+std::pair<unsigned, enum_add_to_group> RTree::pick_next( RTreeNode* group_a,
+ RTreeNode* group_b,
+ RTreeNode* s,
+ std::vector<bool> &assigned_v )
+{
+ double max_increase_difference = std::numeric_limits<double>::min();
+ unsigned max_increase_difference_node = 0;
+ double current_increase_difference = 0;
+
+ enum_add_to_group group_to_add = ADD_TO_GROUP_A;
+
+ /*
+ bounding boxes of the 2 new groups. This info isn't available, since they
+ have no parent nodes (so that the parent node knows the bounding box).
+ */
+ Rect bounding_box_a;
+ Rect bounding_box_b;
+
+ double increase_area_a = 0;
+ double increase_area_b = 0;
+
+ _RTREE_PRINT(" pick_next, assigned_v.size:" << assigned_v.size() );
+
+ // if non leaf node: one of the 2 groups (both groups are the same, either leaf/nonleaf)
+ if( group_a->children_nodes.size() > 0 ){
+ _RTREE_PRINT(" non leaf");
+
+ // calculate the bounding boxes of the 2 new groups.
+ bounding_box_a = Rect( group_a->children_nodes[0].bounding_box );
+ for( unsigned i = 1; i < group_a->children_nodes.size(); i++ ){
+ bounding_box_a.unionWith( group_a->children_nodes[i].bounding_box );
+ }
+
+ bounding_box_b = Rect( group_b->children_nodes[0].bounding_box );
+ for( unsigned i = 1; i < group_b->children_nodes.size(); i++ ){
+ bounding_box_b.unionWith( group_b->children_nodes[i].bounding_box );
+ }
+
+
+ _RTREE_PRINT(" PN1"); // PN1
+ for( unsigned i = 0; i < assigned_v.size(); i++ ){
+ _RTREE_PRINT(" i:" << i << " assigned:" << assigned_v[i]);
+ if( assigned_v[i] == false ){
+
+ increase_area_a = find_enlargement( bounding_box_a, s->children_nodes[i].bounding_box );
+ increase_area_b = find_enlargement( bounding_box_b, s->children_nodes[i].bounding_box );
+
+ current_increase_difference = std::abs( increase_area_a - increase_area_b );
+ _RTREE_PRINT(" PN2 " << i << ": " << current_increase_difference ); // PN2
+ if( current_increase_difference > max_increase_difference ){
+ max_increase_difference = current_increase_difference;
+ max_increase_difference_node = i;
+
+ // TODO tie not solved!
+ if( increase_area_a < increase_area_b ){
+ group_to_add = ADD_TO_GROUP_A;
+ }
+ else{
+ group_to_add = ADD_TO_GROUP_B;
+ }
+ }
+ }
+ }
+ //assert(max_increase_difference_node >= 0);
+ assert(max_increase_difference_node < assigned_v.size());
+ assigned_v[max_increase_difference_node] = true;
+ _RTREE_PRINT(" ... i:" << max_increase_difference_node << " assigned:" << assigned_v[max_increase_difference_node] );
+ }
+ else{ // else leaf node
+ _RTREE_PRINT(" leaf");
+
+ // calculate the bounding boxes of the 2 new groups
+ bounding_box_a = Rect( group_a->children_leaves[0].bounding_box );
+ for( unsigned i = 1; i < group_a->children_leaves.size(); i++ ){
+ std::cout<< " lala";
+ bounding_box_a.unionWith( group_a->children_leaves[i].bounding_box );
+ }
+
+ bounding_box_b = Rect( group_b->children_leaves[0].bounding_box );
+ for( unsigned i = 1; i < group_b->children_leaves.size(); i++ ){
+ std::cout<< " lala";
+ bounding_box_b.unionWith( group_b->children_leaves[i].bounding_box );
+ }
+ std::cout<< "" << std::endl;
+
+ _RTREE_PRINT(" PN1"); // PN1
+ for( unsigned i = 0; i < assigned_v.size(); i++ ){
+ _RTREE_PRINT(" i:" << i << " assigned:" << assigned_v[i]);
+ if( assigned_v[i] == false ){
+ increase_area_a = find_enlargement( bounding_box_a, s->children_leaves[i].bounding_box );
+ increase_area_b = find_enlargement( bounding_box_b, s->children_leaves[i].bounding_box );
+
+ current_increase_difference = std::abs( increase_area_a - increase_area_b );
+ _RTREE_PRINT(" PN2 " << i << ": " << current_increase_difference ); // PN2
+
+ if( current_increase_difference > max_increase_difference ){
+ max_increase_difference = current_increase_difference;
+ max_increase_difference_node = i;
+
+ // TODO tie not solved!
+ if( increase_area_a < increase_area_b ){
+ group_to_add = ADD_TO_GROUP_A;
+ }
+ else{
+ group_to_add = ADD_TO_GROUP_B;
+ }
+ }
+ }
+ }
+ assert(max_increase_difference_node < assigned_v.size());
+ assigned_v[max_increase_difference_node] = true;
+ _RTREE_PRINT(" ... i:" << max_increase_difference_node << " assigned:" << assigned_v[max_increase_difference_node] );
+ }
+
+ _RTREE_PRINT(" node:" << max_increase_difference_node << " added:" << group_to_add );
+ return std::make_pair( max_increase_difference_node, group_to_add );
+}
+
+/* I3 =========================================================================
+
+adjust_tree:
+Ascend from a leaf node L to root, adjusting covering rectangles and propagating node splits as
+necessary
+
+We modified this one from the source in the step AT1 and AT5
+
+AT1) Initialize:
+ Set N=L
+ IF L was spilt previously, set NN to be the resulting second node AND
+ (not mentioned in the original source but that's what it should mean)
+ Assign all entries of first node to L
+AT2) check if done:
+ IF N is root stop
+AT3) adjust covering rectangle in parent entry
+ 1) Let P be the parent of N
+ 2) Let EN be the N's entry in P
+ 3) Adjust EN bounding box so that it tightly enclosses all entry rectangles in N
+AT4) Propagate node split upward
+ IF N has a partner NN resulting from an earlier split
+ create a new entry ENN with ENN "p" pointing to NN and ENN bounding box enclosing all
+ rectangles in NN
+
+ IF there is room in P add ENN
+ ELSE invoke split_node to produce P an PP containing ENN and all P's old entries.
+AT5) Move up to next level
+ Set N=P,
+ IF a split occurred, set NN=PP
+ goto AT1 (was goto AT2)
+*/
+
+bool RTree::adjust_tree( RTreeNode* position,
+ // modified: it holds the last split group
+ std::pair<RTreeNode*, RTreeNode*> &node_division,
+ bool initial_split_performed)
+{
+ RTreeNode* parent;
+ unsigned child_in_parent; // the element in parent node that points to current posistion
+ std::pair< RTreeNode*, bool > find_result;
+ bool split_performed = initial_split_performed;
+ bool root_split_performed = false;
+
+ _RTREE_PRINT(" adjust_tree");
+ _RTREE_PRINT(" AT1");
+
+ while( true ){
+ _RTREE_PRINT(" ------- current tree status:");
+ _RTREE_PRINT_TREE_INS(root, 2, true);
+
+ // check for loop BREAK
+ if( position == root ){
+ _RTREE_PRINT(" AT2: found root");
+ if( split_performed ){
+ root_split_performed = true;
+ }
+ break;
+ }
+
+ if( split_performed ){
+ copy_group_a_to_existing_node( position, node_division.first );
+ }
+
+ /*
+ pick randomly, let's say the 1st entry of the current node.
+ Search for this spatial area in the tree, and stop to the parent node.
+ Then find position of current node pointer, in the parent node.
+ */
+ _RTREE_PRINT(" AT3.1"); // AT3.1 Let P be the parent of N
+ if( position->children_nodes.size() > 0 ){
+ find_result = find_parent( root, position->children_nodes[0].bounding_box, position);
+ }
+ else{
+ find_result = find_parent( root, position->children_leaves[0].bounding_box, position);
+ }
+ parent = find_result.first;
+
+ // parent is a non-leaf, by definition
+ _RTREE_PRINT(" AT3.2"); // AT3.2 Let EN be the N's entry in P
+ for( child_in_parent = 0; child_in_parent < parent->children_nodes.size(); child_in_parent++ ){
+ if( parent->children_nodes[ child_in_parent ].data == position){
+ _RTREE_PRINT(" child_in_parent: " << child_in_parent);
+ break;
+ }
+ }
+
+ _RTREE_PRINT(" AT3.3");
+ // AT3.2 Adjust EN bounding box so that it tightly enclosses all entry rectangles in N
+ recalculate_bounding_box( parent, position, child_in_parent );
+
+
+ _RTREE_PRINT(" AT4"); // AT4
+ if( split_performed ){
+ // create new record (from group_b)
+ //RTreeNode* new_node = new RTreeNode();
+ Rect new_record_bounding;
+
+ RTreeRecord_NonLeaf new_record = create_nonleaf_record_from_rtreenode( new_record_bounding, node_division.second );
+
+ // install new entry (group_b)
+ if( parent->children_nodes.size() < max_records ){
+ parent->children_nodes.push_back( new_record );
+ split_performed = false;
+ }
+ else{
+ parent->children_nodes.push_back( new_record );
+ node_division = quadratic_split( parent ); // AT5
+ split_performed = true;
+ }
+
+ }
+ _RTREE_PRINT(" AT5"); // AT5
+ position = parent;
+ }
+
+ return root_split_performed;
+}
+
+/*
+find_parent:
+The source only mentions that we should "find" the parent. But it doesn't seay how. So we made a
+modification of search.
+
+Initially we take the root, a rect of the node, of which the parent we look for and the node we seek
+
+We do a spatial search for this rect. If we find get an intersecttion with the rect we check if the
+child is the one we look for.
+If not we call find_parent again recursively
+*/
+
+std::pair< RTreeNode*, bool > RTree::find_parent( RTreeNode* subtree_root,
+ Rect search_area,
+ RTreeNode* wanted) const
+{
+ _RTREE_PRINT("find_parent");
+
+ std::pair< RTreeNode*, bool > result;
+ if( subtree_root->children_nodes.size() > 0 ){
+
+ for( unsigned i=0; i < subtree_root->children_nodes.size(); i++ ){
+ if( subtree_root->children_nodes[i].data == wanted){
+ _RTREE_PRINT("FOUND!!"); // non leaf node
+ return std::make_pair( subtree_root, true );
+ }
+
+ if( subtree_root->children_nodes[i].bounding_box.intersects( search_area ) ){
+ result = find_parent( subtree_root->children_nodes[i].data, search_area, wanted);
+ if ( result.second ){
+ break;
+ }
+ }
+ }
+ }
+ return result;
+}
+
+
+void RTree::copy_group_a_to_existing_node( RTreeNode *position, RTreeNode* group_a ){
+ // clear position (the one that was split) and put there all the nodes of group_a
+ if( position->children_nodes.size() > 0 ){
+ _RTREE_PRINT(" copy_group...(): install group A to existing non-leaf node");
+ // non leaf-node: position
+ position->children_nodes.clear();
+ for(auto & children_node : group_a->children_nodes){
+ position->children_nodes.push_back( children_node );
+ }
+ }
+ else{
+ _RTREE_PRINT(" copy_group...(): install group A to existing leaf node");
+ // leaf-node: positions
+ position->children_leaves.clear();
+ for(auto & children_leave : group_a->children_leaves){
+ position->children_leaves.push_back( children_leave );
+ }
+ }
+}
+
+
+
+RTreeRecord_NonLeaf RTree::create_nonleaf_record_from_rtreenode( Rect &new_entry_bounding, RTreeNode* rtreenode ){
+
+ if( rtreenode->children_nodes.size() > 0 ){
+ // found bounding box of new entry
+ new_entry_bounding = Rect( rtreenode->children_nodes[0].bounding_box );
+ for(unsigned i = 1; i < rtreenode->children_nodes.size(); i++ ){
+ new_entry_bounding.unionWith( rtreenode->children_nodes[ i ].bounding_box );
+ }
+ }
+ else{ // non leaf: rtreenode
+ // found bounding box of new entry
+ new_entry_bounding = Rect( rtreenode->children_leaves[0].bounding_box );
+ for(unsigned i = 1; i < rtreenode->children_leaves.size(); i++ ){
+ new_entry_bounding.unionWith( rtreenode->children_leaves[ i ].bounding_box );
+ }
+ }
+ return RTreeRecord_NonLeaf( new_entry_bounding, rtreenode );
+}
+
+
+
+/*
+ print the elements of the tree
+ based on ordered tree walking
+*/
+void RTree::print_tree(RTreeNode* subtree_root, int depth ) const{
+
+ if( subtree_root->children_nodes.size() > 0 ){
+
+ // descend in each one of the elements and call print_tree
+ for( unsigned i=0; i < subtree_root->children_nodes.size(); i++ ){
+ //print spaces for indentation
+ for(int j=0; j < depth; j++){
+ std::cout << " " ;
+ }
+
+ std::cout << subtree_root->children_nodes[i].bounding_box << ", " << subtree_root->children_nodes.size() << std::endl ;
+ _RTREE_PRINT_TREE_INS( subtree_root->children_nodes[i].data, depth+1, used_during_insert);
+ }
+
+ }
+ else{
+ for(int j=0; j < depth; j++){
+ std::cout << " " ;
+ }
+ std::cout << subtree_root->children_leaves.size() << ": " ;
+
+ // print all the elements of the leaf node
+ for(auto & children_leave : subtree_root->children_leaves){
+ std::cout << children_leave.data << ", " ;
+ }
+ std::cout << std::endl ;
+
+ }
+}
+
+
+void RTree::sanity_check(RTreeNode* subtree_root, int depth, bool used_during_insert ) const{
+
+ if( subtree_root->children_nodes.size() > 0 ){
+ // descend in each one of the elements and call sanity_check
+ for(auto & children_node : subtree_root->children_nodes){
+ sanity_check( children_node.data, depth+1, used_during_insert);
+ }
+
+
+ // sanity check
+ if( subtree_root != root ){
+ assert( subtree_root->children_nodes.size() >= min_records);
+ }
+/*
+ else{
+ assert( subtree_root->children_nodes.size() >= 1);
+ }
+*/
+
+ if( used_during_insert ){
+ // allow one more during for insert
+ assert( subtree_root->children_nodes.size() <= max_records + 1 );
+ }
+ else{
+ assert( subtree_root->children_nodes.size() <= max_records );
+ }
+
+ }
+ else{
+ // sanity check
+ if( subtree_root != root ){
+ assert( subtree_root->children_leaves.size() >= min_records);
+ }
+/*
+ else{
+ assert( subtree_root->children_leaves.size() >= 1);
+ }
+*/
+
+ if( used_during_insert ){
+ // allow one more during for insert
+ assert( subtree_root->children_leaves.size() <= max_records + 1 );
+ }
+ else{
+ assert( subtree_root->children_nodes.size() <= max_records );
+ }
+ }
+}
+
+
+
+/*=============================================================================
+ search
+===============================================================================
+Given an RTree whose root node is T find all index records whose rects overlap search rect S
+S1) Search subtrees:
+ IF T isn't a leaf, check every entry E to determine whether E I overlaps S
+ FOR all overlapping entries invoke Search on the tree whose root node is pointed by E P
+S2) ELSE T is leaf
+ check all entries E to determine whether E I overlaps S
+ IF so E is a qualifying record
+*/
+
+
+void RTree::search( const Rect &search_area, std::vector< int >* result, const RTreeNode* subtree ) const {
+ // S1
+ if( subtree->children_nodes.size() > 0 ){ // non-leaf: subtree
+ for(const auto & children_node : subtree->children_nodes){
+ if( children_node.bounding_box.intersects( search_area ) ){
+ search( search_area, result, children_node.data );
+ }
+ }
+ }
+ // S2
+ else{ // leaf: subtree
+ for(const auto & children_leave : subtree->children_leaves){
+ if( children_leave.bounding_box.intersects( search_area ) ){
+ result->push_back( children_leave.data );
+ }
+ }
+ }
+}
+
+
+/*=============================================================================
+ erase
+===============================================================================
+we changed steps D2)
+D1) Find node containing record
+ Invoke find_leaf() to locate the leaf node L containing E
+ IF record is found stop
+D2) Delete record
+ Remove E from L (it happened in find_leaf step FL2)
+D3) Propagate changes
+ Invoke condense_tree, passing L
+D4) Shorten tree
+ If root node has only one child, after the tree was adjusted, make the child new root
+
+return
+0 on success
+1 in case no entry was found
+
+*/
+//int RTree::erase( const RTreeRecord_Leaf & record_to_erase ){
+int RTree::erase( const Rect &search_area, const int shape_to_delete ){
+ _RTREE_PRINT("\n=====================================");
+ _RTREE_PRINT("erase element: " << shape_to_delete);
+ // D1 + D2: entry is deleted in find_leaf
+ _RTREE_PRINT("D1 & D2 : find and delete the leaf");
+ RTreeNode* contains_record = find_leaf( root, search_area, shape_to_delete );
+ if( !contains_record ){ // no entry returned from find_leaf
+ return 1; // no entry found
+ }
+
+ // D3
+ //bool root_elimination_performed = condense_tree( contains_record );
+
+ // D4
+
+ //if( root_elimination_performed ){
+ if( root->children_nodes.size() > 0 ){ // non leaf: root
+ // D4
+ if( root->children_nodes.size() == 1 ){
+ _RTREE_PRINT("D4 : non leaf: ");
+ tree_height--;
+ RTreeNode* t = root;
+ root = root->children_nodes[0].data;
+ delete t;
+ }
+
+ }
+ else { // leaf: root
+ // D4
+ // do nothing
+ }
+ sanity_check( root, 0 );
+ return 0; // success
+}
+
+
+/*
+ find_leaf()
+Given an RTree whose root node is T find the leaf node containing index entry E
+
+FL1) Search subtrees
+ IF T is non leaf, check each entry F in T to determine if F I overlaps E I
+ foreach such entry invoke find_leaf on the tree whose root is pointed to by F P until E is
+ found or all entries have been checked
+FL2) search leaf node for record
+ IF T is leaf, check each entry to see if it matches E
+ IF E is found return T
+ AND delete element E (step D2)
+*/
+
+RTreeNode* RTree::find_leaf( RTreeNode* subtree, const Rect &search_area, const int shape_to_delete ) const {
+ // FL1
+ if( subtree->children_nodes.size() > 0 ){ // non-leaf: subtree
+ for(auto & children_node : subtree->children_nodes){
+ if( children_node.bounding_box.intersects( search_area ) ){
+ RTreeNode* t = find_leaf( children_node.data, search_area, shape_to_delete );
+ if( t ){ // if search was successful terminate
+ return t;
+ }
+ }
+ }
+ }
+ // FL2
+ else{ // leaf: subtree
+ for( std::vector< RTreeRecord_Leaf >::iterator it = subtree->children_leaves.begin(); it!=subtree->children_leaves.end(); ++it ){
+ if( it->data == shape_to_delete ){
+ // delete element: implement step D2)
+ subtree->children_leaves.erase( it );
+ return subtree;
+ }
+ }
+ }
+ return 0;
+}
+
+
+/*
+ condense_tree()
+Given a Leaf node L from which an entry has been deleted eliminate the node if it has too few entries and reallocate its entries
+Propagate node elimination upwards as necessary
+Adjust all covering recsts n the path to the root making them smaller if possible
+
+CT1) Initialize
+ Set N=L
+ Set Q the set of eliminated nodes to be empty
+CT2) // Find parent entry (was there but doesn't make sense)
+ IF N is the root
+ goto CT6
+ ELSE
+ 1) Find parent entry
+ 2) let P be the parent of N
+ 3) and let EN be N's entry in P
+CT3) IF N has fewer than m entries
+ Eliminate underfull node
+ 1) delete EN from P
+ 2) and add N to set Q
+CT4) ELSE
+ adjust EN I to tightly contain all entries in N
+CT5) move up one level in tree
+ set N=P and repeat from CT2
+
+CT6) Re insert orphaned entries
+ Re-inser all entreis of nodes in set Q
+ Entries from eliminated leaf nodes are re-inserted in tree leaves (like in insert)
+ BUT non-leaf nodes must be placed higher in the tree so that leaves of their dependent subtrees
+ will be on the same level as leaves of the main tree. (on the same height they originally were)
+ (not mentioned in the source description: the criteria for placing the node should be
+ again TODO ??? least enlargement)
+
+*/
+// TODO this can be merged with adjust_tree or refactor to reutilize some parts. less readable
+bool RTree::condense_tree( RTreeNode* position )
+{
+ RTreeNode* parent;
+ unsigned child_in_parent = 0; // the element in parent node that points to current posistion
+
+ std::pair< RTreeNode*, bool > find_result;
+ bool elimination_performed = false;
+ bool root_elimination_performed = false;
+ unsigned current_height = tree_height+1;
+ Rect special_case_bounding_box;
+ _RTREE_PRINT(" condense_tree");
+ _RTREE_PRINT(" CT1");
+ // leaf records that were eliminated due to under-full node
+ std::vector< RTreeRecord_Leaf > Q_leaf_records( 0 );
+
+ // < non-leaf records, their height > that were eliminated due to under-full node
+ std::vector< std::pair< RTreeRecord_NonLeaf, unsigned > > Q_nonleaf_records( 0 );
+
+
+ while( true ){
+
+ // check for loop BREAK
+ if( position == root ){
+ _RTREE_PRINT(" CT2 position is root");
+ if( elimination_performed ){
+ root_elimination_performed = true;
+ }
+ break;
+ }
+
+ /*
+ pick randomly, let's say the 1st entry of the current node.
+ Search for this spatial area in the tree, and stop to the parent node.
+ Then find position of current node pointer, in the parent node.
+ */
+ /*
+ special case. if elimination due to children being underfull was performed AND
+ AND parent had only 1 record ,then this one record was removed.
+ */
+ if( position->children_nodes.size() > 0 ){
+ _RTREE_PRINT(" CT2.1 - 2 non leaf: find parent, P is parent");
+ // CT2.1 find parent. By definition it's nonleaf
+ find_result = find_parent( root, position->children_nodes[0].bounding_box, position);
+ }
+ else if( position->children_nodes.size() == 0
+ && position->children_leaves.size() == 0
+ && elimination_performed )
+ { // special case
+ _RTREE_PRINT(" CT2.1 - 2 special case: find parent, P is parent");
+ // CT2.1 find parent. By definition it's nonleaf
+ find_result = find_parent( root, special_case_bounding_box, position);
+ }
+ else{
+ _RTREE_PRINT(" CT2.1 - 2 leaf: find parent, P is parent");
+ // CT2.1 find parent. By definition it's nonleaf
+ find_result = find_parent( root, position->children_leaves[0].bounding_box, position);
+ }
+ // CT2.2 Let P be the parent of N
+ parent = find_result.first;
+
+
+ // parent is a non-leaf, by definition. Calculate "child_in_parent"
+ _RTREE_PRINT(" CT2.3 find in parent, position's record EN");
+ // CT2.3 Let EN be the N's entry in P
+ for( child_in_parent = 0; child_in_parent < parent->children_nodes.size(); child_in_parent++ ){
+ if( parent->children_nodes[ child_in_parent ].data == position){
+ _RTREE_PRINT(" child_in_parent: " << child_in_parent << " out of " << parent->children_nodes.size() << " (size)" );
+ break;
+ }
+ }
+
+ if( position->children_nodes.size() > 0 ){ // non leaf node: position
+ _RTREE_PRINT(" CT3 nonleaf: eliminate underfull node");
+ // CT3 Eliminate underfull node
+ if( position->children_nodes.size() < min_records ){
+ _RTREE_PRINT(" CT3.2 add N to Q");
+ // CT3.2 add N to set Q ( EN the record that points to N )
+ for(auto & children_node : position->children_nodes){
+ _RTREE_PRINT(" i " << i );
+ std::pair< RTreeRecord_NonLeaf, unsigned > t = std::make_pair( children_node, current_height-1);
+ Q_nonleaf_records.push_back( t );
+
+ }
+ special_case_bounding_box = parent->children_nodes[ child_in_parent ].bounding_box;
+
+ _RTREE_PRINT(" CT3.1 delete in parent, position's record EN");
+ // CT3.1 delete EN from P ( parent is by definition nonleaf )
+ if( remove_record_from_parent( parent, position ) ){ // TODO does erase, delete to the pointer ???
+ _RTREE_PRINT(" remove_record_from_parent error ");
+ }
+ elimination_performed = true;
+ }
+ else{
+ _RTREE_PRINT(" CT4 "); /// CT4) if not underfull
+ recalculate_bounding_box( parent, position, child_in_parent );
+ elimination_performed = false;
+ }
+
+ }
+ else{ // leaf node: position
+ _RTREE_PRINT(" CT3 leaf: eliminate underfull node");
+ // CT3 Eliminate underfull node
+ if( position->children_leaves.size() < min_records ){
+ _RTREE_PRINT(" CT3.2 add N to Q " << position->children_leaves.size() );
+ // CT3.2 add N to set Q
+ for(auto & children_leave : position->children_leaves){
+ _RTREE_PRINT(" i " << i );
+ Q_leaf_records.push_back( children_leave ); // TODO problem here
+ special_case_bounding_box = children_leave.bounding_box;
+ }
+
+ _RTREE_PRINT(" CT3.1 delete in parent, position's record EN");
+ // CT3.1 delete EN from P ( parent is by definition nonleaf )
+ if( remove_record_from_parent( parent, position ) ){
+ _RTREE_PRINT(" remove_record_from_parent error ");
+ }
+ elimination_performed = true;
+ }
+ else{
+ _RTREE_PRINT(" CT4 "); /// CT4) if not underfull
+ recalculate_bounding_box( parent, position, child_in_parent );
+ elimination_performed = false;
+ }
+ }
+ _RTREE_PRINT(" CT5 ");// CT5) move up one level in tree
+ position = parent;
+
+ current_height--;
+ }
+ // CT6: reinsert
+ _RTREE_PRINT(" ------ Q_leaf");
+ for( std::vector< RTreeRecord_Leaf >::iterator it = Q_leaf_records.begin(); it != Q_leaf_records.end(); ++it ){
+ _RTREE_PRINT(" leaf:" << (*it).data);
+ }
+ _RTREE_PRINT(" ------ Q_nonleaf");
+ for( std::vector< std::pair< RTreeRecord_NonLeaf, unsigned > >::iterator it = Q_nonleaf_records.begin(); it != Q_nonleaf_records.end(); ++it ){
+ _RTREE_PRINT(" ------- " << it->second );
+ _RTREE_PRINT_TREE( it->first.data, 0);
+ }
+
+ _RTREE_PRINT(" CT6 ");
+ for(auto & Q_leaf_record : Q_leaf_records){
+ insert( Q_leaf_record );
+ _RTREE_PRINT(" inserted leaf:" << (*it).data << " ------------");
+ _RTREE_PRINT_TREE( root, 0);
+ }
+
+
+ for(auto & Q_nonleaf_record : Q_nonleaf_records){
+ insert( RTreeRecord_Leaf() , true, Q_nonleaf_record.second, Q_nonleaf_record.first );
+ _RTREE_PRINT(" inserted nonleaf------------");
+ _RTREE_PRINT_TREE( root, 0);
+ // TODO this fake RTreeRecord_Leaf() looks stupid. find better way to do this ???
+ }
+
+ return root_elimination_performed;
+}
+
+
+/*
+given:
+- a parent
+- a child node
+- and the position of the child node in the parent
+recalculate the parent record's bounding box of the child, in order to ightly contain all entries of child
+
+NOTE! child must be indeed child of the parent, otherwise it screws up things. So find parent and child
+before calling this function
+*/
+void RTree::recalculate_bounding_box( RTreeNode* parent, RTreeNode* child, unsigned &child_in_parent ) {
+ if( child->children_nodes.size() > 0 ){
+ _RTREE_PRINT(" non-leaf: recalculate bounding box of parent "); // non leaf-node: child
+ parent->children_nodes[ child_in_parent ].bounding_box = Rect( child->children_nodes[0].bounding_box );
+ for( unsigned i=1; i < child->children_nodes.size(); i++ ){
+ parent->children_nodes[ child_in_parent ].bounding_box.unionWith( child->children_nodes[i].bounding_box );
+ }
+ }
+ else{
+ _RTREE_PRINT(" leaf: recalculate bounding box of parent "); // leaf-node: child
+ parent->children_nodes[ child_in_parent ].bounding_box = Rect( child->children_leaves[0].bounding_box );
+
+ for( unsigned i=1; i < child->children_leaves.size(); i++ ){
+ parent->children_nodes[ child_in_parent ].bounding_box.unionWith( child->children_leaves[i].bounding_box );
+ }
+ }
+}
+
+/*
+given:
+- a parent
+- a child node
+it removes the child record from the parent
+
+NOTE! child must be indeed child of the parent, otherwise it screws up things.
+So find parent and child before calling this function
+*/
+int RTree::remove_record_from_parent( RTreeNode* parent, RTreeNode* child ) {
+ _RTREE_PRINT( "remove_record_from_parent)" );
+ for( std::vector< RTreeRecord_NonLeaf >::iterator it = parent->children_nodes.begin(); it!=parent->children_nodes.end(); ++it ){
+ if( it->data == child ){
+ // delete element: implement step D2)
+ parent->children_nodes.erase( it );
+ return 0; // success
+ }
+ }
+ return 1; // failure
+}
+
+/*=============================================================================
+TODO update
+===============================================================================
+*/
+
+
+};
+
+/*
+ Local Variables:
+ mode:c++
+ c-file-style:"stroustrup"
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+ indent-tabs-mode:nil
+ fill-column:99
+ End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-29 21:57:03 +0000