src/2geom/sbasis-2d.h

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