src/2geom/sbasis.cpp

   1 /*
   2  *  sbasis.cpp - S-power basis function class + supporting classes
   3  *
   4  *  Authors:
   5  *   Nathan Hurst <njh@mail.csse.monash.edu.au>
   6  *   Michael Sloan <mgsloan@gmail.com>
   7  *
   8  * Copyright (C) 2006-2007 authors
   9  *
  10  * This library is free software; you can redistribute it and/or
  11  * modify it either under the terms of the GNU Lesser General Public
  12  * License version 2.1 as published by the Free Software Foundation
  13  * (the "LGPL") or, at your option, under the terms of the Mozilla
  14  * Public License Version 1.1 (the "MPL"). If you do not alter this
  15  * notice, a recipient may use your version of this file under either
  16  * the MPL or the LGPL.
  17  *
  18  * You should have received a copy of the LGPL along with this library
  19  * in the file COPYING-LGPL-2.1; if not, write to the Free Software
  20  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
  21  * You should have received a copy of the MPL along with this library
  22  * in the file COPYING-MPL-1.1
  23  *
  24  * The contents of this file are subject to the Mozilla Public License
  25  * Version 1.1 (the "License"); you may not use this file except in
  26  * compliance with the License. You may obtain a copy of the License at
  27  * http://www.mozilla.org/MPL/
  28  *
  29  * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
  30  * OF ANY KIND, either express or implied. See the LGPL or the MPL for
  31  * the specific language governing rights and limitations.
  32  */
  33
  34 #include <cmath>
  35
  36 #include <2geom/sbasis.h>
  37 #include <2geom/isnan.h>
  38
  39 namespace Geom{
  40
  41 /*** At some point we should work on tighter bounds for the error.  It is clear that the error is
  42  * bounded by the L1 norm over the tail of the series, but this is very loose, leading to far too
  43  * many cubic beziers.  I've changed this to be \sum _i=tail ^\infty |hat a_i| 2^-i but I have no
  44  * evidence that this is correct.
  45  */
  46
  47 /*
  48 double SBasis::tail_error(unsigned tail) const {
  49     double err = 0, s = 1./(1<<(2*tail)); // rough
  50     for(unsigned i = tail; i < size(); i++) {
  51         err += (fabs((*this)[i][0]) + fabs((*this)[i][1]))*s;
  52         s /= 4;
  53     }
  54     return err;
  55 }
  56 */
  57
  58 double SBasis::tailError(unsigned tail) const {
  59   Interval bs = bounds_fast(*this, tail);
  60   return std::max(fabs(bs.min()),fabs(bs.max()));
  61 }
  62
  63 bool SBasis::isFinite() const {
  64     for(unsigned i = 0; i < size(); i++) {
  65         if(!(*this)[i].isFinite())
  66             return false;
  67     }
  68     return true;
  69 }
  70
  71 std::vector<double> SBasis::valueAndDerivatives(double t, unsigned n) const {
  72     std::vector<double> ret(n+1);
  73     ret.push_back(valueAt(t));
  74     SBasis tmp = *this;
  75     for(unsigned i = 0; i < n; i++) {
  76         tmp.derive();
  77         ret[i] = tmp.valueAt(t);
  78     }
  79     return ret;
  80 }
  81
  82
  83 SBasis operator+(const SBasis& a, const SBasis& b) {
  84     SBasis result;
  85     const unsigned out_size = std::max(a.size(), b.size());
  86     const unsigned min_size = std::min(a.size(), b.size());
  87     result.reserve(out_size);
  88
  89     for(unsigned i = 0; i < min_size; i++) {
  90         result.push_back(a[i] + b[i]);
  91     }
  92     for(unsigned i = min_size; i < a.size(); i++)
  93         result.push_back(a[i]);
  94     for(unsigned i = min_size; i < b.size(); i++)
  95         result.push_back(b[i]);
  96
  97     assert(result.size() == out_size);
  98     return result;
  99 }
 100
 101 SBasis operator-(const SBasis& a, const SBasis& b) {
 102     SBasis result;
 103     const unsigned out_size = std::max(a.size(), b.size());
 104     const unsigned min_size = std::min(a.size(), b.size());
 105     result.reserve(out_size);
 106
 107     for(unsigned i = 0; i < min_size; i++) {
 108         result.push_back(a[i] - b[i]);
 109     }
 110     for(unsigned i = min_size; i < a.size(); i++)
 111         result.push_back(a[i]);
 112     for(unsigned i = min_size; i < b.size(); i++)
 113         result.push_back(-b[i]);
 114
 115     assert(result.size() == out_size);
 116     return result;
 117 }
 118
 119 SBasis& operator+=(SBasis& a, const SBasis& b) {
 120     const unsigned out_size = std::max(a.size(), b.size());
 121     const unsigned min_size = std::min(a.size(), b.size());
 122     a.reserve(out_size);
 123
 124     for(unsigned i = 0; i < min_size; i++)
 125         a[i] += b[i];
 126     for(unsigned i = min_size; i < b.size(); i++)
 127         a.push_back(b[i]);
 128
 129     assert(a.size() == out_size);
 130     return a;
 131 }
 132
 133 SBasis& operator-=(SBasis& a, const SBasis& b) {
 134     const unsigned out_size = std::max(a.size(), b.size());
 135     const unsigned min_size = std::min(a.size(), b.size());
 136     a.reserve(out_size);
 137
 138     for(unsigned i = 0; i < min_size; i++)
 139         a[i] -= b[i];
 140     for(unsigned i = min_size; i < b.size(); i++)
 141         a.push_back(-b[i]);
 142
 143     assert(a.size() == out_size);
 144     return a;
 145 }
 146
 147 SBasis operator*(SBasis const &a, double k) {
 148     SBasis c;
 149     c.reserve(a.size());
 150     for(unsigned i = 0; i < a.size(); i++)
 151         c.push_back(a[i] * k);
 152     return c;
 153 }
 154
 155 SBasis& operator*=(SBasis& a, double b) {
 156     if (a.isZero()) return a;
 157     if (b == 0)
 158         a.clear();
 159     else
 160         for(unsigned i = 0; i < a.size(); i++)
 161             a[i] *= b;
 162     return a;
 163 }
 164
 165 SBasis shift(SBasis const &a, int sh) {
 166     SBasis c = a;
 167     if(sh > 0) {
 168         c.insert(c.begin(), sh, Linear(0,0));
 169     } else {
 170         //TODO: truncate
 171     }
 172     return c;
 173 }
 174
 175 SBasis shift(Linear const &a, int sh) {
 176     SBasis c;
 177     if(sh > 0) {
 178         c.insert(c.begin(), sh, Linear(0,0));
 179         c.push_back(a);
 180     }
 181     return c;
 182 }
 183
 184 #if 0
 185 SBasis multiply(SBasis const &a, SBasis const &b) {
 186     // c = {a0*b0 - shift(1, a.Tri*b.Tri), a1*b1 - shift(1, a.Tri*b.Tri)}
 187
 188     // shift(1, a.Tri*b.Tri)
 189     SBasis c;
 190     if(a.isZero() || b.isZero())
 191         return c;
 192     c.resize(a.size() + b.size(), Linear(0,0));
 193     for(unsigned j = 0; j < b.size(); j++) {
 194         for(unsigned i = j; i < a.size()+j; i++) {
 195             double tri = Tri(b[j])*Tri(a[i-j]);
 196             c[i+1/*shift*/] += Linear(Hat(-tri));
 197         }
 198     }
 199     for(unsigned j = 0; j < b.size(); j++) {
 200         for(unsigned i = j; i < a.size()+j; i++) {
 201             for(unsigned dim = 0; dim < 2; dim++)
 202                 c[i][dim] += b[j][dim]*a[i-j][dim];
 203         }
 204     }
 205     c.normalize();
 206     //assert(!(0 == c.back()[0] && 0 == c.back()[1]));
 207     return c;
 208 }
 209 #else
 210
 211 SBasis multiply_add(SBasis const &a, SBasis const &b, SBasis c) {
 212     if(a.isZero() || b.isZero())
 213         return c;
 214     c.resize(a.size() + b.size(), Linear(0,0));
 215     for(unsigned j = 0; j < b.size(); j++) {
 216         for(unsigned i = j; i < a.size()+j; i++) {
 217             double tri = Tri(b[j])*Tri(a[i-j]);
 218             c[i+1/*shift*/] += Linear(Hat(-tri));
 219         }
 220     }
 221     for(unsigned j = 0; j < b.size(); j++) {
 222         for(unsigned i = j; i < a.size()+j; i++) {
 223             for(unsigned dim = 0; dim < 2; dim++)
 224                 c[i][dim] += b[j][dim]*a[i-j][dim];
 225         }
 226     }
 227     c.normalize();
 228     //assert(!(0 == c.back()[0] && 0 == c.back()[1]));
 229     return c;
 230 }
 231
 232 SBasis multiply(SBasis const &a, SBasis const &b) {
 233     SBasis c;
 234     if(a.isZero() || b.isZero())
 235         return c;
 236     return multiply_add(a, b, c);
 237 }
 238 #endif
 239 SBasis integral(SBasis const &c) {
 240     SBasis a;
 241     a.resize(c.size() + 1, Linear(0,0));
 242     a[0] = Linear(0,0);
 243
 244     for(unsigned k = 1; k < c.size() + 1; k++) {
 245         double ahat = -Tri(c[k-1])/(2*k);
 246         a[k] = Hat(ahat);
 247     }
 248     double aTri = 0;
 249     for(int k = c.size()-1; k >= 0; k--) {
 250         aTri = (Hat(c[k]).d + (k+1)*aTri/2)/(2*k+1);
 251         a[k][0] -= aTri/2;
 252         a[k][1] += aTri/2;
 253     }
 254     a.normalize();
 255     return a;
 256 }
 257
 258 SBasis derivative(SBasis const &a) {
 259     SBasis c;
 260     c.resize(a.size(), Linear(0,0));
 261     if(a.isZero())
 262         return c;
 263
 264     for(unsigned k = 0; k < a.size()-1; k++) {
 265         double d = (2*k+1)*(a[k][1] - a[k][0]);
 266
 267         c[k][0] = d + (k+1)*a[k+1][0];
 268         c[k][1] = d - (k+1)*a[k+1][1];
 269     }
 270     int k = a.size()-1;
 271     double d = (2*k+1)*(a[k][1] - a[k][0]);
 272     if(d == 0)
 273         c.pop_back();
 274     else {
 275         c[k][0] = d;
 276         c[k][1] = d;
 277     }
 278
 279     return c;
 280 }
 281
 282 void SBasis::derive() { // in place version
 283     if(isZero()) return;
 284     for(unsigned k = 0; k < size()-1; k++) {
 285         double d = (2*k+1)*((*this)[k][1] - (*this)[k][0]);
 286
 287         (*this)[k][0] = d + (k+1)*(*this)[k+1][0];
 288         (*this)[k][1] = d - (k+1)*(*this)[k+1][1];
 289     }
 290     int k = size()-1;
 291     double d = (2*k+1)*((*this)[k][1] - (*this)[k][0]);
 292     if(d == 0)
 293         pop_back();
 294     else {
 295         (*this)[k][0] = d;
 296         (*this)[k][1] = d;
 297     }
 298 }
 299
 300 //TODO: convert int k to unsigned k, and remove cast
 301 SBasis sqrt(SBasis const &a, int k) {
 302     SBasis c;
 303     if(a.isZero() || k == 0)
 304         return c;
 305     c.resize(k, Linear(0,0));
 306     c[0] = Linear(std::sqrt(a[0][0]), std::sqrt(a[0][1]));
 307     SBasis r = a - multiply(c, c); // remainder
 308
 309     for(unsigned i = 1; i <= (unsigned)k and i<r.size(); i++) {
 310         Linear ci(r[i][0]/(2*c[0][0]), r[i][1]/(2*c[0][1]));
 311         SBasis cisi = shift(ci, i);
 312         r -= multiply(shift((c*2 + cisi), i), SBasis(ci));
 313         r.truncate(k+1);
 314         c += cisi;
 315         if(r.tailError(i) == 0) // if exact
 316             break;
 317     }
 318
 319     return c;
 320 }
 321
 322 // return a kth order approx to 1/a)
 323 SBasis reciprocal(Linear const &a, int k) {
 324     SBasis c;
 325     assert(!a.isZero());
 326     c.resize(k, Linear(0,0));
 327     double r_s0 = (Tri(a)*Tri(a))/(-a[0]*a[1]);
 328     double r_s0k = 1;
 329     for(unsigned i = 0; i < (unsigned)k; i++) {
 330         c[i] = Linear(r_s0k/a[0], r_s0k/a[1]);
 331         r_s0k *= r_s0;
 332     }
 333     return c;
 334 }
 335
 336 SBasis divide(SBasis const &a, SBasis const &b, int k) {
 337     SBasis c;
 338     assert(!a.isZero());
 339     SBasis r = a; // remainder
 340
 341     k++;
 342     r.resize(k, Linear(0,0));
 343     c.resize(k, Linear(0,0));
 344
 345     for(unsigned i = 0; i < (unsigned)k; i++) {
 346         Linear ci(r[i][0]/b[0][0], r[i][1]/b[0][1]); //H0
 347         c[i] += ci;
 348         r -= shift(multiply(ci,b), i);
 349         r.truncate(k+1);
 350         if(r.tailError(i) == 0) // if exact
 351             break;
 352     }
 353
 354     return c;
 355 }
 356
 357 // a(b)
 358 // return a0 + s(a1 + s(a2 +...  where s = (1-u)u; ak =(1 - u)a^0_k + ua^1_k
 359 SBasis compose(SBasis const &a, SBasis const &b) {
 360     SBasis s = multiply((SBasis(Linear(1,1))-b), b);
 361     SBasis r;
 362
 363     for(int i = a.size()-1; i >= 0; i--) {
 364         r = multiply_add(r, s, SBasis(Linear(Hat(a[i][0]))) - b*a[i][0] + b*a[i][1]);
 365     }
 366     return r;
 367 }
 368
 369 // a(b)
 370 // return a0 + s(a1 + s(a2 +...  where s = (1-u)u; ak =(1 - u)a^0_k + ua^1_k
 371 SBasis compose(SBasis const &a, SBasis const &b, unsigned k) {
 372     SBasis s = multiply((SBasis(Linear(1,1))-b), b);
 373     SBasis r;
 374
 375     for(int i = a.size()-1; i >= 0; i--) {
 376         r = multiply_add(r, s, SBasis(Linear(Hat(a[i][0]))) - b*a[i][0] + b*a[i][1]);
 377     }
 378     r.truncate(k);
 379     return r;
 380 }
 381
 382 /*
 383 Inversion algorithm. The notation is certainly very misleading. The
 384 pseudocode should say:
 385
 386 c(v) := 0
 387 r(u) := r_0(u) := u
 388 for i:=0 to k do
 389   c_i(v) := H_0(r_i(u)/(t_1)^i; u)
 390   c(v) := c(v) + c_i(v)*t^i
 391   r(u) := r(u) ? c_i(u)*(t(u))^i
 392 endfor
 393 */
 394
 395 //#define DEBUG_INVERSION 1
 396
 397 SBasis inverse(SBasis a, int k) {
 398     assert(a.size() > 0);
 399 // the function should have 'unit range'("a00 = 0 and a01 = 1") and be monotonic.
 400     double a0 = a[0][0];
 401     if(a0 != 0) {
 402         a -= a0;
 403     }
 404     double a1 = a[0][1];
 405     assert(a1 != 0); // not invertable.
 406
 407     if(a1 != 1) {
 408         a /= a1;
 409     }
 410     SBasis c;                           // c(v) := 0
 411     if(a.size() >= 2 && k == 2) {
 412         c.push_back(Linear(0,1));
 413         Linear t1(1+a[1][0], 1-a[1][1]);    // t_1
 414         c.push_back(Linear(-a[1][0]/t1[0], -a[1][1]/t1[1]));
 415     } else if(a.size() >= 2) {                      // non linear
 416         SBasis r = Linear(0,1);             // r(u) := r_0(u) := u
 417         Linear t1(1./(1+a[1][0]), 1./(1-a[1][1]));    // 1./t_1
 418         Linear one(1,1);
 419         Linear t1i = one;                   // t_1^0
 420         SBasis one_minus_a = SBasis(one) - a;
 421         SBasis t = multiply(one_minus_a, a); // t(u)
 422         SBasis ti(one);                     // t(u)^0
 423 #ifdef DEBUG_INVERSION
 424         std::cout << "a=" << a << std::endl;
 425         std::cout << "1-a=" << one_minus_a << std::endl;
 426         std::cout << "t1=" << t1 << std::endl;
 427         //assert(t1 == t[1]);
 428 #endif
 429
 430         c.resize(k+1, Linear(0,0));
 431         for(unsigned i = 0; i < (unsigned)k; i++) {   // for i:=0 to k do
 432 #ifdef DEBUG_INVERSION
 433             std::cout << "-------" << i << ": ---------" <<std::endl;
 434             std::cout << "r=" << r << std::endl
 435                       << "c=" << c << std::endl
 436                       << "ti=" << ti << std::endl
 437                       << std::endl;
 438 #endif
 439             if(r.size() <= i)                // ensure enough space in the remainder, probably not needed
 440                 r.resize(i+1, Linear(0,0));
 441             Linear ci(r[i][0]*t1i[0], r[i][1]*t1i[1]); // c_i(v) := H_0(r_i(u)/(t_1)^i; u)
 442 #ifdef DEBUG_INVERSION
 443             std::cout << "t1i=" << t1i << std::endl;
 444             std::cout << "ci=" << ci << std::endl;
 445 #endif
 446             for(int dim = 0; dim < 2; dim++) // t1^-i *= 1./t1
 447                 t1i[dim] *= t1[dim];
 448             c[i] = ci; // c(v) := c(v) + c_i(v)*t^i
 449             // change from v to u parameterisation
 450             SBasis civ = one_minus_a*ci[0] + a*ci[1];
 451             // r(u) := r(u) - c_i(u)*(t(u))^i
 452             // We can truncate this to the number of final terms, as no following terms can
 453             // contribute to the result.
 454             r -= multiply(civ,ti);
 455             r.truncate(k);
 456             if(r.tailError(i) == 0)
 457                 break; // yay!
 458             ti = multiply(ti,t);
 459         }
 460 #ifdef DEBUG_INVERSION
 461         std::cout << "##########################" << std::endl;
 462 #endif
 463     } else
 464         c = Linear(0,1); // linear
 465     c -= a0; // invert the offset
 466     c /= a1; // invert the slope
 467     return c;
 468 }
 469
 470 SBasis sin(Linear b, int k) {
 471     SBasis s = Linear(std::sin(b[0]), std::sin(b[1]));
 472     Tri tr(s[0]);
 473     double t2 = Tri(b);
 474     s.push_back(Linear(std::cos(b[0])*t2 - tr, -std::cos(b[1])*t2 + tr));
 475
 476     t2 *= t2;
 477     for(int i = 0; i < k; i++) {
 478         Linear bo(4*(i+1)*s[i+1][0] - 2*s[i+1][1],
 479                   -2*s[i+1][0] + 4*(i+1)*s[i+1][1]);
 480         bo -= s[i]*(t2/(i+1));
 481
 482
 483         s.push_back(bo/double(i+2));
 484     }
 485
 486     return s;
 487 }
 488
 489 SBasis cos(Linear bo, int k) {
 490     return sin(Linear(bo[0] + M_PI/2,
 491                       bo[1] + M_PI/2),
 492                k);
 493 }
 494
 495 //compute fog^-1. ("zero" = double comparison threshold. *!*we might divide by "zero"*!*)
 496 //TODO: compute order according to tol?
 497 //TODO: requires g(0)=0 & g(1)=1 atm... adaptation to other cases should be obvious!
 498 SBasis compose_inverse(SBasis const &f, SBasis const &g, unsigned order, double zero){
 499     SBasis result; //result
 500     SBasis r=f; //remainder
 501     SBasis Pk=Linear(1)-g,Qk=g,sg=Pk*Qk;
 502     Pk.truncate(order);
 503     Qk.truncate(order);
 504     Pk.resize(order,Linear(0.));
 505     Qk.resize(order,Linear(0.));
 506     r.resize(order,Linear(0.));
 507
 508     int vs= valuation(sg,zero);
 509
 510     for (unsigned k=0; k<order; k+=vs){
 511         double p10 = Pk.at(k)[0];// we have to solve the linear system:
 512         double p01 = Pk.at(k)[1];//
 513         double q10 = Qk.at(k)[0];//   p10*a + q10*b = r10
 514         double q01 = Qk.at(k)[1];// &
 515         double r10 =  r.at(k)[0];//   p01*a + q01*b = r01
 516         double r01 =  r.at(k)[1];//
 517         double a,b;
 518         double det = p10*q01-p01*q10;
 519
 520         //TODO: handle det~0!!
 521         if (fabs(det)<zero){
 522             det = zero;
 523             a=b=0;
 524         }else{
 525             a=( q01*r10-q10*r01)/det;
 526             b=(-p01*r10+p10*r01)/det;
 527         }
 528         result.push_back(Linear(a,b));
 529         r=r-Pk*a-Qk*b;
 530
 531         Pk=Pk*sg;
 532         Qk=Qk*sg;
 533         Pk.truncate(order);
 534         Qk.truncate(order);
 535         r.truncate(order);
 536     }
 537     result.normalize();
 538     return result;
 539 }
 540
 541 }
 542
 543 /*
 544   Local Variables:
 545   mode:c++
 546   c-file-style:"stroustrup"
 547   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 548   indent-tabs-mode:nil
 549   fill-column:99
 550   End:
 551 */
 552 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :