Interpolated Peak Localisation in 3x3 Window

I'm using C++ to create a program to find the sub-pixel location of peaks in an array. Currently I find maxima within a 3x3 window, such that the central pixel in each window is greater than each of it's 8 neighbours.

Is there a well known method for determining the location of the peak to sub-pixel accuracy?

I've read about representing the array by a Taylor expansion up to quadratic terms and taking its derivative at zero to find the offset, but it seems a little heavyweight...

If this seems "heavyweight" for you everything is heavyweight. Generally speaking you need a interpolation algorithm to get from discrete representation to some continuous representation and find the peak. Using "image-processing" implies functions in 2D. I can suggest to use some basic interpolation (linear, bi-linear, cubic, etc.) and find the peaks where derivatives go to 0.

Thanks @Ross. Here's a code snippet of what I wrote to do the job in case anyone else is looking for the same thing.

//
// By approximating with Taylor expansion to quadratic terms, the peak should
// lie at offset [ix,iy] from as calculated by:
//
//   [ix]  =  - [d2I/dx2  d2I/dxy]^-1 . [dI/dx]
//   [iy]       [d2I/dxy  d2I/dy2]      [dI/dy]
//
//
// Assume 'arr' is our array of values (i.e. image) and [x,y] is the location of 
// of a peak pixel in the array.  The interpolated location of the peak is given 
// by the point [x+ix][y+iy].
//

float dx = (arr[x+1][y] - arr[x-1][y]) / 2.f;
float dy = (arr[x][y+1] - arr[x][y-1]) / 2.f;
float dxx = (arr[x+1][y] + arr[x-1][y] - 2 * arr[x][y]);
float dyy = (arr[x][y+1] + arr[x][y-1] - 2 * arr[x][y]);
float dxy = (arr[x+1][y+1] - arr[x+1][y-1] - arr[x-1][y+1] + arr[x-1][y-1]) / 4.f;

float det = 1.f/(dxx*dyy - dxy*dxy);

float ix = x - (dyy*dx - dxy*dy) * det;
float iy = y - (dxx*dy - dxy*dx) * det;