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