src/2geom/bezier.h

   1 /**
   2  * \file bezier.h
   3  * \brief \todo brief description
   4  *
   5  * Copyright 2007  MenTaLguY <mental@rydia.net>
   6  * Copyright 2007  Michael Sloan <mgsloan@gmail.com>
   7  * Copyright 2007  Nathan Hurst <njh@njhurst.com>
   8  *
   9  * This library is free software; you can redistribute it and/or
  10  * modify it either under the terms of the GNU Lesser General Public
  11  * License version 2.1 as published by the Free Software Foundation
  12  * (the "LGPL") or, at your option, under the terms of the Mozilla
  13  * Public License Version 1.1 (the "MPL"). If you do not alter this
  14  * notice, a recipient may use your version of this file under either
  15  * the MPL or the LGPL.
  16  *
  17  * You should have received a copy of the LGPL along with this library
  18  * in the file COPYING-LGPL-2.1; if not, write to the Free Software
  19  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
  20  * You should have received a copy of the MPL along with this library
  21  * in the file COPYING-MPL-1.1
  22  *
  23  * The contents of this file are subject to the Mozilla Public License
  24  * Version 1.1 (the "License"); you may not use this file except in
  25  * compliance with the License. You may obtain a copy of the License at
  26  * http://www.mozilla.org/MPL/
  27  *
  28  * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
  29  * OF ANY KIND, either express or implied. See the LGPL or the MPL for
  30  * the specific language governing rights and limitations.
  31  *
  32  */
  33
  34 #ifndef SEEN_BEZIER_H
  35 #define SEEN_BEZIER_H
  36
  37 #include <2geom/coord.h>
  38 #include <valarray>
  39 #include <2geom/isnan.h>
  40 #include <2geom/d2.h>
  41 #include <2geom/solver.h>
  42 #include <boost/optional/optional.hpp>
  43
  44 namespace Geom {
  45
  46 inline Coord subdivideArr(Coord t, Coord const *v, Coord *left, Coord *right, unsigned order) {
  47 /*
  48  *  Bernstein :
  49  *      Evaluate a Bernstein function at a particular parameter value
  50  *      Fill in control points for resulting sub-curves.
  51  *
  52  */
  53
  54     unsigned N = order+1;
  55     std::valarray<Coord> row(N);
  56     for (unsigned i = 0; i < N; i++)
  57         row[i] = v[i];
  58
  59     // Triangle computation
  60     const double omt = (1-t);
  61     if(left)
  62         left[0] = row[0];
  63     if(right)
  64         right[order] = row[order];
  65     for (unsigned i = 1; i < N; i++) {
  66         for (unsigned j = 0; j < N - i; j++) {
  67             row[j] = omt*row[j] + t*row[j+1];
  68         }
  69         if(left)
  70             left[i] = row[0];
  71         if(right)
  72             right[order-i] = row[order-i];
  73     }
  74     return (row[0]);
  75 /*
  76     Coord vtemp[order+1][order+1];
  77
  78     // Copy control points
  79     std::copy(v, v+order+1, vtemp[0]);
  80
  81     // Triangle computation
  82     for (unsigned i = 1; i <= order; i++) {
  83         for (unsigned j = 0; j <= order - i; j++) {
  84             vtemp[i][j] = lerp(t, vtemp[i-1][j], vtemp[i-1][j+1]);
  85         }
  86     }
  87     if(left != NULL)
  88         for (unsigned j = 0; j <= order; j++)
  89             left[j]  = vtemp[j][0];
  90     if(right != NULL)
  91         for (unsigned j = 0; j <= order; j++)
  92             right[j] = vtemp[order-j][j];
  93
  94             return (vtemp[order][0]);*/
  95 }
  96
  97 template <typename T>
  98 inline T bernsteinValueAt(double t, T const *c_, unsigned n) {
  99     double u = 1.0 - t;
 100     double bc = 1;
 101     double tn = 1;
 102     T tmp = c_[0]*u;
 103     for(unsigned i=1; i<n; i++){
 104         tn = tn*t;
 105         bc = bc*(n-i+1)/i;
 106         tmp = (tmp + tn*bc*c_[i])*u;
 107     }
 108     return (tmp + tn*t*c_[n]);
 109 }
 110
 111
 112 class Bezier {
 113 private:
 114     std::valarray<Coord> c_;
 115
 116     friend Bezier portion(const Bezier & a, Coord from, Coord to);
 117
 118     friend OptInterval bounds_fast(Bezier const & b);
 119
 120     friend Bezier derivative(const Bezier & a);
 121
 122 protected:
 123     Bezier(Coord const c[], unsigned ord) : c_(c, ord+1){
 124         //std::copy(c, c+order()+1, &c_[0]);
 125     }
 126
 127 public:
 128     unsigned int order() const { return c_.size()-1;}
 129     unsigned int size() const { return c_.size();}
 130
 131     Bezier() {}
 132     Bezier(const Bezier& b) :c_(b.c_) {}
 133     Bezier &operator=(Bezier const &other) {
 134         if ( c_.size() != other.c_.size() ) {
 135             c_.resize(other.c_.size());
 136         }
 137         c_ = other.c_;
 138         return *this;
 139     }
 140
 141     struct Order {
 142         unsigned order;
 143         explicit Order(Bezier const &b) : order(b.order()) {}
 144         explicit Order(unsigned o) : order(o) {}
 145         operator unsigned() const { return order; }
 146     };
 147
 148     //Construct an arbitrary order bezier
 149     Bezier(Order ord) : c_(0., ord.order+1) {
 150         assert(ord.order ==  order());
 151     }
 152
 153     explicit Bezier(Coord c0) : c_(0., 1) {
 154         c_[0] = c0;
 155     }
 156
 157     //Construct an order-1 bezier (linear Bézier)
 158     Bezier(Coord c0, Coord c1) : c_(0., 2) {
 159         c_[0] = c0; c_[1] = c1;
 160     }
 161
 162     //Construct an order-2 bezier (quadratic Bézier)
 163     Bezier(Coord c0, Coord c1, Coord c2) : c_(0., 3) {
 164         c_[0] = c0; c_[1] = c1; c_[2] = c2;
 165     }
 166
 167     //Construct an order-3 bezier (cubic Bézier)
 168     Bezier(Coord c0, Coord c1, Coord c2, Coord c3) : c_(0., 4) {
 169         c_[0] = c0; c_[1] = c1; c_[2] = c2; c_[3] = c3;
 170     }
 171
 172     void resize (unsigned int n, Coord v = 0)
 173     {
 174         c_.resize (n, v);
 175     }
 176
 177     void clear()
 178     {
 179         c_.resize(0);
 180     }
 181
 182     inline unsigned degree() const { return order(); }
 183
 184     //IMPL: FragmentConcept
 185     typedef Coord output_type;
 186     inline bool isZero() const {
 187         for(unsigned i = 0; i <= order(); i++) {
 188             if(c_[i] != 0) return false;
 189         }
 190         return true;
 191     }
 192     inline bool isConstant() const {
 193         for(unsigned i = 1; i <= order(); i++) {
 194             if(c_[i] != c_[0]) return false;
 195         }
 196         return true;
 197     }
 198     inline bool isFinite() const {
 199         for(unsigned i = 0; i <= order(); i++) {
 200             if(!IS_FINITE(c_[i])) return false;
 201         }
 202         return true;
 203     }
 204     inline Coord at0() const { return c_[0]; }
 205     inline Coord at1() const { return c_[order()]; }
 206
 207     inline Coord valueAt(double t) const {
 208         int n = order();
 209         double u, bc, tn, tmp;
 210         int i;
 211         u = 1.0 - t;
 212         bc = 1;
 213         tn = 1;
 214         tmp = c_[0]*u;
 215         for(i=1; i<n; i++){
 216             tn = tn*t;
 217             bc = bc*(n-i+1)/i;
 218             tmp = (tmp + tn*bc*c_[i])*u;
 219         }
 220         return (tmp + tn*t*c_[n]);
 221         //return subdivideArr(t, &c_[0], NULL, NULL, order());
 222     }
 223     inline Coord operator()(double t) const { return valueAt(t); }
 224
 225     SBasis toSBasis() const;
 226 //    inline SBasis toSBasis() const {
 227 //        SBasis sb;
 228 //        bezier_to_sbasis(sb, (*this));
 229 //        return sb;
 230 //        //return bezier_to_sbasis(&c_[0], order());
 231 //    }
 232
 233     //Only mutator
 234     inline Coord &operator[](unsigned ix) { return c_[ix]; }
 235     inline Coord const &operator[](unsigned ix) const { return const_cast<std::valarray<Coord>&>(c_)[ix]; }
 236     //inline Coord const &operator[](unsigned ix) const { return c_[ix]; }
 237     inline void setPoint(unsigned ix, double val) { c_[ix] = val; }
 238
 239     /**
 240     *  The size of the returned vector equals n_derivs+1.
 241     */
 242     std::vector<Coord> valueAndDerivatives(Coord t, unsigned n_derivs) const {
 243         /* This is inelegant, as it uses several extra stores.  I think there might be a way to
 244          * evaluate roughly in situ. */
 245
 246          // initialize return vector with zeroes, such that we only need to replace the non-zero derivs
 247         std::vector<Coord> val_n_der(n_derivs + 1, Coord(0.0));
 248
 249         // initialize temp storage variables
 250         std::valarray<Coord> d_(order()+1);
 251         for (unsigned i = 0; i < size(); i++) {
 252             d_[i] = c_[i];
 253         }
 254
 255         unsigned nn = n_derivs + 1;
 256         if(n_derivs > order()) {
 257             nn = order()+1; // only calculate the non zero derivs
 258         }
 259         for (unsigned di = 0; di < nn; di++) {
 260             //val_n_der[di] = (subdivideArr(t, &d_[0], NULL, NULL, order() - di));
 261             val_n_der[di] = bernsteinValueAt(t, &d_[0], order() - di);
 262             for (unsigned i = 0; i < order() - di; i++) {
 263                 d_[i] = (order()-di)*(d_[i+1] - d_[i]);
 264             }
 265         }
 266
 267         return val_n_der;
 268     }
 269
 270     std::pair<Bezier, Bezier > subdivide(Coord t) const {
 271         Bezier a(Bezier::Order(*this)), b(Bezier::Order(*this));
 272         subdivideArr(t, &const_cast<std::valarray<Coord>&>(c_)[0], &a.c_[0], &b.c_[0], order());
 273         return std::pair<Bezier, Bezier >(a, b);
 274     }
 275
 276     std::vector<double> roots() const {
 277         std::vector<double> solutions;
 278         find_bernstein_roots(&const_cast<std::valarray<Coord>&>(c_)[0], order(), solutions, 0, 0.0, 1.0);
 279         return solutions;
 280     }
 281     std::vector<double> roots(Interval const ivl) const {
 282         std::vector<double> solutions;
 283         find_bernstein_roots(&const_cast<std::valarray<Coord>&>(c_)[0], order(), solutions, 0, ivl[0], ivl[1]);
 284         return solutions;
 285     }
 286 };
 287
 288
 289 void bezier_to_sbasis (SBasis & sb, Bezier const& bz);
 290
 291
 292 inline
 293 SBasis Bezier::toSBasis() const {
 294     SBasis sb;
 295     bezier_to_sbasis(sb, (*this));
 296     return sb;
 297     //return bezier_to_sbasis(&c_[0], order());
 298 }
 299
 300 //TODO: implement others
 301 inline Bezier operator+(const Bezier & a, double v) {
 302     Bezier result = Bezier(Bezier::Order(a));
 303     for(unsigned i = 0; i <= a.order(); i++)
 304         result[i] = a[i] + v;
 305     return result;
 306 }
 307
 308 inline Bezier operator-(const Bezier & a, double v) {
 309     Bezier result = Bezier(Bezier::Order(a));
 310     for(unsigned i = 0; i <= a.order(); i++)
 311         result[i] = a[i] - v;
 312     return result;
 313 }
 314
 315 inline Bezier operator*(const Bezier & a, double v) {
 316     Bezier result = Bezier(Bezier::Order(a));
 317     for(unsigned i = 0; i <= a.order(); i++)
 318         result[i] = a[i] * v;
 319     return result;
 320 }
 321
 322 inline Bezier operator/(const Bezier & a, double v) {
 323     Bezier result = Bezier(Bezier::Order(a));
 324     for(unsigned i = 0; i <= a.order(); i++)
 325         result[i] = a[i] / v;
 326     return result;
 327 }
 328
 329 inline Bezier reverse(const Bezier & a) {
 330     Bezier result = Bezier(Bezier::Order(a));
 331     for(unsigned i = 0; i <= a.order(); i++)
 332         result[i] = a[a.order() - i];
 333     return result;
 334 }
 335
 336 inline Bezier portion(const Bezier & a, double from, double to) {
 337     //TODO: implement better?
 338     std::valarray<Coord> res(a.order() + 1);
 339     if(from == 0) {
 340         if(to == 1) { return Bezier(a); }
 341         subdivideArr(to, &const_cast<Bezier&>(a).c_[0], &res[0], NULL, a.order());
 342         return Bezier(&res[0], a.order());
 343     }
 344     subdivideArr(from, &const_cast<Bezier&>(a).c_[0], NULL, &res[0], a.order());
 345     if(to == 1) return Bezier(&res[0], a.order());
 346     std::valarray<Coord> res2(a.order()+1);
 347     subdivideArr((to - from)/(1 - from), &res[0], &res2[0], NULL, a.order());
 348     return Bezier(&res2[0], a.order());
 349 }
 350
 351 // XXX Todo: how to handle differing orders
 352 inline std::vector<Point> bezier_points(const D2<Bezier > & a) {
 353     std::vector<Point> result;
 354     for(unsigned i = 0; i <= a[0].order(); i++) {
 355         Point p;
 356         for(unsigned d = 0; d < 2; d++) p[d] = a[d][i];
 357         result.push_back(p);
 358     }
 359     return result;
 360 }
 361
 362 inline Bezier derivative(const Bezier & a) {
 363     //if(a.order() == 1) return Bezier(0.0);
 364     if(a.order() == 1) return Bezier(a.c_[1]-a.c_[0]);
 365     Bezier der(Bezier::Order(a.order()-1));
 366
 367     for(unsigned i = 0; i < a.order(); i++) {
 368         der.c_[i] = a.order()*(a.c_[i+1] - a.c_[i]);
 369     }
 370     return der;
 371 }
 372
 373 inline Bezier integral(const Bezier & a) {
 374     Bezier inte(Bezier::Order(a.order()+1));
 375
 376     inte[0] = 0;
 377     for(unsigned i = 0; i < inte.order(); i++) {
 378         inte[i+1] = inte[i] + a[i]/(inte.order());
 379     }
 380     return inte;
 381 }
 382
 383 inline OptInterval bounds_fast(Bezier const & b) {
 384     return Interval::fromArray(&const_cast<Bezier&>(b).c_[0], b.size());
 385 }
 386
 387 //TODO: better bounds exact
 388 inline OptInterval bounds_exact(Bezier const & b) {
 389     return bounds_exact(b.toSBasis());
 390 }
 391
 392 inline OptInterval bounds_local(Bezier const & b, OptInterval i) {
 393     //return bounds_local(b.toSBasis(), i);
 394     if (i) {
 395         return bounds_fast(portion(b, i->min(), i->max()));
 396     } else {
 397         return OptInterval();
 398     }
 399 }
 400
 401 inline std::ostream &operator<< (std::ostream &out_file, const Bezier & b) {
 402     for(unsigned i = 0; i < b.size(); i++) {
 403         out_file << b[i] << ", ";
 404     }
 405     return out_file;
 406 }
 407
 408 }
 409 #endif //SEEN_BEZIER_H
 410
 411 /*
 412   Local Variables:
 413   mode:c++
 414   c-file-style:"stroustrup"
 415   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 416   indent-tabs-mode:nil
 417   fill-column:99
 418   End:
 419 */
 420 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :