source/interpolate/spline.h

   1 #ifndef MSP_INTERPOLATE_SPLINE_H_
   2 #define MSP_INTERPOLATE_SPLINE_H_
   3
   4 #include <algorithm>
   5 #include "knot.h"
   6 #include "polynomial.h"
   7
   8 namespace Msp {
   9 namespace Interpolate {
  10
  11 template<typename T, unsigned N>
  12 struct SplineValue
  13 {
  14         typedef LinAl::Vector<T, N> Type;
  15         static T get(const Type &v, unsigned i) { return v[i]; }
  16         static Type make(const T *v) { return Type(v); }
  17 };
  18
  19 template<typename T>
  20 struct SplineValue<T, 1>
  21 {
  22         typedef T Type;
  23         static T get(const Type &v, unsigned) { return v; }
  24         static Type make(const T *v) { return *v; }
  25 };
  26
  27 template<typename T, unsigned N>
  28 struct SplineKnot
  29 {
  30         typedef typename SplineValue<T, N>::Type Value;
  31         T x;
  32         Value y;
  33
  34         SplineKnot(): x(T()) { }
  35         SplineKnot(T x_, const Value &y_): x(x_), y(y_) { }
  36 };
  37
  38 /**
  39 Stores a spline of degree D.  It is a piecewise polynomial function with value
  40 continuity.  Derivatives are not guaranteed to be continuous.
  41
  42 The template parameter N controls the dimensionality of the spline's values.
  43 Multi-dimensional splines can be used to implement parametric curves.
  44
  45 While this class contains everything needed to store and evaluate splines, it
  46 cannot be used to construct non-empty splines.  See also HermiteSpline.
  47 */
  48 template<typename T, unsigned D, unsigned N = 1>
  49 class Spline
  50 {
  51 public:
  52         typedef typename SplineValue<T, N>::Type Value;
  53         typedef SplineKnot<T, N> Knot;
  54
  55         struct Segment
  56         {
  57                 Knot start_knot;
  58                 Polynomial<T, D> polynomials[N];
  59                 Knot end_knot;
  60
  61                 Value operator()(T) const;
  62         };
  63
  64 private:
  65         enum { STRIDE = sizeof(Segment)-sizeof(Knot) };
  66
  67         unsigned short n_segments;
  68         unsigned short capacity;
  69         char *segments;
  70
  71 public:
  72         Spline();
  73 protected:
  74         Spline(const Knot &);
  75 public:
  76         Spline(const Spline &);
  77         Spline &operator=(const Spline &);
  78         ~Spline();
  79
  80 private:
  81         static unsigned data_size(unsigned n) { return sizeof(Knot)+n*STRIDE; }
  82 protected:
  83         void reserve(unsigned);
  84         void add_segment(const Polynomial<T, D> [N], T);
  85
  86 public:
  87         unsigned degree() const { return D; }
  88         unsigned size() const { return n_segments; }
  89         const Segment &segment(unsigned i) const;
  90         const Knot &knot(unsigned i) const { return i==0 ? segment(0).start_knot : segment(i-1).end_knot; }
  91         const Polynomial<T, D> &polynomial(unsigned i, unsigned j = 0) const { return segment(i).polynomials[j]; }
  92
  93         Value operator()(T) const;
  94 };
  95
  96 template<typename T, unsigned D, unsigned N>
  97 inline Spline<T, D, N>::Spline():
  98         n_segments(0),
  99         capacity(0),
 100         segments(0)
 101 { }
 102
 103 template<typename T, unsigned D, unsigned N>
 104 inline Spline<T, D, N>::Spline(const Knot &k):
 105         n_segments(0),
 106         capacity(0),
 107         segments(new char[sizeof(Knot)])
 108 {
 109         new(segments) Knot(k);
 110 }
 111
 112 template<typename T, unsigned D, unsigned N>
 113 inline Spline<T, D, N>::Spline(const Spline &s):
 114         n_segments(0),
 115         capacity(0),
 116         segments(0)
 117 {
 118         *this = s;
 119 }
 120
 121 template<typename T, unsigned D, unsigned N>
 122 inline Spline<T, D, N> &Spline<T, D, N>::operator=(const Spline &s)
 123 {
 124         if(s.segments)
 125         {
 126                 reserve(s.n_segments);
 127                 std::copy(s.segments, s.segments+data_size(s.n_segments), segments);
 128                 n_segments = s.n_segments;
 129         }
 130         else
 131         {
 132                 delete[] segments;
 133                 n_segments = 0;
 134                 capacity = 0;
 135                 segments = 0;
 136         }
 137
 138         return *this;
 139 }
 140
 141 template<typename T, unsigned D, unsigned N>
 142 inline Spline<T, D, N>::~Spline()
 143 {
 144         delete[] segments;
 145 }
 146
 147 template<typename T, unsigned D, unsigned N>
 148 inline void Spline<T, D, N>::reserve(unsigned n)
 149 {
 150         if(n<capacity)
 151                 return;
 152
 153         char *new_segs = new char[data_size(n)];
 154         if(segments)
 155         {
 156                 std::copy(segments, segments+data_size(n_segments), new_segs);
 157                 delete[] segments;
 158         }
 159         segments = new_segs;
 160         capacity = n;
 161 }
 162
 163 template<typename T, unsigned D, unsigned N>
 164 inline void Spline<T, D, N>::add_segment(const Polynomial<T, D> polys[N], T end_x)
 165 {
 166         if(!segments)
 167                 throw std::logic_error("Spline::add_segment");
 168
 169         if(n_segments==capacity)
 170                 reserve(n_segments+1);
 171         ++n_segments;
 172
 173         Segment &seg = *reinterpret_cast<Segment *>(segments+(n_segments-1)*STRIDE);
 174         std::copy(polys, polys+N, seg.polynomials);
 175         seg.end_knot.x = end_x;
 176         seg.end_knot.y = seg(end_x);
 177 }
 178
 179 template<typename T, unsigned D, unsigned N>
 180 inline const typename Spline<T, D, N>::Segment &Spline<T, D, N>::segment(unsigned i) const
 181 {
 182         return *reinterpret_cast<const Segment *>(segments+i*STRIDE);
 183 }
 184
 185 template<typename T, unsigned D, unsigned N>
 186 inline typename Spline<T, D, N>::Value Spline<T, D, N>::operator()(T x) const
 187 {
 188         if(!n_segments)
 189                 return Value();
 190
 191         const Segment *seg = &segment(0);
 192         for(unsigned i=1; (i<n_segments && x>seg->end_knot.x); ++i)
 193                 seg = &segment(i);
 194
 195         return (*seg)(x);
 196 }
 197
 198 template<typename T, unsigned D, unsigned N>
 199 inline typename Spline<T, D, N>::Value Spline<T, D, N>::Segment::operator()(T x) const
 200 {
 201         T v[N];
 202         for(unsigned i=0; i<N; ++i)
 203                 v[i] = polynomials[i](x);
 204
 205         return SplineValue<T, N>::make(v);
 206 }
 207
 208 } // namespace Interpolate
 209 } // namespace Msp
 210
 211 #endif