source/interpolate/spline.h

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