Is there a generalization of the Lagrange polynomial to 3D?

While ignoring most of the OP's wishes, by far the simplest interpolation between four neighbouring pixels is the bilinear one: $$ P(x,y) = a + b\,x + c\,y + d\,xy $$ Substitute the four neighbouring pixel positions: $$ P(i,j) = a + b\,i + c\,j + d\,ij \\ P(i+1,j) = a + b\,(i+1) + c\,j + d\,(i+1)j \\ P(i,j+1) = a + b\,i + c\,(j+1) + d\,i(j+1) \\ P(i+1,j+1) = a + b\,(i+1) + c\,(j+1) + d\,(i+1)(j+1) $$ Solve for the coefficients $(a,b,c,d)$ and substitute. Or make an educated guess.
Then for $i \le x \le i+1$ and $j \le y \le j+1$ we find: $$ P(x,y) = (x-i)(y-j)P(i+1,j+1) + (x-i)(j+1-y)P(i+1,j) \\ + (i+1-x)(y-j)P(i,j+1) + (i+1-x)(j+1-y)P(i,j) $$ I would most probably stick to this if it were my research (but it isn't).