95fd0cf3ba8e45c19da4b878537ad7f84ea5576e
1 /** root finding for sbasis functions.
2 * Copyright 2006 N Hurst
3 * Copyright 2007 JF Barraud
4 *
5 * It is more efficient to find roots of f(t) = c_0, c_1, ... all at once, rather than iterating.
6 *
7 * Todo/think about:
8 * multi-roots using bernstein method, one approach would be:
9 sort c
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.
42 **/
44 #include <cmath>
45 #include <map>
47 #include <2geom/sbasis.h>
48 #include <2geom/sbasis-to-bezier.h>
49 #include <2geom/solver.h>
51 using namespace std;
53 namespace Geom{
55 /** Find the smallest interval that bounds a
56 \param a sbasis function
57 \returns inteval
59 */
61 #ifdef USE_SBASIS_OF
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]));
68 }
69 return result;
70 }
71 #else
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]));
78 }
79 return result;
80 }
81 #endif
83 /** Find a small interval that bounds a
84 \param a sbasis function
85 \returns inteval
87 */
88 // I have no idea how this works, some clever bounding argument by jfb.
89 #ifdef USE_SBASIS_OF
90 OptInterval bounds_fast(const SBasisOf<double> &sb, int order) {
91 #else
92 OptInterval bounds_fast(const SBasis &sb, int order) {
93 #endif
94 Interval res(0,0); // an empty sbasis is 0.
96 for(int j = sb.size()-1; j>=order; j--) {
97 double a=sb[j][0];
98 double b=sb[j][1];
100 double v, t = 0;
101 v = res[0];
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);
105 }else{
106 res[0]=lerp(t, a+v*t, b);
107 }
109 v = res[1];
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);
113 }else{
114 res[1]=lerp(t, a+v*t, b);
115 }
116 }
117 if (order>0) res*=pow(.25,order);
118 return res;
119 }
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
125 \return interval
127 */
128 #ifdef USE_SBASIS_OF
129 OptInterval bounds_local(const SBasisOf<double> &sb, const OptInterval &i, int order) {
130 #else
131 OptInterval bounds_local(const SBasis &sb, const OptInterval &i, int order) {
132 #endif
133 double t0=i->min(), t1=i->max(), lo=0., hi=0.;
134 for(int j = sb.size()-1; j>=order; j--) {
135 double a=sb[j][0];
136 double b=sb[j][1];
138 double t = 0;
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));
142 }else{
143 lo = lerp(t, a+lo*t, b);
144 }
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));
149 }else{
150 hi = lerp(t, a+hi*t, b);
151 }
152 }
153 Interval res = Interval(lo,hi);
154 if (order>0) res*=pow(.25,order);
155 return res;
156 }
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...
170 */
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());
176 }
178 #ifdef USE_SBASIS_OF
179 static void multi_roots_internal(SBasis const &f,
180 SBasis const &df,
181 #else
182 static void multi_roots_internal(SBasis const &f,
183 SBasis const &df,
184 #endif
185 std::vector<double> const &levels,
186 std::vector<std::vector<double> > &roots,
187 double htol,
188 double vtol,
189 double a,
190 double fa,
191 double b,
192 double fb){
194 if (f.size()==0){
195 int idx;
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);
200 }
201 return;
202 }
203 ////usefull?
204 // if (f.size()==1){
205 // int idxa=upper_level(levels,fa);
206 // int idxb=upper_level(levels,fb);
207 // if (fa==fb){
208 // if (fa==levels[idxa]){
209 // roots[a]=idxa;
210 // roots[b]=idxa;
211 // }
212 // return;
213 // }
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;
220 // }
221 // return;
222 // }
223 if ((b-a)<htol){
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);
231 }
232 return;
233 }
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.
246 //ta_hi=ta_lo=a;
247 roots[idxa].push_back(a);
248 ta_hi=ta_lo=a+htol;
249 }else{
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();
254 }
255 if (idxb<(int)levels.size() && fabs(fb-levels.at(idxb))<vtol){//b can be considered a root.
256 //tb_hi=tb_lo=b;
257 roots[idxb].push_back(b);
258 tb_hi=tb_lo=b-htol;
259 }else{
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();
264 }
266 double t0,t1;
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));
274 }else{
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)
282 t_left =t-htol/2;
283 t_right=t+htol/2;
284 ft_left =f(t_left);
285 ft_right=f(t_right);
286 }
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) );
289 }
290 }
292 /** Solve f(t)=c for several c at once.
293 \param f sbasis function
294 \param levels vector of 'y' values
295 \param htol, vtol
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!
314 */
315 std::vector<std::vector<double> > multi_roots(SBasis const &f,
316 std::vector<double> const &levels,
317 double htol,
318 double vtol,
319 double a,
320 double b){
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));
327 return(roots);
328 }
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);
341 return;
342 }
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);
346 }
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];
353 if(d != 0) {
354 double r = s[0][0] / d;
355 if(0 <= r && r <= 1)
356 res.push_back(r);
357 }
358 return res;
359 }
361 /** Find all t s.t s(t) = 0
362 \param a sbasis function
363 \returns vector of zeros (roots)
365 */
366 std::vector<double> roots(SBasis const & s) {
367 switch(s.size()) {
368 case 0:
369 return std::vector<double>();
370 case 1:
371 return roots1(s);
372 default:
373 {
374 Bezier bz;
375 sbasis_to_bezier(bz, s);
376 return bz.roots();
377 }
378 }
379 }
381 };
383 /*
384 Local Variables:
385 mode:c++
386 c-file-style:"stroustrup"
387 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
388 indent-tabs-mode:nil
389 fill-column:99
390 End:
391 */
392 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :