Is there a generalization of the Lagrange polynomial to 3D? What is a way to construct a smooth polynomial surface ($\mathbb{R}^2 \rightarrow \mathbb{R}$) with Lagrange-polynomial properties in every partial derivative? I want to try this for image interpolation.

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).
