diff options
Diffstat (limited to 'src/2geom/orphan-code')
| -rw-r--r-- | src/2geom/orphan-code/arc-length.cpp | 292 | ||||
| -rw-r--r-- | src/2geom/orphan-code/chebyshev.cpp | 126 | ||||
| -rw-r--r-- | src/2geom/orphan-code/intersection-by-bezier-clipping.cpp | 560 | ||||
| -rw-r--r-- | src/2geom/orphan-code/intersection-by-smashing.cpp | 349 | ||||
| -rw-r--r-- | src/2geom/orphan-code/nearestpoint.cpp | 405 | ||||
| -rw-r--r-- | src/2geom/orphan-code/redblack-toy.cpp | 327 | ||||
| -rw-r--r-- | src/2geom/orphan-code/redblacktree.cpp | 575 | ||||
| -rw-r--r-- | src/2geom/orphan-code/rtree.cpp | 1350 |
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 : |