Some bezier spline interpolation code

bezier.hpp

/// this is mix in glsl terms
template <typename T, typename S>
static inline float mix(T a, T b, S t) {
   return (1 - t) * a + t * b;
}

template <typename T, typename S>
static inline T quadratic_bezier(T a, T b, T c, S t) {
   T d = mix(a, b, t);
   T e = mix(b, c, t);
   return mix(d, e, t);
}

template <typename T, typename S>
static inline T cubic_bezier(T a, T b, T c, T d, S t) {
   T e = mix(a, b, t);
   T f = mix(b, c, t);
   T g = mix(c, d, t);
   return quadratic_bezier(e, f, g, t);
}

template <typename T, typename S>
static inline T quintic_bezier(T a, T b, T c, T d, T e, S t) {
   T f = mix(a, b, t);
   T g = mix(b, c, t);
   T h = mix(c, d, t);
   T i = mix(d, e, t);
   return cubic_bezier(f, g, h, i, t);
}

template <typename T, typename S>
static inline T cubic_hermite_spline(T const *p, unsigned np, S t) {
   // every other point in p is a velocity specifying the next
   // and possibly previous control point, that is for N values
   // in p, we actually have N/2 points on the spline, and N/2*2+2
   // control points (including the points in the spline)
   //
   // e.g. we have 8 values in p
   // each p[i] for odd i is an actual point on the spline
   // each p[i] for even i is a velocity that specifies 1 or 2 control
   // points.
   //
   // if p contains 4 values, the whole thing is a cubic spline from
   // p[0] to p[2] with control points
   // p'[1] = p[0]+p[1] and p'[2]=p[2]-p[3]
   //
   // so for one cubic spline segment, one actually needs
   // the two points we are in between plus their velocities
   // the first point's velocity is added to eh position of the point
   // the second point's velocity needs to be subtracted from the
   // position of the point
   assert(np >= 4);
   assert(np % 2 == 0);

   unsigned actual_np = np / 2;
   S f = t * (actual_np - 1);
   int i = int(floor(f));
   S s = fract(f);

   if (t < 0)
      return p[0];

   if (t >= 1)
      return p[np - 2];

   T a = p[i * 2];
   T b = a + p[i * 2 + 1];
   T d = p[i * 2 + 2];
   T c = d - p[i * 2 + 3];

   return cubic_bezier(a, b, c, d, s);
}