source/interpolate/hermitespline.h

   1 #ifndef MSP_INTERPOLATE_HERMITESPLINE_H_
   2 #define MSP_INTERPOLATE_HERMITESPLINE_H_
   3
   4 #include <vector>
   5 #include "spline.h"
   6
   7 namespace Msp {
   8 namespace Interpolate {
   9
  10 template<typename T>
  11 struct Hermite
  12 {
  13         static const Polynomial<T, 3> &h00();
  14         static const Polynomial<T, 3> &h10();
  15         static const Polynomial<T, 3> &h01();
  16         static const Polynomial<T, 3> &h11();
  17         static Polynomial<T, 3> make(T, T, T, T);
  18 };
  19
  20 /**
  21 Implements a cubic Hermite spline.  It is parametrized by the value and slope
  22 at each knot.  The first derivative is guaranteed to be continuous, making it
  23 particularly useful as an interpolator.
  24 */
  25 template<typename T, unsigned N = 1>
  26 class HermiteSpline: public Spline<T, 3, N>
  27 {
  28 public:
  29         using typename Spline<T, 3, N>::Value;
  30
  31         struct Knot: Spline<T, 3, N>::Knot
  32         {
  33                 Value dy;
  34
  35                 Knot(T x_, Value y_, Value d): Spline<T, 3, N>::Knot(x_, y_), dy(d) { }
  36         };
  37
  38         HermiteSpline(const std::vector<Knot> &);
  39         HermiteSpline(Value, Value, Value, Value);
  40 };
  41
  42 template<typename T>
  43 inline const Polynomial<T, 3> &Hermite<T>::h00()
  44 {
  45         static Polynomial<T, 3> p(LinAl::Vector<T, 4>(2, -3, 0, 1));
  46         return p;
  47 }
  48
  49 template<typename T>
  50 inline const Polynomial<T, 3> &Hermite<T>::h10()
  51 {
  52         static Polynomial<T, 3> p(LinAl::Vector<T, 4>(1, -2, 1, 0));
  53         return p;
  54 }
  55
  56 template<typename T>
  57 inline const Polynomial<T, 3> &Hermite<T>::h01()
  58 {
  59         static Polynomial<T, 3> p(LinAl::Vector<T, 4>(-2, 3, 0, 0));
  60         return p;
  61 }
  62
  63 template<typename T>
  64 inline const Polynomial<T, 3> &Hermite<T>::h11()
  65 {
  66         static Polynomial<T, 3> p(LinAl::Vector<T, 4>(1, -1, 0, 0));
  67         return p;
  68 }
  69
  70 template<typename T>
  71 inline Polynomial<T, 3> Hermite<T>::make(T p0, T m0, T p1, T m1)
  72 {
  73         return h00()*p0 + h10()*m0 + h01()*p1 + h11()*m1;
  74 }
  75
  76
  77 template<typename T, unsigned N>
  78 inline HermiteSpline<T, N>::HermiteSpline(const std::vector<Knot> &k):
  79         Spline<T, 3, N>(k.front())
  80 {
  81         typedef SplineValue<T, N> SV;
  82
  83         this->reserve(k.size()-1);
  84         for(unsigned i=1; i<k.size(); ++i)
  85         {
  86                 T dx = k[i].x-k[i-1].x;
  87                 Polynomial<T, 1> t(LinAl::Vector<T, 2>(T(1)/dx, -k[i-1].x/dx));
  88                 Polynomial<T, 3> p[N];
  89                 for(unsigned j=0; j<N; ++j)
  90                         p[j] = Hermite<T>::make(SV::get(k[i-1].y, j), SV::get(k[i-1].dy, j)*dx, SV::get(k[i].y, j), SV::get(k[i].dy, j)*dx)(t);
  91                 this->add_segment(p, k[i].x);
  92         }
  93 }
  94
  95 template<typename T, unsigned N>
  96 inline HermiteSpline<T, N>::HermiteSpline(Value p0, Value m0, Value p1, Value m1):
  97         Spline<T, 3, N>(Knot(0, p0, m0))
  98 {
  99         typedef SplineValue<T, N> SV;
 100
 101         this->reserve(1);
 102         Polynomial<T, 3> p[N];
 103         for(unsigned i=0; i<N; ++i)
 104                 p[i] = Hermite<T>::make(SV::get(p0, i), SV::get(m0, i), SV::get(p1, i), SV::get(m1, i));
 105         this->add_segment(p, T(1));
 106 }
 107 } // namespace Interpolate
 108 } // namespace Msp
 109
 110 #endif