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     Coord vtemp[order+1][order+1];
  48
  49     /* Copy control points      */
  50     std::copy(v, v+order+1, vtemp[0]);
  51
  52     /* Triangle computation     */
  53     for (unsigned i = 1; i <= order; i++) {
  54         for (unsigned j = 0; j <= order - i; j++) {
  55             vtemp[i][j] = lerp(t, vtemp[i-1][j], vtemp[i-1][j+1]);
  56         }
  57     }
  58     if(left != NULL)
  59         for (unsigned j = 0; j <= order; j++)
  60             left[j]  = vtemp[j][0];
  61     if(right != NULL)
  62         for (unsigned j = 0; j <= order; j++)
  63             right[j] = vtemp[order-j][j];
  64
  65     return (vtemp[order][0]);
  66 }
  67
  68
  69 class Bezier {
  70 private:
  71     std::valarray<Coord> c_;
  72
  73     friend Bezier portion(const Bezier & a, Coord from, Coord to);
  74
  75     friend OptInterval bounds_fast(Bezier const & b);
  76
  77     friend Bezier derivative(const Bezier & a);
  78
  79 protected:
  80     Bezier(Coord const c[], unsigned ord) : c_(c, ord+1){
  81         //std::copy(c, c+order()+1, &c_[0]);
  82     }
  83
  84 public:
  85     unsigned int order() const { return c_.size()-1;}
  86     unsigned int size() const { return c_.size();}
  87
  88     Bezier() {}
  89     Bezier(const Bezier& b) :c_(b.c_) {}
  90     Bezier &operator=(Bezier const &other) {
  91         if ( c_.size() != other.c_.size() ) {
  92             c_.resize(other.c_.size());
  93         }
  94         c_ = other.c_;
  95         return *this;
  96     }
  97
  98     struct Order {
  99         unsigned order;
 100         explicit Order(Bezier const &b) : order(b.order()) {}
 101         explicit Order(unsigned o) : order(o) {}
 102         operator unsigned() const { return order; }
 103     };
 104
 105     //Construct an arbitrary order bezier
 106     Bezier(Order ord) : c_(0., ord.order+1) {
 107         assert(ord.order ==  order());
 108     }
 109
 110     explicit Bezier(Coord c0) : c_(0., 1) {
 111         c_[0] = c0;
 112     }
 113
 114     //Construct an order-1 bezier (linear Bézier)
 115     Bezier(Coord c0, Coord c1) : c_(0., 2) {
 116         c_[0] = c0; c_[1] = c1;
 117     }
 118
 119     //Construct an order-2 bezier (quadratic Bézier)
 120     Bezier(Coord c0, Coord c1, Coord c2) : c_(0., 3) {
 121         c_[0] = c0; c_[1] = c1; c_[2] = c2;
 122     }
 123
 124     //Construct an order-3 bezier (cubic Bézier)
 125     Bezier(Coord c0, Coord c1, Coord c2, Coord c3) : c_(0., 4) {
 126         c_[0] = c0; c_[1] = c1; c_[2] = c2; c_[3] = c3;
 127     }
 128
 129     void resize (unsigned int n, Coord v = 0)
 130     {
 131         c_.resize (n, v);
 132     }
 133
 134     void clear()
 135     {
 136         c_.resize(0);
 137     }
 138
 139     inline unsigned degree() const { return order(); }
 140
 141     //IMPL: FragmentConcept
 142     typedef Coord output_type;
 143     inline bool isZero() const {
 144         for(unsigned i = 0; i <= order(); i++) {
 145             if(c_[i] != 0) return false;
 146         }
 147         return true;
 148     }
 149     inline bool isConstant() const {
 150         for(unsigned i = 1; i <= order(); i++) {
 151             if(c_[i] != c_[0]) return false;
 152         }
 153         return true;
 154     }
 155     inline bool isFinite() const {
 156         for(unsigned i = 0; i <= order(); i++) {
 157             if(!IS_FINITE(c_[i])) return false;
 158         }
 159         return true;
 160     }
 161     inline Coord at0() const { return c_[0]; }
 162     inline Coord at1() const { return c_[order()]; }
 163
 164     inline Coord valueAt(double t) const {
 165         return subdivideArr(t, &c_[0], NULL, NULL, order());
 166     }
 167     inline Coord operator()(double t) const { return valueAt(t); }
 168
 169     SBasis toSBasis() const;
 170 //    inline SBasis toSBasis() const {
 171 //        SBasis sb;
 172 //        bezier_to_sbasis(sb, (*this));
 173 //        return sb;
 174 //        //return bezier_to_sbasis(&c_[0], order());
 175 //    }
 176
 177     //Only mutator
 178     inline Coord &operator[](unsigned ix) { return c_[ix]; }
 179     inline Coord const &operator[](unsigned ix) const { return c_[ix]; }
 180     inline void setPoint(unsigned ix, double val) { c_[ix] = val; }
 181
 182     /* This is inelegant, as it uses several extra stores.  I think there might be a way to
 183      * evaluate roughly in situ. */
 184
 185     std::vector<Coord> valueAndDerivatives(Coord t, unsigned n_derivs) const {
 186         std::vector<Coord> val_n_der;
 187         Coord d_[order()+1];
 188         unsigned nn = n_derivs + 1;     // the size of the result vector equals n_derivs+1 ...
 189         if(nn > order())
 190             nn = order()+1;                     // .. but with a maximum of order() + 1!
 191         for(unsigned i = 0; i < size(); i++)
 192             d_[i] = c_[i];
 193         for(unsigned di = 0; di < nn; di++) {
 194             val_n_der.push_back(subdivideArr(t, d_, NULL, NULL, order() - di));
 195             for(unsigned i = 0; i < order() - di; i++) {
 196                 d_[i] = (order()-di)*(d_[i+1] - d_[i]);
 197             }
 198         }
 199         val_n_der.resize(n_derivs);
 200         return val_n_der;
 201     }
 202
 203     std::pair<Bezier, Bezier > subdivide(Coord t) const {
 204         Bezier a(Bezier::Order(*this)), b(Bezier::Order(*this));
 205         subdivideArr(t, &c_[0], &a.c_[0], &b.c_[0], order());
 206         return std::pair<Bezier, Bezier >(a, b);
 207     }
 208
 209     std::vector<double> roots() const {
 210         std::vector<double> solutions;
 211         find_bernstein_roots(&c_[0], order(), solutions, 0, 0.0, 1.0);
 212         return solutions;
 213     }
 214 };
 215
 216
 217 void bezier_to_sbasis (SBasis & sb, Bezier const& bz);
 218
 219
 220 inline
 221 SBasis Bezier::toSBasis() const {
 222     SBasis sb;
 223     bezier_to_sbasis(sb, (*this));
 224     return sb;
 225     //return bezier_to_sbasis(&c_[0], order());
 226 }
 227
 228 //TODO: implement others
 229 inline Bezier operator+(const Bezier & a, double v) {
 230     Bezier result = Bezier(Bezier::Order(a));
 231     for(unsigned i = 0; i <= a.order(); i++)
 232         result[i] = a[i] + v;
 233     return result;
 234 }
 235
 236 inline Bezier operator-(const Bezier & a, double v) {
 237     Bezier result = Bezier(Bezier::Order(a));
 238     for(unsigned i = 0; i <= a.order(); i++)
 239         result[i] = a[i] - v;
 240     return result;
 241 }
 242
 243 inline Bezier operator*(const Bezier & a, double v) {
 244     Bezier result = Bezier(Bezier::Order(a));
 245     for(unsigned i = 0; i <= a.order(); i++)
 246         result[i] = a[i] * v;
 247     return result;
 248 }
 249
 250 inline Bezier operator/(const Bezier & a, double v) {
 251     Bezier result = Bezier(Bezier::Order(a));
 252     for(unsigned i = 0; i <= a.order(); i++)
 253         result[i] = a[i] / v;
 254     return result;
 255 }
 256
 257 inline Bezier reverse(const Bezier & a) {
 258     Bezier result = Bezier(Bezier::Order(a));
 259     for(unsigned i = 0; i <= a.order(); i++)
 260         result[i] = a[a.order() - i];
 261     return result;
 262 }
 263
 264 inline Bezier portion(const Bezier & a, double from, double to) {
 265     //TODO: implement better?
 266     Coord res[a.order()+1];
 267     if(from == 0) {
 268         if(to == 1) { return Bezier(a); }
 269         subdivideArr(to, &a.c_[0], res, NULL, a.order());
 270         return Bezier(res, a.order());
 271     }
 272     subdivideArr(from, &a.c_[0], NULL, res, a.order());
 273     if(to == 1) return Bezier(res, a.order());
 274     Coord res2[a.order()+1];
 275     subdivideArr((to - from)/(1 - from), res, res2, NULL, a.order());
 276     return Bezier(res2, a.order());
 277 }
 278
 279 // XXX Todo: how to handle differing orders
 280 inline std::vector<Point> bezier_points(const D2<Bezier > & a) {
 281     std::vector<Point> result;
 282     for(unsigned i = 0; i <= a[0].order(); i++) {
 283         Point p;
 284         for(unsigned d = 0; d < 2; d++) p[d] = a[d][i];
 285         result.push_back(p);
 286     }
 287     return result;
 288 }
 289
 290 inline Bezier derivative(const Bezier & a) {
 291     //if(a.order() == 1) return Bezier(0.0);
 292     if(a.order() == 1) return Bezier(a.c_[1]-a.c_[0]);
 293     Bezier der(Bezier::Order(a.order()-1));
 294
 295     for(unsigned i = 0; i < a.order(); i++) {
 296         der.c_[i] = a.order()*(a.c_[i+1] - a.c_[i]);
 297     }
 298     return der;
 299 }
 300
 301 inline Bezier integral(const Bezier & a) {
 302     Bezier inte(Bezier::Order(a.order()+1));
 303
 304     inte[0] = 0;
 305     for(unsigned i = 0; i < inte.order(); i++) {
 306         inte[i+1] = inte[i] + a[i]/(inte.order());
 307     }
 308     return inte;
 309 }
 310
 311 inline OptInterval bounds_fast(Bezier const & b) {
 312     return Interval::fromArray(&b.c_[0], b.size());
 313 }
 314
 315 //TODO: better bounds exact
 316 inline OptInterval bounds_exact(Bezier const & b) {
 317     return bounds_exact(b.toSBasis());
 318 }
 319
 320 inline OptInterval bounds_local(Bezier const & b, OptInterval i) {
 321     //return bounds_local(b.toSBasis(), i);
 322     if (i) {
 323         return bounds_fast(portion(b, i->min(), i->max()));
 324     } else {
 325         return OptInterval();
 326     }
 327 }
 328
 329 inline std::ostream &operator<< (std::ostream &out_file, const Bezier & b) {
 330     for(unsigned i = 0; i < b.size(); i++) {
 331         out_file << b[i] << ", ";
 332     }
 333     return out_file;
 334 }
 335
 336 }
 337 #endif //SEEN_BEZIER_H
 338
 339 /*
 340   Local Variables:
 341   mode:c++
 342   c-file-style:"stroustrup"
 343   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 344   indent-tabs-mode:nil
 345   fill-column:99
 346   End:
 347 */
 348 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :