src/2geom/sbasis-roots.cpp

   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.
  13
  14     in essence we are implementing quicksort on a continuous function
  15
  16  *  the gsl poly roots finder is faster than bernstein too, but we don't use it for 3 reasons:
  17
  18  a) it requires convertion to poly, which is numerically unstable
  19
  20  b) it requires gsl (which is currently not a dependency, and would bring in a whole slew of unrelated stuff)
  21
  22  c) it finds all roots, even complex ones.  We don't want to accidently treat a nearly real root as a real root
  23
  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
  27
  28 jfbarraud: eigenvalue root finding could be done directly in sbasis space ?
  29
  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.
  34
  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).
  38
  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.
  41
  42  **/
  43
  44 #include <cmath>
  45 #include <map>
  46
  47 #include <2geom/sbasis.h>
  48 #include <2geom/sbasis-to-bezier.h>
  49 #include <2geom/solver.h>
  50
  51 using namespace std;
  52
  53 namespace Geom{
  54
  55 /** Find the smallest interval that bounds a
  56  \param a sbasis function
  57  \returns inteval
  58
  59 */
  60
  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
  82
  83 /** Find a small interval that bounds a
  84  \param a sbasis function
  85  \returns inteval
  86
  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.
  95
  96     for(int j = sb.size()-1; j>=order; j--) {
  97         double a=sb[j][0];
  98         double b=sb[j][1];
  99
 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         }
 108
 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 }
 120
 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
 126
 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];
 137
 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         }
 145
 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 }
 157
 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!
 172
 173
 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 }
 177
 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){
 193
 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     }
 234
 235     int idxa=upper_level(levels,fa,vtol);
 236     int idxb=upper_level(levels,fb,vtol);
 237
 238     Interval bs = *bounds_local(df,Interval(a,b));
 239
 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.
 244
 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     }
 265
 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.
 271
 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 }
 291
 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
 298
 299 Effectively computes:
 300 results = roots(f(y_i)) for all y_i
 301
 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...
 312
 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){
 321
 322     std::vector<std::vector<double> > roots(levels.size(), std::vector<double>());
 323
 324     SBasis df=derivative(f);
 325     multi_roots_internal(f,df,levels,roots,htol,vtol,a,f(a),b,f(b));
 326
 327     return(roots);
 328 }
 329 //-------------------------------------
 330
 331
 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 }
 347
 348 // It is faster to use the bernstein root finder for small degree polynomials (<100?.
 349
 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 }
 360
 361 /** Find all t s.t s(t) = 0
 362  \param a sbasis function
 363  \returns vector of zeros (roots)
 364
 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 }
 380
 381 };
 382
 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 :