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     size_t size() const {return d.size();}
  73     Linear operator[](unsigned i) const {
  74         assert(i < size());
  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
  95
  96     SBasis() {}
  97     explicit SBasis(double a) {
  98         push_back(Linear(a,a));
  99     }
 100     SBasis(SBasis const & a) :
 101         d(a.d)
 102     {}
 103     SBasis(Linear const & bo) {
 104         push_back(bo);
 105     }
 106     SBasis(Linear* bo) {
 107         push_back(*bo);
 108     }
 109     explicit SBasis(size_t n, Linear const&l) : d(n, l) {}
 110
 111     //IMPL: FragmentConcept
 112     typedef double output_type;
 113     inline bool isZero() const {
 114         if(empty()) return true;
 115         for(unsigned i = 0; i < size(); i++) {
 116             if(!(*this)[i].isZero()) return false;
 117         }
 118         return true;
 119     }
 120     inline bool isConstant() const {
 121         if (empty()) return true;
 122         for (unsigned i = 0; i < size(); i++) {
 123             if(!(*this)[i].isConstant()) return false;
 124         }
 125         return true;
 126     }
 127
 128     bool isFinite() const;
 129     inline double at0() const {
 130         if(empty()) return 0; else return (*this)[0][0];
 131     }
 132     inline double at1() const{
 133         if(empty()) return 0; else return (*this)[0][1];
 134     }
 135
 136     int degreesOfFreedom() const { return size()*2;}
 137
 138     double valueAt(double t) const {
 139         double s = t*(1-t);
 140         double p0 = 0, p1 = 0;
 141         for(unsigned k = size(); k > 0; k--) {
 142             const Linear &lin = (*this)[k-1];
 143             p0 = p0*s + lin[0];
 144             p1 = p1*s + lin[1];
 145         }
 146         return (1-t)*p0 + t*p1;
 147     }
 148     //double valueAndDerivative(double t, double &der) const {
 149     //}
 150     double operator()(double t) const {
 151         return valueAt(t);
 152     }
 153
 154     std::vector<double> valueAndDerivatives(double t, unsigned n) const;
 155
 156     SBasis toSBasis() const { return SBasis(*this); }
 157
 158     double tailError(unsigned tail) const;
 159
 160 // compute f(g)
 161     SBasis operator()(SBasis const & g) const;
 162
 163 //MUTATOR PRISON
 164     //remove extra zeros
 165     void normalize() {
 166         while(!empty() && 0 == back()[0] && 0 == back()[1])
 167             pop_back();
 168     }
 169
 170     void truncate(unsigned k) { if(k < size()) resize(k); }
 171 private:
 172     void derive(); // in place version
 173 };
 174
 175 //TODO: figure out how to stick this in linear, while not adding an sbasis dep
 176 inline SBasis Linear::toSBasis() const { return SBasis(*this); }
 177
 178 //implemented in sbasis-roots.cpp
 179 OptInterval bounds_exact(SBasis const &a);
 180 OptInterval bounds_fast(SBasis const &a, int order = 0);
 181 OptInterval bounds_local(SBasis const &a, const OptInterval &t, int order = 0);
 182
 183 /** Returns a function which reverses the domain of a.
 184  \param a sbasis function
 185
 186 useful for reversing a parameteric curve.
 187 */
 188 inline SBasis reverse(SBasis const &a) {
 189     SBasis result(a.size(), Linear());
 190
 191     for(unsigned k = 0; k < a.size(); k++)
 192         result[k] = reverse(a[k]);
 193     return result;
 194 }
 195
 196 //IMPL: ScalableConcept
 197 inline SBasis operator-(const SBasis& p) {
 198     if(p.isZero()) return SBasis();
 199     SBasis result(p.size(), Linear());
 200
 201     for(unsigned i = 0; i < p.size(); i++) {
 202         result[i] = -p[i];
 203     }
 204     return result;
 205 }
 206 SBasis operator*(SBasis const &a, double k);
 207 inline SBasis operator*(double k, SBasis const &a) { return a*k; }
 208 inline SBasis operator/(SBasis const &a, double k) { return a*(1./k); }
 209 SBasis& operator*=(SBasis& a, double b);
 210 inline SBasis& operator/=(SBasis& a, double b) { return (a*=(1./b)); }
 211
 212 //IMPL: AddableConcept
 213 SBasis operator+(const SBasis& a, const SBasis& b);
 214 SBasis operator-(const SBasis& a, const SBasis& b);
 215 SBasis& operator+=(SBasis& a, const SBasis& b);
 216 SBasis& operator-=(SBasis& a, const SBasis& b);
 217
 218 //TODO: remove?
 219 /*inline SBasis operator+(const SBasis & a, Linear const & b) {
 220     if(b.isZero()) return a;
 221     if(a.isZero()) return b;
 222     SBasis result(a);
 223     result[0] += b;
 224     return result;
 225 }
 226 inline SBasis operator-(const SBasis & a, Linear const & b) {
 227     if(b.isZero()) return a;
 228     SBasis result(a);
 229     result[0] -= b;
 230     return result;
 231 }
 232 inline SBasis& operator+=(SBasis& a, const Linear& b) {
 233     if(a.isZero())
 234         a.push_back(b);
 235     else
 236         a[0] += b;
 237     return a;
 238 }
 239 inline SBasis& operator-=(SBasis& a, const Linear& b) {
 240     if(a.isZero())
 241         a.push_back(-b);
 242     else
 243         a[0] -= b;
 244     return a;
 245     }*/
 246
 247 //IMPL: OffsetableConcept
 248 inline SBasis operator+(const SBasis & a, double b) {
 249     if(a.isZero()) return Linear(b, b);
 250     SBasis result(a);
 251     result[0] += b;
 252     return result;
 253 }
 254 inline SBasis operator-(const SBasis & a, double b) {
 255     if(a.isZero()) return Linear(-b, -b);
 256     SBasis result(a);
 257     result[0] -= b;
 258     return result;
 259 }
 260 inline SBasis& operator+=(SBasis& a, double b) {
 261     if(a.isZero())
 262         a = SBasis(Linear(b,b));
 263     else
 264         a[0] += b;
 265     return a;
 266 }
 267 inline SBasis& operator-=(SBasis& a, double b) {
 268     if(a.isZero())
 269         a = SBasis(Linear(-b,-b));
 270     else
 271         a[0] -= b;
 272     return a;
 273 }
 274
 275 SBasis shift(SBasis const &a, int sh);
 276 SBasis shift(Linear const &a, int sh);
 277
 278 inline SBasis truncate(SBasis const &a, unsigned terms) {
 279     SBasis c;
 280     c.insert(c.begin(), a.begin(), a.begin() + std::min(terms, (unsigned)a.size()));
 281     return c;
 282 }
 283
 284 SBasis multiply(SBasis const &a, SBasis const &b);
 285 // This performs a multiply and accumulate operation in about the same time as multiply.  return a*b + c
 286 SBasis multiply_add(SBasis const &a, SBasis const &b, SBasis c);
 287
 288 SBasis integral(SBasis const &c);
 289 SBasis derivative(SBasis const &a);
 290
 291 SBasis sqrt(SBasis const &a, int k);
 292
 293 // return a kth order approx to 1/a)
 294 SBasis reciprocal(Linear const &a, int k);
 295 SBasis divide(SBasis const &a, SBasis const &b, int k);
 296
 297 inline SBasis operator*(SBasis const & a, SBasis const & b) {
 298     return multiply(a, b);
 299 }
 300
 301 inline SBasis& operator*=(SBasis& a, SBasis const & b) {
 302     a = multiply(a, b);
 303     return a;
 304 }
 305
 306 /** Returns the degree of the first non zero coefficient.
 307  \param a sbasis function
 308  \param tol largest abs val considered 0
 309  \returns first non zero coefficient
 310 */
 311 inline unsigned
 312 valuation(SBasis const &a, double tol=0){
 313     unsigned val=0;
 314     while( val<a.size() &&
 315            fabs(a[val][0])<tol &&
 316            fabs(a[val][1])<tol )
 317         val++;
 318     return val;
 319 }
 320
 321 // a(b(t))
 322 SBasis compose(SBasis const &a, SBasis const &b);
 323 SBasis compose(SBasis const &a, SBasis const &b, unsigned k);
 324 SBasis inverse(SBasis a, int k);
 325 //compose_inverse(f,g)=compose(f,inverse(g)), but is numerically more stable in some good cases...
 326 //TODO: requires g(0)=0 & g(1)=1 atm. generalization should be obvious.
 327 SBasis compose_inverse(SBasis const &f, SBasis const &g, unsigned order=2, double tol=1e-3);
 328
 329 /** Returns the sbasis on domain [0,1] that was t on [from, to]
 330  \param a sbasis function
 331  \param from,to interval
 332  \returns sbasis
 333
 334 */
 335 inline SBasis portion(const SBasis &t, double from, double to) { return compose(t, Linear(from, to)); }
 336 inline SBasis portion(const SBasis &t, Interval ivl) { return compose(t, Linear(ivl[0], ivl[1])); }
 337
 338 // compute f(g)
 339 inline SBasis
 340 SBasis::operator()(SBasis const & g) const {
 341     return compose(*this, g);
 342 }
 343
 344 inline std::ostream &operator<< (std::ostream &out_file, const Linear &bo) {
 345     out_file << "{" << bo[0] << ", " << bo[1] << "}";
 346     return out_file;
 347 }
 348
 349 inline std::ostream &operator<< (std::ostream &out_file, const SBasis & p) {
 350     for(unsigned i = 0; i < p.size(); i++) {
 351         out_file << p[i] << "s^" << i << " + ";
 352     }
 353     return out_file;
 354 }
 355
 356 // These are deprecated, use sbasis-math.h versions if possible
 357 SBasis sin(Linear bo, int k);
 358 SBasis cos(Linear bo, int k);
 359
 360 std::vector<double> roots(SBasis const & s);
 361 std::vector<std::vector<double> > multi_roots(SBasis const &f,
 362                                  std::vector<double> const &levels,
 363                                  double htol=1e-7,
 364                                  double vtol=1e-7,
 365                                  double a=0,
 366                                  double b=1);
 367
 368 }
 369 #endif
 370
 371 /*
 372   Local Variables:
 373   mode:c++
 374   c-file-style:"stroustrup"
 375   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 376   indent-tabs-mode:nil
 377   fill-column:99
 378   End:
 379 */
 380 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
 381 #endif