summaryrefslogtreecommitdiffstats
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/2geom/path.h3
-rw-r--r--src/2geom/poly-dk-solve.cpp128
-rw-r--r--src/2geom/poly-laguerre-solve.cpp294
-rw-r--r--src/2geom/sbasis-to-bezier.cpp10
-rw-r--r--src/2geom/sweep.cpp208
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<D2<SBasis> > 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<d2<>> 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 <iterator>
-
-/*** implementation of the Durand-Kerner method. seems buggy*/
-
-std::complex<double> evalu(Poly const & p, std::complex<double> x) {
- std::complex<double> result = 0;
- std::complex<double> xx = 1;
-
- for(unsigned i = 0; i < p.size(); i++) {
- result += p[i]*xx;
- xx *= x;
- }
- return result;
-}
-
-std::vector<std::complex<double> > DK(Poly const & ply, const double tol) {
- std::vector<std::complex<double> > roots;
- const int N = ply.degree();
-
- std::complex<double> b(0.4, 0.9);
- std::complex<double> 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<double> denom = 1;
- std::complex<double> 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<double> dr = evalu(ply, R)/denom;
- error += norm(dr);
- roots[r_i] = R - dr;
- }
- /*std::copy(roots.begin(), roots.end(), std::ostream_iterator<std::complex<double> >(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 <iterator>
+
+/*** implementation of the Durand-Kerner method. seems buggy*/
+
+std::complex<double> evalu(Poly const & p, std::complex<double> x) {
+ std::complex<double> result = 0;
+ std::complex<double> xx = 1;
+
+ for(unsigned i = 0; i < p.size(); i++) {
+ result += p[i]*xx;
+ xx *= x;
+ }
+ return result;
+}
+
+std::vector<std::complex<double> > DK(Poly const & ply, const double tol) {
+ std::vector<std::complex<double> > roots;
+ const int N = ply.degree();
+
+ std::complex<double> b(0.4, 0.9);
+ std::complex<double> 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<double> denom = 1;
+ std::complex<double> 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<double> dr = evalu(ply, R)/denom;
+ error += norm(dr);
+ roots[r_i] = R - dr;
+ }
+ /*std::copy(roots.begin(), roots.end(), std::ostream_iterator<std::complex<double> >(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 <iterator>
-
-typedef std::complex<double> 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<cdouble >
-laguerre(Poly p, const double tol) {
- std::vector<cdouble > 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<double>
-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 <iterator>
+
+typedef std::complex<double> 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<cdouble >
+laguerre(Poly p, const double tol) {
+ std::vector<cdouble > 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<double>
+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<Geom::D2<Geom::SBasis> > 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 <algorithm>
-
-namespace Geom {
-
-std::vector<std::vector<unsigned> > sweep_bounds(std::vector<Rect> rs) {
- std::vector<Event> events; events.reserve(rs.size()*2);
- std::vector<std::vector<unsigned> > 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<unsigned> open;
- for(unsigned i = 0; i < events.size(); i++) {
- unsigned ix = events[i].ix;
- if(events[i].closing) {
- std::vector<unsigned>::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<std::vector<unsigned> > sweep_bounds(std::vector<Rect> a, std::vector<Rect> b) {
- std::vector<std::vector<unsigned> > pairs(a.size());
- if(a.empty() || b.empty()) return pairs;
- std::vector<Event> 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<unsigned> 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<std::vector<unsigned> > fake_cull(unsigned a, unsigned b) {
- std::vector<std::vector<unsigned> > ret;
-
- std::vector<unsigned> 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 <algorithm>
+
+namespace Geom {
+
+std::vector<std::vector<unsigned> > sweep_bounds(std::vector<Rect> rs) {
+ std::vector<Event> events; events.reserve(rs.size()*2);
+ std::vector<std::vector<unsigned> > 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<unsigned> open;
+ for(unsigned i = 0; i < events.size(); i++) {
+ unsigned ix = events[i].ix;
+ if(events[i].closing) {
+ std::vector<unsigned>::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<std::vector<unsigned> > sweep_bounds(std::vector<Rect> a, std::vector<Rect> b) {
+ std::vector<std::vector<unsigned> > pairs(a.size());
+ if(a.empty() || b.empty()) return pairs;
+ std::vector<Event> 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<unsigned> 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<std::vector<unsigned> > fake_cull(unsigned a, unsigned b) {
+ std::vector<std::vector<unsigned> > ret;
+
+ std::vector<unsigned> 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;
+}
+
+}
generated by cgit v1.2.3 (git 2.25.1) at 2026-05-19 06:31:19 +0000