source/interpolate/hermitespline.h

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