source/interpolate/bezierspline.h

   1 #ifndef MSP_INTERPOLATE_BEZIERSPLINE_H_
   2 #define MSP_INTERPOLATE_BEZIERSPLINE_H_
   3
   4 #include <vector>
   5 #include "spline.h"
   6
   7 namespace Msp {
   8 namespace Interpolate {
   9
  10 template<typename T, unsigned D>
  11 struct Bernstein
  12 {
  13         static Polynomial<T, D> b(unsigned);
  14 };
  15
  16 /**
  17 A spline composed of bezier curves.  Each segment has D-1 control points
  18 between the knots which determine the shape of the curve.  Bezier splines
  19 with 2- or 3-dimensional values are commonly used in graphics.
  20 */
  21 template<typename T, unsigned D, unsigned N = 1>
  22 class BezierSpline: public Spline<T, D, N>
  23 {
  24 public:
  25         using typename Spline<T, D, N>::Knot;
  26
  27         BezierSpline(const std::vector<Knot> &);
  28 };
  29
  30 template<typename T, unsigned D>
  31 inline Polynomial<T, D> Bernstein<T, D>::b(unsigned k)
  32 {
  33         if(k>D)
  34                 throw std::invalid_argument("Bernstein::b");
  35
  36         LinAl::Vector<T, D+1> coeff;
  37
  38         int c = ((D-k)%2 ? -1 : 1);
  39         for(unsigned i=0; i<=D-k; ++i)
  40         {
  41                 coeff[i] = c;
  42                 c = c*-1*static_cast<int>(D-k-i)/static_cast<int>(i+1);
  43         }
  44
  45         unsigned n = 1;
  46         for(unsigned i=0; i<k; ++i)
  47                 n *= D-i;
  48         for(unsigned i=2; i<=k; ++i)
  49                 n /= i;
  50
  51         for(unsigned i=0; i<=D; ++i)
  52                 coeff[i] *= n;
  53
  54         return Polynomial<T, D>(coeff);
  55 }
  56
  57 template<typename T, unsigned D, unsigned N>
  58 inline BezierSpline<T, D, N>::BezierSpline(const std::vector<Knot> &k):
  59         Spline<T, D, N>(k.front())
  60 {
  61         typedef SplineValue<T, N> SV;
  62
  63         if((k.size()-1)%D)
  64                 throw std::invalid_argument("BezierSpline::BezierSpline");
  65
  66         Polynomial<T, D> b[D+1];
  67         for(unsigned i=0; i<=D; ++i)
  68                 b[i] = Bernstein<T, D>::b(i);
  69
  70         for(unsigned i=0; i+D<k.size(); i+=D)
  71         {
  72                 T dx = k[i+D].x-k[i].x;
  73                 Polynomial<T, 1> t(LinAl::Vector<T, 2>(T(1)/dx, -k[i].x/dx));
  74                 Polynomial<T, D> p[N];
  75                 for(unsigned j=0; j<N; ++j)
  76                 {
  77                         for(unsigned c=0; c<=D; ++c)
  78                                 p[j] += b[c]*SV::get(k[i+c].y, j);
  79                         p[j] = p[j](t);
  80                 }
  81                 this->add_segment(p, k[i+D].x);
  82         }
  83 }
  84
  85 } // namespace Interpolate
  86 } // namespace Msp
  87
  88 #endif