trunk/src/2geom/sbasis.h

   1 /**
   2  * \file
   3  * \brief Defines S-power basis function class
   4  *
   5  *  Authors:
   6  *   Nathan Hurst <njh@mail.csse.monash.edu.au>
   7  *   Michael Sloan <mgsloan@gmail.com>
   8  *
   9  * Copyright (C) 2006-2007 authors
  10  *
  11  * This library is free software; you can redistribute it and/or
  12  * modify it either under the terms of the GNU Lesser General Public
  13  * License version 2.1 as published by the Free Software Foundation
  14  * (the "LGPL") or, at your option, under the terms of the Mozilla
  15  * Public License Version 1.1 (the "MPL"). If you do not alter this
  16  * notice, a recipient may use your version of this file under either
  17  * the MPL or the LGPL.
  18  *
  19  * You should have received a copy of the LGPL along with this library
  20  * in the file COPYING-LGPL-2.1; if not, write to the Free Software
  21  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
  22  * You should have received a copy of the MPL along with this library
  23  * in the file COPYING-MPL-1.1
  24  *
  25  * The contents of this file are subject to the Mozilla Public License
  26  * Version 1.1 (the "License"); you may not use this file except in
  27  * compliance with the License. You may obtain a copy of the License at
  28  * http://www.mozilla.org/MPL/
  29  *
  30  * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
  31  * OF ANY KIND, either express or implied. See the LGPL or the MPL for
  32  * the specific language governing rights and limitations.
  33  */
  34
  35 #ifndef SEEN_SBASIS_H
  36 #define SEEN_SBASIS_H
  37 #include <vector>
  38 #include <cassert>
  39 #include <iostream>
  40
  41 #include <2geom/linear.h>
  42 #include <2geom/interval.h>
  43 #include <2geom/utils.h>
  44 #include <2geom/exception.h>
  45
  46 //#define USE_SBASISN 1
  47
  48
  49 #if defined(USE_SBASIS_OF)
  50
  51 #include "sbasis-of.h"
  52
  53 #elif defined(USE_SBASISN)
  54
  55 #include "sbasisN.h"
  56 namespace Geom{
  57
  58 /*** An empty SBasis is identically 0. */
  59 class SBasis : public SBasisN<1>;
  60
  61 };
  62 #else
  63
  64 namespace Geom{
  65
  66 /*** An empty SBasis is identically 0. */
  67 class SBasis{
  68     std::vector<Linear> d;
  69     void push_back(Linear const&l) { d.push_back(l); }
  70
  71 public:
  72     // As part of our migration away from SBasis isa vector we provide this minimal set of vector interface methods.
  73     size_t size() const {return d.size();}
  74     Linear operator[](unsigned i) const {
  75         return d[i];
  76     }
  77     Linear& operator[](unsigned i) { return d.at(i); }
  78     Linear const* begin() const { return (Linear const*)&*d.begin();}
  79     Linear const* end() const { return (Linear const*)&*d.end();}
  80     Linear* begin() { return (Linear*)&*d.begin();}
  81     Linear* end() { return (Linear*)&*d.end();}
  82     bool empty() const {return d.empty();}
  83     Linear &back() {return d.back();}
  84     Linear const &back() const {return d.back();}
  85     void pop_back() { d.pop_back();}
  86     void resize(unsigned n) { d.resize(n);}
  87     void resize(unsigned n, Linear const& l) { d.resize(n, l);}
  88     void reserve(unsigned n) { d.reserve(n);}
  89     void clear() {d.clear();}
  90     void insert(Linear* before, const Linear* src_begin, const Linear* src_end) { d.insert(std::vector<Linear>::iterator(before), src_begin, src_end);}
  91     //void insert(Linear* aa, Linear* bb, Linear* cc} { d.insert(aa, bb, cc);}
  92     Linear& at(unsigned i) { return d.at(i);}
  93     //void insert(Linear* before, int& n, Linear const &l) { d.insert(std::vector<Linear>::iterator(before), n, l);}
  94     bool operator==(SBasis const&B) { return d == B.d;}
  95     operator std::vector<Linear>() { return d;}
  96
  97
  98     SBasis() {}
  99     explicit SBasis(double a) {
 100         push_back(Linear(a,a));
 101     }
 102     SBasis(SBasis const & a) :
 103         d(a.d)
 104     {}
 105     SBasis(Linear const & bo) {
 106         push_back(bo);
 107     }
 108     SBasis(Linear* bo) {
 109         push_back(*bo);
 110     }
 111     explicit SBasis(size_t n, Linear const&l) : d(n, l) {}
 112
 113     //IMPL: FragmentConcept
 114     typedef double output_type;
 115     inline bool isZero() const {
 116         if(empty()) return true;
 117         for(unsigned i = 0; i < size(); i++) {
 118             if(!(*this)[i].isZero()) return false;
 119         }
 120         return true;
 121     }
 122     inline bool isConstant() const {
 123         if (empty()) return true;
 124         for (unsigned i = 0; i < size(); i++) {
 125             if(!(*this)[i].isConstant()) return false;
 126         }
 127         return true;
 128     }
 129
 130     bool isFinite() const;
 131     inline double at0() const {
 132         if(empty()) return 0; else return (*this)[0][0];
 133     }
 134     inline double at1() const{
 135         if(empty()) return 0; else return (*this)[0][1];
 136     }
 137
 138     int degreesOfFreedom() const { return size()*2;}
 139
 140     double valueAt(double t) const {
 141         double s = t*(1-t);
 142         double p0 = 0, p1 = 0;
 143         for(unsigned k = size(); k > 0; k--) {
 144             const Linear &lin = (*this)[k-1];
 145             p0 = p0*s + lin[0];
 146             p1 = p1*s + lin[1];
 147         }
 148         return (1-t)*p0 + t*p1;
 149     }
 150     //double valueAndDerivative(double t, double &der) const {
 151     //}
 152     double operator()(double t) const {
 153         return valueAt(t);
 154     }
 155
 156     std::vector<double> valueAndDerivatives(double t, unsigned n) const;
 157
 158     SBasis toSBasis() const { return SBasis(*this); }
 159
 160     double tailError(unsigned tail) const;
 161
 162 // compute f(g)
 163     SBasis operator()(SBasis const & g) const;
 164
 165 //MUTATOR PRISON
 166     //remove extra zeros
 167     void normalize() {
 168         while(!empty() && 0 == back()[0] && 0 == back()[1])
 169             pop_back();
 170     }
 171
 172     void truncate(unsigned k) { if(k < size()) resize(k); }
 173 private:
 174     void derive(); // in place version
 175 };
 176
 177 //TODO: figure out how to stick this in linear, while not adding an sbasis dep
 178 inline SBasis Linear::toSBasis() const { return SBasis(*this); }
 179
 180 //implemented in sbasis-roots.cpp
 181 OptInterval bounds_exact(SBasis const &a);
 182 OptInterval bounds_fast(SBasis const &a, int order = 0);
 183 OptInterval bounds_local(SBasis const &a, const OptInterval &t, int order = 0);
 184
 185 /** Returns a function which reverses the domain of a.
 186  \param a sbasis function
 187
 188 useful for reversing a parameteric curve.
 189 */
 190 inline SBasis reverse(SBasis const &a) {
 191     SBasis result(a.size(), Linear());
 192
 193     for(unsigned k = 0; k < a.size(); k++)
 194         result[k] = reverse(a[k]);
 195     return result;
 196 }
 197
 198 //IMPL: ScalableConcept
 199 inline SBasis operator-(const SBasis& p) {
 200     if(p.isZero()) return SBasis();
 201     SBasis result(p.size(), Linear());
 202
 203     for(unsigned i = 0; i < p.size(); i++) {
 204         result[i] = -p[i];
 205     }
 206     return result;
 207 }
 208 SBasis operator*(SBasis const &a, double k);
 209 inline SBasis operator*(double k, SBasis const &a) { return a*k; }
 210 inline SBasis operator/(SBasis const &a, double k) { return a*(1./k); }
 211 SBasis& operator*=(SBasis& a, double b);
 212 inline SBasis& operator/=(SBasis& a, double b) { return (a*=(1./b)); }
 213
 214 //IMPL: AddableConcept
 215 SBasis operator+(const SBasis& a, const SBasis& b);
 216 SBasis operator-(const SBasis& a, const SBasis& b);
 217 SBasis& operator+=(SBasis& a, const SBasis& b);
 218 SBasis& operator-=(SBasis& a, const SBasis& b);
 219
 220 //TODO: remove?
 221 /*inline SBasis operator+(const SBasis & a, Linear const & b) {
 222     if(b.isZero()) return a;
 223     if(a.isZero()) return b;
 224     SBasis result(a);
 225     result[0] += b;
 226     return result;
 227 }
 228 inline SBasis operator-(const SBasis & a, Linear const & b) {
 229     if(b.isZero()) return a;
 230     SBasis result(a);
 231     result[0] -= b;
 232     return result;
 233 }
 234 inline SBasis& operator+=(SBasis& a, const Linear& b) {
 235     if(a.isZero())
 236         a.push_back(b);
 237     else
 238         a[0] += b;
 239     return a;
 240 }
 241 inline SBasis& operator-=(SBasis& a, const Linear& b) {
 242     if(a.isZero())
 243         a.push_back(-b);
 244     else
 245         a[0] -= b;
 246     return a;
 247     }*/
 248
 249 //IMPL: OffsetableConcept
 250 inline SBasis operator+(const SBasis & a, double b) {
 251     if(a.isZero()) return Linear(b, b);
 252     SBasis result(a);
 253     result[0] += b;
 254     return result;
 255 }
 256 inline SBasis operator-(const SBasis & a, double b) {
 257     if(a.isZero()) return Linear(-b, -b);
 258     SBasis result(a);
 259     result[0] -= b;
 260     return result;
 261 }
 262 inline SBasis& operator+=(SBasis& a, double b) {
 263     if(a.isZero())
 264         a = SBasis(Linear(b,b));
 265     else
 266         a[0] += b;
 267     return a;
 268 }
 269 inline SBasis& operator-=(SBasis& a, double b) {
 270     if(a.isZero())
 271         a = SBasis(Linear(-b,-b));
 272     else
 273         a[0] -= b;
 274     return a;
 275 }
 276
 277 SBasis shift(SBasis const &a, int sh);
 278 SBasis shift(Linear const &a, int sh);
 279
 280 inline SBasis truncate(SBasis const &a, unsigned terms) {
 281     SBasis c;
 282     c.insert(c.begin(), a.begin(), a.begin() + std::min(terms, (unsigned)a.size()));
 283     return c;
 284 }
 285
 286 SBasis multiply(SBasis const &a, SBasis const &b);
 287 // This performs a multiply and accumulate operation in about the same time as multiply.  return a*b + c
 288 SBasis multiply_add(SBasis const &a, SBasis const &b, SBasis c);
 289
 290 SBasis integral(SBasis const &c);
 291 SBasis derivative(SBasis const &a);
 292
 293 SBasis sqrt(SBasis const &a, int k);
 294
 295 // return a kth order approx to 1/a)
 296 SBasis reciprocal(Linear const &a, int k);
 297 SBasis divide(SBasis const &a, SBasis const &b, int k);
 298
 299 inline SBasis operator*(SBasis const & a, SBasis const & b) {
 300     return multiply(a, b);
 301 }
 302
 303 inline SBasis& operator*=(SBasis& a, SBasis const & b) {
 304     a = multiply(a, b);
 305     return a;
 306 }
 307
 308 /** Returns the degree of the first non zero coefficient.
 309  \param a sbasis function
 310  \param tol largest abs val considered 0
 311  \returns first non zero coefficient
 312 */
 313 inline unsigned
 314 valuation(SBasis const &a, double tol=0){
 315     unsigned val=0;
 316     while( val<a.size() &&
 317            fabs(a[val][0])<tol &&
 318            fabs(a[val][1])<tol )
 319         val++;
 320     return val;
 321 }
 322
 323 // a(b(t))
 324 SBasis compose(SBasis const &a, SBasis const &b);
 325 SBasis compose(SBasis const &a, SBasis const &b, unsigned k);
 326 SBasis inverse(SBasis a, int k);
 327 //compose_inverse(f,g)=compose(f,inverse(g)), but is numerically more stable in some good cases...
 328 //TODO: requires g(0)=0 & g(1)=1 atm. generalization should be obvious.
 329 SBasis compose_inverse(SBasis const &f, SBasis const &g, unsigned order=2, double tol=1e-3);
 330
 331 /** Returns the sbasis on domain [0,1] that was t on [from, to]
 332  \param a sbasis function
 333  \param from,to interval
 334  \returns sbasis
 335
 336 */
 337 inline SBasis portion(const SBasis &t, double from, double to) { return compose(t, Linear(from, to)); }
 338 inline SBasis portion(const SBasis &t, Interval ivl) { return compose(t, Linear(ivl[0], ivl[1])); }
 339
 340 // compute f(g)
 341 inline SBasis
 342 SBasis::operator()(SBasis const & g) const {
 343     return compose(*this, g);
 344 }
 345
 346 inline std::ostream &operator<< (std::ostream &out_file, const Linear &bo) {
 347     out_file << "{" << bo[0] << ", " << bo[1] << "}";
 348     return out_file;
 349 }
 350
 351 inline std::ostream &operator<< (std::ostream &out_file, const SBasis & p) {
 352     for(unsigned i = 0; i < p.size(); i++) {
 353         out_file << p[i] << "s^" << i << " + ";
 354     }
 355     return out_file;
 356 }
 357
 358 // These are deprecated, use sbasis-math.h versions if possible
 359 SBasis sin(Linear bo, int k);
 360 SBasis cos(Linear bo, int k);
 361
 362 std::vector<double> roots(SBasis const & s);
 363 std::vector<std::vector<double> > multi_roots(SBasis const &f,
 364                                  std::vector<double> const &levels,
 365                                  double htol=1e-7,
 366                                  double vtol=1e-7,
 367                                  double a=0,
 368                                  double b=1);
 369
 370 }
 371 #endif
 372
 373 /*
 374   Local Variables:
 375   mode:c++
 376   c-file-style:"stroustrup"
 377   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 378   indent-tabs-mode:nil
 379   fill-column:99
 380   End:
 381 */
 382 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
 383 #endif