Heightfield Interpolation using Cubic Hermite Function

A frequently asked question involves how to read a height value from between the grid points (pixels, samples) in a height field. A height field can be tought of as a square mesh, which has been tesselated and displaced along the square normal according to some input data, typically a grayscale image where white is tall peaks and black is deep canyons. You do this by performing interpolation over the values of the height field.
You need to decide on some interpolator order. The higher order, the smoother the result, but the more expensive the interpolation. Linear or cubic are commonly used. Linear is very cheap, but looks pretty bad (all angles and sharp creases). Cubic is a little more expensive, but looks much better (soft and smooth). Cubic allows you to extract more detail by using multiple samples from within the same "cell" in the height field.
The trick is to first interpolate along the integer grid lines (say, integer X values), to get reference points along the other grid line (say, the Y line) to interpolate between.
For linear interpolation, the pseudo-code looks like so (in pseudo-code):

I(F,x,y) {
  xi = floor(x);
  yi = floor(y);
  xt = x - xi;
  yt = y - yi;
  y0 = F[xi][yi]*(1-xt)+F[xi][yi+1]*xt;
  y1 = F[xi+1][yi]*(1-xt)+F[xi+1][yi+1]*xt;
  return y0*(1-yt)+y1*yt;
}

For cubic, it looks something like (taken from some old terrain code I wrote):

inline float hermite4( float frac, float * ptr )
{
  const float    c     = (ptr[2] - ptr[0]) * 0.5f;
  const float    v     = ptr[1] - ptr[2];
  const float    w     = c + v;
  const float    a     = w + v + (ptr[3] - ptr[1]) * 0.5f;
  const float    b     = w + a;
  return ((((a * frac) - b) * frac + c) * frac + ptr[1]);
}

static float do_sample( unsigned char const * data, uint pitch, 
    uint x, uint z, float xp, float zp )
{
  assert( xp < 1 && zp < 1 );
  assert( xp >= 0 && zp >= 0 );
  assert( x > 0 && z > 0 );
  float xx[ 4 ];
  float yy[ 4 ];
  for( int iz = 0; iz < 4; ++iz ) {
    unsigned char const * data2 = data + x - 1 + (iz + z - 1) * pitch;
    for( int ix = 0; ix < 4; ++ix ) {
      xx[ ix ] = data2[ ix ];
    }
    yy[ iz ] = hermite4( xp, xx );
  }
  float sv = (hermite4( zp, yy ) - 128.f) * (1.f / 128.f);
  return sv;
}