# how to interpolate a 2d function with 6 points?

I'm implementing an algorithm that uses a so called 6-point interpolation, which I never heard before. In the article I'm reading it's described like this:

$\phi(p\Delta x, q\Delta y)=[q(q-1)/2] \phi_{0,-1}+[p(p-1)/2]\phi _{-1,0}$ $+(1-pq-p^2-q^2)\phi_{0,0}+[p(p-2q+1)/2]\phi_{1,0}+[q(q-2p+1)/2]\phi_{0,1}$ $+pq\phi_{1,1}+O(max(\Delta x,\Delta y)^2\phi_{0,0})$

Also, this image should clarify who are $q$ and $p$.

My implementation is running into some problems and I would like to at least exclude this obscure interpolation algorithm as a source of problems. Does anyone know how to obtain it?

-
Rarely does one see an interpolation scheme that's actually discontinuous! Let all $\phi_{i,j}$ be zero except $\phi_{0,-1} = 1$. Then near $(1,\frac12)$, $\phi$ is $\frac18$ on the left side and zero on the right. –  Rahul Oct 13 '11 at 0:04
Considering that the point to be interpolated is always within the upper right square, this doesn't seem like a problem. Right?.. –  Diego Oct 13 '11 at 0:16
Usually one uses such a scheme to interpolate a larger grid of points, in which case I was pointing out that the interpolated $\phi$ won't be continuous from one grid square to the next. If those six points are all you have, then I apologize for misunderstanding the question; but I don't see why you would need a specialized interpolation scheme in that case. –  Rahul Oct 13 '11 at 0:35
oh! It is for a larger grid! I think I see your point. For every point to be interpolated I find the 4 surrounding points and then the last two, according to the current quadrant. –  Diego Oct 13 '11 at 0:58

$+(1+pq-p^2-q^2)\phi_{0,0}$
Notice the sign of $pq$. (If you set $p=q=0$, this term should be $\phi_{0,0}$, for example.)