1
/** root finding for sbasis functions.
2
* Copyright 2006 N Hurst
3
* Copyright 2007 JF Barraud
5
* It is more efficient to find roots of f(t) = c_0, c_1, ... all at once, rather than iterating.
8
* multi-roots using bernstein method, one approach would be:
10
take median and find roots of that
11
whenever a segment lies entirely on one side of the median,
12
find the median of the half and recurse.
14
in essence we are implementing quicksort on a continuous function
16
* the gsl poly roots finder is faster than bernstein too, but we don't use it for 3 reasons:
18
a) it requires convertion to poly, which is numerically unstable
20
b) it requires gsl (which is currently not a dependency, and would bring in a whole slew of unrelated stuff)
22
c) it finds all roots, even complex ones. We don't want to accidently treat a nearly real root as a real root
24
From memory gsl poly roots was about 10 times faster than bernstein in the case where all the roots
25
are in [0,1] for polys of order 5. I spent some time working out whether eigenvalue root finding
26
could be done directly in sbasis space, but the maths was too hard for me. -- njh
28
jfbarraud: eigenvalue root finding could be done directly in sbasis space ?
30
njh: I don't know, I think it should. You would make a matrix whose characteristic polynomial was
31
correct, but do it by putting the sbasis terms in the right spots in the matrix. normal eigenvalue
32
root finding makes a matrix that is a diagonal + a row along the top. This matrix has the property
33
that its characteristic poly is just the poly whose coefficients are along the top row.
35
Now an sbasis function is a linear combination of the poly coeffs. So it seems to me that you
36
should be able to put the sbasis coeffs directly into a matrix in the right spots so that the
37
characteristic poly is the sbasis. You'll still have problems b) and c).
39
We might be able to lift an eigenvalue solver and include that directly into 2geom. Eigenvalues
40
also allow you to find intersections of multiple curves but require solving n*m x n*m matrices.
47
#include <2geom/sbasis.h>
48
#include <2geom/sbasis-to-bezier.h>
49
#include <2geom/solver.h>
55
/** Find the smallest interval that bounds a
56
\param a sbasis function
62
OptInterval bounds_exact(SBasisOf<double> const &a) {
63
Interval result = Interval(a.at0(), a.at1());
64
SBasisOf<double> df = derivative(a);
65
vector<double>extrema = roots(df);
66
for (unsigned i=0; i<extrema.size(); i++){
67
result.extendTo(a(extrema[i]));
72
OptInterval bounds_exact(SBasis const &a) {
73
Interval result = Interval(a.at0(), a.at1());
74
SBasis df = derivative(a);
75
vector<double>extrema = roots(df);
76
for (unsigned i=0; i<extrema.size(); i++){
77
result.extendTo(a(extrema[i]));
83
/** Find a small interval that bounds a
84
\param a sbasis function
88
// I have no idea how this works, some clever bounding argument by jfb.
90
OptInterval bounds_fast(const SBasisOf<double> &sb, int order) {
92
OptInterval bounds_fast(const SBasis &sb, int order) {
94
Interval res(0,0); // an empty sbasis is 0.
96
for(int j = sb.size()-1; j>=order; j--) {
102
if (v<0) t = ((b-a)/v+1)*0.5;
103
if (v>=0 || t<0 || t>1) {
104
res[0] = std::min(a,b);
106
res[0]=lerp(t, a+v*t, b);
110
if (v>0) t = ((b-a)/v+1)*0.5;
111
if (v<=0 || t<0 || t>1) {
112
res[1] = std::max(a,b);
114
res[1]=lerp(t, a+v*t, b);
117
if (order>0) res*=pow(.25,order);
121
/** Find a small interval that bounds a(t) for t in i to order order
122
\param sb sbasis function
123
\param i domain interval
124
\param order number of terms
129
OptInterval bounds_local(const SBasisOf<double> &sb, const OptInterval &i, int order) {
131
OptInterval bounds_local(const SBasis &sb, const OptInterval &i, int order) {
133
double t0=i->min(), t1=i->max(), lo=0., hi=0.;
134
for(int j = sb.size()-1; j>=order; j--) {
139
if (lo<0) t = ((b-a)/lo+1)*0.5;
140
if (lo>=0 || t<t0 || t>t1) {
141
lo = std::min(a*(1-t0)+b*t0+lo*t0*(1-t0),a*(1-t1)+b*t1+lo*t1*(1-t1));
143
lo = lerp(t, a+lo*t, b);
146
if (hi>0) t = ((b-a)/hi+1)*0.5;
147
if (hi<=0 || t<t0 || t>t1) {
148
hi = std::max(a*(1-t0)+b*t0+hi*t0*(1-t0),a*(1-t1)+b*t1+hi*t1*(1-t1));
150
hi = lerp(t, a+hi*t, b);
153
Interval res = Interval(lo,hi);
154
if (order>0) res*=pow(.25,order);
158
//-- multi_roots ------------------------------------
159
// goal: solve f(t)=c for several c at once.
160
/* algo: -compute f at both ends of the given segment [a,b].
161
-compute bounds m<df(t)<M for df on the segment.
162
let c and C be the levels below and above f(a):
163
going from f(a) down to c with slope m takes at least time (f(a)-c)/m
164
going from f(a) up to C with slope M takes at least time (C-f(a))/M
165
From this we conclude there are no roots before a'=a+min((f(a)-c)/m,(C-f(a))/M).
166
Do the same for b: compute some b' such that there are no roots in (b',b].
167
-if [a',b'] is not empty, repeat the process with [a',(a'+b')/2] and [(a'+b')/2,b'].
168
unfortunately, extra care is needed about rounding errors, and also to avoid the repetition of roots,
169
making things tricky and unpleasant...
171
//TODO: Make sure the code is "rounding-errors proof" and take care about repetition of roots!
174
static int upper_level(vector<double> const &levels,double x,double tol=0.){
175
return(upper_bound(levels.begin(),levels.end(),x-tol)-levels.begin());
179
static void multi_roots_internal(SBasis const &f,
182
static void multi_roots_internal(SBasis const &f,
185
std::vector<double> const &levels,
186
std::vector<std::vector<double> > &roots,
196
idx=upper_level(levels,0,vtol);
197
if (idx<(int)levels.size()&&fabs(levels.at(idx))<=vtol){
198
roots[idx].push_back(a);
199
roots[idx].push_back(b);
205
// int idxa=upper_level(levels,fa);
206
// int idxb=upper_level(levels,fb);
208
// if (fa==levels[idxa]){
214
// int idx_min=std::min(idxa,idxb);
215
// int idx_max=std::max(idxa,idxb);
216
// if (idx_max==levels.size()) idx_max-=1;
217
// for(int i=idx_min;i<=idx_max; i++){
218
// double t=a+(b-a)*(levels[i]-fa)/(fb-fa);
219
// if(a<t&&t<b) roots[t]=i;
224
//TODO: use different tol for t and f ?
225
//TODO: unsigned idx ? (remove int casts when fixed)
226
int idx=std::min(upper_level(levels,fa,vtol),upper_level(levels,fb,vtol));
227
if (idx==(int)levels.size()) idx-=1;
228
double c=levels.at(idx);
229
if((fa-c)*(fb-c)<=0||fabs(fa-c)<vtol||fabs(fb-c)<vtol){
230
roots[idx].push_back((a+b)/2);
235
int idxa=upper_level(levels,fa,vtol);
236
int idxb=upper_level(levels,fb,vtol);
238
Interval bs = *bounds_local(df,Interval(a,b));
240
//first times when a level (higher or lower) can be reached from a or b.
241
double ta_hi,tb_hi,ta_lo,tb_lo;
242
ta_hi=ta_lo=b+1;//default values => no root there.
243
tb_hi=tb_lo=a-1;//default values => no root there.
245
if (idxa<(int)levels.size() && fabs(fa-levels.at(idxa))<vtol){//a can be considered a root.
247
roots[idxa].push_back(a);
250
if (bs.max()>0 && idxa<(int)levels.size())
251
ta_hi=a+(levels.at(idxa )-fa)/bs.max();
252
if (bs.min()<0 && idxa>0)
253
ta_lo=a+(levels.at(idxa-1)-fa)/bs.min();
255
if (idxb<(int)levels.size() && fabs(fb-levels.at(idxb))<vtol){//b can be considered a root.
257
roots[idxb].push_back(b);
260
if (bs.min()<0 && idxb<(int)levels.size())
261
tb_hi=b+(levels.at(idxb )-fb)/bs.min();
262
if (bs.max()>0 && idxb>0)
263
tb_lo=b+(levels.at(idxb-1)-fb)/bs.max();
267
t0=std::min(ta_hi,ta_lo);
268
t1=std::max(tb_hi,tb_lo);
269
//hum, rounding errors frighten me! so I add this +tol...
270
if (t0>t1+htol) return;//no root here.
272
if (fabs(t1-t0)<htol){
273
multi_roots_internal(f,df,levels,roots,htol,vtol,t0,f(t0),t1,f(t1));
275
double t,t_left,t_right,ft,ft_left,ft_right;
276
t_left =t_right =t =(t0+t1)/2;
277
ft_left=ft_right=ft=f(t);
278
int idx=upper_level(levels,ft,vtol);
279
if (idx<(int)levels.size() && fabs(ft-levels.at(idx))<vtol){//t can be considered a root.
280
roots[idx].push_back(t);
281
//we do not want to count it twice (from the left and from the right)
287
multi_roots_internal(f,df,levels,roots,htol,vtol,t0 ,f(t0) ,t_left,ft_left);
288
multi_roots_internal(f,df,levels,roots,htol,vtol,t_right,ft_right,t1 ,f(t1) );
292
/** Solve f(t)=c for several c at once.
293
\param f sbasis function
294
\param levels vector of 'y' values
296
\param a, b left and right bounds
297
\returns a vector of vectors, one for each y giving roots
299
Effectively computes:
300
results = roots(f(y_i)) for all y_i
302
* algo: -compute f at both ends of the given segment [a,b].
303
-compute bounds m<df(t)<M for df on the segment.
304
let c and C be the levels below and above f(a):
305
going from f(a) down to c with slope m takes at least time (f(a)-c)/m
306
going from f(a) up to C with slope M takes at least time (C-f(a))/M
307
From this we conclude there are no roots before a'=a+min((f(a)-c)/m,(C-f(a))/M).
308
Do the same for b: compute some b' such that there are no roots in (b',b].
309
-if [a',b'] is not empty, repeat the process with [a',(a'+b')/2] and [(a'+b')/2,b'].
310
unfortunately, extra care is needed about rounding errors, and also to avoid the repetition of roots,
311
making things tricky and unpleasant...
313
TODO: Make sure the code is "rounding-errors proof" and take care about repetition of roots!
315
std::vector<std::vector<double> > multi_roots(SBasis const &f,
316
std::vector<double> const &levels,
322
std::vector<std::vector<double> > roots(levels.size(), std::vector<double>());
324
SBasis df=derivative(f);
325
multi_roots_internal(f,df,levels,roots,htol,vtol,a,f(a),b,f(b));
329
//-------------------------------------
332
void subdiv_sbasis(SBasis const & s,
333
std::vector<double> & roots,
334
double left, double right) {
335
OptInterval bs = bounds_fast(s);
336
if(!bs || bs->min() > 0 || bs->max() < 0)
337
return; // no roots here
338
if(s.tailError(1) < 1e-7) {
339
double t = s[0][0] / (s[0][0] - s[0][1]);
340
roots.push_back(left*(1-t) + t*right);
343
double middle = (left + right)/2;
344
subdiv_sbasis(compose(s, Linear(0, 0.5)), roots, left, middle);
345
subdiv_sbasis(compose(s, Linear(0.5, 1.)), roots, middle, right);
348
// It is faster to use the bernstein root finder for small degree polynomials (<100?.
350
std::vector<double> roots1(SBasis const & s) {
351
std::vector<double> res;
352
double d = s[0][0] - s[0][1];
354
double r = s[0][0] / d;
361
/** Find all t s.t s(t) = 0
362
\param a sbasis function
363
\returns vector of zeros (roots)
366
std::vector<double> roots(SBasis const & s) {
369
return std::vector<double>();
375
sbasis_to_bezier(bz, s);
386
c-file-style:"stroustrup"
387
c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
392
// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
added
removed
inkscape-0.47pre1/src/2geom/sbasis-roots.cpp
