1 #ifndef MSP_INTERPOLATE_SPLINE_H_
2 #define MSP_INTERPOLATE_SPLINE_H_
7 #include "polynomial.h"
10 namespace Interpolate {
12 template<typename T, unsigned N>
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); }
21 struct SplineValue<T, 1>
24 static T get(const Type &v, unsigned) { return v; }
25 static Type make(const T *v) { return *v; }
28 template<typename T, unsigned N>
31 typedef typename SplineValue<T, N>::Type Value;
35 SplineKnot(): x(T()) { }
36 SplineKnot(T x_, const Value &y_): x(x_), y(y_) { }
40 Stores a spline of degree D. It is a piecewise polynomial function with value
41 continuity. Derivatives are not guaranteed to be continuous.
43 The template parameter N controls the dimensionality of the spline's values.
44 Multi-dimensional splines can be used to implement parametric curves.
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.
49 template<typename T, unsigned D, unsigned N = 1>
53 typedef typename SplineValue<T, N>::Type Value;
54 typedef SplineKnot<T, N> Knot;
59 Polynomial<T, D> polynomials[N];
62 Value operator()(T) const;
66 enum { STRIDE = sizeof(Segment)-sizeof(Knot) };
68 unsigned short n_segments;
69 unsigned short capacity;
77 Spline(const Spline &);
78 Spline &operator=(const Spline &);
82 static unsigned data_size(unsigned n) { return sizeof(Knot)+n*STRIDE; }
84 void reserve(unsigned);
85 void add_segment(const Polynomial<T, D> [N], T);
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]; }
94 Value operator()(T) const;
97 template<typename T, unsigned D, unsigned N>
98 inline Spline<T, D, N>::Spline():
104 template<typename T, unsigned D, unsigned N>
105 inline Spline<T, D, N>::Spline(const Knot &k):
108 segments(new char[sizeof(Knot)])
110 new(segments) Knot(k);
113 template<typename T, unsigned D, unsigned N>
114 inline Spline<T, D, N>::Spline(const Spline &s):
122 template<typename T, unsigned D, unsigned N>
123 inline Spline<T, D, N> &Spline<T, D, N>::operator=(const Spline &s)
127 reserve(s.n_segments);
128 std::copy(s.segments, s.segments+data_size(s.n_segments), segments);
129 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