From b6d9822b29b878acf90e7900b73c56d32ce75ae4 Mon Sep 17 00:00:00 2001 From: johanengelen Date: Sun, 23 Dec 2007 19:34:02 +0000 Subject: [PATCH] update to latest 2geom (fixed 2 bugs :)) --- src/2geom/path.h | 3 +- src/2geom/poly-dk-solve.cpp | 128 ++++++------- src/2geom/poly-laguerre-solve.cpp | 294 +++++++++++++++--------------- src/2geom/sbasis-to-bezier.cpp | 10 +- src/2geom/sweep.cpp | 208 ++++++++++----------- 5 files changed, 322 insertions(+), 321 deletions(-) diff --git a/src/2geom/path.h b/src/2geom/path.h index 988babe3e..3c536c1d8 100644 --- a/src/2geom/path.h +++ b/src/2geom/path.h @@ -492,7 +492,8 @@ public: Piecewise > ret; ret.push_cut(0); unsigned i = 1; - for(const_iterator it = begin(); it != end_default(); ++it, i++) { + // ignore that path is closed or open. pw> is always open. + for(const_iterator it = begin(); it != end(); ++it, i++) { ret.push(it->toSBasis(), i); } return ret; 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/sbasis-to-bezier.cpp b/src/2geom/sbasis-to-bezier.cpp index 7283dd87d..415b30d89 100644 --- a/src/2geom/sbasis-to-bezier.cpp +++ b/src/2geom/sbasis-to-bezier.cpp @@ -212,15 +212,15 @@ path_from_piecewise(Geom::Piecewise > const &B, double to for(unsigned i = 0; ; i++) { if(i+1 == B.size() || !are_near(B[i+1].at0(), B[i].at1(), tol)) { //start of a new path + if(are_near(start, B[i].at1()) && sbasis_size(B[i]) <= 1) { + //last line seg already there + goto no_add; + } + build_from_sbasis(pb, B[i], tol); if(are_near(start, B[i].at1())) { //it's closed pb.closePath(); - if(sbasis_size(B[i]) <= 1) { - //last line seg already there - goto no_add; - } } - build_from_sbasis(pb, B[i], tol); no_add: if(i+1 >= B.size()) break; start = B[i+1].at0(); diff --git a/src/2geom/sweep.cpp b/src/2geom/sweep.cpp index 08674ab2f..96104ea90 100644 --- a/src/2geom/sweep.cpp +++ b/src/2geom/sweep.cpp @@ -1,104 +1,104 @@ -#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]++; - 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]++; + 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