src/2geom/sbasis-2d.h

   1 #ifndef SEEN_SBASIS_2D_H
   2 #define SEEN_SBASIS_2D_H
   3 #include <vector>
   4 #include <cassert>
   5 #include <algorithm>
   6 #include <2geom/d2.h>
   7 #include <2geom/sbasis.h>
   8 #include <iostream>
   9
  10 namespace Geom{
  11
  12 class Linear2d{
  13 public:
  14     /*
  15         u 0,1
  16         v 0,2
  17     */
  18     double a[4];
  19     Linear2d() {}
  20     Linear2d(double aa) {
  21         for(unsigned i = 0 ; i < 4; i ++)
  22             a[i] = aa;
  23     }
  24     Linear2d(double a00, double a01, double a10, double a11)
  25     {
  26         a[0] = a00;
  27         a[1] = a01;
  28         a[2] = a10;
  29         a[3] = a11;
  30     }
  31
  32     double operator[](const int i) const {
  33         assert(i >= 0);
  34         assert(i < 4);
  35         return a[i];
  36     }
  37     double& operator[](const int i) {
  38         assert(i >= 0);
  39         assert(i < 4);
  40         return a[i];
  41     }
  42     double apply(double u, double v) {
  43         return (a[0]*(1-u)*(1-v) +
  44                 a[1]*u*(1-v) +
  45                 a[2]*(1-u)*v +
  46                 a[3]*u*v);
  47     }
  48 };
  49
  50 inline Linear extract_u(Linear2d const &a, double u) {
  51     return Linear(a[0]*(1-u) +
  52                   a[1]*u,
  53                   a[2]*(1-u) +
  54                   a[3]*u);
  55 }
  56 inline Linear extract_v(Linear2d const &a, double v) {
  57     return Linear(a[0]*(1-v) +
  58                   a[2]*v,
  59                   a[1]*(1-v) +
  60                   a[3]*v);
  61 }
  62 inline Linear2d operator-(Linear2d const &a) {
  63     return Linear2d(-a.a[0], -a.a[1],
  64                     -a.a[2], -a.a[3]);
  65 }
  66 inline Linear2d operator+(Linear2d const & a, Linear2d const & b) {
  67     return Linear2d(a[0] + b[0],
  68                   a[1] + b[1],
  69                   a[2] + b[2],
  70                   a[3] + b[3]);
  71 }
  72 inline Linear2d operator-(Linear2d const & a, Linear2d const & b) {
  73     return Linear2d(a[0] - b[0],
  74                   a[1] - b[1],
  75                   a[2] - b[2],
  76                   a[3] - b[3]);
  77 }
  78 inline Linear2d& operator+=(Linear2d & a, Linear2d const & b) {
  79     for(unsigned i = 0; i < 4; i++)
  80         a[i] += b[i];
  81     return a;
  82 }
  83 inline Linear2d& operator-=(Linear2d & a, Linear2d const & b) {
  84     for(unsigned i = 0; i < 4; i++)
  85         a[i] -= b[i];
  86     return a;
  87 }
  88 inline Linear2d& operator*=(Linear2d & a, double b) {
  89     for(unsigned i = 0; i < 4; i++)
  90         a[i] *= b;
  91     return a;
  92 }
  93
  94 inline bool operator==(Linear2d const & a, Linear2d const & b) {
  95     for(unsigned i = 0; i < 4; i++)
  96         if(a[i] != b[i])
  97             return false;
  98     return true;
  99 }
 100 inline bool operator!=(Linear2d const & a, Linear2d const & b) {
 101     for(unsigned i = 0; i < 4; i++)
 102         if(a[i] == b[i])
 103             return false;
 104     return true;
 105 }
 106 inline Linear2d operator*(double const a, Linear2d const & b) {
 107     return Linear2d(a*b[0], a*b[1],
 108                   a*b[2], a*b[3]);
 109 }
 110
 111 class SBasis2d : public std::vector<Linear2d>{
 112 public:
 113     // vector in u,v
 114     unsigned us, vs; // number of u terms, v terms
 115     SBasis2d() {}
 116     SBasis2d(Linear2d const & bo)
 117         : us(1), vs(1) {
 118         push_back(bo);
 119     }
 120     SBasis2d(SBasis2d const & a)
 121         : std::vector<Linear2d>(a), us(a.us), vs(a.vs) {}
 122
 123     Linear2d& index(unsigned ui, unsigned vi) {
 124         assert(ui < us);
 125         assert(vi < vs);
 126         return (*this)[ui + vi*us];
 127     }
 128
 129     Linear2d index(unsigned ui, unsigned vi) const {
 130         if(ui >= us)
 131             return Linear2d(0);
 132         if(vi >= vs)
 133             return Linear2d(0);
 134         return (*this)[ui + vi*us];
 135     }
 136
 137     double apply(double u, double v) const {
 138         double s = u*(1-u);
 139         double t = v*(1-v);
 140         Linear2d p;
 141         double tk = 1;
 142 // XXX rewrite as horner
 143         for(unsigned vi = 0; vi < vs; vi++) {
 144             double sk = 1;
 145             for(unsigned ui = 0; ui < us; ui++) {
 146                 p += (sk*tk)*index(ui, vi);
 147                 sk *= s;
 148             }
 149             tk *= t;
 150         }
 151         return p.apply(u,v);
 152     }
 153
 154     void clear() {
 155         fill(begin(), end(), Linear2d(0));
 156     }
 157
 158     void normalize(); // remove extra zeros
 159
 160     double tail_error(unsigned tail) const;
 161
 162     void truncate(unsigned k);
 163 };
 164
 165 inline SBasis2d operator-(const SBasis2d& p) {
 166     SBasis2d result;
 167     result.reserve(p.size());
 168
 169     for(unsigned i = 0; i < p.size(); i++) {
 170         result.push_back(-p[i]);
 171     }
 172     return result;
 173 }
 174
 175 inline SBasis2d operator+(const SBasis2d& a, const SBasis2d& b) {
 176     SBasis2d result;
 177     result.us = std::max(a.us, b.us);
 178     result.vs = std::max(a.vs, b.vs);
 179     const unsigned out_size = result.us*result.vs;
 180     result.resize(out_size);
 181
 182     for(unsigned vi = 0; vi < result.vs; vi++) {
 183         for(unsigned ui = 0; ui < result.us; ui++) {
 184             Linear2d bo;
 185             if(ui < a.us && vi < a.vs)
 186                 bo += a.index(ui, vi);
 187             if(ui < b.us && vi < b.vs)
 188                 bo += b.index(ui, vi);
 189             result.index(ui, vi) = bo;
 190         }
 191     }
 192     return result;
 193 }
 194
 195 inline SBasis2d operator-(const SBasis2d& a, const SBasis2d& b) {
 196     SBasis2d result;
 197     result.us = std::max(a.us, b.us);
 198     result.vs = std::max(a.vs, b.vs);
 199     const unsigned out_size = result.us*result.vs;
 200     result.resize(out_size);
 201
 202     for(unsigned vi = 0; vi < result.vs; vi++) {
 203         for(unsigned ui = 0; ui < result.us; ui++) {
 204             Linear2d bo;
 205             if(ui < a.us && vi < a.vs)
 206                 bo += a.index(ui, vi);
 207             if(ui < b.us && vi < b.vs)
 208                 bo -= b.index(ui, vi);
 209             result.index(ui, vi) = bo;
 210         }
 211     }
 212     return result;
 213 }
 214
 215
 216 inline SBasis2d& operator+=(SBasis2d& a, const Linear2d& b) {
 217     if(a.size() < 1)
 218         a.push_back(b);
 219     else
 220         a[0] += b;
 221     return a;
 222 }
 223
 224 inline SBasis2d& operator-=(SBasis2d& a, const Linear2d& b) {
 225     if(a.size() < 1)
 226         a.push_back(-b);
 227     else
 228         a[0] -= b;
 229     return a;
 230 }
 231
 232 inline SBasis2d& operator+=(SBasis2d& a, double b) {
 233     if(a.size() < 1)
 234         a.push_back(Linear2d(b));
 235     else {
 236         for(unsigned i = 0; i < 4; i++)
 237             a[0] += double(b);
 238     }
 239     return a;
 240 }
 241
 242 inline SBasis2d& operator-=(SBasis2d& a, double b) {
 243     if(a.size() < 1)
 244         a.push_back(Linear2d(-b));
 245     else {
 246         a[0] -= b;
 247     }
 248     return a;
 249 }
 250
 251 inline SBasis2d& operator*=(SBasis2d& a, double b) {
 252     for(unsigned i = 0; i < a.size(); i++)
 253         a[i] *= b;
 254     return a;
 255 }
 256
 257 inline SBasis2d& operator/=(SBasis2d& a, double b) {
 258     for(unsigned i = 0; i < a.size(); i++)
 259         a[i] *= (1./b);
 260     return a;
 261 }
 262
 263 SBasis2d operator*(double k, SBasis2d const &a);
 264 SBasis2d operator*(SBasis2d const &a, SBasis2d const &b);
 265
 266 SBasis2d shift(SBasis2d const &a, int sh);
 267
 268 SBasis2d shift(Linear2d const &a, int sh);
 269
 270 SBasis2d truncate(SBasis2d const &a, unsigned terms);
 271
 272 SBasis2d multiply(SBasis2d const &a, SBasis2d const &b);
 273
 274 SBasis2d integral(SBasis2d const &c);
 275
 276 SBasis2d derivative(SBasis2d const &a);
 277
 278 SBasis2d sqrt(SBasis2d const &a, int k);
 279
 280 // return a kth order approx to 1/a)
 281 SBasis2d reciprocal(Linear2d const &a, int k);
 282
 283 SBasis2d divide(SBasis2d const &a, SBasis2d const &b, int k);
 284
 285 // a(b(t))
 286 SBasis2d compose(SBasis2d const &a, SBasis2d const &b);
 287 SBasis2d compose(SBasis2d const &a, SBasis2d const &b, unsigned k);
 288 SBasis2d inverse(SBasis2d const &a, int k);
 289
 290 // these two should probably be replaced with compose
 291 SBasis extract_u(SBasis2d const &a, double u);
 292 SBasis extract_v(SBasis2d const &a, double v);
 293
 294 SBasis compose(Linear2d const &a, D2<SBasis> const &p);
 295
 296 SBasis compose(SBasis2d const &fg, D2<SBasis> const &p);
 297
 298 D2<SBasis> compose_each(D2<SBasis2d> const &fg, D2<SBasis> const &p);
 299
 300 inline std::ostream &operator<< (std::ostream &out_file, const Linear2d &bo) {
 301     out_file << "{" << bo[0] << ", " << bo[1] << "}, ";
 302     out_file << "{" << bo[2] << ", " << bo[3] << "}";
 303     return out_file;
 304 }
 305
 306 inline std::ostream &operator<< (std::ostream &out_file, const SBasis2d & p) {
 307     for(unsigned i = 0; i < p.size(); i++) {
 308         out_file << p[i] << "s^" << i << " + ";
 309     }
 310     return out_file;
 311 }
 312
 313 };
 314
 315 /*
 316   Local Variables:
 317   mode:c++
 318   c-file-style:"stroustrup"
 319   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 320   indent-tabs-mode:nil
 321   fill-column:99
 322   End:
 323 */
 324 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
 325 #endif