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.

If you fix the value $x$ one variable you have the Lagrange interpolation formula
$$ f(x, y) = \sum_{k=1}^N f(x, y_k)\prod_{i \neq k} \frac{y_i - y}{y_i - y_j} $$
For each fixed value $y = y_k$ you can construct a Lagrange polynomial for the function $f(x, y_k)$.
$$ f(x, y_k) = \sum_{k=1}^N f(x_k, y_k)\prod_{i \neq k} \frac{x_i - x}{x_i - x_j} $$
The result is expansion in some kind of bilinear Lagrange basis.
$$ \prod_{i \neq k} \frac{y_i - y}{y_i - y_j}  \prod_{i \neq k} \frac{x_i - x}{x_i - x_j} = \delta(i = k)\delta(i = k)$$
We have constructed one of many possible B-splines.

In theory the Weistestrass approximation theorem says continuous functions on the interval can be uniformly approximated by polynomials.  Using the argument above it makes sense that 2-variable polynomials are dense in continuous functions on the square.  This is a special case of the Stone-Weierstrass Theorem:
$$ \overline{ \mathbb{R}[x,y]} = C^0\big([0,1]^2\big) $$
Using linear interpolation we can always draw some kind of surface connecting your control points and having the value specified, $(x_k, y_k, f(x_k,y_k))$.  The Stone Weierstrass theorem says it can always be a polynomial and the difference between your approximation $f_1(x,y) = \sum a x^m y^n$ to your original function $f$ can be made arbitrarily small (as small as you want).
$$  \sup_{(x,y)\in [0,1]^2} |f_1(x,y) - f(x,y)| < \epsilon $$
It practice it would be really great to know the degree of the polynomials you find and how to build them (find the coefficients).   Splines let you do just that.
