1 #ifndef MSP_INTERPOLATE_SPLINE_H_
2 #define MSP_INTERPOLATE_SPLINE_H_
5 #include "polynomial.h"
8 namespace Interpolate {
10 template<typename T, unsigned N>
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); }
19 struct SplineValue<T, 1>
22 static T get(const Type &v, unsigned) { return v; }
23 static Type make(const T *v) { return *v; }
26 template<typename T, unsigned N>
29 typedef typename SplineValue<T, N>::Type Value;
33 SplineKnot(): x(T()) { }
34 SplineKnot(T x_, const Value &y_): x(x_), y(y_) { }
38 Stores a spline of degree D. It is a piecewise polynomial function with value
39 continuity. Derivatives are not guaranteed to be continuous.
41 The template parameter N controls the dimensionality of the spline's values.
42 Multi-dimensional splines can be used to implement parametric curves.
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.
47 template<typename T, unsigned D, unsigned N = 1>
51 typedef typename SplineValue<T, N>::Type Value;
52 typedef SplineKnot<T, N> Knot;
57 Polynomial<T, D> polynomials[N];
60 Value operator()(T) const;
64 enum { STRIDE = sizeof(Segment)-sizeof(Knot) };
66 unsigned short n_segments;
67 unsigned short capacity;
75 Spline(const Spline &);
76 Spline &operator=(const Spline &);
80 static unsigned data_size(unsigned n) { return sizeof(Knot)+n*STRIDE; }
82 void reserve(unsigned);
83 void add_segment(const Polynomial<T, D> [N], T);
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]; }
92 Value operator()(T) const;
95 template<typename T, unsigned D, unsigned N>
96 inline Spline<T, D, N>::Spline():
102 template<typename T, unsigned D, unsigned N>
103 inline Spline<T, D, N>::Spline(const Knot &k):
106 segments(new char[sizeof(Knot)])
108 new(segments) Knot(k);
111 template<typename T, unsigned D, unsigned N>
112 inline Spline<T, D, N>::Spline(const Spline &s):
120 template<typename T, unsigned D, unsigned N>
121 inline Spline<T, D, N> &Spline<T, D, N>::operator=(const Spline &s)
125 reserve(s.n_segments);
126 std::copy(s.segments, s.segments+data_size(s.n_segments), segments);
127 n_segments = s.n_segments;
140 template<typename T, unsigned D, unsigned N>
141 inline Spline<T, D, N>::~Spline()
146 template<typename T, unsigned D, unsigned N>
147 inline void Spline<T, D, N>::reserve(unsigned n)
152 char *new_segs = new char[data_size(n)];
155 std::copy(segments, segments+data_size(n_segments), new_segs);
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)
166 throw std::logic_error("Spline::add_segment");
168 if(n_segments==capacity)
169 reserve(n_segments+1);
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);
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
181 return *reinterpret_cast<const Segment *>(segments+i*STRIDE);
184 template<typename T, unsigned D, unsigned N>
185 inline typename Spline<T, D, N>::Value Spline<T, D, N>::operator()(T x) const
190 const Segment *seg = &segment(0);
191 for(unsigned i=1; (i<n_segments && x>seg->end_knot.x); ++i)
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
201 for(unsigned i=0; i<N; ++i)
202 v[i] = polynomials[i](x);
204 return SplineValue<T, N>::make(v);
207 } // namespace Interpolate