source/interpolate/spline.h

   1 #ifndef MSP_INTERPOLATE_SPLINE_H_
   2 #define MSP_INTERPOLATE_SPLINE_H_
   3
   4 #include <algorithm>
   5 #include <vector>
   6 #include "knot.h"
   7 #include "polynomial.h"
   8
   9 namespace Msp {
  10 namespace Interpolate {
  11
  12 template<typename T, unsigned N>
  13 struct SplineValue
  14 {
  15         typedef LinAl::Vector<T, N> Type;
  16         static T get(const Type &v, unsigned i) { return v[i]; }
  17         static Type make(const T *v) { return Type(v); }
  18 };
  19
  20 template<typename T>
  21 struct SplineValue<T, 1>
  22 {
  23         typedef T Type;
  24         static T get(const Type &v, unsigned) { return v; }
  25         static Type make(const T *v) { return *v; }
  26 };
  27
  28 template<typename T, unsigned N>
  29 struct SplineKnot
  30 {
  31         typedef typename SplineValue<T, N>::Type Value;
  32         T x;
  33         Value y;
  34
  35         SplineKnot(): x(T()) { }
  36         SplineKnot(T x_, const Value &y_): x(x_), y(y_) { }
  37 };
  38
  39 /**
  40 Stores a spline of degree D.  It is a piecewise polynomial function with value
  41 continuity.  Derivatives are not guaranteed to be continuous.
  42
  43 The template parameter N controls the dimensionality of the spline's values.
  44 Multi-dimensional splines can be used to implement parametric curves.
  45
  46 While this class contains everything needed to store and evaluate splines, it
  47 cannot be used to construct non-empty splines.  See also HermiteSpline.
  48 */
  49 template<typename T, unsigned D, unsigned N = 1>
  50 class Spline
  51 {
  52 public:
  53         typedef typename SplineValue<T, N>::Type Value;
  54         typedef SplineKnot<T, N> Knot;
  55
  56         struct Segment
  57         {
  58                 Knot start_knot;
  59                 Polynomial<T, D> polynomials[N];
  60                 Knot end_knot;
  61
  62                 Value operator()(T) const;
  63         };
  64
  65 private:
  66         enum { STRIDE = sizeof(Segment)-sizeof(Knot) };
  67
  68         unsigned short n_segments;
  69         unsigned short capacity;
  70         char *segments;
  71
  72 public:
  73         Spline();
  74 protected:
  75         Spline(const Knot &);
  76 public:
  77         Spline(const Spline &);
  78         Spline &operator=(const Spline &);
  79         ~Spline();
  80
  81 private:
  82         static unsigned data_size(unsigned n) { return sizeof(Knot)+n*STRIDE; }
  83 protected:
  84         void reserve(unsigned);
  85         void add_segment(const Polynomial<T, D> [N], T);
  86
  87 public:
  88         unsigned degree() const { return D; }
  89         unsigned size() const { return n_segments; }
  90         const Segment &segment(unsigned i) const;
  91         const Knot &knot(unsigned i) const { return i==0 ? segment(0).start_knot : segment(i-1).end_knot; }
  92         const Polynomial<T, D> &polynomial(unsigned i, unsigned j = 0) const { return segment(i).polynomials[j]; }
  93
  94         Value operator()(T) const;
  95 };
  96
  97 template<typename T, unsigned D, unsigned N>
  98 inline Spline<T, D, N>::Spline():
  99         n_segments(0),
 100         capacity(0),
 101         segments(0)
 102 { }
 103
 104 template<typename T, unsigned D, unsigned N>
 105 inline Spline<T, D, N>::Spline(const Knot &k):
 106         n_segments(0),
 107         capacity(0),
 108         segments(new char[sizeof(Knot)])
 109 {
 110         new(segments) Knot(k);
 111 }
 112
 113 template<typename T, unsigned D, unsigned N>
 114 inline Spline<T, D, N>::Spline(const Spline &s):
 115         n_segments(0),
 116         capacity(0),
 117         segments(0)
 118 {
 119         *this = s;
 120 }
 121
 122 template<typename T, unsigned D, unsigned N>
 123 inline Spline<T, D, N> &Spline<T, D, N>::operator=(const Spline &s)
 124 {
 125         if(s.segments)
 126         {
 127                 reserve(s.n_segments);
 128                 std::copy(s.segments, s.segments+data_size(s.n_segments), segments);
 129                 n_segments = s.n_segments;
 130         }
 131         else
 132         {
 133                 delete[] segments;
 134                 n_segments = 0;
 135                 capacity = 0;
 136                 segments = 0;
 137         }
 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