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 "sbasis.h"
  48 #include "sbasis-to-bezier.h"
  49 #include "solver.h"
  50
  51 using namespace std;
  52
  53 namespace Geom{
  54
  55 Interval bounds_exact(SBasis const &a) {
  56     Interval result = Interval(a.at0(), a.at1());
  57     SBasis df = derivative(a);
  58     vector<double>extrema = roots(df);
  59     for (unsigned i=0; i<extrema.size(); i++){
  60         result.extendTo(a(extrema[i]));
  61     }
  62     return result;
  63 }
  64
  65 Interval bounds_fast(const SBasis &sb, int order) {
  66     Interval res;
  67     for(int j = sb.size()-1; j>=order; j--) {
  68         double a=sb[j][0];
  69         double b=sb[j][1];
  70
  71         double v, t = 0;
  72         v = res[0];
  73         if (v<0) t = ((b-a)/v+1)*0.5;
  74         if (v>=0 || t<0 || t>1) {
  75             res[0] = std::min(a,b);
  76         }else{
  77             res[0]=lerp(t, a+v*t, b);
  78         }
  79
  80         v = res[1];
  81         if (v>0) t = ((b-a)/v+1)*0.5;
  82         if (v<=0 || t<0 || t>1) {
  83             res[1] = std::max(a,b);
  84         }else{
  85             res[1]=lerp(t, a+v*t, b);
  86         }
  87     }
  88     if (order>0) res*=pow(.25,order);
  89     return res;
  90 }
  91
  92 Interval bounds_local(const SBasis &sb, const Interval &i, int order) {
  93     double t0=i.min(), t1=i.max(), lo=0., hi=0.;
  94     for(int j = sb.size()-1; j>=order; j--) {
  95         double a=sb[j][0];
  96         double b=sb[j][1];
  97
  98         double t = 0;
  99         if (lo<0) t = ((b-a)/lo+1)*0.5;
 100         if (lo>=0 || t<t0 || t>t1) {
 101             lo = std::min(a*(1-t0)+b*t0+lo*t0*(1-t0),a*(1-t1)+b*t1+lo*t1*(1-t1));
 102         }else{
 103             lo = lerp(t, a+lo*t, b);
 104         }
 105
 106         if (hi>0) t = ((b-a)/hi+1)*0.5;
 107         if (hi<=0 || t<t0 || t>t1) {
 108             hi = std::max(a*(1-t0)+b*t0+hi*t0*(1-t0),a*(1-t1)+b*t1+hi*t1*(1-t1));
 109         }else{
 110             hi = lerp(t, a+hi*t, b);
 111         }
 112     }
 113     Interval res = Interval(lo,hi);
 114     if (order>0) res*=pow(.25,order);
 115     return res;
 116 }
 117
 118 //-- multi_roots ------------------------------------
 119 // goal: solve f(t)=c for several c at once.
 120 /* algo: -compute f at both ends of the given segment [a,b].
 121          -compute bounds m<df(t)<M for df on the segment.
 122          let c and C be the levels below and above f(a):
 123          going from f(a) down to c with slope m takes at least time (f(a)-c)/m
 124          going from f(a)  up  to C with slope M takes at least time (C-f(a))/M
 125          From this we conclude there are no roots before a'=a+min((f(a)-c)/m,(C-f(a))/M).
 126          Do the same for b: compute some b' such that there are no roots in (b',b].
 127          -if [a',b'] is not empty, repeat the process with [a',(a'+b')/2] and [(a'+b')/2,b'].
 128   unfortunately, extra care is needed about rounding errors, and also to avoid the repetition of roots,
 129   making things tricky and unpleasant...
 130 */
 131 //TODO: Make sure the code is "rounding-errors proof" and take care about repetition of roots!
 132
 133
 134 static int upper_level(vector<double> const &levels,double x,double tol=0.){
 135     return(upper_bound(levels.begin(),levels.end(),x-tol)-levels.begin());
 136 }
 137
 138 static void multi_roots_internal(SBasis const &f,
 139                                  SBasis const &df,
 140                                  std::vector<double> const &levels,
 141                                  std::vector<std::vector<double> > &roots,
 142                                  double htol,
 143                                  double vtol,
 144                                  double a,
 145                                  double fa,
 146                                  double b,
 147                                  double fb){
 148
 149     if (f.size()==0){
 150         int idx;
 151         idx=upper_level(levels,0,vtol);
 152         if (idx<(int)levels.size()&&fabs(levels.at(idx))<=vtol){
 153             roots[idx].push_back(a);
 154             roots[idx].push_back(b);
 155         }
 156         return;
 157     }
 158 ////usefull?
 159 //     if (f.size()==1){
 160 //         int idxa=upper_level(levels,fa);
 161 //         int idxb=upper_level(levels,fb);
 162 //         if (fa==fb){
 163 //             if (fa==levels[idxa]){
 164 //                 roots[a]=idxa;
 165 //                 roots[b]=idxa;
 166 //             }
 167 //             return;
 168 //         }
 169 //         int idx_min=std::min(idxa,idxb);
 170 //         int idx_max=std::max(idxa,idxb);
 171 //         if (idx_max==levels.size()) idx_max-=1;
 172 //         for(int i=idx_min;i<=idx_max; i++){
 173 //             double t=a+(b-a)*(levels[i]-fa)/(fb-fa);
 174 //             if(a<t&&t<b) roots[t]=i;
 175 //         }
 176 //         return;
 177 //     }
 178     if ((b-a)<htol){
 179         //TODO: use different tol for t and f ?
 180         //TODO: unsigned idx ? (remove int casts when fixed)
 181         int idx=std::min(upper_level(levels,fa,vtol),upper_level(levels,fb,vtol));
 182         if (idx==(int)levels.size()) idx-=1;
 183         double c=levels.at(idx);
 184         if((fa-c)*(fb-c)<=0||fabs(fa-c)<vtol||fabs(fb-c)<vtol){
 185             roots[idx].push_back((a+b)/2);
 186         }
 187         return;
 188     }
 189
 190     int idxa=upper_level(levels,fa,vtol);
 191     int idxb=upper_level(levels,fb,vtol);
 192
 193     Interval bs = bounds_local(df,Interval(a,b));
 194
 195     //first times when a level (higher or lower) can be reached from a or b.
 196     double ta_hi,tb_hi,ta_lo,tb_lo;
 197     ta_hi=ta_lo=b+1;//default values => no root there.
 198     tb_hi=tb_lo=a-1;//default values => no root there.
 199
 200     if (idxa<(int)levels.size() && fabs(fa-levels.at(idxa))<vtol){//a can be considered a root.
 201         //ta_hi=ta_lo=a;
 202         roots[idxa].push_back(a);
 203         ta_hi=ta_lo=a+htol;
 204     }else{
 205         if (bs.max()>0 && idxa<(int)levels.size())
 206             ta_hi=a+(levels.at(idxa  )-fa)/bs.max();
 207         if (bs.min()<0 && idxa>0)
 208             ta_lo=a+(levels.at(idxa-1)-fa)/bs.min();
 209     }
 210     if (idxb<(int)levels.size() && fabs(fb-levels.at(idxb))<vtol){//b can be considered a root.
 211         //tb_hi=tb_lo=b;
 212         roots[idxb].push_back(b);
 213         tb_hi=tb_lo=b-htol;
 214     }else{
 215         if (bs.min()<0 && idxb<(int)levels.size())
 216             tb_hi=b+(levels.at(idxb  )-fb)/bs.min();
 217         if (bs.max()>0 && idxb>0)
 218             tb_lo=b+(levels.at(idxb-1)-fb)/bs.max();
 219     }
 220
 221     double t0,t1;
 222     t0=std::min(ta_hi,ta_lo);
 223     t1=std::max(tb_hi,tb_lo);
 224     //hum, rounding errors frighten me! so I add this +tol...
 225     if (t0>t1+htol) return;//no root here.
 226
 227     if (fabs(t1-t0)<htol){
 228         multi_roots_internal(f,df,levels,roots,htol,vtol,t0,f(t0),t1,f(t1));
 229     }else{
 230         double t,t_left,t_right,ft,ft_left,ft_right;
 231         t_left =t_right =t =(t0+t1)/2;
 232         ft_left=ft_right=ft=f(t);
 233         int idx=upper_level(levels,ft,vtol);
 234         if (idx<(int)levels.size() && fabs(ft-levels.at(idx))<vtol){//t can be considered a root.
 235             roots[idx].push_back(t);
 236             //we do not want to count it twice (from the left and from the right)
 237             t_left =t-htol/2;
 238             t_right=t+htol/2;
 239             ft_left =f(t_left);
 240             ft_right=f(t_right);
 241         }
 242         multi_roots_internal(f,df,levels,roots,htol,vtol,t0     ,f(t0)   ,t_left,ft_left);
 243         multi_roots_internal(f,df,levels,roots,htol,vtol,t_right,ft_right,t1    ,f(t1)  );
 244     }
 245 }
 246
 247 std::vector<std::vector<double> > multi_roots(SBasis const &f,
 248                                       std::vector<double> const &levels,
 249                                       double htol,
 250                                       double vtol,
 251                                       double a,
 252                                       double b){
 253
 254     std::vector<std::vector<double> > roots(levels.size(), std::vector<double>());
 255
 256     SBasis df=derivative(f);
 257     multi_roots_internal(f,df,levels,roots,htol,vtol,a,f(a),b,f(b));
 258
 259     return(roots);
 260 }
 261 //-------------------------------------
 262
 263 #if 0
 264 double Laguerre_internal(SBasis const & p,
 265                          double x0,
 266                          double tol,
 267                          bool & quad_root) {
 268     double a = 2*tol;
 269     double xk = x0;
 270     double n = p.size();
 271     quad_root = false;
 272     while(a > tol) {
 273         //std::cout << "xk = " << xk << std::endl;
 274         Linear b = p.back();
 275         Linear d(0), f(0);
 276         double err = fabs(b);
 277         double abx = fabs(xk);
 278         for(int j = p.size()-2; j >= 0; j--) {
 279             f = xk*f + d;
 280             d = xk*d + b;
 281             b = xk*b + p[j];
 282             err = fabs(b) + abx*err;
 283         }
 284
 285         err *= 1e-7; // magic epsilon for convergence, should be computed from tol
 286
 287         double px = b;
 288         if(fabs(b) < err)
 289             return xk;
 290         //if(std::norm(px) < tol*tol)
 291         //    return xk;
 292         double G = d / px;
 293         double H = G*G - f / px;
 294
 295         //std::cout << "G = " << G << "H = " << H;
 296         double radicand = (n - 1)*(n*H-G*G);
 297         //assert(radicand.real() > 0);
 298         if(radicand < 0)
 299             quad_root = true;
 300         //std::cout << "radicand = " << radicand << std::endl;
 301         if(G.real() < 0) // here we try to maximise the denominator avoiding cancellation
 302             a = - std::sqrt(radicand);
 303         else
 304             a = std::sqrt(radicand);
 305         //std::cout << "a = " << a << std::endl;
 306         a = n / (a + G);
 307         //std::cout << "a = " << a << std::endl;
 308         xk -= a;
 309     }
 310     //std::cout << "xk = " << xk << std::endl;
 311     return xk;
 312 }
 313 #endif
 314
 315 void subdiv_sbasis(SBasis const & s,
 316                    std::vector<double> & roots,
 317                    double left, double right) {
 318     Interval bs = bounds_fast(s);
 319     if(bs.min() > 0 || bs.max() < 0)
 320         return; // no roots here
 321     if(s.tailError(1) < 1e-7) {
 322         double t = s[0][0] / (s[0][0] - s[0][1]);
 323         roots.push_back(left*(1-t) + t*right);
 324         return;
 325     }
 326     double middle = (left + right)/2;
 327     subdiv_sbasis(compose(s, Linear(0, 0.5)), roots, left, middle);
 328     subdiv_sbasis(compose(s, Linear(0.5, 1.)), roots, middle, right);
 329 }
 330
 331 // It is faster to use the bernstein root finder for small degree polynomials (<100?.
 332
 333 std::vector<double> roots1(SBasis const & s) {
 334     std::vector<double> res;
 335     double d = s[0][0] - s[0][1];
 336     if(d != 0) {
 337         double r = s[0][0] / d;
 338         if(0 <= r and r <= 1)
 339             res.push_back(r);
 340     }
 341     return res;
 342 }
 343
 344 std::vector<double> roots(SBasis const & s) {
 345     switch(s.size()) {
 346     case 0:
 347         return std::vector<double>();
 348     case 1:
 349         return roots1(s);
 350     default:
 351         return sbasis_to_bezier(s).roots();
 352     }
 353 }
 354
 355 };
 356
 357 /*
 358   Local Variables:
 359   mode:c++
 360   c-file-style:"stroustrup"
 361   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 362   indent-tabs-mode:nil
 363   fill-column:99
 364   End:
 365 */
 366 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :