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