From 62b84873d9ce273a3b7a59631c97a9b5bfe2dc06 Mon Sep 17 00:00:00 2001 From: johanengelen Date: Sat, 22 Mar 2008 21:18:07 +0000 Subject: [PATCH] update to latest 2geom --- src/2geom/angle.h | 100 +++++----- src/2geom/matrix.cpp | 36 ++-- src/2geom/matrix.h | 1 + src/2geom/poly-dk-solve.cpp | 128 ++++++------- src/2geom/poly-laguerre-solve.cpp | 294 +++++++++++++++--------------- src/2geom/sweep.cpp | 210 ++++++++++----------- 6 files changed, 387 insertions(+), 382 deletions(-) diff --git a/src/2geom/angle.h b/src/2geom/angle.h index 05364fb22..4d548ab49 100644 --- a/src/2geom/angle.h +++ b/src/2geom/angle.h @@ -1,50 +1,50 @@ -/** - * \file angle.h - * \brief Various trigoniometric helper functions - * - * Authors: - * Johan Engelen - * - * Copyright (C) 2007 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. - * - */ - -#ifndef LIB2GEOM_SEEN_ANGLE_H -#define LIB2GEOM_SEEN_ANGLE_H - -namespace Geom { - -#ifndef M_PI -# define M_PI 3.14159265358979323846 -#endif - -inline double deg_to_rad(double deg) { return deg*M_PI/180.0;} - -inline double rad_to_deg(double rad) { return rad*180.0/M_PI;} - -} - -#endif +/** + * \file angle.h + * \brief Various trigoniometric helper functions + * + * Authors: + * Johan Engelen + * + * Copyright (C) 2007 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. + * + */ + +#ifndef LIB2GEOM_SEEN_ANGLE_H +#define LIB2GEOM_SEEN_ANGLE_H + +namespace Geom { + +#ifndef M_PI +# define M_PI 3.14159265358979323846 +#endif + +inline double deg_to_rad(double deg) { return deg*M_PI/180.0;} + +inline double rad_to_deg(double rad) { return rad*180.0/M_PI;} + +} + +#endif diff --git a/src/2geom/matrix.cpp b/src/2geom/matrix.cpp index 6d7c772c0..f90bb6d42 100644 --- a/src/2geom/matrix.cpp +++ b/src/2geom/matrix.cpp @@ -108,9 +108,9 @@ void Matrix::setIdentity() { //TODO: use eps bool Matrix::isIdentity(Coord const eps) const { - return are_near(_c[0], 1.0) && are_near(_c[1], 0.0) && - are_near(_c[2], 0.0) && are_near(_c[3], 1.0) && - are_near(_c[4], 0.0) && are_near(_c[5], 0.0); + return are_near(_c[0], 1.0, eps) && are_near(_c[1], 0.0, eps) && + are_near(_c[2], 0.0, eps) && are_near(_c[3], 1.0, eps) && + are_near(_c[4], 0.0, eps) && are_near(_c[5], 0.0, eps); } /** Answers the question "Does this matrix perform a translation, and \em{only} a translation?" @@ -118,9 +118,9 @@ bool Matrix::isIdentity(Coord const eps) const { \return A bool representing yes/no. */ bool Matrix::isTranslation(Coord const eps) const { - return are_near(_c[0], 1.0) && are_near(_c[1], 0.0) && - are_near(_c[2], 0.0) && are_near(_c[3], 1.0) && - (!are_near(_c[4], 0.0) || !are_near(_c[5], 0.0)); + return are_near(_c[0], 1.0, eps) && are_near(_c[1], 0.0, eps) && + are_near(_c[2], 0.0, eps) && are_near(_c[3], 1.0, eps) && + (!are_near(_c[4], 0.0, eps) || !are_near(_c[5], 0.0, eps)); } /** Answers the question "Does this matrix perform a scale, and \em{only} a Scale?" @@ -128,9 +128,9 @@ bool Matrix::isTranslation(Coord const eps) const { \return A bool representing yes/no. */ bool Matrix::isScale(Coord const eps) const { - return !are_near(_c[0], 1.0) || !are_near(_c[3], 1.0) && //NOTE: these are the diags, and the next line opposite diags - are_near(_c[1], 0.0) && are_near(_c[2], 0.0) && - are_near(_c[4], 0.0) && are_near(_c[5], 0.0); + return !are_near(_c[0], 1.0, eps) || !are_near(_c[3], 1.0, eps) && //NOTE: these are the diags, and the next line opposite diags + are_near(_c[1], 0.0, eps) && are_near(_c[2], 0.0, eps) && + are_near(_c[4], 0.0, eps) && are_near(_c[5], 0.0, eps); } /** Answers the question "Does this matrix perform a uniform scale, and \em{only} a uniform scale?" @@ -138,9 +138,9 @@ bool Matrix::isScale(Coord const eps) const { \return A bool representing yes/no. */ bool Matrix::isUniformScale(Coord const eps) const { - return !are_near(_c[0], 1.0) && are_near(_c[0], _c[3]) && - are_near(_c[1], 0.0) && are_near(_c[2], 0.0) && - are_near(_c[4], 0.0) && are_near(_c[5], 0.0); + return !are_near(_c[0], 1.0, eps) && are_near(_c[0], _c[3], eps) && + are_near(_c[1], 0.0, eps) && are_near(_c[2], 0.0, eps) && + are_near(_c[4], 0.0, eps) && are_near(_c[5], 0.0, eps); } /** Answers the question "Does this matrix perform a rotation, and \em{only} a rotation?" @@ -148,13 +148,17 @@ bool Matrix::isUniformScale(Coord const eps) const { \return A bool representing yes/no. */ bool Matrix::isRotation(Coord const eps) const { - return !are_near(_c[0], _c[3]) && are_near(_c[1], -_c[2]) && - are_near(_c[4], 0.0) && are_near(_c[5], 0.0) && - are_near(_c[0]*_c[0] + _c[1]*_c[1], 1.0); + return !are_near(_c[0], _c[3], eps) && are_near(_c[1], -_c[2], eps) && + are_near(_c[4], 0.0, eps) && are_near(_c[5], 0.0, eps) && + are_near(_c[0]*_c[0] + _c[1]*_c[1], 1.0, eps); } bool Matrix::onlyScaleAndTranslation(Coord const eps) const { - return are_near(_c[0], _c[3]) && are_near(_c[1], 0) && are_near(_c[2], 0); + return are_near(_c[0], _c[3], eps) && are_near(_c[1], 0, eps) && are_near(_c[2], 0, eps); +} + +bool Matrix::isSingular(Coord const eps) const { + return are_near(det(), 0.0, eps); } bool Matrix::flips() const { diff --git a/src/2geom/matrix.h b/src/2geom/matrix.h index c39c99716..c9f244d62 100644 --- a/src/2geom/matrix.h +++ b/src/2geom/matrix.h @@ -91,6 +91,7 @@ class Matrix { bool isScale(double eps = EPSILON) const; bool isUniformScale(double eps = EPSILON) const; bool onlyScaleAndTranslation(double eps = EPSILON) const; + bool isSingular(double eps = EPSILON) const; bool flips() const; diff --git a/src/2geom/poly-dk-solve.cpp b/src/2geom/poly-dk-solve.cpp index fdc1cefe5..87d238f14 100644 --- a/src/2geom/poly-dk-solve.cpp +++ b/src/2geom/poly-dk-solve.cpp @@ -1,64 +1,64 @@ -#include "poly-dk-solve.h" -#include - -/*** implementation of the Durand-Kerner method. seems buggy*/ - -std::complex evalu(Poly const & p, std::complex x) { - std::complex result = 0; - std::complex xx = 1; - - for(unsigned i = 0; i < p.size(); i++) { - result += p[i]*xx; - xx *= x; - } - return result; -} - -std::vector > DK(Poly const & ply, const double tol) { - std::vector > roots; - const int N = ply.degree(); - - std::complex b(0.4, 0.9); - std::complex p = 1; - for(int i = 0; i < N; i++) { - roots.push_back(p); - p *= b; - } - assert(roots.size() == ply.degree()); - - double error = 0; - int i; - for( i = 0; i < 30; i++) { - error = 0; - for(int r_i = 0; r_i < N; r_i++) { - std::complex denom = 1; - std::complex R = roots[r_i]; - for(int d_i = 0; d_i < N; d_i++) { - if(r_i != d_i) - denom *= R-roots[d_i]; - } - assert(norm(denom) != 0); - std::complex dr = evalu(ply, R)/denom; - error += norm(dr); - roots[r_i] = R - dr; - } - /*std::copy(roots.begin(), roots.end(), std::ostream_iterator >(std::cout, ",\t")); - std::cout << std::endl;*/ - if(error < tol) - break; - } - //std::cout << error << ", " << i<< std::endl; - return roots; -} - - -/* - 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:encoding=utf-8:textwidth=99 : +#include "poly-dk-solve.h" +#include + +/*** implementation of the Durand-Kerner method. seems buggy*/ + +std::complex evalu(Poly const & p, std::complex x) { + std::complex result = 0; + std::complex xx = 1; + + for(unsigned i = 0; i < p.size(); i++) { + result += p[i]*xx; + xx *= x; + } + return result; +} + +std::vector > DK(Poly const & ply, const double tol) { + std::vector > roots; + const int N = ply.degree(); + + std::complex b(0.4, 0.9); + std::complex p = 1; + for(int i = 0; i < N; i++) { + roots.push_back(p); + p *= b; + } + assert(roots.size() == ply.degree()); + + double error = 0; + int i; + for( i = 0; i < 30; i++) { + error = 0; + for(int r_i = 0; r_i < N; r_i++) { + std::complex denom = 1; + std::complex R = roots[r_i]; + for(int d_i = 0; d_i < N; d_i++) { + if(r_i != d_i) + denom *= R-roots[d_i]; + } + assert(norm(denom) != 0); + std::complex dr = evalu(ply, R)/denom; + error += norm(dr); + roots[r_i] = R - dr; + } + /*std::copy(roots.begin(), roots.end(), std::ostream_iterator >(std::cout, ",\t")); + std::cout << std::endl;*/ + if(error < tol) + break; + } + //std::cout << error << ", " << i<< std::endl; + return roots; +} + + +/* + 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:encoding=utf-8:textwidth=99 : diff --git a/src/2geom/poly-laguerre-solve.cpp b/src/2geom/poly-laguerre-solve.cpp index 766f16eda..921ec3e3b 100644 --- a/src/2geom/poly-laguerre-solve.cpp +++ b/src/2geom/poly-laguerre-solve.cpp @@ -1,147 +1,147 @@ -#include "poly-laguerre-solve.h" -#include - -typedef std::complex cdouble; - -cdouble laguerre_internal_complex(Poly const & p, - double x0, - double tol, - bool & quad_root) { - cdouble a = 2*tol; - cdouble xk = x0; - double n = p.degree(); - quad_root = false; - const unsigned shuffle_rate = 10; - static double shuffle[] = {0, 0.5, 0.25, 0.75, 0.125, 0.375, 0.625, 0.875, 1.0}; - unsigned shuffle_counter = 0; - while(std::norm(a) > (tol*tol)) { - //std::cout << "xk = " << xk << std::endl; - cdouble b = p.back(); - cdouble d = 0, f = 0; - double err = abs(b); - double abx = abs(xk); - for(int j = p.size()-2; j >= 0; j--) { - f = xk*f + d; - d = xk*d + b; - b = xk*b + p[j]; - err = abs(b) + abx*err; - } - - err *= 1e-7; // magic epsilon for convergence, should be computed from tol - - cdouble px = b; - if(abs(b) < err) - return xk; - //if(std::norm(px) < tol*tol) - // return xk; - cdouble G = d / px; - cdouble H = G*G - f / px; - - //std::cout << "G = " << G << "H = " << H; - cdouble radicand = (n - 1)*(n*H-G*G); - //assert(radicand.real() > 0); - if(radicand.real() < 0) - quad_root = true; - //std::cout << "radicand = " << radicand << std::endl; - if(G.real() < 0) // here we try to maximise the denominator avoiding cancellation - a = - sqrt(radicand); - else - a = sqrt(radicand); - //std::cout << "a = " << a << std::endl; - a = n / (a + G); - //std::cout << "a = " << a << std::endl; - if(shuffle_counter % shuffle_rate == 0) - ;//a *= shuffle[shuffle_counter / shuffle_rate]; - xk -= a; - shuffle_counter++; - if(shuffle_counter >= 90) - break; - } - //std::cout << "xk = " << xk << std::endl; - return xk; -} - -double laguerre_internal(Poly const & p, - Poly const & pp, - Poly const & ppp, - double x0, - double tol, - bool & quad_root) { - double a = 2*tol; - double xk = x0; - double n = p.degree(); - quad_root = false; - while(a*a > (tol*tol)) { - //std::cout << "xk = " << xk << std::endl; - double px = p(xk); - if(px*px < tol*tol) - return xk; - double G = pp(xk) / px; - double H = G*G - ppp(xk) / px; - - //std::cout << "G = " << G << "H = " << H; - double radicand = (n - 1)*(n*H-G*G); - assert(radicand > 0); - //std::cout << "radicand = " << radicand << std::endl; - if(G < 0) // here we try to maximise the denominator avoiding cancellation - a = - sqrt(radicand); - else - a = sqrt(radicand); - //std::cout << "a = " << a << std::endl; - a = n / (a + G); - //std::cout << "a = " << a << std::endl; - xk -= a; - } - //std::cout << "xk = " << xk << std::endl; - return xk; -} - - -std::vector -laguerre(Poly p, const double tol) { - std::vector solutions; - //std::cout << "p = " << p << " = "; - while(p.size() > 1) - { - double x0 = 0; - bool quad_root = false; - cdouble sol = laguerre_internal_complex(p, x0, tol, quad_root); - //if(abs(sol) > 1) break; - Poly dvs; - if(quad_root) { - dvs.push_back((sol*conj(sol)).real()); - dvs.push_back(-(sol + conj(sol)).real()); - dvs.push_back(1.0); - //std::cout << "(" << dvs << ")"; - //solutions.push_back(sol); - //solutions.push_back(conj(sol)); - } else { - //std::cout << sol << std::endl; - dvs.push_back(-sol.real()); - dvs.push_back(1.0); - solutions.push_back(sol); - //std::cout << "(" << dvs << ")"; - } - Poly r; - p = divide(p, dvs, r); - //std::cout << r << std::endl; - } - return solutions; -} - -std::vector -laguerre_real_interval(Poly const & ply, - const double lo, const double hi, - const double tol) { -} - -/* - 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:encoding=utf-8:textwidth=99 : +#include "poly-laguerre-solve.h" +#include + +typedef std::complex cdouble; + +cdouble laguerre_internal_complex(Poly const & p, + double x0, + double tol, + bool & quad_root) { + cdouble a = 2*tol; + cdouble xk = x0; + double n = p.degree(); + quad_root = false; + const unsigned shuffle_rate = 10; + static double shuffle[] = {0, 0.5, 0.25, 0.75, 0.125, 0.375, 0.625, 0.875, 1.0}; + unsigned shuffle_counter = 0; + while(std::norm(a) > (tol*tol)) { + //std::cout << "xk = " << xk << std::endl; + cdouble b = p.back(); + cdouble d = 0, f = 0; + double err = abs(b); + double abx = abs(xk); + for(int j = p.size()-2; j >= 0; j--) { + f = xk*f + d; + d = xk*d + b; + b = xk*b + p[j]; + err = abs(b) + abx*err; + } + + err *= 1e-7; // magic epsilon for convergence, should be computed from tol + + cdouble px = b; + if(abs(b) < err) + return xk; + //if(std::norm(px) < tol*tol) + // return xk; + cdouble G = d / px; + cdouble H = G*G - f / px; + + //std::cout << "G = " << G << "H = " << H; + cdouble radicand = (n - 1)*(n*H-G*G); + //assert(radicand.real() > 0); + if(radicand.real() < 0) + quad_root = true; + //std::cout << "radicand = " << radicand << std::endl; + if(G.real() < 0) // here we try to maximise the denominator avoiding cancellation + a = - sqrt(radicand); + else + a = sqrt(radicand); + //std::cout << "a = " << a << std::endl; + a = n / (a + G); + //std::cout << "a = " << a << std::endl; + if(shuffle_counter % shuffle_rate == 0) + ;//a *= shuffle[shuffle_counter / shuffle_rate]; + xk -= a; + shuffle_counter++; + if(shuffle_counter >= 90) + break; + } + //std::cout << "xk = " << xk << std::endl; + return xk; +} + +double laguerre_internal(Poly const & p, + Poly const & pp, + Poly const & ppp, + double x0, + double tol, + bool & quad_root) { + double a = 2*tol; + double xk = x0; + double n = p.degree(); + quad_root = false; + while(a*a > (tol*tol)) { + //std::cout << "xk = " << xk << std::endl; + double px = p(xk); + if(px*px < tol*tol) + return xk; + double G = pp(xk) / px; + double H = G*G - ppp(xk) / px; + + //std::cout << "G = " << G << "H = " << H; + double radicand = (n - 1)*(n*H-G*G); + assert(radicand > 0); + //std::cout << "radicand = " << radicand << std::endl; + if(G < 0) // here we try to maximise the denominator avoiding cancellation + a = - sqrt(radicand); + else + a = sqrt(radicand); + //std::cout << "a = " << a << std::endl; + a = n / (a + G); + //std::cout << "a = " << a << std::endl; + xk -= a; + } + //std::cout << "xk = " << xk << std::endl; + return xk; +} + + +std::vector +laguerre(Poly p, const double tol) { + std::vector solutions; + //std::cout << "p = " << p << " = "; + while(p.size() > 1) + { + double x0 = 0; + bool quad_root = false; + cdouble sol = laguerre_internal_complex(p, x0, tol, quad_root); + //if(abs(sol) > 1) break; + Poly dvs; + if(quad_root) { + dvs.push_back((sol*conj(sol)).real()); + dvs.push_back(-(sol + conj(sol)).real()); + dvs.push_back(1.0); + //std::cout << "(" << dvs << ")"; + //solutions.push_back(sol); + //solutions.push_back(conj(sol)); + } else { + //std::cout << sol << std::endl; + dvs.push_back(-sol.real()); + dvs.push_back(1.0); + solutions.push_back(sol); + //std::cout << "(" << dvs << ")"; + } + Poly r; + p = divide(p, dvs, r); + //std::cout << r << std::endl; + } + return solutions; +} + +std::vector +laguerre_real_interval(Poly const & ply, + const double lo, const double hi, + const double tol) { +} + +/* + 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:encoding=utf-8:textwidth=99 : diff --git a/src/2geom/sweep.cpp b/src/2geom/sweep.cpp index b9ef71b5c..c0e1c3997 100644 --- a/src/2geom/sweep.cpp +++ b/src/2geom/sweep.cpp @@ -1,105 +1,105 @@ -#include "sweep.h" - -#include - -namespace Geom { - -std::vector > sweep_bounds(std::vector rs) { - std::vector events; events.reserve(rs.size()*2); - std::vector > pairs(rs.size()); - - for(unsigned i = 0; i < rs.size(); i++) { - events.push_back(Event(rs[i].left(), i, false)); - events.push_back(Event(rs[i].right(), i, true)); - } - std::sort(events.begin(), events.end()); - - std::vector open; - for(unsigned i = 0; i < events.size(); i++) { - unsigned ix = events[i].ix; - if(events[i].closing) { - std::vector::iterator iter = std::find(open.begin(), open.end(), ix); - //if(iter != open.end()) - open.erase(iter); - } else { - for(unsigned j = 0; j < open.size(); j++) { - unsigned jx = open[j]; - if(rs[jx][Y].intersects(rs[ix][Y])) { - pairs[jx].push_back(ix); - } - } - open.push_back(ix); - } - } - return pairs; -} - -std::vector > sweep_bounds(std::vector a, std::vector b) { - std::vector > pairs(a.size()); - if(a.empty() || b.empty()) return pairs; - std::vector events[2]; - events[0].reserve(a.size()*2); - events[1].reserve(b.size()*2); - - for(unsigned n = 0; n < 2; n++) { - unsigned sz = n ? b.size() : a.size(); - events[n].reserve(sz*2); - for(unsigned i = 0; i < sz; i++) { - events[n].push_back(Event(n ? b[i].left() : a[i].left(), i, false)); - events[n].push_back(Event(n ? b[i].right() : a[i].right(), i, true)); - } - std::sort(events[n].begin(), events[n].end()); - } - - std::vector open[2]; - bool n = events[1].front() < events[0].front(); - for(unsigned i[] = {0,0}; i[n] < events[n].size();) { - unsigned ix = events[n][i[n]].ix; - bool closing = events[n][i[n]].closing; - //std::cout << n << "[" << ix << "] - " << (closing ? "closer" : "opener") << "\n"; - if(closing) { - open[n].erase(std::find(open[n].begin(), open[n].end(), ix)); - } else { - if(n) { - //n = 1 - //opening a B, add to all open a - for(unsigned j = 0; j < open[0].size(); j++) { - unsigned jx = open[0][j]; - if(a[jx][Y].intersects(b[ix][Y])) { - pairs[jx].push_back(ix); - } - } - } else { - //n = 0 - //opening an A, add all open b - for(unsigned j = 0; j < open[1].size(); j++) { - unsigned jx = open[1][j]; - if(b[jx][Y].intersects(a[ix][Y])) { - pairs[ix].push_back(jx); - } - } - } - open[n].push_back(ix); - } - i[n]++; - if(i[n]>=events[n].size()) {break;} - n = (events[!n][i[!n]] < events[n][i[n]]) ? !n : n; - } - return pairs; -} - -//Fake cull, until the switch to the real sweep is made. -std::vector > fake_cull(unsigned a, unsigned b) { - std::vector > ret; - - std::vector all; - for(unsigned j = 0; j < b; j++) - all.push_back(j); - - for(unsigned i = 0; i < a; i++) - ret.push_back(all); - - return ret; -} - -} +#include "sweep.h" + +#include + +namespace Geom { + +std::vector > sweep_bounds(std::vector rs) { + std::vector events; events.reserve(rs.size()*2); + std::vector > pairs(rs.size()); + + for(unsigned i = 0; i < rs.size(); i++) { + events.push_back(Event(rs[i].left(), i, false)); + events.push_back(Event(rs[i].right(), i, true)); + } + std::sort(events.begin(), events.end()); + + std::vector open; + for(unsigned i = 0; i < events.size(); i++) { + unsigned ix = events[i].ix; + if(events[i].closing) { + std::vector::iterator iter = std::find(open.begin(), open.end(), ix); + //if(iter != open.end()) + open.erase(iter); + } else { + for(unsigned j = 0; j < open.size(); j++) { + unsigned jx = open[j]; + if(rs[jx][Y].intersects(rs[ix][Y])) { + pairs[jx].push_back(ix); + } + } + open.push_back(ix); + } + } + return pairs; +} + +std::vector > sweep_bounds(std::vector a, std::vector b) { + std::vector > pairs(a.size()); + if(a.empty() || b.empty()) return pairs; + std::vector events[2]; + events[0].reserve(a.size()*2); + events[1].reserve(b.size()*2); + + for(unsigned n = 0; n < 2; n++) { + unsigned sz = n ? b.size() : a.size(); + events[n].reserve(sz*2); + for(unsigned i = 0; i < sz; i++) { + events[n].push_back(Event(n ? b[i].left() : a[i].left(), i, false)); + events[n].push_back(Event(n ? b[i].right() : a[i].right(), i, true)); + } + std::sort(events[n].begin(), events[n].end()); + } + + std::vector open[2]; + bool n = events[1].front() < events[0].front(); + for(unsigned i[] = {0,0}; i[n] < events[n].size();) { + unsigned ix = events[n][i[n]].ix; + bool closing = events[n][i[n]].closing; + //std::cout << n << "[" << ix << "] - " << (closing ? "closer" : "opener") << "\n"; + if(closing) { + open[n].erase(std::find(open[n].begin(), open[n].end(), ix)); + } else { + if(n) { + //n = 1 + //opening a B, add to all open a + for(unsigned j = 0; j < open[0].size(); j++) { + unsigned jx = open[0][j]; + if(a[jx][Y].intersects(b[ix][Y])) { + pairs[jx].push_back(ix); + } + } + } else { + //n = 0 + //opening an A, add all open b + for(unsigned j = 0; j < open[1].size(); j++) { + unsigned jx = open[1][j]; + if(b[jx][Y].intersects(a[ix][Y])) { + pairs[ix].push_back(jx); + } + } + } + open[n].push_back(ix); + } + i[n]++; + if(i[n]>=events[n].size()) {break;} + n = (events[!n][i[!n]] < events[n][i[n]]) ? !n : n; + } + return pairs; +} + +//Fake cull, until the switch to the real sweep is made. +std::vector > fake_cull(unsigned a, unsigned b) { + std::vector > ret; + + std::vector all; + for(unsigned j = 0; j < b; j++) + all.push_back(j); + + for(unsigned i = 0; i < a; i++) + ret.push_back(all); + + return ret; +} + +} -- 2.30.2